Gaussian FFT
6 次查看(过去 30 天)
显示 更早的评论
I am using the Matlab fft() function to get FFT from a Gaussian, y(t)=exp(-a*t^2) and compare to the continuous Fourier transform. With other words, I use FFT to approximate CFT. CFT of a Gaussian is also a Gaussian. Surprisingly, the FFT resulted in a spectrum that oscillates between positive and negative values. Why FFT does not approximate the CFT ? I would appreciate your comments. Below is a small code I wrote. %========================================================================== % % FFT of gaussian: Comparison of Fourier transform with DFT % % FFT Fourier pair: y(t)=exp(-a*t^2) = Y(f)=(pi/a)^(1/2)*exp[-(pi*f)^2/a] % E. O. Brigham, "The FFT and its applications", Prentice-Hall, Inc., 1988 %========================================================================== % ----------------------- set up signal parameters ------------------------ a=4; % set parameter "a" n=256; % number of grid nodes Tmax=16.0; % set Tmax % --------------------------- create a signal ----------------------------- dt=Tmax/(n-1);t=-Tmax/2:dt:Tmax/2; % set a grid in the time domain y=exp(-a*t.^2); % map signal on grid nodes % ---------------------------- frequency grid ----------------------------- fo=1.0/Tmax;f_max=1.0/(2*dt); % sampling and maximum frequency f(1:n+1)=-(n/2)*fo:fo:n/2*fo; % set frequency grid % ----------------- get FFT of the discretized function ------------------- % + frequencies: f(1)=0,f(2)=fo,...f(n/2+1)=f_max, % - frequencies: f(n/2+1)=-f_max,f(n/2+2)=-f_max+fo,...,f(n-1)=-2*fo,f(n)=-fo Y=fft(y)*dt; % fft(y)*dt Y1(1:n/2)=Y(n/2+1:n); % set Y to negative frequencies Y1(n/2+1:n+1)=Y(1:n/2+1); % set Y to positive frequencies % ---------- create a discrete version of the Fourier transform ----------- t_anal=-Tmax/2:Tmax/1000:Tmax/2; % set a grid for analytical "y" y_anal=exp(-a*t_anal.^2); % map signal on grid nodes f_anal=-f_max:f_max/1000:f_max; % set frequency grid Y_anal=(pi/a)^0.5*exp(-f_anal.^2*pi^2/a); % analytical FFT(y_anal) % ------------------------------- figures --------------------------------- figure(1);clf;subplot(2,1,1); % subplot plot(t_anal,y_anal,'b:');hold on; % plot signal y_anal(t) plot(t,y,'r.');hold on; % plot signal y(t) xlabel('t(s)');ylabel('f(t)'); % labels legend('analytical','discretized'); % legend title('Input signal'); % figure title axis([-Tmax/2 Tmax/2 -0.1 1.1]); % axis
subplot(2,1,2); % subplot plot(f_anal,Y_anal,'b:');hold on; % plot Re(Y_anal) plot(f,real(Y1),'r.');hold on; % plot Re(Y) xlabel('f (1/s)');ylabel('FFT(y)'); % labels legend('Re(Y) - analytical','Re(Y) - FFT'); % legend title('Real part of Y'); % figure title axis([-f_max f_max -1.1 1.1]); % axis % ---------------------------------- end ----------------------------------
3 个评论
采纳的回答
Dr. Seis
2012-6-5
Try plotting the absolute value of the frequency domain amplitudes... when your Gaussian is not centered at time 0, the phase may be non-zero.
2 个评论
Dr. Seis
2012-6-7
The code isn't necessary here... that's just the way things are. The only way to get the FFT result to look like the CFT result is to shift the peak of the Gaussian to time = 0. Since the FFT doesn't know what the time is associated with your sample points, the only way this can be accomplished is "wrap-around". What I mean is that your time series (of length N samples) will have its' peak at the first sample (t = 0), then the second sample will be equal to your Nth sample, and your third sample will equal the N-1th sample, and the fourth sample will equal the N-2th sample, and so on...
This will have the effect of making your time series more-or-less symmetrical. A real-valued, symmetrical time series will result in a real-valued, symmetrical frequency spectrum. That's just the way it works. As soon as you shift your Gaussian away from time = 0, you will introduce phase. A signal with non-zero phase means your frequency amplitudes will become complex. It should be expected that your frequency amplitudes will "jump around" the more your Gaussian moves away from time = 0... until you get past the middle of your time series, then it will begin to "jump around" less and less until your Gaussian "wraps around" again.
更多回答(1 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!