Why is spectrum.periodogram not recommended, and how to substitute pwelch in it's place?

A message pops up in Matlab when I use spectrum.periodogram to find the median frequency of a signal, saying that it is not recommended. Why is this? I also get the same message for spectrum.pwelch

I am implementing the solution described by Median Frequency using spectrum.periodogram which is:

psdest = psd(spectrum.periodogram,x,'Fs',1000,'NFFT',length(x)); 
   normcumsumpsd = cumsum(psdest.Data)./sum(psdest.Data);
   Ind = find(normcumsumpsd <=0.5,1,'last');
   fprintf('Median frequency is %2.3f Hz\n',psdest.Frequencies(Ind));

After lots of research I still don't understand the output of the psd well enough so that I can susbstitute spectrum.periodgram. Normally I find my PSD using:

[Pxx,Fx] = pwelch(s,[],[],[],fs);
plot(Fx,10*log10(Pxx))

Ideally I would like to find the median frequency by manipulating the Pxx and Fx values but I am struggling to relate Pxx directly to psdest.Data. I would be grateful if someone could point me in the right direction. Thank you in advance!

Answer 1

[Daniel kiracofe]

Well, first, is your signal a random variable? Or more of a periodic signal (i.e. a sine wave)? If it's periodic, then just a simple FFT will suffice to give you the frequency. You only need to mess with power spectral density if you have a random signal.

As to why periodogram is not recommended... first, let's establish one fact: you can never actual measure power spectral density, because to do that you'd need an infinitely long sample of the data. You can only estimate power spectral density with a finite length sample. And, as it turns out, the periodogram is not a very good estimate. One is problem is that as you take longer and longer data samples, you would expect to get better and better estimates. But this does not happen with periodogram. The estimate remains noisy even for more and more data.

welch's method (pwelch()) is better because as you take more and more data, the estimate gets better. There are other methods, but pwelch is pretty reasonable, and I use it a lot.

Hope this helps. I've also got some tutorials on power spectral density on my website: http://www.mechanicalvibration.com/Random_vibrations.html

[Mary]

Thanks for the quick reply. I have been using pwelch (signal is EMG, so random) for the psd estimate. Thanks for clarifying re periodogram.

My problem now is trying to get the median frequency, even when I use spectrum.welch in the the psd (see above) the message about not being recommended pops up, and says to use pwelch instead.

Is there an easy way to use the output of pwelch (Pxx and Fx in the code above) to calculate the median frequency? I tried substituting psdest.Data with Pxx, but that didn't work and I'm not sure what the format of the psdest.Data is, and haven't been able to figure that out from the matlab help.

[Daniel kiracofe]

I don't get any warning message in version 2011. What version are you using? The documentation in version 2011 does say that the "spectrum" object is obsolete and has been replaced by other functions, such as pwelch(). I would suspect that that warning message is just to let you know that the spectrum object might be removed from future versions of matlab, so if you want to write code that will work on the next version, you should avoid using it. Or at least that is my guess.

psd using spectrum.welch and pwelch() should use the exact same algorithm, although some of the defaults may be different. I was under the impression that psdest.Data and Pxx in your code above should both represent the same thing. I would have thought that normcumsumpsd = cumsum(Pxx)./sum(PXxx) would have worked. Can you elaborate on why that didn't work for you?

[Mary]

I was using version 2014. The message pops up if you hover near the code with the mouse in the editor window.

I got the code to work in the end. Just a silly mistake in translation. Thanks!

[krn99]

Mary can you please send the code for PSD using pwelch and to find Median and Mean Frequency

[Greg Dionne]

@krn99: see medfreq and meanfreq in the Signal Processing Toolbox.

Answer 2

[Sahaj Sandhu]

Hey, have you calculated other features also ? like mean power or total power ?

Answer 3

[bhavya kailkhura]

Is there an implementation of pwelch for 2d data? For example, if I want to plot psd of an image with dc component centered, how can I use pwelch to do that?

Thanks!

[balandong]

Calculate at each dimension then sum the two

Answer 4

[Greg Dionne]

If you have R2015a or later, try medfreq.

The spectrum package is no longer recommended for use.