Normalizing an FFT Vector

Hey all,

I am currently trying to understand how exactly to normalize an fft eigenstate in the context of the 1-d particle in an infinite well. My code is as below:

L = 100;

x = linspace(0,L,L+1);

psi = sqrt(2/L)*sin(pi*x/L); %particle in a box ground state%

norm = trapz(x,psi.*conj(psi)); %normalized state has L-2 norm = 1%

F_T = fft(psi);

norm2 = trapz(2*pi*x/L, F_T.*conj(F_T));

The psi is the discretized array for my particle in position space that, as calculated in "norm", has L-2 norm of 1. Now I wanted to transform this into momentum space via the fft, but when I try to get the L-2 norm of it again as in "norm2", the norm is no longer 1. I have seen suggestions of dividing the fft by online but it doesn't seem to resolve the issue either. Any advice would be appreciated.

Thanks!

Answer 1

[Paul]

Hi bil,

Define psi

syms x real
L = 100;
psi(x) = sqrt(sym(2)/L)*sin(sym(pi)*x/L)*rectangularPulse(0,L,x); %particle in a box ground state%

Its energy computed in the space domain is

Ecs = int(psi(x)*conj(psi(x)),x,-inf,inf)

Ecs = 

Hmm, int can't handle that with the inf limits. Recompute with finite limits

Ecs = int(psi(x)*conj(psi(x)),x,-1e10,1e10)

Ecs = 

1

Its Continuous Fourier Transform (CFT)

syms w real
Psi(w) = simplify(fourier(psi(x),x,w))

Psi(w) = 

figure
fplot(abs(Psi(w)))
hold on
xlabel('rad/sec')

Its energy computed in the frequency domain is

Ecf = int(Psi(w)*conj(Psi(w)),w,-inf,inf)/2/sym(pi)

Ecf = 

Can't find a closed form expression.

Verify the imaginary part of the integrand is 0

imag(Psi(w)*conj(Psi(w)))

ans = 

0

Integrate just the real part

Ecf = int(simplify(real(Psi(w)*conj(Psi(w)))),w,-inf,inf)/2/sym(pi)

Ecf = 

1

As expected, Ecf = Ecs.

Samples of psi. I'm going to change the number of samples to illustrate a point later.

N = 289;
x = linspace(0,L,N);
dx = x(2) - x(1)

dx = 0.3472

psi = double(psi(x)); %particle in a box ground state%

Approximation of the energy integral via trapz yields the exact answer.

Ecs = trapz(x,psi.*conj(psi)) %normalized state has L-2 norm = 1%

Ecs = 1

The energy in the discrete space samples is

Eds = sum(psi.*conj(psi))

Eds = 2.8800

The continuous space energy is approximated by scaling by dx

Ecapproximate = Eds*dx

Ecapproximate = 1.0000

"DTFT" of the samples (quotes becasue the indpendependent variable is space, not time), w in rad/sample

[dtft,w] = freqz(psi,1,linspace(-pi,pi,2048));

Add to the plot. The DTFT has to be scaled by dx to approximate the CFT

plot(w/dx,abs(dtft*dx))

Discrete energy computed from the DTFT

Edtft = trapz(w,dtft.*conj(dtft))/2/pi

Edtft = 2.8800

Approximate Ec after scaling by dx

Ecapproximate = Edtft*dx

Ecapproximate = 1.0000

DFT of the samples of psi

dft = fftshift(fft(psi));

Discrete energy in the frequency domain

Edft = sum(dft.*conj(dft))/N

Edft = 2.8800

Mulitply by dx to approximate Ec

Ecapproximate = Edft*dx

Ecapproximate = 1.0000

Can also approximate Ec via trapezoidal integration of the DFT

Define frequency vector (N is odd), rad/distance

w = ((-(N-1)/2):((N-1)/2))/N*2*pi/dx;

Integrate, divide by 2*pi as required when w is in rad/sample

Ecapproximate = trapz(w,(dft*dx).*conj(dft*dx))/2/pi

Ecapproximate = 1.0000

Add DFT to plot

stem(w,abs(dft)*dx)
xlim([-0.5 0.5])
legend('CFT','DTFT','DFT')

[bil]

Wow you went above and beyond this is great thank you so much!

[bil]

Hi, I'm not sure if you'll still see this, but I transferred the same code to Python and realized that using the Simpson's form of integration caused the calculated norm to be much less than 1. Do you know the cause of this?

[Paul]

