447 views (last 30 days)

Show older comments

Hi, The question is to calculate PSD using FFT function in MATLAB. Ive already done it with pwelch command in MATLAB and now it's time to do it with FFT command and compare the results. If I have file named: file2.Mat which contains 3 columns. first column is time, second Force and the third is acceleration. the sampling is 4000Hz and the number of NFFT is ,let us say, 4444. I know that we need to multiply the window with time column. And then what?

Does anybody know how to do it? Ive been work on this for like 3 hours.

regards

Noel Khan
on 16 Oct 2020

Frantz Bouchereau
on 29 Jul 2021

Wayne King
on 21 Nov 2011

Why don't you just use spectrum.periodogram?

Fs = 1e3;

t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*100*t)+randn(size(t));

plot(psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',length(x)));

If you want to do it simply with fft()

xdft = fft(x);

xdft = xdft(1:length(x)/2+1);

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

psdest = 1/(length(x)*Fs)*abs(xdft).^2;

freq = 0:Fs/length(x):Fs/2;

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

grid on;

Compare the plots.

Chappi
on 17 Dec 2019

Wayne King
on 21 Nov 2011

The PSD is an even function of frequency, so you only need from 0 to the Nyquist, if you want to conserve the total power, you have to multiply all frequencies except 0 and the Nyquist by two if you only keep 1/2 the frequencies. 0 and the Nyquist only occur once in the PSD estimate, all other frequencies occur twice. If you look at the example I gave you, then you see it agrees with the scaling in MATLAB's periodogram.

The answer about multiplying by a window, Hanning, Hamming, Blackman, Tukey. etc. is that it depends. A window reduces the bias in the periodogram, but that comes at the cost of reduced frequency resolution (a broader main lobe).

Chris Schwarz
on 27 Aug 2020

An adjustment to Wayne's code that gives an exact match to the periodogram is:

Fs = 1e3;

t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*100*t)+randn(size(t));

plot(psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',length(x)));

xdft = fft(x);

xdft = xdft(1:length(x)/2+1);

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

freq = 0:Fs/length(x):Fs/2;

df = diff(freq);

df = df(1);

% psdest = 1/(length(x)*Fs)*abs(xdft).^2; % original

psdest = abs(xdft/length(x)).^2/(2*df); % modified

hold on

plot(freq,10*log10(psdest),'r');

grid on;

karinkan
on 7 Mar 2021

close all

Fs = 1e3;

t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*100*t)+randn(size(t));

L =length(x);

NFFT = 2^nextpow2(L);

plot(psd(spectrum.periodogram,x,'Fs',Fs,'NFFT',NFFT));

df = Fs/NFFT;

freq = 0:df:Fs/2;

xdft = fft(x,NFFT);

xdft_s = xdft(1:NFFT/2+1);

amp = abs(xdft_s)/NFFT;

psdest = amp.^2/(df); % original

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

hold on

plot(freq,10*log10(psdest),'r');

grid on;

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

Start Hunting!