Partial Fourier transform of large sparse 3D array

8 views (last 30 days)
I am trying to calculate the Fourier transform of a large 3D array (30,000 x 30,000 x 1500 pixels) with about 30,000 x 15,000 nonzero elements. There is at most 1 nonzero pixel in each x-y plane (first 2 dimensions). I cannot downsample the array because I need enough input resolution to get a specific data span in the output.
To avoid storing/loading multiple 10TB arrays, I'm thinking of doing 1D Fourier transform on each dimensions sequentially. I would appreciate further suggestions to speed up computation and reduce memory load. Moreover, if I only want to compute the first 256 x 256 x 768 pixels of the output, is there a way to further reduce computation?
  1 Comment
David Goodmanson
David Goodmanson on 3 Dec 2023
Hi Linh,
you say that there is at most 1 nonzero pixel in each x-y plane, but there are only 1500 x-y planes so this appears to conflict with the statement that there are 30,000 x 15,000 nonzero elements.

Sign in to comment.

Answers (2)

Matt J
Matt J on 3 Dec 2023
Edited: Matt J on 3 Dec 2023
You can do the partial Fourier transforms by computing a set of partial DFT matrices:
dftPartial=@(N,k)sparse(exp(-(2j*pi/N)*(0:N-1).*(0:k-1)'));
Fxy=dftPartial(3e4,256);
Fz= dftPartial(1500,768);
and combine them into a 3D operator by downloading the KronProd class,
F=KronProd({Fxy,Fz},[1,1,2])
The 3D sparse array, X, that is being transformed can be represented as an ndSparse array,
Then, you shoul be able to do simply,
Y=full(F*X);
  1 Comment
Matt J
Matt J on 3 Dec 2023
This example took about 2 min. on my 16 GB laptop
>> X=ndSparse.sprand(F.domainsizes,1500/size(F,2));
>> tic; Y=full(F*X); toc
Elapsed time is 126.221071 seconds.
>> whos F X Y
Name Size Kilobytes Class Attributes
F 50331648x1350000000000 207247 KronProd
X 30000x30000x1500 36 ndSparse sparse
Y 256x256x768 786432 double complex

Sign in to comment.


Matt J
Matt J on 3 Dec 2023
Edited: Matt J on 3 Dec 2023
There is at most 1 nonzero pixel in each x-y plane (first 2 dimensions) ... I'm thinking of doing 1D Fourier transform on each dimensions sequentially.
That would seem to be more work than necessary. The 2D DFT of an impulse (and any sub-rectangle of it) is something that can be computed anaytically. So, if each slice is an impulse image, it's really only the FFT along the 3rd dimension that you would need to do numerically, which should be very tractable.

Tags

Community Treasure Hunt

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

Start Hunting!