about FFT and wave packet

Hi there, I am trying to verify the fourier transformation on the Gaussian in momentum and position space. The analytical Gaussian in momentum space is

g = sqrt(1/(sqrt(pi)*s))*exp(-p^2/(4*s^2)), where s is the standard deviation. If we integrate g^2, we see that g^2 is normalized to 1.

The inverse fourier transformation of g is the Gaussian in position, which is given analytically as f f = sqrt(2*s/sqrt(pi))*exp(-2s^2*x^2)

Again, if we integrate f^2, we got the normalization to 1. I try the Fourier transformation of f in mathematica and I reproduce g. Now I need to do that numerically in matlab, I got the following code

N  = 2048;                    
dx = 0.1963;               
x  = dx*[-N/2:(N/2-1)];    
p  = 0.0156*[-N/2+1:N/2];  
dp = abs(p(2)-p(1));

s = 0.2;   
f = sqrt(2*s/sqrt(pi))*exp(-2*s^2.*x.^2); 
g = fftshift(fft(f))/sqrt(N); % numerical solution
gt = 1/(sqrt(2)*pi^(1/4)*sqrt(s))*exp(-p.^2/(8*s^2));  % analytical expression
sum(abs(f).^2*dx)     % we see |f|^2 is normalized to one
sum(abs(g).^2*dp)     % |g|^2 is not normalized to one 
subplot(2,1,1); plot(x, abs(f).^2);
subplot(2,1,2); plot(p, abs(g).^2); hold on; plot(p, abs(gt).^2, 'r:');

I don't understand why g is not same as gt, which is the analytical solution.

Answer 1

[Chris Turnes]

You are using different scaling factors for the two expressions; the FFT version is scaled with dx = 0.1963 and the "analytical expression" with dp = 0.0156.

First, you'll notice that while sum(abs(g).^2*dp) is not equal to 1, if you correct the scaling factor to use dx (which appears to define the grid over which g is calculated), you actually will get the proper normalization:

sum(abs(g).^2*dx)
ans =
   1.0000e+00

Second, the analytical grid is defined over -N/2+1:N/2 while the numerical grid is -N/2:(N/2-1). If you're expecting the same results, the grids should be defined over the same point range.

To get the correct normalization then, you just need to re-scale the output of the FFT to be relative to dx, not dp:

g = g*sqrt(dx/dp);

Once you do that, you'll see

sum(abs(g).^2*dp)
ans =
   1

and the plots will look correct.

[Ki]

It looks good now. Just have a quick question. If I start with the function f, let it evolves in time, so during the evolution I need to check how the fourier transform on f (i.e. g) looks like, so I need to multiply the fft(f) with sqrt(dx/dp) each time to get a right normalization?

[Chris Turnes]

I suppose the normalization that you need will depend on exactly how the problem is set up and what you're trying to compare. I think the main issue here is just making sure that the time grid of the analytical solution and the time grid of the numerical solution that you are comparing against are the same.