matlab fft phase response square wave I am new to phase analysis, recently I have been trying compare the phase between the input signal and output signal of a system. From the code below, I simulated a square wave and plotted the phase, from basic theory, I thought that since this is a simulated signal, I should get near zero values other than the values at each input frequency. So my questions are (based on the following code):

1) Why does Matlab give me pi/2 as phase response of actual input signals, I don't understand the calculation behind it (does it have to do with my signal being sine and not cosine?)

2) If I don't filter out low (near zero) values, the phase response is a sloping line that goes toward very large numbers, and seems to be proportional to the number of samples of FFT I have, is this due to how unwrap works?

3) If I am trying to compare an input signal vs output signal, both of which for the sake of argument looks similar to the square wave I simulated, how do I determine if the system is linear phased? (from theory, the phase delay should be the same therefore, after finding the phase response, should I check if the difference of phase for each input signal between output and input is constant?

4) Also, group delay is considered as the derivative of phase response, and phase delay is the normalization of phase with respect to its corresponding frequency, a perfect square wave should have constant values at each harmonic for group delay and phase delay shouldn't it?, if so can I somehow see this result in my simulation?

close all;

clear all;

fs = 16000;

end_time = 6-1/fs;

x = 0:1/fs:end_time;

x = x';

n=2000;

freq = 32;

wave = sin(2*pi* freq*x)/(4*pi);

for i = 3:2:n

    wave = wave + sin(i*2*pi*freq*x)/(4*i*pi);

end

wave = awgn(wave,55,'measured');

wave (1:10000) = 0;

NFFT = length(wave);

fft_result = fft(wave, NFFT);

for i = 1:NFFT

    if abs(fft_result(i)) < 10^-4

        fft_result(i) = 0;

    end

end

fft_mag = db(abs(fft_result));

fft_phase = unwrap(angle(fft_result));

fstep = fs/NFFT;

f_range = 0:fs/NFFT:fs-1/NFFT;

figure

plot(f_range, fft_phase);

Thank you for being patient enough to read all of this and giving me answers.
=)

Re question (1):

Most FFTs, or the way atan2() is commonly used with complex FFT results to compute a phase, produce a zero phase (strictly real) output for strictly symmetric input.  This means you want cosine waves, periodic in aperture, not sine waves starting at zero at the window edge, for FFT results with zero phase.


How to Implement the moving average filter in time domain in MATLAB? I have a sinusoidal x(t) = sin(500πt) signal  is corrupted by random noise. The corrupted signal is sampled with sampling frequency fs = 5 kHz and passed through a 5-point moving average filter to reduce the noise. I have written the following MATLAB Code for the same:-
fs=5e3; t=0:1/fs:0.02;
x=sin(2*pi*250*t);           % uncorrupted signal

r=0.1*randn(1,length(x));    % random noise
XR=x+r;                      % signal + noise (filter input)

y=0;
% Implement the moving average filter in time domain
for n=5:length(XR),          % M = 5
    y(n) =                   % missing code using the command 'sum'
end

figure(12);
subplot(211); plot(t,XR,t,x,'r'); legend('Corrupted','Original');
xlabel('time (sec)'); ylabel('Amplitude'); grid;

subplot(212); plot(t,y,'k',t,x,'r'); legend('After filtering','Original');
xlabel('time (sec)'); ylabel('Amplitude');grid;

But, I am not sure how to implement the moving average filter in time domain using sum command to get the output (y[n]) in the for loop as mentioned above in the code. Any suggestions would be highly appreciated! Thanks a lot for the help

HINT (really an answer but not providing code so its a hint):
In your setup, a $M$-point moving average filter uses the average of the previous $M-1$ samples and the current sample to compute the current output. So for $n=5$, the output will be $\frac{1}{5}(x[1]+x[2]+x[3]+x[4]+x[5])$. In general, for $n=k$ the output of this $M$-point MA filter is:
$$ \frac{1}{M}\sum_{n=k-M+1}^kx[k]$$


Observation Matrix in Kalman Filter I've been trying to understand the concept of the Kalman Filter. I came across this great article which makes the concept sufficiently clear.

However I could not understand the concept of the matrix $\mathbf{H}$ which I think is the observation matrix. In Equation 8 in the article, the sensors are modeled using this matrix. I've been through a few posts in this site, where I found that the $\mathbf{H}$ matrix is always a combination of 1's and 0's.

Is this always true for any system where I am trying to estimate the value of a parameter like velocity or position?

If you can do it that way, it makes it straight forward to implement the Kalman filter. As an example - consider a constant velocity model in Cartesian coordinates, but the only thing you are measuring is the position (not velocity) and you are measuring the position in a Cartesian coordinate system. Then your $\mathbf{H}$ matrix will only pick off the position value as you've described.

Now consider a Sonar or Radar system - which measure range and bearing rather than x, y, z. Now the measurements take place in a Polar/Spherical coordinate system but it is easier to keep the model for your targets as having constant velocities in the cartesian coordinate system. Now your state propagation matrix is the same as in the previous case, but now your $\mathbf{H}$ matrix needs to convert your Cartesion state coordinates to Polar/Speherical coordinates for your measurements. e,g, $R =\sqrt{x^2+y^2+z^2}$. So now your $\mathbf{H}$ matrix is not simply extracting one of the state variables. In fact your $\mathbf{H}$ is not even a constant value - it varies with the values of the state vector. Thus the transformation in the measurement matrix $\mathbf{H}$ is non-linear, so the Kalman filter is not even optimal anymore. In this case you now get into the Extended Kalman Filter to try and handle the non-linearity in the transformation.

This was meant to give you one specific example. Also, note that there are other ways of handling the Cartesian to Spherical conversion issue.


What sonic artifacts are caused by phase vocoders? When using a phase vocoder to alter the spectral envelope of a signal, the input signal is split up into many (often hundreds of) bands. With many phase vocoders, the output levels these individual bands are often controllable by the user.

What artifacts are caused by filtering these individual bands in a phase vocoder?
If and why, phase vocoders are better suited (and cause less artifacts) than FIR filters (or other methods) for surgical changes to such narrow frequency bands.  For example, are phase vocoders going to cause the same, or less, pass band ´ringing´ than using an FIR filter for the same task?

I understand phase vocoders only up to the point that they use an inverse DFT to resynthesise the sound, though Im not sure exactly how the individual bands are resynthesized, so I dont know what the problems/artifacts will be.

The term "phase vocoder" may encompass more than one resynthesis method.  One method uses a pruned frequency reassignment of the DFT basis vectors, and thus resynthesizes sound using a sparser set of basis vectors with slightly different frequencies than were present in the original DFT analysis windows.

This sparse basis may not be able to recreate a full spectrum wide-band signal, such as produced by a transient, as well as the original non-reassigned full DFT basis.

Phase vocoder analysis also makes the assumption that the spectrum is fairly stationary between FFT frames, an assumption that can also be violated by transients.

One term for the artifacts produced by this type of analysis/resynthesis is transient smearing.


limit to possible noise shaping? I want to do noise shaping in a 100kHz, 16 bit application, so as to shift all quantization noise to the 25khz-50kHz band, with minimal noise in the DC-25kHz band.

I set up mathematica to create a 31-sample error filter kernel via reinforcement learning which works well: After a little learning, I can get about ~16dB enhancement of high frequency noise for about the same amount of reduction in the low frequency band (the central line is the unshaped dither noise level). This is in line with the "Gerzon-Craven" noise shaping theorem.

Now to my problem:

I cannot manage to shape the noise even more even after extensive learning, although the Gerzon-Craven theorem doesnt forbid it. For example, it should be possible to achieve 40 dB reduction in the low band and 40 dB enhancement in the high band.

So is there another fundamental limit I am running into ?

I tried looking into the Shannon noise/sampling/information theorems but after fiddling with that a while, I only managed to derive a single limit from it: the Gerzon-Craven theorem, which seems to be a direct outcome of the Shannon theorem.

Any help is appreciated.

EDIT: more info

First off, the filter kernel that produces the above noise shaping, Note that the most recent sample is on the right side. Numeric values of the BarChart rounded to .01: {-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29} (Not exactly the bar char but produces similar curve)

Another note about error feedback implementation:

I tried two different implementations of the error feedback. First I compared the rounded output sample to the desired value and used this deviation as error. Second I compared the rounded output sample to the (input+error feedback). Although both methods produce quite different kernels, both seem to level off at about the same noise shaping intensity. The data posted here uses the second implementation.

Here is the code that is used to calculate the digitized wave samples. step is the stepsize for rounding. wave is the undigitized waveform (typically just zeros when no signal is applied).

TestWave[kernel_?VectorQ] :=
 Module[{k = kernel, nf, dith, signals, twave, deltas},
  nf = Length@k;
  dith = RandomVariate[TriangularDistribution[{-1, 1}*step], l];
  signals = deltas = Table[0, {l}];
  twave = wave;
  Do[
   twave[[i]] -= k.PadLeft[deltas[[;; i - 1]], nf];
   signals[[i]] = Round[twave[[i]] + dith[[i]], step];
   deltas[[i]] = signals[[i]] - twave[[i]];
   , {i, l}];
  signals
  ]

The reinforcement method:

The "score" is computed by looking at the noise power spectrum. The goal is to minimize noise power in the band DC-25kHz. I am not penalizing noise in the high frequency band, so arbitrarily high noise there would not diminish score. I am introducing noise in the kernel weights to learn. Maybe, therefore, I am in a (very broad and deep) local minimum, but I consider this extremely unlikely.

Comparison to standard filter design:

Mathematica allows to generate filters iteratively. These can have much better contrast than 36 dB when their frequency response is plotted; up to 80-100 dB. Numeric values: {0.024, -0.061, -0.048, 0.38, -0.36, -0.808, 2.09, -0.331, -4.796, 6.142, 3.918, -17.773, 11.245, 30.613, -87.072, 113.676, -87.072, 30.613, 11.245, -17.773, 3.918, 6.142, -4.796, -0.331, 2.09, -0.808, -0.36, 0.38, -0.048, -0.061, 0.024}

However, when applying those in the actual noise shaping, they are (a) clamped to the same ~40dB contrast, (b) perform worse than the learned filter doing actually no noise attenuation.

Basic dithering without noise shaping

Basic dithered quantization without noise shaping works like this:

Figure 1. Basic dithered quantization system diagram. Noise is zero-mean triangular dither with a maximum absolute value of 1. Rounding is to nearest integer. Residual error is the difference between output and input, and is calculated for analysis only.

The triangular dither increases the variance of the resulting residual error by a factor of 3 (from $\frac{1}{12}$ to $\frac{1}{4}$) but decouples the mean and variance of the net quantization error from the value of the input signal.  That means that the net error signal is uncorrelated with the input but higher moments are not decoupled, so it's not truly completely independent random error, but no one has determined that people can hear any dependency of the higher moments in the net error signal on the input signal in an audio application.

With independent additive residual error we would have a simpler model of the system:

Figure 2. Approximation of basic dithered quantization. Residual error is white noise.

