How do you do a radial Fourier transform in MATLAB?

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?

Answer 1

[David Goodmanson]

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

[N/A]

Thanks, David. What about multiplying

-imag(fftshift(fft(ifftshift(rFr))))*delr

by q? Also, when I implement your code, exchanging the lines

Frpos = exp(-5*rpos);    % F(r), r >=0
rFrpos = rpos.*Frpos;

with

Frpos = @(rpos) interp1(r,fr,rpos,'spline');    % F(r), r >=0
rFrpos = rpos.*(Frpos(rpos)-1);

because 1) my function is a vector and 2) I subtract 1 from it, I get a function f(Q) that is symmetric about the x-axis, which is not at all what I expect. Am I misunderstanding something? Let me know if it would help to attach example data. Thanks again.

[David Goodmanson]

hey S,

The transform is between the obfect rFr = r*F(r) and the object qFq = q*F(q). So qFq already contains the multiplicative factor of q. If you want F(q) itself, you can just divide qFq by q.

I forgot to mention it, but I was treating f(r)-1 as a new function F(r), which it looks like you are doing above. However you want to accomplish it, you need to define a function rFrpos on the domain 0<=r<=r0, and then create an antisymmetric function rFr on the domain -r0<=r<=r0. The code does that with

rFr = [-fliplr(rFrpos(2:end)) rFrpos]; 

Once rFr is an antisymmetric function, and with the use of fftshift and ifftshift, qFq has to come out as an antisymmetric function.

[N/A]

Thanks, David. I'm still running into difficulty, however. For 'rFr' as I defined it with

Frpos = @(rpos) interp1(r,fr,rpos,'spline');    % F(r), r >=0
rFrpos = rpos.*(Frpos(rpos)-1);
rFr = [-fliplr(rFrpos(2:end)) rFrpos];

I get the following antisymmetric function:

But for 'qFq' I get something symmetric:

Do you know what the issue might be?

I really appreciate your help.

-Sam

[David Goodmanson]

what kind of plot commands are you doing, exactly? If on the plots the horizontal variable is x, then the idea is that f(-x) = -f(x) so that the function is antisymmetric about x = 0. But neither plot extends into negative values of x so there is no way to say whether either function is antisymmetric.

The following plot is antisymmetric about x=0:

x = -10:.01:10;
xpos = x(x>=0);
ypos = xpos.*exp(-xpos);
y = [-fliplr(ypos(2:end)) ypos];
figure(1)
plot(x,y)
grid on

as opposed to this one which is symmetric about x = 0

ypos = xpos.^4.*exp(-xpos);
y = [fliplr(ypos(2:end)) ypos];   % no minus sign
figure(2)
plot(x,y)
grid on

[N/A]

Thanks, David. For the plot shown above of 'rFr' I am using:

plot(rpos,rFr(rpos>=0))

where

Frpos = @(rpos) interp1(r,fr,rpos,'spline');    % F(r), r >=0
rFrpos = rpos.*(Frpos(rpos)-1);
rFr = [-fliplr(rFrpos(2:end)) rFrpos];

Note the '-1' in 'rFrpos'.

For the plot shown above of 'qFq' I am using:

plot(q(q>0),qFq(q>0))

So I am confused about why the plot of 'qFq' has the form that it does.

[David Goodmanson]

Depending on what rFr is, the plot of qFq may well be correct. It looks like when you zoom in there will be a lot of oscillations. What most matters are the plots for positive r and positive q that you have, because those are the physcically significant ranges. Not knowing what your fr is and the number of points there are, I can't reproduce the plots.

[JH]

Thank you for the very helpfull code @David Goodmanson.

It's seems to work very well for cases where is monotonically decreasing. E.g. if then transform is given by and the numerical transform matches the analytical one:

However for decreasing but oscillating function, e.g. the numerical and analytical results don't seem to match (case j=2 in the MWE)

Probably there some silly error I'm making here regarding how q should be mapped in the numerical and analytical case(?)

My best

---------------

Here's MWE based on the previous code

clear
% create a function for positive r
n = 2e2;
delr = 1;
r = (-n:n)*delr;
rpos = r(r>=0);
N = length(r);
delq = 2*pi/(N*delr);   % golden rule for fft
q = (-n:n)*delq;

% function parameters
a = 20*delr; 
b = 1e-1/delr;

% test cases (ground truths)
j = 2;

% real space functions and their analytical Fourier transforms
if j == 1
    Frpos = exp(-rpos/a); % F(r), r >=0 
    Fq_theor = 8*pi*a^3./((1 + a^2*q.^2).^2); 
elseif j == 2 
    Frpos = exp(-rpos/a).*sinc(rpos*b); % F(r), r >=0
    Fq_theor = 8*pi*a^3./(1 + a^4*(b^2 - q.^2).^2 + 2*a^2*(b^2 + q.^2));
end
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 
qFq = -imag(fftshift(fft(ifftshift(rFr))))*delr;
figure(1);clf; 

subplot(1,3,1)
    plot(rpos,Frpos)
    legend('f(r>0)')
    legend('Location','northoutside')
subplot(1,3,2)
    plot(q,qFq./q,'b-',q,Fq_theor/(2*pi),'r--')
    grid on; axis tight
    xlim([-1 1])
    legend('Numeric F(q)','Theoretical F(q)')
    legend('Location','northoutside')

% 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;
subplot(1,3,3)
    plot(r,rFr./r,'--',r,rFr2./r,'o')
    grid on; axis tight
    xlim([-50 50])    
        legend('original f(r)','Back transfered f(r)')
        legend('Location','northoutside')

[David Goodmanson]

Hi JH

The analytic expression I came up with is

