How to draw a graph about the FourierTransform of a 2D Function?

Question

Here is a 2D function, and it cannot be calculated by the FourierTransfrom

n = 2; δ = 0.2;
p[x_, y_] = Cosh[n*Sqrt[2 Pi]*δ*Sqrt[x^2 + y^2]]*Exp[-((δ^2*(x^2 + y^2) + 2 n^2*Pi)/2)];

p[x,y] ContourPlot

and my problem is how to draw a graph about its FourierTransform Function?

Here is an example: if I have a function that can be calculated by FourierTransform like below:

n = 2; δ = 0.2;
p[x_, y_] = Cosh[n*Sqrt[2 Pi]*δ*x]*Cosh[n*Sqrt[2 Pi]*δ*y]*Exp[-((δ^2*x^2+ 2 n^2*Pi)/2)]*Exp[-((δ^2*y^2+ 2 n^2*Pi)/2)];
uu1[tx_, ty_] = FourierTransform[p[x, y], {x, y}, {tx, ty}];
ContourPlot[p[x, y], {x, -40, 40}, {y, -40, 40}]
ContourPlot[uu[tx, ty]^2, {tx, -0.5, 0.5}, {ty, -0.5, 0.5}]

we will get the answer below:

And here, I suppose the function p[x,y] can not use FourierTransform ,I can use the function Fourier

data = Table[p[x, y], {x, -40, 40, 0.4}, {y, -40, 40, 0.4}];
ListDensityPlot[Re[Fourier[data]]]

but the answer is not I wanted! I don't know why? and I hope I can get help from you!

Answer 1

You are missing two things:

• Quadrant swapping, to put the low frequency components of the FFT at the centre of the array (look up Matlab's fftshift function for more information)

• Zero padding, which effectively interpolates between data points in the transformed domain and helps to reduce the "crosshair" effect you can see in Nasser's results. See e.g. this for more details

E.g.

data = Table[p[x, y], {x, -40, 40, 1}, {y, -40, 40, 1}];

(* pad array with zeros - here I am padding to 1024x1024 *)
data = ArrayPad[data, {#, 1024 - Length[data] - #} &[512 - Round[Length[data]/2]]];

(* FFT with quadrant swapping *)
result = Re @ With[{qs = RotateRight[#, {512, 512}] &}, qs @ Fourier @ qs @ data];

(* Cut out the interesting part of the data and plot it *)
ListContourPlot[ArrayPad[result, -400]^2]

Answer 2

ps. I tried to get the symbolic Fourier, but gave up. Too hard to do this symbolically.

expr = p[x, y] /. {n -> 2, del -> 0.5};

FourierTransform[expr, {x, y}, {w1, w2}]

Not working.... So, try numerical FFT. Using meshgrid and pick some spatial x,y frequencies, which I am not sure now if it meets spatial Nyquist or not. Then using FFT (Fourier) on the 2D image.

n = 2; del = 0.5;

p[x_, y_] = Cosh[n*Sqrt[2 Pi]*del*Sqrt[x^2 + y^2]]* Exp[-((del^2*(x^2 + y^2) + 2 n^2*Pi)/2)];

meshgrid[x_List, y_List]:={ConstantArray[x, Length[x]],Transpose@ConstantArray[y,Length[y]]}

{xx, yy} = meshgrid[Range[-4 Pi, 4 Pi, .05], Range[-4 Pi, 4 Pi, .05]];
c = p[xx, yy];

ListPlot3D[Transpose@c, DataRange -> {{-4 Pi, 4 Pi}, {-4 Pi, 4 Pi}}, PlotRange -> All]

{nRow, nCol} = Dimensions[c];
d = c*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
fw = Fourier[d, FourierParameters -> {1, 1}];
(*adjust for better viewing as needed*)
fudgeFactor = 10;
abs = fudgeFactor*Log[1 + Abs@fw];
Labeled[Image[abs/Max[abs], ImageSize -> 300],Style["Magnitude spectrum", 18]]

If you use different sampling in the above meshgrid, resolution of the 2D looking spectrum will change, like this (these are all centered so that center of spectrum is DC frequency).

arg = Arg@fw;
Labeled[Image[arg/Max[arg], ImageSize -> 300],Style["Phase spectrum", 18]]