为何我的FFT输出存在额外谐波？

API 探路者

2026-6-23

FFT实信号频域转换结果异常问题排查

我用C++实现了基于欧拉公式分离实虚部的FFT算法，将实信号从时域转频域，但输出结果异常——除正确谐波外出现额外谐波。

测试输入：16元素数组 {9, -2.82, -5, 2.82, 1, 2.82, -5, -2.82, 9, -2.82, -5, 2.82, 1, 2.82, -5, -2.82}
实际输出幅值序列：{0,0,8,0,40,0,24,0,0,0,24,0,40,0,8,0}
预期幅值序列：{0,0,0,40,0,32,0,0,0,32,0,40,0,0,0}

以下是我的实现代码：

// for example dn = {9, -2.82, -5, 2.82, 1, 2.82, -5, -2.82, 9, -2.82, -5, 2.82, 1, 2.82, -5, -2.82}, kn = 16 is the length of array
printf("\n-----------------------------------------------------------Shifted Data--------------------------------------------------------------------------------------------------------\n");
float* dm = (float*)malloc(kn * sizeof(float));

uint32_t rev = 0;
uint32_t n = 0;
int order = 0;
for (int i = 0; i < kn; i++) {
    n = i;
    if (!(n & (n - 1)))
        order++;
    for (int j = 0; j < log2(kn); j++) {  // shift input data 
        rev = rev << 1;
        rev = rev | (n & 1);
        n = n >> 1;
    }
//  printf("rev = %d, i = %d order = %d\n", rev, i, order);
//  std::cout << "rev = " << std::bitset<9>(rev) << " i = " << std::bitset<9>(i) <<std::endl;
    dm[i] = dn[rev];
    printf("dm[%d] = %f d[%d] = %f, %d\n\n", i, dm[i], i, dn[i], rev);
    rev = 0;
}

int numselect = kn;
int selectcounter = 0;

while (numselect != 1) {
    numselect = numselect / 2;
    selectcounter++;
}

selectcounter--;
int tempval1 = 0;
int tempval2 = 0;
int tempvalcounter = 0;
int counter = 0;
double **ftnreal = new double* [log2(kn)]; // m*N array for real output
for (int i = 0; i < log2(kn); i++)
    ftnreal[i] = new double[kn];

selectcounter = 0;
for (int i = 0; i < log2(kn); i++) {
    for (int j = 0; j < kn; j++)
        ftnreal[i][j] = 0.0;
    if (i == 0) {
        for (int j = 0; j < kn / 2; j++) {
            ftnreal[i][j * 2] = dm[2 * j] + cos(2 * M_PI * i / 2) * dm[2 * j + 1]; //calculating data for the first column of elements
            ftnreal[i][j * 2 + 1] = dm[2 * j] - cos(2 * M_PI * i / 2) * dm[2 * j + 1];
            printf("ftnreal[%d][%d] = dm[%d] + cos(2 * %f*%d/2)*dm[%d]\n", i, j*2,j*2,M_PI,i,2*j + 1);
            printf("ftnreal[%d][%d] = dm[%d] - cos(2 * %f*%d/2)*dm[%d]\n", i, j * 2 + 1, j * 2, M_PI, i, 2 * j + 1);
        }
    }
    else {
        tempval1 = pow(2, i + 1);
        tempval2 = pow(2, i);
        counter = 0;
        tempvalcounter = 0;
        printf("i = %d, tempval1 = %d, tempval2 = %d\n", i, tempval1, tempval2);
        for (int j = 0; j < kn / 2; j++) {
                ftnreal[i][counter] = ftnreal[i - 1][counter] + cos(2 * M_PI * (float(tempvalcounter) / float(tempval1))) * ftnreal[i - 1][counter + tempval2]; //calculating data for the rest of columns
                ftnreal[i][counter + tempval2] = ftnreal[i - 1][counter] - cos(2 * M_PI * (float(tempvalcounter) / float(tempval1))) * ftnreal[i - 1][counter + tempval2];
                printf("ftnreal[%d][%d] = ftnreal[%d][%d] + cos(2*%f*(%f/%f))*ftnreal[%d][%d]       / %f      %f\n", i,counter,i - 1,counter,M_PI,float(tempvalcounter), float(tempval1),i - 1,counter + tempval2, float(tempvalcounter) / float(tempval1), cos(2 * M_PI * (float(tempvalcounter) / float(tempval1))));
                printf("ftnreal[%d][%d] = ftnreal[%d][%d] - cos(2*%f*(%f/%f))*ftnreal[%d][%d]       / %f      %f\n", i, counter + tempval2, i - 1, counter, M_PI, float(tempvalcounter), float(tempval1), i - 1, counter + tempval2, float(tempvalcounter) / float(tempval1), cos(2 * M_PI * (float(tempvalcounter) / float(tempval1))));
                counter++;
                tempvalcounter++;
                if (tempvalcounter % tempval2 == 0) {
                    counter = counter + tempval2;
                    tempvalcounter = 0;
                }
        }
    }
    cout << endl;
}


double** ftnimag = new double* [log2(kn)]; // same goes for imaginary part
for (int i = 0; i < log2(kn); i++)
    ftnimag[i] = new double[kn];

