Fast Creation of Approximate 
Local  Fourier and InverseFourier 
Continuous Transforms

The Package NFourierLocalTransform

(*:Sources:
	A . Nersessian .  
                  Quasiexponential function and Fourier quasitransforms (in Russian) .   Reports of NAS RA, vol .  101, N1, 2001, pp 5 - 11. 
	 *)
	  

	  

BeginPackage["Calculus`FourierTransform`"];
Off[NIntegrate::"slwcon"];
Off[
        NIntegrate::"ncvb"];
Off[NIntegrate::"ploss"]; 
Off[
        General::spell];
Off[
        General::spell1];
Off[
        General::"ovfl"];
Off[
        General::"unfl"];
Off[
        FindRoot::"frsec"];
Off[
        NIntegrate::"inovf"];
Off[Break::"nofwd"];

 NFourierLocalTransform::usage = " NFourierLocalTransform[expr, \
x->{t,tmin,tmax}] gives the symbolicform of the Fourier  transform of expr, a \
function of x, as a function of t on the interval from tmin to tmax.  It is \
defined by 1/Sqrt[2 Pi] NIntegrate[expr Exp[I w t], {t, -Infinity, \
Infinity}].    NFourierLocalTransform[expr, x->{t,tmin,tmax}, \
FourierParameters -> {a, b}] gives the symbolic form of the parameterized \
transform Sqrt[Abs[b]/(2 Pi)^(1-a)] NIntegrate[expr Exp[I b w t], {t, \
-Infinity, Infinity}].";

\(NInverseFourierLocalTransform::usage = "NInverseFourierTransform[expr, w, t] gives the numerical value of the inverse
Fourier transform of expr, a function of w, at t.  It is defined by
1/Sqrt[2 Pi] NIntegrate[expr Exp[-I w t], {w, -Infinity, Infinity}].
NInverseFourierTransform[expr, w, t, FourierParameters -> {a, b}] gives
the numerical value of the parameterized inverse transform
Sqrt[Abs[b]/(2 Pi)^(1+a)] NIntegrate[expr Exp[-I b w t],
{w, -Infinity, Infinity}].";\)


 \(NFourierLocalTransform::small = \*""(Subscript[L, 2]) norm of the function is smaller than the required precision."";\)

 NFourierLocalTransform::longcoef = "You are just calculating \
`1` integrals. It would take a long time for slow computers. Otherwise \
decrease precision or the interval on which you are calculating the  Fourier \
transform.";

 \(NFourierLocalTransform::unknnopt = "Unknown Option(s) `1`.";\)

 NFourierLocalTransform::optx = "Unknown option `1`.";

  NFourierLocalTransform::lngtm = "It would take a long time to \
calculate about-> `1`  integrals.";


 \(NInverseFourierLocalTransform::small = \*""(Subscript[L, 2]) norm of the function is smaller than the required precision."";\)

 NInverseFourierLocalTransform::longcoef = "You are just \
calculating `1` integrals. It would take a long time for slow computers. \
Otherwise decrease precision or the interval on which you are calculating the \
 Fourier transform.";

 NNInverseFourierLocalTransform::unknnopt = "Unknown Option(s) \
`1`.";

 \(NInverseFourierLocalTransform::optx = "Unknown option `1`.";\)

  NInverseFourierLocalTransform::lngtm = "It would take a long \
time to calculate about-> `1`  integrals.";




Options[
            NFourierLocalTransform] = Options[
              NInverseFourierLocalTransform] = 
            Flatten[{InterpolationPoints -> Automatic, 
                InterpolationOrder -> Automatic, 
                FourierParameters -> {0, 1}, Symbolic -> False, 
                Options[NIntegrate]}];

Begin["`Private`"];
<< 
        Utilities`FilterOptions`;
\(flst[l_List] := 
          Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

            RotateRight[Fourier[fy]/\@nn, Floor[nn/2]]];\)


IntrvlAB[f_, ϵ_, δ_, x_, opt___] := 
          Module[{NInt, Sk, LR, RL, FindPnt, St, En, StLimit, EndLimit, 
              StartPoint, 
              optI}, 
{StLimit, 
                EndLimit} = {-∞, ∞}; 
optI = 
              FilterOptions[NIntegrate, opt]; 
StartPoint = 
              0; 
LR = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StartPoint, EndLimit}, 
                      Evaluate[optI]] ≤ ϵ, -1, 
                1]; 
RL = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StLimit, StartPoint}, 
                      Evaluate[optI]] ≤ ϵ, 
                1, -1]; 
Sk = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StLimit, EndLimit}, 
                      Evaluate[optI]] ≤ ϵ, 0, 
                1]; 
