% Gaussian FFT demonstration
% 5-14-2008

clear all; close all;

alpha=1;  % Parameter value in the Gaussian
L=20;  % Define the computational domain [-L/2,L/2]
n=128;   % Define the number of Fourier modes 2^n

x2=linspace(-L/2,L/2,n+1);  % Define the domain discretization
x=x2(1:n);  % Consider only the first n points because of periodicity

u=1/(alpha*sqrt(2*pi))*exp(-x.*x./(2*alpha^2));  % Gaussian function
ut=fft(u);  % Take the FFT of our Gaussian function
utshift=fftshift(ut);  % Shift FFT

subplot(1,3,1), plot(x,u)  % Original function
xlabel('x'); ylabel('u')
subplot(1,3,2), plot(real(ut))  % Unshifted transform
axis([0 128 -15 15]);
xlabel('Fourier Modes'); ylabel('real(ut)')
subplot(1,3,3), plot(-64:63, real(utshift)) %  Shifted transform
axis([-64 64 -15 15]);
xlabel('Fourier Modes'); ylabel('real(utshift)')

figure(2)
k=(2*pi/L)*[0:(n/2-1) (-n/2):-1]; % Vector of wave numbers
ut1=ifft(-i*k.*ut);  % First derivative
ut2=ifft(k.^2.*ut);  % Second derivative
ut3=ifft(-i*k.^3.*ut);  % Third derivative
plot(x,ut1,'r',x,ut2,'b',x,ut3,'k')
xlabel('x'); ylabel('u^(n)')
legend('u^(1)','u^(2)','u^(3)')