for (int i = 0; i < log2(kn); i++) {
    for (int j = 0; j < kn; j++)
        ftnimag[i][j] = 0.0;
    if (i == 0) {
        for (int j = 0; j < kn/2; j++) {
            ftnimag[i][j * 2] = -sin(2 * M_PI * i / 2) * dm[2 * j + 1];
            ftnimag[i][j * 2 + 1] = sin(2 * M_PI * i / 2) * dm[2 * j + 1];
        }
    }
    else {
        tempval1 = pow(2, i + 1);
        tempval2 = pow(2, i);
        counter = 0;
        int tempvalcounter = 0;
        printf("i = %d, tempval1 = %d, tempval2 = %d\n", i, tempval1, tempval2);
        for (int j = 0; j < kn / 2; j++) {
            ftnimag[i][counter] = ftnimag[i - 1][counter] - sin(2 * M_PI * (float(tempvalcounter) / float(tempval1))) * ftnimag[i - 1][counter + tempval2];
            ftnimag[i][counter + tempval2] = ftnimag[i - 1][counter] + sin(2 * M_PI * (float(tempvalcounter) / float(tempval1))) * ftnimag[i - 1][counter + tempval2];
            printf("ftnreal[%d][%d] = ftnimag[%d][%d] - sin(2*%f*(%f/%f))*ftnreal[%d][%d]       / %f     %f\n", i, counter, i - 1, counter, M_PI, float(tempvalcounter), float(tempval1), i - 1, counter + tempval2, float(tempvalcounter) / float(tempval1), -sin(2*M_PI* (float(tempvalcounter) / float(tempval1))));
            printf("ftnreal[%d][%d] = ftnimag[%d][%d] + sin(2*%f*(%f/%f))*ftnreal[%d][%d]       / %f    %f\n", i, counter + tempval2, i - 1, counter,M_PI, float(tempvalcounter), float(tempval1), i - 1, counter + tempval2, float(tempvalcounter) / float(tempval1), sin(2*M_PI* (float(tempvalcounter) / float(tempval1))));
            counter++;
            tempvalcounter++;
            if (tempvalcounter % tempval2 == 0) {
                counter = counter + tempval2;
                tempvalcounter = 0;
            }
        }
    }
    cout << endl;
}
printf("-----------------------------------------------------------Real--------------------------------------------------------------------------------------------------------\n");

//printf("%d", selectcounter);
for (int i = 0; i < log2(kn); i++) {
    for (int j = 0; j < kn; j++) {
        cout << ftnreal[i][j] << " ";
    }
    cout << endl;
    cout << endl;
    cout << endl;
}

printf("-----------------------------------------------------------Imaginary--------------------------------------------------------------------------------------------------------\n");

for (int i = 0; i < log2(kn); i++) {
    for (int j = 0; j < kn; j++) {
        cout << ftnimag[i][j] << " ";
    }
    cout << endl;
    cout << endl;
    cout << endl;
}

printf("-----------------------------------------------------------Abs--------------------------------------------------------------------------------------------------------\n");

tempval1 = log2(kn) - 1;
printf("%d\n", tempval1);
float* ftnabs = (float*)malloc(kn * sizeof(float));
for (int i = 0; i < kn; i++) {
    ftnabs[i] = sqrt(ftnreal[tempval1][i] * ftnreal[tempval1][i] + ftnimag[tempval1][i] * ftnimag[tempval1][i]);
    cout << ftnabs[i] << " ";
    cout << 2 * ftnabs[i] / kn << " ";
}

核心问题点

位逆序置换逻辑错误
代码中order变量完全多余，且逆序索引计算依赖浮点数log2(kn)存在精度风险，位反转逻辑未正确实现。正确的位逆序应针对kn的二进制位数（如16对应4位），逐位反转输入索引的二进制位。
FFT蝶形运算公式错误
实部与虚部的蝶形运算被错误分离，未遵循标准FFT复蝶形公式：
```
X[k] = X0[k] + W_N^k * X1[k]
X[k+N/2] = X0[k] - W_N^k * X1[k]
```
其中旋转因子W_N^k = cos(2πk/N) - i·sin(2πk/N)，代码未结合实虚部联动计算，导致结果偏离预期。
旋转因子角度计算错误
计算旋转因子时使用pow(2, i+1)作为分母，而非当前阶段的块大小或总长度kn，导致旋转因子角度错误，引入额外谐波。
内存管理不规范
使用malloc和new分配的内存未释放，存在内存泄漏；log2(kn)返回浮点数，直接作为数组下标易引发精度问题，应改用整数计算二进制位数。

修正建议

修复位逆序：先计算kn的二进制位数m（如int m = 0; int temp = kn; while(temp>1) { temp>>=1; m++; }），然后对每个索引i，反转其m位二进制得到逆序索引。
合并实虚部蝶形运算：按照标准FFT公式，同时更新实部和虚部，避免分离计算的逻辑错误。
修正旋转因子：根据当前蝶形阶段的块大小计算正确的旋转因子角度，确保W_N^k的参数准确。
完善内存释放：在程序结束前，用free(dm)、free(ftnabs)释放堆内存，用delete[]释放ftnreal和ftnimag的二级数组。

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

为何我的FFT输出存在额外谐波？

核心问题点 #

修正建议 #

核心问题点

修正建议