# Calculate PSD using FFT

조회 수: 441 (최근 30일)
Peppe 2011년 11월 21일
편집: Aditya 2023년 10월 5일
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
##### 댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Noel Khan 2020년 10월 16일
Do you have the pwelch implementation? I would like to see how you got your estimate if you dont mind :)

댓글을 달려면 로그인하십시오.

### 채택된 답변

Frantz Bouchereau 2021년 7월 29일
Here is a popular MATLAB doc page that explains the relationship between FFT and true power spectra: Power Spectral Density Estimates Using FFT.
##### 댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Aditya 2023년 10월 5일
편집: Aditya 2023년 10월 5일
This page appears to have a conceptual error in the first example. Amplitudes must be multiplied by two before taking one-sided PSD (as pointed out in reply by Chris Schwarz on 27 Aug 2020). Can the documentation be fixed?

댓글을 달려면 로그인하십시오.

### 추가 답변 (4개)

Wayne King 2011년 11월 21일
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).
##### 댓글 수: 0이전 댓글 -2개 표시이전 댓글 -2개 숨기기

댓글을 달려면 로그인하십시오.

Wayne King 2011년 11월 21일
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.
##### 댓글 수: 2없음 표시없음 숨기기
Wayne King 2011년 11월 21일
If you want to use a window, like Hamming, etc. You can do:
plot(psd(spectrum.periodogram('Hamming'),x,'Fs',Fs,'NFFT',length(x)));
Chappi 2019년 12월 17일
Hi, can I ask about how you can apply Hamming window using FFT method? I dont have Signal Processing Tool so I can not use periodogram function. Thank you

댓글을 달려면 로그인하십시오.

Peppe 2011년 11월 21일
Why xdft = xdft(1:length(x)/2+1); % what does that mean? xdft(2:end-1) = 2*xdft(2:end-1); %? dubbel side? and why?
I dont understand these two lines psdest = 1/(length(x)*Fs)*abs(xdft).^2; freq = 0:Fs/length(x):Fs/2;
##### 댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Peppe 2011년 11월 21일
Shouldnt I multiply the hanning windows with my time vector?

댓글을 달려면 로그인하십시오.

Chris Schwarz 2020년 8월 27일
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;
##### 댓글 수: 2없음 표시없음 숨기기
PK 2021년 1월 18일
psdest = abs(xdft/length(x)).^2/(2*df); % modified
why '2*df'? For the zero frequency also 2*df?
Thanks
karinkan 2021년 3월 7일
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;

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Parametric Spectral Estimation에 대해 자세히 알아보기

### Community Treasure Hunt

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

Start Hunting!

Translated by