Fast Fourier Transform 
on Nonequispaced Grid

The Package   NUFourier

(*:Summary:This package calculates the list  (∑_(k = 1)^N) (Subscript[f, k]) 
                    (E^(I π (Subscript[t, k]) (Subscript[w, j]))), j = 1,...M  with a given precision for the given complex valued listf={(Subscript[f, k])}, real list t={(Subscript[t, k])},k = 1,2,..,N and real list w={(Subscript[w, k])}, k = 1,2,..,M. *)
  
  
  
  
  

         

	  

  


BeginPackage["NumericalMath`NUFourier`"];
NUFourier::incsh = "Lists have incompatible shapes.";
NUFourier::warn="NUFourier[f,{t,w}] works best if (max{(Subscript[t, k])}-min{(Subscript[t, k])})(max{(Subscript[w, k])}-min{(Subscript[w, k])})≤20 N, where N=Length[f]. Otherwise the program works slow .";
Off[General::"spell"]; 

NUFourier::usage="NUFourier[f,{t,w}] evaluates the sum (∑_(k = 1)^N)(Subscript[f, k])\!\(E\^\(I\\ π\\ \\ \(t\_k\) w\_j\)\), j=1,...M with a given precision for the given complex valued list f={(Subscript[f, k])}, real list t={(Subscript[t, k])}, k=1,2,..,N and real list w={(Subscript[w, k])}, k=1,2,..,M.";
Options[
          NUFourier] = {PrecisionGoal -> 
            Automatic};
Begin["`Private`"];
dt[n_Integer, al_] := 
        Table[ Evaluate[(E^(al (k^2)))], {k, -Floor[n/2], 
            n - Floor[n/2] - 1}];
\(per[l_List, p_Integer] := 
	
        Block[{nn, lb, le}, 
		nn = Length[l]; 
		lb = 
            Take[l, {1, p}]; 
		le = 
            Take[l, {nn - p + 1, nn}]; 
		Flatten[{le, l, 
              lb}]
		];\)
ge[n_Integer, m_] := 
        Table[0, {k, 
            m n - 1 - Floor[ m n/2] - n + Floor[n/2] + 1}];
gb[
          n_Integer, m_] := 
        Table[0, {k, Floor[m n/2] - Floor[n/2]}];
nul[ll_List, 
          m_] := Flatten[{gb[Length[ll], m], ll, 
            ge[Length[ll], m]}];
lim[prec_] := 
        Which[prec < 9, {prec, 70}, prec < 11, {prec, 80}, 
          prec < 14, {prec, 100}, 
          True, {prec, 120}];
in = 
        Interpolation[Table[lim[prec], {prec, 0, 16}], 
          InterpolationOrder -> 1];
