Partial Fourier transform of large sparse 3D array
6 次查看(过去 30 天)
显示 更早的评论
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 个评论
David Goodmanson
2023-12-3
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.
回答(2 个)
Matt J
2023-12-3
编辑:Matt J
2023-12-3
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 个评论
Matt J
2023-12-3
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
Matt J
2023-12-3
编辑:Matt J
2023-12-3
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.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!