Using FFT in for-loop is extremely slow - how to accelerate?

Below is the code that I have set up to run FFT on img_norm, a 4D data set. This code takes nearly 1 hour to run partly because I need "n" in the fft function to be 10,000 for my analysis. Any suggestions or thoughts on how I can speed up this calculation would be greatly appreciated! Thank you.

xdim = 128;
ydim = 44;
zdim = 25;
newspacing = 10000;

ft = zeros(xdim, ydim, zdim, newspacing);
ftshift = zeros(xdim, ydim, zdim, newspacing);
ftshiftmag = zeros(xdim, ydim, zdim, newspacing);
phase = zeros(xdim, ydim, zdim, newspacing);

for x = 1:xdim
    for y = 1:ydim
        for z = 1:zdim
             if mask(x,y,z) == 1
            ft(x,y,z,:) = fft(squeeze(img_norm(x,y,z,:)),newspacing); %calculate fourier transform of time series in each voxel
            ftshift(x,y,z,:) = fftshift(ft(x,y,z,:)); %shift fourier transorm of each time series to be centered around 0
            ftshiftmag(x,y,z,:) = abs(ftshift(x,y,z,:)); %get the magnitude of the shifted fourier transform
            phase(x,y,z,:) = angle(ftshift(x,y,z,:)); %calculate the phase angle
             end
        end
    end
end

[Matt J]

Your results should be consuming 60 GB in double floating point. Do you really have that much RAM?

Answer 1

[Matt J]

To conserve memory, I only store the results for x,y,z inside the mask. It would be straightforward to re-embed them in 4D arrays if you really have sufficient RAM.

 I=img_norm(mask(:),:);
 
   ft=fft(I,newspacing,2);
   ftshift=fftshift(ft,2);
   ftshiftmag=abs(ftshift);
   phase=angle(ftshift);
   

[Apoorva Ayyagari]

Thank you, Matt!