用c语言写重叠相加法

重叠相加法（Overlap-add method）是一种数字信号处理中常用的算法，用于实现长信号的快速卷积。下面是用C语言实现重叠相加法的示例代码：

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

#define PI 3.14159265358979323846

// FFT实现
void fft(double *xr, double *xi, int n, int sign) {
    int i, j, k, m, nv2;
    double tr, ti, ur, ui, sr, si, tvr;

    m = (int)log2(n);
    nv2 = n / 2;

    j = nv2;
    for (i = 1; i < n - 1; i++) {
        if (i < j) {
            tr = xr[j];
            ti = xi[j];
            xr[j] = xr[i];
            xi[j] = xi[i];
            xr[i] = tr;
            xi[i] = ti;
        }
        k = nv2;
        while (k <= j) {
            j -= k;
            k /= 2;
        }
        j += k;
    }

    for (i = 1; i <= m; i++) {
        k = 1 << i;
        for (j = 0; j < nv2; j += k) {
            for (sr = 0.0, si = 0.0, m = j; m < j + k / 2; m++) {
                ur = cos(PI * (m - j) / (double)(k / 2));
                ui = sign * sin(PI * (m - j) / (double)(k / 2));
                tr = xr[m + k / 2] * ur - xi[m + k / 2] * ui;
                ti = xr[m + k / 2] * ui + xi[m + k / 2] * ur;
                xr[m + k / 2] = xr[m] - tr;
                xi[m + k / 2] = xi[m] - ti;
                xr[m] += tr;
                xi[m] += ti;
            }
        }
    }

    if (sign == -1) {
        for (i = 0; i < n; i++) {
            xr[i] /= n;
            xi[i] /= n;
        }
    }
}

// 重叠相加法实现
void overlap_add(double *x, int lx, double *h, int lh, double *y, int ly) {
    int n = 1;
    while (n < lx + lh - 1) {
        n *= 2;
    }

    double *xr = (double*)malloc(n * sizeof(double));
    double *xi = (double*)malloc(n * sizeof(double));
    double *hr = (double*)malloc(n * sizeof(double));
    double *hi = (double*)malloc(n * sizeof(double));
    double *yr = (double*)malloc(n * sizeof(double));
    double *yi = (double*)malloc(n * sizeof(double));
    int i;

    for (i = 0; i < n; i++) {
        xr[i] = 0.0;
        xi[i] = 0.0;
        hr[i] = 0.0;
        hi[i] = 0.0;
        yr[i] = 0.0;
        yi[i] = 0.0;
    }

    for (i = 0; i < lx; i++) {
        xr[i] = x[i];
    }

    for (i = 0; i < lh; i++) {
        hr[i] = h[i];
    }

    fft(xr, xi, n, 1);
    fft(hr, hi, n, 1);

    for (i = 0; i < n; i++) {
        yr[i] = xr[i] * hr[i] - xi[i] * hi[i];
        yi[i] = xr[i] * hi[i] + xi[i] * hr[i];
    }

    fft(yr, yi, n, -1);

    for (i = 0; i < ly; i++) {
        y[i] = yr[i];
    }

    free(xr);
    free(xi);
    free(hr);
    free(hi);
    free(yr);
    free(yi);
}

int main() {
    double x[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double h[] = {1.0, 2.0, 3.0};
    double y[7];
    overlap_add(x, 5, h, 3, y, 7);
    int i;
    for (i = 0; i < 7; i++) {
        printf("%f ", y[i]);
    }
    printf("\n");
    return 0;
}

这段代码实现了重叠相加法，其中用到了FFT算法。输入信号x的长度为lx，卷积核h的长度为lh，输出信号y的长度为ly（ly = lx + lh - 1）。在main函数中，我们给出了一个示例，其中x为{1.0, 2.0, 3.0, 4.0, 5.0}，h为{1.0, 2.0, 3.0}，y的长度为7。输出结果为{1.0, 4.0, 10.0, 16.0, 22.0, 16.0, 9.0}。