MATLAB: How do you do a radial Fourier transform in MATLAB

fftMATLAB

I have 'r' and a function 'f(r)' as vectors of numbers, with r ranging from 0.05 to 17.95. I want to calculate the radial Fourier transform:
which is the equation that I'm trying to calculate, except for some normalizing constants. After loading 'r' and 'fr', the latter being f(r), here's what I try so far:
Q = linspace(0.01,17.95,1000)';
R = linspace(0.01,17.95,1000)';
fR = @(R) interp1(r,fr,R,'spline');
integrand_fft = @(R,Q) (R.^2).*(fR(R)-1)./(Q.*R);
fQ_fft_v = -imag(fft((R.^2).*integrand_fft(R,Q)));
but the result appears symmetric and does not reach and asymptotic value that I expect, so I think that I did something wrong. How do I set up this code correctly?

Best Answer

Hi S,
here is one way. First of all, the transform pair is
q*F(q) =            Int r*F(r) sin(q*r) dr
r*F(r) = (1/(2*pi))*Int q*F(q) sin(q*r) dq
so the basic functions are qFq = q*F(q) and rFr = r*F(r). If you want F(r) itself you can always obtain it from rFr ./ r and similarly with q. The code below takes r*F(r) for the domain (0,rmax) and for mathematical convenience creates an odd (antisymmetric) function on the domain (-rmax,rmax). That means the sine terms in the fft will survive because they are odd, and the cosine terms will not because they are even. The fft produces an antisymmetric function of q in the domain (-qmax,qmax). For a final result you can just ignore the functions for r<0 and q<0.
The code starts with rFr, does a transform to find qFq. Then, to test the inverse transform there is a transform back to find rFr2 and compare it with rFr. It's the same.
% create a function for positive r
n = 2e4;
delr = .001;
r = (-n:n)*delr;
rpos = r(r>=0);
Frpos = exp(-5*rpos);    % F(r), r >=0
rFrpos = rpos.*Frpos;
% create antisymmetric function with r=0 at the center
% see plot 2
rFr = [-fliplr(rFrpos(2:end)) rFrpos];
% transform to q, multiply by delr to approximate the Riemmann integral 
N = length(r);
qFq = -imag(fftshift(fft(ifftshift(rFr))))*delr;
delq = 2*pi/(N*delr);   % golden rule for fft
q = (-n:n)*delq;
figure(1)
plot(q,qFq)
grid on
% transform back to r, multiply by delq for Riemann integral
% factor of N because the ifft divides by N and that must be canceled out
rFr2 = (1/(2*pi))*N*imag(fftshift(ifft(ifftshift(qFq))))*delq;
figure(2)
plot(r,rFr,r,rFr2)
grid on
Related Question