\(flst[l_List] := 
          Block[{fy, nn = Length[l]}, 
	fy = RotateLeft[l, Floor[nn/2]]; 

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

\(tiv[t_List, m_, qq_] := 
          Module[{nn, d, a, aa, aaa, e, b, c, nz, ss, pss}, 

            nn = 2\ m\ \ qq; 
d = Ceiling[Sort[nn\ \ t + nn]]; 

            a = Complement[Range[1, 2\ nn], d]; 

            aa = Table[0. , {Length[a]}]; 
aaa = Transpose[List[a, aa]]; 

            e = Complement[Range[1, 2\ nn], a]; 
b = Split[d]; 

            c = Table[Length[b[\([i]\)]], {i, 1, Length[b]}]; 

            nz = Transpose[List[e, c]]; 

            ss = Transpose[Sort[Join[nz, aaa]]] // Last; 

            pss = Partition[ss, 2]; 

            Table[Apply[Plus, pss[\([i]\)]], {i, 1, Length[pss]}]];\)

\(pe[al_, {f_List, t_List}, q_Integer, m_, qq_] := 

          Module[{nn, s, zero, u, usum, l, 
              ll, α, β}, 
nn = 
              Length[f]; 
s = 
              tiv[t, m, qq]; 
zero = 
              Table[0, {q}]; 
u = Join[zero, s, zero]; 

            usum = ReplacePart[FoldList[Plus[#1, #2] &, 0, u], 0, 1]; 

            l = Range[\(-m\)\ qq - q/2, m\ qq + q/2]; 

            ll = l + \ m\ qq + q/2 + 1; 
α = 
              Map[usum[\([#]\)] + 1\  &, ll]; 
β = 
              Map[usum[\([q + #]\)]\  &, 
                ll]; \ \ \ \ 

\ \(\@π\/\@al\) 
              Table[∑\+\(k = α[\([j + m\ qq + prec + 1]\)]\)\%\(β[\([j + m\ qq + prec + 1]\)]\)f[\([k]\)] 
                    E\^\(-\(\(\(π\^2\) \((2\ m\ qq\ \ t[\([k]\)] - 2\ j)\)\^2\)\/\(4\ al\)\)\), {j, \(-m\)\ qq - prec, 
                  m\ qq + 
                    prec}]


];\)

\(superA2[w_List, f_List, prec_] := 
	
          Module[{h, m, nn, al0, al2, al1, p, ce, al, a, aa, b, ll, lp, lr, 
              lf, ww, le, lu, ls, lk, lq}, 
h = 
              1; 
m = 1.8; 
nn = 
              Length[f]; 
al0 = \(h\ 2.3\ \ m\ \((prec)\)\)\/\((m - 1)\); 
al2 = 
              FindRoot[\ \@\(m\/2\ \ Sinh[\(2  α\)\/m]\)\%4 == 
                  E\^\(\(-\ 1.15\)\ prec + α/2 + Log[α]/4\), {α, al0, al0 + 1}, 
MaxIterations -> 
                  15]; 
al1 = \((\(al2 // First\) // 
                    Last)\) + 1\/\(prec\^2 + 1\); 
p = 
              2\ \((Round[\ \(π\^\(-1\)\) \@al1\ \@\(2.3\ prec + \(0.25\ al1\)\/m\^2 - Log[4  al1/π]/2\)\ \ ])\); 
		ce = 
              Ceiling[p/2]; 
\ \ \ \ \ 	al = 
              al1\/nn\^2; 
		a = π^2/\((4\ al)\); 
		aa = 
              4  a/nn\^2; 
		b = Floor[nn/2]; 
		ll = 
              Table[k, {k, 0, 
                  p - 1}]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ lp = 
              flst[f\ dt[nn, al]]; 
		lr = per[lp, ce]; 
		lf = 
              Map[Evaluate[E^\((\(-aa\)\ #\^2)\)] &, ll]; 
		ww = 
              Map[Evaluate[\(-2\)\ Ceiling[\(-\(p\/2\)\) + b\ #1] + 
                      nn\ #1] &, w]; 
		le = 
              Evaluate[Exp[\(-aa\)\ ww\^2/4]]; 
		lu = 
              Evaluate[Exp[aa\ ww]]; 
		ls = 
              Map[Evaluate[\((Ceiling[b\ #\ \  - p/2] + 1 + ce + b)\)] &, 
                w]; 
		lk = Outer[lr[\([#1 + #2]\)] &, ll, ls]; 
		lq = 
              Map[Evaluate[lu]^#\  &, ll]; 
		\(√π\/√al\) 
              Evaluate[le\ \ Apply[Plus, lf\ lq\ lk]]];\)

\(NUFourier[f_, {w_, ww_}, opt___] := 
        Block[{nn, mm, a, c, fnew, wwnew, wnew, qq, m, α, q, s2, 
            prec}, 
{prec} = \({PrecisionGoal} /. {opt}\) /. 
              Options[NUFourier]; 
If[prec === Automatic, 
            prec = $MachinePrecision - 11]; 
nn = 
            Length[f]; 
\ \ \ \ \ \ \ nnnew = 
            Length[w]; 
If[nn == nnnew, , 
            Message[NUFourier::incsh]; 
            Break[]]; \ 
If[\((Max[w] - 
                    Min[w])\) \((Max[ww] - Min[ww])\) ≤ \ 
              20\ nn, , Message[NUFourier::warn]]; 
mm = 
            Length[
              ww]; 
If[\((\(nn\^\(1 - 5\/8\)\) 
                    mm\^\(5\/8\) < in[prec])\) || \((nn ≤ 40)\), 
            Table[∑\+\(k = 0\)\%\(nn - 1\)f[\([k + 1]\)] 
                  E\^\(I\ π\ \ w[\([k + 1]\)]\ \ ww[\([j]\)]\), {j, 1, 
                mm}], 
a = \(\((w // First)\) + \((w // Last)\)\)\/2; 
c = \(\((ww // First)\) + \((ww // Last)\)\)\/2; 
fnew = 
              Table[f[\([k + 
                        1]\)]\ E\^\(I\ π\ \ w[\([k + 1]\)]\ \ c\), {k, 0, 
                  nn - 1}]; 
wwnew = \((ww - 
                    c)\) \((\((w // Last)\) - \((w // 
                        First)\))\); 
wnew = \((w - 
                    a)\)/\((\((w // Last)\) - \((w // 
                        First)\))\); 
qq = \ 
              Ceiling[\((wwnew // Last)\)]; 
m = 
              Ceiling[1.7 + 
                  4\/\(prec + 1\)]; 
α = \(m \((Log[2] 2\^\(4 - 1.5  prec\) + \ Log[10]\ prec)\)\)\/\(m - 1\); 
q = 2\ prec; 
s2 = 
              pe[α, {fnew, wnew}, q, m, 
                qq]; 
Table[\(E\^\(I\ π\ \ a\ \((ww[\([j]\)] - c)\)\)\) 
                  E\^\(α \((wwnew[\([j]\)]\/\(2\ m\ qq\))\)\^2\), {j, 
                  1, mm}]\ \ superA2[\(\(\ \)\(wwnew\)\)\/\(\(\ \)\(m\ qq\)\), nul[s2, 1.8], 
                prec]]
];\)
End[];
EndPackage[];

Introduction and implementation Notes

The package NUFourier calculates the list

(∑_(k = 1)^N) (Subscript[f, k]) 
            (E^(I π  (Subscript[t, k]) (Subscript[w, j]))), j = 1,2,..,M                                                     (1)

for given complex valued list f={(Subscript[f, k])},real list t={(Subscript[t, k])}(k=1,2,..,N) and real list w={(Subscript[w, k])}(k=1,2,..,M).  The fast calculation of this Fourier Transform on (in general)nonuniform grids is one of the important problems  in applied mathematics.It is now central to many areas, notable spectral analysis in signal processing  when the input data is not uniformly spaced,as well as for mathematical sources of the computer tomography.  There are many papers of known autors devoted to fast solution of this problem (see [1]). In [2] it  was solved  on the base of  the Gaussian Bell function  but the technique of the precision providing was not properly  developed.  The current algorithm  NUFourier[f,{t,w}] is based on the paper[3] which provides both fast calculation  and about exact required precision for the sum (1).
    
        References
[1] A.F.Ware. Fast Fourier Transforms for Irregularly Spaced Data, SIAM Review, vol 40,No 4,1998,pp.838-856.
[2] A.Dutt,V.Rokhlin,Fast Fourier Transforms for Nonequispaced Data,Siam J.Sci.Comput.,vol 14,No 6,1993, pp.1368-1393.
[3] A. Nersessian, Quasiexponential Function and Fourier Quasitransforms,Doklady NAN Armenii (Reports of the National Academy of Sciences of Armenia), vol 101,No 1,2001,pp.5-11.

Implementation notes (see above  formula (1) )

◼NUFourier[expr,{list1,list2}] gives an approximate Discrete Fourier Transform of  the expr  on the general type nonequally spaced real grid list2 when the initial Fourier nodes correspond the general type nonequally spaced real grid list2  for f=expr, t=list1 and w=list2).
◼NUFourier has the option PrecisionGoal->p with default setting  p = $MachinePrecision - 11.
◼The option PrecisionGoal corresponds to relative (Subscript[L, 2])-error of the output list.
◼If N M is not big ( ≤  about (100^2) ) or only N is not big (N≤40) NUFourier works by the direct calculation of the sum (1).  Otherwise the wavelet technique is used.
◼NUFourier works more faster if  M >> N.
◼The theoretical complexity of the direct calculation of the sum (1)  is about 2 N×M arithmetic operations. The cost of the  NUFourier  is proportional to N(LogN+p) for PrecisionGoal->p. 
◼ NUFourier works best if  N >> 1, M >> 1,  M > N  /10 and (max {(Subscript[t, k])} - min {(Subscript[t, k])}) (max {(Subscript[w, k])} - 
              min {(Subscript[w, k])}) ≤ 20 N

Demonstration of the package

\(All\ examples\ were\ solved\ on\ the\ computer

      PENTIUM\ 4, 1.5\ GHz, 512\ RAM\)

Here we calculate the sum (∑_(k = 1)^N) (Subscript[f, k]) 
          (E^(I π (Subscript[t, k]) (Subscript[(w ), j]))), j = 1, ... 2  N  for a random complex list {(Subscript[f, k])} and random real lists {(Subscript[t, k])} and {(Subscript[w, j])} and compare the times of  direct evaluation  with the time of  program  NUFourier[] when the last is working with default precision = 5.

\(Do[Block[{nn = 2^y, mm = 2\ nn, t, f, x, s1, 
          s2}, 
Print["**** M->", mm, "****N->", 
          nn]; 
t = 
          Sort[Table[\ Random[Real, {0, 2}], {nn}]]; 
f = 
          Table[2  Random[
                Complex, {\(-1\) - I, 
                  1 + 2  I}], {nn}]; 
x = 
          Sort[Table[
              Random[Real, {0, mm}], {mm}]]; 
{time1, s2} = 
          NUFourier[f, {t, x}] // Timing; 
\ {time2, s1} = 
          Table[∑\+\(k = 0\)\%\(nn - 1\)f[\([k + 1]\)] 
                  E\^\(I\ π\ \ t[\([k + 1]\)]\ \ x[\([j]\)]\), {j, 1, 
                mm}] // 
            Timing; 
Print["Time->{direct,fast,direct/fast}->", {time2, time1, 
            time2/time1}]; 
Print[\*""Relative (Subscript[L, 2])-error->"", 
          Sqrt[Apply[Plus, Abs[s1 - s2]\^2]]/
            Sqrt[Apply[Plus, Abs[s1]^2]]]
], {y, 4, 
        24}]\)

Here we calculate the previous sum omly by the program NUFourier[] for more big values of N and M

\(Do[Block[{nn = 2^y, mm = 2\ nn, t, f, x, s1, 
          s2}, 
Print["**** M->", mm, "****N->", 
          nn]; 
t = 
          Sort[Table[Random[Real, {0, 2}], {nn}]]; 
f = 
          Table[2  Random[
                Complex, {\(-1\) - I, 
                  1 + 2  I}], {nn}]; 
x = 
          Sort[Table[
              Random[Real, {0, mm}], {mm}]]; 
{time1, s2} = 
          NUFourier[f, {t, x}] // 
            Timing; 
Print["Time->{fast}->", {time1}];
], {y, 4, 24}]\)

Conclusion

From this and several other experiments we can conclude that if N and M are growing together then the time of 
NUFourier (in the last resort for large values of N and M) is much less than the time of direct calculation of the sum (1).