From: Ryan on
All,

I’m doing some signal processing research which requires the use of the periodogram. I coded a periodogram script on my own and then compared the results to the MATLAB’s periodogram.m. I get very good (almost perfect) agreement between my code and the syntax

NFFT=262144;
[Pxx,w] = periodogram(x,[],NFFT);

Where x is a 266655 length signal vector. However if I use the syntax

periodogram(x,[],NFFT)

what's plotted is not scaled the same as log(Pxx). The plot is the correct shape but the location on the Y axis is off by about 35db. This 35db shift is very important to my subsequent calculations and I would like to know which approach is correct. I'm wondering if anyone can explain where or what this scale factor is?!

The help file has this to say about the different syntax:


"[Pxx,w] = periodogram(x) returns the power spectral density (PSD) estimate Pxx of the sequence x using a periodogram. The power spectral density is calculated in units of power per radians per sample.

periodogram(...) with no outputs plots the power spectral density in dB per unit frequency in the current figure window."

I'm not sure how to interpret "units of power per radians per sample" vs. "dB per unit frequency."

I've tracked periodogram.m through computeperiodogram.m to computepsd.m to computeDFT.m to fft.m and it made perfect sense. I think this is why my method and [Pxx,w] = periodogram(x,[],NFFT); compare well. However, when you call periodogram(x,[],NFFT) lines 141-144 of periodogram.m call

hspec = spectrum.periodogram({winName,winParam});
hpsd.Metadata.setsourcespectrum(hspec);

plot(hpsd);

This code produces the plot with "improper" scaling. I don’t know how to track this syntax.

For reference, "my method" of producing the periodogram is as follow

dfftx=fft(x,NFFT)/sqrt(NFFT); %take fft and scale by sqrt(N)

fhatxx=dfftx.'*conj(dfftx); %mult by complex conj

fhatxx=2*fhatxx; % mult by two to look at the half spectrum

fhatxx=fhatxx/2/pi; %scale by 2*pi since no sampling freq was presented

figure
plot(log(fhatxx))

I know I got a little sloppy with multiplying all the terms by 2 and not excluding frequencies at zero and pi but let's just ignore those. As I said, the above algorithm agrees very well with the result of [Pxx,w] = periodogram(x,[],NFFT); but not well at all with periodogram(x,[],NFFT). If someone could explain this to me I would be very grateful.


Thanks,
Ryan