C++实现基3、基5混合基FFT时9点等复合长度计算结果错误求解

问题背景 #

我正在C中实现基5、基3 FFT算法，已经成功完成基2 FFT的编写，但在实现基3或基5逻辑时存在Bug：对3个采样点做FFT可以得到正确结果，但对9个采样点（3*3）做FFT时输出结果不正确，仅有首个计算值符合预期。
我最初参考的是可正常运行的Python代码，尝试直接将其迁移到C中。

参考Python代码 #

import numpy as np
def fft(x):
    """
    radix-2,3,5 FFT algorithm
    """
    N = len(x)
    if N <= 1:
        return x
    elif N % 2 == 0:
        # For multiples of 2 this formula works
        even = fft(x[0::2])
        odd =  fft(x[1::2])
        T = [np.exp(-2j*np.pi*k/N)*odd[k] for k in range(N//2)]
        return [even[k] + T[k] for k in range(N//2)] + \
               [even[k] - T[k] for k in range(N//2)]
    elif N % 3 == 0:
        # Optional, implementing factor 3 decimation
        p0 = fft(x[0::3])
        p1 = fft(x[1::3])
        p2 = fft(x[2::3])
        # 注：原代码缩进有误，已修正为列表推导式写法
        return [p0[k % (N//3)] +
                p1[k % (N//3)] * np.exp(-2j*np.pi*k/N) + 
                p2[k % (N//3)] * np.exp(-4j*np.pi*k/N) for k in range(N)]
    elif N % 5 == 0:
        #factor 5 decimation
        p0 = fft(x[0::5])
        p1 = fft(x[1::5])
        p2 = fft(x[2::5])
        p3 = fft(x[3::5])
        p4 = fft(x[4::5])

        return [p0[k % (N//5)] +
                p1[k % (N//5)] * np.exp(-2j*np.pi*k/N) + 
                p2[k % (N//5)] * np.exp(-4j*np.pi*k/N) + 
                p3[k % (N//5)] * np.exp(-6j*np.pi*k/N) +
                p4[k % (N//5)] * np.exp(-8j*np.pi*k/N)
               for k in range(N)]

x = [1,1.00071,1.00135,1.00193,1.00245,1.0029,1.00329,1.00361,1.00387]
assert(np.allclose(fft(x), np.fft.fft(x)))

存在问题的C++代码 #

fft.hpp #

#define _USE_MATH_DEFINES
#pragma once
#include <cmath>
#include <vector>
#include <complex>


using std::vector;
using std::complex;

vector<complex<float>> slicing(vector<complex<float>> vec, unsigned int X, unsigned int Y, unsigned int stride)
{
    // To store the sliced vector
    vector<complex<float>> result;

    // Copy vector using copy function()
    int i = X;
    while (result.size() < Y)
    {
        result.push_back(vec[i]);
        i = i + stride;
    }
    // Return the final sliced vector
    return result;
}

void fft(vector<complex<float>>& x)
{
    // Check if it is splitted enough

    const size_t N = x.size();
    if (N <= 1)
        return;

    else if (N % 2 == 0)
    {
        //Radix-2
        vector<complex<float>> even = slicing(x, 0, N / 2, 2); //split the inputs in even / odd indices subarrays
        vector<complex<float>>  odd = slicing(x, 1, N / 2, 2);

        // conquer
        fft(even);
        fft(odd);

        // combine
        for (size_t k = 0; k < N / 2; ++k)
        {
            complex<float> t = std::polar<float>(1.0, -2 * M_PI * k / N) * odd[k];
            x[k] = even[k] + t;
            x[k + N / 2] = even[k] - t;
        }
    }
    else if (N % 3 == 0)
    {
        //Radix-3
        //factor 3 decimation
        vector<complex<float>> p0 = slicing(x, 0, N / 3, 3);
        vector<complex<float>> p1 = slicing(x, 1, N / 3, 3);
        vector<complex<float>> p2 = slicing(x, 2, N / 3, 3);

        fft(p0);
        fft(p1);
        fft(p2);

        for (int i = 0; i < N; i++)
        {
            complex<float> temp = p0[i % (int)N / 3];
            temp += (p1[i % (int)N / 3] * std::polar<float>(1.0, -2 * M_PI * i / N));
            temp += (p2[i % (int)N / 3] * std::polar<float>(1.0, -4 * M_PI * i / N));
            x[i] = temp;
        }
    }
    else if (N % 5 == 0)
    {
        //Radix-5
        //factor 5 decimation
        vector<complex<float>> p0 = slicing(x, 0, N / 5, 5);
        vector<complex<float>> p1 = slicing(x, 1, N / 5, 5);
        vector<complex<float>> p2 = slicing(x, 2, N / 5, 5);
        vector<complex<float>> p3 = slicing(x, 3, N / 5, 5);
        vector<complex<float>> p4 = slicing(x, 4, N / 5, 5);

        fft(p0);
        fft(p1);
        fft(p2);
        fft(p3);
        fft(p4);
        for (int i = 0; i < N; i++)
        {
            complex<float> temp = p0[i % (int)N / 5];
            temp += (p1[i % (int)N / 5] * std::polar<float>(1.0, -2 * M_PI * i / N));
            temp += (p2[i % (int)N / 5] * std::polar<float>(1.0, -4 * M_PI * i / N));
            temp += (p3[i % (int)N / 5] * std::polar<float>(1.0, -6 * M_PI * i / N));
            temp += (p4[i % (int)N / 5] * std::polar<float>(1.0, -8 * M_PI * i / N));
            x[i] = temp;
        }
    }
}

main.cpp #

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <iostream>
#include "fft.hpp"


typedef vector<complex<float>> complexSignal;

int main()
{
    complexSignal abit;
    int N = 9;
    abit.push_back({1,0});
    abit.push_back({1.00071 ,0 });
    abit.push_back({1.00135 ,0 });
    abit.push_back({1.00193 ,0 });
    abit.push_back({1.00245 ,0 });
    abit.push_back({1.0029 ,0 });
    abit.push_back({1.00329 ,0 });
    abit.push_back({1.00361 ,0 });
    abit.push_back({1.00387 ,0 });
    std::cout << "Before:" << std::endl;
    for (int i = 0; i < N; i++)
    {
        std::cout << abit[i] << std::endl;
    }
    std::cout << "After:" << std::endl;
    fft(abit);
    for (int i = 0; i < N; i++)
    {
        std::cout << abit[i] << std::endl;
    }
    return 0;
}

结果对比 #

实际输出 #

(9.02011,0)
(5.83089,-4.89513)
(0.700632,-3.98993)
(-0.000289979,0.000502368)
(-0.00218513,0.000362784)
(-0.00179241,0.00139188)
(-0.000289979,-0.000502368)
(0.000175771,-0.00354373)
(-0.003268,-0.00558837)

预期输出 #

(9.020109999999999+0j)
(-0.0032675770104925446+0.005588577982060319j)
(-0.0023772289746976797+0.0024179090499282354j)
(-0.0022250000000012538+0.0011691342951078987j)
(-0.002185194014811494+0.00036271471530890747j)
(-0.0021851940148113033-0.00036271471530980844j)
(-0.0022249999999994774-0.0011691342951105632j)
(-0.002377228974696629-0.0024179090499291786j)
(-0.00326757701049002-0.005588577982061138j)

问题原因 #

Bug出在模运算的括号优先级错误：C++中%和/运算符优先级相同，遵循左结合规则，所以i % (int)N / 3实际执行顺序是(i % (int)N) / 3，而不是预期的i % (N/3)。比如N=9时，i=1的计算结果是1%9/3=0，i=3的计算结果是3%9/3=1，导致从p0/p1/p2中取数的索引完全错误。
基5逻辑存在完全相同的优先级错误。

修复方案 #

修改fft.hpp中基3、基5部分的模运算括号即可：

• 基3部分：将所有i % (int)N / 3替换为i % (int)(N / 3)
• 基5部分：将所有i % (int)N / 5替换为i % (int)(N / 5)
  修改后的基3逻辑示例：

for (int i = 0; i < N; i++)
{
    int k = i % (int)(N / 3);
    complex<float> temp = p0[k];
    temp += (p1[k] * std::polar<float>(1.0, -2 * M_PI * i / N));
    temp += (p2[k] * std::polar<float>(1.0, -4 * M_PI * i / N));
    x[i] = temp;
}

修改后重新编译运行，即可得到和numpy.fft一致的正确结果。

内容的提问来源于stack exchange，提问作者yarin Cohen