Basic question about FFT and PSD
in [Matlab]
|
From: Agus Nugroho on 25 Jun 2010 07:02 Dear everyone, I am sorry to bring such a basic question. I hope you don't mind to answer my questions below: 1. I have calculated the FFT for a sample signal and I got a spectrum with x-axis represented frequency and y-axis represented by amplitude. And the result is correct,, the amplitude and the frequency are correct. However, when I tried to calculate by using the PSD formula (Periodogram), I got a different result. I understand that the PSD has x-axis in which represented by frequency and y-axis is represented by power/frequency. So my question is " what's the different between these two methods, since both of them are used to obtain the frequency spectrum of signal?" In my opinion, the resulted frequency should be same. 2. Is it possible to transform my FFT spectrum directly to PSD spectrum? I mean may I simple by squared the amplitude (Y(f))^2 in FFT to obtain the PSD spectrum? I am sorry for my understanding, and I am grateful for your answer. Regards, Agus
From: Wayne King on 25 Jun 2010 08:04 "Agus Nugroho" <agus_nugroho01(a)yahoo.com> wrote in message <i022bc$nqs$1(a)fred.mathworks.com>... > Dear everyone, > > I am sorry to bring such a basic question. I hope you don't mind to answer my questions below: > > 1. I have calculated the FFT for a sample signal and I got a spectrum with x-axis represented frequency and y-axis represented by amplitude. And the result is correct,, the amplitude and the frequency are correct. However, when I tried to calculate by using the PSD formula (Periodogram), I got a different result. > > I understand that the PSD has x-axis in which represented by frequency and y-axis is represented by power/frequency. So my question is " what's the different between these two methods, since both of them are used to obtain the frequency spectrum of signal?" In my opinion, the resulted frequency should be same. > > 2. Is it possible to transform my FFT spectrum directly to PSD spectrum? I mean may I simple by squared the amplitude (Y(f))^2 in FFT to obtain the PSD spectrum? > > I am sorry for my understanding, and I am grateful for your answer. > > Regards, > > Agus Hi Agus, To compute the periodogram using fft(), you have to scale the absolute value squared DFT coefficients by (1/length(data))*(1/Fs) where Fs is the sampling frequency if your PSD estimate is in power/Hz, or (1/length(data))*(1/(2*pi)) if your PSD estimate is power/(radians/sample). If your PSD estimate is one-sided (say positive frequencies only, which is a very common thing to do), then you have an additional factor of 2 for all frequencies except zero (DC) and the Nyquist. Remember, the PSD estimates are scaled by the reciprocal of the length of the data (not the length of the data after any zero padding). I’ll try to explain with code, which we’ll check against spectrum.periodogram() and psd() in the Signal Processing Toolbox just to show you the result is consistent with the SPT routines. Case 1: Two-sided PSD estimate with sampling frequency. PSD estimate in power/Hz (dB) • The scaling factor is: dt/N where dt is the sampling interval and N is the length of the input data ignoring any padding of course- dt/N is equivalent to 1/(N*Fs). Since the PSD estimate is two-sided, there is no extra factor of two required. reset(RandStream.getDefaultStream); dt = 1/10000; %Fs=10000 t = 0:dt:.1024-dt; %time vector of length 1024 y = cos(2*pi*1000*t)+randn(size(t)); %1 kHz sinewave in noise—signal for all %examples. h = spectrum.periodogram; %spectrum object (periodogram) %two-sided PSD with sampling frequency specified psd_est_2side=psd(h,y,'spectrumtype','twosided','Fs',10000); %now let’s do it by brute force and check agreement y_fft=abs(fft(y)).^2; psd_y=(dt/length(y))*y_fft; %plot the comparison between PSD estimates plot(psd_est_2side.Frequencies./1000,10*log10(psd_y)) xlabel('kHz'); ylabel('Power/Hz (dB)'); grid on; figure; plot(psd_est_2side); Case 2: One-sided PSD with sampling frequency—For real-valued data, the DFT is conjugate symmetric. Therefore, most people only plot the positive frequencies. If you want Parseval’s relation to still hold you have to multiply all the positive frequencies which normally also occur as negative frequencies by a factor of 2. This is because you’ve essentially “lost” half your power by omitting the negative frequencies. • The scaling factor for a one-sided PSD estimate with a known sampling frequency (Power/Hz) is: (2*dt)/length(data) Note: The zero frequency estimate (DC) and the Nyquist frequency estimates are NOT multiplied by 2. DC does not occur as a negative frequency. In a PSD estimate with the DC component centered, Matlab computes estimates over the semi-open interval (-Nyquist,Nyquist], therefore the Nyquist frequency also does not occur twice. %onesided estimate psd_est_1side=psd(h,y,'spectrumtype','onesided','Fs',10000); y_fft=abs(fft(y)).^2; Y_fft=y_fft(1:513); % we only keep the positive frequencies. Y_fft=(dt/1024)*Y_fft; Y_fft(2:end-1)=2*Y_fft(2:end-1); %additional factor of two for all but DC and the %Nyquist %compare plot(psd_est_1side.Frequencies./1000,10*log10(Y_fft)); xlabel('kHz'); ylabel('Power/Hz (dB)'); grid on; figure; plot(psd_est_1side); Case 3: PSD estimate for normalized frequency-- (Power (in dB)/rad/sample) If normalized frequency is used, the scaling factors are: • 1/length(data)*(1/(2*pi)) for a two-sided estimate • 2/length(data)*(1/(2*pi)) for a one-sided estimate (again, see note concerning DC and the Nyquist!) % doing the two-sided case -- not centering DC psd_est_2side_normalized=psd(h,y,'spectrumtype','twosided'); y_fft=abs(fft(y)).^2; Y_fft=(1/length(y))*(1/(2*pi))*y_fft; plot(psd_est_2side_normalized.Frequencies./pi, 10*log10(Y_fft)); xlabel('Normalized Frequency (x \pi rad/sample)'); ylabel('Power/rad/sample (dB)'); grid on; figure; plot(psd_est_2side_normalized); The one-sided case with normalized frequency requires the additional factor of two except at DC and the Nyquist. Hope that helps, Wayne
From: Agus Nugroho on 25 Jun 2010 18:11 Dear Dr. Wayne, Thank you so much for your answer. Its really clear and detail. As soon, I will try to implement and test to my data. Thank you very much for your kindness. Regards, Agus
|
Pages: 1 Prev: LDA ouput as classes seperated well in plot Next: Canon command and state names |