FindPnt[a_, b_, k_, j_] :=  
              Module[{NInt, St, En, Pnt}, 
Pnt = 
                  If[j == 1, ((a + b)/2), a]; 
St = 
                  If[k == 1, Pnt, StLimit]; 
En = 
                  If[k == 1, EndLimit, 
                    Pnt]; 
NInt = √NIntegrate[
                      Evaluate[Abs[f]^2], {x, St // N, En // N}, 
                      Evaluate[
                        optI]]; 
Which[
Abs[a - b] < δ && 
                    j == 1, 
                  If[((LR + RL + k LR - k RL)/2) == k, b, 
                    a], 
k (( LR + RL + k LR - k RL)/2) NInt < k (( LR + RL + k LR - k RL)/2) ϵ, 
                  FindPnt[(3/4) a + 
                      j  (a/4) - (1/2) Abs[Abs[Sign[a - 1]] - 1], Pnt, k, 
                    1], 
k (( LR + RL + k LR - k RL)/2) NInt > k (( LR + RL + k LR - k RL)/2) ϵ, 
                  FindPnt[2  a (( 1 - j)/2) + ((a + b)/2) ((1 + j)/2), 
                    b, k, j], 
NInt == ϵ && 
                    j == -1, 
                  a
]
]; 
Which[
Sk == 1  &&  
                StLimit == -∞  &&  
                EndLimit == ∞, {FindPnt[RL, 
                    0, -1, -1] // N, 
                FindPnt[LR, 0, 1, -1] // 
                  N}, 
Sk == 1  &&  
                StLimit ≠ -∞  &&  
                EndLimit ≠ ∞, 
{LR = -1; RL = 1; FindPnt[StLimit, EndLimit, -1, 1] // N, LR = 1; 
                RL = -1; 
                FindPnt[StLimit, EndLimit, 1, 1] // 
                  N}, 
Sk == 0, 
              Message[ 
                NFourier (NInverseFourier) 
                  LocalTransform::small]
]
];


IntrvlABI[f_, ϵ_, δ_, x_, opt___] := 
          Module[{NInt, Sk, LR, RL, FindPnt, St, En, StLimit, EndLimit, 
              StartPoint, 
              optI}, 
{StLimit, 
                EndLimit} = {-∞, ∞}; 
optI = 
              FilterOptions[NIntegrate, opt]; 
StartPoint = 
              0; 
LR = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StartPoint, EndLimit}, 
                      Evaluate[optI]] ≤ ϵ, -1, 
                1]; 
RL = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StLimit, StartPoint}, 
                      Evaluate[optI]] ≤ ϵ, 
                1, -1]; 
Sk = 
              If[√NIntegrate[
                      Evaluate[Abs[f]^2], {x, StLimit, EndLimit}, 
                      Evaluate[optI]] ≤ ϵ, 0, 
                1]; 
FindPnt[a_, b_, k_, j_] :=  
              Module[{NInt, St, En, Pnt}, 
Pnt = 
                  If[j == 1, ((a + b)/2), a]; 
St = 
                  If[k == 1, Pnt, StLimit]; 
En = 
                  If[k == 1, EndLimit, 
                    Pnt]; 
NInt = √NIntegrate[
                      Evaluate[Abs[f]^2], {x, St // N, En // N}, 
                      Evaluate[
                        optI]]; 
Which[
Abs[a - b] < δ && 
                    j == 1, 
                  If[((LR + RL + k LR - k RL)/2) == k, b, 
                    a], 
k (( LR + RL + k LR - k RL)/2) NInt < k (( LR + RL + k LR - k RL)/2) ϵ, 
                  FindPnt[(3/4) a + 
                      j  (a/4) - (1/2) Abs[Abs[Sign[a - 1]] - 1], Pnt, k, 
                    1], 
k (( LR + RL + k LR - k RL)/2) NInt > k (( LR + RL + k LR - k RL)/2) ϵ, 
                  FindPnt[2  a (( 1 - j)/2) + ((a + b)/2) ((1 + j)/2), 
                    b, k, j], 
NInt == ϵ && 
                    j == -1, 
                  a
]
]; 
Which[
Sk == 1  &&  
                StLimit == -∞  &&  
                EndLimit == ∞, {FindPnt[RL, 
                    0, -1, -1] // N, 
                FindPnt[LR, 0, 1, -1] // 
                  N}, 
Sk == 1  &&  
                StLimit ≠ -∞  &&  
                EndLimit ≠ ∞, 
{LR = -1; RL = 1; FindPnt[StLimit, EndLimit, -1, 1] // N, LR = 1; 
                RL = -1; 
                FindPnt[StLimit, EndLimit, 1, 1] // 
                  N}, 
Sk == 0, 
              Message[ 
                NInverseFourierLocalTransform::small]
]
];


