What is the practical meaning of the variance, covariance, mean value? What is the practical meaning of variance, covariance and mean for a signal?

Mean of a signal can be practically visualized as the dc average value present in the signal (for a complete sinusoidal period), for e.g

Variance of a signal is the difference between the normalized squared sum of instantaneous values with the mean value. In other words it provides you with the deviation of the signal from its mean value. It gives you the spread of your signal's data set.

Co-variance of a signal consists of unnormalized relation with the signal correlated with itself and also with other signals. It provides you the degree of similarity between the signals of interest.

https://en.wikipedia.org/wiki/Variance

http://www.dspguide.com/ch2/2.htm

https://en.wikipedia.org/wiki/Covariance

EDIT : FOR COMPLETE INSIGHT PLEASE REFER ALSO TO @Dan Boschen's ANSWER.


Sample & Hold: Estimate jitter delay of an ADC I've got a question regarding the effects of sample and hold.
My input signal is $x$, and my output signal $x_{SH}$. The error between the signals is $e_{SH} = x_{SH}-x$.

The sampling frequency is set to $f_s = 1000kHz$ and the jitter delay to $t_{jitter}=4\cdot 10^{-4}s$.

The influence of the jitter is implemented in Matlab as the following:

t = (0:N-1)/f_S;
x = sin(2*pi*f*t);
t_SH = t+Tjitter*(2*rand(size(t))-1);
x_SH = sin(2*pi*f*t_SH);
e_SH = x_SH-x;
rxx = xcorr(x,x,'biased');
ree = xcorr(e_SH,e_SH,'biased');

Now I'd like to estimate $t_{jitter}$ using the formula:

$SNR = 20 \log_{10}\left[\frac{1/\sqrt{2}}{\Delta v_{rms}}\right]=20 \log_{10} \left[\frac{1}{2\pi f t_{jitter}}\right]\\$
$\Rightarrow t_{jitter} = \frac{1}{ 2\pi f \left(\sqrt{ \frac{P_{signal}}{P_{noise}} }\right) }=\frac{1}{ 2\pi f \left(\sqrt{ \frac{max(r_{xx})}{max(r_{ee})}} \right) }$

Now I only achieve the same result of $t_{jitter}$ if I add a factor of $\sqrt{3}$, which I can't explain why.

$\Rightarrow t_{jitter} = \frac{1}{ 2\pi f \left(\sqrt{ \frac{P_{signal}}{P_{noise}\cdot 3} }\right) }$

I'd be happy for some help! Thanks in advance!

I think the jitter to SNR formula is based on a gaussian jitter.

You use "rand" which yields uniformly distributed numbers. You should use randn() instead which yields numbers distributed according to a gaussian distribution.


Puzzled by CIC filters I have been trying to wrap my head around CIC filters.

Say that I want to do filtered decimation by a factor D.

Option#1: boxcar filter. Each output sample is formed by taking the sum of D input samples, and multiplying/shifting to get a proper average. D-1 adds per output sample.

Option#2: 1-stage CIC filter with an accumulator running at the input rate, and a comb filter running at the output rate. D adds per output sample.

So what is the fuzz about? Does CIC only make sense for higher orders? If order N CIC is like convolving the rect() function N-1 times, that would be a B-spline, right? I have seen authors praising B-splines for their functional benefits. Being able to do them efficiently might be handy.

Thanks for any input
Knut

Edit: pseudo matlab code showing the alternatives as I understand them.

N=100;
x = randn(N,1);
D = 5;

%% decimated moving average, ~100 adds
for k = 1:N/D
  y(k) = (1/D)*(x(k*D) + x(k*D-1) + ... + x(k*D-D+1));
end

 %% recursive moving average, ~200 adds
y(1) = x(1);

for k = 2:N
  y(k) = x(k) + y(k-1);
end

for k = D:N
    z(k) = y(k) - y(k-D);
end

z = z(1:D:end);

%% CIC, ~120 adds
y(1) = x(1);

for k = 2:N
  y(k) = x(k) + y(k-1);
end

z = y(1:D:end);

for k = 2:N/D
  w(k) = z(k) - z(k-1);
end

It seems you mixed up CIC filters and recursive moving-average

A recursive moving-average uses a comb-filter and an accumulator.

$$ H(z) = \frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}} $$

The numerator corresponds to a comb filter (i.e a delay line of N elements and 1 adder).
The denominator corresponds to an accumulator. The 1/N factor is to get a moving-average, otherwise you'd have a moving-accumulator. You only need 2 additions, a delay line, an accumulator and a divider. With a non-recursive moving-average filter, you'd need N additions and you'd need a delay line anyway. A recursive moving-average is much
less computionally expansive than a non-recursive moving-average filter.

Now the CIC filter combines the recursive moving-average filter with decimation and interpolation.

For a CIC decimation filter, you insert the rate transition block between the accumulator and the comb filter. Say you have a 16-sample moving average-filter (N=16) and you want to decimate by 8 (R = 8). You'd need a delay line of 16 elements. However, by inserting the decimation between the accumulator and the comb filter, you can reduce the number of elements in the delay line. You'd only need N/R elements, i.e. 2. You'd still need to subtract 2 numbers, however you would do it at sampling rate of fs/8.

I suggest you consult this site

https://www.researchgate.net/publication/228648728_Understanding_cascaded_integrator-comb_filters

Update: see this nice write-up by Richard Lyons-
https://www.dsprelated.com/showarticle/1337.php


Why is the digital frequency highest at $\pi$ and what does it mean in terms of analog frequency? I can't seem to comprehend why the digital frequency is highest at $\pi$. I know the relationship between the digital and analog frequency is given by:
$$\omega = 2\pi\frac{F}{F_s}=2\pi f$$
Where $F$ is the analog frequency and $F_s$ is the sampling frequency.
Why does the digital frequency increase from $0$ to $\pi$ and then decrease from $\pi$ to $2\pi$?

Check this figure.

Consider the signal $x[n]=A$ which is constant but can be written as $x[n]=A\cos(0n)$, so its frequency is zero. There is no oscillation.

Let's increase the frequency. For the signal $x[n]=A\cos(\pi n/4)$, it oscillates at frequency $\omega_0=\pi/4$. The oscillation frequency increases.

For the frequency $\omega_0=\pi$, the signal is $x[n]=A\cos(\pi n)=A(-1)^n$, thus the signal oscillates between $A$ and $-A$ for each pair of consequtive samples. You can see that the signal cannot oscillate at a higher frequency, there is no faster oscillation than this.

For frequency $\omega_0=7\pi/4$, the signal can be written as $$x[n]=A\cos\Big(\frac{7\pi n}{4}\Big) = A\cos\Big(2\pi n-\frac{\pi n}{4}\Big)=A\cos\Big(\frac{\pi n}{4}\Big)$$ which oscillates at a lower frequency.

Finally, for $\omega_0=2\pi$, you have a constant signal again, no oscillation.

It is the nature of discrete time that introduces this effect.


Transform minimum phase FIR into linear phase FIR I've seen examples of transforming a linear phase FIR into a minimum phase FIR, but is there a simple process to transform a minimum phase FIR into a linear phase FIR?

I would like to end up with a single FIR of linear phase. If the magnitude ends up being squared during the process (or the number of taps double), then that would be OK.

It is generally impossible to transform a given minimum-phase FIR system into a linear phase FIR system with the same magnitude response. There is one special case for which this is possible, and that is if the zeros of the minimum phase system inside the unit circle haven even multiplicity. Because in that case you can, for each zero location, mirror half of the multiple zeros across the unit circle, which gives you a linear phase filter without changing the magnitude. But that's of course a quite artificial situation, so in general there's no way to obtain a linear phase filter from a minimum phase filter without changing the magnitude.

A way to transform a minimum-phase FIR system into a linear-phase FIR system with a squared magnitude is by simply doubling all zeros on the unit circle, and for each zero inside the unit circle adding a mirrored zero outside the unit circle. If $H(z)$ is a given $N^{th}$ order causal minimum-phase transfer function, the corresponding linear phase transfer function is given by

$$G(z)=H(z)H(1/z^*)z^{-N}\tag{1}$$

Note that the number of coefficients is (approximately) doubled and that

$$|G(e^{j\omega})|=|H(e^{j\omega})|^2\tag{2}$$


Mathematics of Harris corner point detection This is the mathematical expression for Harris corner detection:

But I have the following doubts:

What is the physical significance of $u$ and $v$? Many references say it is the magnitude by which the window $w$ shifted. So how much is the window shifted? One pixel or two pixels?
Is the summation over the pixel positions covered by the window?
Assuming simply $w(x,y) = 1$ , $I(x,y)$ is intensity of the single pixel at $(x,y)$ or the summation of the intensities inside the window with center at $(x,y)$?
According to wiki they say the image is 2D , denoted by I and then asks to consider an image patch over the area $(x,y)$, then uses the notation $I(x,y)$

I am finding it confusing to grasp the mathematical explanation. Anyone has an idea?

The meaning of that formula is really quite simple. Imagine you take two same-sized small areas of an image, the blue one and the red one:

