How do you do a radial Fourier transform in MATLAB?

조회 수: 29 (최근 30일)
N/A
N/A 2020년 5월 7일
편집: David Goodmanson 2023년 2월 3일
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?

채택된 답변

David Goodmanson
David Goodmanson 2020년 5월 7일
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
  댓글 수: 10
JH
JH 2023년 2월 2일
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
David Goodmanson 2023년 2월 3일
편집: David Goodmanson 2023년 2월 3일
Hi JH, I see they're equivalent and modified my comment accordingly.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

N/A
N/A 2020년 5월 8일
편집: N/A 2020년 5월 8일
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.
  댓글 수: 5
N/A
N/A 2020년 5월 9일
편집: N/A 2020년 5월 9일
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
David Goodmanson 2020년 5월 10일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Frequency Transformations에 대해 자세히 알아보기

태그

제품


릴리스

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by