Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: Re: Numerical Differentiation using Fourier Transform

• To: mathgroup at smc.vnet.net
• Subject: [mg32982] Re: [mg32948] Re: Numerical Differentiation using Fourier Transform
• From: "Philippe Dumas" <dumasphi at noos.fr>
• Date: Sat, 23 Feb 2002 02:37:58 -0500 (EST)
• References: <a4vgmt$r4b$1@smc.vnet.net> <200202210707.CAA02160@smc.vnet.net>
• Reply-to: "Philippe Dumas" <dumasphi at noos.fr>
• Sender: owner-wri-mathgroup at wolfram.com

Hi Jens-Peer

I don't think you're right about the previous question by Mike Croucher.
His practical problem is to use an approximation for a function given a
finite set of sampling points. Tha'ts quite common thing !
And he clearly wants to use a Fourier approximation (instead of
interpolation polynomials, or whatever). This has decisive advantages if you
want to remove noise for example from an experimental array of values.
One point, however, is very important to have in mind. This is the so called
'Gibbs' phenomenon arising because of discontinuity between the first point
and the last point of the array. If the function does not take the same
value, then you will get strange things due to non absolute convergence of
the Fourier approximation. Just use the following code by replacing the
function f[x] by a function taking different values at f[0] and f[xmax]. You
will see what I mean by very strange (and yet perfectly well understood)
phenomenon §

**************************
Clear[test, xt, Foutest, n, nHalf, xmax, ismooth]
xmax = 6.
nHalf = 200
n = 2*nHalf + 1
xt = Table[i*xmax/n , {i, 0, n - 1}];
f[x_] := E^(-(x - xmax/2.)^2)
test = Table[f[i*xmax/n] , {i, 0, n - 1}];
Foutest = Chop[Fourier[test]];
FouDer = Foutest;
FouDer[[1]] = 0.;
ismooth = 1000000;
freq1 = -I*Table[i*Exp[-(i/ismooth)^2], {i, 1, nHalf}]; (* Here is one half
of  the frequency terms *)
freq2 = I*Table[i*Exp[-(i/ismooth)^2], {i, 1, nHalf}]; (* Here is the other
half . Notice the Reverse order below*)

freq = Join[{0}, freq1, Reverse[freq2]]; (* here is the full array of
I*omega terms to get the derivative  *)

FouDer = FouDer*freq;
der = InverseFourier[FouDer];
ListPlot[Transpose[{xt, test}], PlotJoined -> True]
dplot = ListPlot[Transpose[{xt, der}], PlotJoined -> True,
    PlotStyle -> {Thickness[0.001]}]
dplottheo =
  Plot[f'[x], {x, 0, xmax},
    PlotStyle -> {Thickness[0.007], RGBColor[1, 0, 0]}]
Show[dplot, dplottheo]
*********************************


Philippe Dumas
99, route du polygone
03 88 84 67 80
67100 Strasbourg

----- Original Message -----
From: "Jens-Peer Kuska" <kuska at informatik.uni-leipzig.de>
To: mathgroup at smc.vnet.net
Subject: [mg32982] [mg32948] Re: Numerical Differentiation using Fourier Transform


> Hi,
>
> you really know that you dont mix up continuous and discrete
> Fourier transforms ? For a discrete to discrete Fourier transform
> you have sampled the values of  the 2Pi periodic function f with
> N points x[k]=2Pi*k, k=0,1,..,N-1 and f[k]=f[x[k]] and with the
> DFT you calculate a[j] with
>
> f[k]=Sum[a[j]*Exp[2Pi*I*j*k/N],{j,0,N-1}]/N
>
> Can yo tell me what you mean with the continuous derivative ?
> of the discrete f[k] ? This ist simply not defined because k can
> have only discrete values.
>
> You can use a discrete approximation to a derivative (f[i+1]-f[i-1])/2
> and interprete this as a convolution
>
> test = Table[Sin[2Pi*x/20], {x, 0, 20, 20/511}];;
>
> fkernel = Fourier[
>            RotateLeft[
>              Drop[Table[
>                If[Abs[i] <= 1, -i, 0]/2, {i, -Length[test]/2,
> Length[test]/2}],
>                -1],
>              Length[test]/2]
>           ];
>
> dfftest = InverseFourier[Chop[fkernel*Fourier[test]]];
>
>
> But is is up to you, to use a higher order approximation of the
> first derivative (f[i-2] - 8*f[i-1] + 8*f[i+1] - f[i+2])/12
>
> Regards
>   Jens
>
>
> Mike wrote:
> >
> > Hi
> >
> > I need to differentiate a function that i only know numerically and
> > for various reasons I have been asked to do this via a Fourier
> > Transform Method.  Before doing it on the horrible function I am
> > dealing with, I thought I would get to grips with the method using
> > much easier functions.
> >
> > So say I want to differentiate
> >
> > E^(-x^2)
> >
> > w.r.t x using Fourier Transforms : analytically i could proceed as
> > follows
> >
> > FourierTransform[E^(-x^2), x, w]] = 1/(Sqrt[2]*E^(w^2/4))
> >
> > multiply this by I*w and taking the inverse Transform yields
> >
> > (-2*x)/E^x^2
> >
> > as expected so i am along the right lines here.  Now I assume that I
> > only know this function numerically and see if i can reproduce this
> > result.  First I made a table of values :
> >
> > test = Table[E^(-x^2), {x, 0.01, 20, 0.01}];
> >
> > and find it's DFT using
> >
> > Fourier[test]
> >
> > now in order to find the derivative I need to multiply this by I*w and
> > take the inverse transform.  My problem is how to do this?  what are
> > my values of w?  The DFT looks very different to the one I found
> > analytically and to be honest I don't have a clue of whats going on.
> > Any help would be appreciated - Thanks
> >
> > Mike
>
>

• References:
  • Re: Numerical Differentiation using Fourier Transform
    • From: Jens-Peer Kuska <kuska@informatik.uni-leipzig.de>