快速 卷积 （用 fft 快速傅里叶变换 做）c 源码

#include "FFT_Tool.h" #include "Converse.h" #include <math.h> double * conv_FFT_2(double xin1[],double xin2[],int length1,int length2){ //卷积定理 int L = length1; int M = length2; int length = 0; for (int i = 0; i < 100; i++) { if (((int) pow(2, i)) >= (L + M - 1)) { length = (int) pow(2, i); break; } } double * xout = new double[length]; Complex * temp1=new Complex[length]; Complex * temp2=new Complex[length]; Complex * temp3 = new Complex[length]; for (int k = 0; k < length; k++) { for (int i = 0; i < L; i++,k++) { Complex c; c.real=xin1[i]; c.imag=0; *(temp1+i) =c; } for (i = 0; i < length - L; i++,k++) { Complex c; c.real=0; c.imag=0; *(temp1+i + L) = c; } } for ( i = 0; i < M; i++) { Complex c; c.real=xin2[i]; c.imag=0; *(temp2+i) = c; } for ( i = 0; i < length - M; i++) { Complex c; c.real=0; c.imag=0; *(temp2+i + M) = c; } Complex *temp4 = FFT_T_1(temp1,length); Complex *temp5 = FFT_T_1(temp2,length); for ( i = 0; i < length; i++) { *(temp3+i) = complexEE(temp4[i], temp5[i]); } temp1 = IDFT(temp3,length); for ( i = 0; i < L + M - 1; i++) { *(xout+i) = (temp1+i)->real; } return xout; } double * conv_FFT_1(double xin1[],double xin2[],int length1,int length2){ int L = length1; int M = length2; int length = 0; int i,k; for (i = 0; i < 100; i++) { if (((int) pow(2, i)) >= (L + M - 1)) { length = (int) pow(2, i); break; } } double * xout = new double[length]; double * temp1=new double[length]; double * temp2=new double[length]; Complex * temp3=new Complex[length]; for(i=0;i<L;i++){ *(temp1+i)=xin1[i]; } for(i=L;i<length;i++){ *(temp1+i)=0; } for(i=0;i<M;i++){ *(temp2+i)=xin2[i]; } for(i=M;i<length;i++){ *(temp2+i)=0; } for(i=0;i<length;i++){ //temp3[i]=new Complex(temp1[i],temp2[i]); temp3[i].real=temp1[i]; temp3[i].imag=temp2[i]; } Complex *xk=FFT_T_1(temp3,length); Complex *xk1=new Complex[length]; Complex *xk2=new Complex[length]; xk1->real=complexAdd(*(xk+0),complexG(*(xk+0))).real/2; xk1->imag=complexAdd(*(xk+0),complexG(*(xk+0))).imag/2; xk2->real=complexDec(*(xk+0),complexG(*(xk+0))).imag/2; xk2->imag=-complexDec(*(xk+0),complexG(*(xk+0))).real/2; for( k=1;k<length;k++){ (xk1+k)->real=complexAdd(*(xk+k),complexG(*(xk+length-k))).real/2; (xk1+k)->imag=complexAdd(*(xk+k),complexG(*(xk+length-k))).imag/2; (xk2+k)->real=complexDec(*(xk+k),complexG(*(xk+length-k))).imag/2; (xk2+k)->imag=-complexDec(*(xk+k),complexG(*(xk+length-k))).real/2; } Complex *xk3=new Complex[length]; for(i=0;i<length;i++){ *(xk3+i)=complexEE(xk1[i], xk2[i]); } Complex *temp4=IDFT(xk3,length); for(i=0;i<length;i++){ *(xout+i)=(temp4+i)->real; } // return xout; } double * conv_Partion(double xin1[], double xin2[],int length1,int length2){ // xin2为长序列 ,重叠相加法 int length = length1 + length2 ; double * xout = new double[length]; int M = length1; int L = M; //初始值 double * temp1 = new double[L]; for (int j = 0; j < L; j++) { *(temp1+j) = xin2[j]; } double *t1 = new double[L + M - 1]; t1 = conv_FFT_1(xin1, temp1,M,L); for (int k = 0; k < L + M - 1; k++) { *(xout+k) = *(t1+k); } for (int l = L; l < length2; l = l + L) { if (length2 - l >= L) { double *temp = new double[L]; for (int j = 0; j < L; j++) { *(temp+j) = xin2[j + l]; } double *t = new double[L + M - 1]; t = conv_FFT_1(xin1, temp,M,L); //重叠相加 for (int k = 0; k <= M - 1; k++) { *(xout+l + k) = *(t+k) + *(xout+l + k); } for ( k = 0; k < L; k++) { *(xout+l + M - 1 + k) = *(t+M - 1 + k); } } else { double *temp = new double[length2 - l]; for (int j = 0; j < length2 - l; j++) { *(temp+j) = xin2[j + l]; } double *t = new double[length2 - l + M - 1]; t = conv_FFT_1(xin1, temp,M,length2 - l); //重叠相加 for (int k = 0; k < M - 1; k++) { *(xout+l + k) = *(t+k) + *(xout+l + k); } for ( k = 0; k < length2 - l; k++) { //int m = 1 + M + k; *(xout+l + M - 1 + k) = *(t+M - 1 + k); } } } return xout; }