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.
online but it doesn't seem to resolve the issue either. Any advice would be appreciated.Thanks!
采纳的回答
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)
Hmm, int can't handle that with the inf limits. Recompute with finite limits
Ecs = int(psi(x)*conj(psi(x)),x,-1e10,1e10)
Its Continuous Fourier Transform (CFT)
syms w real
Psi(w) = simplify(fourier(psi(x),x,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)
Can't find a closed form expression.
Verify the imaginary part of the integrand is 0
imag(Psi(w)*conj(Psi(w)))
Integrate just the real part
Ecf = int(simplify(real(Psi(w)*conj(Psi(w)))),w,-inf,inf)/2/sym(pi)
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)
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%
The energy in the discrete space samples is
Eds = sum(psi.*conj(psi))
The continuous space energy is approximated by scaling by dx
Ecapproximate = Eds*dx
"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
Approximate Ec after scaling by dx
Ecapproximate = Edtft*dx
DFT of the samples of psi
dft = fftshift(fft(psi));
Discrete energy in the frequency domain
Edft = sum(dft.*conj(dft))/N
Mulitply by dx to approximate Ec
Ecapproximate = Edft*dx
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
Add DFT to plot
stem(w,abs(dft)*dx)
xlim([-0.5 0.5])
legend('CFT','DTFT','DFT')
7 个评论
bil
2023-4-1
bil
2023-4-2
Paul
2023-4-2
I'm afraid I don't. Did you get Simpson's rule to work in Matlab?
bil
2023-4-2
David Goodmanson
2023-4-4
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
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))
I = Psi(w)*conj(Psi(w))
I = simplify(real(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))
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)
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.
更多回答（1 个）
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%
F_T = (fftshift(psi));
F_T = norm(F_T)
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.
so then I should integrate the wavefunction using that spacing instead.
Paul
2023-4-1
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?
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
