MATLAB: Normalization for the fft of a wavefunction in momentum space

fftnormalization

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

Best Answer

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
Related Question