The window function equals 0 outside the red rectangle (for simplicity, we can assume the window is simply constant within the red rectangle). So the window function selects which pixels you want to look at and assigns relative weights to each pixel. (Most common is the Gaussian window, because it's rotationally symmetric, efficient to calculate and emphasizes the pixels near the center of the window.) The blue rectangle is shifted by (u,v).

Next you calculate the sum of squared difference between the image parts marked red and blue, i.e. you subtract them pixel by pixel, square the difference and sum up the result (assuming, for simplicity that the window = 1 in the area we're looking at). This gives you one number for every possible (u,v) -> E(u,v).

Let's see what happens if we calculate that for different values of u/v:

First keep v=0:

This should be no surprise: The difference between the image parts is lowest when the offset (u,v) between them is 0. As you increase the distance between the two patches, the sum of squared differences also increases.

Keeping u=0:

The plot looks similar, but the sum of squared differences between the two image parts is a lot smaller when you shift the blue rectangle in the direction of the edge.

A full plot of E(u,v) looks like this:

The plot looks a bit like a "canyon": There's only a small difference if you shift the image in the direction of the canyon. That's because this image patch has a dominant (vertical) orientation.

We can do the same for a different image patch:

Here, the plot of E(u,v) looks different:

No matter which way you shift the patch, it always looks different.

So the shape of the function E(u,v) tells us something about the image patch

if E(u,v) is near 0 everywhere, there is no texture in the image patch you're looking at
if E(u,v) is "canyon-shaped", the patch has a dominant orientation (this could be an edge or a texture)
if E(u,v) is "cone-shaped", the patch has texture, but no dominant orientation. That's the kind of patch a corner-detector is looking for.

  Many references say it is the magnitude by which the window 'w' shifted...so how much is the window shifted?one pixel...two pixels?

Normally, you don't calculate E(u,v) at all. You're only interested in the shape of it in the neighborhood of (u,v)=(0,0). So you just want the Taylor expansion of E(u,v) near (0,0), which completely describes the "shape" of it.

  Is the summation over the pixel positions covered by the window?

Mathematically speaking, it's more elegant to let the summation range over all pixels. Practically speaking there's no point in summing pixels where the window is 0.


Which one is good Image processing language? For my academic project, I need to process image to get pixel color and also needs to change the pixel color.I need pattern recognition. It is a web application. All I am trying is to capture the image of debit card by webcam and read the details. Please suggest me programming language.

I'am using Opencv, you can program in it with C++ and Python, for example.


Which is the most suitable programming language for image processing in scientific researches? We are currently working on a platform (web application for a graduation final project we used mean stack) that able researchers to test default methods recognition (methods of Enhancement, Extraction, and matching ), add his own methods too (a script), compose a set of methods ( build a process of recognition) and visualize the outputs of each step. We want to know what a programming language is used the most for testing image processing in scientific articles.for the moment we are using Matlab for its simplicity.

Regards.

Your question as is has a very simple answer: "the best language is the one I use". Pros and Cons (P&C) would give some more value to it, as commented by Peter and Marcus. So, back to the topic, but concerning with P&C, here is my $\epsilon$...

Your application does not seem to be restrict only to simulation, but also to web interface; moreover, you want a roadmap to it, probably in a asynchronous way dealing with branches, and therefore with some version control.

Matlab can be a choice if your audience has all toolboxes you have for this work (web + image + signal + ... ). Does it?

Forget control versioning if you intend to use Simulink.

I would make thumbs up for Python as well. You have support of a huge community. The scientific packages are numerous. It is a complete language, from web applications to visual to scripting to ... And it is free!


Calculating SPL from pressure signal - Amplitude vs Power method I have a pressure signal from a Fluent FFowcs-Williams Hawkings acoustics analysis.

I converted this pressure signal into the frequency domain in order to get SPL values, using Matlab. I used the amplitude of the FFT:

A = fft(s10);        %Amplitude
B = abs(A);          %Magnitude of amplitude (complex numbers to real)
C = B/L;             %Normalize amplitude by dividing by number of samples
D = C(1:fix(L/2)+1); %Select one half of 2-sided spectrum
E = 2*D;             %Compensate for loss of energy by taking only 1 half of spectrum
F = E/sqrt(2);       %Rms of amplitude

SPL = 20*log10(F/Pref);

As far as I can see, Fluent uses the Power Spectral Density of the FFT to calculate SPL, instead of using the amplitude. The results of this calculation corresponds to what Fluent does in its FFT:

A = fft(s10);                 %Amplitude
G = A(1:fix(L/2)+1);          %Select 1-sided spectrum
PSD = (2*(abs(G)).^2)/(Fs*L); %Compensate for loss of energy by multiplying by 2; scale by
                              %dividing by Fs and L; square
SPL2 = 10*log10(PSD/(2e-5)^2);

When plotting these two different sets of results, there is a constant difference between the SPL values of the two signals. In this case it's a 12.0377 dB shift.

My question is, which of these two methods is the correct one to use in order to calculate SPL, and why? Or have I made a mistake in my calculations?

Help would be much appreciated!!

In amplitude method, you didn't divide by sampling frequency. In PSD method you divided by sampling frequency. This is the only difference between the two methods. I think because of that, there is a difference in the value. This may help. Suggest me, if I'm wrong.


Pilot tone frequency doubling U.S. Broadcast FM stereo includes a $19\textrm{ kHz}$ pilot tone that needs to be frequency doubled to help demodulate the DSB stereo difference signal.

How can this frequency doubling be done digitally?
What signal paths need to be delay matched (etc.) to avoid phase errors causing loss of stereo separation?

You can multiply it (modulate it) with itself, which will give you a signal with double the frequency and a DC component which is trivial to remove.

The other way to do it would be to $h * |x(n)|$ with $h$ being a resonator (and $x$ your pilot tone). That is, having tracked the pilot tone, rectify it and then pass it through a very simple resonator at 38 kHz to get rid of harmonics. In this case, you will only have to phase match the filter's delay.


How exactly is analog Signal to Noise Ration (SNR) defined? Say, my ADC reads noise signal only, not sure how to define a threshold. Let's say my noise is sitting at 0V (zero mean), oscillating between +5 mV to -5 mV, so Vnoise = 10 mVpp. I'm looking for 20 dB SNR with target present. Does it occur when ADC reads either 10^(20/20)(5 mV) = 50 mV or 10^(20/20)(-5 mV) = -50 mV OR when it's 10^(20/20)(10 mV) = 100 mV or 10^(20/20)(-10 mV) = -100 mV?

Now, a different scenario, what if my noise oscillates between -5 mV to +5 mV but there's a noise spike at +7 mV and -6 mV (just 2 values out of 10000 samples collected), do I change my threshold based on maximum noise level, not mean noise level?

It is not clear from your question if you are referring to ADC measurements using volts (typical would be counts) or if you are referring to the ADC input. However your question on SNR can apply to either case. Your reference to a "threshold" is unclear.

The signal to noise ratio is typically referred to the signal power, which occupies a bandwidth, as a ratio to the noise power over that same bandwidth.  Any other noise we have the opportunity to filter out in our subsequent processing.

Also, the noise power with reference to a voltage reading would proportional to the square of the rms voltage (1 sigma) not the peak-peak (which could be 6 sigma or more).

Further, this rms voltage would be representative of the (square root of the) power across the entire measurement bandwidth (in terms of an analog measurement), or the digital Nyquist bandwidth in terms of an ADC output (and again here I would use rms counts, where counts are increments of the lsb of the ADC).

So to compute SNR, if you are using rms quantities, it would be 20 Log the ratio of the rms signal to the rms noise in band. If you are using rms squared (power or variance) quantities, it would be 10 log the ratio of the variance of the signal to the variance of the noise is band (the portion of the noise that occupies the same bandwidth of the signal).


Power Spectral Density of Digital Modulations I want to simulate the power spectral densities of different modulations in MATLAB (without Simulink). Unfortunately I have only very basic understanding of signal processing, hence I struggle with this topic. Currently I have a very simple spectrum like this:
 There are different objects (A and B in the picture above) sending signals on different frequencies. I am not simulating any bits or calculating a sine, I am just assuming that there is something going on on a certain frequency with a certain power level.

Now, I want to simulate different digital modulations, e.g. PSK, BPSK, GMSK, FSK. I understand that the spectrum should look something like this:

There is not a single center frequency, but also side bands. A and B shall illustrate different modulations. My question now is, how can I calculate these different modulations? In text books I may find figures displaying the 'power spectral density' of ASK, but often times there is no formula how to calculate this. I guess the spectral density depends on the actual bit sequence, but I just need the envelope.

In one text book ("Digital Modulation Techniques" by F. Xiong) I found this formula for ASK:
  $$\Phi_s(f) = {A^2T \over 4} \mathrm{sinc}^2\big[T(f - f_c)] + {A^2 \over 4} \delta(f - f_c)\big]$$
where $T$ is bit duration, $A$ is carrier amplitude, $f_c$ is carrier frequency.
At first I found this formula promising, because it does not require a bit sequence, but I do not understand the last addend.

Are there any formulas like this for other modulation techniques?
If not, how could I calculate these graphs?

At first I found this formula promising, because it does not require a bit sequence,

That's the point about PSDs: They are an expectation value of a stochastic process! That is, the randomness (in this case, the data you're transmitting) is being "averaged away".

Deriving the spectral shape of any modulation (including the pulse shaping, don't forget that – you don't want to end up with a result like the formula you cite, as that is hardly band-limited!) is actually not that hard: You find a formula for the autocorrelation of the signal, and fourier transform that, and then you find the expectation. The trick is in finding a suitable form of the autocorrelation formula!

  but I do not understand the last addend.

That's a speciality of ASK: Most constellations are designed so that the mean of the information symbols is 0; that's not the case for ASK.

To make it short: The PSD $\Phi_{vv}$ of the baseband signal $v$, consisting of a sequence of information symbols $I_n$ with variance $\sigma_i^2$ and mean $\mu_i$, modulated at a symbol period $T$ with a pulse $g(t)$ (the pulse's spectrum being $G(f)$) is (from Proakis "Digital Communications", 3rd ed., eq. 4-4-17, on page 207):

$$\Phi_{ii}(f) = \frac{\sigma_i^2}{T} \left\lvert G(f) \right\rvert^2 + \frac{\mu_i^2}{T^2} \sum\limits_{m=-\infty}^{\infty} \left\lvert G\left(\frac mT\right) \right\rvert^2\delta\left(f-\frac mT\right)\tag{4-4-17}$$

where $\delta$ is the dirac delta distribution.

This applies to all linear modulations – PSK, QAM, ASK, OOK, but not to non-linear ones (FSK etc).

In the case of the example you cite, a rectangular pulse shaper $g(t)$ was chosen, which is $1$ for the period $-\frac T2, \frac T2$ around a "symbol time", and $0$ else. That one has a $\DeclareMathOperator{\sinc}{sinc}\sinc$ Fourier transform shape. That's where your $\sinc^2$ comes from!

The second term, the $\frac{A^2}4 \delta$ is what the sum in eq. $(4-4-17)$ collapses to, because the $\sinc$ is zero in every $mT$ but for $m=0$.

As said above, most constellations are designed to have zero mean ($\mu_i = 0$), so that the sum term just disappears, and your spectral shape is only defined by the spectrum of your pulse. (Only ASK and OOK differ, as far as commonly used modulations are concerned, as far as I can tell.)

Your formula uses $(f-f_c)$ instead of $f$, because it's already the bandpass signal's shape. Don't do that to you! Just consider the signal spectra in baseband, because we can mathematically shift those anywhere we like (especially by an $f_c$).

If you don't have access to good literature: Proakis 3rd ed can be had only for less than a fancy coffee at an American coffee house chain at online used book sellers. It's not my favourite edition of the book, as it cuts the math parts where the second edition really excelled, but it's still a very solid book, and as far as I can see, the cheapest edition of the "very established" digital comms books (being Proakis, Sklar, Tse/Viswanath, and, if you're German, Kammeyer). I like Lapidoth's book, too, and it can be had for free on his website. (update: I just found out that you can get Tse online, too, which is awesome! Might be a bit too advanced, though.)


Ornstein Uhlenbeck with drift The Ornstein-Uhlenbeck (OU) process

$dX_t = -\frac{1}{\mu} X_t + \sqrt{\frac{2\sigma^2}{\mu}} dW_t $

generates coloured noise with autocorrelation function $R(t) = \langle X_t,X_{t'}\rangle = \sigma^2 e^{-|t-t'|/\mu}$.

I would be interested in modelling a process $ Y_t = \alpha t + X_t$. After some searching, I came across the terms "OU process with drift" or "trending OU process".

Does anyone know what is the autocorrelation function of the process $Y_t$? Thank you!

