0
$\begingroup$

I'm learning some basic signal processing to use in my research. So I'm trying to understand the pwelch function in MATLAB to calculate the PSD. My understanding was that the PSD is the density of power in the frequency bins, and thus the PSD of a sine wave should be a delta function. However, when I run my simple code I get a distribution.

enter image description here

  • Why is this happening?
  • Is this something to do with windowing?
  • Or is my understanding of PSD wrong?

EDIT:

Here is my MATLAB code:

rng('default');
fs = 1
n = 0:1/fs:200000;
x = sin(2*pi*n);
pwelch(x);
$\endgroup$
  • $\begingroup$ Are you sure you're using pwelch correctly? I just tried and I get a single peak at the expected frequency. $\endgroup$ – MBaz Jun 6 '16 at 18:03
  • $\begingroup$ Hi, @Mbaz. I have added my matlab code to the original post. Am I doing something wrong? $\endgroup$ – elenasto Jun 6 '16 at 18:20
  • $\begingroup$ Well, your code changes your question significantly. $\endgroup$ – Marcus Müller Jun 6 '16 at 18:55
  • $\begingroup$ Now, your y-axis makes sense. You're basically looking at values that are nearly all zero, and all that you see are numerical artefacts $\endgroup$ – Marcus Müller Jun 6 '16 at 19:07
3
$\begingroup$

The signal

sin(2*pi*n)

For $n$ integer will give you zero, or something around machine epsilon. Try setting your signal to

sin(2*pi*n*0.01)

And you should get more like what you expect.

$\endgroup$
  • 1
    $\begingroup$ @ Peter Thanks the worked. But I'm not quite sure why the former didn't work. Whats wrong with getting zero at integer n? $\endgroup$ – elenasto Jun 6 '16 at 18:53
  • $\begingroup$ @SharanBanagiri what's the value of $\sin(2\pi n) \,\forall n\in\mathbb N$? $\endgroup$ – Marcus Müller Jun 6 '16 at 18:55
  • $\begingroup$ Zero. But I'm not quite sure why that would matter. I'm sorry if I'm being a little thick. $\endgroup$ – elenasto Jun 6 '16 at 18:56
  • $\begingroup$ the spectrum of a constantly zero signal must be zero. Since it is only damn near zero, all operations are extremely close to how exact your machine can calculate, and hence all the effects you observe are just "waste". $\endgroup$ – Marcus Müller Jun 6 '16 at 19:09
  • $\begingroup$ @SharanBanagiri Look up Nyquist frequency. This tells you that in this case you can only find the power of frequencies lower than half the sampling frequency. Since your sampling frequency is 1, then the frequency of the signal can be at most a half. $\endgroup$ – fibonatic Jun 6 '16 at 20:15
1
$\begingroup$

The PSD of a sine observerved for infinite time is indeed a dirac impulse.

Is this something to do with windowing?

But: you're not observing the actual PSD of a sine, but you're observing a small excerpt of time. So, yes, this is due to windowing.

Think about it: because you don't do anything with an infinite duration of observation, it's like you take a rectangular window, and use it to cut out a piece of infinity (your observation length). Now, that's multiplication with a rectangle in time domain, which is equivalent to convolution with a sinc in frequency domain.

You can get a very sharp PSD if you just apply the DFT to a sine signal actually sampled exactly for a multiple of its period, but that's not the case in general, either.

Or is my understanding of PSD wrong?

That's hard to tell, because I don't know your understanding of Power Spectral Density!

Being strict, the PSD is actually a property of a random process, namely:

The power spectral density is the Fourier Transform of the autocorrelation function of a wide-sense stationary random process.

As you can see, that's a property of a random process that you'll never be able to observe directly, but can only deduct from other observations, given some assumptions on the process observed.

Now, the "technical" understanding of PSD is "I went, and looked into infinitesimally small parts of the spectrum, and noted down the Energy in there, after observing those for a finite amount of time".

What's implicit in there is that by observing for a finite time, you can make assumptions on the expectation value of the magnitude square of the Fourier Transform of the observed signal (which is a random process); by the Theorem of Wiener-Chintschin, this is the same as the Fourier Transform of the autocorrelation function – but that, again, makes the assumptions

  • that your random process is wide-sense stationary (which it isn't if your observation window is not a multiple of a period) and
  • that your observation's is an efficient estimator of the signals properties.
$\endgroup$
  • $\begingroup$ @ Marcus Thank you for your reply. What you said about the sampling for a multiple of period makes sense. I added my code above and it goes to a multiple of a period. However what I found after playing with it a little is that it for some reason depends on the sampling frequency. Whatever he length of my time interval may be, for fs<2 the sharp peak goes away and I get a graph like above. But for fs>2, I get a peak. Do you know why that is? $\endgroup$ – elenasto Jun 6 '16 at 18:48
  • $\begingroup$ I don't really understand your question, but this sounds like a Shannon-Nyquist problem. $\endgroup$ – Marcus Müller Jun 6 '16 at 18:56
  • 1
    $\begingroup$ I just realized that we wasted our time. with fs == 1 you just produce an array of values that are very close to zero, and hence numerically unstable, before trying to estimate the spectrum of that. $\endgroup$ – Marcus Müller Jun 6 '16 at 19:08
  • $\begingroup$ Yes, I get it now. It was an issue with the sampling freq being so small that the it wasn't able to sample the sine wave properly. The problem goes away with a higher sampling freq. Thanks a lot for your help :) $\endgroup$ – elenasto Jun 6 '16 at 19:09

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.