Gaussian FFT

11 次查看（过去 30 天）
George Petrov
George Petrov 2012-6-4
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 个评论
显示 1更早的评论
Oleg Komarov
Oleg Komarov 2012-6-5
Most importantly i think "t_anal" and similar are NOT good names.
George Petrov
George Petrov 2012-6-7
Hello Oleg,
Sorry about the name. "_anal" comes from "analytical". I can replace it with something else and
put the code back, if you would like to take a look and run it.
George

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Dr. Seis
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 个评论
George Petrov
George Petrov 2012-6-7
Hello Elige,
Thanks for your suggestion. When I plot the absolute value of the fouirer transform in the frequency domain it is a Gaussian. Without the absolute value the signal alternates between positive and negative. It jumps from + to -, then + and - again. The code I attached looks awful. I can attach it again, if you would like. This FFT example bothers me because it is supposed to be really straightforward.
Thanks: George
Dr. Seis
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 个）

George Petrov
George Petrov 2012-6-11
Hello Elige, Thanks a lot. I think I got it: using the symmetry property will ensure a real-valued symmetrical spectrum. Indeed, I see a small imaginary component, probably due to the phase I introduce. I will give it a try with a fully symmetrical profile. Thanks again. George

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by