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$
4
  • $\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, 2016 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, 2016 at 18:20
  • $\begingroup$ Well, your code changes your question significantly. $\endgroup$ Jun 6, 2016 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$ Jun 6, 2016 at 19:07

3 Answers 3

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$
5
  • 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, 2016 at 18:53
  • $\begingroup$ @SharanBanagiri what's the value of $\sin(2\pi n) \,\forall n\in\mathbb N$? $\endgroup$ Jun 6, 2016 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, 2016 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$ Jun 6, 2016 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, 2016 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$
4
  • $\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, 2016 at 18:48
  • $\begingroup$ I don't really understand your question, but this sounds like a Shannon-Nyquist problem. $\endgroup$ Jun 6, 2016 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$ Jun 6, 2016 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, 2016 at 19:09
0
$\begingroup$

Looking at the power density spectrum you get, it seems evident that you are seen a too wide spectrum portion. Nyquist theorem says that after n*1/2 sampling frequency, spectrum is repeated (n=1,2,3...) and centered on multiple sampling frequency value, to infinite. The spectrum shown says that your sinus is summed to a constant and that sampling is about 0.54 (the constant comes from lobe at 0 Hz). Your graph repeats itself after about 0.54/2=0.27. Delta function spectrum for a sinus is a theoretical outcome for an infinite signal. This is not the case. It is multiplied with a window in time, that generates a convolution of relative spectra.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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