1
$\begingroup$

I am trying to calculate this simple FT on a grid {x,y} but it tasks too long and gives me some error message even though the results are Ok. Are there other efficient ways to do this on MMA?

this is the Fourier transform

FT[xp_, yp_] := 
 1/\[Pi] NIntegrate[
   1/Sqrt[2 \[Pi] Sqrt[ (x^2 + y^2)]]
     Exp[-I xp x] Exp[-I yp y], {x, -1, 1}, {y, -Sqrt[1 - x^2], Sqrt[
    1 - x^2]}]

I would like to calculate this on a square grid of step=0.1 but it takes too long, here I used step=1 and it takes 80 sec is this reasonable?

step = 1.;
res = ParallelTable[{x1, y1, 
     Abs[FT[x1, y1] Conjugate[FT[x1, y1]]]}, {x1, -10., 10., 
     step}, {y1, -10., 10., step}]; // AbsoluteTiming   

NIntegrate::slwcon: 
   Numerical integration converging too slowly; suspect one of the following:
    singularity, value of the integration is 0, highly oscillatory integrand,
    or WorkingPrecision too small.    

{81.4518, Null}   

ListPlot3D[Flatten[res, 1], PlotRange -> All, BoxRatios -> {1, 1, 1}, 
 ColorFunction -> "Rainbow"]   

enter image description here

Update#1

the answer by @Daniel Lichtblau gives the output very fast but the results are not exactly the same as in the question. For example, if we set the plot range as follows we can see there is a wave-like profile at the base of the peak which is absent in the answer?

ListPlot3D[Flatten[res, 1], PlotRange -> {0,0.03}, BoxRatios -> {1, 1, 1}, 
 ColorFunction -> "Rainbow"]   

enter image description here

Here are the results from the answer

enter image description here

Update#2

As @yarchik mentioned the 2D integral in the question is not FT as I claimed which is a stupid mistake, the Best way to speed it up is by transforming it to polar coordinate, and results come in no time and looks like this

enter image description here

$\endgroup$
8
  • 1
    $\begingroup$ Have you tried to sample your function on a regular grid and to apply Fourier (which applies FFT)? $\endgroup$ Jan 23 at 13:09
  • $\begingroup$ @HenrikSchumacher, what do you mean by regular grid? isn't the same what I am doing? $\endgroup$ Jan 23 at 13:58
  • 1
    $\begingroup$ Please, have a look at the documentation of Fourier first and see whether it suits your needs. $\endgroup$ Jan 23 at 14:21
  • 1
    $\begingroup$ What you are computing is not the 2D Fourier transform. You introduce artifacts, such as the wave-like profile at the base of the peak, by integrating in a small disk-shape domain. Also, it is known analytically that the Fourier transform of a gaussian peak is again a gaussian peak. The best way to verify is to perform a transform twice. You should get the original peak in the case of FT. $\endgroup$
    – yarchik
    Jan 24 at 9:42
  • 2
    $\begingroup$ No need to remove. Someone may find the whole discussion useful. $\endgroup$
    – yarchik
    Jan 24 at 11:09

1 Answer 1

4
$\begingroup$

This is more or less the suggestion from @HenrikSchumacher. At least that's the intent. Note that I make the step slightly offset from unity in order to avoid a singularity at the origin.

step = 1.01;
pts = Table[{x, y, 1/Sqrt[2 \[Pi] Sqrt[(x^2 + y^2)]]}, {x, -10., 10., 
    step}, {y, -10., 10., step}];
ft = Fourier[pts[[All, All, 3]]];
dims = Dimensions[ft]

(* Out[38]= {20, 20} *)

Now center the transform by rotating halfway. This pusts the DC component in the middle, modulo any off-by-one error in the code.

shift = Ceiling[dims[[1]]/2]; shiftedft = 
 RotateLeft[Map[RotateLeft[#, shift] &, ft], shift];

Recreate the proper coordinates.

newvals = 
  Table[{10*(i - shift)/shift, 10*(j - shift)/shift, 
    Abs[shiftedft[[i, j]]*Conjugate[shiftedft[[i, j]]]]}, {i, 
    dims[[1]]}, {j, dims[[2]]}];

Now we can (re)create the plot.

ListPlot3D[Flatten[newvals, 1], BoxRatios -> {1, 1, 1}, 
 PlotRange -> All]

enter image description here

$\endgroup$
3
  • $\begingroup$ that is very fast but did not give the desired results. Kindly see my update $\endgroup$ Jan 23 at 19:16
  • $\begingroup$ @valarmorghulis, you should play with discretization parameters (see first Table[..] in Daniel's answer) to make the substructure more visible. $\endgroup$
    – Rom38
    Jan 24 at 10:34
  • $\begingroup$ @Rom38, I did a mistake by calling that FT which is just an ordinary 2D integral, see my update 2 $\endgroup$ Jan 24 at 10:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.