Differentiation and Integration

Needs["StableM`"];
SetOptions[FourierTransform,FourierParameters-> {1,1}];
SetOptions[InverseFourierTransform,FourierParameters-> {1,1}];
SetOptions[Fourier, FourierParameters -> {-1,1}];
SetOptions[InverseFourier, FourierParameters -> {-1,1}];
ReOrder::usage = "ReOrder[list] reorders a frequency domain list for the time domain and vice versa.  If Length[list] is odd the output list will be truncated to the next smallest integer.";
XTable::usage="XTable[SR, N] puts out a list of x values that correspond to the stable x-axis in the time domain for a given length, N, and SR of the DFFT list.";
DFTCF::usage= "DFTCF[cf,SR,N] puts out a discrete FT of a distribution specified by the characteric function cf in the frequency domain compatible with FourierParameters -> {-1, 1}.  SR is the sampling rate; SR is greater than 0.  N is the number of points and must be even.  The output is scaled ready for inversion.";
CF1c::usage= "CF1c[u, α, β, γ, δ] outputs a compiled version of the S(α, β, γ, δ; 1) characteristic function.  All the inputs must be real numbers.";
ReOrder=With[{n=Length[#]},Join[Drop[#,Floor[n/2]],Drop[#,-Floor[n/2]]]]&;
XTable  = Compile[{SR, {N,_Integer}},Module[{n},
n=2 Floor[N/2];
Range[-n/2,n/2-1,1]/(SR)]];
DFTCF [cf_,SR_,n_Integer]:=Module[{rlist,n2},
n2=Floor[n/2];
rlist=2. Pi SR/n Join[Range[0,n2-1],Range[-n2,-1]];
Map[cf[ #]&,rlist]SR/n
	];
CF1c=Compile[{#1,#2,#3,#4,#5},If[#1==0,1,If[#2==1,E^(I*#5*#1-#4*Abs[#1]*((2*I*#3*Log[Abs[#1]]*Sign[#1])/Pi+1)),E^(I*#1*#5-Abs[#1]^#2*#4^#2*(1-I*Sign[#1]*#3*Tan[(Pi*#2)/2]))]]];

Differentiation in the transform domain can be done with multiplication of each term by -i t.

differentiate = Compile[{t} , -I t];
scale = 8192/16;
sdft = DFTCF[CF1c[#1, 1.5, 0.5, 1, 0] &, 16, 8192];
df = scale DFTCF[differentiate[#1] &, 16, 8192];(*this isn't scaled so it is multiplied by scale*)
conv = Chop[InverseFourier[sdft df]];
g1 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv]}], PlotRange -> {{-20, 20}, All}, Joined -> True, PlotStyle -> Blue]

differentiate=Compile[{t} ,-I t];
scale=8192/16;
sdft=DFTCF[CF1c[#1,1.5,0.5,1,0]&,16,8192];
df=scale DFTCF[differentiate[#1]&,16,8192];(*this isn't scaled so it is multiplied by scale*)
conv=Chop[InverseFourier[sdft df]];
g1=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]}],PlotRange->{{-20,20},All},Joined->True,PlotStyle->Blue]

Differentiate twice

conv = Chop[InverseFourier[ sdft df df]];
g1 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv]}], PlotRange -> {{-20, 20}, All}, Joined -> True, PlotStyle -> Blue]

conv=Chop[InverseFourier[ sdft df df]];
g1=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]}],PlotRange->{{-20,20},All},Joined->True,PlotStyle->Blue]

The following creates a distribution function.  This is also multipication operation for each term but the first element must be replaced with π.  The constant of integration has to be added; in this case it is 1/2.

integrate = Compile[{t} , If[t == 0, Pi, -1/(I t)]];
scale = 8192/16;
int = scale DFTCF[integrate[#1] &, 16, 8192];
conv = Chop[InverseFourier[int sdft]];
g1 = ListPlot[Transpose[{XTable[16, 8192], ReOrder[conv] + 1/2}], PlotRange -> {{-20, 20}, All}, Joined -> True, PlotStyle -> Blue]

integrate=Compile[{t} ,If[t==0,Pi,-1/(I t)]];
scale=8192/16;
int=scale DFTCF[integrate[#1]&,16,8192];
conv=Chop[InverseFourier[int sdft]];
g1=ListPlot[Transpose[{XTable[16,8192],ReOrder[conv]+1/2}],PlotRange->{{-20,20},All},Joined->True,PlotStyle->Blue]

Clearly this is not correct.  The FFT is expecting a periodic function that begins and ends at the same point.  So here is a trick, make it look like a periodic function: join the pdf with an inverted reverse of itself.  Put it back in the Fourier domain, use double the integral routine, but for the display only use the left half of the result.

conv = Chop[InverseFourier[ sdft]];
sdft2 = Fourier[Join[ReOrder[conv], -Reverse[ReOrder[conv]]]];
int = 2 scale DFTCF[integrate[#1] &, 2 16, 2 8192];
conv = Chop[InverseFourier[int sdft2]];
g1 = ListPlot[Transpose[{XTable[16, 8192], Drop[conv, -8192] + 1/2}], PlotRange -> {{-20, 20}, All}, Joined -> True, PlotStyle -> Blue]

conv=Chop[InverseFourier[ sdft]];
sdft2=Fourier[Join[ReOrder[conv],-Reverse[ReOrder[conv]]]];
int=2 scale DFTCF[integrate[#1]&,2 16,2 8192];
conv=Chop[InverseFourier[int sdft2]];
g1=ListPlot[Transpose[{XTable[16,8192],Drop[conv,-8192]+1/2}],PlotRange->{{-20,20},All},Joined->True,PlotStyle->Blue]

This is what both sides of it looks like.

g1 = ListPlot[Transpose[{XTable[2 4, 2 8192], conv + 1/2}], PlotRange -> {{-2 scale, 2 scale}, All}, Joined -> True, PlotStyle -> Blue]

g1=ListPlot[Transpose[{XTable[2 4,2 8192],conv+1/2}],PlotRange->{{-2 scale,2 scale},All},Joined->True,PlotStyle->Blue]

txcolor=RGBColor[181/255,18/255,27/255];
t1 = Graphics[Text[Hyperlink[Framed[Style["Financial Home", {txcolor, Bold, Larger}]], "../map/Overview.html"], {-5, 0}]];
t2 = Graphics[Text[Hyperlink[Framed[Style["Next", {txcolor, Bold, Larger}]], "../ft05/Convolution.html"], {0, -0.5}]];
t3 = Graphics[Text[Hyperlink[Framed[Style["Download Notebook", {txcolor, Bold, Larger}]], "../ft04/Differentiation.nb"], {5, 0}]];
t4 = Graphics[Text[Hyperlink[Framed[Style["Previous", {txcolor, Bold, Larger}]], "../ft03/ShiftScale.html"], {0, 0.5}]];
Show[ t1, t2, t3, t4, ImageSize -> 500]
CellPrint[TextCell["\n\n© Copyright " <> ToString[DateList[][[1]]] <> " mathestate    " <> DateString[DateList[], "DateShort"], "Text"]]

© Copyright 2008 mathestate    Fri 29 Feb 2008