3D Fourier Series fftn?

24 次查看（过去 30 天）
sbr
sbr 2024-1-23
评论： David Goodmanson 2024-1-26
I need to reproduce a formula in matlab according to the following screenshot:
is a known matrix, are the 3 dimensions of a 3D object. Does the above equation represent a 3D Fourier Series? If so, I could use the built in fftn-function to decompose the 3D-data into its components via
Zx = fftn(G3Dx);
amp = abs(Zx);
ang = angle(Zx);
seems to me like the normalization term. But shouldn't that depend on the gridsize and not on the physical dimensions?

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

采纳的回答

David Goodmanson
David Goodmanson 2024-1-24
编辑：David Goodmanson 2024-1-24
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 个评论
sbr
sbr 2024-1-24
Thank you for your explanation!
So in the 3D case, if my discretized x,y,z grid has N samples in every direction, the normalization expands to 1/N^3, or in this case 1/(8*N^3). I also tried to confirm this with Parseval's theorem.
Lets say I got a 3D matrix P for the 3D fftn, then I did it as following:
Zx = fftn(P); %line1
a = abs(Zx)./(8*(samples_x^3)); % Magnitude %line2
theta_x = unwrap(angle(Zx)); % Phase %line3
%check Parseval Theorem
time3D = sum(sum(sum(abs(P.^2))));
freq3D = sum(sum(sum(abs(Zx.^2))))./(8*(samples_x^3));
figure;
bar([time3D, freq3D], 'BarWidth', 0.5);
xticks(1:2);
xticklabels({'Time Domain', 'Frequency Domain'});
ylabel('Energy');
title('Energy Comparison');
Energy is equal in time and freq. domain so the normalization seems to be correct now.
But one question still comes up for me: Do I apply the normalization direction when calling the fftn command in line1 and therefore affect the amplitude and frequency components, or apply it only on the magnitude calculation in line2?
I understand the math in all you code lines but am not 100% sure about the theory behind that line:
delf = 1/(N*delx); % fft golden rule; also delf = 1/(2*Lx)
Does it define the smallest frequency step I can evalute with the fft in the chosen grid?
Also, you used fftshift. Is that necessary or is it just for convenient visualisation?
Thank you!
David Goodmanson
David Goodmanson 2024-1-26
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by