converse.cpp

用快速傅立叶变换做的快速卷积,计算速度快.

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