Normalization for the fft of a wavefunction in momentum space

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

Answer 1

[David Goodmanson]

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

[PETER SWOVICK]

Thanks for your help! With your adjustments everything works! I am just a little bit curious, though, why do you multiply by the dp to the response from the fft? I am guessing that because in the fourier transform above,

pVec(p) = 1/sqrt(2*pi)) * Int xVec(x)*exp(-i*p*x) dx

the fft algorithm would only does the steps of xVec(x)*exp(-i*p*x), right? Thanks for the help!

[David Goodmanson]

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.