7
$\begingroup$

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?

$\endgroup$

migrated from stackoverflow.com Nov 24 '11 at 21:07

This question came from our site for professional and enthusiast programmers.

5
$\begingroup$

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

$\endgroup$
12
$\begingroup$

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

$\endgroup$
  • $\begingroup$ +1 Oh, I knew about the windowing but couldn't make the link. Thanks! $\endgroup$ – petrichor Nov 24 '11 at 14:54
  • 1
    $\begingroup$ A single spike appears regardless of window used, if the frequency is an exact fit for the window length. gist.github.com/236567 $\endgroup$ – endolith Nov 28 '11 at 17:16
  • 1
    $\begingroup$ 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. $\endgroup$ – Itamar Katz Nov 29 '11 at 7:01
  • $\begingroup$ see example git://gist.github.com/1403819.git $\endgroup$ – Itamar Katz Nov 29 '11 at 7:10
  • $\begingroup$ @ItamarKatz: Yeah you're right. I meant "without a window". $\endgroup$ – endolith Dec 1 '11 at 21:04

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.