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

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;
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?

• Just for curiosity, in which situation will this technique make the frequency resolution higher? I just tested with a Sin[Range[50.]] list, but only find some small spurious peak generated, the shape of main peak seems to remain unchanged. – xzczd Jun 5 '18 at 8:55
• @xzczd I was looking at Table[Sin[2*Pi*f*t], {t, 0, 1, inc}] where inc is very small (say ~0.005). The table looks nearly continuous and the frequency can be determined by eye. However, for f=1.5, for example, or any non-integer f, Fourier does not fully capture the peak. – Questino Jun 5 '18 at 13:28
• So, this technique is useful when the original frequency resolution isn't that good? – xzczd Jun 5 '18 at 15:56
• Yes, I think so. It might be useful given a data set which is short (in the time domain, for example), but which has a high sampling rate. Fourier gives low resolution plots, even though more information is contained in the data, because Fourier only samples integer multiples of the fundamental frequency 1/T, where T is the length of the data in time. – Questino Jun 5 '18 at 18:45

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= {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= {0., Null}

Out= 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).

• Because of the improved scaling with n, this reduces my runtime from 8hrs to 0.9s, so thanks! The key philosophy change from my approach and attempts with Compile is that you simply apply the usual Fourier but modify it by shifting the exponential. Very smooth. – Questino Jun 5 '18 at 21:35
• A heterodyne Fourier transform. I used this method, with a big frequency shift and decimation, to compute fragments of long transforms on a minicomputer for my MIT Ph.D. thesis in 1978. I think the padding method is simpler for interpolation. – John Doty Jun 6 '18 at 14:56
• @JohnDoty (1) Some testing also indicates the zero padding method is a hair faster, if you use machine double zeros. (2) I have an unrelated question for you which I will raise off-forum. (3) I do recall creating beats from near-equal signal pairs, but the term "heterodyne" was new to me. (4) Not many of us still around who remember minicomputers (let alone worked with them). – Daniel Lichtblau Jun 6 '18 at 19:08

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.

• If you use powers of 2^2, wouldn't that make it "fourier"? (with apologies to Neal Sloane) – Daniel Lichtblau Jun 5 '18 at 3:40