doc

FFT算法原理与实现

  • 1星
  • 日期: 2018-10-09
  • 大小: 333KB
  • 所需积分:1分
  • 下载次数:5
  • favicon收藏
  • rep举报
  • 分享
  • free评论
标签: FFT算法

FFT算法

数字信号处理

数字信号处理,缩写为DSP,是面向电子信息学科的专业基础课,它的基本概念、基本分析方法已经渗透到了信息与通信工程,电路与系统,集成电路工程,生物医学工程,物理电子学,导航、制导与控制,电磁场与微波技术,水声工程,电气工程,动力工程,航空工程,环境工程等领域。

在数字信号处理中常常需要用到离散傅立叶变换(DFT),以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征,但是算法计算量大,耗时长,不利于计算机实时对信号进行处理。因此至DFT被发现以来,在很长的一段时间内都不能被应用到实际的工程项目中,直到一种快速的离散傅立叶计算方法——FFT,被发现,离散傅立叶变换才在实际的工程中得到广泛应用。需要强调的是,FFT并不是一种新的频域特征获取方式,而是DFT的一种快速实现算法。本文就FFT的原理以及具体实现过程进行详尽讲解。

FFT原理与实现   2010-10-07 21:10:09|  分类: 数字信号处理 |  标签:fft  dft  |举报|字号 订阅 在数字信号处理中常常需要用到离散傅立叶变换(DFT),以获取信号的频域特征。尽 管传统的DFT算法能够获取信号频域特征,但是算法计算量大,耗时长,不利于计算机实 时对信号进行处理。因此至DFT被发现以来,在很长的一段时间内都不能被应用到实际的 工程项目中,直到一种快速的离散傅立叶计算方法——FFT,被发现,离散傅立叶变换才在 实际的工程中得到广泛应用。需要强调的是,FFT并不是一种新的频域特征获取方式,而 是DFT的一种快速实现算法。本文就FFT的原理以及具体实现过程进行详尽讲解。 DFT计算公式 本文不加推导地直接给出DFT的计算公式: [pic] 其中x(n)表示输入的离散数字信号序列,WN为旋转因子,X(k)为输入序列x(n)对应的 N个离散频率点的相对幅度。一般情况下,假设x(n)来自于低通采样,采样频率为fs,那 么X(k)表示了从-fs/2率开始,频率间隔为fs/N,到fs/2- fs/N截至的N个频率点的相对幅度。因为DFT计算得到的一组离散频率幅度值实际上是在 频率轴上从成周期变化的,即X(k+N)=X(k)。因此任意取连续的N个点均可以表示DFT的计 算效果,负频率成分比较抽象,难于理解,根据X(k)的周期特性,于是我们又可以认为 X(k)表示了从零频率开始,频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。 N点DFT的计算量 根据(1)式给出的DFT计算公式,我们可以知道每计算一个频率点X(k)均需要进行N次 复数乘法和N-1次复数加法,计算N各点的X(k)共需要N^2次复数乘法和N*(N- 1)次复数加法。当x(n)为实数的情况下,计算N点的DFT需要2*N^2次实数乘法,2*N*(N- 1)次实数加法。 旋转因子WN的特性   1.WN的对称性 [pic]        2.WN的周期性 [pic]    3.WN的可约性 [pic]根据以上这些性质,我们可以得到式(5)的一系列有用结果 [pic] 基-2 FFT算法推导 假设采样序列点数为N=2^L,L为整数,如果不满足这个条件可以人为地添加若干个0 以使采样序列点数满足这一要求。首先我们将序列x(n)按照奇偶分为两组如下: [pic]于是根据DFT计算公式(1)有: [pic]至此,我们将一个N点的DFT转化为了式(7)的形式,此时k的取值为0到N- 1,现在分为两段来讨论,当k为0~N/2- 1的时候,因为x1(r),x2(r)为N/2点的序列,因此,此时式(7)可以写为: [pic]而当 k取值为N/2~N-1时,k用k’+N/2取代,k’取值为0~N/2- 1。对式(7)化简可得: [pic] 综合以上推导我们可以得到如下结论:一个N点的DFT变换过程可以用两个N/2点的DF T变换过程来表示,其具体公式如式(10)所示DFT快速算法的迭代公式: [pic]上式中X'(k’)为偶数项分支的离散傅立叶变换,X''(k’’)为奇数项分支的离散 傅立叶变换。 式(10)的计算过程可以用图1的蝶形算法流图直观地表示出来。 [pic]图1 时间抽取法蝶形运算流图 在图1中,输入为两个N/2点的DFT输出为一个N点的DFT结果,输入输出点数一致。运用这 种表示方法,8点的DFT可以用图2来表示: [pic] 图2 8点DFT的4点分解 根据公式(10),一个N点的DFT可以由两个N/2点的DFT运算构成,再结合图1的蝶形信号流 图可以得到图2的8点DFT的第一次分解。该分解可以用以下几个步骤来描述:   1.将N点的输入序列按奇偶分为2组分别为N/2点的序列   2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2   3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果 根据式(10)我们可以对图2中的4点DFT进一步分解,得到图3的结果,分解步骤和前面一 致。 [pic] 图3 8点DFT的2点分解 最后对2点DFT进一步分解得到最终的8点FFT信号计算流图: [pic] 图4 8点DFT的全分解 从图2到图4的过程中关于旋转系数的变化规律需要说明一下。看起来似乎向前推一级, 在奇数分组部分的旋转系数因子增量似乎就要变大,其实不是这样。事实上奇数分组部 分的旋转因子指数每次增量固定为1,只是因为每向前推进一次,该分组序列的数据个数 变少了,为了统一使用以原数据N为基的旋转因子就进行了变换导致的。每一次分组奇数 部分的系数WN,这里的N均为本次分组前的序列点数。以上边的8点DFT为例,第一次分组 N=8,第二次分组N为4,为了统一根据式(4)进行了变换将N变为了8,但指数相应的需要乘 以2。 N点基-2 FFT算法的计算量 从图4可以看到N点DFT的FFT变换可以转为log2(N)级级联的蝶形运算,每一级均包含有N /2次蝶形计算。而每一个蝶形运算包含了1次复数乘法,2次复数加法。因此N点FFT计算 的总计算量为:复数乘法——N/2×log2(N)    复数加法——N×log2(N)。假设被采样的序列为实数序列,那么也只有第一级的计算为实数 与复数的混合计算,经过一次迭代后来的计算均变为复数计算,在这一点上和直接的DF T计算不一致。因此对于输入序列是复数还是实数对FFT算法的效率影响较小。一次复数 乘法包含了4次实数乘法,2次实数加法,一次复数加法包含了2次复数加法。因此对于N 点的FFT计算需要总共的实数乘法数量为:2×N×log2(N);总的复数加法次数为:2xNxlo g2(N)。 N点基-2 FFT算法的实现方法 从图4我们可以总结出对于点数为N=2^L的DFT快速计算方法的流程:    1.对于输入数据序列进行倒位序变换。 该变换的目的是使输出能够得到X(0)~X(N- 1)的顺序序列,同样以8点DFT为例,该变换将顺序输入序列x(0)~x(7)变为如图4的x(0) ,x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其实现方法是:假设顺序输入序列一次村 在A(0)~A(N- 1)的数组元素中,首先我们将数组下标进行二进制化(例:对于点数为8的序列只需要LO G2(8) = 3位二进制序列表示,序号6就表示为110)。二进制化以后就是将二进制序列进行倒位, 倒位的过程就是将原序列从右到左书写一次构成新的序列,例如序号为6的二进制表示为 110,倒位后变为了011,即使十进制的3。第三步就是将倒位前和倒位后的序号对应的数 据互换。依然以序号6为例,其互换过程如下: temp = A(6); A(6) = A(3); A(3) = A(6); 实际上考虑到执行效率,如果对于每一次输入的数据都需要这个处理过程是非常浪费时 间的。我们可以采用指向指针的指针来实现该过程,或者是采用指针数组来实现该过程 。   2.蝶形运算的循环结构。  从图4中我们可以看到对于点数为N = 2^L的fft运算,可以分解为L阶蝶形图级联,每一阶蝶形图内又分为M个蝶形组,每个蝶 形组内包含K个蝶形。根据这一点我们就可以构造三阶循环来实现蝶形运算。编程过程需 要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。   3.浮点到定点转换需要注意的关键问题  上边的分析都是基于浮点运算来得到的结论,事实上大多数嵌入式系统对浮点运算支持 甚微,因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点方式。对于简单的 DFT运算从浮点到定点显得非常容易。根据式(1),假设输入x(n)是经过AD采样的数字序 列,AD位数为12位,则输入信号范围为0~4096。为了进行定点运算我们将旋转因子实部 虚部同时扩大2^12倍,取整数部分代表旋转因子。之后,我们可以按照(1)式计算,得到 的结果与原结果成比例关系,新的结果比原结果的2^12倍。但是,对于使用蝶形运算的 fft我们不能采用这种简单的放大旋转因子转为整数计算的方式。因为fft是一个非对称 迭代过程,假设我们对旋转因子进行了放大,根据蝶形流图我们可以发现其最终的结果 是,不同的输入被放大了不同的倍数,对于第一个输入x(0)永远也不会放大。举一个更 加形象的例子,还是以图4为例。从图中可以看出右侧的X(0)可以直接用下式表示:  [pic] 从上式我们可以看到不同输入项所乘的旋转因子个数(注意这里是个数,就算是wn^0,也 被考虑进去了,因为在没有放大时wn^0等于1,放大后所有旋转因子指数模均不为1,因 此需要考虑)。这就导致输入不平衡,运算结果不正确。经查阅相关资料,比较妥善的做 法是,首先对所有旋转因子都放大2^Q倍,Q必须要大于等于L,以保证不同旋转因子的差 异化。旋转因子放大,为了保证其模为1,在每一次蝶形运算的乘积运算中我们需要将结 果右移Q位来抵消这个放大,从而得到正确的结果。之所以采用放大倍数必须是2的整数 次幂的原因也在于此,我们之后可以通过简单的右移位运算将之前的放大抵消,而右移 位又代替了除法运算,大大节省了时间。          4.计算过程中的溢出问题  最后需要注意的一个问题就是计算过程中的溢出问题。在实际应用中,AD虽然有12位的 位宽,但是采样得到的信号可能较小,例如可能在0~8之间波动,也就是说实际可能只有 3位的情况。这种情况下为了在计算过程中不丢失信息,一般都需要先将输入数据左移P 位进行放大处理,数据放大可能会导致溢出,从而使计算错误,而溢出的极限情况是这 样:假设我们数据位宽为D位(不包括符号位),AD采样位数B位,数字放大倍数P位,旋转 因此放大倍数Q位,FFT级联运算带来的最大累加倍数L位。我们得到:  [pic]  假设AD位宽12,数据位宽32,符号位1位,因此有效位宽31位,采样点数N,那么我们可 以得到log2(N)+P+Q<=19,假设点数128,又Q>=L可以得到放大倍数P<=5。  FFT代码示例  #define N 128 #define POWER 6//该值代表了输入数据首先被放大的倍数,放大倍数为2^POWER #define P_COEF 8//该值代表了旋转因子被放大的倍数,放大倍数为2^POWER #if (N == 4) #define L 2//L的定义满足L = log2(N) #elif (N == 8) #define L 3//L的定义满足L = log2(N) #elif (N == 16) #define L 4//L的定义满足L = log2(N) #elif (N == 32) #define L 5//L的定义满足L = log2(N) #elif (N == 64) #define L 6//L的定义满足L = log2(N) #elif (N == 128) #define L 7//L的定义满足L = log2(N) #endif //AD采样位数为12位,本可以采用s16 x[N]定义数据空间,但是为了节省存储空间,fft结果也将存储在该变量空间。由于受 //fft影响变量需要重新定义,定义的方式及具体原因如下: //fft过程会乘以较大旋转因子,因此需要32位定义 //fft过程会产生负数,因此需要有符号定义 //fft过程会出现复数,因此需要定义二维数组,[][0]存放实部,[][1]存放虚部 s32 x[N][2] = {}; //定义*p[N]是为了在第一次指针初始化以后,该数组指针按照位倒序指向数据变量x //p[i][0]代表了指向数据的实部,p[i][1]代表了指向数据的虚部 s32 *p[N]; //定义旋转因子矩阵 //旋转因子矩阵存储了wn^0,wn^1,wn^2...wn^(N/2- 1)这N/2个复数旋转因子 s16 w[N>>1][2] = {256,0,256,-13,255,-25,253,-38,251,- 50,248,-62,245,-74,241,-86,237,-98,231,-109,226,                   -121,220,-132,213,-142,206,-152,198,-162,190,-172,181,- 181,172,-190,162,-198,152,                   -206,142,-213,132,-220,121,-226,109,-231,98,-237,86,- 241,74,-245,62,-248,50,-251,38,                   -253,25,-255,13,-256,0,-256,-13,-256,-25,-255,-38,-253,- 50,-251,-62,-248,-74,-245,-86,                   -241,-98,-237,-109,-231,-121,-226,-132,-220,-142,-213,- 152,-206,-162,-198,-172,-190,-181,                   -181,-190,-172,-198,-162,-206,-152,-213,-142,-220,-132,- 226,-121,-231,-109,-237,-98,-241,                   -86,-245,-74,-248,-62,-251,-50,-253,-38,-255,-25,-256,- 13}; void init_pointer(void) { unsigned char i = 0; unsigned char j = 0; unsigned char k = 0; for(i = 0; i < N; i++) { j = 0; for(k = 0; k < L; k++) { j |= (((i >> k)&0x01)<<(L-k-1)); } p[i] = &x[j][0]; } } /* *description:基2fft主函数,该函数将借助指针数组p对全局变量数组x进行fft计算 * 并将结果存储在数组x中 *global var:rw->x; r- >p,w。(r表示读,w表示写,rw表示读写) */ void fft2(void) { u8 i;//i用于表示蝶形图级联的阶数 u8 j;//表示蝶形分组起始点序列,蝶形分组跨度为2^i u8 k;//k表示蝶形组内偶数分支序列点号 u8 gp_distance = 1;//蝶形分组跨度 u8 temp;//temp用于记录当前组间距离的一半,同时也表示了同一碟形两分支间的距离 u8 gp_hf = 0;//记录当前组的中间点序号 u8 delta = N;//旋转因子下标增量,本来下标初始值应该为N/2,但是。。 s16 *pw = &(w[0][0]); int tp_result[2]; //用于临时存放旋转因子和奇数分组的乘积 //输入信号序列放大 for(i = 0; i < N; i++) { x[i][0] <<= POWER; x[i][1] <<= POWER; } for(i = 0; i < L; i++) { temp = gp_distance; gp_distance <<= 1; for(j = 0; j < N; j+=gp_distance) { gp_hf = temp + j; pw = &(w[0][0]); for(k = j; k < gp_hf; k++)//完成一组内的所有蝶形运算 { //蝶形运算中的一组复数乘法过程 tp_result[0] = pw[0] * (p[k+temp][0]) - pw[1] * (p[k+temp][1]); tp_result[1] = pw[0] * (p[k+temp][1]) + pw[1] * (p[k+temp][0]); tp_result[0] >>= P_COEF; tp_result[1] >>= P_COEF; //蝶形运算中的2组复数加法过程 p[k+temp][0] = p[k][0] - tp_result[0]; p[k+temp][1] = p[k][1] - tp_result[1]; p[k][0] += tp_result[0]; p[k][1] += tp_result[1]; pw += delta; } } delta >>= 1; } }
更多简介内容

推荐帖子

评论

combat
一篇技术性文章,得研究一下才能弄明白。
2020-10-24 14:41:38回复
登录/注册

意见反馈

求资源

回顶部

datasheet推荐 换一换

About Us 关于我们 客户服务 联系方式 器件索引 网站地图 最新更新 手机版 版权声明

北京市海淀区知春路23号集成电路设计园量子银座1305 电话:(010)82350740 邮编:100191

电子工程世界版权所有 京ICP证060456号 京ICP备10001474号 电信业务审批[2006]字第258号函 京公海网安备110108001534 Copyright © 2005-2020 EEWORLD.com.cn, Inc. All rights reserved
$(function(){ var appid = $(".select li a").data("channel"); $(".select li a").click(function(){ var appid = $(this).data("channel"); $('.select dt').html($(this).html()); $('#channel').val(appid); }) })