c - 在VBA Excel中实现Rosetta Code FFT

我试图将FFT Rosetta Code实现到VBA excel中我无法完全按照Rosetta代码中所写的那样重建相同的输出数据起初我以为这是类型转换不匹配，但将输入值缩放1.1也将输出值缩放1.1我认为我的代码唯一有问题的地方是我如何在代码中实现引用的数组，而不是Rosetta所写的我发现Rosetta通过在其奇数递归调用中写出+step和buf+step来改变数组的地址我不知道VBA中有类似的构造，所以我只是传递了一个额外的递归参数shift，它跟踪在传递到下一个递归调用时地址的移动。
我的轮班执行有什么问题，或者是其他什么问题？
Rosetta声称它的FFT是内存到位的，但是他们制作了数据的额外副本并将其存储到out[]如果你能解释一下为什么会这样。
C语言中的FFT-Rosetta码

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);
}

我在VBA中尝试FFT

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

输出比较

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)

最佳答案

您已声明c1为Long：

 "Dim c1(1 To 2) As Long"

将其更改为Double并且它可以工作：

 "Dim c1(1 To 2) As Double"