In the approximate model the output is simply input plus independent white noise residual error.

Dithering with noise shaping

I can't read Mathematica very well so instead of your system I'll analyze the system from Lipshitz et al. "Minimally audible noise shaping" J. Audio Eng. Soc., Vol.39, No.11, November 1991:

Figure 3. Lipshitz et al. 1991 system diagram (adapted from their Fig. 1). Filter (italicized in the text) includes in it a one sample delay so that it can be used as an error feedback filter. Noise is triangular dither.

If the residual error is independent from current and past values of signal A, we have a simpler system:

Figure 4. An approximate model of the Lipshitz et al. 1991 system. Filter is the same as in Fig. 3 and includes in it a one sample delay. It is no longer used as a feedback filter. Residual error is white noise.

In this answer I will work with the more easily analyzed approximate model (Fig. 4). In the original Lipshitz et al. 1991 system, Filter has a generic infinite impulse response (IIR) filter form that covers both IIR and finite impulse response (FIR) filters. In the following we will assume that Filter is a FIR filter, as I believe based on my experiments with your coefficients that that is what you have in your system. The transfer function of Filter is:

$$H_\mathit{Filter}(z) = -b_1z^{-1} - b_2z^{-2} - b_3z^{-3} - \ldots$$

The factor $z^{-1}$ represents a one-sample delay. In the approximate model there is also a direct summation path to output from residual error. This gets summed with the negated output of Filter, forming the full noise shaping filter transfer function:

$$H(z) = 1 - H_\mathit{Filter}(z) = 1 + b_1z^{-1} + b_2z^{-2} + b_3z^{-3} + \ldots.$$

To go from your Filter coefficients, which you list in order $\ldots, -b_3, -b_2, -b_1$, to the full noise shaping filter transfer function polynomial coefficients $1, b_1, b_2, b_3, \ldots$, the sign of the coefficients is changed to account for the negation of Filter output in the system diagram, and the coefficient $b_0 = 1$ is appended to the end (by horzcat in the Octave script below), and finally the list is reversed (by flip):

pkg load signal
b = [-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29];
c = flip(horzcat(-b, 1));
freqz(c)
zplane(c)

The script plots the magnitude frequency response and the zero locations of the full noise shaping filter:

Figure 5. Magnitude frequency response of the full noise-shaping filter.

Figure 6. Z-plane plot of poles ($\times$) and zeros ($\circ$) of the filter. All the zeros are inside the unit circle, so the full noise-shaping filter is minimum-phase.

I think the problem of finding the filter coefficients can be reformulated as the problem of designing a minimum-phase filter with a leading coefficient of 1. If there are inherent limitations to the frequency response of such filters, then these limitations are transferred to equivalent limitations in noise shaping that uses such filters.

Conversion from all-pole design to minimum-phase FIR

A procedure for design of different but in many ways equivalent filters are described in Stojanović et al., "All-Pole Recursive Digital Filters Design Based on Ultraspherical Polynomials", Radioengineering, vol 23, no 3, September 2014. They calculate denominator coefficients of the transfer function of an IIR all-pole low-pass filter. Those always have a leading denominator coefficient of 1 and have all poles inside the unit circle, a requirement of stable IIR filters. If those coefficients are used as the coefficients of the minimum-phase FIR noise shaping filter, they will give an inverted high-pass frequency response compared to the low-pass IIR filter (transfer function denominator coefficients becoming numerator coefficients). In your notation one set of coefficients from that article is {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, which could be tested for the noise shaping application although it isn't exactly to the specification:

Figure 7. Magnitude frequency response of the FIR filter using coefficients from Stojanović et al. 2014.

Figure 8. Pole-zero plot of the FIR filter using coefficients from Stojanović et al. 2014.

The all-pole transfer function is:

$$H(z) = \frac{1}{1 + a_1z^{-1} + a_2z^{-2} + a_3z^{-3} + \ldots}$$

So, you can design a stable all-pole IIR low-pass filter and use the $a$ coefficients as $b$ coefficients to get a minimum-phase high-pass FIR filter with a leading coefficient of 1.

To design an all-pole filter and to convert that to a minimum-phase FIR filter, you will not be able to use IIR filter design methods that start from an analog prototype filter and map the poles and zeros into the digital domain using bilinear transform. That includes cheby1, cheby2, and ellip in Octave and Python's SciPy. These methods will give zeros away from z-plane origin so the filter will not be of the required all-pole type.

Answer to the theoretical question

If you don't care how much noise there will be at frequencies above quarter of the sampling frequency, then Lipshitz et al. 1991 addresses your question directly:

  For such weighting functions, which go to zero over part of the band, there is no theoretical limit to the weighted noise-power reduction obtainable from the circuit of Fig. 1. This would be the case if, for example, one assumes that the ear has zero sensitivity between, say, 20 kHz and the Nyquist Frequency, and chooses the weighting function to reflect this fact.

Their Fig 1. shows a noise shaper with a generic IIR filter structure with both poles and zeros, so different to the FIR structure that you have at the moment, but what they say applies also to that, because a FIR filter impulse response can be made arbitrarily close to the impulse response of any given stable IIR filter.

Octave script for filter design

Here is an Octave script for coefficient calculation by another method that I think is equivalent to the Stojanovici et al. 2014 method parameterized as $\nu=0$ with the right choice of my dip parameter.

pkg load signal
N = 14; #number of taps including leading tap with coefficient 1
att = 97.5; #dB attenuation of Dolph-Chebyshev window, must be positive
dip = 2; #spectrum lift-up multiplier, must be above 1
c = chebwin(N, att);
c = conv(c, c);
c /= sum(c);
c(N) += dip*10^(-att/10);
r = roots(c);
j = (abs(r(:)) <= 1);
r = r(j);
c = real(poly(r));
c .*= (-1).^(0:(N-1)); #if this complains, then root finding has probably failed
freqz(c)
zplane(c)
printf('%f, ', flip(-c(2:end))), printf('\n'); #tobalt's format

It starts with a Dolph-Chebyshev window as the coefficients, convolves it with itself to double the transfer function zeros, adds to the middle tap a number that "lifts up" the frequency response (considering the middle tap as being at zero time) so that it is everywhere positive, finds the zeros, removes zeros that are outside the unit circle, converts the zeros back to coefficients (the leading coefficient from poly is always 1), and flips the sign of every second coefficient to make the filter high-pass. Results from (an older but almost equivalent version of) the script look promising:

Figure 9. Magnitude frequency response of the filter from (an older but almost equivalent version of) the above script.

Figure 10. Pole-zero plot of the filter from (an older but almost equivalent version of) the above script.

The coefficients from (an older but almost equivalent version of) the above script in your notation: {0.357662, -2.588396, 9.931419, -26.205448, 52.450624, -83.531276, 108.508775, -116.272581, 102.875781, -74.473956, 43.140431, -19.131434, 5.923468}. The numbers are large which could lead to numerical problems.

Octave implementation of noise shaping

Finally, I did my own implementation of noise shaping in Octave and don't get problems like you did. Based on our discussion in comments, I think the limitation in your implementation was that the noise spectrum was evaluated using a rectangular window a.k.a. "no windowing", which spilled the high-frequency spectrum to the low frequencies.

pkg load signal
N = length(c);
M = 16384; #signal length
input = zeros(M, 1);#sin(0.01*(1:M))*127;
er = zeros(M, 1);
output = zeros(M, 1);
for i = 1:M
  A = input(i) + er(i);
  output(i) = round(A + rand() - rand());
  for j = 2:N
    if (i + j - 1 <= M)
      er(i + j - 1) += (output(i) - A)*c(j);
    endif
  endfor
endfor
pwelch(output, max(nuttallwin(1024), 0), 'semilogy');

Figure 11. Quantization noise spectral analysis from the above Octave implementation of noise shaping for constant zero input signal. Horizontal axis: Normalized frequency. Black: no noise shaping (c = [1];), red: your original filter, blue: the filter from section "Octave script for filter design".

Figure 12. Time domain output from the above Octave implementation of noise shaping for constant zero input signal. Horizontal axis: sample number, vertical axis: sample value. Red: your original filter, blue: the filter from section "Octave script for filter design".

The more extreme noise shaping filter (blue) results in very large quantized output sample values even for zero input.


Books for learning signals and systems? Please recommend books for learning/studying signals and systems
I know about oppenhiem,please write about other books that are easy to understand

Linear Signals and systems by Lathi

Check out MIT lectures of Prof. Oppenheim too. Those were super helpful.


The selection of the amplitude of the signal, consisting of rectangular pulses, at different frequencies I already asked a similar question, but did not receive an answer and did not find a solution, I want to make another attempt to get to the bottom of the truth.

There is a signal consisting of rectangular pulses at different frequencies and other noises. The amplitude of the rectangular pulses may vary. The duration of the pulses can be either 100 μs, or vary by 150-200 pulses per 200 μs.

This is a real signal, which consists of pulses at a frequency of 1 kHz, pulses at a frequency of 1.4 kHz and other noises.

How can one distinguish the amplitude of these pulses and find out the duration? If more details, I need to find:

The amplitude of the pulses is 1 kHz. (the pulses themselves are not important, I just need the average amplitude of each pulse)
The pulse amplitude of 1.4 kHz. (the pulses themselves are not important, I just need the average amplitude of each pulse)
The duration of 1 kHz impulses (the exact time is not required, I just need “yes” if the duration is 100 µs, “no” if the duration is 200 µs.)
The duration of 1.4 kHz impulses (the exact time is not needed, I need “yes” if the duration is 100 µs, “no” if the duration is 200 µs.)

I thought that you can use two band IIR filters, for example 6 orders of magnitude at a frequency of 1 kHz and a frequency of 1.4 kHz, respectively. But after the filter, the amplitude and duration of the initial rectangular signal are lost, and I cannot calculate the initial amplitude and duration of the pulses.

I saw a device that somehow solves this problem on a simple microcontroller, when applying such a signal to the input, at the output I see the amplitudes of both signals, so I concluded that somehow it can be done programmatically. For a month I have been trying to find a solution to this problem on the Internet in both my native language and internationally, but my searches have come to nothing. Perhaps someone solved a similar problem and simply did not share the solution, or no one came across a similar task ... I re-read a lot of pages of “stackexchange” in search of a solution to the problem, but I couldn’t find the answer. I will be glad of any help in solving this problem!

Here is an attempt to simulate this signal in a matlab:

Fs = 100000;
t = 0:1/Fs:0.01;
F = 1.5;
f0 = 1000;
tau = 0.0001;
s = (square(2*pi*t*f0, f0*tau*100) + 1)* F/2;
s = s+randn(size(t))/50;

F2 = 2;
f02 = 1400;
tau2 = 0.0001;
s2 = (square(2*pi*t*f02, f02*tau2*100) + 1)* F2/2;
s2 = s2+randn(size(t))/50;

a=s+s2;

Next, I created mask “b” and tried to use circular cross correlation for verification:

b=0;
for i=1:10
    b(i)=1;
end
for i=11:100
    b(i)=0;
end
b=[b,b,b,b,b];

xcorr = fft(a).*conj(fft(b,length(a)));

At the output, I got this picture.

Here is the result of circular cross correlation with no bias of the original signal, with bias of 30:

Next, on the threshold, I find the beginning of the impulse, and did the integration over 100 μs and a dump:

temp=0;
i=1;
corrf=[];
for g=1:1000
   if (corr(g)>42) || (temp==1)
      temp=1;
      corrf(i)=corr(g);
      i=i+1;
    end
end
corrf=corrf(1:900);
intdump(corrf,10);

The following is what I believe to be an optimized approach for performance in the presence of additive white noise when no other information is known about quantity of pulses or their amplitude distributions, beyond that they are 100 or 200 us long rectangular pulses and repeat at the 1KHz and 1.4KHz rates. This can be even further improved if any other further information is known (including learning from the data to see if anything else about the waveform is repetitive and therefore can be predicted once identified).

Note that the composite waveform consisting of the two pulse types combined  (pulse type meaning 1KHz and 1.4KHz pulses) repeats every 5 ms.

Create a mask to establish the timing of the pulse types. If the starting time between the two pulses is fixed and known then you can create a single reference waveform that is 0 everywhere except where a pulse could be located, and 1 in those positions using 100 us assumed durations. The longer the mask the better the result up to the maximum expected duration of a burst of pulses but keep it in 5 ms durations to take advantage of the cyclical property for using circular cross correlation.

If the 100 and 200 us pulses are equiprobable then in would actually be slightly better to use 200 us pulses in the mask (the break-even is when the ratio of 200 us pulses to 100 us pulses is $\sqrt{2}$, since 200 us masks would have 3 dB more noise which is offset by the signal increase if the mask allows the accumulation of $\sqrt{2}$ more pulses.)

Perform a circular cross correlation of the mask with a capture of the data. The result will have cyclical peaks showing the starting position of the 5ms composite pattern.

The circular cross correlation is given by the following which could be used as a quick test of the procedure:

$$XCORR = ifft(fft(b)conj(fft(a))$$

Where fft(a) is the fft of your mask and fft(b) is the fft of your data.

Here is a very simple example of that to demonstrate how this is done properly and the expected result since the OP is having trouble with this step:

>> mask = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];
>> // create three cases of test data for example
>> dataAlligned = mask;
>> dataShiftOne = circshift(mask,1,2);
>> dataShiftTwo = circshift(mask,2,2);

Thus each would be:

mask = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];

