3
$\begingroup$

I am trying to understand leakage phenomena when using DFT. In the following I am reffering to Lyons: Understanding Digital Signal Processing (2nd ed.), p.69-71.

Having real valued cosine with $N$ samples as an input, the amplitude response of N-point DFT can be approximated by (in terms of bin index $m$)

$|Y(m)|\approx\frac{N}{2}\operatorname{sinc}(k-m)$,

where $k$ is number of cycles inside $N$ samples.

If $k$ is a whole number, there is only 1 bin non-zero. Otherwise spectral leakage occurs. So far so good.

Now I am want to replicate this in Matlab using following code.

N = 32;                % No of samples
fs = 32000;            % sampling frequency
m = linspace(0,N-1,N);
k = 8.5;               % No of cycles within N samples
X = sin(2*pi*k/N*m);   % cosine signal
Y = fft(X);
f = linspace(0,fs*(N-1)/N,N);
stem(f,abs(Y));hold on

With $k$ being a whole number, no leakage is obtained.

However, using $k=8.5$ produces peaks at 8 kHz and 9 kHz which differ in magnitude, althought slightly ($|Y((8 \text{kHz})/f_\text{s})|=10.15$ and $|Y((9 \text{kHz})/f_\text{s})|=10.25$ ). In my reference book number 10.17 is given.

My questions are:

  1. Shouldn't the two peaks be exactly the same? Why they are not?

  2. Additionaly, I have plotted function $|Y(m)|$ with ''continious'' $m$. The bins of original spectra should be exactly on the resulting sinc function. However, farther away we go from 8.5 kHz, bigger the discrepancy. Aren't the DFT bins just ''uncovered'' parts of underlying continious spectra?

$\endgroup$
  • 1
    $\begingroup$ the answer to Q1 is that there is a circular nature to the DFT and that a real sinusoid at, say 9 kHz, has both a positive and negative frequency components at $\pm$9 kHz. some of that at 9 kHz is leaking from -9 kHz, which due to the nature of the DFT is at some positive-indexed bin above N/2 . $\endgroup$ – robert bristow-johnson Feb 27 '17 at 21:30
3
$\begingroup$

No, the two peaks don't need to be the same. The discrete-time Fourier transform (DTFT) (samples of which are computed by the DFT) is also not symmetric with respect to the sinusoid's frequency (as explained in Robert's comment). If you compute the DTFT of a truncated sinusoid

$$x[n]=\sin(\omega_0n)(u[n]-u[n-N])$$

you get

$$X(e^{j\omega})=\frac{1}{2j}\left(e^{-j(\omega-\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega-\omega_0)N}{2}\right)}{\sin\left(\frac{\omega-\omega_0}{2}\right)}-e^{-j(\omega+\omega_0)\frac{N-1}{2}}\frac{\sin\left(\frac{(\omega+\omega_0)N}{2}\right)}{\sin\left(\frac{\omega+\omega_0}{2}\right)}\right)$$

The first term in the parentheses has its main contribution around $\omega=\omega_0$, whereas the second term has it around $\omega=-\omega_0$. Note that the two functions are no sinc functions; only for $\omega$ very close to $\omega_0$ (or $-\omega_0$) can they be approximated by a sinc function. And that still wouldn't take into account their overlap.

So your results are correct, and the plot below shows your DFT on top of the DTFT as computed above:

enter image description here

Here's the Matlab/Octave code for the DTFT plot:

N = 32;
k = 8.5;
fs = 32000;
w = 2 * pi * linspace(0,1,500);
w0 = 2*pi*k/N;
X = exp(-1i*(w-w0)*(N-1)/2) .* sin(N/2*(w-w0))./sin((w-w0)/2);
X = X - exp(-1i*(w+w0)*(N-1)/2) .* sin(N/2*(w+w0))./sin((w+w0)/2);
X = X/(2*1i);
plot(w*fs/(2*pi),abs(X))
$\endgroup$
  • $\begingroup$ That explains why sinc function is merely an approximation. Thank you. I tried to replicate DTFT plot in above figure but to no avail. Can you provide me with the code? $\endgroup$ – KlemenD Feb 28 '17 at 8:38
  • $\begingroup$ @KlemenD: I've added the code. You just need to implement the given formula. $\endgroup$ – Matt L. Feb 28 '17 at 8:47
  • $\begingroup$ I was setting w = 2*pi * linspace(0,fs,500) and w0 = 2*pi*k/N*fs which gave me wrong result. I assume the DTFT formula expects w and w0 being relative frequencies? $\endgroup$ – KlemenD Feb 28 '17 at 9:33
  • $\begingroup$ @KlemenD: Yes, that's usually what is meant by writing $X(e^{j\omega})$. You also sometimes see $X(e^{j\omega T})$, where $T$ is the sampling interval. In that case, $\omega$ is not normalized. $\endgroup$ – Matt L. Feb 28 '17 at 9:38
1
$\begingroup$

For sampled strictly real data, the Sinc function will be repeated periodically across an infinite spectrum (this is what allows undersampling to work) (sum called a Dirichlet kernel), and that infinite series of Sinc functions will be conjugate mirrored (needed to cancel out any imaginary components in the real data).

The "interference" from conjugate mirrored Sinc (with the shape of the assumed symmetric Sinc) is especially large for frequency peaks near DC and Fs/2, where higher portions of the 2 Sincs overlap. (Thus, the conjugate mirror interference will be smaller near Fs/4, use that frequency if you want to see symmetry).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.