FFT原理及实现

netcfln 发表于 2011-12-22 08:51

<DIV>
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198789111a1.jpg" target=_blank></A><A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_1311988002W998.jpg" target=_blank></A>两点DFT简化
假设输入为x,x;输出为X,X. 伪代码如下（简单的单个2点蝶形算法） :
// ------------------------------------------------------------------
#define   N     2
#define   PI    3.1415926
 
// ------------------------------------------------------------------
int i, j
for(i=0, X=0.0; i<N; i++)
   for(j=0; j<N; j++)
      X += x * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198776322q2.jpg" target=_blank></A> 
X =  x*(1-0) + x*(1-0) = x + 1*x;
X =  x*(1-0) + x*(-1-0) = x - 1*x;
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198766966qu.jpg" target=_blank></A> 
FFT实现流程图分析(N=8,以8点信号为例)
FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-point DFTs
 
<SPAN style="FONT-FAMILY: 'Times New Roman'; FONT-SIZE: 10.5pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-bidi-font-size: 10.0pt; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US> 
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198766966qu.jpg" target=_blank><IMG border=0 src="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198766966qu.jpg" ; .load="imgResize(this, 650);"></A>

 

8点FFT流程图(Layer表示层, gr表示当前层的颗粒)
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198776322q2.jpg" target=_blank></A> 
 
 
 
下面以LayerI为例:
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198776322q2.jpg" target=_blank><IMG border=0 src="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198776322q2.jpg" ; .load="imgResize(this, 650);"></A>
LayerI部分,具有4个颗粒,每个颗粒2个输入
(注意2个输入的来源,由时域信号友情提供,感谢感谢J)
<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">我们将输入<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>x<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">分为两部分<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>x_r, x_i.<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">具有实部和虚部<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>,<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">时域信号本没有虚部的<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>,<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">因此可以让<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>x_i<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">为<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>0.<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">那么为什么还要画蛇添足分为实部和虚部呢<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>?<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">这是因为<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>LayerII, LayerIII<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">的输入是复数<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>,<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">为了编码统一而强行分的<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>.<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">当然你编码时可以判断当前层是否为<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>1<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">来决定是否分<SPAN style="FONT-FAMILY: 'Times New Roman'; COLOR: black; FONT-SIZE: 10pt; mso-fareast-font-family: 宋体; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA" lang=EN-US>.<SPAN style="FONT-FAMILY: 宋体; COLOR: black; FONT-SIZE: 10pt; mso-fareast-language: ZH-CN; mso-font-kerning: 1.0pt; mso-ansi-language: EN-US; mso-bidi-language: AR-SA; mso-bidi-font-family: 'Times New Roman'">但是我想每个
人最后都会倾向分的.
旋转因子 tw = cos(2*PI*k/N)-j*sin(2*PI*k/N);也可以分为实部和虚部,令其为tw_r, tw_i;
则tw = tw_r - j*tw_i;
 
X = (x_r + j*x_i) + (tw_r–j*tw_i) * (x_r+j*x_i)
则
         X_R = x_r + tw_r*x_r + tw_i*x_i;
          X_I = x_i - tw_i*x_r + tw_r*x_i;
 
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198789111a1.jpg" target=_blank><IMG border=0 src="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198789111a1.jpg" ; .load="imgResize(this, 650);"></A>
LayerII部分,具有2个颗粒,每个颗粒4个输入
(注意4个输入的来源,由LayerI友情提供,感谢感谢J)
<A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_131198789111a1.jpg" target=_blank></A><A href="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_1311988002W998.jpg" target=_blank><IMG border=0 src="http://blog.chinaunix.nethttp://blog.chinaunix.net/attachment/201107/30/25977911_1311988002W998.jpg" ; .load="imgResize(this, 650);"></A>
LayerIII部分,具有1个颗粒,每个颗粒8个输入
(注意8个输入的来源,由LayerII友情提供,感谢感谢J)
LayerI, LayerII, LayerIII从左往右,蝶形信号运算流非常明显!
 
 
  
 
 
假令输入为x, x,输出为X, X. x分解为x_r, x_i部分
则该蝶形运算为
X
 = (x_r-j*x_i) + (x_r-j*x_i)*(cos(2*PI*k/N)-j*sin(2*PI*k/N));
再令cos(2*PI*k/N)为tw1, sin(2*PI*k/N)为tw2则
X = (x_r-j*x_i) + (x_r-j*x_i)*(tw1-j*tw2);
 