dataAlligned = [1 0 0 0 0 1 0 0 0 0 1 0 0 0 0];

dataShiftOne = [0 1 0 0 0 0 1 0 0 0 0 1 0 0 0];

dataShiftTwo = [0 0 1 0 0 0 0 1 0 0 0 0 1 0 0];

The circular cross correlation is done for each case above as:

>> corr = @(data,mask) ifft(fft(data).*conj(fft(mask)));

Below shows the result of this for each case:

>> corr(dataAllgined, mask)
ans =
   3   0   0   0   0   3   0   0   0   0   3   0   0   0   0
>> corr(dataShiftOne, mask)
ans =
   0   3  -0   0   0   0   3  -0   0   0   0   3  -0   0   0
>> corr(dataShiftTwo, mask)
ans =
   0   0   3   0  -0   0   0   3   0  -0   0   0   3   0  -0

Given the sparse nature of the mask, the implementation for the circular cross correlation would likely be more efficient in the time domain by repeated multiply and sums at all circular rotations of the mask. Given the mask only has values of 0 and 1 this would be quite simple to compute even in a microcontroller, as the multiply operation would just be a selection of which samples to sum.

If the starting delay between the two pulse types is not known, then do the above with two masks; one for the 1KHz pulses and one for the 1.4KHz pulses.

Once timing is established for the above techniques, then we know when a pulse of each type should appear if a pulse is to be present in that interval. So at this point to detect each pulse use an integrate and dump that would integrate over 100us in two successive durations at all the expected start times. (This can be implemented with a moving average filter that has a width equal to the pulse duration, and taking the samples of the moving average output at each expected end of pulse location.)  A decision metric is established using a noise optimized threshold on the integrate and dump output to decide pulse present=1 or pulse not present = 0. The use of two successive durations is to determine if the pulse was 100us long or 200 is long (decoding 1,1 or 1,0 for that interval).

The magnitude of each pulse is the magnitude of the integrate and dump output scaled by the number of samples—- keep a running average of the detected magnitude for each pulse type and adjust the threshold based on the averaged magnitude to optimize the decision for pulse present or not present.


How do I compare two resultant images from different binarization algorithms? I have a binarized image to which i add noise and then try to filter out the noise using various thresholding algorithms such as otsu and niblack.How can i compare the resultant image with the original image so as to find the percentage error that exists between the two??

the original image is as such:

and the resultant image is:

I need a way to find the percentage error that is present.

Well i think you need to compare the effectiveness of various methods of noise removal that you would be using to get back the original image :

I would suggest to take the following parameters for measure

Mean Square Error (MSE)
Peak Signal to Noise Ratio (PSNR in dB)
Structural Content (SC)

and many others

You can find the MATLAB codes out here


Identify fundamental frequencies in short-period waveform Maybe I'm being stupid and failing to see something here...

I need to analyse a signal over a relatively short time period (~1-2 cycles for the lowest frequency, perhaps 20-30 cycles for the highest frequency). My aim is to identify the frequencies which are present - usually two or three, never more than six. So naturally, I perform a Fourier transform, applying a Blackman-Harris window to get a reasonable compromise between peak breadth and sidelobe suppression.

A priori I expected that when there were two frequencies I'd see two peaks, with three frequencies I'd see three peaks etc. What actually happens is that with two frequencies I get three peaks ($f_{1}$, $f_{2}$, $f_{1}+f_{2}$), with three frequencies I get up to seven peaks ($f_{1}$, $f_{2}$, $f_{3}$, $f_{1}+f_{2}$, $f_{2}+f_{3}$, $f_{1}+f_{3}$, $f_{1}+f_{2}+f_{3}$) and so on up to 63 peaks with six frequencies. This isn't terribly useful. Because the peaks are relatively broad because of the short sampling time, it is almost always the case that some of them overlap and are indistinguishable from one another, and even if (in the six-frequency case) I could see all 63 peaks, there are $ ^{63}C_{6} = 67,945,521$ possible ways of guessing which six peaks represent the actual frequencies I want!

Yet I can't imagine that what I'm trying to do is that unusual, so I can't help thinking that I'm either doing the transform wrong, or missing some other, much simpler way of achieving my objective. Which is it?

Finding the sum (and/or difference) of source frequencies in an FFT result usually means that there is an inter-modulation distortion in the source mixing (mixer).  Likely a multiplication instead of a strictly linear summing of source signals.

Regarding your last bit of thinking, trying to do blind source separation using just a short FFT is trying to do something unusual.  Any DSP solution is extremely unlikely to be simpler.


Bilinear Interpolation on A Subsection Of An Image Is it possible to do bilinear interpolation on each 8 by 8 block of an image one at a time until you cover the entire image rather than doing it on the entire image? What are the consequences of this?

Bilinear interpolation can be viewed as traditional sample rate conversion (insert 0s between existing samples, convolve with a triangular kernel in each dimension). Somewhat unusually, the size of the (continous-time) triangular filter kernel is not scaled relative to resampling ratio before it is sampled; each sample is always the result of 4 neighbouring input sample (unless it falls exactly on one of them). This makes it quirky seen with dsp goggles but maps very efficiently to compute hw.

When doing convolution you need overlap between separate chunks, or you will get some kind of discontinuity. For bilinear interpolation the overlap need only be one sample wide.


BER for OFDM-QPSK over Rayleigh Channel The rayleigh channel is a single tap.
   1.
    h=1/sqrt(2)*randn(1,length(txSig_cp))+1i*randn(1,length(txSig_cp));
    fading=h.*txSig_cp;

   2.
    h=1/sqrt(2)*(randn(1,1)+1i*randn(1,1));
    fading=conv(h,txSig_cp);

Which one is correct? or both also wrong. Because I have obtained different results when compared to the berfading function in MATLAB (berfading(EbNoVec,'psk',4,1)).

w = conv(u,v) returns the convolution of vectors u and v. If u and v
are vectors of polynomial coefficients, convolving them is equivalent
to multiplying the two polynomials. (from MATLAB documentation)

Both expressions are correct. In 2., whenever an addend of a convolution operation calls randn in a convolution computation procedure, the callee returns new values for both real and imaginary parts of h.


Inverse Fast Fourier with overlap I had to work with sound(denoising algorithms), when I read about an algorithm it said that we do DTSTFT(Discrete-time short-time fourier transform) we will create an overlap between each "chunk" of data, for example if I have $$s=[1,2,3,4,5,6]$$ I divide it into $$s^{(1)}=[1,2,3,4]\\s^{(2)}=[3,4,5,6]$$ and then I use FFT on $s^{(1)}$ and $s^{(2)}$(this for $50\%$ overlap)

My problem is how can I do IFFT for this after I touched $S^{(1)}$ and $S^{(2)}$?  Everywhere I looked said to do the overlap but had not say a thing about the Inverse

There are several ways.

Start with taking the iFFT of $ S^{(1)} $ and $ S^{(2)} $.  Let's call the results $ m^{(1)} $ and $ m^{(2)} $, since they are modified.  They are back in the time domain and you wish to recombine them into one signal again.  Ideally, they would line up perfectly and you could just truncate one and keep the other.  This won't happen, so the trick is to fade the latter into the former over the overlap region.  At the beginning of the region the weighting of $ m^{(1)} $ should be one and the weighting of $ m^{(2)} $ should be zero.  At the end of the region, the reverse should be true.  The simplest solution is a linear transition, but that may not be the best.

For a 50% overlap, there is a nice property of a VonHann window in that the appropriate weighting is implicitely correct in the transition region and you can simply sum the two modified signals over the overlap region.  For other overlap sizes the same function can still be used, but it needs to be resized appropriately.  The window function can be applied in the time domain before you take the initial FFT, or after the iFFT is done.  Of course, which you do greatly affects the bin values, so it depends on the nature of your "touching" on which way is better.  The equivalent of a VonHann window can also be done in the frequency domain by calculating a new value for each bin by subtracting the average value of its neighbors and taking half that value.  The equation would be $Y[k] = -.25 * X[k-1] + .5 * X[k] - .25 * X[k+1]$ where the $X$s are the original bin values.  Doing this before or after your "touching" is the same as if you were windowing the signal before going into the frequency domain or after coming out.

The less touching you do, the more similar the two modified signals are going to be so the smoother the fading transition.  For too much touching, the mismatch between the two modified signals may be too great and the results won't be so good.  However, with a fading over the overlap region, the results are always going to be a smooth signal.

The VonHann window is good for spectrogram displays because it spreads the value of a sharp peak (near integer frequencies) to the neighboring bins and reduces the leakage of in between frequencies, so when you do a tone sweep through the frequency range, the display remains somewhat consistent and less frame size dependent.  Chances are that this same property will make it desirable for you to do the window function before you do your touching.

