Fourier DFT scaling problem

[Arthur]

Hi all.

I hope someone can help me with what seems to be a fairly elementary problem, nevertheless I have spent a good deal of time trying to come to a solution without any success. I am fairly new to the Fourier transform.

My problem is that the output of the Fourier[] command has a scaling dependency on my sampling parameters which I cannot figure out how to remove.

In general I have a numerical function Psi[x] (a wavefunction) and I would like to compute an approximation to the Fourier transform F[Psi[x]] ~~ Psihat[k].

I am assuming the best way to achieve this is using the built in Fourier[] function.

Now although I expect to have a numerical function Psi[x] my problem is apparent even when I define Psi[x] analytically.

If I have



Psi[x_] := (Pi^(-1/4)) (E^(-x^2/2))

and take the continuous FT

FourierTransform[Psi[x], x, k]

OUT: (Pi^(-1/4)) (E^(-k^2/2))

The output is as expected an identical function of k (since Psi[x] was a normalized eigenfunction of the FT).

If I try to take the DFT of Psi[x] however:

/******CODE******/

n1 = 2^8.; (* sampling points *)
L = 32.; (* domain [-L/2,L/2] *)
dx = L/n1; (*step size*)

X = Table[-dx n1/2 + (n - 1) dx, {n, n1}]; (* lattice *)

rot[X_] := RotateLeft[X, n1/2]; (* rotate input/output of fourier *)
F[X_] := rot@Fourier@rot@X; (* DFT *)

SetOptions[ListLinePlot, PlotRange -> {{-L/2, L/2}, {0, 1}},
DataRange -> {-L/2, L/2}, Filling -> Axis];

ListLinePlot[Abs[Psi[X]]^2]
ListLinePlot[Abs[F@Psi[X]]^2]


/******END******/
The second graph of the squared modulus of F@Psi[X] is evidently incorrectly scaled. Moreover this scaling error is not improved by increasing the sampling frequency, and is in fact exacerbated by it.

It seems reasonable that the scaling be dependent on the value dx, though I have done a lot of experimentation with the FourierParameters command to remove this dependence though I can't get it to work in a way that is independent of the choice of Psi[x].

Any help would be *greatly* appreciated!

Thanks and best regards,

Arthur

[Kevin J. McCann]

Fourier will do the job, but I do not understand why you are doing the
rotates.

Here is a slightly different take on your problem with a Parseval check
at the end. (Just copy and drop it into an empty notebook.)

Cheers,

Kevin


Psi[x_] := \[Pi]^(-1/4) E^(-(x^2/2))

\[ScriptCapitalN] = 2^8;
L = 32.0;
\[CapitalDelta]x = L/\[ScriptCapitalN];
x = -16 + Table[i, {i, 0, \[ScriptCapitalN] - 1}]*\[CapitalDelta]x;
\[CapitalDelta]k = 1/L;

f = Psi[x];
k = Table[n \[CapitalDelta]k, {n, 0, \[ScriptCapitalN] - 1}];

F = \[CapitalDelta]x Fourier[f, FourierParameters -> {1, -1}];
ListPlot[Transpose[{k, Abs[F]^2}], PlotRange -> All, Joined -> True]
(* This part rotates the k-spectrum so it "looks right" *)
krot = Table[
n \[CapitalDelta]k, {n, -(\[ScriptCapitalN]/2), \[ScriptCapitalN]/
2 - 1}];
Frot = RotateRight[F, \[ScriptCapitalN]/2];
ListPlot[Transpose[{krot, Abs[Frot]^2}], PlotRange -> All,
Joined -> True]
(* Check that Parseval is right *)
{Total[Abs[f]^2] \[CapitalDelta]x, Total[Abs[F]^2] \[CapitalDelta]k}

[W Craig Carter]

Hello Kevin,
The snippet you provided was useful and timely for me. Would you happen
to have something similar for spectral derivatives using Fourier[],
FourierDCT[], and FourierDST[]?
Thanks in any case,

Craig Carter

[Joseph Gwinn]

In article <jpflom$ktd$1...@smc.vnet.net>, Arthur <art...@ymail.com>
wrote:

In the documentation on Fourier[], the actual equation is given. The
coefficients depend on the length of the list given to Fourier. You
will have to divide by Sqrt[length] to have coefficients independent of
length, if I recall. There are also a number of alternative
normalizations one can choose using options. But you will need to
figure out the correct normalization and options for your application.

Joe Gwinn

[Dana DeLouis]

> FourierTransform[Psi[x], x, k]
>
> OUT: (Pi^(-1/4)) (E^(-k^2/2))

Hi. I don't have a solution, but does anything here help?
Does changing the default settings help in some way?

xData={-1,+1};
xSignal={+1,-1};
xDefault = {0,1};

Psi[x_] := (Pi^(-1/4)) (E^(-x^2/2))

SetOptions[{FourierTransform,InverseFourierTransform},FourierParameters->xDefault];

FourierTransform[Psi[x],x,k]
E^(-(k^2/2))/\[Pi]^(1/4)

SetOptions[{FourierTransform,InverseFourierTransform},FourierParameters->xSignal];

FourierTransform[Psi[x],x,k]
Sqrt[2] E^(-(k^2/2)) \[Pi]^(1/4)

SetOptions[{FourierTransform,InverseFourierTransform},FourierParameters->xData];

FourierTransform[Psi[x],x,k]
E^(-(k^2/2))/(Sqrt[2] \[Pi]^(3/4))

= = = = =
Dana DeLouis
Mac & Math 8
= = = = =