X_R = x_r + x_r*tw1 - x_i*tw2;
X_I = x_i
 
 
 
 
 
 
 
 
 
             x_r = x_r + x_r*tw1 + x_i*tw2;
             x_i = x_i - x_r*tw2 + x_i*tw1;
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
譬如8点输入x
1. 先分割成2部分: x, x, x, x和 x, x, x, x
2. 信号x, x, x, x再分割成x, x和 x, x
   信号x, x, x, x再分割成x, x和 x, x
3.无法分割了,已经分割成2点了J.
 
 
 
如上图:
在LayerI的时候,我们是对2点进行DFT.(一共4次DFT )
输入为    x&x;   x&x;    x&x;    x&x
输出为    y,y; Y,y; Y,y; Y,y;
 
流程:
I.希望将输入直接转换为x, x, x, x, x, x, x, x的顺序
II.对转换顺序后的信号进行4次DFT
 
步骤I代码实现
/**
 *反转算法. 这个算法效率比较低!先用起来在说,之后需要进行优化.
 */
static void bitrev(void )
{
   int   p=1, q, i;
   int   bit_rev[ N ];
   float  xx_r[ N ];
   
   bit_rev[ 0 ] = 0;
   while( p < N )
   {
      for(q=0; q<p; q++) 
      {
          bit_rev[ q ]    = bit_rev[ q ] * 2;
          bit_rev[ q + p ] = bit_rev[ q ] + 1;
      }
      p *= 2;
   }
   for(i=0; i<N; i++)  xx_r[ i ] = x_r[ i ];   
   for(i=0; i<N; i++)  x_r = xx_r[ bit_rev ];
}
// ------------------------此刻序列x重排完毕------------------------
 
步骤II代码实现
int j;
float TR;    //临时变量
float tw1; //旋转因子
/*两点DFT */
for(k=0; k<N; k+=2)
{
   //<A href="http://student.csdn.net/link.php?url=http://writeblog.csdn.net/Editor/FCKeditor/editor/fckeditor.html?InstanceName=ctl00_ContentPlaceHolder1_EntryEditor1_richTextEditor_richTextEditor&amp;Toolbar=Default#_top" target=_blank>两点DFT简化</A>告诉我们tw1=1
   TR = x_r; // TR就是A, x_r就是B.
   x_r  = TR + tw1*x_r;
   x_r = TR - tw1*x_r;
}
 
 
 
 
 
 
在LayerII的时候,我们希望得到z,就需要对y进行DFT.
y,y; y,y; y,y; y,y;
z,         z; z,z; z,z; z,z;
 
在LayerIII的时候,我们希望得到v,就需要对z进行DFT.
z,z; z,z; z,z; z,z;
v,v; v,v; v,v; v,v;
 
准备
 
 
令输入为x, x,输出为y, y
这个N绝对不是上面的8,这个N是当前颗粒的输入样本总量
对于LayerI而言N是2;对于LayerII而言N是4;对于LayerIII而言N是8
 
复数乘法:(a+j*b) * (c+j*d)
实部= a*c – bd;
虚部= ad + bc;
 
旋转因子:
 
 
 
 
 
 
 
 
实现(C描述)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#include "complex.h"
 
// --------------------------------------------------------------------------
#define  N    8 //64 
#define  M    3 //6    //2^m=N 
#define  PI   3.1415926
// --------------------------------------------------------------------------
 
float  twiddle = {1.0, 0.707, 0.0, -0.707};
float  x_r = {1, 1, 1, 1, 0, 0, 0, 0}; 
float  x_i;                        //N=8 
/*
float  twiddle = {1,      0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733, 
                       0.7075, 0.6349, 0.5561,   0.4721,   0.3835,   0.2912,   0.1961,   0.0991, 
                       0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349, 
                      -0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951};      //N=64 
float  x_r={1,1,1,1,1,1,1,1,
               1,1,1,1,1,1,1,1,
               1,1,1,1,1,1,1,1,
                             1,1,1,1,1,1,1,1, 
                             0,0,0,0,0,0,0,0, 
                             0,0,0,0,0,0,0,0, 
                             0,0,0,0,0,0,0,0, 
                             0,0,0,0,0,0,0,0,};
float  x_i;
*/
FILE *fp;
 
// ----------------------------------- func -----------------------------------
/**
 *初始化输出虚部
 */
static void fft_init(void )
{
   int i;
   for(i=0; i<N; i++)  x_i = 0.0; 
}
 
/**
 *反转算法.将时域信号重新排序.
 *这个算法有改进的空间
 */
