%fft and pwelch方法求取功率谱
load x.mat
Fs = 1;
t = (0:1/Fs:1-1/Fs).';
Nx = length(x);
% Window data
w = hanning(Nx);
xw = x.*w;
% Calculate power
nfft = Nx;
X = fft(xw,nfft);
mx = abs(X).^2;
% Normalize by window power. Multiply by 2 (except DC & Nyquist)
% to calculate one-sided spectrum. Divide by Fs to calculate
% spectral density.
mx = mx/(w'*w);
NumUniquePts = nfft/2+1;
mx = mx(1:NumUniquePts);
mx(2:end-1) = mx(2:end-1)*2;
Pxx1 = mx/Fs;
Fx1 = (0:NumUniquePts-1)*Fs/nfft;
[Pxx2,Fx2] = pwelch(x,[],[],[],Fs);
plot(1./Fx1,Pxx1,1./Fx2,Pxx2,'r:');
legend('PSD via FFT','PSD v ia pwelch')