8

I am playing with the periodogram of MATLAB. I created a simple script to observe how it behaves:

rng(1);  %# initialize the random number generator

Fs = 1000;  %# Sampling frequency
duration = 0.1; %# seconds

A = 1; %# Sinusoid amplitude
f = 150; %# Sinusoid frequency
eps = 0.01;

t = 0:1/Fs:duration;
x = A * sin(2*pi*f*t) + eps * randn(size(t));

periodogram(x,[],1024,Fs);

enter image description here

I have no problem with the code and can write my own periodogram function using the algorithms given in the documentation but I wonder the theoretical reason behind the comb-like hills which are not 150 Hz. What do I get those instead of getting a single spike over 150 Hz? Is there anything special in the distances of the peaks of these hills?

Lorem Ipsum
  • 5,944
  • 3
  • 34
  • 42
petrichor
  • 235
  • 1
  • 3

2 Answers2

13

A single spike (as you call it) appears theoretically only for a infinite-length sinusoid. Since your signal is 100 samples length, it is not infinite. You actually multiplied your infinite signal with a window which has a value of 1 over 100 samples, and 0 elsewhere. Since multiplication in time domain is equivalent to convolution in the frequency domain, your spectrum is a convolution of the single spike and the frequency response of the window (btw it is called rectangular window). This is the function you got.

I suggest you read about windowing: http://en.wikipedia.org/wiki/Window_function

Itamar Katz
  • 239
  • 1
  • 3
  • +1 Oh, I knew about the windowing but couldn't make the link. Thanks! – petrichor Nov 24 '11 at 14:54
  • 1
    A single spike appears regardless of window used, if the frequency is an exact fit for the window length. https://gist.github.com/236567 – endolith Nov 28 '11 at 17:16
  • 1
    That is not correct. For a rectangular window this is true, since you sample the window function in the frequency domain exactly in its zeros, so you are 'blind' to the side lobes. However it is not true for a general window function. – Itamar Katz Nov 29 '11 at 07:01
  • see example git://gist.github.com/1403819.git – Itamar Katz Nov 29 '11 at 07:10
  • @ItamarKatz: Yeah you're right. I meant "without a window". – endolith Dec 01 '11 at 21:04
5

I'm not entirely satisfied by Itamar Katz's answer, so here's my explanation.

The DFT of an $N$ length complex signal, $x[n]=e^{\imath 2\pi f n/N}$ is

$$X[k]=\mathcal{F}\{x[n]\}=\frac{e^{\imath 2\pi (f-k)}-1}{e^{\imath 2\pi (f-k)/N}-1}$$

So, the power or the magnitude squared response is given by

$$\left\vert X[k]\right\vert^2 = \left(\frac{\sin\left(\pi(f-k)\right)}{\sin\left(\pi(f-k)/N\right)}\right)^2$$

As you can see, the above expression is zero whenever $f-k$ is an integer. You can convince yourself that the denominator is zero at only one point, and at this point, taking limits gives you a value $N^2$ for the ratio. Hence, there is no point at which the expression blows up.

Now when you take the log of the above expression, $log_{10}(0)$ is $-\infty$ (or for that matter, in any base) and hence you get nulls everywhere you had a zero. This is what results in the "comb like hills" in your plot.

Here's a short illustration in Mathematica:

Clear@X
X[f_, n_] := (Sin[π (f - #)]/Sin[π (f - #)/n])^2 &
Plot[X[3, 10][k], {k, -5, 5}, PlotRange -> All]

enter image description here

Frequency is on the x-axis and power (linear) is on the y-axis. You can see that the zeros occur at integer values and the peak is at 3, which is the frequency I had chosen. Now taking $log_{10}$ of the above, you get nulls which give rise to the comb like structure

enter image description here

Here's another example with a larger $N$, showing more nulls.

enter image description here

Lorem Ipsum
  • 5,944
  • 3
  • 34
  • 42