static void bitrev(void )
{
   int   p=1, q, i;
   int   bit_rev[ N ]; //
   float  xx_r[ N ];   //
   
   bit_rev[ 0 ] = 0;
   while( p < N )
   {
      for(q=0; q<p; q++) 
      {
          bit_rev[ q ]    = bit_rev[ q ] * 2;
          bit_rev[ q + p ] = bit_rev[ q ] + 1;
      }
      p *= 2;
   }
   
   for(i=0; i<N; i++)  xx_r[ i ] = x_r[ i ];   
   
   for(i=0; i<N; i++)  x_r = xx_r[ bit_rev ];
}
 
/* ------------ add by sshc625 ------------ */
static void bitrev2(void )
{
   return ;
}
 
 
/* */ 
void display(void )
{
   printf("\n\n");
   int  i; 
   for(i=0; i<N; i++)
      printf("%f\t%f\n", x_r, x_i);
}
 
/**
 *
 */
void fft1(void ) 
{ fp = fopen("log1.txt","a+");
   int    L, i, b, j, p, k, tx1, tx2;
   float  TR, TI, temp; //临时变量
   float  tw1, tw2;
   
   /*深M.对层进行循环. L为当前层,总层数为M. */
   for(L=1; L<=M; L++)
   {
      fprintf(fp,"----------Layer=%d----------\n", L);
      /* b的意义非常重大,b表示当前层的颗粒具有的输入样本点数 */
      b = 1;
      i = L - 1;
      while(i > 0)
      {
          b *= 2;
          i--;
      }
      
// --------------是否外层对颗粒循环,内层对样本点循环逻辑性更强一些呢! --------------
      /*
       * outter对参与DFT的样本点进行循环
       * L=1,循环了1次(4个颗粒,每个颗粒2个样本点)
       * L=2,循环了2次(2个颗粒,每个颗粒4个样本点)
       * L=3,循环了4次(1个颗粒,每个颗粒8个样本点)
       */
      for(j=0; j<b; j++)
      {
          /*求旋转因子tw1 */
          p = 1;
          i = M - L; // M是为总层数, L为当前层.
          while(i > 0)
          {
             p = p*2;
             i--;
          }
          p  = p * j;
          tx1 = p % N;
          tx2 = tx1 + 3*N/4;
          tx2 = tx2 % N;
          // tw1是cos部分,实部; tw2是sin部分,虚数部分.
          tw1 = ( tx1>=N/2)? -twiddle  : twiddle[ tx1 ];
          tw2 = ( tx2>=N/2)? -twiddle : twiddle;
          
          /*
           * inner对颗粒进行循环
           * L=1,循环了4次(4个颗粒,每个颗粒2个输入)
           * L=2,循环了2次(2个颗粒,每个颗粒4个输入)
           * L=3,循环了1次(1个颗粒,每个颗粒8个输入)
           */
          for(k=j; k<N; k=k+2*b)
          {
             TR = x_r; // TR就是A, x_r就是B.
             TI = x_i;
             temp = x_r;
             /*
              *如果复习一下 (a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么
              *就能理解为什么会如下运算了,只有在L=1时候输入才是实数,之后层的
              *输入都是复数,为了让所有的层的输入都是复数,我们只好让L=1时候的
              *输入虚部为0
              * x_i*tw2是两个虚数相乘 
              */
             fprintf(fp,"tw1=%f, tw2=%f\n", tw1, tw2);
             x_r  = TR + x_r*tw1 + x_i*tw2;
             x_i  = TI - x_r*tw2 + x_i*tw1;
             
             x_r = TR - x_r*tw1 - x_i*tw2;
             x_i = TI + temp*tw2    - x_i*tw1;
             
             fprintf(fp,"k=%d, x_r=%f, x_i=%f\n", k, x_r, x_i);
             fprintf(fp,"k=%d, x_r=%f, x_i=%f\n", k+b, x_r, x_i);
          }//
      }//
   }//
}
 
/**
 * ------------ add by sshc625 ------------
 *该实现的流程为
 * for( Layer )
 *    for( Granule )
 *        for( Sample )
 * 
 *
 *
 *
 */
