Chinaunix首页 | 论坛 | 博客

steven_miao的ChinaUnix博客miaoliangjiuhon.blog.chinaunix.net

首页　| 　博文目录　| 　关于我

steven_miao

  • 博客访问： 1829422
  • 博文数量： 195
  • 博客积分： 4227
  • 博客等级： 上校
  • 技术积分： 2835
  • 用 户 组： 普通用户
  • 注册时间： 2010-09-04 10:39
文章分类

全部博文（195）

文章存档

2013年（1）

2012年（26）

2011年（168）

我的朋友
最近访客
推荐博文
相关博文
FFT初解

分类： LINUX

2011-05-26 10:58:09

初解

谨以此文纪念过往的岁月

一．前言

首先申明俺不是一个算法工程师，俺是一个底层驱动工程师，有人会发问一个底层驱动工程师需要这个吗？但是我不幸的告诉你，确实是需要的，不过我们不要像算法工程师那样搞得很精通，但是还是需要去了解这是个什么东西。说实话，这个东西在大学时候学过，还好好的去理解了一样，不过到现在忘的差不多了，这愈发的让我明白一句话，好记性不如烂笔头，如果以前有好好记录的好习惯，那现在只要把以前的东西拿出来看看再印证一下就可以了。不过历史不可以如果。为了不让明天继续懊悔今天，在这里记录下本人学习的一些记录。

二．FFT的数学基础

FFT是什么，额，说白了就是一种时域和频域变换的手段。什么是时域什么是频域，额，你google吧，鄙视百度！

这里说明一下，相比较FFT而言，我们更加关心离散傅里叶变换(DTF),因为我们底层驱动工程师面对的是一个一个离散的数据。也许大家还是对FFT的变换的概念还不理解，那我们举个例子来说吧：

1 AD在一个1Hz正弦波周期采集了1024个数据(即你的采样频率为1024Hz)，从时域上看1s内采集了1024数据，那1024个数据做1024点的FFT变换。我们以数据的下标作为横轴，而数据的值作为纵轴，你会发现第一个点的值最大，我们可以从该值计算出我们被采样的频率为1Hz，值是这样算出来的1024(Hz)/1024*1 = 1Hz,如果最大的值点是第8个点，那被采样频率是8Hz。不过要满足乃奎斯定律，即采样频率大于被采样频率的两倍。

为什么会有东西出现呢，主要是时域的信号好采样，至于频域的信号难以采样，于是傅里叶这个天才就发明了这个东东，不过这可是现代数字信号处理的基础。

闲话不说了，还是来看DFT的数学基础(下面的内容不少拷贝wiki的)：

对于复数序列x_{0},. x_{1},. ...,. x_{n-1}，离散傅里叶变换公式为：

y_j = .sum_{k=0}^{n-1} e^{-{2.pi.imath .over n} j k}x_k .qquad j = 0,1,.dots,n-1.

下面，我们用N次单位根WN来表示e^{-j.frac{2.pi}{N}}。

WN的性质：

  1. 周期性，WN具有周期N。
  2. 对称性：W_{N}^{k+.frac{N}{2}}=-W_{N}^{k}。
  3. W_{N}^{ikn}=W_{.frac{N}{i}}^{kn}

为了简单起见，我们下面设待变换序列长度n = 2r。 根据上面单位根的对称性，求级数y_k=.sum_{n=0}^{N-1} W_{N}^{kn}x_n时，可以将求和区间分为两部分：

.begin{matrix}y_k=.sum_{n=2i} W_{N}^{kn} x_n + .sum_{n=2i+1} W_{N}^{kn}x_n..= .sum_{i} W_{.frac{N}{2}}^{ki}x_{2i} + W_{N}^{k}.sum_{i} W_{.frac{N}{2}}^{ki}x_{2i+1}..= F_{even}(k) + W_{N}^{k}F_{odd}(k)&&&&&&(i.in.mathbb{Z}).end{matrix}

Fodd(k) 和 Feven(k)是两个分别关于序列.left.{x_n.right.}_0^{N-1}奇数号和偶数号序列N/2点变换。由此式只能计算出yk的前N/2个点，对于后N/2个点，注意 Fodd(k) 和 Feven(k) 都是周期为N/2的函数，由单位根的对称性，于是有以下变换公式：

  • y_{k+.frac{N}{2}} = F_{even}(k) - W_{N}^{k}F_{odd}(k)
  • y_k = F_{even}(k) + W_{N}^{k}F_{odd}(k)。

这样，一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据不难分析出此时算法的时间复杂度为O(NlogN)

三．8点快速傅里叶变换

主要是蝶形算法

FFT算法图



