I attempted to implement the FFT Rosetta Code into VBA excel. I was unable to reconstruct the same output data exactly as written in the Rosetta Code. At first I thought it was type conversion mismatch, but scaling the input values by 1.1 scaled the output values by 1.1 as well. The only aspect of my code that I think could be in question is how I implemented the referenced arrays in code versus what Rosetta writes. I discovered Rosetta shifts the addresses of the arrays by writing out + step and buf + step in its odd recursion calls. I am unaware of a similar construction in VBA, so I simply passed an additional recursion parameter shift, which keeps track of the shifting of addresses as its passed into the next recursion call.
What is wrong with my shift implementation or is it something else?
Rosetta claims its FFT is memory in place, however they make an extra copy of the data and store it into out[]. I would appreciate an explanation of why this is so.
FFT Rosetta Code in C
void _fft(cplx buf[], cplx out[], int n, int step)
{
if (step < n) {
_fft(out, buf, n, step * 2);
_fft(out + step, buf + step, n, step * 2);
for (int i = 0; i < n; i += 2 * step) {
cplx t = cexp(-I * PI * i / n) * out[i + step];
buf[i / 2] = out[i] + t;
buf[(i + n)/2] = out[i] - t;
}
}
}
void fft(cplx buf[], int n)
{
cplx out[n];
for (int i = 0; i < n; i++) out[i] = buf[i];
_fft(buf, out, n, 1);
}
My attempted FFT in VBA
Private Sub rec_fft(ByRef buf, ByRef out, ByVal n, ByVal step, ByVal shift)
If step < n Then
Dim ii As Long
Dim pi As Double
pi = 3.14159265358979 + 3.23846264338328E-15
Dim c1(1 To 2) As Long
Dim c2(1 To 2) As Double
Call rec_fft(out, buf, n, step * 2, shift)
Call rec_fft(out, buf, n, step * 2, shift + step)
For i = 1 To n / (2 * step)
ii = 2 * step * (i - 1)
c1(1) = Cos(-1 * pi * CDbl(ii) / CDbl(n))
c1(2) = Sin(-1 * pi * CDbl(ii) / CDbl(n))
c2(1) = c1(1) * out(shift + 1 + ii + step, 1) - c1(2) * out(shift + 1 + ii + step, 2)
c2(2) = c1(1) * out(shift + 1 + ii + step, 2) + c1(2) * out(shift + 1 + ii + step, 1)
buf(shift + 1 + ii / 2, 1) = out(shift + 1 + ii, 1) + c2(1)
buf(shift + 1 + ii / 2, 2) = out(shift + 1 + ii, 2) + c2(2)
buf(shift + 1 + (ii + n) / 2, 1) = out(shift + 1 + ii, 1) - c2(1)
buf(shift + 1 + (ii + n) / 2, 2) = out(shift + 1 + ii, 2) - c2(2)
Next i
End If
End Sub
Private Sub fft(ByRef buf, ByVal n)
Dim out() As Double
ReDim out(1 To n, 1 To 2)
For i = 1 To n
out(i, 1) = buf(i, 1)
out(i, 2) = buf(i, 2)
Next i
Call rec_fft(buf, out, n, 1, 0)
End Sub
Output Comparison
Input Data: 1 1 1 1 0 0 0 0
Rosetta FFT : 4 (1, -2.41421) 0 (1, -0.414214) 0 (1, 0.414214) 0 (1, 2.41421)
My FFt : 4 (1, -3) 0 (1, -1) 0 (1, 1) 0 (1, 3)
For i = 1 To n Step (2 * step)
(now you see a potential problem use a variable named "step"). The third is the integer division in the array indexing; instead of (/2) use (\2). There may be other issues, but those are the ones that I saw.