void fft2(void )
{  fp = fopen("log2.txt","a+");
   int    cur_layer, gr_num, i, k, p;
   float  tmp_real, tmp_imag, temp;  //临时变量,记录实部
   float  tw1, tw2;//旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.
      
   int   step;     //步进
   int   sample_num;  //颗粒的样本总数(各层不同,因为各层颗粒的输入不同)
   
   /*对层循环 */
   for(cur_layer=1; cur_layer<=M; cur_layer++)
   {     
      /*求当前层拥有多少个颗粒(gr_num) */
      gr_num = 1;
      i = M - cur_layer;
      while(i > 0)
      {
          i--;
          gr_num *= 2;
      }
      
      /*每个颗粒的输入样本数N' */
      sample_num   = (int)pow(2, cur_layer);
      /*步进.步进是N'/2 */
      step      = sample_num/2;
      
      /* */
      k = 0;
      
      /*对颗粒进行循环 */
      for(i=0; i<gr_num; i++)
      {
          /*
           *对样本点进行循环,注意上限和步进 
           */
          for(p=0; p<sample_num/2; p++)
          {  
             //旋转因子,需要优化...   
             tw1 = cos(2*PI*p/pow(2, cur_layer));
             tw2 = -sin(2*PI*p/pow(2, cur_layer));
             
             tmp_real = x_r;
             tmp_imag = x_i;
             temp = x_r;
             
             /*(tw1+jtw2)(x_r+jx_i)
              *
              * real : tw1*x_r - tw2*x_i
              * imag : tw1*x_i + tw2*x_r
              *我想不抽象出一个
              * typedef struct {
              * double real; //实部
              * double imag; //虚部
              * } complex;以及针对complex的操作
              *来简化复数运算是否是因为效率上的考虑!
              */
              
             /*蝶形算法 */
             x_r  = tmp_real + ( tw1*x_r - tw2*x_i );
             x_i  = tmp_imag + ( tw2*x_r + tw1*x_i );
             /* X = A(k)+WB(k)
              * X = A(k)-WB(k)的性质可以优化这里*/
             //旋转因子,需要优化...
             tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
             tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
             x_r = tmp_real + ( tw1*temp - tw2*x_i );
             x_i = tmp_imag + ( tw2*temp + tw1*x_i );
             
             printf("k=%d, x_r=%f, x_i=%f\n", k+p, x_r, x_i);
             printf("k=%d, x_r=%f, x_i=%f\n", k+p+step, x_r, x_i);
          }
          /*开跳!:) */
          k += 2*step;
      }  
   }
}
/*
 *后记:
 *究竟是颗粒在外层循环还是样本输入在外层,好象也差不多,复杂度完全一样.
 *但以我资质愚钝花费了不少时间才弄明白这数十行代码.
 *从中我发现一个于我非常有帮助的教训,很久以前我写过一部分算法,其中绝大多数都是递归.
 *将数据量减少,减少再减少,用归纳的方式来找出数据量加大代码的规律
 *比如FFT
 * 1.先写死LayerI的代码;然后再把LayerI的输出作为LayerII的输入,又写死代码; ......
 *   大约3层就可以统计出规律来.这和递归也是一样,先写死一两层,自然就出来了!
 * 2.有的功能可以写伪代码,不急于求出结果,降低复杂性,把逻辑结果定出来后再添加.
 *   比如旋转因子就可以写死,就写1.0.流程出来后再写旋转因子.
 *寥寥数语,我可真是流了不少汗! Happy!
 */
 
  
void dft(void )
{
   int   i, n, k, tx1, tx2;
   float  tw1,tw2;
   float  xx_r,xx_i;
   
   /*
    * clear any data in Real and Imaginary result arrays prior to DFT
    */
   for(k=0; k<=N-1; k++)
      xx_r = xx_i = x_i = 0.0;
      
   // caculate the DFT
   for(k=0; k<=(N-1); k++)
   {
      for(n=0; n<=(N-1); n++)
      {
          tx1 = (n*k);
          tx2 = tx1+(3*N)/4;
          tx1 = tx1%(N);
          tx2 = tx2%(N);
          if(tx1 >= (N/2))
             tw1 = -twiddle; 
          else 
             tw1 = twiddle; 
          if(tx2 >= (N/2)) 
             tw2 = -twiddle; 
          else 
             tw2 = twiddle;
          xx_r = xx_r+x_r*tw1; 
          xx_i = xx_i+x_r*tw2; 
      } 
      xx_i = -xx_i;
   }
   // display
   for(i=0; i<N; i++) 
      printf("%f\t%f\n", xx_r, xx_i); 
}
 
// ---------------------------------------------------------------------------
int main(void )
{
   fft_init( );
   bitrev( );
   // bitrev2( );
   //fft1( );
   fft2( );
   display( );
   
   system("pause" );
   // dft();               
   return 1;
}
 
 </DIV>

abc3214 发表于 2012-08-25 17:52

为什么还要画蛇添足分为实部和虚部呢

页: [1]

Chinaunix's Archiver

FFT原理及实现