I'm interested in taking some partial derivatives of a 3 dimensional array in Matlab - say $A(i,j,k)$ approximates $f(x_i,y_j,z_k)$. I need to approximate things like $\partial_{xy}f$, $\partial_{yz}f$, etc. I've already built the finite difference matrices for a few of these derivatives, but ultimately this is going to be slow and inaccurate (I'd eventually like to operate on $256^3$ or even $512^3$ arrays).
Does anyone have a good suggestion for computing these partials with FFT? I've tried the following naive thing and it doesn't seem to work:
N=128;
L=2*pi;
x=L/N*(-N/2:N/2-1);y=x;z=x;
[X,Y,Z]=meshgrid(x,y,z);
f=sin(X).*sin(Y).*sin(Z); % A 2-pi periodic function to test
fp_exact=cos(X).*sin(Y).*sin(Z); % Its exact x-partial
ik=1i*(2*pi/L)*[0:nx/2-1 0 -nx/2+1:-1]; % The spectral differentiation vector
fhat=fftn(Z); % Compute the 3-D FFT
C=mat2cell(fhat,nx,ones(1,nx),ones(1,nx)); % So I can act on each column
result=cellfun(@(x) x.*ik',C,'UniformOutput',false); % Multiply each column of Z by ik
fprimehat=cell2mat(result); % Convert back to 3D array
fprime=ifftn(fprimehat); % IFFT
The above code seems hopeless - I think I'm missing something big and haven't really thought about this hard enough. With the above definition of ik, I can compute a 1D spectral derivative with just
fprime=ifft(ik.*fft(f));
Thanks in advance for any tips.