Hope this helps,

Ced

=================================================

Followup:

The process works like this:

1) Divide your signal into overlapping chunks

2) For each chunk

A) Apply Window function

B) Do the FFT

C) Apply touches

D) Do the iFFT

E) Reverse Window function

3) Recombine chunks into a new signal

Generally speaking, when you recombine the signal you will want a weighting function that goes from 1 to 0 across your overlap region so you are fading out of one chunk as you are fading into the next chuck.

If you use an overlap of greater than 50%, you will have regions where more than two chunks exist for the same time domain and you'll have to figure out an appropriate recombining strategy.  My understanding is that overlaps over 50% are in applications where the end result is the frequency domain and your interest is tracking frequency components.

You can use any window function you want, and any overlap amount.  If you follow the steps above, if you don't do any "touching" in the frequency domain, you should get your original signal back.  The advantage of the VonHann window using a 50% overlap is that the window function actually works as a weighting function so steps 2E and 3 are significantly simplified.

Check out this article:  "Overlap_Add_OLA_STFT_Processing"

The Hamming window is covered below figure 8.10.


SNR estimation in segmented speech signals I was called to estimate the signal-to-noise ratio in a collection of recordings that were manually annotated, meaning that the speech and noise segments are known.

The noise is uncorrelated with the signal and is additive noise: $x(i)=d(i)+n(i)$ ('d' stands for desired). I have no information regarding the sound intensity (no sound meter was employed) neither the voltage-sound-intensity relationship in devices.

I found a formula here that seems to be a possible solution (authors explain the annotated way before diving into the no boundaries problem):

$$\text{SNR}=10 \log_{10}{\left(\frac{P(x)-P(n)}{P(n)}\right)}$$

Where $P(x)$ and $P(n)$ are the power of the contamined signal 'x' and the noise signal 'n' ('n' is basically a copy of 'x' but with zero energy in the speech segments). The mean was subtracted from both signals as requested in the paper.

When I apply the formula in two separate recordings (one clean and the other containing much noise, different signals and different noise) I am getting negative measurements: -6.250953 vs. -7.793706, and a difference of just 1.5 dB between them (I was expecting to see a 20dB diference).

Does it make sense to have negative numbers? How can I interpret the 1.5 dB difference in a meaningful way? Thanks in advance for any light on this!

Below are two copies of the same magnitude spectrum of the example of the clean recording I mentioned earlier (I did not include the noisy one). The formula estimated a SNR of -6.250953 dB. I did not include the oscillogram but it oscillates in the [-1, 1] range:

This is what I am taking as the $x(i)=d(i)+n(i)$ signal:

And this is what I am taking as the noise $n(i)$ signal (all energy except the black regions which represent the desired signal):

Therefore the formula can be applied very straightforward

The given equation should be used with caution. Not always we can separate powers of two mixed signals just like that. It will require them to be orthogonal (or independent and at least one of them zero-mean).

In your formula the power of signal $y[n]$ is defined as $$P=\frac{1}{N}\sum_{i=1}^{N}|y[i]|^2$$
where $N$ is the length of $y$. You shouldn't forget the normalizing factor $\frac{1}{N}$.

Just to make it clear, $x$ should be a recording of signal plus noise, and $n$ should be just noise (silence/ no speech). The oscillogram of $n$ should look like the background of oscillogram of $x$ for the formula to make sense.

Regarding the negative values, logarithm of a number less than one is negative. It implies that the numerator is smaller than the denumerator, or the signal power is less than that of noise. A negative SNR can be seen in practice, where the signal is weak in nature. GPS signal can be a good example. But I really think that your noise is polluted with some sort of harmonics that don't appear in $x$, that is why the overall result is less than one and the SNR becomes negative.


IDFT of H(z) sampled in N values If a have a causal IIR filter described by $H(z)$ and I sample it in $N$ equispaced values around the unit circle, I get a DFT of $N$ points. That DFT corresponds to $h[n]$ truncated in $n=N-1$ or to the DFT of $h[n]$ periodized with period $N$ (i.e. with copies of itself overlapped every $N$ points)?

The fact that sampling in one domain corresponds to aliasing in the other domain is an important property of the Fourier transform. Of course, this can also be shown formally. If $X[k]$ is the length $N$ DFT of $x[n]$, we have

$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{1}$$

Furthermore, $H(z)$ is related to the impulse response by

$$H(z)=\sum_{m=-\infty}^{\infty}h[m]z^{-m}\tag{2}$$

If we define

$$X[k]=H(e^{j2\pi k/N})=\sum_{m=-\infty}^{\infty}h[m]e^{-j2\pi mk/N}\tag{3}$$

we can determine $x[n]$ as follows. From $(1)$ and $(3)$ we have

$$\begin{align}x[n]&=\frac{1}{N}\sum_{k=0}^{N-1}\sum_{m=-\infty}^{\infty}h[m]e^{-j2\pi mk/N}e^{j2\pi nk/N}\\&=\frac{1}{N}\sum_{m=-\infty}^{\infty}h[m]\sum_{k=0}^{N-1}e^{j2\pi(n-m)k/N}\tag{4}\end{align}$$

The last sum in $(4)$ equals $N$ for $n-m=lN$ for integer $l$, and zero otherwise. Consequently, this last sum can be rewritten as

$$\sum_{k=0}^{N-1}e^{j2\pi(n-m)k/N}=N\sum_{l=-\infty}^{\infty}\delta[n-m-lN]$$

and $(4)$ becomes

$$\begin{align}x[n]&=\sum_{m=-\infty}^{\infty}h[m]\sum_{l=-\infty}^{\infty}\delta[n-m-lN]\\&=\sum_{l=-\infty}^{\infty}\sum_{m=-\infty}^{\infty}h[m]\delta[n-m-lN]\\&=\sum_{l=-\infty}^{\infty}h[n-lN]\tag{5}\end{align}$$

where I've used the fact that the sum over $m$ is the convolution of $h[n]$ with $\delta[n-lN]$, and consequently equals $h[n-lN]$.

Eq. $(5)$ establishes the result that the IDFT of an equidistantly sampled discrete-time Fourier transform (DTFT) equals an aliased version of the sequence corresponding to the given DTFT.


Rationally related frequencies and the Fourier Series representation Suppose that we have the signal $$x(t) = e^{j\omega t} + e^{j\frac{3}{2} \omega t},$$ and we want to find a Fourier Series representation for that signal. Is this possible? According to my understanding, since the ratio of the frequencies of the two complex exponentials above is a rational number, the signal is periodic. Since it is periodic, it is possible to represent the signal as $$x(t) = \sum_{k = -\infty}^{\infty} a_k e^{jk\omega_0t},$$
where $\omega_0$ is the fundamental frequency of the signal?

Yes, that's easily possible. The fundamental frequency of the given signal is $\omega_0=\omega/2$, and its Fourier series coefficients are

$$a_k=\begin{cases}1,\quad k\in\{2,3\}\\0,\quad\text{otherwise}\end{cases}$$


The peak in frequency spectrum from sampled data is higher than the true amplitude I am given a task to do sampling on the function
$$f(t) = 1.3 \sin(2\pi \cdot 10 \cdot t)$$
with sampling frequency 4 Hz and start at time t_start = 0.0877s, so there should be aliasing and in the frequency spectrum, the amplitude should be at the perceived frequency.
The perceived frequency is 2 Hz as the folding frequency.
I chose to take $N=128$ data points.

Then after completing FFT with sampled $f(t)$, I wanted to plot frequency spectrum, however I have a peak of 1.815 at 2 Hz, however the height of peak should be around 1.3. Here I attached this spectrum.

Is there any explanation to this? Usually due to the leakage, the peak should be less than the actual amplitude, however I got higher peak. It seems quite strange for me.
Here is also a Matlab code, if you find it more convenient to check.
t = 0.0877:0.25:32.087;
signal = 1.3.*sin(2.*pi.*10.*t);

y = fft(signal);

%65th data point corresponds to fft at 2 Hz
abs(y(65))/(numel(t)/2)

Thank you in advance.

10 Hz aliasing down to two Hz which is the Nyquist frequency.
When you sample a Nyquist frequency signal, you only get two samples per period and the two samples are exactly 180 degree apart, i.e. you always get the same numbers with alternating sines.
Which number you get depends on the relative phase fo your sample clock and the signal.If you sample a consine you samle at the peaks and you get [1 -1 1 -1 ...] .  If you sample a sine wave , you sample at the zeros and you get [0, 0, 0, 0 ...]. For any other phase, you get something in between.
Keep in mind that the sampling theorem requires the signal the be bandlimit to "less than half the sample rate" not to "equal or less than half the sample rate".


Seeking advice on how to denoise a video from poor CCTV lighting I have grainy CCTV footage that will need some amount of clean up. I request for ideas on the best way to clean out the noise. I am attaching a frame extracted from the video here for hints on the specific noise cleaning I will need. Just for context, the guy in the image robbed and brutalized a vulnerable female a couple of days back. The guy wore a COVID mask thorough out the ordeal.  In the brief spell when he put the mask off, he was in a region where lighting was too poor. Hence the image that would give us a clean shot is way too grainy.

Something pretty basic to do with good generalization would be using one of the following:

The Bilateral Filter.
The Non Local Means Filter.

They are easy to implement and have a decent performance on many models.
Good resource about denoising is - Which Noise Reduction Algorithms Are Used in Commercial RAW Image Processors?


How can I change parameters into a text string in MATLAB? I am trying to compute the Morlet wavelet transform (cmor) of a signal, but I have to do it several times with different parameters Fb and Fc. Is there a way to automate it? For example, I have to do the transform with: cmor1-1, cmor1-2, ..., cmor1-10.
I created a vector Fc=1:10, but I'm not able to replace Fc parameter into cmor1-Fc.
I tried this way:

morlet='cmor1-Fc'
for i=1:10
morlet(1,8:9)=i;
end

but it doesn't work (there appear squares instead of the number, and it would properly work only for 2-digit numbers.
Thank you

You can use this code

text = 'cmor1-';
for i=1:10
     s{i}=[text num2str(i)];
end

where s{i} contains the string cmor1-i


Controllability's dependence on B matrix Assuming the LTI system:

$\dot{x}=Ax+Bu\\$ , where $x\in R^n$

I understand that

the $B$ matrix is usually given, and it is up to us to derive the optimal $u$ given an initial state of the system and a desired state where we want the system to arrive, in a finite time
a full-diagonal $B$ matrix could guarantee global controllability for the system
any $B$ matrix has to be dependent on the system and the application

In case the $B$ matrix is not a priori given, how could we cook up the "ideal" $B$ matrix, or at least any possible matrix that can render the system controllable?

potential ways, that might or might not be correct:

calculating the controllability matrix for any possible combination and keep those combinations that offer ideally global controllability, or at least the maximum rank of the controllability matrix. But this is computationally too demanding
through Cayley-Hamilton theorem, finding the eigenvalues of $A$ matrix and accordingly deciding for the minimum and maximum number of columns of $B$ matrix. But this implies that we know or can find, the values of the components of the $B$ matrix

