Why am I getting wrong spectrum applying FFT comparing with correct spectrum?

8 次查看(过去 30 天)
Hello! I am trying to calculate the spectrum of an exponentially dampened oscillating field but I got the wrong result comparing with correct spectrum (different peak amplitude and FWHM). However, I tried with sine function and got a correct result. So I am not sure which part goes wrong (fft part or definition of correct spectrum with MATLAB). Thank you very much!
clear all
clc
Fs = 100; % Sampling frequency
T = 1/Fs; % Sampling period
duration= 50; %second
L=duration*Fs; %sampling length
t = 0:1/Fs:duration-1/Fs; % Time vector
S = 0.7*exp(-t/5).*exp(-i*2*pi*5*t); % Original signal
figure (1)
plot(t,S) % plot original signal in time domain
Y = fft(real(S)); % fast fourier transform
P2 = abs(Y/L); % return it back to correct amplitude
P1 = P2(1:L/2+1); % one side spectrum
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L; % Frequency vector
figure (2)
plot(f,P1) % plot spectrum in frequency domain
w = 2*pi*Fs*(0:(L/2))/L; % angular Frequency vector
w1=(0:0.001:100);
FFT=(0.7/(2*pi))*(1./((1/5)-i*(w1-2*pi*5))); % define correct spectrum
figure (4)
plot(w1,real(FFT)) % plot correct spectrum in angular frequency domain
figure (5)
plot(w,P1) % plot spectrum in angular frequency domain

采纳的回答

David Goodmanson
David Goodmanson 2022-11-1
Hi Jiang,
There are a couple of things going on here. The expression for FFT is not correct. Since you are taking the real part of S, the transform is on the cosine, (1/2)(exp(+...)+exp(-...)) . The answer, Integral (......) dt, is
0.7*(1/2)* (1./((1/t0)-i*(w1-w0)) + 1./((1/t0)-i*(w1+w0))); % with t0 = 5, w0 = 2*pi*5
There is no factor of 2pi in the denominator of this expression.
The fft approximates this integral. Taking the fft does the sum, and that result is multiplied by T, the width of each interval. So Y gets a factor of T instead of (1/L). Dividing by L is used in lots of situations, but not this one.
With the changes, the plots agree.
Fs = 100; % Sampling frequency
T = 1/Fs; % Sampling period
duration= 50; %second
L=duration*Fs; %sampling length
t = 0:1/Fs:duration-1/Fs; % Time vector
S = 0.7*exp(-t/5).*exp(-i*2*pi*5*t); % Original signal
figure (1)
plot(t,S) % plot original signal in time domain
Y = fft(real(S)); % fast fourier transform
% P2 = abs(Y/L); % return it back to correct amplitude
P2 = abs(Y)*T;
P1 = P2(1:L/2+1); % one side spectrum
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L; % Frequency vector
figure (2)
plot(f,P1) % plot spectrum in frequency domain
w = 2*pi*Fs*(0:(L/2))/L; % angular Frequency vector
% one sided spectrum
w1=(0:0.001:100);
w0 = 2*pi*5;
t0 = 5;
FFT = 0.7*(1/2)*( 1./((1/t0)-i*(w1-w0)) + 1./((1/t0)-i*(w1+w0)) );
FFT = 2*FFT; % one sided
% plot abs(FFT) instead of real(FFT) since it is a fairer comparison to abs(Y).
figure(4)
plot(w1,abs(FFT)) % plot correct spectrum in angular frequency domain
figure (5)
plot(w,P1) % plot spectrum in angular frequency domain

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Spectral Measurements 的更多信息

产品

Community Treasure Hunt

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

Start Hunting!

Translated by