How does Fourier work with arbitrary FourierParameters?

Question

Mathematica's Fourier function allows you to insert an arbitrary real number in the exponent of the discrete Fourier transform, via FourierParameters, so that the transform becomes something like $$ \sum_{j=1}^n u_j e^{2\pi i b (j-1)(k-1)/n} $$ with $b$ arbitrary.

I find this very useful in my work, but I don't understand how it works. The usual description of the FFT algorithm (e.g. as outlined in "Numerical Recipes") depends in an essential way on the $n$-periodicity of exponential factor. When I insert an arbitrary $b$ and call Fourier, am I still getting a fast Fourier transform? Or is it secretly doing a slow DFT? In the former case, how does Mathematica get around the need for a periodic exponential factor?

Answer 1

Here is an example of how to use FourierParameters for interpolating a spectrum. First we generate some data, take the ordinary Fourier transform and plot.

data = Table[Sin[\[Pi] k/33], {k, 256}] // N;
ListLinePlot[data]
ft = Fourier[data, FourierParameters -> {-1, -1}];
ListLinePlot[Abs@ft[[1 ;; 128]], PlotRange -> All]
ListPlot[Abs@ft[[1 ;; 20]], PlotRange -> All]

The last, expanded view of the spectrum shows that there are rather few points defining the peak.

The following Dynamic enables the interpolation to be seen when FourierParameters are altered. The actual frequency and the magnitude of the Fourier peak should be

 (*   {0.0151515, 0.511146}   *)

DynamicModule[{b = 1, pos, a = 0, ft, freqs},
 ft = Abs@Fourier[data, FourierParameters -> {-1, b}];
 pos = Position[ft, Max[ft]][[1, 1]];
 freqs = Table[(n - 1)/Length[data], {n, Length[data]}];
 Column[{
   Dynamic[
    Row[{"Frequency (Hz) = ", a[[1]], "     Magnitude = ", a[[2]] }]],
   Dynamic@ListLinePlot[Transpose[{freqs, ft}], PlotRange -> All,
     Epilog -> {Red, PointSize[0.01], 
       Point[a = Transpose[{freqs, ft}][[pos]]]},
     ImageSize -> 12 72, Frame -> True, Mesh -> All],
   Slider[Dynamic[b, {b = #;
       ft = Abs@Fourier[data, FourierParameters -> {-1, b}];
       freqs = Table[b (n - 1)/Length[data], {n, Length[data]}];
       pos = Position[ft, Max[ft]][[1, 1]]
       } &], {0, 1}, Appearance -> "Labeled"]
   }]
 ]

As you can see we can get many more points into the peak and have thus interpolated the spectrum. Thus we can more accurately work out the frequency.

Hope that helps.

Answer 2

Here is an implementation based on Bailey's paper.

ClearAll[ffrft$set] ;
ffrft$set[length_, span_] := Block[
{trig, work},
trig =Exp[I*Pi*span/length*Range[0.0, length-1]^2] ;
work = Fourier[Flatten[{Conjugate[trig] , Exp[-I*Pi*span/length*Range[-length, -1.0]^2]}] , FourierParameters -> {1, 1} ];
{trig, work}
] ;

ClearAll[ffrft$get] ;
ffrft$get[length_, signal_, trig_, work_] := Block[
{ffrft},
ffrft = ConstantArray[0.0, 2*length]  ;
ffrft[[1;;length]] = trig*signal ;
ffrft = work*Fourier[ffrft, FourierParameters -> {1, 1} ] ;
ffrft =Fourier[ffrft, FourierParameters -> {1, -1} ][[1;;length]] ;
ffrft*trig/(2*length)
] ;

Example:

length = 1024 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
{trig, work} = ffrft$set[length, b]  ;
ref = Fourier[signal, FourierParameters -> {1, b}] ;
res = ffrft$get[length, signal, trig, work] ;
ref-res // Chop // DeleteDuplicates
(* {0} *)

Timing:

length = 8192 ;
signal = Sin[2*Pi*0.12345*Range[length]] ;
b = 1/Sqrt[2.0]  ;
t$ref = First@RepeatedTiming[Fourier[signal, FourierParameters -> {1, b}] ];
t$res = First@RepeatedTiming[{trig, work} = ffrft$set[length, b]  ; ffrft$get[length, signal, trig, work] ];
t$ref/t$res

For FourierParameters -> {1, 1} WM is most likely uses FFT from MKL, while for FourierParameters -> {a, b} it is not clear. While it is possible to compute with n*log(n), I'm not sure it is the case in WM.