Hi,
This question is regards to calculating the coefficients of a 2D-FFT.
It's a problem that's been bugging me for quite some time and I would
really appreciate some help! I am willing to tip handsomly for a quick
and correct answer. (up to double of price)
The equation is a 2nd kind spherical hankel function of the zeroth order.
It is described here,
http://mathworld.wolfram.com/SphericalHankelFunctionoftheSecondKind.html.
so basically, h(x) = i/x * exp(-ix).
However, here's a few twists, instead of x we have ko*sqrt(Z^2 + C^2),
and from these three terms ko is a constant of about 150. Z is from
0-->a where a is a real number of about 0.03 and C is from 0 --> b
where b is a real number of about 0.02. And because this function has
Z^2 and C^2 terms it is obviously even, and thus the fourier series
would be an expansion of cosines.
In short,
hankel ( ko*sqrt(Z^2 + C^2)) = SUMrs * { hrs * cos (r*pi*Z/a) * cos (s*pi*C/b) }
where hrs are the coefficients and is the usual double integral of the
hankel function multiplied by the two cosine terms. But the double
integral cannot be evaluated in maple in a short period of time so it
is no good, so a sampling algorithm such as the 2D fft must be used.
The resulting coefficient matrix should be 16 x 16, which means r =
0..15 and s = 0..15. Therefore not many samples need to be taken. If
16 x 16 is not possible then another small sized matrix is ok as well.
In conclusion, the coefficients hrs should be computed for say the
values ko = 150, a = 0.03 and b = 0.02. And the way to check would be
to substitute an arbitrary value into the hankel function and
evaluating it and using the fourier series expansion to calculate it
and making sure they are close to 3 digits of accuracy. Also, since
the matlab, maple or mathematica code could be reused for
arbitrary values of ko, a and b that would also be great to have. If
any other methods to calculate the coefficients are better, that's
also adequate!
What I think the process should be:
Forcing the function to be even so that the expansion would be in
cosines as required.
Avoid the singularity at 0. (must)
Taking samples at adequate intervals not exceeding 20 samples.
Taking the 2D fft to compute coefficients.
Check and verify.
I have been stuck on the first 3 :(
The following are what I have so far:
A)
Matlab code to calculated coefficients for 1 variable h(x) and
checking it at x = 3:
(As can be seen, matlab doesn't have a function for the spherical
hankels and I had to convert the normal hankel to the spherical hankel
by the sqx term multiplied by the 0.5th order of the Hankel function
H)
function result = shankeltest (T)
Ns = 25;
first = 0.00000001;
last = T;
amount = (first:(last-first)/Ns:(last-((last-first)/(Ns))))
sqx = sqrt (0.5*pi./amount);
takesample = sqx .* besselh (0.5, 2, amount);
coefficients = fft (takesample)
%result = coefficients;
n = (0:1:Ns-1);
result = 1/Ns * sum(coefficients.*(exp((2*pi*i/T)*n*3)));
B)
Matlab code attempt at making the function even and evaluating the coefficients...
(Note: twiddle is used to avoid the singularity at 0 for the imaginary
part of the function...)
(Note: the plots look even! but the coefficients still retain their
imaginary parts which means the expansion will be interms of both
cosines and sines!)
(Note: what I know is that if the function is even the coefficients
will be real and when converted back to the trig form of the fourier
expansion the sine terms disappear!)
Ns = 20;
last = T;
twiddle = 0.00001;
amount = (-last:2*last/(Ns-1):last);
amountoriginal = amount + twiddle
amount = amountoriginal.^2.^0.5
sqx = sqrt (0.5*pi./amount);
takesample = sqx .* besselh (0.5, 2, amount)
plot (real(takesample))
%plot (imag(takesample))
coefficients = fft (takesample)
%result = coefficients;
%n = (0:1:Ns-1);
n = (0-ceil(Ns/2):1:ceil(Ns/2)-1)
%A = sum (coefficients.*exp(i*n*3))
result = 1/Ns * sum(coefficients.*(exp((2*pi*i/(2*T))*n*3)));
C)
Maple: failed attempt at evaluating coefficients explicitly!
evalf(Int(Int(exp(-I*150*sqrt(X^2+Z^2))*I/(150*sqrt(X^2+Z^2)),
X=0.0001..0.03), Z=0.0001..0.02)); |