Fq_theor = 8*pi*a^3./((1 + a^2*(b+q).^2).*(1 + a^2*(b-q).^2));

but as you pointed out below this is equivalent to the expression you have.

To make the traces in the middle plot above agree, you have to make a change to the sinc function. You have sin(br)/br but unfortunately the Matlab convention for sinc is not sin(x)/x but rather the other fairly common one, sin(pi x)/(pi x). I would have preferred the first expression myself but it was not to be. Rather than having to think about whether to multiply or divide the argument by pi when using the Matlab function I just made another one which you can put onto the bottom of the script you have, or perhaps save it separately.

function y = sinc1(x)
% sinc function, sin(x)/(x), no checks on the input
%
% y = sinc1(x)
y = sin(x)./x;
y(x==0) = 1;
end

[JH]

Much obliged, David! =)

You reply is very helpful and appreciated :). I'll probably be bit lazy and use an anonymous function based on your suggestion

sinc2 = @(x) sinc(x/pi);

I think to match the analytical expression with the numerical results then, one has to divide the former with due to different conventions with FFT(?)

% create a function for positive r
n = 1e3;
delr = 1;
r = (-n:n)*delr;
rpos = r(r>=0);

% sinc(x) in Matlab (and numpy) is sin(pi x)/(pi x)
sinc2 = @(x) sinc(x/pi);

a = 50*delr; 
b = 3e-1/delr;

N = length(r);
delq = 2*pi/(N*delr);   % golden rule for fft
q = (-n:n)*delq;

% f(r) and analytical F(q)
Frpos = exp(-rpos/a).*sinc2(rpos*b); % F(r), r >=0    
Fq_theor = 1/(2*pi)*(8*pi*a^3./((1 + a^2*(b+q).^2).*(1 + a^2*(b-q).^2)));

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 
qFq = -imag(fftshift(fft(ifftshift(rFr))))*delr;

figure(1);clf; pause(0.01)

% analytical solution to the fourier transform is 
subplot(1,3,1)
    plot(rpos,Frpos)
    legend('f(r>0)')
    legend('Location','northoutside')

subplot(1,3,2)
    plot(q,qFq./q,'b-',q,Fq_theor,'r--')
    grid on; axis tight
    xlim([-1 1])
    legend('Numeric F(q)','Theoretical F(q)')
    legend('Location','northoutside')

% 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;
subplot(1,3,3)

% visualization discarding step size (adaptive)
    s = round(length(r)/500);
    plot(r,rFr./r,'--',r(1:s:end),rFr2(1:s:end)./r(1:s:end),'o')
    grid on; axis tight
    xlim([-50 50])    
        legend('original f(r)','Back transfered f(r)')
        legend('Location','northoutside')

And the results seem to match perfectly - many thanks! =)

P.S: I forgot this thread for few days, and then also found out in another context that sinc(x) has different convention in Matlab (same as issue in numpy).

I think expressions are equivalent (using Mathematica shamelessly here)

[David Goodmanson]

Hi JH, I see they're equivalent and modified my comment accordingly.

Answer 2

[N/A]

Thanks, David. I think I'm getting closer to resolving this. Can you please explain the following?

delq = 2*pi/(N*delr);   % golden rule for fft
q = (-n:n)*delq;

Also, when I change n to 17000 (from 2e4) and delR to 0.00001, it radically alters the result. Can you please explain what the optimal choice should be? The function f(r) only goes out to 17.95, which is why I changed the value of n.

[David Goodmanson]

again, not knowing what your f(r) looks like makes it hard to comment.

[N/A]

Below is an example f(r). As mentioned, it goes out to only about 18, and I am expecting f(Q) to be interesting only in about the same range (out to about 18 or so in Q-space). So I think the issue is making sure there are enough data points in the x-axis of the FFT over this range. That's why I asked my questions above about the x-axis for both f(r) and f(Q). (When I run the code and zoom into the relevant region, it does seem ok, but there are far too few data points to be able to say for sure.)

[David Goodmanson]

This function is very spead out and slow moving in the r window. Nothing wrong with that, but it does mean that in the q window, qFq will be scrunched up against the origin. Keeping array spacing delr the same and extending the r window out further, by a factor of 5 for example (consequently the number of points increases by a factor of 5) will improve the frequency resolution in q [see last paragraph]. In any case you will have to zoom in toward the origin in the q plot as you are doing.

It makes sense that there are oscillations in f(q) since the f(r) function does not get going until you are well away from the origin. That resembles a shift in f(r) which leads to an oscillatory phase factor in f(q).

the golden rule is just a necessary consequence of the fft. In terms of f and t, suppose you have a time array, N points with spacing delt. The time window is T = delt*N. With fourier transform by fft, an even number of oscillations have to occur in the time window. So f1*T = 1, f2*T = 2, f3*T = 3 etc. That means that (f_n+1 - f_n)*T = delf*T = 1. So delf = 1/(N*delt) and delf*delt = 1/N. Since w = 2*pi*f, then delw*delt = 2*pi/N. That's the basic situation with r and q.

Note that if delr*delq = 2*pi/N, then keeping delr the same and increasing the width of the r window by increasing N means that delq gets smaller. So there is Increased frequency resolution.

[N/A]

Thank you, David. Finally, can you just please comment on why you used 'iffshift' and 'ffshift' here, as well as why you need to multiply the FFT by 'delR'? Naively, I would think that that multiplying it would be unnecessary here as it would be incorporated into the FFT itself. -Sam

[David Goodmanson]

The fft operates on functions with an array index. It doesn't know anything about what the distance between array points is, so to create an integral you have to supply that yourself.

The r and q arrays have zero at the center, but the fft wants zero as the first element. fftshift and ifftshift accomplish that.