本文共 3326 字，大约阅读时间需要 11 分钟。

 

1)

蝶形变换：普通的FFT算法称为基2的FFT算法，这种算法的核心是蝶形变换  长度为n=2^k1的变换共需要做   k1   *   n/2   次蝶形变换，(如上图所示)若需变换数据表示为一个复数数组c[],则每次蝶形变换有2个输入   c[i],c[i+s],两个输出:c[i],c[i+s]，s成为翅间距。   每个变换的基本算法是：

        t=wr   *   c[i+s];       c[i+s]=c[i]-t;     c[i]=c[i]+t;       前面说过，长度为n=2^k1的变换共需要做   k1   *   n/2次变换，实际的程序是一个3层循环,共需要k1*k2*(k3/2)次变换(k2*k3/2=n/2)。前面的wr是w的整数次方，w=e^(-   2*PI/k3)   (k3=2,4,8,16...n,PI是圆周率)，也成为旋转因子，例如n=32的变换需要log2(32)=5趟变换:       第1趟变换需要16*1次变换，翅间距是1,     若w=e^(-2*PI/2)，   则wr=w^1       第2趟变换需要8*2次变换，   翅间距是2,     若w=e^(-2*PI/4)，   则wr=w^1,w^2       第3趟变换需要4*2次变换，   翅间距是4,     若w=e^(-2*PI/8)，   则wr=w^1,w^2,w^3,w^4         第4趟变换需要2*8次变换，   翅间距是8,     若w=e^(-2*PI/16)，则wr=w^1,w^2,w^3,w^4,w^5,w^6,w^7,w^8       第5趟变换需要16*1次变换，翅间距是16,   若w=e^(-2*PI/32)，则wr=w^1,w^2,w^3,w^4,w^5...w^15,w^16

void  fft(struct complex c_IO[],int m){	int L, B, j, k,  p, q;	for(k = 0; k < SAMPLECOUNT; k++)	{	  g_out[k] = c_IO[k];	}	for(L = 1.0; L <= m; L++)	{		B=(int)pow(2.0,	L-1);		for(j=0;j<=B-1;j++)		{			p=j*(int)pow(2.0, m-L);			q=(int)pow(2.0, L);			for(k=j;k<=SAMPLECOUNT-1;k=k+q)			{				temp[k]=complex_add(g_out[k],complex_mult(W_n_k(SAMPLECOUNT,p) ,g_out[k+B] ) );				temp[k+B]=complex_remove(g_out[k],complex_mult(W_n_k(SAMPLECOUNT,p) ,g_out[k+B] ) );				g_out[k]=temp[k];				g_out[k+B]=temp[k+B];			}		}	}}

 

//*************************************************************************// *// * 函数名称：// *   FFT()// *// * 参数:// *   complex
   
     * TD- 指向时域数组的指针// *   complex
    
      * FD- 指向频域数组的指针// *   r－2的幂数，即迭代次数// *// * 返回值:// *   无。// *// * 说明:// *   该函数用来实现快速付立叶变换。// *// ************************************************************************/void FFT(	complex
     
       *TD,  complex
      
        *FD, 			complex
       
         *X1,  complex
        
          *X2, complex
         
           *Omega, int r){ // 付立叶变换点数 long count; // 循环变量 int i,j,k; // 中间变量 int bfsize,p; // 角度 double angle; complex
          
            *X; // 计算付立叶变换点数 count = 1 << r; // 分配运算所需存储器 //Omega = new complex
           
            [count / 2]; //X1 = new complex
            
             [count]; //X2 = new complex
             
              [count]; // 计算加权系数 for(i = 0; i < count / 2; i++) { angle = -i * 3.1415926 * 2 / count; Omega[i] = complex
              
               (cos(angle), sin(angle)); } // 将时域点写入X1 memcpy(X1, TD, sizeof(complex
               
                ) * count); // 采用蝶形算法进行快速付立叶变换 for(k = 0; k < r; k++) { for(j = 0; j < 1 << k; j++) { bfsize = 1 << (r-k); for(i = 0; i < bfsize / 2; i++) { p = j * bfsize; X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2]; X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * Omega[i * (1<
               
              
             
            
           
          
         
        
       
      
     
    
   

2)复数数组排序，在基2的蝶形变换中，复数数组需要重新排序，c[i]   要放置到数组c的第   reverse(c[i])   的位置，m=reverse(n)   函数的算法是这样的:

若   n的   k位2进制的为b[],   b[k-1],B[k-2],...b[2],b[1],b[0],(   b[i]   等于1或者0,b[0]为最低bit).   则m=reverse(n)的2进制的为   b[0],b[1],b[2],b[3],...b[k-1]   (b[k-1]为最低bit).

2.更复杂的变换算法：基2的蝶形变换算法不止一种，它可分为2类，一类为基2时分傅立叶变换，另一类为基2频分傅立叶变换。上例的变为基2时分算法，在每一趟变换中，翅间距依次变大，第一趟为2，最后一趟为n/2,数组重排在变换之前进行，基2频分算法正好相反，翅间距依次缩小，起于n/2，止于2，数组重排在蝶形变换之后进行。   在＜傅立叶变换＞一书中，提到3种基2时分变换，3种基2频分变换。上述算法称为基2时分FFT第二种算法。我在看你的这个程序的同时，还看到朱志刚写的一个FFT程序，这个程序的算法是基2时分FFT第一种算法，它比经典的算法更复杂，需要对wr数组进行逆序排列。 3.更复杂的FFT算法，除了基2   的FFT算法外，还有更加复杂的基4算法，基8算法，甚至基3,基5算法，纯粹的基4算法只能计算长度为4^k的变换，但它比基2的算法速度更高。为了提高速度，很多FFT算法使用混合基算法。如我看到的2个效率很高程序均使用了混合基算法。第一个程序是来自：http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html，它使用了基2，基4,甚至基8混合算法，共提供了6同样功能的算法。可对长度为2^k的序列做FFT，这个程序的写的很复杂，我现在尚不能完全看懂。另一个程序来自：http://hjem.get2net.dk/jjn/fft.htm。相对于前者，这个程序相对简单一点。它使用了基2，基3,基4,基5,基8,基10   混合算法，几乎可以计算任意长度的FFT。具体的，当序列长度n为2，3,5,7,11,13,17,19,23,29,31,37等小素数时，或者n的最大素因数小于等于37时，可计算这个序列的FFT。

转载地址：http://mzyob.baihongyu.com/