3D Fourier Series fftn?
14 views (last 30 days)
Show older comments
sbr
on 23 Jan 2024
Commented: David Goodmanson
on 26 Jan 2024
I need to reproduce a formula in matlab according to the following screenshot:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1597786/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1597791/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1597796/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1597801/image.png)
Zx = fftn(G3Dx);
amp = abs(Zx);
ang = angle(Zx);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1597806/image.png)
0 Comments
Accepted Answer
David Goodmanson
on 24 Jan 2024
Edited: David Goodmanson
on 24 Jan 2024
Hi sbr
Absolutely this is a 3d fourier series. The reason that the normalization 1/8*LxLyLz does not depend on grid size is that x,y,z are continuous variables, whereas m,n,k are discrete integers. Doing the fft by numerical computation involves discretization of x y and z, (e.g. Nx points in x from 0 to 2*Lx etc.) which is where the grid size that the user must choose will come in. Consider a 1d case:
Lx = 5;
N = 1e5;
n0 = 12;
x = ((0:N-1)/N)*2*Lx; % N points from 0 to 2*Lx
y = cos(2*pi*x*n0/(2*Lx));
delx = (2*Lx)/N;
delf = 1/(N*delx); % fft golden rule; also delf = 1/(2*Lx)
f = (-N/2:N/2-1)*delf;
figure(1); grid on
plot(x,y)
z = fftshift(fft(y))/N; % 1/N corrresponds to 1/(2*L) in continuous case
figure(2); grid on
plot(f,real(z),'o')
xlim([-2 2])
The smallest possible (nonzero) frequency is delf = 1/10. The cosine wave has peaks of amplitude 1/2 at frequencies +-(1/10)*n0 = +- 1.2. To emulate the integral, you multiply by delx to stand in for dx, and divide by 2*Lx as in posted the expression. This gives((2*Lx)/N)/(2*Lx) = 1/N, and that factor gives the correct height peaks in the frequency domain.
2 Comments
David Goodmanson
on 26 Jan 2024
You can divide either on line 1 or line 2, although if you nomalize Zx on line 1 then you don't have to remember to normalize it in the parseval calculation. However you do it, the output of the fft is going to need that factor somewhere along the line.
fftshift is just for convenience.
delf from what I call the golden rule is indeed the frequency step in the frequency grid. Consider a total time T with N points (the fft is circular so N points means N intervals, not N-1, but let's ignore that detail for now), so that the time spacing is
delt = T/N.
The smallest nonzero frequency that can fit into the interval T with an integral number of oscillations has frequency 1/T, so that with cos(2*pi*f*T), the argument is 2*pi*(1/T)*T = 2*pi. Multiples of that frequency have argument 2*pi*n which is correct. So to create the frequency grid, the spacing is
delf = 1/T.
Combining those two you get
delt*delf = 1/N
which is the way I tend to think of it because it can be used to create one grid spacing from the other, and for time and frequency grids coming from some other source it provides a check.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!