Correct way to evaluate the PSD

Dear All,

In order to evaluate the psd of a signal I'm using the following code:

%%Generating the signal in time domain
A=46;                   % [V] amplitude of the signal
fc=1E9;                 % [Hz] signal frequency
T=50E3;                 % [pts] number of points in time domain
fs=5E9;                 % [Hz] sampling frequency in time domain
time=0:1/fs:T/fs-1/fs;  % time vector
Sig=A*sin(2*pi*fc*time);% signal in time domain

%%Converting the signal in frequency domain
NFFT=2^nextpow2(T);     % adding necessary zeros

% Frequency transform & scaling the result so that it is not a function of the length of Sig
FSig=fft(Sig, NFFT)/T;  
% Taking only half of spectrum
FSig=FSig(1:NFFT/2);
% Multiply with sqrt(2) to keep the same energy (because it's computed as one side); The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2. 
if rem(NFFT, 2)     
    FSig(2:end) = FSig(2:end)*sqrt(2);      % odd NFFT excludes Nyquist point 
else
    FSig(2:end -1) = FSig(2:end-1)*sqrt(2);  % even NFFT
end
% Correcting with sqrt(N/NFFT) so it will not depend on zero padding length
FSig=sqrt(T/NFFT)*FSig;

Magnitude=abs(FSig);                        % signal magnitude in freqeuncy domain
f=fs/2*linspace(0, 1, NFFT/2);              % frequency vector

%%Estimating the PSD
PSD=FSig.*conj(FSig)*T/fs;

Now, evaluating the psd of the same time signal using the psd function found in Matlab:

h = spectrum.welch;       % Instantiate a welch object.
    psd(h,Sig,'fs',fs);   % Plot the one-sided PSD.

the result is smoother. Reading on Matlab functions I found that the window in psd.m results from psdchk.m which in turn is created by a Hanning window (whose size is variable).

My question is how is correct to evaluate the psd ? I need to create a signal in time domain having a smooth psd.

Regards,

AndMJ

Answer 1

[Wayne King]

You can't compare your code to a Welch estimate. A Welch estimate does not simply use a window. A Welch estimate breaks up the data into (usually overlapped) windowed segments, estimates the PSD for each segment, and then averages those segments together.

The Welch estimate greatly reduces the variance of the PSD estimate (by increasing the degrees of freedom) compared to the periodogram (or modified periodogram). Your code using fft() essentially produces a periodogram.

To get a smooth PSD you have a few choices:

1.) Compute a parametric PSD using pburg() for example. These estimates are very smooth, however (big however), you have to specify a model for the PSD and it is very sensitive to model misspecification.

2.) Use pwelch() or pmtm() -- these are nonparametric methods that estimate the PSD. They smooth the PSD estimate in fundamentally different ways, but I think you can use either. I prefer pmtm() but pwelch() is also good.

I would recommend option #2