I'm afraid I don't. Did you get Simpson's rule to work in Matlab?

[bil]

Sorry, disregard my question I have figured out the error. Thanks again!

[David Goodmanson]

Hi Paul, it's kind of interesting to work out the integral for Ecf, although I could not find a truly simple way. If you insert w = x/100, dw = dx/100 you get

  pi Int (1+e^(ix))(1+e^(-ix)) / (x^2-pi^2)^2 dx
= pi Int (2+e^(ix)+e^(-ix))    / (x^2-pi^2)^2 dx
= pi Int (2+e^(ix)+e^(-ix))   (1/2api^2) d/da 1/(x^2-a^2pi^2) dx      @(a=1)
= pi Int (2+e^(ix)+e^(-ix))   (1/2api^2) d/da (1/2api) [1/(x-api) - 1/(x+api)] dx
          ^   ^       ^
Integrate first, differentiate second. 
the 2 term multiplies factors like 1/(x-api) and the principal value integrates to 0.
drop that term.  The next two terms are complex conjugates so the net result is

(1/4api^2) d/da (1/a) Int 2 Re e^(ix)  [1/(x-api) - 1/(x+api)] dx
                                            ^            ^
let x--> x+api in the first term and x--> x-api in the second term, obtain

= 2/(4api^2) d/da (1/a) Re 2i sin(api) Int e^(ix) /x  dx

The integral is Int (cos(x)+isin(x))/x dx. cos part is an odd function, integrates to 0.
Int sin(x)/x dx = pi.  this leaves

(2 2i ipi)/(4api^2) d/da (1/a) sin(api) = (-1/api) d/da (1/a) sin(api)    @ a=1

differentiation of the 1/a factor still leaves sin(api) = 0   @ a=1
differentiation of sine factor gives

(-1/api) (1/a) pi cos(api)   @ a=1    = 1    ok

[Paul]

Hi David,

Matlab had trouble too. If you can go back through the history of this Answer you'll see that I was originally using vpaintegral because I couldn't wrangle Matlab into finding a closed form expression. It wasn't until I forced Matlab to only integrate the real part of the integrand that it could come up with a closed form solution (I'm pretty sure I tried rewriting it in terms of complex exponentials first). But even that result is interesting.

Repeating the code from above.

syms x real
L = 100;
psi(x) = sqrt(sym(2)/L)*sin(sym(pi)*x/L)*rectangularPulse(0,L,x); 
syms w real
Psi(w) = simplify(fourier(psi(x),x,w))

Psi(w) = 

I = Psi(w)*conj(Psi(w))

I = 

I = simplify(real(I))

I = 

So the integrand has a form that at least looks potentially tractable. Maybe. But, I wouldn't want to integrate it by hand. The anti-derivative looks very complicated to me

E(w) = simplify(int(I,w))

E(w) = 

But taking the limits results in the miracle, even with those log functions that (I think) only cancel each other in the limit

limit(E(x),x,inf) - limit(E(x),x,-inf)

ans = 

[David Goodmanson]

Hi Paul,

Although this analytical solution aspect is not as important as the sum total of what you showed in your answer, I have a couple of comments. In this case with the limits of integration at +-inf, it does work to find the indefinite integral as a set of real functions and evauate it at the end points, but I don't think it's preferred. And for those functions, it appears that you might have to take someone's else's word for it on the values at +-inf.

The method I mentioned reduced everyhing to solving [ Int sin(x)/x = pi ]. That's dead easy to prove with complex variables. Even better, one can also go with complex integration of the original integrand from the get-go and deal with second order poles, which gives a shorter and more direct derivation than the one I posted.

Answer 2

[VBBV]

Try with fftshift function ,

L = 100;
x = linspace(0,L,L+1);
psi = sqrt(2/L)*sin(pi*x/L); %particle in a box ground state%

%<<norm is standard builtin function in matlab
Norm = trapz(x,psi.*conj(psi)) %normalized state has L-2 norm = 1%

Norm = 1

F_T = (fftshift(psi));
F_T = norm(F_T)

F_T = 1

[bil]

Hi, thanks for the response! Perhaps this is getting a little technical, but could you explain the "norm2" in my code is wrong (or how exactly the fftshift resolves things)? I was under the impression that for a discretized array in position space (my "x" array), in momentum space, it would become so then I should integrate the wavefunction using that spacing instead.

[Paul]

I thought the point in the Question is to transform to the frequency domain. But I don't see anything in the frequency domain in this Answer?