Normalization for the fft of a wavefunction in momentum space
10 views (last 30 days)
Show older comments
PETER SWOVICK
on 8 Apr 2020
Commented: David Goodmanson
on 9 Apr 2020
I am currently working on an fft that takes a guassian wavefunction from position space to momentum space by doing an fft. But when I do this, for some reason, I do not
seem to get a normalized state using the typical dividing by the length of sampling as seen in many documentation examples. If I force a state to be normalized through
pVecFFTnorm=pVecFFTTrial./(sqrt(sum(abs(pVecFFTTrial).^2)*dp));
I get a result for a quantity that I do not expect (I get a pExpectSquare, shown in the code below, that deviates significantly from 0.5) What could be done to ensure that my norm will be 1 for my wavefunction? The code is provided below:
%Note: the parameters of a=1, hbar=1 p0=0, and x0=0 have been numerically assumed
%for wave pack and all other functions below.
N=2048;
clf
wavePack=@(x)((1/pi)^(1/4)*exp((-(x).^2)/2));
xDomFFT=linspace(-20,20,N);
xVecFFT=wavePack(xDomFFT);
dx=(xDomFFT(2)-xDomFFT(1));
normBeforeFFT=sum(xVecFFT.^2)*dx;
%xPoints = 2^nextpow2(xPoints);
pVecFFTpreShift=fft(xVecFFT,N);
pVecFFT=fftshift(fft(xVecFFT));
%the fft may not necesarily be normalized. This line is to make sure that
%the fft is normalized.
pDom=linspace(-N/2+1,N/2,(N));
dp=abs(pDom(2)-pDom(1));
pVecFFTTrial=pVecFFT/N;
pVecFFTnorm=pVecFFTTrial./(sqrt(sum(abs(pVecFFTTrial).^2)*dp));
plot(pDom,pVecFFTnorm);
xlim([-25,25]);
xlabel("p")
ylabel("\psi(p)")
title("\psi(p) on a grid")
pPopFFT=abs(pVecFFTTrial).^2;
%plot(pDom,pPopFFT)
IsNorm=sum(pPopFFT)*dp %This should be 1, I am getting this to be some other number.
%NewNorm=std(pVecFFTTrial)
pExpect=pPopFFT.*pDom;
expectP=sum(pExpect)*dp %This should be close to 0
%plot(pDom,pExpect)
pExpectSquare=pPopFFT.*(pDom.^2); %This should be close to 0.5
0 Comments
Accepted Answer
David Goodmanson
on 8 Apr 2020
Hi Peter,
the factor of 1/N is used to find the amplitudes of continuous wave sines and cosines. Here you are approximating a fourier integral with the fft, which is a different matter. If you want to see .5 for the second moment in the momentum domain, then the appropriate transform pair is
pVec(p) = 1/sqrt(2*pi)) * Int xVec(x)*exp(-i*p*x) dx
xVec(x) = 1/sqrt(2*pi)) * Int pVec(p)*exp( i*p*x) dp
The only requirement is that the product of the two factors in front equals 1/(2*pi), and with the symmetric scaling here you will get .5 in the p domain.
You are approximating this with fft arrays of discrete spacing. If you were using for example exp(2*pi*i*f*t) in the transform, the golden rule for an N-point fft with array spacings dt and df is
dt*df = 1/N [1]
Since you are using exp(i*k*x) the 2*pi is absorbed into k to give
dx*dk = 2*pi/N
The code below just reproduces what is above. I changed the number of points to 2000 since I think the 2^n preference is mostly useless. WIth k and x it doesn't matter as much, but from [1] you can see that if N is a number like 2000 then both dt and df could have convenient spacings.
N = 2000;
wavePack = @(x)((1/pi)^(1/4)*exp((-(x).^2)/2));
x = linspace(-20,20,N+1);
x(end) = []; % N points, of the form (-N/2:N/2-1)*dx
xVec = wavePack(x);
figure(1)
plot(x,xVec)
grid on
dx=(x(2)-x(1));
normBeforeFFT = sum(xVec.^2)*dx
% approximate the fourier transform integral
pVec_nosh = (1/sqrt(2*pi))*fft(xVec)*dx;
pVec = fftshift(pVec_nosh);
dp = 2*pi/(dx*N);
p = (-N/2:N/2-1)*dp;
fig(2)
plot(p,abs(pVec).^2)
grid on
Int0 = sum(abs(pVec).^2)*dp
Int1 = sum(p.*abs(pVec).^2)*dp
Int2 = sum(p.^2.*abs(pVec).^2)*dp
normBeforeFFT = 1.0000
Int0 = 1.0000
Int1 = -1.2465e-17
Int2 = 0.5000
2 Comments
David Goodmanson
on 9 Apr 2020
I think maybe you meant multiplying by dx. And that's right, the fft takes an array with indices, all its calculations in terms of indices and it doesn't really know anything about dx. You are reproducing a Riemann integral so you have to multiply by the step width, just like multplying by h in simpson's rule.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!