looped 3D FFT of 4D data

16 Ansichten (letzte 30 Tage)
sbr
sbr am 30 Jan. 2024
Kommentiert: William Rose am 1 Feb. 2024
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!

Antworten (2)

William Rose
William Rose am 30 Jan. 2024
Bearbeitet: William Rose am 30 Jan. 2024
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.
  4 Kommentare
William Rose
William Rose am 31 Jan. 2024
Bearbeitet: William Rose am 31 Jan. 2024
[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
William Rose am 1 Feb. 2024
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.

Melden Sie sich an, um zu kommentieren.


sbr
sbr am 31 Jan. 2024
"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.
  1 Kommentar
William Rose
William Rose am 31 Jan. 2024
@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.

Melden Sie sich an, um zu kommentieren.

Produkte


Version

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by