I tried converting some C code I found on the internet to D.
When I run the code below:
import std.math;
import std.stdio;
static double[2][] fft(in double[2][] timeDomain)
{
static void fft_rec(sizediff_t n, sizediff_t offset, sizediff_t
delta, in double[2][] timeDomain, double[2][] freqDomain, double[2][]
scratch)
{
if (n > 2) /* Perform recursive step. */
{
const N2 = n / 2; /* half the number of points in FFT */
/* Calculate two (n/2)-point DFT's. */
fft_rec(N2, offset, 2 * delta, timeDomain, scratch,
freqDomain);
fft_rec(N2, offset + delta, 2 * delta, timeDomain,
scratch, freqDomain);
/* Combine the two (n/2)-point DFT's into one n-point
DFT. */
foreach (k; 0 .. N2)
{
const double TWO_PI =
6.2831853071795864769252867665590057683943L;
const
k00 = offset + k * delta,
k01 = k00 + N2 * delta,
k10 = offset + 2 * k * delta,
k11 = k10 + delta,
cs = cos(k * TWO_PI / n),
sn = sin(k * TWO_PI / n),
tmp0 = cs * scratch[k11][0] + sn * scratch[k11][1],
tmp1 = cs * scratch[k11][1] - sn * scratch[k11][0];
freqDomain[k01][0] = scratch[k10][0] - tmp0;
freqDomain[k01][1] = scratch[k10][1] - tmp1;
freqDomain[k00][0] = scratch[k10][0] + tmp0;
freqDomain[k00][1] = scratch[k10][1] + tmp1;
}
}
else /* Perform 2-point DFT. */
{
const k00 = offset, k01 = k00 + delta;
freqDomain[k01][0] = timeDomain[k00][0] -
timeDomain[k01][0];
freqDomain[k01][1] = timeDomain[k00][1] -
timeDomain[k01][1];
freqDomain[k00][0] = timeDomain[k00][0] +
timeDomain[k01][0];
freqDomain[k00][1] = timeDomain[k00][1] +
timeDomain[k01][1];
}
}
auto freqDomain = new double[2][timeDomain.length];
fft_rec(timeDomain.length, 0, 1, timeDomain, freqDomain, new
double[2][2 * timeDomain.length]);
return freqDomain;
}
void main()
{
const result = fft([[0, 1], [1, -1]]);
pragma(msg, result);
writeln(result);
}
I get the following results:
Compile time: [[1,0],[1,0]]
Run time: [[1, 0], [-1, 2]]
Why is this?
