looped 3D FFT of 4D data

I try setting up a fft of 3d meassurement data I obtained over some timesteps t, which then combined gives 4D data.

For that, I used matlab's fftn command, but my problem here is, that I get zero amplitude over all timesteps for the fft analysis.

I followed some multidimensional fft posts I found here on the forum to build my code, nevertheless it seems to have an error inside.

%P: (x*y*z*t)-matrix. P(:,:,:,1) == 3D-matrix of equally spaced gridpoints

% and their meass. value at timestep t=t1

load P.mat

%https://de.mathworks.com/matlabcentral/answers/13896-fftshift-of-an-image

%recorded datapoints per dimension

[Ny,Nx,Nz,L] = size(PQR_4D{1}); % pixels; Nx,y,z == samples_x,y,z;

% L = 3D meassurement data obtained on L evenly spaced timesteps

%dimensions of cube

Lx = 2*30e-3; Ly = 2*5e-3; Lz = 2*90e-3; % meters

%pixel sampling rates for X,Y,Z direction

Fs_x = Nx/Lx; % pixels per meter

Fs_y = Ny/Ly;

Fs_z = Nz/Lz;

%compute the pixel sizes in each direction as the reciprocal of the pixel sampling rates

dx = 1/Fs_x; % meters per pixel

dy = 1/Fs_y;

dz = 1/Fs_z;

%compute a linear distance scale for the X,Y,Z-axes of the cropped image

x = dx*(0:Nx-1)'; % meters

y = dy*(0:Ny-1)';

z = dz*(0:Nz-1)';

%compute the frequency increments for X,Y,Z

dFx = Fs_x/Nx; % cycles per centimeter

dFy = Fs_y/Ny;

dFz = Fs_z/Nz;

%create a frequency domain array for X,Y,Z

Fx = (-Fs_x/2:dFx:Fs_x/2-dFx)'; % cycles per centimeter

Fy = (-Fs_y/2:dFy:Fs_y/2-dFy)';

Fz = (-Fs_z/2:dFz:Fs_z/2-dFz)';

figure;

title("Single-Sided Amplitude Spectrum of X(t)");

xlabel("f (Hz)"); ylabel("|P1(f)|"); hold on; grid on;

for timestep=1:size(P,4)

Y = fftshift( fftn(P(:,:,:,1))/(Nx*Ny*Nz) );

% 1/N corrresponds to 1/(2*L) in continuous case

P2 = abs(Y);

P1 = P2(1:fix(Nx/2)+1);

P1(2:end-1) = 2*P1(2:end-1);

f = Fs_x/Nx*(0:(Nx/2));

plot(f,P1,"LineWidth",3);

legend();

end

Lx, Ly, Lz (unit: meters) define a 3D cube in which voltage-data is meassured on each gridpoint. This 3D data is sampled 'L'-times, with each sample beeing taken 5.5e-5 seconds apart. I also include my datafile P.mat, but it only contains 4 out of 60 timesteps since it would be too big.

I would appreiciate if someone could show me my fault in the code, thanks anyway!

Answer 1

[William Rose]

@sbr,

P2 is 64x64x64. P1 is 1x33. Your plot shows that one column along one corner-edge of the cube P2 is zero. But P2 is not all zeros. For example, after running your code (I had to make a couple of changes for it to run), I plotted a few slices through the cube:

figure;
for i=1:20, plot(1:64,P2(:,3*i,10),'-r.'); hold on; end

and I got the result below:

Since I am not famliar with your data, I don't know what this means, but it is not all zeros.

[William Rose]

@sbr,

Here's another way top look at the P2 data.

Make 64 slices through the 64x64x64 cube. Each slice is 64x64. Normalize each slice by its maximum value, so that the max pixel value in each slice is one. Then arrange the 64 slices into an 8x8 grid and show the result. I already have P2 from running your code, when I do the following

P2nn=zeros(64,64,64);

for k=1:64,P2nn(:,:,k)=P2(:,:,k)/max(max(P2(:,:,k))); end

out2=imtile(P2nn,GridSize=[8,8]);

imshow(out2)

And I get this:

You could try slicing it along the other directions.

[sbr]

@William Rose,

Thank you for your answers!

"I had to make a couple of changes for it to run": Were there other bugs than the fauly plot or is the FFT in general set up correctly?

"Since I am not famliar with your data, I don't know what this means, [...] ": I wanna try and explain the dataset a bit more. In the original post I said it is meassured voltage values in a cube, but that was maybe simplified too much. In fact, 'P' is a 4D matrix containing the derivatives of a magnetic field. One can access the first timestep P(t=t1) with 'P(:,:,:,1)'. The 4th dimension just sweeps through all timesteps. And then each entry of this 3D matrix just contains the rate of change of the magnetic field at the respective 3D grid at time t. The grid is evenly spaced with

[x_mag,y_mag,z_mag] = meshgrid(linspace(0,Lx,Nx),linspace(0,Ly,Ny),linspace(0,Lz,Nz));

The data is known to be periodic, that's the basic motivation to make an fft: "[...] Hence the sources P can be observed to be repeating itself at every Lx, Ly, Lz. This allows to compute the source frequency components within a cube."