There are different ways of expressing a control system. For example,

Observable Canonical Form
Controllable Canonical Form

If you express your system in controllable canonical form, then $B$ take the form of a column matrix with all zeros except one in the top position.

Link: https://en.wikibooks.org/wiki/Control_Systems/Standard_Forms


Can anyone suggests some basic readings for video processing? Actually, I am new to this field and would like to work on a project in Video processing. I have a little background in Image processing but would like to study further in video processing. So, If anyone could suggest some readings as pre-requisites for video processing.

I recommend Video Tracking: Theory and Practice  as comprehensive overview of what have been done in this area for last 15 years. Also, 3D Computer Vision: Efficient Methods and Applications for quick understanding useful algorithms.

This books will be a good start, but when you have enough knowledge, you should start looking into research papers, because information there is more updated than in books.


Is matlab fft faster for real input? I realize that matlab fft is twice as fast for real input. But does matlab fft take the fact that the transform will be conjugate symmetric into consideration? Since the input is real, it would have half as much information as a complex input, so this can explain the factor of 2. I was thinking that it perhaps should be a factor of 4. 2 for conjugate symmetry, 2 for having no complex component. Am I thinking wrong?

Matlab, like most computational packages, implement a pretty neat set of fast Fourier transform algorithms referred to as FFTW: fastest Fourier transform in the west. You can find more out about it here.

FFTW actually runs a preflight on your machine when it is first installed. At that time, it runs a bunch of test computations to figure out which prime factorization would be fastest given your processor speeds and cache sizes. Then, it uses a symbolic program to generate a bunch of FFT algorithms specific to your machine for a whole bunch of cases, like real valued, symmetric, etc. So, it's hard to say exactly what's going on under the hood except that if there's an efficiency to be had, FFTW will take advantage.


Why only take 205 sample when using a sample rate of 10khz or 8khz? I've often seen i my DSP book, that when they record a sound, they use a sample rate at 10khz or 8 khz, which makes sense since the fs is 2 times bigger than the sampling signals max freq.
What i don't understand is why they only take 205 samples, rather than 10000 samples or 8000 samples, wouldn't each sample be "unperiodic"??..
What is the magic behind the number 205 ???

EDIT:
I know they aren't able to sample 8000 times here, but why not take 320 samples then?? than just 205??

Example:

  Estimating DTMF Tones with the Goertzel Algorithm

  The minimum duration of a DTMF signal defined by the ITU standard is
  40 ms. Therefore, there are at most 0.04 x 8000 = 320 samples
  available for estimation and detection. The DTMF decoder needs to
  estimate the frequencies contained in these short signals.

  One common approach to this estimation problem is to compute the
  Discrete-Time Fourier Transform (DFT) samples close to the seven
  fundamental tones. For a DFT-based solution, it has been shown that
  using 205 samples in the frequency domain minimizes the error between
  the original frequencies and the points at which the DFT is estimated.

  Nt = 205; original_f = [lfg(:);hfg(:)]  % Original frequencies
  original_f =

     697
     770
     852
     941
    1209
    1336
    1477

  k = round(original_f/Fs*Nt);  % Indices of the DFT estim_f =
  round(k*Fs/Nt)      % Frequencies at which the DFT is estimated
  estim_f =

     702
     780
     859
     937
    1210
    1327
    1483

  To minimize the error between the original frequencies and the points
  at which the DFT is estimated, we truncate the tones, keeping only 205
  samples or 25.6 ms for further processing.

  tones = tones(1:205,:);

The frequency of DFT bin centers is i*Fs/N, for integer i.  Thus different N will create different arrays of frequencies.  For some N, this will produce an array that has the lowest nearest bin center distance from the frequency set that you wish to measure (in thus case, a vector of DFT frequencies).  According to the referenced document, for some standard sample rate, N=205 will produce that vector with the closest DFT bin centers to DFT frequencies.

Since each DFT bin has a Sinc shaped roll-off, the closer a bin frequency is to a desired measurement frequency, the lower the magnitude measurement error will be without needing further interpolation computation.

For N greater than 320, the DFT will have greater frequency resolution, but the time resolution gets worse and will overlap multiple tones, which creates greater ambiguity in which tone is being detected.


External/Internal Description of a System When an LTI system is NOT completely controllable/observable, does the linear differential equation governing the system give the internal or the external description of a system?

External. That's the whole point: The behavior of any LTI system can be completely characterized by it's transfer function, impulse response, differential equation, etc. (which can all be converted into each other).

That's a "black box" description and gives no information on the internal details of the system. There many different physical realizations of the same system, but you can't tell from just the outside behavior


How to detect letters in an image? I have a project that requires detecting and processing letters in a given input image. I am also new to subject of image processing. There are lots of subsubjects in image processing. I believe I don't have to learn all the concepts of image processing to achive my purpose.

Question above; how can I detect letters in a given input image?
Note: I already know matlab. I need to do this by implementing algorithm, not by a built-in function.

I am gonna tell my thoughts. I want to write an application which takes an input image as an argument. If there is are English letters on that image (photo), it will detect them by circling each of the letters. I don't know which subject in digital image processing I must look at.

The field is knows as Optical Character Recognition (OCR) and it's a broad field.

Since it's clearly an academic exercise, you can keep it simple. Start by image segmentation : find shapes where a shape is just a set of connected pixels. [1] Then compare the shapes to known shapes. Don't expect an exact match, do an approximate match. What overall properties can you reliably determine for the shape? Aspect ratio? Holes in the middle?

[1] At this point, you can already draw circles around each shape that's not too big nor too small.


low frequency issue with noise spectrum (of derivative) I have hard time in understanding the results I get when I plot the spectrum of the time-derivative of a noise signal:

$|DFT[\frac{dx(t)}{dt}](\omega)|^2$

as a function of frequency $\omega$ (see image below), where $DFT$ is the discrete Fourier Transform, $x(t)$ is a discrete time serie of random samples extracted from a gaussian distribution.

In python:

import numpy as np
import matplotlib.pyplot as plt
x = np.random.randn(100000)
F = np.abs(fft.fft(np.diff(x)))**2
plt.loglog(F[:50000])

The F.transform of a derivative is $i\omega$ times the F.transform. So I'd expect the result to be  $\omega^2 |DFT[x]|^2$.

This happens, but what are the features at low frequency? If I run the script several times, a plateau shows up at low frequency, but not always with the same amplitude and not always at the same low frequencies (see in the fig. the differences between various runs of the same script). Probably something strange happens also at high frequency? I heard "leakage" can explain things like this, but I'n not an expert. Can you help me to understand, and correct it?

I think your python code is wrong. Below I put a modified and corrected MATLAb/OCTAVE code instead. I hope it's clear from the code about the difference.

 N = 32*256;                       % Length of the time domain signal
 x = randn(1,N);                   % Normal type random vector
 wk = linspace(-pi,pi-2*pi/N,N);   % DFT frequencies from -pi to pi
 X = fftshift(fft(x,N));           % DFT X[k] of x[n] FROM -pi to pi
 Px = (wk.^2).*(abs(X).^2);        % Simulated |DFT|^2 of the derivative signal
 figure,plot(linspace(-1,1,N),Px)  % plot the $DFT|^2 from -pi to pi
 title('Simplified Periodogram');  % QUADRATIC behaviour is recognised.
 xlabel('Normalized Frequency (x \pi rad/sample)'); % No DC plateau ?

 % LOG-LOG plot.
 figure,loglog(10E3*wk(N/2+1:end),Px(N/2+1:end))

 % Alternatively: compute PSD based on time-domain signal
 pX = fftshift(j*wk.*X);            % DFT of the derivative
 dx = real(ifft(pX,N));                % time-domain signal from I-DFT
 figure,periodogram(dx,hamming(N));    % compute & display its Periodogram
 figure,pwelch(dx);                    % compute & display its Welch

The resulting plots are:

The loglog plot of MATLAB is the following:

there is no DC plateau here. I believe it's your diff function which adds the DC content to your differentiated data. You should apply the DSP standard diagnostic procedure to investigate your case and find out which step deviates from the expected behaviour...


How to take acoustic measurements of a clasroom I'm an electrical engineering major and I have a question regarding a project I'm working on.

I'm supposed to design 3 band-stop filters on MATLAB for echo cancellation/suppression.

My problem is not filter designing, but with taking acoustic measurements of a classroom since I have little or no information in that area.

The supervisor suggested taking coefficients using impulse response invariant! I could not find something solid or a method while searching online.

  So my question is, what is the simplest way of taking acoustic
  measurements of a normal size room/classroom? How to find the
  coefficients through impulse response invariant method?

All suggestions are appreciated, current tools I have (dynamic mic, speakers, audio interface (Focusrite Scarlet 2i2), RoomEQ wizard, MATLAB/Simulink, FPGA DE2 board with cyclone)

