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