Is it possible to correctly perform a multi-dimensional FFT on a 1D linearised version of a 3D array?

I have a 3D array indexed by A_3D(ind_y,ind_x,ind_z). This has then been converted to a single 1D linear vector B_1D(ind) such that the coordinates cycle through the z coorindinates fastest, followed by y, and finally by x which is the slowest to increment.

I am able to take the FFT of the 3D matrix A_3D no problem, using fftn(). I would like to know if it is possible in Matlab to successfully take a FFT of B_1D directly. In other words, is it possible to do a multi-dimensional FFT when working with a linearised 3D array?

My attempt below fails, because when calling fft() it simply takes a 1D FFT of the entire vector, not caring that there are blocks of points which correspond to different dimensions:

Nx = 3;
Ny = 2;
Nz = 4;
N = Nx*Ny*Nz;

[X , Y , Z] = meshgrid(1:Nx,1:Ny,1:Nz);

%%% Do the FFT on a 3D matrix, and reshape the output to 1D vector
A_3D = rand(size(X));   % Create random 3D matrix, indexed by A_3D(indy,indx,indz)
A_3D_FFT = fftn(A_3D);  % 3D FFT
A_1D_FFT = reshape( permute(A_3D_FFT,[3,1,2]),[1 N] ); % Reshape output of FFT so cycle through z fastest, then y, then x

%%% Do FFT correctly on a reshaped, 1D version of the original 3D array?
B_1D = reshape( permute(A_3D,[3,1,2]),[1 N] ); % Instead of reshaping output of FFT, reshape original matrix instead
B_1D_FFT = fftn(B_1D, [Ny,Nx,Nz]); % Is there a way to do this, even though B_1D has been flattened?

max( abs(A_1D_FFT - B_1D_FFT) ) % Check if equivalent

I know that this is possible in C++, because in the FFTW library help docs you can specify a "FFT plan" for fftwnd using

my_plan = fftw3d_create_plan_specific(Nx,Ny,Nz,...

which gives information about how the initial 3D matrix was unwrapped to 1D. Then you can do something like

B_1D_FFT = fftwnd(my_plan , 1 , B_1D , 1 , Nx*Ny*Nz , ... );

and in this case the transform is done correctly. I would like to do the same thing in Matlab (which I believe is based on the same fftw library)

[Thomas Barrett]

Can anyone help with this?

[Matt J]

I don't believe there is a way to do it with the same efficiency as in 3D form.

Regardless, I think you should explain the advantage you hope to gain by doing this. There is virtually no computational or memory cost to reshaping your 1D vector back to 3D for the purposes of doing the fft.

[Thomas Barrett]

Hi @Matt J. In the following example on my system, the time spent doing the 1D-3D and 3D-1D reshaping is about 44% of the total time (about the same as doing the fft itself)

Nx = 256;
Ny = 256;
Nz = 128;
N = Nx*Ny*Nz;
A_1D = rand(N,1);

tic
for k = 1:20
    A_3D = permute( reshape( A_1D, [Nz,Ny,Nx] ) , [2,3,1] );
    A_3D = fftn(A_3D, [Ny,Nx,Nz]);
    A_1D = reshape( permute(A_3D,[3,1,2]) , [N 1] );
end
toc

Is this what you mean?

[Matt J]

That's almost what I mean. It's not the reshape() that is costing you time. It is the permute(). There is no need to permute that I can see. You could just as easily do,

Nx = 256;
Ny = 256;
Nz = 128;
N = Nx*Ny*Nz;
A0 = rand(N,1);


tic
for k = 1:20
    B = reshape( A0, [Nz,Ny,Nx] ) ;
    A_3D = fftn(B);
    A_1D = reshape( A_3D, N,1);
end
toc

Elapsed time is 3.544917 seconds.

and we can see below that the time spent is essentially the same time as if only the fftn() operation is included:

B= reshape( A0, [Nz,Ny,Nx] ) ;
tic
for k = 1:20    
    A_3D = fftn(B); 
end
toc

Elapsed time is 3.483294 seconds.

[Thomas Barrett]

Great, thanks @Matt J you are totally right - there was no need to permute, since it was being undone again. I didn't realise that permute was so intensive, and that reshape was so efficient.

I was originally interested in whether Matlab allows access to some of the other fftw features (such as "plans", etc), but doing the reshape is sufficient. Do you want to add it as an answer, then I can accept it for completeness?

Answer 1

[Matt J]

Reshaping to and from 3D format should not add significant cost:

Nx = 256;
Ny = 256;
Nz = 128;
N = Nx*Ny*Nz;
A0 = rand(N,1);

tic
for k = 1:20
    B = reshape( A0, [Nz,Ny,Nx] ) ;
    A_3D = fftn(B);
    A_1D = reshape( A_3D, N,1); %Strangely, A_1D=A_3D(:) is much slower ????
end
toc

Elapsed time is 3.575770 seconds.

When compared to processing without reshaping, it is almost the same:

B= reshape( A0, [Nz,Ny,Nx] ) ;
tic
for k = 1:20    
    A_3D = fftn(B); 
end
toc

Elapsed time is 3.422608 seconds.

Answer 2

[rajmouli jujjavarapu]

F_d = fft(t_d, [], n);

At the n term if you could specify the dimension, in your case it is 3. I guess this computes each layer individually.

[Thomas Barrett]

This is not correct. It is not as simple as that.

Specifying n as you say is only when t_d is multi-dimensional, in which case n then tells fft along which dimension to operate.

Instead, my vector is 1D, and I need to somehow give the information of Nx, Ny, Nz sizes to fft.

[rajmouli jujjavarapu]

So if you reshaped the multi-dimentional vector to a signal dimention and compute FFT, if this is your question then I dont think I know of any function that can do it. As far as I see it, for the fft function it will be a new matrix with single dimention.

You can try to make your own fuction, which can reshape it back to the multidimension vector and compute FFT

[Thomas Barrett]

For sure the output of the FFT of the 1D vector will be another 1D vector - this I am happy with. But it should not simply be a direct 1D FFT of that vector (it should somehow account for the fact that this is a linearised 3D matrix). And yes, it can certainly be done by reshaping before using fftn and then reshape back after, but this is computationally costly, and I would like to avoid it.

[rajmouli jujjavarapu]

I understand. I didnot find any function for this, but lets wait and see. somebody will have a good idea and will reply.

Thank you.