How do I interpret the output from fftn? (performed on a 3D matrix: time x spatial dimension 1 x spatial dimension 2)

45 visualizzazioni (ultimi 30 giorni)
I have a 3D dataset, comprised of time-series signals recorded at a range of locations in two spatial dimensions (front/back and side to side). My intention is to measure the amplitude/power of waves travelling in those two spatial dimensions, at each temporal frequency.
To achieve this, I have performed a 3D fast Fourier transform on my matrix (7 rows reflecting front/back spatial locations x 5 columns reflecting left to right spatial locations x 512 timepoints). I have then shifted the zero frequency to the centre of the output, and performed an absolute transform on that output.
This was implemented using the following line of code:
3DFFT=abs(fftshift(fftn(data)));
The output of my analysis shows significant differences between two conditions at specific locations in the two spatial dimensions, within a specific frequency. However, I'm having trouble understanding what the values at each spatial location represent. If I had performed a 2D FFT (with fft2), then my undrstanding is that the output would represent the strength of waves travelling in one direction or the other within the 1 spatial dimension.
However, with 2 spatial dimensions, does the output at each individual location reflect the strength of waves travelling to or from that point? And if so, does that include: A) forwards and backwards travelling waves, B) side to side travelling waves, and also C) diagonal travelling waves?
My confusion arose from running tests on the same data, but with fft2 (performing a 2D FFT) on each of the columns (left to right spatial locations) separately, and noticing the outputs were different from the 3D FFT, prompting me to think that each cell in the fftn output is affected by both the strength of waves travelling front/back and left to right.
Thank you in advance for any help that can be provided!

Risposta accettata

William Rose
William Rose il 28 Nov 2024 alle 7:50
When you compare the 2d FFT to the 3d FFT, you must view slices through the 3d FFT that correspond to the plane of the 2d FFT. You will know you have got the right slice orientation if the dimensions of the slice through the 3d FFT match the dimensions of the 2d FFT. You must also analyze many slices of the 3d FFT to see which slices, if any, look like the 2d FFT.
I understand this stuff better when I analyze simulated data. With only 7x5 spatial points, you are extremely limited in your ability to identify spatial waves, travelling or otherwise. Therefore I use 100 time points and 24x20 spatial points in the three examples below. Examples 1 and 2 show spatial waves that do not travel. Example 3 shows a travelling wave.
% This code is used in all three examples.
% Define the time and space independent variables.
nt=100; nx=24; ny=20; % Assume delta t=delta x=delta y=1
t=0:(nt-1); x=0:(nx-1); y=0:(ny-1); % vectors t, x, y
[T,X,Y]=meshgrid(t,x,y); % 3D arrays for T, X, Y
Example 1: Oscillation along the x-direction only, with spatial frequency fx=1/8:
ft1=0; fx1=1/8; fy1=0; % frequency of oscillation along t, x, y
a=cos(2*pi*ft1*T+2*pi*fx1*X+2*pi*fy1*Y);
A=abs(fftshift(fftn(a)));
The plots of a(time and space) and A (3d FFT of a) are below. Since a and A are 3D arrays, we must plot slices through them. In this example, a(t,x,y) does not change with time. The top panel in the figure below shows a(x,y) at t=0. If we chose a later time, it would look exactly the same.
The bottom left panel shows abs(A(fx,fy)) at ft=0. Remember that the temporal frequency of the oscillation is 0 in this example. The plot has spikes at fx,fy=-1/8, 0 and at fx,fy=+1/8, 0.
The bottom right panel shows abs(A(fx,fy)) at ft=0.04 Hz. There is no power in the 3d signal at this temporal frequency. Therefore abs(A)=0 at all spatial frequencies, for this slice.
Example 2: Oscillation along the x- and y-direction, with spatial frequencies fx=1/8, fy=1/5. The signal is the same at all times, i.e. ft=0:
ft1=0; fx1=1/8; fy1=1/5; % frequency of oscillation along t, x, y
a=cos(2*pi*ft1*T+2*pi*fx1*X+2*pi*fy1*Y);
A=abs(fftshift(fftn(a)));
The plots of a(time and space) and A (3d FFT of a) are below. Since a and A are 3D arrays, we can only plot slices through them. In this example, a(t,x,y) does not change with time. The top panel in the figure below shows a(x,y) at t=0. If we chose a later time, it would look exactly the same. We rotated the view to look down on the surface.
The bottom left panel shows abs(A(fx,fy)) at ft=0. Remember that the temporal frequency of the oscillation is 0 in this example. The plot has spikes at fx,fy=-1/8,-1/5 and at fx,fy=+1/8,+1/5. We rotated the view to look down on the surface.
The bottom right panel shows abs(A(fx,fy)) at ft=0.04 Hz. There is no power in the 3d signal at this temporal frequency. Therefore abs(A)=0 at all spatial frequencies, for this slice.
Example 3: Travelling wave, with spatial frequencies fx1=1/8, fy1=1/5 (same as example 2), and temporal frequency ft1=0.1 Hz.
ft1=0.1; fx1=1/8; fy1=1/5; % frequency of oscillation along t, x, y
a=cos(2*pi*ft1*T+2*pi*fx1*X+2*pi*fy1*Y);
A=abs(fftshift(fftn(a)));
The figure shows plots of a(t,x,y) at three times (top row), and it shows plots of A (3d FFT of a) at four frequencies (bottom row). In some cases, the view is rotated to look down on the surface.
The top row shows a(x,y) at t=0, t=1, t=2. The plots show that the wave is travelling down and to the left. The direction of travel in the x,y plane is given by the vector -(fx1,fy1). The speed is proportional to ft1.
The bottom row shows abs(A(fx,fy)) at ft=0 Hz, ft=+0.1 Hz, ft=-0.1 Hz, and ft=0.2 Hz. Note the differnt vertical axis scales in the four plots. At f=0 Hz and f=0.2 Hz, there is no power in the 3d signal, at any of the spatial frequencies. At f=+-0.1 Hz, there IS power in the signal, at the expected spatial frequencies: fx=+-1/8, fy=+-1/5.
The scripts for the three examples are attached.

Più risposte (0)

Categorie

Scopri di più su Fourier Analysis and Filtering in Help Center e File Exchange

Tag

Prodotti


Release

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by