fft_tool.cpp

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

#include "FFT_Tool.h"
#include "iostream.h"
#include <math.h>
int converse(int t,int length) {
	  int s=0;
        for(int i=0;i<32;i++){
           if((int)pow(2, i)>=length){
              s=i;
              break;
           }
        }
		
       int * a=new int[s];
	  for(int j=0;j<s;j++)
		*(a+j)=0;
        
        int n=0;
        while(t!=0){
          *(a+n)=t%2;
          n++;
          t=t/2;
        }
        int sum=0;
        n=0;
        while(s>0){
           sum=sum+(int)(*(a+s-1)*pow(2, n));
           n++;
           s--;
        }
        return sum;
}



Complex * FFT_T_1(Complex xin[],int length){
	  //时间抽取1维fft
        
		Complex *ptr=new Complex[length];
		int i, j, k=0;
        Complex temp; 
        
		bool *flag=new bool[length] ;
        for (i = 0; i < length; i++) {
            *(flag+i) = false;//标志 

        }
        //求逆序 
        for (i = 0; i < length; i++) {

            j = i;
            k = converse(i,length);
          
            if ((j != k) && (*(flag+j) == false)) {
                
                temp = xin[k];
                xin[k] = xin[j];
                xin[j] = temp;
                *(flag+k) = true;
            }
        }

        for (i = 0; i < length; i = i + 2) {
            //两点DFT
            *(ptr+i) = complexAdd(xin[i], xin[i + 1]);
            *(ptr+i + 1) = complexDec(xin[i], xin[i + 1]);
        }

        for (i = 0; i < length; i++) {
            //xin存xout副本
            xin[i] = *(ptr+i);
        }

        int m = 2;
        int t;
        while (m != length) {
            //蝶形算法
            for (j = 0; j < length; j = j + 2 * m) {
                for (i = 0; i < m; i++) {
                    //k = i; 
					Complex c;
					c.real=cos(i * 2 * PI / (2 * m));
					c.imag=sin(-i * 2 * PI / (2 * m));
                    *(ptr+i + j) = complexAdd(xin[i + j], complexEE(c, xin[i + j + m]));
                    *(ptr+i + j + m) = complexDec(xin[i + j], complexEE(c, xin[i + j + m]));

                }
            }
            for (t = 0; t < length; t++) {
                xin[t] = *(ptr+t);//暂存 

            }
            m = m * 2;
        }



return ptr;

}


Complex * IDFT(Complex xk[],int length){
	//fft反变换
		Complex *ptr=new Complex[length];

     	double N=(double)1/length;
       
        for(int i=0;i<length;i++)
            *(ptr+i)=complexG(xk[i]);
        ptr=FFT_T_1(ptr,length);
        for(int j=0;j<length;j++){
			xk[j]=complexRC(N, complexG(*(ptr+j)));
        }
       ptr=xk;

        return ptr;

}