periodogram function
EUIX Wang
"periodogram uses an nfft-point FFT of the windowed data (x.*window) to compute the periodogram. If the value you specify for nfft is less than the length of x, then x.*window is wrapped modulo nfft. If the value you specify for nfft is greater than the length of x, then x.*window is zero-padded to compute the FFT."
Now what exactly does the first part mean? I understand how zero padding works what is wrapped modulo?
Also does this work if your value of nfft is not a integer multiple of your number of data points? For example if I have 1000 data points and I pick nfft=210. Then does matlab zero pad the data with 50 extra zeros when executing periodogram?
Wayne King
Hi, what it means is that it divides the data into nonoverlapping segments of length NFFT, possibly padding the last one with zeros if a case like you describe holds. Then it stacks those segments and adds them together to form a length NFFT segment. So for example:
X=1:20;
NFFT=10;
Y = buffer(X,NFFT);
Y=sum(Y,2);
Hope that helps,
Wayne
TideMan
> "EUIX Wang" <kuh...@gmail.com> wrote in message <hnbksd$jc...@fred.mathworks.com>...
Wayne
Are you sure about this?
I thought it transferred all the segments to the freq domain (using
FFT), squared (to get PSD), then took the mean.
If you do this in the time domain, as you suggest, you could
conceivably cancel the signal before it gets transformed.
Wayne King
Hi Tideman, I'm pretty sure it does it that way, but you're definitely making me think, which is always good :) Anyway, here's what I think it's doing with an example
reset(RandStream.getDefaultStream);
x = randn(10,1);
nfft = 5;
win = hamming(10);
y = x.*win;
%Periodogram with nfft less than segment length
[Pxx,Fxx] = periodogram(x,hamming(10),nfft);
% Pxx is 0.0024 0.4686 0.4955
% now the brute force way
% wrapping the windowed time data with datawrap()
% same as sum(buffer(y,5),2)
Y = datawrap(y,nfft);
ydft = fft(Y,5);
% Window normalization
U = win'*win;
ydft = (1/(2*pi))*ydft.*conj(ydft)/U;
% select the bins for one-sided estimate
ydft = ydft(1:(nfft+1)/2);
% scale for one-sided estimate 2*everything but DC
ydft(2:end) = 2*ydft(2:end);
Wayne
TideMan
> TideMan <mul...@gmail.com> wrote in message <3a040376-fd14-418e-9574-31ae51d0a...@q2g2000pre.googlegroups.com>...
Yes, this is what my personal fspect function does as well.
Except that I cheat when it comes to scaling and simply set the area
under the spectrum to the variance.
EUIX Wang
Hey thanks for the reply but I think I am still missing something. Tell me if I am understanding this correctly. If for example I have a dataset of 1000 points and using nfft=210. Is matlab using datawrap on the 1000 points, thus creating 210 points and then applying the periodogram function to that 210 point dataset (that is periodogram(data210,[],210,Fs))?
Because I tried that and compared it to taking the periodogram directly (periodogram(data1000,[],210,Fs)) and while the shape seems to be the same the two data sets seem to differ by some factor. I am posting this pretty late but perhaps it has to do some with the sampling frequency Fs?
EUIX Wang
Wayne King
> I did a little bit of tinker and I think the factor is N/nfft, where N is the total number of points. I say I think because the exact decimal places don't always match up. Although I am not sure if that is because of precision issues in matlab. However I am not where to come up with this extra factor just yet. It's getting late so I will think about it tomorrow.
Hi, It sounds like it's definitely just the scaling factor from your description. In your code are you using 1/Fs as part of the scaling factor? I think you'll find that the scaling is
1/(Fs*N) where N is the length of the input series (before wrapping). Then, for the one-sided PSD estimate, you scale everything except DC by 2. For example:
reset(RandStream.getDefaultStream);
x=randn(10,1);
N = length(x);
nfft=5;
Fs=1000;
%Periodogram with nfft less than segment length
[Pxx,Fxx]=periodogram(x,[],nfft,Fs);
% Pxx is 0.0039 0.0061 0.0077
% now the brute force way
% wrapping the windowed time data with datawrap()
% same as sum(buffer(x,5),2)
X=datawrap(x,nfft);
xdft=fft(X,5);
Pwr = xdft.*conj(xdft);
% Scaling
Pwr = (1/(Fs*N))*Pwr;
% select the bins for one-sided estimate
Pwr=Pwr(1:(nfft+1)/2);
% scale for one-sided estimate 2*everything but DC
Pwr(2:end)=2*Pwr(2:end);
% compare Pwr and Pxx
Hope that helps,
Wayne