利用FFT进行快速卷积c程序例子 [复制链接]

/*当序列长度>=128时，这种方法就比普通线性卷积速度快*/
#include<stdio.h>
#include<stdlib.h>
#include"FFT.h"
const int L1=4;
const int L2=5;
//普通线性卷积
void conv(float x[],int m,float y[],int n,float z[],int l)
{
    int i,j;
    for(i=0; i<l; i++)
    {
           z=0.0;
           for(j=0; j<n; j++)
                  if(i-j>=0&&i-j<=m-1)
        z+=y[j]*x[i-j];
    }
}
void main()
{
    int i, L;
    L=L1+L2-1;//保证L为2的幂次
    float x[L1]={1,1,1,1};
    float y[L2]={1,1,1,1,1};
    float *z1, *z2;
    z1=(float *)calloc(L,sizeof(float));
    z2=(float *)calloc(L,sizeof(float));

    //直接线性卷积
    conv(x,L1,y,L2,z1,L);

    //利用FFT
    //1.填充序列长度到L
    float *x1, *y1;
    float *xi, *yi, *zi;
    x1=(float *)calloc(L,sizeof(float));
    y1=(float *)calloc(L,sizeof(float));
                    xi=(float *)calloc(L,sizeof(float));
    yi=(float *)calloc(L,sizeof(float));
    zi=(float *)calloc(L,sizeof(float));
    for(i=0; i<L1; i++)
          x1=x;
    for(i=0; i<L2; i++)
          y1=y;
    //2.进行傅里叶正变换
    fft(x1,xi,3,1);
    fft(y1,yi,3,1);
    //3.频谱相乘
    for(i=0; i<L; i++)
    {
        z2=x1*y1-xi*yi;
        zi=x1*yi+xi*y1;
    }
    //4.傅里叶反变换
    fft(z2,zi,3,-1);

    //输出结果
    printf("普通线性卷积:        快速卷积:\n");
    for(i=0; i<L; i++)
          printf("z1[%d]=%f        z2[%d]=%f\n",i,z1,i,z2);

    free(z1);free(z2);free(x1);free(y1);
    free(xi);free(yi);free(zi);
}