how to do fft to a gaussian function

I have a Gaussian wave function that is psi = exp(-x.^2/sigma^2) with sigma = 1e-5 and x range x = -3e-5:1e-7:3e-5. I can get a perfect Gaussian shape by plotting this function. But when I do fft to this equation, I always get a delta function. Should I get a Gaussian function in momentum space? Thanks very much for answering my question. This is my code.

close all; clear all; clc;

sigma = 1e-5;

x = -3e-5:1e-7:3e-5;

psi = exp(-x.^2/sigma^2);

A = 1/sqrt(sum(abs(psi).^2*1e-7));

psi = A.*psi;

figure

plot(x,psi);

xlim([-3e-5,3e-5]);

xlabel('x axis (m)')

ylabel('psi')

title('gaussian distribution sigma = 10 um')

phi = fft(psi);

phi = fftshift(phi);

k = linspace(-1e-100,1e-100,length(phi));

figure

plot(k,abs(phi))

<<

>>

[Rahul]

I have noticed something very similar, did you find a resolution to the problem ?

Answer 1

[David Goodmanson]

Hello Shijie,

This is happening because of basic behavior of ffts. For an N-point fft, suppose f(k) is the fourier transform of f(x) and let

delx, delk  = spacing of x array, k array      respectively
Wx, Wk  = characteristic width of f(x), f(k)        "
nx, nk  = number of points across Wx, Wk            "

In your case Wx = sigma/2, approximately. These quantities obey

nx = Wx/delx;   nk = Wk/delk
delx*delk = 2*pi/N                 % fundamental rule for fft
Wx*Wk ~~ 1/2,                      % Heisenberg uncertainty principle

nx*nk = N/(4*pi)                   % another fft rule, which follows from lines above

To plot each function well, one needs a significant number of points in both nx and nk. In your case nx = 50 by construction, N = 601 so nk is 1. You get almost a delta function for f(k) as you saw. To make nk larger, you have to either decrease nx (case 2 below) or increase N (case 3). I used individual markers on the plots and removed the limit on the x axis so in some cases you will have to zoom in to see what is happening.

As for the k axis, I don't know what is going on with the 1e-100 scaling. Your gaussian width is on the order of 1e-5, so your k values are going to be around 1e5, which is what the code below does.

sigma = 1e-5;
delx = 1e-7; x = -3e-5:delx:3e-5;       % case 1, original
%delx = 1e-6; x = -3e-4:delx:3e-4;       % case 2
%delx = 1e-7; x = -6e-4:delx:6e-4;       % case 3
N = length(x);
psi = exp(-x.^2/sigma^2);
A = 1/sqrt(sum(abs(psi).^2*delx));
psi = A.*psi;
figure(1)
plot(x,psi,'o');
%xlim([-3e-5,3e-5]);
xlabel('x axis (m)')
ylabel('psi')
title('gaussian distribution sigma = 10 um')

Int_x = sum(psi.^2)*delx  % should be 1

phi = fft(psi)*delx;    % multiply by delx to approximate the fourier integral
phi = fftshift(phi);
delk = 2*pi/(N*delx);
k = (-N/2:N/2-1)*delk;  % new k
figure(2)
plot(k,abs(phi),'o')

Int_k = sum(abs(phi.^2))*delk   % should be 2 pi