%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')