在 VBA Excel 中实现 Rosseta 代码 FFT

Question

我试图在 VBA excel 中实现 FFT Rosetta Code。我无法重建与 Rosetta Code 中所写完全相同的输出数据。起初我认为这是类型转换不匹配，但将输入值缩放 1.1 也会将输出值缩放 1.1。我认为我的代码中唯一可能存在问题的方面是我如何在代码中实现引用的数组，而不是 Rosetta 编写的内容。我发现 Rosetta 通过在其奇数递归调用中写出 + step 和 buf + step 来移动数组的地址。我不知道 VBA 中有类似的结构，所以我只是传递了一个额外的递归参数 shift，它在传递到下一个递归调用时跟踪地址的移动。

我的班次实施有什么问题还是其他原因？

Rosetta 声称它的 FFT 是内存到位，但是他们制作了一份额外的数据副本并将其存储到 out[] 中。我希望能解释一下为什么会这样。

C 中的 FFT 罗塞塔代码

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)

Answer 1

您已将 c1 声明为Long：

 "Dim c1(1 To 2) As Long"

把它改成Double就可以了：

 "Dim c1(1 To 2) As Double"