3D Fourier transform in cylindrical coordinate

I have a numerical function, f(rho,phi,z), in a 50x50x50 array, and trying to get ft(k_rho,k_phi,k_z).

What I got about the Fourier transform in cylindrical coordinates is:

f(k_rho,k_phi,k_z) = int_0^{infty} drho int_0^{2pi} dphi int_{-infty}^{infty} dz f(rho,phi,z)*exp(-i2pik_rho*rho)*cos(k_phi-phi)*exp(-i2pi*k_z*z)

and perform Fourier transform on z-coordinate first, i.e.

f(rho,phi,k_z) = int_{-infty}^{infty} dz f(rho,phi,z)*exp(-i2pi*k_z*z).

and perform Hankel transform for $\rho$ and $\phi$ coordinates.

What I do not understand is how to numerically perform Fourier transfrom in z-coordinate only in Matlab.

I have been trying to use fft function is Matlab, but no idea how to apply fft on z-coordinates only. I am pretty sure that neither ft=fft(f) nor ft=fft(f,[],3) would give me the right answer.

I might even be wrong mathematically. Could anyone explain what I should be doing? Thank you.