Gaussian FFT

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 ----------------------------------

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[Oleg Komarov]

Most importantly i think "t_anal" and similar are NOT good names.

[George Petrov]

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

Answer 1

[Dr. Seis]

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.

[George Petrov]

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]

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.

Answer 2

[George Petrov]

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