Hi: This is a statistics viewpoint and not a DSP one. I just say that because there are termimology issues that have caused confusion with my answers in the past.  The autocorrelation is only defined for constant mean processes. A process with drift doesn't have a constant mean so you would need to detrend the process first which would bring you back to the zero mean case. One way to do that would be run a regression with time as your "x", in order to estimate $\alpha$. Then, once you have the estimate, you can subtract the trend off. I'm not sure if this is the answer you wanted to hear but the trend makes it not possible to estimate correlations because they can't be calculated because the mean of the process is not constant.


Calculate lines of the Inverse of Autocorrelation Matrix I need to calculate the inverse of a autocorrelation matrix

$$\mathbf R_{xx} = E\left\{\mathbf x \mathbf x^T\right\}$$

Where the samples $\mathbf x$ are $266000\times 1$ vectors, which means I'll have a $266000\times 266000$ matrix and need to invert it. But I only need the first 729 lines. The problem is, I have to find a classifier $\mathbf h$, where
$$\mathbf h = \mathbf R_{xx}^{-1}\mathbf u$$

where $\mathbf R_{xx}$ is the autocorrelation matrix and $\mathbf u$ is a $266000\times 1$ vector.

I have $\mathbf u$ and I can calculate $\mathbf R_{xx}$, but I can't save it on memory (it's a matrix larger than 400GB and I only have 380GB available). I'm only interested in the first 729 lines of $\mathbf h$

Use a rank one Cholesky update approach.


Confusion in deriving formula for fourier tansform of impulse train I was trying to derive fourier transform for impulse train :

I know how to solve for this using using properties of fourier transform. But now I wanted to use a brute force approach to it so I did the following

What makes you think it's wrong? You're just not done yet. In your final expression, if $\omega n T $ is a multiple of $2\pi$, you'll sum "infinitely many ones" which gives you "infinite", whereas for other values, you're going to some numbers "equally spaced" over the complex unit circle, which gives you zero. This is exactly how a train of deltas behaves.

Of course this argument is very handwavy and somewhat dangerous as these intuitions may mislead us when making statements about infinite sums.

As a suggestion: try with a pulse train of the form $\sum_{n=-N}^N \delta(t-nT)$. Compute its Fourier transform. Then consider the limit $N\rightarrow \infty$. You should get something like $\sum_{n=-N}^N {\rm e}^{-j\omega n T}$, which you can simplify using the geometric series. Then, forming the limit over $N$ should be easier.


Deconvolution and Polynomial division I found this comment inside MATLAB's deconv.m function

  Deconvolution and polynomial division are the same operations as a digital filter's impulse response $B(z)/A(z)$.

What does this statement mean?

Does this mean that one can recover the original signal, $x[n]$, given an impulse response $h[n]$ and output $y[n]$? If so, how?

Is it correct that if I subject an impulse signal $\delta[n]$ as input to a filter with the coefficients of $y[n]$ as the denominator and the coefficients of $h[n]$ as numerator I will get the original $x[n]$?

In addition to the answer from Matt L., I want to emphasize the fact why the comment mentions the polynomial division. Consider to finite-length sequences $h[n]$ and $x[n]$, with length $N$ (one sequence can be zero-padded to length $N$).

Now, look at the convolution of both:

$$y[n]=h[n]*x[n]=\sum_{n'=0}^{N-1}h[n']x[n-n']$$

As an example, consider the two sequences $h[n]=[1, 2, 3]$ and $x[n]=[4, 6, 5]$. Their convolution is (in Matlab)

>> conv([1 2 3], [4 6 5])

ans =

     4    14    29    28    15

Now, consider the product of the two polynomials $h(x)=1+2x+3x^2$ and $x(x)=4+6x+5x^2$ (in Mathematica):

Expand[(1 + 2 x + 3 x^2) (4 + 6 x + 5 x^2)]
>>> 4+14 x+29 x^2+28 x^3+15 x^4

As you can see, the coefficients of the polynomial of the product are equal to the values of the convolution of both sequences.

So, now to get $x[n]$ from $y[n]$, you need to perform the polynomial division of $y(x)/h(x)$:

PolynomialQuotient[4 + 14 x + 29 x^2 + 28 x^3 + 15 x^4, (1 + 2*x + 3*x^2), x]
>>> 4+6 x+5 x^2

Hence, the coefficients of the polynomial yield the original sequence.

This technique can be extended to infinite sequences.

To get the connection to the Z-Transform here: The Z-Transform of the finite sequences $x[n]$ and $h[n]$ are given by

$$H(z)=1+2z^{-1}+3z^{-2}\\ X(z)=4+6z^{-1}+5z^{-2}$$

And the Z-Transform of the output $y[n]$ is given by

$$Y(z)=H(z)X(z).$$

Here, you again have the polynomial multiplication. Also, here it becomes even easier to spot, why it also works with infinite sequences: Then both $H(z)$ and $X(z)$ become fractions of polynomials, but the principle of the polynomial multipication remains.


Determining what filter (if any) has been applied to an image I'm trying to determine what "filter" has been applied to an image.  I have the original, raw image and the processed image.  The processed image has been put through a pipeline (black box) that does various things to the image in an iterative manner.

The filter, as far as I can tell, is similar to an unsharp-mask, but with "negative bowls" mitigated.  However, the processed image cannot be reproduced by a simple unsharp mask.

I've attempted to measure the "filter function" by taking the ratio of the power spectra of the raw and unprocessed images, but the best I've been able to do with that is unsharp mask with a non-gaussian smoothing function.

What methods are available for determining the filter applied to an image?

(not related to Determining the type of filter on an image)

In case your filter is linear, this should be quite easy. You have a linear model with unknown parameters (9 if your support is 3x3, 16 if 4x4, etc...). You have a lot of constraints - the original image and result image. You can use linear least squares. If you can generate any input, it is even easier. Recall that a linear time-invariant system is fully determined by its response to delta impulse, which is its kernel.
In case your filter is not linear (like Paul says), there is almost nothing that you can do to find it out, except using your common sense, or attempting to optimize another model, given your inputs and outputs. Image processing logic can be very non-linear, including ad-hoc solutions, so it looks quite hard to me.


MSK signal recoverd though FM Discriminator? I am attempting to write a FFSK (differentially encoded MSK) modulator.  I have found and implemented a demodulator that works with sample data taken from a radio of a FFSK signal.  The demodulator uses a FM Discriminator to transform the IQ data to a FFSK signal that can be decoded.  I am using the method described in the article "Minimum Shift Keying: A Spectrally Efficient Modulation" by Subbarayan Pasupathy to modulate a FFSK signal.  When a FM Discriminator is applied to the output of the modulator a square wave alternating between +/- 4.5 deg is produced.

Below is a snippet of my matlab code used for modulation:

% differentially encode to produce a 1-to-1 mapping between input bits and output frequency for FFSK
diffBits = [bits(1) zeros(1, length(bits) - 1)];
for i=2:length(diffBits)
    diffBits(i) = xor(diffBits(i - 1), bits(i));
end
iBits = diffBits(1:2:end);
qBits = diffBits(2:2:end);

% rectangular pulses of duration 2*T w/ Q delayed from I by T
iPulses = [kron(iBits * 2 - 1, ones(1, 2 * SYMBOL_SAMPLES)) zeros(1, SYMBOL_SAMPLES)];
qPulses = [zeros(1, SYMBOL_SAMPLES) kron(qBits * 2 - 1, ones(1,2 * SYMBOL_SAMPLES))];
t = -T:SAMPLE_PERIOD:(MESSAGE_DURATION - SAMPLE_PERIOD);
plotTime = t/T;

% baseband symbols (Parallel generation of Type I MSK, see Couch Figure 5-36)
iSymbols = iPulses .* cos(pi * t / (2 * T));
qSymbols = qPulses .* sin(pi * t / (2 * T));

Here is a graph of the output of the FM discriminator from the sample data:

Here is a graph of the output of the FM discriminator from the modulator data:

My question is what is different between the real data and my modulator that could produce this different output?  I am looking for a way to make my modulator produce data similar to the sample data.

If this helps, here is the FFT of the baseband of the sample data (derotated and resampled down):

% samples are complex int16 values read in from a file
num_samples_fft = 4096*2;
Fs = 781250;
freqOffset = -25200
shifted_samples = samples .* exp(2*j*pi*freqOffset/Fs*(0:length(samples)-1));
resampled_samples = resample(shifted_samples, 96, 3125);
resampled_samples_fft = fft(resampled_samples,num_samples_fft);
resampled_samples_fft = resampled_samples_fft([num_samples_fft/2+1:end 1:num_samples_fft/2]);
figure;
plot(f,2*abs(resampled_samples_fft(1:num_samples_fft)));
title('Sample Data FFT');

Whereas if I do the same process to data that I modulate this is the resulting FFT graph:

It is possible that the sample data has a pulse-shaping filter before the frequency modulator. That would explain the first plot.

Look at the power spectrum of the sample data and check the sidelobe levels as compared to textbook MSK and your modulator. A pulse shaping filter will reduce the sidelobe levels.

Square waves are expected in the second plot, because the data entering the modulator has no pulse shaping on it. The spectrum should be textbook MSK.


Signal enhancing algorithm (deconvolution??) I need an algorithm that will detect the correct peak location in a "signal" (actually a series of data points).
The signal is generated by a machine. The signal is nothing similar to an electrical signal (audio/radio frequency).

The data points could look like this: [1, 2, 1, 3, 4, 5, 4, 3, 2, 1, 13, 14, 15, 18, 23, 19, 17, 15, 15, 15, 14, 11, 9, 4, 1, 1, 2, 2, 1, 2]

I think we could imagine that the signal was accidentally passed through a low-pass filter, and now we want to restore it.

NOTES:

After a peak, if the signal is not coming down fast enough we can consider that there are multiple peaks (you can best see that in the 'red' signal at the end of the series).

One advantage on my side is that the whole data series is available (save on disk). So, the signal does not have to be processed in real time.

I can see in your Signal on the right, the "picks" are separated equidistantly. So this is my proposition: Subsample your signal and add a amplified Gaussian signal centered at each sample.


Receiving multiple signals at the same frequency at the same time I'm doing a project for my school it's a communication platform for different mobile devices I'm gonna build and I need to send data from the devices to station am making, I figured how to use signal spreading to encode multiple signals and send it to multiple devices like a cell tower.

But what I wanna know is how the tower I'm building is gonna receive multiple signals and not cancel each other out?

One basic concept is that if multiple signals are mutually orthogonal in some sense (time domain, frequency domain, encoding, spreading function, phase relationship, etc.) they won't interfere in that domain, and thus can't cancel.


Transmission Bandwidth and Sampling frequency In LTE standards, for a transmission bandwith of 2.5 MHz, the sampling frequency is 3.84 MHz. I was wondering with this choice, are we violating Nyquist theorem?

I found the details in the following link (page 15)

http://www.freescale.com/files/wireless_comm/doc/white_paper/3GPPEVOLUTIONWP.pdf

I would have assumed that with a bandwith of 2.5 MHz, then the channel sampling frequency (signal input sampling frequency) should be at least 5 MHz?

Thanks.

Complex (I/Q) signal processing is used in this case. With a complex-valued signal sampled at a rate $B \text{ Hz}$, you can unambiguously represent a total bandwidth of $B \text{ Hz}$.

What you're thinking of is the real-valued signal case, where a signal sampled at $B \text{ Hz}$ can unambiguously represent a total bandwidth of $\frac{B}{2} \text{ Hz}$. This is the most common statement of the Nyquist theorem, but it only holds for real-valued signals.

Note that the total number of "values" per second is the same, as each complex sample consists of a real and imaginary part. So, there's no "information gain" to be made by using complex sampling, it's just more convenient for many applications (including digital communications signal processing).


Example image that is degraded by shot noise and easily misinterpreted? For illustrating the effects of shot noise (more precisely, emitter failure in fluorescence imaging) I am looking for a good example image. The image should ideally be easily misinterpreted, especially by non-aided human analysts, when a fair amount of shot noise is applied to it.

Does anyone know of such an image?

Would something like Figure 4 in this paper by Dumitriu, Rodriguez, and Morrison (2011) work? They imaged the spines on a neuron using six different gain settings, which produce different amounts of shot noise. One can see how the putative shape of the spine changes based on the amount of noise present. This example might be a little esoteric for a general audience, but it does sound like you're doing some kind of microscopy, so perhaps it's a good fit.

Dumitriu, D, Rodriguez A, and Morrison JH. High-throughput, detailed, cell-specific neuroanatomy of dendritic spines using microinjection and confocal microscopy. Nat Protoc. 2011 Aug 25;6(9):1391-411. doi: 10.1038/nprot.2011.389.
Article at Nature Protocols


Determining SINAD from measured data I would like to determine the SINAD of a distorted signal from a measurement and have troubles to get consistent data. The problem is (at least to my understanding) that I only have a very limited number of signal periods and I'm not sure, whether I have enough data to provide this measurement at all.
In order to determine the SINAD I need the power of my signal as well as the power of "all the rest", i.e. the power of noise and distortion. To obtain these powers, I took my oscilloscope data and calculated the spectrum via an FFT. I didn't use any window function (i.e. I used a rectangular window) to calculate the FFT since I know the fundamental frequency I expect from my signal and can choose an integer multiple of the fundamental period as my window width.
From the frequency bin at my desired frequency I determined the fundamental amplitude. The noise + distortion power was then calculated by subtracting the signal's RMS value from the total RMS value according to

From these two RMS values I then calculated the SINAD according to:

The problem is now that my results vary depending on the number of fundamental periods (I took between 1 and 33 periods, I don't have more) from the measurement I take to calculate the FFT, where increasing the number of fundamental periods degrades the SINAD.
This makes sense to me to a certain extent since if I extend the window size of my input data, the frequency resolution of my spectrum increases. This means that parts of the noise + distortion in the vicinity of the fundamental bin which were combined into this bin for a course frequency resolution are now contained in the neighboring bins of the fundamental and treated as noise/distortion and not as part of the signal anymore.
I couldn't really find a lot on how to eliminate this issue. Googling on how to measure SINAD mostly leads to the definitions themselves. Is there a way to determine the SINAD for the limited amount of signal periods? How would I have to calculate it? I also tried the matlab function sinad(), however, it has similar issues (although with this function, the dependence on the number of signal periods is less compared to my own calculation).

Similar to the determination of the confidence intervals when estimating the population variance from a smaller sample, with a limited amount of data you will have a large variability in your estimates. There really is no way around reducing this error or variability other than taking more data (to the extent it is stationary).
For an optimized approach to measuring the SINAD, I suggest using correlation of your distorted signal with the reference waveform before any noise is added. I detail this approach further at these post:
How can I find SNR, PEAQ, and ODG values by comparing two audios?
Noise detection


Picking out a signal that appears to be noise inside a large signal I have a signal from a photodiode sensor that has two types of noise. One type of noise is ambient light white noise that gets introduced just from the surrounding environemnt. The other type of noise is due to the motion of the sensor.

Currently i don't seem to have a problem filtering out the white noise when the sensor is held still. But when the sensor is moving i get very large noise from motion, this noise is much larger than the signal i am interested in and as a result makes the signal I want look like noise.

The signal i am interested is anywhere between 0.75Hz and 4Hz, the motion signal can be anywhere from 0.25Hz to 5Hz.

I am currently using a butterworth filter between 0.5 and 4 Hz and then a moving average to smooth it out. But this causes the signal to be lost in the motion noise.

I have looked into adaptive filters but am not sure whether i should try to implement them or not because i do not have enough experience with them to tell whether they will help my problem or not.

Are there techniques that are used to pick out signals that share frequency with noise?

Adaptive filters could be used in situations where you have a reference signal, what adaptive filters basically do is to try modelling the signal that is present in captured signal using reference signal. So, we can subtract and/or extract the required signal.

In the situation mentioned in question, you can use adaptive filter if you have reference for either the signal you are interested in or the noise you want to remove. If you have the noise signal as reference adaptive filters will remove the required signal(signal of interest) from input and only noise will be remaining, now this noise can be subtracted from original input signal to get the required signal(signal of interest).

Same could be done if signal of interest is present as a reference, but in this case the adaptive filters output will be the required output.


When input data into a low pass filter, are the first several data unusable? For example, I want to feed a DC into a low pass filter. The the filter's first several output will not have the same DC value I want.

Does it mean that I need to dispose first several results?
How many shall I dispose?
What is the technical term describing this phenomenon, step response?
Does it only happen to DC or to other frequency signals as well?

Here comes sample code and diagram:

Fs = 64e3;
f = 8e3;
T = 1/Fs;
t = 0:T:0.002;
input = 0.002*cos(2*pi*f*t);
input_with_DC = input+0.5;

[b,a]=cheby1(1,2,0.4);
output = filter(b,a,input_with_DC);

plot(input_with_DC)
hold
plot(output,'r')

It can be seen from above diagram that the DC rises at first and then approach to the actual DC value.

You have to judge yourself if for your purposes you need to get rid of the first few output values. The phenomenon you observe is determined by two factors. The first is the delay of the filter (which is usually frequency dependent). This delay is a consequence of the causality of the filter. For a linear phase FIR filter that delay is independent of frequency and equals $(N-1)/2$ samples, where $N$ is the filter length (as pointed out in a comment by MBaz). The other factor lies in the nature of low pass filters and has nothing to do with a non-ideal (causal) implementation. It is the non-zero rise time that any low pass filter exhibits, even an ideal (non-causal) brick-wall filter. The rise time is proportional to the inverse of the cut-off frequency.

For minimum-phase filters with a very low order (as in your example), the first factor plays only a small role, so the delay you observe is mainly caused by the non-zero rise time of the low pass filter. If you decrease the cut-off frequency you will observe an increase in rise time and vice versa. Note that the input signal that your filter sees is actually a jump from zero (no input signal) to a DC value (plus some small signal), so you pretty much see the system's step response.

In general the response to a change in the input signal is called transient response, and it is not restricted to a change in DC value.


OFDM: How to do the carrier recovery of the quadrature(I/Q) demodulation For an OFDM system, after the OFDM modulation, the quadrature modulation is used to move the signal from the baseband to the pass-band.
But when it comes to the receivers, I can't find any way to make the carrier recovery, not only the frequency but also the phase.

I can't find any way to make the carrier recovery,

I don't think you can or should really talk "carrier recovery" in OFDM systems; after all, OFDM is a multi-carrier system. So, you'd need to recover some 2 to a couple thousand carriers, depending on the OFDM system's number of subcarriers.

  but also the phase

Uh-oh. Have you considered what this means? When do we use OFDM?

We use OFDM when we have a frequency-selective channel. Which means that the phase is different for different frequencies within your band. Which means you can't say "it's a single phase that needs to be recovered". You need to synchronize all phases of all subcarriers, if you want phase synchronization.

Essentially, you try to use single-carrier terminology on a multi-carrier system. That won't work.

What we should do, here, is talk about frequency and phase synchronization.

In OFDM systems, there's different ways to do that, but the one you'll most likely encounter (because it's very commonly used in bursty/packeted systems like Wifi) is a pretty clever thing called "Schmidl&Cox synchronization" (read the excellent original paper [1]!).

Basically, what that does is a three-step process:

Use the autocorrelation properties of the specifically crafted synchronization OFDM symbol to estimate the beginning of a frame
Use the phase that an autocorrelation coefficient of a symmetrical part has to determine what is called the "fine frequency offset", i.e. the frequency error modulo subcarrier spacing
Determine the "coarse frequency offset" simply by counting the "energy-carrying" subcarriers.

Note that the above is one way to synchronize OFDM; there are others, and they depend on other mechanics.

Phase synchronization in OFDM systems is most commonly done using phase reference signals ("pilot symbols") spread across the time / subcarrier plane. You just define "the symbol on that subcarrier at that time is not a data symbol but always has the value $X$", and then you look what you receive instead of $X$, correct the phase of that so that you do get $X$, and use that phase correction factor.

But that's again system-specific – there's enough OFDM systems that don't do phase sync at all, but use a differential modulation inside (DAB / DAB+ for example; that one oesn't do S&C, either).

To conclude:

There's no "carrier recovery" in OFDM, since there's not a single carrier
The equivalent "frequency synchronization" and "phase synchronization" is something that is very system-specific, and there's a lot of different way it's done – and all depend on properties of the transmit signal and the constraints of the receiver. There can't be a single answer to the question "how to sync OFDM".
If you want to learn one way to sync OFDM, it's probably Schmidl&Cox[1].

[1] T. M. Schmidl and D. C. Cox, "Robust Frequency and Timing - Synchronization for OFDM," available online


Advice on filtering a CW signal for processing This question is related to my other recent questions.  In the last question, I learned that I should filter about 4 * WPM (words per minute) around the carrier frequency if I want to recover Morse code.  One problem is, learning the words per minute.  I filter because I want to differentiate times when the system is keyed, to times that are just noise.  Filtering requires WPM, WPM seems to require filtering.

Now I am simply going through the whole audio, and if a window has an element above a certain threshold, the whole window is marked as being 'on' or 'keyed'.  The window is big enough for a whole alternation of the carrier to be examined (sampling frequency over carrier freq, or fs/fc).

Are there better ways to approach this?  I have been working on this for awhile now.  I am doing this in Python.

The recommendation to use 4* WPM just means that you should estimate the upper limit of the WPM that you expect to receive.  High speed CW starts around 40 WPM (usually too fast to write or type the decoded signal.) So, 50WPM would be a good upper limit which means you could use a bandwidth of 200Hz.

As an alternative to what you are doing, you could use the goertzel algorithm.  That algorithm works like a narrow bandpass with a signal level detector.  You run your signal through it, and it delivers a number that tells you how strong the signal is at the center frequency that you gave the goertzel method.

This web site  gives a clear example of how to write a goertzel method to do what you are trying to accomplish.  Use the tandem goertzel to work in real time.  The example is in C, but it should be clear enough that you can translate it to python.

Varying the RESETSAMPLES value will allow you to change how narrow the filter is.  The higher RESETSAMPLES is, the narrower the filter will be.  The article suggests 200 samples, so I would start there.  If your signal is too noisy and you get a lot of bad morse code symbols, then you can make the number larger.  What I usually do is make it some multiple of the wave length of the frequency I am trying to detect in samples.  So, to detect 1kHz using 44100Hz sampling, I would first try 44 (44100/1000.) If that is too "loose" then 88, then 132 and so forth until I get a filter that is narrow enough.  If you go too high, though, you will start "smearing" the morse code symbols and won't be able to tell when one ends and the next begins.


Difference Equation - Time Response analysis of Discrete Systems We know that that the difference equation formula for computing an output sample at time $ n$ is based upon its past/present input samples and past output samples in time domain:

$$
\begin{aligned}
y(n) &= b_0 \,x(n) + b_1 \,x(n - 1) + \cdots + b_M \,x(n - M) \\
&\qquad- a_1 \,y(n - 1) - \cdots - a_N \,y(n - N) \\
&=\sum_{i=0}^M b_i \,x(n-i) - \sum_{j=1}^N a_j \,y(n-j)
\end{aligned}
$$

My question is why $j$ begins from 1 rather than from value 0 like for $i$ and also how can we write the below output equation in difference equation form:

$$
\begin{aligned}
y(n) &= 0.01\ x(n - 5) + 0.002\ x(n - 1) + 0.99\ y(n - 1)\\
y(n) &= 0.01\ x(n) + 0.002\ x(n - 1) \\
y(5) &= 0.01\ x(3) + x(1) \\
\end{aligned}
$$

Secondly why $y(n)$ sequence values are calculated as a sum of a series of input / output, cause I have been studying that $$y(n) = \{2,4,6,8\}$$ when $$x(n) = \{1,2,3,4\}$$ and $$n =\{1,2,3,4\}$$ (example per say) i.e. $y(n) = 2\cdot x(n)$.

There was no mention that for $y(n)$ ($y(2)$ say) value we need to add values of $x(n)$ and $y(n)$ i.e. (take values of $x(1)$ or $y(1)$) or $y(n-1)$)?

The answer to your first question (why $j$ starts at $1$) is contained in your first sentence: the output is computed from the current input value ($i=0$), from past input values ($i\ge 1$), and from past output values (i.e. $j\ge 1$). Imagine that you had the following difference equation with $j$ starting at $0$:

$$y(n)=\sum_{i=0}^Mb_ix(n-i)-\sum_{j=0}^Na_jy(n-j)\tag{1}$$

Then you could rewrite (1) as

$$y(n)=\sum_{i=0}^Mb_ix(n-i)-\sum_{j=1}^Na_jy(n-j)-a_0y(n)\tag{2}$$

and consequently

$$y(n)(1+a_0)=\sum_{i=0}^Mb_ix(n-i)-\sum_{j=1}^Na_jy(n-j)\tag{3}$$

Now you can divide (3) by $1+a_0$ to arrive at an equation like the one given in your question (with $j$ starting at $1$) (and with all coefficients scaled by the constant $1+a_0$). So you can always use the difference equation with $j$ starting at $1$ without loss of generality.

Concerning the 3 examples of difference equations you gave I'll just discuss the first one, the others are totally analogous. Here you have $M=5$ and $N=1$ with coefficients $b_1=0.002$, $b_5=0.01$, and $b_0=b_2=b_3=b_4=0$, and $a_1=-0.99$.

Finally, concerning your last question about the system $y(n)=2x(n)$. That's again a special case of your general difference equation with all $a_i=0$ (i.e. no recursion) and with $b_0=1$ (and all other $b_i=0$, $i\ge 1$). So all your examples are special cases of the general difference equation.


Self study question about Padé Approximation and transfer functions Given a signal $x[n] = [1, 2, 3, 4, 5]$. How many transfer functions can be found with Padé approximation which have a causal impulse response and start with those five samples. How many of them are stable?

Found this question and I'm not sure if I understood it right. I'd say none of them is stable since Padé can't guarantee stable solutions. The first five samples will be perfectly recreated in the impulse response. Therefor $h(z)$ between $0-4$ looks like $x[n]$. But what about how many? Are there multiple solutions? Is this question badly formulated?

A further question is then: Find the Padé model with exaclty one pole. One pole? Isn't that an all pole model? Or is there some kind of rule like always the same amount of poles and zeros? Can't solve it with the Padé euqations.

The following might help. Given the data of $5$ samples,

$$ x[n] = [1, 2, 3, 4, 5]$$

You will have possible choices for the orders $(p,q)$ of $a[k]$ and $b[k]$, such as $\{(4,0),(3,1),(2,2),(1,3),(0,4)\}$.

You can solve or each of the possibilities and check whether they yield stable systems or not, by looking at the roots of $a[k]$... afaik they should all be causal as this is Pade modeling assumption (unless equations modified to handle otherwise).

Note that for certain data sets (including this one) some tail coefficients might turn to be $0$ and the actual order can be less than that of indicated by $p$ or $q$.

Also note that a model with one pole is not an all-pole model. All-pole model requires that $q=0$, and does not depend on $p$.


Using a low pass filter to interpolate signal In my DSP university textbook, the interpolation process is described as follows:

  In order to represent a baseband signal $x[k]$ at an increased sampling rate with the same shapes of its time-domain and frequency-domain representations, the insertion of zeros must be followed by low pass filtering with a total bandwidth of $\frac{2 \pi}{R}$ to perform the actual interpolation, with $R$ being the upsampling factor.

The process is depicted in the following figure:

I understand well how we insert zeros between samples of the original signal (i.e. how to go from the first to the second plot of the figure), but I don’t understand how applying a low pass filter to the signal we obtain in the second plot results in the third one. To me, it looks like those zero samples magically get a value, so I would appreciate any explanation of this “magic”.

Another view which results in superior interpolation filter design is revealed by reviewing the spectrum of the signal with zeros inserted (resulting in an increase of the sampling rate by $I$ when you insert $I-1$ zeroes) and you will see that the original spectrum is replicated at integer multiples of the original sampling rate. The ideal low pass filter is the one that can pass the original spectrum with no distortion and reject all the new alias copies perfectly. (Ideal is not achievable but this defines the filter design targets to approach that).  The ideal filter will interpolate the non zero values to grow the zeros to the value for the final interpolated waveform.

What is occurring, to explain the OP's question conceptually, is that the signal requires high frequency content when going from zeros to a large value in just one sample (and inserting zeros causes this to occur as evidenced by the additional high frequency content in the spectrum). This makes complete sense since the frequency content would be related to a change in magnitude versus a change in time (a big change in a short time requires high frequencies). A filter that removes the high frequencies results in the inability for the signal to change rapidly from one sample to the next (due to the memory of the previous samples which is what a filter does). A very easy way to see this occurring is with a simple moving average filter over $M$ samples when you interpolate by $M$ by inserting $M-1$ zeros; such a filter will hold the last sample as a zero-order hold, growing each zero to the last non-zero sample. This isn't a good way to do interpolation since only one non-zero sample is in the memory of the filter, but very easy to visualize. The "magic" happens with better designed filters as described in this answer which then consider many more samples to perform what is essentially interpolation with higher order polynomials.

With regards to Sinc reconstruction that is commonly used to explain "perfect" reconstruction; it is worth mentioning that Sinc reconstruction is only perfect when you can use a Sinc, but unfortunately the Sinc function extends to infinity which is impossible to do in practical realization (it is equivalent to filtering with a brickwall filter). Therefore the Sinc response is truncated in practice, which on its own results in a very poor filter, which is then typically windowed to improve performance. Ultimately, its use when truncated and windowed results in inferior reconstruction filters. Given the approach of paying attention to where the aliases actually reside in frequency, you can realize optimized filters (superior to simply truncating and windowing a Sinc function) using the multi-band filter algorithms such as least squares using firls() function that is available in MATLAB, Python and Octave.  I demonstrate this in the plots below from an X4 interpolation example with $f_s = 10$ KHz. The upper plot compares options for interpolation filters with the same number of taps, showing the poor performance of a truncated Sinc in green, and then the improved windowed Sinc filter in green along with the least squares multiband filter in red. All have the same number of taps, but the multiband approach results in 10 dB better rejection of the higher frequency images (and not visible from the scale of the plot, but lower passband ripple distortion as well).

The filter magnitude for the multiband filter was shifted down to overlay with the spectrum after zero-insert showing how the multiband filter maximizes rejection where needed thus resulting in the least amount of distortion in the resulting interpolated signal.

So the spectrum above with it's replication at the higher frequencies (images) has the zeros inserted in the time domain. When those higher frequencies are removed through filtering, the spectrum represents the original signal, at the higher sampling rate, as if it was sampled directly.

A truncated and properly windowed Sinc filter will not achieve the rejection or passband performance of this filter with the same number of taps since it would provide rejection at all locations in between the aliases, while the multi-band filter maximizes the rejection where it is needed. The resulting distortion is predictable and can be traded with the desired filter complexity. This is a standard interpolation design approach for high performance and efficient interpolation (especially when converted to polyphase structures!).


Practical Gardner timing error implementation I'm trying to demodulate a QPSK communication.
For the moment I do:

  Signal detection by computing the ratio between the power of the carrier wave compared to the signal power.
  Initial clock synchronization with a preamble.
  Initial phase recovery with the same preamble.

This works great with my signal but now I want to go further and take into consideration a clock frequency offset betweem the TX Clock and the RX Clock.

I made some research and a good way to do it seem to be Gardner Timing Error Detector which should be able to compensate for my clock offset.

To implement it I do:
$$\epsilon_k = Real[(signal(t_k)-signal(t_{k-1}))*conj(signal(t_{k-1/2})]$$
$$t_k = t_{k-1}+T+\gamma*\epsilon_k$$

My problem is that I don't know how to chose the proportional parameter.
$$\gamma$$

I've sometime seen it with something called the normalized narrowband noise BlT but I have no idea what it is. If you could point me in the right direction any help would be very appreciated!

If you have access to Mengali's book 'Synchronization Techniques for Digital Receivers', you will find the details how to choose it for a tracking loop, be it for frequency, phase or timing. Otherwise, just read about a standard discrete-time phase locked loop and see how the choice of paramters affect the loop bandwidth, which is what you are after.


Split audio file into audio phrases I have couple of  audio recordings of length 5 -7 min. I also have corresponding text to each audio file.
I am not subject matter expert on speech to text but I am vaguely aware that much smaller audio files (10 sec) are used to fine-tune speech recognition model .
So I want to split audio files into meaningful partitions .
Correct me if I am wrong in assumption that I cannot train model with large audio files.

Some toolkits like Kaldi can handle longer inputs. Others indeed require 10 seconds or so. It depends on the tools you gonna use.


What is spatial correlation? What does the spatial correlation of a signal say about it?How is it different from autocorrelation?

`

Spatial correlation is relevant in sensor arrays and MIMO wireless communications. In the case of MIMO communications where multiple spatially separate antennas receive different observations of a multipath signal, spatial correlation can be used to evaluate the 'antennae diversity'. Higher the diversity, better the mitigation of signal fading effects introduced by the wireless medium. As a result the link quality is improved.
http://en.wikipedia.org/wiki/Antenna_diversity


matlab How to design lpf\bpf\hpf without builtin functions related to Filtering voice signal from randomly added frequencies

I'm trying to create a simple low pass filter but i'm not sure what is wrong.

my code:

[y,Fs]=wavread('sound.wav');
m = length(y); %30586
Fc = 500;               %% frequency cut
n = [0:m-1];            %% kernel length / samples
Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
Hm(1) = Tc/pi;
figure,plot(Hm);

1.

LPF: How to get to the 0-500 range?

We used this code in the class but with Fc=45,m=100,Fs=200 and it looks ok.
but when I try this on this wav file it's not working as expected. I thought it will pass from 0-500 but instead it returns some crowled plot.

only when I set Fc to be a low number I can see the sin wave. so i've try to manipulate by setting Fc=real(ceil(m/500)). it better than nothing but it's not it.

2.

hpf = 1-lpf, right? well is seems not to be. my guess is that my lpf is wrong

3.

I will deal the ripple after i'll manage to do this

You seem to be confused with the time and frequency domain. Your code for an LP filter is the inverse Fourier transform of an Ideal low-pass filter. So by getting a frequency response for the function you wrote you'll get an ideal LP response you want.

However, there are two mistakes:

1) Your frequency for the ideal LP filter needs to be $[-\frac{m-1}{2} : \frac{m-1}{2}]$

2) You need to change the value for the 0-frequency only if M is an odd number.

Now, you can use the function freqz to get your freqency response, so your final code would be like this:

%[y,Fs]=wavread('sound.wav');
%m = length(y); %30586
m = 30586;
Fc = 500;               %% frequency cut
n = -(m-1)/2 : (m-1)/2;            %% kernel length / samples
%Fs = Fs;    %22000      %% sample frequency / (nyquist - %2*m)
Fs = 22000;
Tc = 2*pi*(Fc/Fs);      %% theta c
Hm = sin(Tc*n) ./ (pi*n);
if (rem(m, 2) == 1)
    index = (m+1)/2;
    Hm(index) = Tc/pi;
end

[H w] = freqz(Hm, 1, 1024); % 1024 - or how many points you want for your Fourier Transform

% Now your H is your frequency response, and w is digital frequency in [0,pi].
% You can now plot the amplitude response:
figure,plot(w, abs(H));

And you have your Low-Pass filter. However, I'd like to suggest using the many ways to create filters in Matlab, they were made for a reason, and allow for much easier filter design :)

One more remark, I have to suggest that you know the difference between analog and digital frequency, as they have different notation. Analog frequencies you write in uppercase, e.g. $F, \Omega$, and digital frequencies in lower-case, $f, \omega$. Here's a version of the code which is a little prettier:

M = 30586;
Fc = 500;
Fs = 22000;
wc = 2*pi*Fc/Fs;

n = -(M-1)/2:(M-1)/2;
h = (sin(n*wc))./(n*pi);
if ( M-2*fix(M/2) ) > 0
        index = (M+1)/2;
        h(index) = wc/pi;
        display(index);
end

[H w] = freqz(h, 1, 1024);
plot(w, abs(H));

Cheers :)

EDIT:

When you use freqz, the digital frequency w can be somewhat hard to control, because it's not created in the same way, or more clearly, the points are not selected like a normal frequency we create with, for example, F = 0 : Fs/(M-1) : Fs, even though it should be the "same" frequency.
The result is that it's not really possible to get the output filtered signal with Y = H .* X, and  y = ifft(Y).

However, I would suggest using the function filter(b, a, x) for filtering your signal, which works good with freqz. You can proceed to do just that:

% ------
% Create a random signal x
n = 0:M-1;
F = -Fs/2:Fs/(M-1):Fs/2;
Fx = [2000 7000 10000];
x = cos(2*pi*Fx(1)/Fs*n) + cos(2*pi*Fx(2)/Fs*n) + cos(2*pi*Fx(3)/Fs*n);
X = abs(fftshift(fft(x, Nfft))); % fftshift because we need domain [-Fs/2, Fs/2]
% ---------

% -------
% Now filter the signal
y = filter(H, 1, x);
Y = abs(fftshift(fft(y)));
% -------

That's it :)


How do I interpret the DCT step in the MFCC extraction process? In most audio processing tasks, one of the most used transformations is MFCC (Mel-frequency cepstral coefficients).

I mostly know the math that's behind the MFCC: I understand both the filterbank step and the Mel frequency scaling.

What I don't get is the DCT (Discrete Cosine Transform) step: What kind of information do I get in this step? What is a good visual representation of this step?

You can think of the DCT as a compression step. Typically with MFCCs, you will take the DCT and then keep only the first few coefficients. This is basically the same reason that the DCT is used in JPEG compression. DCTs are chosen because their boundary conditions work better on these types of signals.

Let's contrast the DCT with the Fourier transform. The Fourier transform is made up of sinusoids that have an integer number of cycles. This means, all of the Fourier basis functions start and end at the same value -- they do not do a good job of representing signals that start and end at different values. Remember that the Fourier transform assumes a periodic extension: If you imagine your signal on a sheet of paper, the Fourier transform wants to roll that sheet into a cylinder so that the left and right sides meet.

Think of a spectrum that is shaped roughly like a line with negative slope (which is pretty typical). The Fourier transform will have to use a lot of different coefficients to fit this shape. On the other hand, the DCT has cosines with half-integer numbers of cycles. There is, for example, a DCT basis function that looks vaguely like that line with negative slope. It does not assume a period extension (instead, an even extension), so it will do a better job of fitting that shape.

So, let's put this together. Once you've computed the Mel-frequency spectrum, you have a representation of the spectrum that is sensitive in a way similar to how human hearing works. Some aspects of this shape are more relevant than others. Usually, the larger more overarching spectral shape is more important than the noisy fine details in the spectrum. You can imagine drawing a smooth line to follow the spectral shape, and that the smooth line you draw might tell you just about as much about the signal.

When you take the DCT and discard the higher coefficients, you are taking this spectral shape, and only keeping the parts that are more important for representing this smooth shape. If you used the Fourier transform, it wouldn't do such a good job of keeping the important information in the low coefficients.

If you think about feeding the MFCCs as features to a machine learning algorithm, these lower-order coefficients will make good features, since they represent some simple aspects of the spectral shape, while the higher-order coefficients that you discard are more noise-like and are not important to train on. Additionally, training on the Mel spectrum magnitudes themselves would probably not be as good because the particular amplitude at different frequencies are less important than the general shape of the spectrum.


How does Jpeg encoding affect the final file size of 24BPP RGB and 8BPP grayscale image? The reason I'd like to know this is because of a file size problem. I am trying to reduce an image's file size from converting it from TIFF to JPEG. A user challenge me that he notice a RGB TIFF image is three times the size of the same image but in grayscale TIFF version, but the post-convert grayscale version of the JPEG image is almost the same size as the RGB JPEG version.

He insist that a 8BPP grayscale jpeg should have ~1/3 of the size of the 24BPP RGB jpeg.

Then I have been told by some technical adviser afterward that JPEG will ever only exist in the form of 24BPP file, and that the grayscale jpg images I found online are all lying in their image meta data saying bit depth is 8. But so far I haven't found any proof online that confirm nor deny that claim.

At this point I have no idea if my conversion method is flawed and inflated the file size, or the fact that the so call "8BPP JPEG that takes only 1/3 of the file size" doesn't exist from the beginning.

Is there any documentation regarding this issues online?

Take a standard R-G-B coded 24-bit image. First, suppose that all channels are coded independently. If by chance, R, G, and B are equal, each one will be converted to the same JPEG file with size $s$.

But such an image would be visually grayscale. The corresponding JPEG compression of the grayscale would be about 3 times smaller (taking headers, etc. into account).

But in JPEG (and most image compression formats), color channels are not treated independently. They are first decorrelated into one luminance (L) and two chrominance images (U, V). If all channels were equal, ie if the image were grayscale as above, Y is the natural grayscale, U and V are zero. The latter two are compressed with almost zero-size, while the luminance is just JPEG compressed.

So, an all gray 24-bit image and the corresponding 8-bit grayscale are JPEG compressed almost the same. In natural images, RGB channels slightly differ:

So if you convert them to YUV (luminance-chrominance), you get:

Typically, the luminance one remains almost 8-bit in dynamics, while the two chrominance ones  are flatter, on 4-5 bits. The two latter can be compressed more heavily, as they can be downsampled (they look more low-pass) and amplitude-quantized.

Hence, to the first order, the most of the information is concentrated on the luminance part. So, if the channels do not differ too much, compressing the RGB or the grayscale with JPEG, for sufficiently high rations, is about the same.


Wavelet Fidelity Filter (video compression) I have read a specification of Dirac (video compression format). A key element of this format is the use of the discrete wavelet transform. Section 15.6.3 of the specification defines individual wavelets by making use of a lifting scheme. I am interested in the Fidelity filter for improved down-conversion and anti-aliasing. Where did it come from? Is there a paper that describes it?

For convenience, the lifting filters are defined as:

[ −2, 10, −25, 81, 81, −25, 10, −2 ]
[ −8, 21, −46, 161, 161, −46, 21, −8 ]

Eventually, I found the answer in this white paper. I quote:

  The Fidelity filter was specifically developed for VC-2 using an
  exhaustive search process to select a filter set with simple integer
  filter coefficients, limited to 8-bit resolution, and offering good
  alias rejection in both analysis and synthesis operations.


Measure of intermittency/continuousness of a signal I have three signals (below) each having the same standard deviation, however, are clearly very different temporally.

Is there some such metric that could be calculated for each of these signals to give an indication of how they vary over time, such as a measure of continuousness or intermittency if such a metric exists.

If possible, I would like to do this by calculating a single number for each of these signals, rather than windowing into smaller chunks.

Thanks

Any metric measuring the relative importance of high-frequency content vs low-frequency content will be correlated with the amount of abrupt changes in the signal.

You can try spectral roll-off (the frequency under which x% of the energy is concentrated) or spectral centroid.


Coherence vs. Magnitude Squared Coherence currently I am writing my masther thesis. The theory part is about the turbulent wind field generation, where the coherence (not magnitude squared) is used:

$$\text{coh}(f) = \frac{|P_{xy}|}{\sqrt{P_{xx} P_{yy}}}$$

I am wondering what is the general advantage of using the magnitude squared coherence instead:

$$\text{coh}^2(f) = \frac{|P_{xy}|^2}{P_{xx} P_{yy}}$$

I think this relationship is similar to that of the correlation coefficient $R$ and coefficient of determination $R^2$. I tried to google the difference regarding the coherence, but unfortunately I don't find any (good) explanations.

Hopefully someone can help me. Thanks in advance!

Apollo3zehn

I believe that the calculation of the magnitude squared coherence is faster for computers (multiplication is faster than square root calculations). And that's pretty much it. :P


Separating excitation from filtering for sound sythesis From my days at college I rememeber that LPC splits a signal into a series of slowly changing filter coefficents and an audio-rate residual error signal, such that the original signal gets reproduced when you run the residual signal as "excitation" through a time-varying filter which is controlled by those coefficients.

The cool thing is that you can replace the residual signal by another signal (e.g. a pulse train) and many of the characteristics of the original signal will be preserved. The excitation signal simply does not carry much information (except the pitch) and thus can be replaced.

I find this concept absolutlely intriguing.

I was hoping that this method could be used to sythesize real-world instruments. I could e.g. replace the residual signal of a real Guitar, by one which was created by a Karplus-Strong algorithm. I would have full control over that signal and if I feed it through the filter coefficients taken from a real guitar, this should take aways much of the artificial feel of Karplus-Strong.

The same reasoning led me to believe that one could build a fine vocal pitch correction this way. You pitch corect the residual and most of the unwanted artifacts will get filtered away in the filtering stage.

But all my attempts failed somewhat. For one, there does not seem to be much support for LPC in the open-source community(Luckily "lpanal" is still around). I played with csound and supercollider and I am getting the feeling that I am on the wrong track.

Does the above make any sense or where did I go wrong?
Is this idea of separating excitation from filtering widely used in sound synthesis? (I believe there is some audio compression algorithm, which uses LPC)
Are there other methods which separate excitation from filtering?

Does the above make any sense or where did I go wrong?

Yes.

  Is this idea of separating excitation from filtering widely used in sound synthesis? (I believe there is some audio compression algorithm, which uses LPC)

Yes it does.

For decades it was the strongest approach for speech - and by extension singing voice - synthesis.
As for musical applications, in the framework of subtractive synthesis, filters or resonator banks (Moog 914) were used to model the coloration brought by the body of an instrument, or to create vocal sounds.
The vocoder uses this principle to apply the formants (filters) extracted from a vocal signal onto a synthetic excitation signal coming from a synth.
It is not uncommon to find a filter simulating the body of an instrument as the final step in a physical model (for example a simple bowed instrument model).
Modal synthesis is another sound synthesis technique based on the representation of sound as an excitation shaped through a filter.
One could also consider that the virtual amp/speaker modelling techniques used to build today's electric guitar sound are inspired by this approach - the signal coming straight from the guitar being considered as a raw excitation shaped by the footprint of an amp & speaker.
Actually, any convolution processor (such as a reverberator) might be seen as an application of this technique - when you sing in a room, your voice is the excitation, the room the resonator.
This idea is also relevant to pitch-shifting. When a speaker speaks at a higher or lower pitch, the fundamental is modified, but the formants do not move (ideally...) - so a good pitch-shifting technique should preserve formants. To pitch-shift vocals or musical instruments sounds in a meaningful way, it might be a good idea to pitch-shift only the excitation part of the signal, while keeping the filter representation unchanged. In practice, though, we use techniques that do not explicitly make the actual separation (PSOLA in the time-domain, or spectral envelope manipulations in the frequency-domain when using a phase vocoder).

  Are there other methods which separate excitation from filtering?

Cepstral analysis. In some situations, adaptive filtering.


Single-pole IIR low-pass filter - which is the correct formula for the decay coefficient? A single-pole IIR low-pass filter can be defined in discrete time as y += a * (x - y), where y is the output sample, x is the input sample and a is the decay coefficient.

However, the definition of a varies.
On Wikipedia, it's defined as 2πfc/(2πfc+1) (where fc is the cutoff frequency).

But here, a is defined as follows: 1 - e^-2πfc.

Their graphs look similar, but which one is more accurate?

The blue graph is the Wikipedia formula, the green is the second one, fc is the x-axis.

The given single-pole IIR filter is also called exponentially weighted moving average (EWMA) filter, and it is defined by the following difference equation:

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1],\qquad 0<\alpha<1\tag{1}$$

Its transfer function is

$$H(z)=\frac{\alpha}{1-(1-\alpha)z^{-1}}\tag{2}$$

The exact formula for the required value of $\alpha$ that results in a desired $3\;\textrm{dB}$ cut-off frequency $\omega_c$ was derived in this answer:

$$\alpha=-y+\sqrt{y^2+2y},\qquad y=1-\cos(\omega_c)\tag{3}$$

Even though it should be easy enough to compute $\alpha$ from $(3)$, there are several approximative formulas floating around the internet. One of them is

$$\alpha\approx 1-e^{-\omega_c}\tag{4}$$

In this answer I've explained how this formula is derived (namely, via the impulse invariant transformation of the corresponding continuous-time low pass filter). This answer compares the approximation $(4)$ with the exact formula, and it is shown that $(4)$ is only useful for relatively small cut-off frequencies (of course, small compared to the sampling frequency).

In the wikipedia link given in the question, there is yet another approximative formula:

$$\alpha\approx\frac{\omega_c}{1+\omega_c}\tag{5}$$

[Note that in all formulas of this answer $\omega_c$ is normalized by the sampling frequency.] This approximation is also derived from discretizing the corresponding analog transfer function, this time not via the impulse invariant method, but by replacing the derivative by a backward difference:

$$\frac{dy(t)}{dt}{\huge |}_{t=nT}\approx\frac{y(nT)-y((n-1)T)}{T}\tag{6}$$

This is equivalent to replacing $s$ by $(1-z^{-1})/T$ in the continuous-time transfer function

$$H(s)=\frac{1}{1+s\tau}\tag{7}$$

which results in

$$H(z)=\frac{1}{1+(1-z^{-1})\tau/T}=\frac{\frac{1}{1+\tau/T}}{1-\frac{\tau/T}{1+\tau/T}z^{-1}}\tag{8}$$

Comparing $(8)$ to $(2)$ we see that

$$\alpha=\frac{1}{1+\tau/T}\tag{9}$$

Since the (continuous-time) $3\;\textrm{dB}$ cut-off frequency is $\Omega_c=1/\tau$, we obtain from $(9)$

$$\alpha=\frac{\Omega_cT}{1+\Omega_cT}\tag{10}$$

Equating the discrete-time cut-off frequency $\omega_c$ with $\Omega_cT$ in $(10)$ gives the approximation $(5)$.

The figure below shows the actually achieved cut-off frequency for a given desired cut-off frequency for the two approximations $(4)$ and $(5)$. Clearly, both approximations become useless for larger cut-off frequencies, and I would suggest that approximation $(5)$ is generally useless.


How do I show an impulse response leads to a zero-phase frequency response? I'm trying to understand how to show that with real coefficients, the phase response of a filter is 0. Here is the impulse response

$h[n] = b_1d[n+1] + b_0d[n] + b_1d[n-1]$

How should I approach showing this?

1) Define $d[n]$ as a sinusoidal

2) Apply angle addition formulas

3) Simplify

4) Interpret the results

Hope this isn't too much.   Sounds like a homework problem.

Ced


Understanding the resulting image matrix when differentiating image Let $A$ be a image matrix and let $P(i,j)$ be the gray level of pixel $i,j$. Let $0$ be black and $255$ be white Assume I want to differentiate this image with respect to the columns $(x)$ as in I want $P(i,j)=P(i,j+1)-P(i,j)$ in the new image.

I get that I can achieve this with convolution of $A$ with the filter $f = [1 \ -1]$ since you flip the filter in the $x-$direction (so pixel $i,j$ would be multiplied by $-1$ and $i,j+1$ with $1$ which gives me what I want.

My understanding is that this process is supposed to highlight differences in the $x-$direction.

Two questions:

How does the gray-level scale in the resulting image work? If one pixel had gray level $240$ and the next pixel had gray level $5$, this would result in the new pixel having the gray level $-235$, what does that mean in the scale of $0$ to $255$? Does my scale change from $0...255$ to $-255...255$?
Does it matter if I convolve $A$ with the filter $[-1 \ 1]$ instead of $[1 \ -1]$?. What would be the result?

Your interpretation is correct: directional derivation operators highlight variation in a given direction. Here, you use the $2$-point discrete derivative in the $x$-direction (along image rows). It may emphasize vertical features.

First, such operators indeed extend the initial image range. However, one often uses the absolute value of the derivative to compute a gradient. With the modulus, the derivative values fall back in the initial image range.

Second, as $[1 \;-1] = -[-1\; 1]$, and since convolution is a linear operator, you will get the opposite result as in the previous case. And it won't matter if you take the absolute value.


Analog to digital conversion using Python I'm a beginner. I try emulate analog signal conversion to digital (including sampling by time and quantizing by level) using Python. Here is my code:

import numpy as np
import matplotlib.pyplot as plt

time_of_view        = 1.; # s.
analog_time         = np.linspace (0, time_of_view, 10e5); # s.

sampling_rate       = 20.; # Hz
sampling_period     = 1. / sampling_rate; # s
sample_number       = time_of_view / sampling_period;
sampling_time       = np.linspace (0, time_of_view, sample_number);

carrier_frequency   = 9.;
amplitude           = 1;
phase               = 0;

quantizing_bits     = 4;
quantizing_levels   = 2 ** quantizing_bits / 2;
quantizing_step     = 1. / quantizing_levels;

def analog_signal (time_point):
    return amplitude * np.cos (2 * np.pi * carrier_frequency * time_point + phase);
sampling_signal     = analog_signal (sampling_time);
quantizing_signal   = np.round (sampling_signal / quantizing_step) * quantizing_step;

fig = plt.figure ()
plt.plot (analog_time,   analog_signal (analog_time) );
#plt.stem (sampling_time, sampling_signal);
plt.stem (sampling_time, quantizing_signal, linefmt='r-', markerfmt='rs', basefmt='r-');
plt.title("Analog to digital signal conversion")
plt.xlabel("Time")
plt.ylabel("Amplitude")

plt.show()

In general I satisfied with it: if I play with variables the result is expected for me. Except one thing: according to the listing sampling_rate = 20 and carrier_frequency = 9 since sampling_rate / carrier_frequency > 2 I expect successful conversion but Instead I got distorted signal (look at picture: envelope of the output signal has clearly expressed amplitude modulation). Where is my mistake? Is my implementation correct?

You see this because you're sampling too close to Nyquist rate, but do not worry, although it seems to be aliased however, your sampling system works well and the signal is not aliased. To be more sure, just take a look at sampled signal in frequency domain:
As you see, only 9th bin is active, all others are zero, so no aliasing is occurred in sampling.The code (MATLAB) for FFT calculation is:

Fs = 20;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = length(sampling_signal);  % Length of signal
t = (0:L-1)*T;                % Time vector
y = sampling_signal;     % Sinusoids plus quantization noise
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

Note! I've changed time_of_view to 10 secs, to avoid FFT leakage (the side lobes you see around 9th bin).

It is O.K to increase signal rate this close to Nyquist rate, as far as you're not intended to return the signal back to the analog domain,
 because then your DAC requires a very sharp Low-pass filter in output and its realization is too hard (if not impossible). This paper intuitively expalins sampling problems, Sampling: What Nyquist Didn't Say, and What to Do About It


Trivial echo separation Say we have a signal S, and we composite it with a 0.1-volume 30-sample-delayed version of itself, so:

T(k) = S(k) + 0.1*S(k-30)

How would one rearrange this equation to make S the subject?

i.e given the resultant signal T, and the fact that we know the parameters of the echo (0.1 and 30), how to reproduce the original S?

The way the signal $T(k)$ is generated is by applying an FIR (finite impulse response) filter to the signal $S(k)$. The transfer function of this FIR filter is

$$H(z)=1+0.1z^{-30}$$

If you want to compensate for such a filter, you need a filter with a transfer function which is the inverse of $H(z)$:

$$G(z)=1/H(z)=\frac{1}{1+0.1z^{-30}}$$

This is a recursive filter with an infinitely long impulse response (IIR). If you filter the signal $T(k)$ with this IIR filter, the result will be $S(k)$.


Should I (and if yes, how can I) assure channel-dependant temporal accuracy in a filterbank? Imagine a filterbank with N channels. A possible implementation is a sliding short-time FFT. If you want to capture frequencies as low as 100 Hz, the length of your sliding window should be at least 10 ms. However, because of this long window length (or in general filter length), you cannot track fast temporal fluctuations in the energy at high frequencies.

Is this ever considered as a problem? How would it be adressed? A possible solution would be to different filter lengths for different channels (i.e. shorter filters for longer channels). However,

I guess that in this case you lose the efficiency of the FFT algorithm?
Does it even make sense to consider fast energy fluctuations at higher frequencies, while assuming that the energy in lower frequencies keeps constant?

Thanks in advance!

Using different filter lengths for different channels is basically what the Wavelet Transform does. I think it basically addresses this issue. Wavelet transforms divide the spectrum up into spectral bins of unequal width. There are efficient options for computing the Wavelet Transform (a tree structured filter bank approach usually).

In the linked Wikipedia Article, you can see a picture comparing the time-frequency tiling of the STFT to that of the Wavelet Transform. It is directly applicable to your question.

I realize that I addressed the "How can I" part of your question, but not the "Should I" part... I'm not sure I can answer the "Should I" part unless I am given more information regarding the target application. So, I will briefly say that Wavelets have found widespread use in compression and image recognition. They are often helpful for multiresolution analysis which (in rough terms) is when one wishes to find something of a certain shape but has no idea how large it will be (for example a pulse of a certain shape that may have been compressed or expanded significantly). If you don't require these types of properties, then a uniform filter bank is often the way to go.


Need help in understanding minimum variance estimator and CRLB concept using an example I am reading the book, Fundamentals of Statistical Signal processing, Estimation Theory Volume 1 by Steven M. Kay. In Chapter 2, there  is a Fig 2.5 which illustrates which estimator to select based on the CRLB. Basically, the estimator with less variance is preferred since its pdf is more concentrated about the true value. I have certain doubts regarding this statement and the term "inconsistent estimator".

Considering 2 estimators for a parameter $\theta$ found using two estimation technique. Estimation technique A is the ML estimator and B is the ordinary least square having the CRLB denoted as CRLB_A and CRLB_B for the unknown parameter $\theta$ respectively.

Doubt 1: At high SNR, the CRLB_A > CRLB_B. So, should I select CRLB_B at high SNR? For low SNR, it is the reverse.

Doubt 2: Let the $mse(A) = E[(\Delta \theta \Delta \theta^T)] < CRLB_A$

What does this imply? Does it mean that the ML estimator of A is inconsistent?

You may want an estimator with minimun variance because each time you calculate it, it is likely to be closer to its expectation, and in the case that the estimator is unbiased, closer to the true parameter. An estimator is inconsisten if increasing the number of samples of the dataset to calculate it doesn't reduces its variance.

The sample mean estimator of X:
$$
\bar X = \frac{1}{N}\sum_{i=0}^{N-1}x_i
$$
is unbiased and its variance is:
$$
\mathbb var(X)=\frac{\sigma^2}{N}
$$
It is clear that as N increases, the variance gets smaller, so it is consistent.

Regarding question 1, NO. For instance, one of the two estimators can have a higher CRLB always, for both high or low SNR. The criteria for choosing an estimator depends on the problem, the prior knowledge about it, the computation constraints and so on. Normally it is desired a minumn variance and a non-biased estimator, i.e. an efficient estimator. But, if your decision is only based on the CRLB, then choose that one with the smaller one.

Take this 'silly' example: for a dataset you decide that the estimator for the mean is always (arbitrarily) 3.5 even if you know nothing about the data and you do no calculations... It is clear that this estimator has variance 0, because it will be always 3.5, but it is strongly biased, it is randomly biased indeed... is it a good estimator? No, although it has 0 variance. The key point here is that yo can't have an unbiased estimator with a lower variance than its CRLB (if the CRLB exists). If the variance of your estimator is lower than the CRLB, it is biased.

If you have some experience with computer simulation, Matlab or similars, it is a good exercise running Montecarlo simulations and plotting the variance and the CRLB of each realisation versus the SNR for a fixed number of samples, or/and for a fixed SNR plotting the variance and the CRLB of the estimator versus different dataset lengths.


Why are there so many windowing functions? Many windowing functions are listed here in the Mathematica documentation. I tried using a few to reduce leakage when computing a Discrete Fourier Transform. From what I could tell it made little difference which windowing function is used. Two of them that are very dissimilar are the BartlettHannWindow and BlackmanHarrisWindow. Could someone provide an example where the BarttletHannWindow is a very good choice, and another where the BlackmanHarrisWindow is a very good choice. That might help me understand why there are so many to choose from.

Aside from reduction of spectral leakage, there is a one major trade-off to be made when choosing a window function. Below you can see a figure with various parameters. Two of them are most important:

Main Lobe Width
Peak Side Lobe Level

Width of the main lobe affects resolution of your analysis. I am sure you know that multiplication in time domain is convolution in frequency domain. By default you always use rectangular window, which has frequency response in a form of $\mathrm{sinc}$ function.

For example it might be the case that two pure sinusoids (Fourier Transform being Dirac deltas) with frequencies being very close to each other, peaks can be very hard to distinguish and thus become one. On figure below I compared 4 window functions for two sinusoids separated exactly by 3 frequency bins. You can see how wide is the main lobe of Nuttall window and how it affects the separation, whereas Rectangular window is the window with narrowest possible width of main lobe (i.e. lowest Normalized Equivalent Noise Bandwidth - NENBW)

On the other hand level of first side lobe affects the ability to recognize adjacent frequencies when spectral leakage is present and side lobes of frequency with high amplitude are almost "covering" other frequency next to it with much lower amplitude:

And here is the trade-off... You can't get very low level of side lobes and very narrow main lobe, i.e. it's not possible to have window function with main peak like in rectangular window, and very low level of side lobes as in Nuttal window. Here is the comparison of some windows:

Apart from parameters mentioned above, there is more features such as:

Side Lobe Roll-off rate (how fast side lobes decrease - for example see Bartlett and Nuttal above)
Existence of analytical frequency response equation (useful when performing theoretical calculations)
Summing of window samples to constant value in time domain when doing overlap (for example Hann window with $33\%$ of overlap)

This is extremely useful when you want to perform STFT and then reconstruct time domain signal.

So now to answer your question. Why do we have so many window functions?

Depending on task whether it is simple FFT analysis of: closely spaced sinusoids, white noise, impulsive sounds, or even FIR filter design, different type of window is recommended. Obviously in 90% of cases people are using Hamming/Hanning window or no window at all. Nevertheless for some specific signals you might find that extracted features depending on window choice tend to be more discriminative.

For more info about Window Functions please refer to:

  F. J. Harris - On the Use of windows for Harmonic Analysis with the
  Discrete Fourier
  Transform

  G. Heinzel, et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows

  M. Cerna, A. F. Harvey - The Fundamentals of FFT-Based Signal
  Analysis and
  Measurement

  H. A. Gaberson - A Comprehensive Windows
  Tutorial


Rotated object detection by Viola jones Algorithm In my project I need to detect rotated and non-rotated faces in real time. I read original articles and try to use common implementation in matlab and openCV. And as I understood Haar features are not rotation invariant. What should I do for successes detection of rotated faces ? I have some ideas about it. Which is more 'real'?

Try to rotate original image with some angle step
Train own classifier on rotated examples
Use any modified Viola Jones Algorithm from research articles
Use another algorithm (like SURF, KLT, Skin based face detection)

You are correct, Haar features, as well as LBP and HOG are not rotation invariant. I have experimented with vision.CascadeObjectDetector in MATLAB's Computer Vision System Toolbox, and found that the face detector model it comes with can tolerate about 15 degrees of in-plane rotation.

So if you want to handle in-plane rotation, the easiest thing is to rotate the image by 15 degree increments and run the detector each time.  If you are talking about out-of-plane rotation, then you would need multiple detectors, such as the frontal and profile detectors that come with MATLAB.

There is another possibility. If you need to detect the faces in video, then you can handle rotations and even partial occlusions by tracking them using KLT.  See this example.


Reconstructed output mismatch for LTI system I have a system with measured input (u) and output (y).

I assume that this is an linear time-invariant (LTI) system and I want to find its impulse response (ir).

This is the Matlab code I use to find ir:

ir = real(ifft(fft(y)./fft(u)));

And this is how I reconstruct the system output:

Y = conv(ir, u);

However, as can be seen below, the reconstructed output differs greatly from the true output. What am I missing here?

The full matlab code including data can be found here

The problem is that the formula $$H(f) = \frac{Y(f)}{X(f)},$$ where $H(f)$ is the frequency response, and $X(f)$, $Y(f)$ are the input and output, is valid only for frquencies $f$ where $|X(f)| \neq 0$.

Furthermore, if $|X(f_0)| \approx 0$, you're going to run into numerical errors, since $|H(f_0)| \rightarrow \infty$.

Your input signal has a bunch of places where it's zero, or close to zero. I wrote a little function that takes care of those cases:

function h = impres(x,y)

X = fft(x);
Y = fft(y);

for i = 1:length(X)
    if abs(X(i)) > 10
        H(i) = Y(i)/X(i);
    else
        H(i) = 0;
    end
end

h = ifft(H);

end

While not perfect, and could be tweaked, I get:

where the red line is the estimated output.

One more thing to be aware of: the equivalent in time of the IDFT is the circular convolution, not the linear convolution. In this case there is not too much difference, but you want to be as correct as you can in your code.