2

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)
2
  • 1
    I see 3 issues. The first is the use of the variable "step". Step is a VB keyword used in For-Next loops. The second is the For-Next declaration; it should be something like 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.
    – TnTinMn
    Commented Jul 1, 2015 at 2:15
  • @TnTinMn Wow, thanks. I did not concern these options. I will try to implement them. Commented Jul 1, 2015 at 10:44

1 Answer 1

1

You have declared c1 as Long:

 "Dim c1(1 To 2) As Long"

Change it to Double and it works:

 "Dim c1(1 To 2) As Double"
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.