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>