加载中…
正文 字体大小:

Cooley-Tukey

(2012-12-02 15:46:48)
标签:

杂谈

Cooley-Tukey是FFT的一种快速算法模型,如下图中,左下角为幅值图,右下角为相角图,将其按奇偶分为两个序列,再递归依次计算子序列,其变换的结果与图中心点对称:

Cooley-Tukey

具体路线很像蝶形,首先要将初始序列按序列号逆序,在8个序列的情况下,如001变为100,则序列1与序列4交换,然后对逆序后的序列进行log(N)次迭代计算,其迭代就是交换序列对的值,以基2时间域为例,具体蝶形变换流程图如下:

Cooley-Tukey

FFT以k为基只能计算k^n的尺寸的图像,采用补零的方法保证k^n的尺寸并不能计算到准确的结果,因为尺寸改变后频域中心和系数变化了,当采用混合基的方法可以涵盖大部分尺寸的图像,如采用基2,基3,基5,基8的混合基。以基2为例,C++代码如下:

//Created by pritry

void fft(CComplex* pData, int n, int nFlag)
{
 if(pData == NULL || n < 0) return;

 int i, j, k, hn, M, L, B, p, J;

 CComplex* Wn = NULL;
 CComplex* x = NULL;

 //强制为2的幂次方
 int d = 1;
 while(d < n) d <<= 1;
 int N = d;

 x = new CComplex[N];

 memset(x, 0, sizeof(CComplex)*N);
 memcpy(x, pData, sizeof(CComplex)*n);

 //初始序列位逆序
 hn = N >> 1;
 j = 0;

 for(i = 0; i < N - 1; ++i)
 {
  if(i < j)
  {
   CComplex tmp = x[i];
   x[i] = x[j];
   x[j] = tmp;
  }

  k = hn;

  while(j >= k)
  {
   j -= k;
   k >>= 1;
  }

  j += k;
 }

 //蝶形运算
 Wn = new CComplex[hn];
 M = (int)(log(N+1.0f)/log(2.0f));

 for(i = 0; i < hn; ++i)
 {
  double angle = nFlag*2*M_PI*i/n;
  Wn[i] = CComplex(cos(angle), sin(angle));
 }

 for(L = 1; L <= M; ++L)
 {
  B = 1 << (L-1);
  
  for(J = 0; J <= B-1; ++J)
  {
   p = J*(1<<(M-L));

   for(k = J; k <= n-1; k += (1 << L))
   {
    CComplex tmp1 = x[k] + x[k+B]*Wn[p];
    CComplex tmp2 = x[k] - x[k+B]*Wn[p];

    x[k]   = tmp1;
    x[k+B] = tmp2;
   }
  }
 }

 if(nFlag == 1)
 {
  for(int k = 0; k < N; ++k)
  {
   x[k] /= n;
  }
 }

 memcpy(pData, x, sizeof(CComplex)*n);

 delete [] Wn;
 delete [] x;
}

 

0

阅读 评论 收藏 转载 喜欢 打印举报
已投稿到:
  • 评论加载中,请稍候...
发评论

       

    发评论

    以上网友发言只代表其个人观点,不代表新浪网的观点或立场。

      

    新浪BLOG意见反馈留言板 不良信息反馈 电话:4006900000 提示音后按1键(按当地市话标准计费) 欢迎批评指正

    新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 会员注册 | 产品答疑

    新浪公司 版权所有