"P2 is 64x64x64. P1 is 1x33": Seems like I extracted the single-sided spectrum as if my data was 1D

[William Rose]

@sbr,

[edit: Change "90 degree phase change" to "180 degree phase change". Add comment about the direction of X, Y in the images.]

You're welcome.

"Were there other bugs..." I had to change a variable name. See attached code. I also changed the plot so that the plot of the one-sided FFTs would be separate for each of the 4 times in the file. I made the plots other plots by running commands at the Matlab prompt, as posted above.

"Seems like I extracted the single-sided spectrum as if my data was 1D." True.

"Hence the sources P can be observed to be repeating itself at every Lx, Ly, Lz. This allows to compute the source frequency components within a cube." Do you mean by this that the spatial period of the oscillating field is Lx,Ly,Lz? If so, then you only get one period of the oscillation, in each dimension, in your sample, since your sample dimensions are Lx by Ly by Lz. If the oscillation were sinusoidal with those frequencies or periods, then you;d expect the FFT to be zero except at the very lowest frequency in each dimension.

Let's examine the raw data.

load P
P1=P(:,:,:,1);;         % P at time 1
Pn1=P1/max(P1,[],'all');  % normalize P by its overall max
out1=imtile(Pn1,GridSize=[8,8]);
figure; imshow(out1); 
title('Slices of P1, perp to z, norm. by overall max');

Each of the 64 images in the montage above is a slice showing the X-Y plane, at a particular Z-value. The z-axis position increases across the columns of row 1, then across the column of row 2, etc. The first 32 images are roughly identical. Then there is an abrupt 180 degree phase change, which is maintained for the remaining 32 slices. Below are enlarged images of slice 32 and slice 33, made with the command figure; imshow(P1(:,:,32)); title('k=32'), etc. In the images below, X increases down,and Y increases from left to right.

This is not what I expected, based on your description of the data. Is this what you expect for dB/dx, or whatever is being measured?

This reminds me of a summer project I had, working on a particle accelerator at the Los Alamos National Laboratory: map the magnetic field strength in 3 dimensions between the poles of a beam-bending magnet. The scale was a bit different than yours; the current in the magnet was a few hundred amps, and the separation between measuring points was on the order of cm, not um.

I realize the montage plot above will be quite small when posted online, because Matlab Answers is shrinking the graphical output, for some users' posts. That has been discussed elsewhere, and I think people at Mathworks are lookng into it. If you run the code yourself, I think you will see an image that is not so tiny.

[William Rose]

@sbr,

It appears that a comment I made earlier did not post, so I am trying again. Sorry if it's a duplicate.

Views of slices of the raw data (not the FFT): X-Y, Y-Z, and Z-X slices. ThisThe other dimension, within each montage, inreases along the columns of row 1, then along the columns of row 2, etc.

These views show that something changes abruptly, half way along each axis.

Answer 2

[sbr]

"Do you mean by this that the spatial period of the oscillating field is Lx,Ly,Lz?": Yes, at least that's what I read out of the paper my script refers to: "The consecutive reflections of the sequence form a 2-D pattern that has periodicities of Lx and Ly in the x- and y-directions, respectively." That later also applies for the 3rd dimension. Do you have IEEE access by any chance? Then I could share the paper with you.

"This is not what I expected, based on your description of the data. Is this what you expect for dB/dx, or whatever is being measured: Absolutely true, there must be a sign-bug somewhere, I am gonna look that up!

"This reminds me of a summer project I had, working on a particle accelerator at the Los Alamos National Laboratory: map the magnetic field strength in 3 dimensions between the poles of a beam-bending magnet.": Yeah, I am onto something similar. My problem is to calculate eddy-current losses in a Permanent Magnet based on 2D flux density data over an electrical period. And the paper I am currently trying to understand does that with some fancy Fourier expansions.

[William Rose]

@sbr,

I love how you have lots of clear comments in your code, and you state the units frequently. Everbody should follow your example.

I have not kept up my IEEE membership, but I can access IEEE journals via my institution.

It is easy to mix up dimensions when working with multidimensional arrays, especially if the array dimesions are the same. Therefore, whenever possible, I make the dimensions unequal in the different directions, to avoid the possibility that dimensions get switched without me knowing it. If you must use 64x64x64, then experiment with arrays with unequal sizes to make sure things are working as you expect. Here is a simple example of what can go wrong.

x=.02:.02:1; y=x;
z=(1-cos(2*pi*(0.5*x)'))*(1-cos(2*pi*(3*y)))/4;
surf(x,y,z,'EdgeColor','none'); grid on;
xlabel('X'); ylabel('Y'); view(45,35);

Things look fine, but they're not, and there's no error, so I woudn't know. In this case, since I made the data, I know the frequency should be one half cycle along X, and 3 cycles along Y. The plot is the opposite. What I need to do is plot the transpose of z:

figure; surf(x,y,z','EdgeColor','none'); grid on;
xlabel('X'); ylabel('Y'); view(45,35);

That is correct. It demonstrates that plotting errors can happen without one realizing it, if the array dimensions are equal.

This is getting complicated. You may email me securely by clicking on the envelope icon that appears in the pop-up when you click on the WR next to my posts.