psd via fft and pwelch

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