different way to compute power spectral density

Question

Renaud Dufour 2015-1-6

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/169078-different-way-to-compute-power-spectral-density

评论： Soumak Bhattacharjee 2019-11-28

Dear everyone,

I am writing a piece of code to compute power spectral density estimate of a signal and wanted to compare two approaches :

compute the FFT of the signal and square its amplitude
compute the biased Autocorrelation function and take the FFT

I expected the two methods to give a similar result, however for the second one the Parseval sum in frequency space is wrong.

% Define spatial signal
dx   = 1e-3;
Fsx  = 1/dx; % sampling frequency
X    = 0:dx:2;
L    = length(X);     % = 2001
Lfft = 2^nextpow2(L); % = 2048
dfx  = 1/(Lfft*dx);
T   = 0.1;
Z   = sin(2*pi*X/T);
% Compute PSD from FFT
fourier_image  = fftshift(fft(X,Lfft)/Lfft);
PSDfromFFT     = abs(fourier_image).^2;
freq           = -Fsx/2:dfx:Fsx/2-dfx; % Frequency axis
% Compute PSD from Autocorrelation
Rxx = xcorr(X,'biased');
lag = (-(L-1):1:(L-1))*dx;
PSDfromRxx = abs(fftshift(fft(Rxx,Lfft)/Lfft));
% Parseval sums
fprintf('\nParseval sums : \n\n');
fprintf('Spatial sum = %f\n', sum(X(:).^2));
fprintf('Frequency sum (from FFT) = %f\n', sum(PSDfromFFT(:))*Lfft);
fprintf('Frequency sum (from ACF) = %f\n', sum(PSDfromRxx(:))*Lfft);
% Plots
figure;
subplot(2,1,1); plot(X,Z,'-b'); title('Signal z(x)');
subplot(2,1,2); plot(lag,Rxx,'-r'); title('Autocorrelation');
figure;
plot(freq,PSDfromFFT,'-or'); title('PSD based on FFT'); hold on;
plot(freq, PSDfromRxx, '-ob'); title('PSD based on autocorrelation');
legend('PSD based on FFT','PSD based on autocorrelation'); hold off;

I get the following result for the Parseval sums :

Parseval sums : 
Spatial sum = 2668.667000
Frequency sum (from FFT) = 2668.667000
Frequency sum (from ACF) = 7853.268426

So the first method is ok but there is a factor of ~3 for the second one, and I cannot figure out what is the problem. Note that I scale the fft with 1/N which is a convention for the project I'm working on.

Thanks in advance for the inputs