segment[f_, ϵ_, x_, opt___] := 
          Block[{δ = 100, ab}, 
ab = 
              IntrvlAB[f, (10^(-ϵ)), δ, x, 
                opt]; 
If[((ab // Last) - (ab // 
                        First)) ≥ 200 δ, ab, 
              IntrvlAB[f, (10^(-ϵ)), 
                Max[(((ab // Last) - (ab // First))/300.) , 1], x, 
                opt]]
];

segmentI[
          f_, ϵ_, x_, opt___] := 
        Block[{δ = 100, ab}, 
ab = 
            IntrvlABI[f, (10^(-ϵ)), δ, x, 
              opt]; 
If[((ab // Last) - (ab // 
                      First)) ≥ 200 δ, ab, 
            IntrvlAB[f, (10^(-ϵ)), 
              Max[(((ab // Last) - (ab // First))/300.) , 1], x, 
              opt]]
];
roG[x_, 
          al_] := (((Sqrt[π]) (E^(-(((π^2) (x^2))/(4 al)))))/(Sqrt[al]));
al[m_, 
            prec_] :=   ((m (Log[2] (2^(4 - 1.5  prec)) +  Log[10] prec))/(m - 1));

\(nFT[f_, {x_ -> λ_}, {m_, nn_}, seg_List, 
            opt___] := \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Module[{d0, α0, α, h, p1, p2, cof, p, optI, ss, mm, ord, flong, symbsol, four, sint, 
              maingrid}, 

{prec, {a, b}, 
                symbsol, ord, 
                mm} = \({PrecisionGoal, FourierParameters, Symbolic, 
                    InterpolationOrder, InterpolationPoints} /. {opt}\) /. 
                Options[\ NFourierLocalTransform]; 
If[
              prec === Automatic, 
              prec = $MachinePrecision - 11]; 
If[
              ord === Automatic, ord = 3]; 
If[
              mm === Automatic, mm = 20]; 
optI = 
              FilterOptions[NIntegrate, 
                opt]; 
α0 = \ \ \(1\/\((m - 
                        1)\)\) \((m \((Log[5\/2 - Round[prec]\/10] + \ 
                        Log[10]\ prec + \(75\/100\) 
                          Log[\@\(seg[\([2]\)] - seg[\([1]\)]\)])\))\); 
α = 
              al /. FindRoot[
                  al == \ \ SetPrecision[\(1\/\((m - 
                                1)\)\) \((m\ \((Log[
                                  5\/2 - Round[prec]\/10] + \ 
                                Log[10]\ prec + \(75\/100\) 
                                  Log[\@\(seg[\([2]\)] - seg[\([1]\)]\)] + 
                                Log[m/\((\ al)\)]\/2)\))\), 
                      100], {al,  .9  α0,  .91  α0}, 
                  MaxIterations -> 3, 
                  WorkingPrecision -> 100]; 
h = 
              hh /. FindRoot[
                  SetPrecision[\(E\^\(\(α\/4\) \((1\/\(2  m\))\)\^2\)\) \((\(-\(\(\@\(π\/2\)\ \((\(-1\) + Erf[\(\@2\ hh\ π\)\/\@α])\)\)\/\(nn\ \@α\)\)\))\)\^\(1/8\), 100] == 
                    SetPrecision[10\^\(\(-prec\)/4\), 100], {hh, \(61\/100\) 
                      prec, \(62\/100\) prec}, MaxIterations -> 3, 
                  WorkingPrecision -> 
                    100]; 
h = \(115\/100\) 
                h; 
p2 = 
              Round[seg[\([2]\)]\ nn/2 + h]; 
p1 = 
              Round[seg[\([1]\)]\ nn/2 - h]; 
If[
              Abs[p2 - p1] > 30000, 
              Message[\ NFourierLocalTransform::lngtm, \ 
                Abs[p2 - p1]\ ]]; 
cof = 
              Table[NIntegrate[
                  Evaluate[roG[nn\ x - 2  s, α] f], {x, 
                    Max[seg[\([1]\)], \(2  s - 2  h\)\/nn], 
                    Min[\(2  s + 2  h\)\/nn, seg[\([2]\)]]}, 
                  AccuracyGoal -> prec + 1, Evaluate[optI]], {\ s, p1, 
                  p2}]; 


If[
              symbsol, 
\(E\^\(\(2  
                          I\ π\ \((p1 - 
                              1)\)\ λ\)\/nn\)\) \(E\^\(α \((\(\(\ \)\(λ\)\)\/nn)\)\^2\)\) \(∑\+\(s = 1\)\%\(p2 - p1 + 1\)\(E\^\(\(2  I\ π\ \ s\ λ\)\/nn\)\) 
                    cof[\([s]\)]\)
, 
maingrid = 
                Table[\((s - 1)\)/mm, {s, 1, 
                    mm\ \ nn\ }]; 
flong = 
                PadRight[cof, nn\ \ mm]; 
four = 
                Sqrt[mm\ nn]\ Fourier[flong]; 
sint = 
                Interpolation[
                  Transpose[{Append[maingrid, nn], 
                      Append[four, four[\([1]\)]]}], 
                  InterpolationOrder -> ord, 
                  PeriodicInterpolation -> 
                    True]\ ; 
\(E\^\(α\ \((λ/nn)\)\^2\)\) 
                E\^\(2\ I\ π\ \ λ\ \ \((\ p1)\)/nn\)\ \ sint[λ]]
	];\)

\(nFTI[f_, {x_ -> λ_}, {m_, nn_}, seg_List, 
            opt___] := \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Module[{d0, α0, α, h, p1, p2, cof, p, optI, ss, mm, ord, flong, symbsol, four, sint, 
              maingrid}, 

{prec, {a, b}, 
                symbsol, ord, 
                mm} = \({PrecisionGoal, FourierParameters, Symbolic, 
                    InterpolationOrder, InterpolationPoints} /. {opt}\) /. 
                Options[\ NFourierLocalTransform]; 
If[
              prec === Automatic, 
              prec = $MachinePrecision - 11]; 
If[
              ord === Automatic, ord = 3]; 
If[
              mm === Automatic, mm = 20]; 
optI = 
              FilterOptions[NIntegrate, 
                opt]; 
α0 = \ \ \(1\/\((m - 
                        1)\)\) \((m \((Log[5\/2 - Round[prec]\/10] + \ 
                        Log[10]\ prec + \(75\/100\) 
                          Log[\@\(seg[\([2]\)] - seg[\([1]\)]\)])\))\); 
α = 
              al /. FindRoot[
                  al == \ \ SetPrecision[\(1\/\((m - 
                                1)\)\) \((m\ \((Log[
                                  5\/2 - Round[prec]\/10] + \ 
                                Log[10]\ prec + \(75\/100\) 
                                  Log[\@\(seg[\([2]\)] - seg[\([1]\)]\)] + 
                                Log[m/\((\ al)\)]\/2)\))\), 
                      100], {al,  .9  α0,  .91  α0}, 
                  MaxIterations -> 3, 
                  WorkingPrecision -> 100]; 
h = 
              hh /. FindRoot[
                  SetPrecision[\(E\^\(\(α\/4\) \((1\/\(2  m\))\)\^2\)\) \((\(-\(\(\@\(π\/2\)\ \((\(-1\) + Erf[\(\@2\ hh\ π\)\/\@α])\)\)\/\(nn\ \@α\)\)\))\)\^\(1/8\), 100] == 
                    SetPrecision[10\^\(\(-prec\)/4\), 100], {hh, \(61\/100\) 
                      prec, \(62\/100\) prec}, MaxIterations -> 3, 
                  WorkingPrecision -> 
                    100]; 
h = \(115\/100\) 
                h; 
p2 = 
              Round[seg[\([2]\)]\ nn/2 + h]; 
p1 = 
              Round[seg[\([1]\)]\ nn/2 - h]; 
If[
              Abs[p2 - p1] > 30000, 
              Message[\ NInverseFourierLocalTransform::lngtm, \ 
                Abs[p2 - p1]\ ]]; 
cof = 
              Table[NIntegrate[
                  Evaluate[roG[nn\ x - 2  s, α] f], {x, 
                    Max[seg[\([1]\)], \(2  s - 2  h\)\/nn], 
                    Min[\(2  s + 2  h\)\/nn, seg[\([2]\)]]}, 
                  AccuracyGoal -> prec + 1, Evaluate[optI]], {\ s, p1, 
                  p2}]; 


If[
              symbsol, 
\(E\^\(\(2  
                          I\ π\ \((p1 - 
                              1)\)\ λ\)\/nn\)\) \(E\^\(α \((\(\(\ \)\(λ\)\)\/nn)\)\^2\)\) \(∑\+\(s = 1\)\%\(p2 - p1 + 1\)\(E\^\(\(2  I\ π\ \ s\ λ\)\/nn\)\) 
                    cof[\([s]\)]\)
, 
maingrid = 
                Table[\((s - 1)\)/mm, {s, 1, 
                    mm\ \ nn\ }]; 
flong = 
                PadRight[cof, nn\ \ mm]; 
four = 
                Sqrt[mm\ nn]\ Fourier[flong]; 
sint = 
                Interpolation[
                  Transpose[{Append[maingrid, nn], 
                      Append[four, four[\([1]\)]]}], 
                  InterpolationOrder -> ord, 
                  PeriodicInterpolation -> 
                    True]\ ; 
\(E\^\(α\ \((λ/nn)\)\^2\)\) 
                E\^\(2\ I\ π\ \ λ\ \ \((\ p1)\)/nn\)\ \ sint[λ]]
	];\)

 \(NFourierLocalTransform[ff_, x_ -> {λ_, aa_, bb_}, 
            opt___] := 
Module[{norm, normnew, f, s, a, b, 
              nn, m, prec, ss, opt1, opt2, 
              optlist}, 

optI = 
              FilterOptions[NIntegrate, 
                opt]; 
opt1 = {InterpolationPoints, 
                InterpolationOrder, FourierParameters, Symbolic, 
                AccuracyGoal, Compiled, GaussPoints, MaxPoints, MaxRecursion, 
                Method, MinRecursion, PrecisionGoal, SingularityDepth, 
                WorkingPrecision}; 
optlist = 
              Table[{opt}[\([i]\)] // First, {i, 1, 
                  Length[{opt}]}]; 
opt2 = 
              Complement[optlist, opt1]; 
If[
              Length[opt2] > 0, 
              Table[Message[\ NFourierLocalTransform::optx, 
                  opt2[\([i]\)]], {i, 1, Length[opt2]}]; 
              Break[]]; 

{prec, {a, b}, 
                symbsol} = \({PrecisionGoal, FourierParameters, 
                    Symbolic} /. {opt}\) /. 
                Options[\ NFourierLocalTransform]; 
If[
              prec === Automatic, 
              prec = $MachinePrecision - 11]; 
norm = 
              Sqrt[NIntegrate[
                  Evaluate[
                    Abs[ff]\^2], {x, \(-∞\), \(-100\), \(-50\), 0, 
                    50, 100, \(+∞\)}, 
                  Evaluate[optI]]]; 
f = 
              ff/norm; 
s = 
              segment[f, \((\(-Log[10, 4.5]\) + prec)\), x, 
                opt]; 
m = 
              Rationalize[
                Max[23/10, 
                  400\/\(\((bb - aa)\) \((s[\([2]\)] - s[\([1]\)])\)\)]]; 
nn = 
              Ceiling[m \((bb - 
                      aa)\)]; 
ss = \@\(Abs[b]\/\((2  π)\)\^\(1 - a\)\)\ \ \ \ nFT[
                  f\ Exp[I\ \(\(aa + bb\)\/2\) 
                        x], {x -> \((λ - \(aa + bb\)\/2)\)/\((π\ b)\)}, {m, nn}, s, opt]; 
normnew = 
              Abs[1 - Sqrt[
                    NIntegrate[
                      Evaluate[Abs[ss]\^2], {λ, 
                        aa, \((\(aa + bb\)\/2 + aa)\)/
                          2, \(aa + bb\)\/2, \((\(aa + bb\)\/2 + bb)\)/2, 
                        bb}, Evaluate[
                        optI]]]]; 
Print["Warning! If you extend the output function by zero out of the interval"]
\ \ \ \ \ \ \ \ \ Print["(", aa, ",", 
                bb, ") then  you will derive a global relative error (by L_2 norm) about ->", \ 
                Max[normnew, 10\^\(-prec\) // N]]; 
ssnew[
                t_] = \((ss /. λ -> 
                      t)\)\ norm; 
If[
              symbsol, 
Format[ssnew] := 
                SequenceForm["Symbolic[", {{aa, 
                      bb}}, ","<>"]"], 
Format[
                  ssnew] := 
                SequenceForm["FunctionInterpolation[", {{aa, 
                      bb}}, ","<>"]"]
]; 
ssnew
];\)

 \(NInverseFourierLocalTransform[ff_, 
            x_ -> {λ_, aa_, bb_}, 
            opt___] := 
Module[{norm, normnew, f, s, a, b, 
              nn, m, prec, ss, opt1, opt2, 
              optlist}, 

optI = 
              FilterOptions[NIntegrate, 
                opt]; 
opt1 = {InterpolationPoints, 
                InterpolationOrder, FourierParameters, Symbolic, 
                AccuracyGoal, Compiled, GaussPoints, MaxPoints, MaxRecursion, 
                Method, MinRecursion, PrecisionGoal, SingularityDepth, 
                WorkingPrecision}; 
optlist = 
              Table[{opt}[\([i]\)] // First, {i, 1, 
                  Length[{opt}]}]; 
opt2 = 
              Complement[optlist, opt1]; 
If[
              Length[opt2] > 0, 
              Table[Message[\ NInverseFourierLocalTransform::optx, 
                  opt2[\([i]\)]], {i, 1, Length[opt2]}]; 
              Break[]]; 

{prec, {a, b}, 
                symbsol} = \({PrecisionGoal, FourierParameters, 
                    Symbolic} /. {opt}\) /. 
                Options[\ NFourierLocalTransform]; 
If[
              prec === Automatic, 
              prec = $MachinePrecision - 11]; 
norm = 
              Sqrt[NIntegrate[
                  Evaluate[
                    Abs[ff]\^2], {x, \(-∞\), \(-100\), \(-50\), 0, 
                    50, 100, \(+∞\)}, 
                  Evaluate[optI]]]; 
f = 
              ff/norm; 
s = 
              segmentI[f, \((\(-Log[10, 4.5]\) + prec)\), x, 
                opt]; 
m = 
              Rationalize[
                Max[25/10, 
                  400\/\(\((bb - aa)\) \((s[\([2]\)] - s[\([1]\)])\)\)]]; 
nn = 
              Ceiling[m \((bb - 
                      aa)\)]; 
ss = \@\(Abs[b]\/\((2  π)\)\^\(1 + a\)\)\ \ \ \ nFT[
                  f\ Exp[I\ \(\(aa + bb\)\/2\) 
                        x], {x -> \((λ - \(aa + bb\)\/2)\)/\((π\ b)\)}, {m, nn}, s, opt]; 
normnew = 
              Abs[1 - Sqrt[
                    NIntegrate[
                      Evaluate[Abs[ss]\^2], {λ, 
                        aa, \((\(aa + bb\)\/2 + aa)\)/
                          2, \(aa + bb\)\/2, \((\(aa + bb\)\/2 + bb)\)/2, 
                        bb}, Evaluate[
                        optI]]]]; 
Print["Warning! If you extend the output function by zero out of the interval"]
\ \ \ \ \ \ \ \ \ Print["(", aa, ",", 
                bb, ") then  you will derive a global relative error (by L_2 norm) about ->", \ 
                Max[normnew, 10\^\(-prec\) // N]]; 
ssnew[
                t_] = \((ss /. λ -> \(-t\))\)\ norm; 
If[
              symbsol, 
Format[ssnew] := 
                SequenceForm["Symbolic[", {{aa, 
                      bb}}, ","<>"]"], 
Format[
                  ssnew] := 
                SequenceForm["FunctionInterpolation[", {{aa, 
                      bb}}, ","<>"]"]
]; 
ssnew
];\)
End[];
EndPackage[]

Introduction and implementation Notes

The package  NFourierLocalTransform is intended for calculating  approximate Fourier Transform which is one of the main mathematical tools both in theoretical investigations and in applications to applied mathematics, physics, engineering etc. 
	 It is a well known  problem to find an efficient algorithm for Fourier Transform and up to now it was solved only on the base of calculation  on a finite set of a real line and consequent interpolation. But this way is not efficient especially for oscillating inputs. Our algorithm works on the base of wavelet technique and generates rather good results for wide class of  functions (see examples below).
	It is incredible to find numerically all values on a whole real line for  Fourier Transform of a general type function .  So our algorithm works only for the functions which  values can be considered zero out of  some interval (making a given precision by (Subscript[L, 2])norm) and output gives on a reqired finite interval with  requested precision.  However it informs a user  about (Subscript[L, 2])- norm of a global error for generated output if we extend the output function by zero out of the requested interval.

◼  NFourierLocalTransform[expr, x->{t,tmin,tmax}] gives an approximate continuous Fourier transform of  expr  of  independent variable x,  as a function from t on the interval from tmin to  tmax with required precision.
◼ The Fourier transform of a function f(x) is by default defined to be \(TraditionalForm\`f\&∧ = \(1\/\@\(2  π\)\) \(∫\_\(-∞\)\%∞\( f(
                x)\)\ \(e\^\(i\ t\ x\)\) dx\)\).  Other definitions are used in some scientific and technical fields.  Different choices of definitions can be specified using the option FourierParameters. With the setting FourierParameters->{a, b } the Fourier transform computed by  NFourierLocalTransform is \(TraditionalForm\`\(\@\(\(\(|\)\(b\)\(|\)\)\/\((2  π)\)\^\(1 - a\)\)\) \(∫\_\(-∞\)\%∞\( f(
              x)\)\ \(e\^\(i\ b\ t\ x\)\) dx\)\).  Some common choices for {a, b }are {0, 1} (default; modern physics), {1, -1} (pure mathematics; systems engineering), {-1, 1} (classical physics), {0, -2 Pi} (signal processing). 
◼  NFourierLocalTransform[expr, x->{t,tmin,tmax}] yields an expression depending on the continuous variable t∈[tmin,tmax]] that represents the continuous  Fourier transform of expr with respect to the continuous variable x.  The output is InterpolatingFunction (default) or Symbolic object. 
◼  By (Subscript[L, 2])norm of a function g(t) on the interval (a, b) we consider the following integral \(\@\(\(1\/\(b - a\)\) ∫\_a\%b | g \((t)\)\( | \^2\)dt\)\).
◼ The following options can be given in addition to options of  NIntegrate:

FourierParameters                       {0, 1}                              form of  Fourier transform
Symbolic                                             False                               form of output 
(SuperStar[InterpolationPoints])                  Automatic                      n  umber of   points that the Interpolation  function use in the program
(SuperStar[InterpolationOrder])                    Automatic                      the degree of the polynomial curves of the output

\(\(\^*\) This\ option\ doesn' 
      t\ \ available\ when\ \ the\ value\ of\ the\ option\ Symbolic\ is\ \(\(True\)\(.\)\)\)

◼ The degree of the polynomial curves of the output  is specified by the option InterpolationOrder. The default setting is InterpolationOrder -> 3. 
◼ The default setting for InterpolationPoints is 20.

Demonstration of the package

All examples were solved on the computer
 PENTIUM  4,  1.5 GHz, 512 RAM

Two examples with Direct and Inverse Fourier Transforms

1.1 (Oscillating function)

Here is a function whose Fourier Transform is seeking

f[x_] = ((1/2)/(1 + x^6)) (E^(4  I  x)) + 
          (I/(2  Cosh[x - 6]));

Here are the plots of the functions  Re[f[x]] and  Im[f[x]]

Plot[Evaluate[{Re[f[x]], Im[f[x]]}], {x, -7, 12}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

Here is the Fourier transform of  f

ff[
          x_] = (1/2) 
            g[x + 4] + (I/2) (Sqrt[(π/2)]) Sech[((π x)/2)] 
            (E^(I x 6));

where

g[ω_] := 
        (1/3) (E^(-Abs[ω])) (Sqrt[(π/2)]) + 
          (1/3) (((I (-3 I + (Sqrt[3])) (E^(-(1/2) I (I - (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[-ω])/(-I + (Sqrt[3]))) - ((I (3 I + (Sqrt[3])) (E^(-(1/2) I (I + (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[-ω])/(I + (Sqrt[3]))) - ((I (3 I + (Sqrt[3])) (E^(-(1/2) I (-I - (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[ω])/(I + (Sqrt[3]))) + ((I (-3 I + (Sqrt[3])) (E^(-(1/2) I (-I + (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[ω])/(-I + (Sqrt[3]))));

Here are the plots of the functions Re[ff[x]] and  Im[ff[x]]

Plot[Evaluate[{Re[ff[x]], Im[ff[x]]}], {x, -12, 5}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

Both functions f and ff  are rather sharp changed and oscillated.  This calculates the Fourier transform of the function f  on the interval (-10,5) with the precision (10^(-4)).  The output is in Symbolic form. The message informs about the global (Subscript[L, 2])error if we extend the output function by zero out of the interval (-10,5)

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -10, 5}, 
          PrecisionGoal -> 4, Symbolic -> True] // 
        Timing

Here is the plot of the absolute error

Plot[Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -10, 5}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

Here is the actual (Subscript[L, 2])error of the output

Sqrt[
      NIntegrate[
          Evaluate[Abs[sf[[2]][x] - ff[x]]^2], {x, -10, 5}]/
        15]

The output can be given also (default) as FunctionInterpolation object

tf1 =  
      NFourierLocalTransform[f[x], x -> {λ, -10, 5}, 
          PrecisionGoal -> 4] // Timing

Here is the plot of the  absolute error. The precision of the output is unsatisfactory as the function is high-oscillating and the values 
of the options InterpolationPoints and InterpolationOrder are small

Plot[Abs[tf1[[2]][x] - ff[x]], {x, -10, 5}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

We can correct the  precision of the output by the choice of the corresponding options

tf2 =  
      NFourierLocalTransform[f[x], x -> {λ, -10, 5}, 
          PrecisionGoal -> 4, InterpolationPoints -> 60, 
          InterpolationOrder -> 7] // Timing

Here is the plot of the absolute error  which is similar to the symbolic one

Plot[Abs[tf2[[2]][x] - ff[x]], {x, -10, 5}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

Now we continue  the function sf[[2]][x] by zero out of the interval  (-10,5)

sfnew[x_] = 
        Which[x < -10, 0, x > 5, 
          0, -10 ≤  x ≤ 5, sf[[2]][
            x]];

Now we are calculating the inverse Fourier transform of  sfnew[x]. Thus we are returned  to initial image f using approximate values of ff generated above. Here is the graph of the error

stf =  
      NInverseFourierLocalTransform[sfnew[x], 
          x -> {λ, -6, 1}, PrecisionGoal -> 4, 
          Symbolic -> True] // Timing

Plot[Evaluate[Abs[stf[[2]][x] - f[x]]], {x, -6, 1}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

1.2 (Simple example)

Here is a function whose Fourier Transform is seeking

f[x_] = Sinh[((π x)/4)]/Cosh[π x];

ff[
          x_] = (((((-1))^(1/8)) ((-1 + I) + (Sqrt[2])) Sinh[(x/2)])/(((Sqrt[2 π]) + 2 (Sqrt[π]) Cosh[x])));

This calculates the Fourier transform of the function f[x] on the interval (-30,30) with the precision (10^(-13))

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -30, 30}, 
          PrecisionGoal -> 13, Symbolic -> True] // 
        Timing

Here is the plot of the absolute error

Plot[Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -30, 30}, 
      PlotRange -> All, PlotStyle -> Hue[ .6]]

Now show that our output is the global Fourier Transform of the function f with rather high precision. We continue  the function sf by zero out of the interval  (-30,30)

sfnew[x_] = 
        If[Abs[x] ≤ 30, Evaluate[sf[[2]][x]], 
          0];

Now we are calculating  the inverse of the Fourier transform for the output function

stf =  
      NInverseFourierLocalTransform[sfnew[x], 
          x -> {λ, -2, 1}, PrecisionGoal -> 8, 
          Symbolic -> True] // Timing

Here is the graph of the actual  error

Plot[Evaluate[Abs[stf[[2]][x] - f[x]]], {x, -2, 1}, 
      PlotRange -> All, PlotStyle -> Hue[ .6]]

Here is the actual (Subscript[L, 2])error

Sqrt[
      NIntegrate[
        Evaluate[Abs[stf[[2]][x] - f[x]]^2], {x, -2, 
          1}]]

Non smooth inputs

2.1 (The derivative of the function f has jumps at  the points x=±1)

Here is a function whose Fourier Transform is seeking

f[x_] = 
        If[Abs[x] ≤  1, (2 - (x^2)) (E^x), 0];

This function is continuous but the first derivatie has jumps on the points x=±1.Here is an exact Fourier transform of  f

ff[
          t_] = -((I (E^(-1 - I t)) (-3 - 
                    (t^2) + (E^(2 + 2 I t)) (-1 - 4 I t + 
                          (t^2))))/((Sqrt[2 π]) (((-I + t))^3)));

Here are the plots of the functions Re[ff[x]], Im[ff[x]]

Plot[{Re[ff[x]], Im[ff[x]]}, {x, -22, 22}, PlotRange -> All, 
      PlotStyle -> {RGBColor[0.996109,  0,  0], 
          RGBColor[0,  0,  0.996109]}]

This calculates the Fourier transform of the function f[x] on the interval (-5,11) with the precision (10^(-15))

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -5, 11}, 
          PrecisionGoal -> 15, Symbolic -> True] // 
        Timing

Here is the plot of the absolute error

Plot[Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -5, 11}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

Here is the actual (Subscript[L, 2 ])error

Sqrt[
      NIntegrate[
          Evaluate[Abs[sf[[2]][x] - ff[x]]^2], {x, -5, 11}]/
        16]

2.2 (The function has jumps at the points x=±1 and x=0)

Here is the function whose Fourier Transform is seeking

f[x_] = 
        Which[x ≤  -1,  0, -1 < x < 0, 2, 
          0 ≤  x, (E^(-x))];

Here is the plot of f

Plot[f[x], {x, -2, 4}, PlotRange -> All, 
      PlotStyle -> Hue[0.6]]

Here are an exact Fourier transform of f with plot of the functions Re[ff[x]], Im[ff[x]]

ff[
          t_] = -((I (E^(-I t)) (-2 (I + 
                          t) + (E^(I t)) (2 I + t)))/((Sqrt[2 π]) t (I + 
                    t)));

Plot[{Re[ff[x]], Im[ff[x]]}, {x, -22, 22}, PlotRange -> All, 
      PlotStyle -> {RGBColor[0.996109,  0,  0], 
          RGBColor[0,  0,  0.996109]}]

This calculates the Fourier transform of the function f[x] on the interval (-11,1) with the precision (10^(-11))

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -11, 1}, 
          PrecisionGoal -> 11, Symbolic -> True] // 
        Timing

Here is the plot of the absolute error. The actual error is less from the required due to the jumps of the function. It is impossible to achieve more precision for this example

Plot[Evaluate[Im[sf[[2]][x] - ff[x]]], {x, -11, 1}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

2.3 (The first derivative of the function has jumps at the points x=3 and x=0)

Here is a function whose Fourier Transform is seeking

f[x_] = Abs[x - 3] (E^(-x Sign[x]/3));

Here are the plots of the functions Re[f[x]] and  Im[f[x]]

Plot[Evaluate[{Re[f[x]], Im[f[x]]}], {x, -20, 26}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

Here is the Fourier transform of  f

ff[t_] := (1/(Sqrt[2 π])) (-((9 (E^(I (I + 3 t))))/(((I + 3 t))^2)) + ((27 I t)/(((I + 3 t))^2)) - ((9 (E^(-1 + 3 I t)))/(((I + 3 t))^2)) + ((9 (Cos[2 ArcTan[3 t]] - I (I + 3 t + Sin[2 ArcTan[3 t]])))/(1 + 9 (t^2))))

Here are the plots of the functions Re[ff[x]] and  Im[ff[x]]

Plot[Evaluate[{Re[ff[x]], Im[ff[x]]}], {x, -4, 8}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

This calculates the Fourier transform of the function f[x] on the interval (-24,24) with the precision (10^(-8)). The output is in Symbolic form

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -24, 24}, 
          PrecisionGoal -> 8, Symbolic -> True] // 
        Timing

Here is the plot of the absolute error

Plot[Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -24, 24}, 
      PlotRange -> All, PlotStyle -> {Hue[ .6], Hue[1]}]

Here is the actual (Subscript[L, 2 ])error

Sqrt[
      NIntegrate[
          Evaluate[Abs[sf[[2]][x] - ff[x]]^2], {x, -24, 24}]/
        48]

Comparison with MATHEMATICA current program NFouriertransform

3.1 Calculation of the Fourier Transform in a very small interval

Consider the example from item 1.1 again.

f[x_] = ((1/2)/(1 + x^6)) (E^(4  I  x)) + 
          (I/(2  Cosh[x - 6]));

ff[
          x_] = (1/2) 
            g[x + 4] + (I/2) (Sqrt[(π/2)]) Sech[((π x)/2)] 
            (E^(I x 6));

where

g[ω_] := 
        (1/3) (E^(-Abs[ω])) (Sqrt[(π/2)]) + 
          (1/3) (((I (-3 I + (Sqrt[3])) (E^(-(1/2) I (I - (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[-ω])/(-I + (Sqrt[3]))) - ((I (3 I + (Sqrt[3])) (E^(-(1/2) I (I + (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[-ω])/(I + (Sqrt[3]))) - ((I (3 I + (Sqrt[3])) (E^(-(1/2) I (-I - (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[ω])/(I + (Sqrt[3]))) + ((I (-3 I + (Sqrt[3])) (E^(-(1/2) I (-I + (Sqrt[3])) ω)) (Sqrt[(π/6)]) UnitStep[ω])/(-I + (Sqrt[3]))));

This calculates the Fourier transform of the function f  on the interval  (-8.000001,-8)  with the precision (10^(-14)) (actually we are calculating in the  asinle point).  The output is in Symbolic form

sf =  
      NFourierLocalTransform[f[x], 
          x -> {λ, -8.000001, -8}, 
          PrecisionGoal -> 13, Symbolic -> True] // 
        Timing

The message say that the output is only  0.0004% of the gloal Fourier Transform (abot a point). Here is the plot of the absolute error

Plot[
      Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -8.000001, -8}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

Now we are calculating the Fourier transform of f in the single point x=-8 by the program  NFourierTransform.

<< Calculus`FourierTransform`

The same settings do not work in required precision(10^(-14)).

NFourierTransform[Evaluate[f[x]], x, -8, 
          PrecisionGoal -> 13] - ff[-8. ] // Timing

We choose the value of the option WorkingPrecision to  acheive the required precision . The time is a bit more.

NFourierTransform[Evaluate[f[x]], x, -8, PrecisionGoal -> 13, 
          WorkingPrecision -> 85] - ff[-8. ] // Timing

Conclusion: NFourierLocalTransform works for single values of Fourier Transform not 
more worse then NFourierTransform

3.2 Calculation of the Fourier Transform in a long interval

Consider the following oscillating function

f[x_] = (1/Cosh)[x - 17];

Here is the exact Fourier transform

ff[x_] = (Sqrt[(π/2)]) Sech[((π x)/2)] 
          (E^(I x  17));

This calculates Fourier transform on the interval (-35,19) with default options

sf =  
      NFourierLocalTransform[f[x], x -> {λ, -75, 9}, 
          Symbolic -> True] // Timing

Here is the plot of the graph of the  error

Plot[Evaluate[Abs[sf[[2]][x] - ff[x]]], {x, -75, 9}, 
      PlotRange -> All, PlotStyle -> {Hue[1]}]

So actually we have the global output with the error about  (10^(-5)).
    Now we are calculating the Fourier transform by the program  NFourierTransform on a given list of outputs and then generating a function by  an  Interpolation.  After  9 seconds of work (our program  worked 3 seconds)  the precision is only (10^(-1)) (our program give precision (10^(-5))).

l = 
        Table[-75.  + 84  (k/850.) , {k, 0, 
            850}];
{t1, lf} = 
        Map[NFourierTransform[f[x], x, #] &, l] // 
          Timing;
{t2, num} = 
        Interpolation[Transpose[{l, lf}], InterpolationOrder -> 11] // 
          Timing;
Print["time->", t1 + t2];

Here is the plot of the graph of the  error

Plot[Abs[num[x] - ff[x]]^2, {x, -75, 9}, 
      PlotRange -> All]

Here is the plot of the graph of the  error which is more accurate on the interval (-15,9)

Plot[Abs[num[x] - ff[x]]^2, {x, -15, 9}, 
      PlotRange -> All]

Conclusion: NFourierLocalTransform  hase  explicite advantage over NFourierTransform 
               even for very easy input functions(when input is both  from (C^∞) and decreasing as (E^(-\[LeftBracketingBar]x\[RightBracketingBar])))