I tried using RoomEQ wizard and followed their tutorial but did not work properly or maybe I`m not getting the whole picture.

It doesn't matter where your microphone is in the room when you take the measurement. What you need to do is find where the peaks are in the cross correlation (each peak corresponds to an echo at a certain lag/delay), figure out what the corresponding inverse filter needs to be to cancel out the sound.

If you're trying to do echo cancellation in MATLAB and have the signal processing toolbox, take a look at the example in xcorr by typing:

  doc xcorr

for a simple example.

Edit 1:To find the inverse filter required, you need to calculate the transfer function between the source and the mic. (edit: sentence in here was misleading) Once your filters are set up, as long as nothing is changing you should be able to cancel the echoes produced by the room.

Edit 2: To get the room impulse response you typically use a long LFM chirp (some people pop a balloon) and a microphone, then calculate the transfer function between the input signal and the recorded signal. Here's an example paper of what you can do: http://www.artalabs.hr/papers/im-aaaa2003.pdf

Edit 3:

  You do realise that reverberation (or echoes as you described) are dependent on location and even time?

you said you are making an "impulse response invariant" filter. So as long as neither your source or microphone are moving, the IIR filter created from your Room Impulse Response measurement will cancel out your echoes...just make sure the mic and speaker stay in the same spot.


Difference between slow, medium, fast AGC In digital receivers, at the RF level (not audio), I often read about "slow AGC" and "fast AGC", and sometimes even "medium AGC" (but I think the latter is for audio). From my understanding, the slow AGC reacts to changes in noise level and variations in signal level over a relatively long period of time (fading for example). And from what I understand, the fast AGC is used more on symbol level, trying to keep them at a relatively stable power level, but not too fast so that amplitude-dependent symbols would not be affected (like M-QAM or M-PSK with M > 4).

Is there more to it? Also, would the "medium AGC" have any use in the RF chain?

I don't believe there is much more to it than that, especially noting that you are aware not to track so fast that you track out amplitude information for modulated signals where the amplitude is used as part of the modulation (so more so for M-QAM than M-PSK). You can also design AGC's with different rise/fall times so can have "fast-attack" and "slow-decay" for example.

A reason you may want to have a dynamic AGC is distinguishing operation between acquisition and tracking. During acquisition it would be likely to use a much faster AGC at the expense of noise performance for the benefit of acquiring the signal faster. Then once acquired the loops (AGC, carrier tracking, symbol tracking) can run slower at the benefit of better noise performance but fast enough to track all expected dynamics.


Can this be considered wide sense stationary? I was discussing this problem with one of my classmates.
The picture shows a recording of the heart rate during before and after sleep.

Can the whole process be considered wide sense stationary? (I say yes because the mean is approximately constant, but how can I estimate the autocorrelation?)
Can the whole process be considered ergodic in autocorrelation and mean? (I'm not sure at all about this)

for 1, a random process is WSS if its autocovariance and mean ensemble average don't vary as a function of specific time instance, just lag.  As mentioned before, it's hard to conclude this from simply one realization of the process about the ensemble average.

However for 2, you can at least estimate what the autocovariance and mean are if the process is ergodic (which, unfortunately, you can't conclude with just a single realization).  If the process is indeed ergodic, the ensemble average is equal to the time average in the limit as your interval of average approaches infinity, which is kind of nice with regards to computability.


Bandpass general equation to difference equation I know this is a very basic question and I am coming out from a quarter of DSP. I want to create a function in Java which can taken in two parameters, either (double centerFrequency, double bandWidth) or (double startFrequency, double cutoffFrequency).

Based on the information I've gathered, I want a non-ideal IIR filter because this is for a project I am doing with an accelerometer that needs to be updated real-time. Hopefully I am not going too much into programming, but the general form of my filter function would be like this:

double BPFilter(double x, double centerFrequency, double bandWidth) {
    double filteredY;
    static double[] y = {0, 0, 0};
    //y[0] represents y[n]
    //y[1] represents y[n-1]
    //y[2] represents y[n-2]

    y[0] = A*y[2] + B*y[1] + C*x;
    y[2] = y[1];
    y[1] = y[0];
    return y[0];
}

What would be the best way to calculate A, B, and C? Currently I know how to do continuous transfer functions, $H(s)$, so if that is translatable to A, B, and C I could do it that way. The sampling rate is at $100\textrm{ Hz}$.

TL;DR: Discretize a continuous transfer function by hand (no MATLAB)

The transfer function of a continuous-time second-order band-pass filter is given by

$$H(s)=\frac{\frac{\omega_0}{Q}s}{s^2+\frac{\omega_0}{Q}s+\omega_0^2}\tag{1}$$

where $\omega_0$ is the center frequency in radians per second, and $Q$ is the quality factor. For $Q\gg 1$, the term $\omega_0/Q$ closely approximates the $3$ dB bandwidth $W$ (in radians per second). Note that the condition $Q\gg 1$ is equivalent to $W\ll\omega_0$.

If you want to transform the transfer function $(1)$ to the discrete-time domain, there are several possibilities. The two most commonly used ones are the bilinear transform and the impulse invariance method. Both of them are covered extensively in the literature and on the web (including this site). Note that the bilinear transform introduces frequency warping, so you have to pre-warp the parameters of the continuous-time filter in such a way that after the transformation the discrete-time filter has the desired properties. A very good explanation of IIR filter design using the bilinear transform can be found in chapter 11 of the (freely available) book Introduction to Signal Processing by S.J. Orfanidis.

You can find specific formulas for second-order discrete-time filters derived from the bilinear transform in the Audio-EQ-Cookbook. Also take a look at this question and its answers.


What happens when the poles of this z-transform function are outside the ROC for a signal? I am given a z-transform function for a signal $h[n]$. It's $H(z)=\frac{2z^2-0.75z}{(z-0.25)(z-0.5)}$. I am supposed to find $h[n]$ and check the stability of the system for these cases:
a)ROC: $|z|>0.5$
a)ROC: $|z|<0.25$
a)ROC: $0.25<|z|<0.5$

I found $h[n]$ for the general case to be:
$H(z)=0.5z^{-1}\frac{z}{z-0.5}+0.25z^{-1}\frac{z}{z-0.25}$
$h[n]=(0.5)^nu[n-1] + 0.25^nu[n-1]$
The poles of a function are $0.5$ and $0.25$
The first case is when ROC: $|z|>0.5$. That would mean that i have a circle with radius $0.5$ on the graph on which the causal signals are the ones outside the circle because those values of $z$ make the signal (or a function) converge and anti-causal ones inside the circle. Yet the answer in the book is:
a)$h[n]=(0.5)^nu[n] + 0.25^nu[n]$ and the system is stable. Shouldn't it diverge for the $0.25$ part?
Also when b) ROC: $|z|<0.25$ the answer in the book is:
b) $h[n]=-(0.5)^nu[-n-1] + -0.25^nu[-n-1]$ and the system is unstable.
I could understand this as now the inside of a $0.25$ radius circle is causal and since pole 0.25 is on the circle and 0.5 is outside the circle the $u[n-1]$ is transformed to $u[-n-1]$ to make it causal, but why does the factor $0.25$ and $0.5$ change signes? I will also write down the solution for the c) part:
c) $h[n]=-(0.5)^nu[-n-1] + 0.25^nu[n]$
Can anyone explain these solutions as i'm confused about the whole ROC involvement in these solutions particular.

Each ROC has it’s own impulse response  .  There isn’t a general impulse response or a single ROC in most cases.

Your problem can have 3 ROC with 3 impulse responses and 3 stability conditions

One ROC is the causal ROC, Another is the anti-causal ( running the filter in reverse time) and the last, which isn’t common has a name that requires me to look it up, but this isn’t my homework problem so i’m going to stay put.

I know that in the first edition of Oppenheim and Schaefer, they explicitly say you can’t write down an impulse response from the Z transform without stating the ROC


How is frequency related to data rate? How does (if it does) frequency affect/relates to Data Rate.

Why for example using the same modulation, LoRa (in the physical layer),
you get somewhere around 250kbps at 2.4Ghz and only 20kbps at 868Hz ?

There's simply no inherent relation between carrier frequency and data rate.

I assume a classical comms engineering education here, but since that most often introduces the concept of equivalent baseband early, a good explanation would be:

  This becomes obvious when you think about the complex baseband representation of a passband being the same, regardless of the passband's center frequency, and that the complex baseband carries identical information to the passband signal.

  Thus, carrier frequency cannot have direct influence on the information transport capabilities of any system, given that you mix the same signal to these carriers.

Indirectly, there might be influences: Bandwidth is proportional to capacity; it's impossible to get 1 GHz bandwidth centered around 433 MHz, but cheap around 60 GHz.

Noise and interference are the other limit to capacity; these are different for different bands and technologies.


Finding the deterministic autocorrelation function (ACF) from its power spectrum The power spectrum of a stationary discrete-time random signal is

  $$\Phi_{xx}(e^{j\omega})=\begin{cases} 1 & |\omega|<\pi/2 \\ 0 & \pi/2
  <|\omega| \le\pi \end{cases} $$

  (a) Determine the autocorrelation function $\phi_{xx}[m]$

I know that the power spectrum is the DTFT of the ACF at 0, but how can I find the ACF from the power spectrum, given that we need to find the ACF for all possible n values?

The power spectral density is the fourier transformed ACF. So what you have to to is to do the inverse fourier transform of the spectrum. And since the spectrum is a rectangle in frequency, it will be a $\frac{\sin(x)}{x}$ in the time domain.

Additionally:

  I know that the power spectrum is the DTFT of the ACF at 0

The signals power is given by $\phi_{xx}(\tau = 0)$:
$$
\phi_{example}(\tau) = \int_{-\infty}^{\infty} x_{1}(t)* x_{2}(t-\tau)\mathrm{d}x
$$
If you autocorrelate (correlate a signal with itself), your subscripts vanish ($ x_1, x_2$ become $x$) and if you calculate the ACF at $\tau = 0$ you will get the squared power of the signal:
$$
\phi_{example}(\tau = 0) = \int_{-\infty}^{\infty} x(t)* x(t)\mathrm{d}x = \int_{-\infty}^{\infty} x(t)^2\mathrm{d}x
$$

If your signal does not have a finite amount of energy, you will have to use the following definition of the ACF:
$$
\phi_{example}(\tau) = \lim_{T \to \infty} \frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}} x(t)* x(t-\tau)\mathrm{d}x
$$


Image Processing: Flatten a 3D Ball There are many descriptions of how to turn a 2D image into a 3D one, however I want to do the opposite, in particular to a ball. As an example, consider the following ball:

If only the (relatively) slow change in intensity due to the light could be removed, the ball would look 2D.

My first approach was to try get a low frequency version of the original, and divide through by that image. A simple way to do that is to blur the image. Here is a sample result using that approach:

This does the job of flattening the light, but there are some problems:

There is a light 'halo' around sharp edges.
Large flat patches (e.g. the black patches) lose their uniformity in the middle.
The edges of the ball are affected by the background during the blurring (there's a dark ring around the ball, particularly noticeable at the bottom right corner).

I've tried mucking about with filters (e.g. bilinear) but haven't really found anything that works. I've also tried to use wavelets (a frequency based approach makes sense to me), but those approaches tended to distort the fine detail, which is very important.

Any ideas on how to do this are greatly appreciated!

My best attempt at this so far is kind of inspired by the difference of Gaussians. I'm putting it as an answer, even though I'm not totally happy with it.

Basically I made a blurred version of the image using two different Gaussian kernels. One captures the fine detail (small kernel) and the other captures the broad detail (large kernel). These two were then merged. This new blurred image is then subtracted from the original.

UPDATE: Made the process iterative, added contrast equalisation and changed blur divide to blur subtract.

output_image = get_input_image()

for i in xrange(4):
    # Magic numbers
    sigma_broad = 301
    sigma_fine  = 5
    merge_ratio = 0.5

    # Get blur
    broad_detail = cv2.GaussianBlur(output_image, (sigma_broad, sigma_broad), 0)
    fine_detail  = cv2.GaussianBlur(output_image, (sigma_fine, sigma_fine), 0)
    blurred = merge_ratio*broad_detail + (1-merge_ratio)*fine_detail

    # Subtract blur from original (and correct for average drop in intensity)
    output_image = (output_image + np.mean(blurred)) - blurred

# Contrast equalisation
clahe = cv2.createCLAHE(clipLimit=1.0, tileGridSize=(4,4))
output_image = clahe.apply(output_image.astype(np.uint8))

Result:

The reasons that I'm not totally satisfied with this result:

The information I want is altered. The objects in the ball (e.g. the dark patches) are disproportionally reduced.
There are some artefacts.
There are lots of magic numbers and parameters to fiddle with. Ideally the algorithm should be robust for a variety of balls.

ANOTHER UPDATE:
I had a go at a Laplacian Pyramid filter, which I adapted from here. I got decent results, posting mainly for interest. This strategy allows for targeting frequency bands, but it's very fiddly.


How to calculate critical damping of a system with two springs and a damper (or two springs and two dampers)? Background
For a simple system where you have a mass attached to a spring and damper in parallel:

We can calculate the critical damping from the equation of motion:
$mx_{tt} + cx_t + kx = 0$
$ms^2 + cs + k = 0$
$s= \frac{-c ± \sqrt{c^2-4mk}}{2m}$
There are then three conditions:

$c^2 <4mk$ (under damping)
$c^2 >4mk$ (over damping)
$c^2 =4mk$ (critical damping)

The damping ratio is then expressed by $\frac{c}{\sqrt{4mk}}$.
Question
I'm wondering if or how this can be extended to more complex systems.
Let's say you add a second spring to the parallel system:

This is a viscoelastic model where the constitutive relationship is in terms of stress (σ) and strain (Ɛ):
$σ = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ}$
Is it possible to calculate the critical damping or damping ratio in the same manner? If so how would it work?
What about with two springs and two dampers like this?

$σ = (η_1+η_2)\dot{ϵ} + \frac{η_1η_2(E_1+E_2)}{E_1E_2} - (\frac{η_1}{E_1} + \frac{η_2}{E_2})\dot{σ} - \frac{η_1η_2}{E_1E_2}\ddot{σ}$
Is it possible to do the same and if so how?
My Guess
My guess is strain is equivalent to $x$ and stress is equal to force, so we can reword the 3 element system:
$0 = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ} - σ$
$0 = E_1x + \frac{η(E_1+E_2)}{E_2}\dot{x} - m\frac{η}{E_2}\dddot{x} - m\ddot{x}$
$0 = E_1 + \frac{η(E_1+E_2)}{E_2}s - ms^2 - m\frac{η}{E_2}s^3$
Which according to Wolfram gives the following roots:
$s = \frac{-E_1E_2}{η(E_1+E_2)}$
$s = ±\sqrt{\frac{E_1}{m}}$
But that is not useful for solving any critical damping scenario that I can figure out.
Is there some other approach or what did I screw up?
Any help or guidance is appreciated. Thanks.

Re-writing the strain-stress equation
$$
0 = E_1Ɛ + \frac{η(E_1+E_2)}{E_2}\dot{Ɛ} - \frac{η}{E_2}\dot{σ} - σ
$$
for displacement/restoring force variables
$$
0 = E_1x + \frac{η(E_1+E_2)}{E_2}\dot{x} + m\frac{η}{E_2}\dddot{x} + m\ddot{x}
$$
and applying the Laplace transform, we arrive at the equation
$$
0 = E_1 + \frac{η(E_1+E_2)}{E_2}s + ms^2 + m\frac{η}{E_2}s^3 \tag1
$$
The formula for root $s = \frac{-E_1E_2}{η(E_1+E_2)} $ is only valid when $m=0$, and for root $s = ±\sqrt{-\frac{E_1}{m}}$ is only valid when $η=0$. In a general case of nonzero $m, η$, the equation (1) has three roots. The general formula for cubic equation roots is complicated. To write it down, the Wikipedia article introduces two first-level "determinants" $Δ_0 = b^2-3ac$ and $Δ_1 = 2b^3-9abc+27a^2d$, where $a, b, c, d$ are the cubic equation coefficients: $ax^3+bx^2+cx+d=0$, Then it defines an "upper-level" determinant $C$, built with the "first-level" determinants,
$$
C=\sqrt[3]{{Δ_1 ± \sqrt{Δ_1^2 - 4Δ_0^3}}\over{2}}
$$
the sign ± is selected to be plus or minus in order to guarantee that $C$ is never zero, and the three roots are
$$
x_k = -{1\over{3a}}\left(b + ξ^kC+{{Δ_0}\over{ξ^kC}}\right), \, k \in \{0,1,2\}, \, ξ={{-1+\sqrt{-3}}\over2} \tag2
$$
For a parallel-serial-connected spring-damper-mass system of your question, $a = m\frac{η}{E_2}, b = m, c = \frac{η(E_1+E_2)}{E_2}, d = E_1$.
The equation (1) reduces to a second-order ODE describing the simple damped harmonic oscillator:

for a trivial case of $η=0$, where there is no damping at all, and
the spring E2 does not participate in the motion equations;
the other non-dissipative case is $η\to{\infty}$, the springs act
in parallel with the total stiffness of $E_1+E_2$;
with $E_2\to\infty$, the bottom spring is "turned off", the system
has one spring and one damper in parallel, this system you consider in the
beginning of your question;
with $E_1=0$, the upper spring is "turned off", the system has one
spring and one damper in series, the system also has overdamped,
underdamped and critically damped solutions that can be found with the
formulas from your link.

For non-zero, finite values of $η, E_0, E_1, m$, iff the system has degenerated roots for certain combinations of $η, E_0, E_1, m$ values, then one of its linearly independent solutions has the form $t·exp(x_kt)$, the root $x_k$ has a negative real part. This solution may be called the critically damped solution. To find these combinations of $η, E_0, E_1, m$ values, one must search for degenerate solutions of Eq. 1 with multiple roots. The formulas are certainly not that simple as the formulas for a damped harmonic oscillator, but the ODE for the mechanical system of eq. 1, by virtue of an electrical-mechanical analogy, is identical to the ODE of a circuit with a Cauer topology of the generic three-pole unbalanced T-section LC low pass filter, the latter being often discussed on electronics forums.


Detection of sharp peaks in an audio signal I have an audio signal and I'm interested to detect the very sharp peaks.

If I zoom, here is something that should NOT be detected (this is voice actually):

NO GOOD

Here is something that should be detected:

GOOD

Ideally, a CNN with some training would be the best solution but I'm too limited in terms of computation to run it.

So, I'm looking for an easy solution with a simple signal processing algorithm (filters and others things).

I have already developed a simple algorithm based on the maximum of a 10ms window and if I see that the maximum is low, then suddenly very high, and then low again, I have a detection.

But it's not outstanding. Any better idea? Is there any filter that can find me "sharp" peaks?

"Sharp peaks" in the time domain translates to high frequency content. Your simple algorithm of min and max values would be improved by first passing your data through a high pass filter - otherwise slow (low frequency) rumbles could degrade your detection methods.

To take your solution to the next level, you should use an "onset detection" algorithm. These are well studied. A recommended paper to start is:

A Tutorial on Onset Detection in Music Signals - Juan Pablo Bello et al, September 2005.


Next steps for analyzing this PWM algorithm and choosing a low pass filtering technique? The Broadcom BCM2711 GPIO controller used by the Raspberry Pi has a hardware PWM implementation option using an algorithm that provides an average value defined by a rational fraction $N/M$ defined as follows in section 8.3. PWM Implementation starting on page 156:

  Set context = 0
  context = context + N
  if (context >= M)
  context = context - M
  send 1
  else
  send 0
  Repeat from step 2

  where context is a register which stores the result of the additions/subtractions.

I've implemented this in a Python script to try to start to get an idea how I'll need to filter it in order to get a low-noise, stable voltage. Here with $M=200$ I've incremented through $0 \le N \le M$ for over 1000 time points each, then applied a gaussian smoothing with a $\sigma$ of 10 time points just to get a quick look.

I took the Fourier transform of the middle 1024 points, ignoring a $3 \sigma$ buffer on each side.

I was surprised to see that some values of $N$ produce signals that are orders of magnitude more noisy than other nearby values of $N$.

Questions:

This kind of analysis is new territory for me, if I've made any major errors in the analysis I would appreciate knowing about it.
If I need near-DC voltages I can just run the BCM2711 at a moderately high frequency and filter the heck out of the output with a simple RC filter.  But if I wanted to work closer to the frequency that the algorithm was running, what additional analysis should I do on this algorithm before deciding on a more aggressive low-pass filter and simulating it's performance?

Here's the Python script for the plots and analysis:

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage
from numpy.fft import rfft, rfftfreq

sigma = 10
ignore = 3 * sigma
M = 200
big = []
for N in range(M+1):
    context = 0
    bob = []
    for i in range(1024 + 2 * ignore):
        context += N
        if (context >= M):
            context -= M
            bob.append(1)
        else:
            bob.append(0)
    big.append(bob)

big = np.array(big, dtype=float)
fbig = ndimage.gaussian_filter1d(big, sigma=sigma, mode='mirror', order=0, axis=-1)
fbig_trunc = fbig[:, ignore:-ignore]
print('fbig_trunc.shape: ', fbig_trunc.shape)
ft = rfft(fbig_trunc, axis=-1)
power = np.abs(ft)**2
power /= power.max()
dB_rel = 10 * np.log10(power)
print('power.shape: ', power.shape)
print('power.max(): ', power.max())
print('dB_rel.shape: ', dB_rel.shape)
print('dB_rel.max(): ', dB_rel.max())

if True:
    plt.figure()
    for thing in fbig:
        plt.plot(thing, linewidth=0.8)
    plt.ylim(-0.02, 1.02)
    plt.title('sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
    plt.xlabel('time steps', fontsize=16)
    plt.show()

if True:
    plt.figure()
    maxi = dB_rel[:, 1:201].max()
    plt.imshow(dB_rel[:, :201], vmax=maxi, vmin=-100)
    plt.colorbar()
    plt.gca().set_aspect(2.0)
    plt.title('FT relative power (dB), sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
    plt.xlabel('frequency', fontsize=16)
    plt.ylabel('N', fontsize=16)
    plt.show()

if True:
    plt.figure()
    std = fbig_trunc.std(axis=-1)
    plt.plot(10*np.log10(std**2))
    plt.title('stdev relative power, sigma={:2d}, M={:3d}'.format(sigma, M), fontsize=16)
    plt.xlabel('N', fontsize=16)
    plt.ylabel('stdev (dB rel)', fontsize=16)
    plt.show()

The pulse generation algorithm sounds like it is more like Pulse Density Modulation (PDM), not Pulse Width Modulation (PWM), because PWM would pack the ones into the beginning of counting, and leave the zeroes in the end.

Indeed some values of N and M cause the output pattern to have a certain repetitive bit stream, for example the 100/200 means the output is 50% 0 and 50% 1, so the output is a bitstream of 1,0,1,0... which is easy to filter. If you set N and M to be exactly ratio of 1/3, the output pattern will also be repetitive, 0,0,1,0,0,1... So, in light of this this also explains the spectrogram of N=50 M=200, which is 25%, 0,0,0,1,0,0,0,1.

These patterns that have a short repetition rate show up as pure tone peaks in the spectrogram. Choosing the values so that there is no obvious short repeating patterns basically dithers the bitstream and widens the spectrum to have multiple lower peaks.

Depending on what is the usage of the output will determine how to filter it. For many cases a simple RC filter is enough, as long as the PDM bitstream output rate is high enough for the RC filter cutoff frequency.


Calculating fourier transform I have just recently started doing fourier transforms and I'm a little confused.

Can someone walk me through in detail how to calculate the Fourier transform of

I'm not looking for answer, just an explanation of how to go about doing it.

Sinusoidals are one of the easier signals to analyze given how you can manipulate the expression by algebraic tricks.

for this function particularly, I would convert sine to its euler form first:

Then use the Fourier Transform definition:


Sample Rate: How to determine sample rate I have this ADC and I have this 169p micro-controller. The system oscillator clock frequency (FOSC) is 8 MHz.

According to the ADC documentation the conversion process takes 16 cycles. So if I have FOSC of 8MHz would I get a sampling rate of 8000000/16 = 500,000 samples per second? The ADC has 4 channels so can I divide the sample rate by 4? 500000/4 = 125,000 samples each 16 bits (2 leading zeros, 2 for channel, 12 for conversion).

This question isn't particularly on-topic, but the answer is pretty simple. Look at the second paragraph of the General Description on the first page of the ADC datasheet:

  The conversion process and data acquisition are controlled using
  CS and the serial clock signal, allowing the device to easily interface
  with microprocessors or DSPs. The input signal is sampled on the falling edge of CS and conversion is initiated at this point. There are no pipeline delays associated with the part.

CS in this case refers to the chip select pin on the ADC's SPI-compatible 4-wire serial bus. That signal is asserted each time the host processor initiates a transaction to the converter. The mechanics of an SPI read go like this:

The SPI master (i.e. the host processor) asserts the CS signal.
If the SPI master has any data to write to the device, it presents it in serial fashion on the MOSI (master out, slave in) signal. From the ADC's perspective, this would be the DIN pin.
The SPI master begins to apply a clock signal to the ADC on the SCLK pin. This is used to simultaneously shift data in and out of the ADC serially.
For this ADC device, the transaction continues for 16 bits. After that, the master has received 16 serial bits from the ADC, which contain the conversion result.
The master deasserts the CS pin to end the transaction.

In the above timeline, the speed of the conversion is therefore limited by two factors:

The SPI clock rate (as it determines how long it takes to extract 16 bits).
The maximum conversion rate that the ADC supports.

#2 in the above list is very clearly stated to be 1 MHz. #1 is dictated by the capabilities of your host processor, and can be found in the datasheet. Specifically, look at the bottom of page 166:

  Bit 0 – SPI2X: Double SPI Speed Bit

  When this bit is written logic one the SPI speed (SCK Frequency) will be doubled when the SPI is in Master mode (see Table 18-5
  ). This means that the minimum SCK period will be two CPU clock periods.

It is stated elsewhere in the datasheet that a "CPU clock period" as referenced here is equivalent to $\frac{1}{f_{OSC}}$. Therefore, if your processor has a CPU clock frequency of 8 MHz, its maximum SPI clock rate is:

$$
f_{max,SPI} = \frac{f_{OSC}}{2} = 4\text{ MHz}
$$

Remember though that this is a serial bus with 16 bits per conversion transaction. Therefore, the maximum ADC sample rate is:

$$
f_{max,ADC} = \frac{f_{max,SPI}}{16} = 0.25\text{ MHz}
$$

Thus, the processor is the limiting factor here; it can only support a maximum of a 250 kHz sample rate. Furthermore, this is the maximum throughput that is available to the entire ADC device, so if you're using more than one channel on the ADC, this 250 kHz would have to be split between the channels you're using.


Use Scale Space Representation to Filter Single Image Currently I hope to use scale space representation to filter one image. Features in one image can be filtered using an Gaussian smooth filter with one optimal sigma. It means different features in one image can be expressed best in different scale under scale space representation.

For example, I have one image with one tree in it. In the scale space representation, three sigma values are used and they are represented as sigma0, sigma1 and sigma2. The ground is best expressed in the smoothed image with sigma0 because it contains textures mainly. The branches are best expressed in the smoother image with sigma1 and the trunk is with the smoother image with sigma2. If I hope to filter the image, I hope that the filtered pixels for the group is from the smoothed image with sigma0.
The filtered pixels for the branches are from the smoothed image with sigma1. The filtered pixels for the trunk are from the smoothed image with sigma2.

It requires that I need to determine in which smoothed image one pixel is expressed best. Is this idea plausible?
I am trying to use differece-of-Gaussian of two successive smoothed images to perform the above task. Is there any other way to combine the three smoothed image?

I use MATLAB to implement the idea. The values of the three sigmas is 1.0, 2.0 and 3.0. The corresponding size of Gaussian kernel is 3, 5 and 7. I use the function fspecial to generate the kernel. Are the parameter reasonable? Please share your experience with the scale space representation to help me. You can provide some links to useful papers.

I used factor 5-6 going from the Standard Deviation (Std) of the Gaussian to the Kernel:
radius = ceil(6 * STD);

Though it means more computation power is required.


What is localizability in Computer Vision? Please consider the following excerpt (from Rodrigo R, Zouqi M, Chen Z, Samarabandu J.: Robust and efficient feature tracking for indoor navigation):

  Robust feature tracking is a requirement for many computer vision tasks such as indoor robot navigation. However, indoor scenes are characterized by poorly localizable features.

I am trying to understand what localizability is, and would appreciate help.  (PS: I know there is a dsp tag called "local-features", but the tag currently has no definition.)

In the robot navigation problem, the localization problem refers to the real time estimation of its position and orientation under various backgrounds. This is usually achieved by some natural landmark selection (laser points, camera views, etc.), and the features in the image (corners, tiny lines with different orientations, etc.). So the localizability reflects how easy, robust and consistent the robot's position and orientation can be estimated with the choice of feature points.

Hugues Sert et al's Mobile robot algebraic localizability includes more strict definition with details, and its relation to observability.


linear phase notch filter matlab I need to filter 9 Hz from a signal sampled by 256 Hz using a linear phase filter.

If someone can bring an explanation or a code example I would greatly appreciate it.

I have tested his code

function [Hd1,Hd2]=NotchFIR(FS,F0,FG,FGAINdB) % Design
    if (F0 > FS/2)
        warnstr =['WARNING! The Sampling frequency is less than twice the notch frequency. Response will be inaccurate.']
    end
    if (FG == F0)
        warnstr =['WARNING! User should not set the gain of the notch frequency. Response will be inaccurate.']
    end
    FGAIN = 10^(FGAINdB/20);
    w0 = 2*pi*F0/FS;
    wg = 2*pi*FG/FS;
    GH = abs(1 - 2*cos(w0)*exp(-j*wg) + exp(-j*2*wg));
    G = FGAIN/GH;
    h = [1 -2*cos(w0) 1];
    freqz(G*h,100,[], FS);
    Hd1=G*h;
    Hd2=100; end

     t=0:1/256:1;
     y=sin(2*pi*9*t);    [Hd1,Hd2]=NotchFIR(256,8.9,9.1,1)
     x=filter(Hd1,Hd2,y);
     figure
     plot(t,x);

But the frequency response of the digital filter shows that the selection of the notch frequency wasn't sharp and many frequencies were attenuated.

@omar: I don't know how you obtained your expression for the 'h' coefficients but it's wildly incorrect.  Study the MATLAB documentation regarding the 'firpm()' and 'fir1()' filter design commands.  Experiment, experiment, and experiment again until you understand those commands. Then when you try to use one of those commands to design your desired filter you'll find that your filter needs a truly HUGE number of taps.


Discrete Cosine Transform (DCT) I need your help. My lecturer told me to do a study to finish the task that she gave us yesterday but the thing is there is no specific explaination on how to solve the question. Through the task, she gave us only this information. I went through youtubes in case there'll be explaination on those but so far i didnt find anything. The task is due in the next 4 days.
I don't need answer, i just need to know how it works and so on and i'll do the calculation.
[UPDATE:]
Here is my calculation. So the answer for H(k) will be in vector right? If it is, i proceed with the second question but only 1/2 were multiply to the original formula for F(k) right? What can i conclude for the comparisons of both signal? Please correct me if i'm wrong.

I can't recall which type of DCT your formula refers to, but assuming the formula is right, then the summation in the first equation computes $N$ coefficients of the DCT sequence $F[k]$ of the $N$-point input sequence $x[n]$. A similar interpretation applies to the inverse formula.
The following Matlab/Octave code implements the first formula:
N = 6;              % length of input sequence x[n], and output sequence F[k]
x = [1,0,2,5,3,7];  % the input

F = zeros(1,N);     % the output: DCT coeffs

for k = 0 :N-1
    csum = 0;
    for n = 0: N-1
        csum = csum + 2*x(n+1)*cos((pi/(2*N))*k*(2*n+1));
    end
    F(k+1) = csum;
end


Downsampling vs. ADC with lower Sample-Rate I would like to know in general where the advantages and disadvantages are:

Downsampling of a high sampled continuous analog signal in FPGA or µC and
Direct use of a lower sampling rate, i.e. using an analog to digital converter with a lower sampling rate ?
Are there any differences or does it have the same effect ?
Because almost every DSP-System i see uses a way Oversampling followed by Decimation (Downsampling + AAF-Filter) instead of adequate Sampling Rate without Decimation.

Greetings

With downsampling you have complete control over the process and it comes down to what compromise of processing complexity, delay, aliasing and loss of passband you can accept.
With a lower rate A/D you are pretty much at the mercy of someone elses spectral trade-offs and in addition you get the quantization/noise of one analog pass.
-k


What are some easy to detect reference mark shapes and corresponding detection methods? Question Details:

There are a variety of image processing algorithms useful for detecting reference marks (contour tracing, FAST, MSER, Haar Feature-based Cascade Classifier, or even just Harris Edge detector, etc).

There are a variety of simple geometric shapes one might pick for use as reference marks (circle, square, X, etc).

I am wondering if someone has some experience with what combinations work best? Some additional criteria I have:

using an algorithm implemented in openCV would be great
Background invariance: detection algorithm should be robust against
random background features(folds in clothing, random sharp corners, random
colors present in background, speckle noise)
Distance invariance: detection algorithm should be stable with
respect to scale. (i.e. is should work whether the reference mark
takes up a 100x100 pixel square in the image or 30x30 pixel square)
Angle invariance: detection algorithm should not fail if the
reference mark is viewed at an angle
I am using a 640x480 webcam as a source of the images, so image processing at better than 5-10 fps would be nice

Additional Info:

My goal is to print out a fixed number black shapes on white paper. Place them in arbitrary positions in the scene that I am imaging with a webcam, and try to detect the center coordinates of each of the black shapes using a c++ program w/ openCV.

I am wondering if there is any consensus about what type of shape is the easiest to detect and how to do it.

Approaches that I have tried:

Detecting concentric circles using contour tracing: 1) Find the contours that are n times nested 2) compute the center point.  The problem is that it is hard to resolve the rings at far distances (>5 feet) and angles.

Detecting cross shapes using contour tracing: 1) find contours, 2) find convex hull of each contour, 3) find number of convexity defects 4) compute the center if convexity defects == 4. Alternatively I have tried finding the depth of each convexity defect and comparing it to the width of the bounding rectangle around the contour to determine if it is a "real" convexity defect or just due to noise. This method works but sometimes certain background features, like folds in clothing, can make extra objects detected.

These marks (in the electronics manufacturing industry) are usually called fiducials or fiducial markers.

This paper from Accutron suggests:

Though I have no reference (or even measure) of the "optimality" of the disk / circle.