《C语言数值Recipes》four1函数FFT实现运行时错误求助

核心问题是原书代码采用1-based数组索引，而C语言默认是0-based索引，直接复制代码会导致数组访问越界。

具体错误分析 #

• 原书假设data数组从索引1开始，长度为2*nn+1，但你的测试数组是0-based，长度为2*nn（比如测试用的8个元素对应nn=4）。
• 代码中for (i = m; i <= n; i += istep)的i <= n会让j = i + mmax达到n + mmax，远超数组最大索引2*nn-1，直接触发越界访问。
• 第一个循环中j的初始化和更新逻辑都是基于1-base设计，放到0-base数组里会导致索引错位，甚至访问到数组末尾之后的内存。

修正方案 #

将所有数组访问逻辑适配C语言的0-base规则，调整关键循环条件：

1. 初始化j = 0，而非j = 1
2. 第一个循环i从0开始，步长2，终止条件i < n（n=2*nn）
3. 内层循环i的起始值改为m-1，终止条件改为i < n
4. 调整位反转置换中j的判断条件为j >= m

修正后的完整代码 #

#include <iostream>
#include <stdio.h>
using namespace std;

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void four1(double* data, unsigned long nn, int isign)
{
    unsigned long n, mmax, m, j, istep, i;
    double wtemp, wr, wpr, wpi, wi, theta;
    double tempr, tempi;

    n = nn << 1;
    j = 0; // 适配0-base起始索引
    for (i = 0; i < n; i += 2) { // i从0开始遍历
        if (j > i) {
            SWAP(data[j], data[i]);
            SWAP(data[j + 1], data[i + 1]);
        }
        m = n >> 1;
        while (m >= 2 && j >= m) { // 适配0-base的判断逻辑
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax = 2;
    while (n > mmax) {
        istep = mmax << 1;
        theta = isign * (6.28318530717959 / mmax);
        wtemp = sin(0.5 * theta);
        wpr = -2.0 * wtemp * wtemp;
        wpi = sin(theta);
        wr = 1.0;
        wi = 0.0;
        for (m = 1; m < mmax; m += 2) {
            for (i = m-1; i < n; i += istep) { // 起始索引改为m-1，终止条件i < n
                j = i + mmax;

                tempr = wr * data[j] - wi * data[j + 1];
                tempi = wr * data[j + 1] + wi * data[j];
                data[j] = data[i] - tempr;
                data[j + 1] = data[i + 1] - tempi;
                data[i] += tempr;
                data[i + 1] += tempi;
            }
            wr = (wtemp = wr) * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}
#undef SWAP


int main()
{
    // 测试用例
    double data[] = {1, 1, 2, 0, 1, 3, 4, 0};

    four1(data, 4, 1);

    // 遍历所有8个元素
    for (int i = 0; i < 8; i++) {
        cout << data[i] << " ";
    }
    cout << endl;
    return 0;
}

补充说明 #

修正后代码解决了索引越界问题，测试用例可正常运行。注意原书four1函数要求nn必须是2的整数次幂，若使用非2次幂的nn需要额外处理。

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