Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

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:

enter image description here

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]]]

enter image description here

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

share|improve this question

2 Answers 2

up vote 4 down vote accepted

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]

enter image description here

share|improve this answer
    
This code works, Thanks a lot! –  FromBaidu Sep 6 '13 at 1:58
    
but I don't understander qs @ Fourier @ qs @ data –  FromBaidu Sep 6 '13 at 1:59

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};

Mathematica graphics

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]

Mathematica graphics

{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]]

Mathematica graphics

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).

Mathematica graphics

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

Mathematica graphics

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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