Power Spectral Density Estimates Using FFT---MATLAB

% Power Spectral Density Estimates Using FFT
% pr8_2_1 
clear all; clc; close all;

randn('state',0);                    %  Random number initialization 
Fs = 1000;                           %  sampling frequency 
t = 0:1/Fs:1-1/Fs;                   %  Time scale 
f1=50; f2=120;                       %  Frequency of two sinusoidal components 
x=cos(2*pi*f1*t)+3*cos(2*pi*f2*t)+randn(size(t)); %  The signal 
%  Use FFT
N = length(x);                       % x length 
xdft = fft(x);                       % FFT
xdft = xdft(1:N/2+1);                %  Positive frequency 
psdx = (1/(Fs*N)) * abs(xdft).^2;    %  Calculate the power spectral density 
psdx(2:end-1) = 2*psdx(2:end-1);     %  ride 2(2:end-1)
freq = 0:Fs/length(x):Fs/2;          %  Frequency scale 
subplot 211
plot(freq,10*log10(psdx),'k')        %  Take logarithm to plot 
grid on; xlim([0 Fs/2]);
title(' use FFT Periodic graph of ')
xlabel(' frequency /Hz')
ylabel(' Power spectral density /(dB/Hz)')
%  call periodogram function 
[Pxx,f]=periodogram(x,rectwin(length(x)),N,Fs);
subplot 212
plot(freq,10*log10(Pxx),'k');        %  Take logarithm to plot 
grid on; xlim([0 Fs/2]);
title(' call periodogram Periodic graph of function ')
xlabel(' frequency /Hz')
ylabel(' Power spectral density /(dB/Hz)')
mxerr = max(psdx'-Pxx)               %  Find the maximum difference between the two methods 
set(gcf,'color','w');