Normalizing an FFT Vector

Question

bil 2023-4-1

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1939259-normalizing-an-fft-vector

评论： David Goodmanson 2023-4-5

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2023-4-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1939259-normalizing-an-fft-vector#answer_1206494

编辑：Paul 2023-4-2

在 MATLAB Online 中打开

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')

7 个评论
显示 5更早的评论隐藏 5更早的评论

David Goodmanson 2023-4-4

在 MATLAB Online 中打开

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 2023-4-4

在 MATLAB Online 中打开

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 2023-4-5

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 2023-4-1

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1939259-normalizing-an-fft-vector#answer_1206449

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

bil 2023-4-1

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.