about FFT and wave packet
2 views (last 30 days)
Show older comments
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.
0 Comments
Answers (1)
Chris Turnes
on 2 Sep 2015
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.
2 Comments
Chris Turnes
on 4 Sep 2015
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.
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!