__no_init float fft_r[FFT_N]; //FFT输入实数序列，保存变换后的频域实部
__no_init float fft_i[FFT_N]; //保存变换后的频域虚部
__no_init float harm[HARM_M+1]; //各次谐波幅值
const Uint8 FFT_BIT[8]={ //供完成倒位序排列使用的位定义
0x01,0x02,0x04,0x08,0x10,0x20,0x40,0x80};

  1. void fft(void)
  2. {
  3.     Uint16 i,j,k,L;
  4.     Uint16 b,c,p,n;
  5.     Uint8 x,xc;
  6.     float Gr,Gi,Wr,Wi; for(i=0;i<FFT_N;i++)
  7.     {
  8.         x =i;
  9.         xc=0;
  10.         for(j=0;j<FFT_M;j++)
  11.         {
  12.             if((x&0x01)!=0)
  13.             {
  14.                 xc|=FFT_BIT[(FFT_M-1)-j];
  15.             }
  16.             x>>=1;
  17.         }
  18.         fft_i[xc]=fft_r[i];
  19.     }
  20.     for(i=0;i<FFT_N;i++)
  21.     {
  22.         fft_r[i]=fft_i[i];
  23.         fft_i[i]=0;
  24.     }

  26. //以上是倒叙就是将数组中的数据进行倒叙，即将数组序号倒序，如FFT_BIT[5]存储到FFT_BIT[3],为什么呢？5 = 0x110，进行二进制倒转变为0x011,变为3.

  28.     for(L=1;L<=FFT_M;L++)
  29.     {
  30.         b=1;
  31.         p=FFT_N>>1;
  32.         for(i=0;i<L-1;i++)
  33.         {
  34.             b<<=1; //b=2^(L-1)
  35.             p>>=1; //p=N/(2^L)
  36.         }
  37.         c=b<<1; //c=b*2
  38.         for(j=0;j<b;j++)
  39.         {
  40.             n=p*j;
  41.             for(k=j;k<FFT_N;k+=c)
  42.             {
  43.                 Gr = fft_r[k];
  44.                 Gi = fft_i[k];
  45.                 Wr = Sinf[n+FFT_N/4]*fft_r[k+b] + Sinf[n]*fft_i[k+b];
  46.                 Wi = Sinf[n+FFT_N/4]*fft_i[k+b] - Sinf[n]*fft_r[k+b];
  47.                 fft_r[k] = Gr + Wr;
  48.                 fft_i[k] = Gi + Wi;
  49.                 fft_r[k+b] = Gr - Wr;
  50.                 fft_i[k+b] = Gi - Wi;
  51.             }
  52.         }
  53.     }
  54. }

//以上循环分为三个一个根据变换点的阶数循环，即上图中level1->level2->level3。第二个循环是group步进，level1为2  level 2为4。 第三个循环是需要两两相乘数据之间的，如level1的length为1，level2的length为2，level3的length为4。

关于旋转因子，就是cos和sin的值，不过这两者之间是可以互相转换的。

也许到此大家对于快速傅里叶有了大概的理解，下面数C++的code，也是源于网上。这段code更好的描述的FFT的蝶形算法，值得参考。


  1. void FFT(unsigned long & ulN, vector<complex<double> >& vecList)

  3. {

  5.        //得到幂数

  7.        unsigned long ulPower = 0; //幂数

  9.        unsigned long ulN1 = ulN - 1;

  11.        while(ulN1 > 0)

  13.        {

  15.               ulPower++;

  17.               ulN1 /= 2;

  19.        }

  21.        //反序

  23.        bitset<sizeof(unsigned long) * 8> bsIndex; //二进制容器

  25.        unsigned long ulIndex; //反转后的序号

  27.        unsigned long ulK;

  29.        for(unsigned long p = 0; p < ulN; p++)

  31.        {

  33.               ulIndex = 0;

  35.               ulK = 1;

  37.               bsIndex = bitset<sizeof(unsigned long) * 8>(p);

  39.               for(unsigned long j = 0; j < ulPower; j++)

  41.               {

  43.                      ulIndex += bsIndex.test(ulPower - j - 1) ? ulK : 0;

  45.                      ulK *= 2;

  47.               }

  49.  

  51.               if(ulIndex > p)

  53.               {

  55.                      complex<double> c = vecList[p];

  57.                      vecList[p] = vecList[ulIndex];

  59.                      vecList[ulIndex] = c;

  61.               }

  63.        }

  65.        //计算旋转因子

  67.        vector<complex<double> > vecW;

  69.        for(unsigned long i = 0; i < ulN / 2; i++)

  71.        {

  73.               vecW.push_back(complex<double>(cos(2 * i * PI / ulN) , -1 * sin(2 * i * PI / ulN)));

  75.        }

  77.  

  79.        for(unsigned long m = 0; m < ulN / 2; m++)

  81.        {

  83.               cout<< "\nvW[" << m << "]=" << vecW[m];

  85.        }

  87.  

  89.        //计算FFT

  91.        unsigned long ulGroupLength = 1; //段的长度

  93.        unsigned long ulHalfLength = 0; //段长度的一半

  95.        unsigned long ulGroupCount = 0; //段的数量

  97.        complex<double> cw; //WH(x)

  99.        complex<double> c1; //G(x) + WH(x)

  101.        complex<double> c2; //G(x) - WH(x)

  103.        for(unsigned long b = 0; b < ulPower; b++)

  105.        {

  107.               ulHalfLength = ulGroupLength;

  109.               ulGroupLength *= 2;

  111.               for(unsigned long j = 0; j < ulN; j += ulGroupLength)

  113.               {

  115.                      for(unsigned long k = 0; k < ulHalfLength; k++)

  117.                      {

  119.                             cw = vecW[k * ulN / ulGroupLength] * vecList[j + k + ulHalfLength];

  121.                             c1 = vecList[j + k] + cw;

  123.                             c2 = vecList[j + k] - cw;

  125.                             vecList[j + k] = c1;

  127.                             vecList[j + k + ulHalfLength] = c2;

  129.                      }

  131.               }

  133.        }

  135. }

四．总结

FFT数字信号处理的重中之重。

阅读(14753) | 评论(1) | 转发(3) |
1

上一篇：Cache 之解

下一篇：The Discrete Cosine Transform (DCT)

给主人留下些什么吧！~~

小雅贝贝2011-05-26 18:10:01

看不到图片啊。

回复 | 举报
关于我们 | 关于IT168 | 联系方式 | 广告合作 | 法律声明 | 免费注册

Copyright 2001-2010 ChinaUnix.net All Rights Reserved 北京皓辰网域网络信息技术有限公司. 版权所有

感谢所有关心和支持过ChinaUnix的朋友们

16024965号-6