How to make this higher resolution transform as fast as Fourier?

Question

Given a table of data, the Fourier command gives the discrete frequency content of the data, where the resolution in frequency space is equal to 1/N, where N is the number of data points.

Higher resolution can be obtained with the same set of data by writing a custom Fourier command which samples frequencies in between those sampled by the Mathematica command. However, a naive application of this is very slow. For example, the following, for sets of real-valued data parametrized by k and m,

{kmax, mmax} = {2, 2};
Do[data[k][m] = Table[RandomReal[], {i, 1, 200}], {k, 1, kmax}, {m, 1, mmax}];
lenn = Length@data[1][1];
frac = 1/2;
fou = Table[Exp[2*Pi*I*(r - 1)*(s - 1)/lenn], {r, 1, lenn}];
Do[myFourier[k][m]=Table[Abs[Total[data[k][m]*fou]], {s, 1, lenn, frac}], {k, 1, kmax}, {m, 1, mmax}]; //Timing

takes about 1.4s to run on my computer. This is ~1000 times slower than

Do[Fourier[data[k][m]], {k, 1, kmax}, {m, 1, mmax}]; // Timing

which takes 0.0013s.

My understanding is that the Mathematica Fourier command is already optimized, and so I have attempted to Compile this code in various ways, but have not been able to find a successful implementation.

The question is, why is this simple table multiplication so slow, and can Compile be used to increase its efficiency?

Answer 1

Substantial speed improvement can be made by using machine arithmetic throughout. For example, could do as below.

fourierVec[s_, n_] := Exp[2.*Pi*I*Range[0., n - 1]*(s - 1)/n]

n = 2*10^2;
data = RandomReal[1, n];

No interpolation in this experiment, just comparing time and values to Fourier.

Timing[
 four1 = 1/Sqrt[N[n]]*Table[data.fourierVec[s, n], {s, n}];]
Timing[four2 = Fourier[data];]
Max[Abs[four1 - four2]]

(* {0.015625, Null}

{0., Null}

8.04419*10^-14 *)

Now compute with the half-stride.

Timing[
 myFourier1 = 1/Sqrt[N[n]]*Table[data.fourierVec[s, n], {s, 1., n, .5}];]

(* Out[124]= {0.03125, Null} *)

Not so bad. But we are running an O(n^2) algorithm and that is simply not going to scale like Fourier at O(n log n). Can we recover this result using just Fourier? Yes, if we compute an "offset" FT and riffle with the regular FT.

Timing[
 myFourier2 = 
   Riffle[Fourier[data], Most@Fourier[fourierVec[1.5, n]*data]];]
Max[Abs[myFourier1 - myFourier2]]

(* Out[132]= {0., Null}

Out[133]= 8.3098*10^-14 *)

This can be extended in a fairly straightforward manner to other fractions. It appears to be substantially the same as the interpolation variant given by @JohnDoty, other than the issue of dropping the last value (which I did in order to get agreement with the corrected DFT from the original post).

Answer 2

Fourier uses the Fast Fourier Transform (FFT), much faster than a direct method. But you can easily create what you want just by padding the data with zeros, since the delta frequency is inversely related to the array length.

dat = RandomReal[1, 10];
Fourier[dat]
(* {1.53116 + 0. I, -0.0783886 - 0.207404 I, -0.298173 + 
    0.179654 I, -0.0245808 + 0.0339534 I, -0.0390086 + 
    0.501077 I, -0.075875 + 0. I, -0.0390086 - 0.501077 I, -0.0245808 - 
    0.0339534 I, -0.298173 - 0.179654 I, -0.0783886 + 0.207404 I} *)

Interpolate with a padded FFT. n interpolation factor. The Sqrt is to get the normalization the same as the unpadded:

interpolatedFourier[dat_, n_] := 
   Fourier[Join[dat, ConstantArray[0, (n - 1) Length[dat]]]] Sqrt[n]
interpolatedFourier[dat, 2]
(* {1.53116 + 0. I, -0.022339 + 1.03333 I, -0.0783886 - 0.207404 I, 
    0.339664 + 0.429062 I, -0.298173 + 0.179654 I, 
    0.0707628 - 0.0261038 I, -0.0245808 + 0.0339534 I, 
    0.279942 + 0.0442931 I, -0.0390086 + 0.501077 I, -0.380536 + 
    0.0354403 I, -0.075875 - 1.97484*10^-17 I, -0.380536 - 
    0.0354403 I, -0.0390086 - 0.501077 I, 
    0.279942 - 0.0442931 I, -0.0245808 - 0.0339534 I, 
    0.0707628 + 0.0261038 I, -0.298173 - 0.179654 I, 
    0.339664 - 0.429062 I, -0.0783886 + 0.207404 I, -0.022339 - 
    1.03333 I} *)

Note that with n==2, every other point matches the uninterpolated version.

Added note: if you're really after speed, the FFT is most efficient when the length of the array is a power of two. So, it may be useful to pad out to such a length. The resulting transform will have increased resolution, but will not sample the same frequencies the unpadded FFT sampled unless its length was a power of 2. Details are left to the reader.