# Finding PSD from Autocorrelation, FFT, Periodogram and Pwelch

31 views (last 30 days)
Yada Kijsathan on 27 Aug 2022
Answered: David Goodmanson on 3 Sep 2022
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

David Goodmanson on 3 Sep 2022
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');

### Categories

Find more on Parametric Spectral Estimation in Help Center and File Exchange

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!