Basics of leakage phenomena in DFT and its approximation with sinc function

Question

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?

Answer 1

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:

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))

Answer 2

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).