在 C 中使用 FFTW 进行傅里叶回归

Question

我尝试使用 FFTW 库进行傅里叶回归

#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;
}

但回归与输入数据不匹配。我哪里做错了？