使用FFTW库实现傅里叶回归结果不匹配，请求排查问题

问题描述 #

我尝试用FFTW库实现傅里叶回归，但生成的回归曲线与输入数据不匹配，以下是我的C语言代码和部分输出结果，请问哪里出错了？

原始代码 #

#include <stdio.h>
#include <math.h>
#include <fftw3.h>

#define PI 3.14159265358979323846

int main() {
    int n = 100;
    double y[n];
    double dt = 0.01;
    double t[n];
    double a0, a[n/2], b[n/2];
    int i;
// Generate data
for (i = 0; i < n; i++) {
    t[i] = i * dt;
    y[i] = sin(2 * PI * 0.5 * t[i]) + sin(2 * PI * 2 * t[i]);
}

// Create FFTW plan
fftw_complex* in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
fftw_complex* out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);
fftw_plan p = fftw_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

// Load data into input array and perform FFT
for (i = 0; i < n; i++) {
    in[i][0] = y[i];
    in[i][1] = 0.0;
}
fftw_execute(p);

// Extract Fourier coefficients
a0 = out[0][0] / n;
for (i = 1; i < n/2; i++) {
    a[i] = 2.0 * out[i][0] / n;
    b[i] = -2.0 * out[i][1] / n;
}

// Generate regression line
double y_reg[n];
for (i = 0; i < n; i++) {
    y_reg[i] = a0;
    for (int j = 1; j < n/2; j++) {
        y_reg[i] += a[j] * cos(2 * PI * j * t[i] / n) + b[j] * sin(2 * PI * j * t[i] / n);
    }
}

// Print data and regression line
printf("t\ty\ty_reg\n");
for (i = 0; i < n; i++) {
    printf("%f\t%f\t%f\n", t[i], y[i], y_reg[i]);
}

// Free memory and destroy plan
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);

return 0;
}

典型输出 #

t       y       y_reg
0.300000    0.221232    0.041583
0.310000    0.142533    0.043084
0.320000    0.073815    0.044592
0.330000    0.016414    0.046106
0.340000    -0.028520   0.047628
0.350000    -0.060050   0.049156
0.360000    -0.077460   0.050690
0.370000    -0.080272   0.052231
0.380000    -0.068250   0.053778
0.390000    -0.041406   0.055331
0.400000    -0.000000   0.056890
0.410000    0.055467    0.058456
0.420000    0.124255    0.060027
0.430000    0.205404    0.061604
0.440000    0.297740    0.063186
0.450000    0.399903    0.064774

错误分析与修正 #

1. 核心错误：重构时的频率公式错误 #

生成回归曲线的循环中，角频率计算误用了采样点数n作为分母，正确的分母应该是总采样时长T（T = n * dt）。

原始错误代码片段：

y_reg[i] += a[j] * cos(2 * PI * j * t[i] / n) + b[j] * sin(2 * PI * j * t[i] / n);

修正后的代码片段：

double T = n * dt; // 计算总采样时长
y_reg[i] += a[j] * cos(2 * PI * j * t[i] / T) + b[j] * sin(2 * PI * j * t[i] / T);

由于本例中n*dt=100*0.01=1，也可以简化为：

y_reg[i] += a[j] * cos(2 * PI * j * t[i]) + b[j] * sin(2 * PI * j * t[i]);

这个错误导致重构的信号频率被缩小了100倍，最终回归曲线几乎是一条缓慢变化的直线，完全偏离原始信号。

2. 次要问题：频谱泄漏导致的不完全匹配 #

原始信号包含0.5Hz的分量，而当前采样参数（总时长1s）的基频为1Hz，0.5Hz并非基频的整数倍，FFT无法精准捕捉该分量，会将其能量分散到多个频率bin中。即使修正频率公式，回归曲线也无法完全重合于原始信号，但2Hz的分量会被正确还原，整体趋势与原始信号一致。

如果希望完全匹配原始信号，可以调整采样参数：将总时长改为2s（比如n=200，dt=0.01），此时基频为0.5Hz，0.5Hz和2Hz都是基频的整数倍（对应j=1和j=4），FFT就能精准捕捉两个分量，回归曲线会与原始信号完全重合。

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