Finding PSD from Autocorrelation, FFT, Periodogram and Pwelch

Dear all,

I wanted to compare the PSD resulting from using Autocorrelation technique, FFT, Periodogram and Pwelch.

So far, the PSD from FFT, Periodogram and Pwelch are the same. However, I get PSD result from Autocorrelation differently in amplitude.

Moreover, after some certain frequency, why the PSD is constant?

I would be appreciated if someone could help me out.

Here's my matlab script

Sen1_acc_x = stim(:,1); %% A given dataset

fs = 20000;

% Using FFT

N = length(Sen1_acc_x);

xdft = fft(Sen1_acc_x);

xdft = xdft(1:N/2+1);

psdx = (1/(fs*N)) * abs(xdft).^2;

psdx(2:end-1) = 2*psdx(2:end-1);

freq = 0:fs/length(Sen1_acc_x):fs/2;

figure(1);

subplot(4,1,1);

plot(freq,10*log10(psdx));

grid on

title('PSD Using FFT')

xlabel('Frequency (Hz)')

ylabel('Power/Frequency (dB/Hz)')

% Using Periodogram

subplot(4,1,2);

periodogram(Sen1_acc_x,rectwin(length(Sen1_acc_x)),length(Sen1_acc_x),fs)

% Using Pwelch

[Sfs1x,fks1x] = pwelch(Sen1_acc_x,window(@rectwin,length(Sen1_acc_x)),[],length(Sen1_acc_x),fs);

subplot(4,1,3);

plot(fks1x,10*log10(Sfs1x));

grid on

title('PSD Using Pwelch');

% Using corr

acf = xcorr(Sen1_acc_x,Sen1_acc_x);

m = length(acf);

n = pow2(nextpow2(m));

pow = fft(acf,n);

f = (0:n-1)*(fs/n); % frequency vector

subplot(4,1,4);

plot(f(1:floor(n/2)),10*log10(pow1:floor(n/2))));

xlabel('Frequency (Hz)')

ylabel('Power')

title('PSD using acf');

grid on

Answer 1

[David Goodmanson]

Hello YK,

Since you are using ffts, the appropriate correlation is going to be circular correlation. This can be done by using

xcorr(Sen1_acc_x, [Sen1_acc_x   Sen1_acc_x])

(assuming it's a row vector) and pulling out the approprate part, but I used ffts instead. Using ffts to obtain the circular correlation that you are going to plug into an fft leads to kind of a circular argument ha ha but it works all right.

I don't have your data so I used random numbers for the input, and all four plots match. The amplitude issue arose because you need to multiply by the same factor that you used in the fft calculation of subplot 1.

Not having your data it's hard to comment on the fact that the PSD goes to a constant at a certain point. Perhaps the data went through a lowpass filter.

% Using circular correlation
% N and freq are borrowed from 'Using FFT'
acf = ifft(fft(Sen1_acc_x).*fft(flip(Sen1_acc_x)));  % circular corr
acf = circshift(acf,1);        % put zero lag value at the beginning of the array
S = real(fft(acf));            % remove small nuisance imaginary part 
S = (1/(fs*N))*S;              % factor that's needed 
S = S(1:(N/2)+1);
S(2:end-1) = 2*S(2:end-1);
subplot(4,1,4);
plot(freq,10*log10(S))
xlabel('Frequency (Hz)')
ylabel('Power')
title('PSD using acf');