3D Fourier Series fftn?

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?

Answer 1

[David Goodmanson]

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.

[sbr]

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]

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.