FFT快速傅里叶变换(蝶形算法)详解

来源:百度  [  文档由 瑾懿诗快 贡献   ]  责编:杨丽  |  侵权/违法举报

1、二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下:先对各行逐一进行一维FFT,然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示:for(int i=0;i;i+)FFT_1D(ROW[i],N);for(int j=0;j;j+)FFT_1D(COL[j],M);其中,ROW[i]表示矩阵的第i行。注意这只是一个简单的记法,并不能完全照抄。还需要通过一些语句来生成各行的数据。同理,COL[i]是对矩阵的第i列的一种简单表示方法。所以,关键是一维FFT算法的实现。2、例程:includeincludeincludedefine N 1000定义复数类型*/typedef struct{double real;double img;}complex;complex x[N],*W;输入序列,变换核*/int size_x=0;输入序列的大小,在本程序中仅限2的次幂*/double PI;圆周率*/void fft();快速傅里叶变换*/void initW();初始化变换核*/void change();变址*/void add(complex,complex,complex*);复数加法*/void mul(complex,complex,complex*);复数乘法*/void sub(complex,complex,complex*);复数减法*/void output();int main(){int i;输出结果*/system("cls");PI=atan(1)*4;printf("Please input the size of x:\\n");scanf("%d",&size_x);printf("Please input the data in x[N]:\\n");for(i=0;i;i+)scanf("%lf%lf",&x[i].real,&x[i].img);initW();fft();output();return 0;}快速傅里叶变换*/void fft(){int i=0,j=0,k=0,l=0;complex up,down,product;change();for(i=0;i(size_x)/log(2);i+){/*一级蝶形运算*/l=1;for(j=0;j;j+2*l){/*一组蝶形运算*/for(k=0;k;k+){/*一个蝶形运算*/mul(x[j+k+l],W[size_x*k/2/l],&product);add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}初始化变换核*/void initW(){int i;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;i;i+){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}变址计算,将x(n)码位倒置*/void change(){complex temp;unsigned short i=0,j=0,k=0;double t;for(i=0;i;i+){k=i;j=0;t=(log(size_x)/log(2));while((t-)>0){j=j;j|=(k&1);k=k>>1;}if(j>i){temp=x[i];x[i]=x[j];x[j]=temp;}}}输出傅里叶变换的结果*/void output(){int i;printf("The result are as follows\\n");for(i=0;i;i+){printf("%.4f",x[i].real);if(x[i].img>=0.0001)printf("+.4fj\\n",x[i].img);else if(fabs(x[i].img))printf("\\n");else printf("%.4fj\\n",x[i].img);}}void add(complex a,complex b,complex*c){c->real=a.real+b.real;c->img=a.img+b.img;}void mul(complex a,complex b,complex*c){c->real=a.real*b.real-a.img*b.img;c->img=a.real*b.img+a.img*b.real;}void sub(complex a,complex b,complex*c){c->real=a.real-b.real;c->img=a.img-b.img;}www.egvchb.cn防采集请勿采集本网。

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

点击图片看大图

因转码可能存在排版等问题,敬请谅解!以下文字仅供您参考:

第一章快速傅里叶变换(FFT)4.1填空题(1)如果序列是一长度为64点的有限长序列,序列是一长度为128点的有限长序列,记(线性卷积),则为点的序列,如果采用基算法以快速卷积的方式实现线性卷积,则的

第五章 快速傅里叶变换 本章目录? ? ? ? ?

给出了蝶形运算误差分析模型,利用FFT信号流图的特点,针对截断、舍入和收敛舍入3种量化方法,得到了准确的定点和块浮点两种FFT算法的均方误差上下限.最后给出了噪信比结果,并用Matlab对其进行了仿真,结果

直接计算DFT的问题及改进的途径 按时间抽取的基2-FFT算法 按频率抽取的基2-FFT算法 快速傅里叶逆变换(IFFT)算法 Matlab实现2 5.1 引言?

因为频域抽样函数,反变换回来时域就是方波) 序列福利叶变换的关系是特殊的\"离散傅立叶变换\",也就是时域序列被认为是各种方波抽样信号的叠加,认为复数的角度只取0和∏这两种情况,于是你就看到

DFT在实际应用中很重要: 可以计算信号的频

变址计算,库里-图基算法 for(i=0;i;i+){ k=i;j=0;t=ex;while((t-)>0){ j;j|=k&1;k>>=1;} if(j>=i){ y[i]=x[j];y[j]=x[i];} } 用变址后的y向量进行计算 for(i=0;i;i+){ t=1;for(j=0;j;j+t){ for(k=0;k;k+)

谱、功率谱和线性卷积等。?

通常用图1中蝶形算法的信号流图来表示式⑸的离散傅里叶变换运算。例如,N=8=2的抽样点的信号序列x(n)的离散傅里叶变换,可用如图2所示的FET算法的信号流图来计算。① N=2点的离散傅里叶变换的计算全由蝶

直接按DFT变换进行计算,当序列长度N很

大时,计算量非常大,所需时间会很长。

FFT并不是一种与DFT不同的变换,而是 DFT的一种快速计算的算法。?3 5.2 直接计算DFT的问题及改进的途径?

DFT的运算量

设复序列x(n) 长度为N点,其DFT为

nk X (k ) ? ? x(n)WN n ?0 N ?1

k=0,…,N-1

(1)计算一个X(k) 值的运算量

复数乘法次数: N 复数加法次数: N-14 5.2.1 DFT的运算量

(2)计算全部N个X(k) 值的运算量

复数乘法次数: N2 复数加法次数: N(N-1)

(3)对应的实数运算量

X (k ) ? ? x(n)Wn ?0N ?1 n ?0N ?1nk N

nk nk ? ?[Re x(n) ? j Im x(n)][Re WN ? j Im WN ] n ?0N ?1

nk nk ? ?{[Re x(n) ? ReWN ? Im x(n) ? ImWN ]

nk nk ? j[Re x(n) ? ImWN ? Im x(n) ? ReWN ]}5 一次复数乘法: 4次实数乘法 一个X(k) :

+ 2次实数加法4N次实数乘法 + 2N+2(N-1)= 2(2N-1)次实数加法

所以 整个N点DFT运算共需要:

实数乘法次数: 4 N2

实数加法次数: N×2(2N-1)= 2N(2N-1)6 DFT运算量的结论N点DFT的复数乘法次数举例 N 2 4 8 16 32 N2 4 16 64 256 1028 N 64 128 256 512 1024 N2 4049 16384 65 536 262 144 1 048 576

结论:当N很大时,其运算量很大,对实时性很强的信号 处理来说,要求计算速度快,因此需要改进DFT的计算 方法,以大大减少运算次数。7 5.2.2 减少运算工作量的途径

nk 主要原理是利用系数 WN 的以下特性对DFT进行分解:

(1)对称性

k ( N ?n ) nk ? ? nk (WN ) ? WN ? WN

(2)周期性 (3)可约性 另外,

( n? N ) k n( k ? N ) nk WN ? WN ? WN

mnk nk WmN ? WN

nk nk / m WN ? WN /m

( k ? N / 2) k WN ? ?WNN /2 WN ? ?18 5.3 按时间抽取的基2-FFT算法? ?算法原理 按时间抽取基-2FFT算法与直接计算 DFT运算量的比较 按时间抽取的FFT算法的特点 按时间抽取FFT算法的其它形式流程图? ?9 5.3.1 算法原理

设N=2L,将x(n)按 n 的奇偶分为两组:

x(2r ) ? x1 (r )

x(2r ? 1) ? x2 (r )rN =0,1,…, ? 1 2则

X (k ) ? DFT [ x(n)] ? ? x(n)WNnkn ?0N ?1?n ?0 n为偶数? x(n)WN ?1nk N?n ?0 n为奇数

nk x ( n ) W ? NN ?110 ?2 rk ( 2 r ?1) k ? ? x(2r )WN ? ? x(2r ? 1)WN r ?0 r ?0N ?1 2 r ?0 N ?1 2 r ?0n ?0 n为偶数 N ?1 2? x(n)WN ?1nk N?n ?0 n为奇数 N ?1 2

nk x ( n ) W ? NN ?1rk k rk ? ? x1 (r )W N ? WN x ( r ) W ? 2 N 2 2

k ? X 1 (k ) ? WN X 2 (k )

式中,X1(k)和X2(k)分别是x1(n)和x2(n)的N/2的DFT。

另外,式中k的取值范围是:0,1, …,N/2-1 。11 k X ( k ) ? X ( k ) ? W 因此, 1 N X 2 (k ) 只能计算出X(k)的前一半值。

后一半X(k) 值, N/2 , N/2 +1, …,N ?利用可得到rk r ( N 2?k ) WN ? W 2 N 2N X1 ( ? k ) ? 2X2(N 2?1?r ?0

x1 (r )Wr ( N 2? k ) N 2rk ? ? x1 (r )WN 2 ? X 1 (k ) r ?0N 2 ?1同理可得N ? k ) ? X 2 (k ) 212 考虑到 及前半部分X(k)

( N 2? k ) N 2 k k WN ? WN ?WN ? ?WN

k X (k ) ? X 1 (k ) ? WN X 2 (k )

因此可得后半部分X(k)X (k ?

k=0,1, …,N/2-1N N N k?N 2 ) ? X 1 (k ? ) ? WN X 2 (k ? ) 2 2 2

k ? X1 (k ) ? WN X 2 (k )

k=0,1, …,N/2-113 蝶形运算

k X (k ) ? X 1 (k ) ? WN X 2 (k )

k X (k ) ? X1 (k ) ?WN X 2 (k )

蝶形运算式

蝶形运算信 号流图符号

因此,只要求出2个N/2点的DFT,即X1(k)和X2(k),再 经过蝶形运算就可求出全部X(k)的值,运算量大大减少。14 以8点为例第一次按奇偶分解

以N=8为例, 分解为2个4点 的DFT,然后 做8/2=4次蝶形 运算即可求出 所有8点X(k)的 值。WN01 WNWN23 WN15 蝶形运算量比较?N点DFT的运算量

复数乘法次数: N2 复数加法次数: N(N-1)?

分解一次后所需的运算量=2个N/2的DFT+ N/2蝶形:

复数乘法次数: 2*(N/2)2+N/2=N2/2+N/2 复数加法次数: 2*(N/2)(N/2-1)+2*N/2=N2/2?

因此通过一次分解后,运算工作量减少了差 不多一半。16 进一步按奇偶分解

由于N=2L,因而N/2仍是偶数 ,可以进一步把每个N/2点 子序列再按其奇偶部分分解为两个N/4点的子序列。

以N/2点序列x1(r)为例 则有

X 1 (k ) ?N 2 ?1 r ?0

x1 (2l ) ? x3 (l ) ? N l ? 0,1, ? , ?1 ? x1 (2l ? 1) ? x4 (l ) ? 4rk N 2? x (r )W1?N 4 ?1 l ?0? x (2l )W1N 4 ?1 l ?0 k N 22lk N 2?N 4 ?1 l ?0? x (2l ? 1)W1

( 2l ?1) k N 2?N 4 ?1 l ?0? x3 (l )Wlk N 4?W

lk x ( l ) W ? 4 N4

k ? X 3 (k ) ? WN / 2 X 4 (k )

k=0,1,…,N ?1 417 且?N ? k X1 ? ? k ? ? X 3 (k ) ? WN / 2 X 4 (k ) ?4 ?

k=0,1,…,N ?1 4

由此可见,一个N/2点DFT可分解成两个N/4点DFT。

同理,也可对x2(n)进行同样的分解,求出X2(k)。18 以8点为例第二次按奇偶分解19 算法原理

对此例N=8,最后剩下的是4个N/4= 2点的DFT,2点

DFT也可以由蝶形运算来完成。

以X3(k)为例。

X 3 (k ) ?即N / 4 ?1?l ?0

x3 (l )W

lk N /4?

lk x ( l ) W ? 3 N /4 l ?01

k=0, 10 0 X 3 (0) ? x3 (0) ? W20 x3 (1) ? x(0) ? W2 x(4) ? x(0) ? WN x(4)0 X 3 (1) ? x3 (0) ? W21x3 (1) ? x(0) ? W21 x(4) ? x(0) ? WN x(4)

这说明,N=2M的DFT可全部由蝶形运算来完成。20 以8点为例第三次按奇偶分解N=8按时间抽取法FFT信号流图21 5.3.2 按时间抽取基2-FFT算法与直接计算DFT运算量的比较

由按时间抽取法FFT的信号流图可知,当N=2L时,共有 L 级 蝶形运算;

每级都由 N/2 个蝶形运算组成,而每个蝶形有 1 次复乘、 2 次复加,因此每级运算都需 N/2 次复乘和 N 次复加。22 这样 L 级运算总共需要: 复数乘法:N N ? L ? log 2 N 2 2

复数加法:N ? L ? N log2 N 直接DFT算法运算量 复数乘法: N2 复数加法: N(N-1) 直接计算DFT与FFT算法的计算量之比为M N2 2N M? ? N log2 N log2 N 223 FFT算法与直接DFT算法运算量的比较N N2N log 2 N 2

计算量 之比MNN2N log 2 N 2

计算量 之比M24 8 16 32 64416 64 256 1028 404914 12 32 80 1924.04.0 5.4 8.0 12.8 21.4128256 512 1024 204816 38465 536 262 144 1 048 576 4 194 3044481 024 2 304 5 120 11 26436.664.0 113.8 204.8 372.424 5.3.3 按时间抽取的FFT算法的特点? ? ? ?

序列的逆序排列 同址运算(原位运算) 蝶形运算两节点间的距离r WN的确定25 序列的逆序排列?

序列的逆序排列

由于 x(n) 被反复地按奇、偶分组,所以流图输入端的 排列不再是顺序的,但仍有规律可循: 因为 N=2M , 对于任意 n(0≤n ≤N-1),可以用M个 二进制码表示为:

n( DEC) ? (nM ?1nM ?2 ? n2 n1n0 ) ( BIN )?0 nM ?1 , nM ?2 ,?, n2 , n1 , n0 ? ? ?1

n 反复按奇、偶分解时,即按二进制码的“0” “1” 分解。26 倒位序的树状图(N=8)27 码位的倒位序(N=8)

自然顺序 n 0 1 2 3 二进制数 000 001 010 011 倒位序二进制数 000 100 010 110? 倒位序顺序数 n0 4 2 645 6100101 110001101 01115 37111111728 倒位序的变址处理(N=8)29 同址运算(原位运算)?

同址运算(原位运算)

某一列任何两个节点k 和j 的节点变量进行蝶形运算 后,得到结果为下一列k、j两节点的节点变量,而和其他 节点变量无关。

这种原位运算结构可以节省存储单元, 降低设备成本。 例运算前

A(k ) ? X 1 (k ) X (k )运算后? A(k )A(k ?N ) ? X 2 (k ) 2k WNN ? A(k ? N ) X (k ? ) 2 230 观察原位运算规律31 蝶形运算两节点间的距离?

蝶形运算两节点间的距离

以N=8为例: 第一级蝶形,距离为: 第二级蝶形,距离为: 第三级蝶形,距离为: 1 2 4

规律:对于共L级的蝶形而言,其m级蝶形运算的节

点间的距离为 2 m?132 r WN的确定?W 的确定r N

以N=8为例:r m ? 1时,WN ? WNj / 4 ? W2jm ? W20 , j ? 0r m ? 2时,WN ? WNj / 2 ? W2jm ? W4j , j ? 0,1r m ? 3时,WN ? WNj ? W2jm ? W8j , j ? 0,1,2,3N ? 2M , 第L级:r WN ? W2jL , j ? 0,1,2,?,2 L?1 ? 1? 2 L ? 2 M ? 2 L?M ? N ? 2 L?M?W ? Wr N j N ?2 L ? M?e?j2? N ?2 L ? M?j?e?j2? ? j ?2 M ? L N?W

j ?2 M ? L N33 5.4 按频率抽取的基2-FFT算法?算法原理

先把输入按n的顺序分成前后两半 再把输出X(k)按k的奇偶分组 设序列长度为N=2L,L为整数 前半子序列x(n) 后半子序列 x(n ?N ) 20≤n≤ 0≤n≤N ?1 2 N ?1 234 5.4.1 算法原理

由DFT定义得

nk X (k ) ? ? x(n)WN n ?0N / 2 ?1 n ?0 N ?1N ?1?? x(n)WN / 2?1 n ?0nk N?

n? N / 2?

nk x(n)WN?? x(n)Wn ?0nk N?N / 2?1?n ?0N ( n? 2 ) k x(n ? )WN 2N?N / 2 ?1?

k? ? N N nk 2 ? x(n) ? x(n ? )WN ? ? WN 2 ? ?

k=0,1, …,N35 X (k ) ?N / 2 ?1?n ?0

k? ? N N nk 2 x ( n ) ? x ( n ? ) W ? W ? N ? N 2 ? ?由于 所以W ?eWN k 2 NN 2 N?j2? N ? N 2? e? j? ? ?1? (?1)kN / 2 ?1则

X (k ) ??n ?0N ? nk ? k ? x(n) ? (?1) x(n ? 2 )?WN ? ?

k=0,1, …,N36 然后按k的奇偶可将X(k)分为两部分? k ? 2r ? ? k ? 2r ? 1则式

X (k ) ?N / 2 ?1N r=0,1, …, ? 1 2N ? nk ? k x ( n ) ? ( ? 1 ) x ( n ? )?WN ? 2 ? ??n ?0可转化为

X (2r ) ?N / 2 ?1?n ?0N ? 2 nr N / 2 ?1 ? N ? nr ? ? x ( n ) ? x ( n ? )?WN / 2 x ( n ) ? x ( n ? ) W ? ? ? ? N 2 2 ? n ?0 ? ? ?

[ x ( n) ? x ( n ?N / 2 ?1 N n ( 2 r ?1) N n nr )] ? WN ? ? {[ x(n) ? x(n ? )]WN } ? WN 2 2 2 n ?0

X (2r ? 1) ?N / 2 ?1?n ?037 令N ? x1 (n) ? x(n) ? x(n ? ) ? 2 ? ? n ? x2 (n) ? ? x(n) ? x(n ? N )?WN ? ? 2 ? ? ? ?

X (2r ) ?N / 2 ?1N n=0,1, …, ? 1 2代入?n ?0N ? nr ? x ( n ) ? x ( n ? ) ? WN / 2 ? 2 ? ?N n nr )]WN } ?WN 2 2

X (2r ? 1) ?N / 2 ?1 n ?0? {[ x(n) ? x(n ?1 nr N 2可得? ( 2r ) ?N 2 ?1 n ?0? x (n)WN 2 ?1 n ?0 2?(2r ? 1) ?? x (n)WN r=0,1, …, ? 1 2nr N 2

为2个N/2点的DFT,合起来正好是N点X(k)的值。38 蝶形运算将N ? x ( n ) ? x ( n ) ? x ( n ? ) 1 ? 2 ? ? n ? x2 (n) ? ? x(n) ? x(n ? N )?WN ? ? 2 ? ? ? ?

称为蝶形运算

与时间抽选基2FFT算法中的蝶形运算符号略有不同。39 例 按频率抽取(N=8)

例 按频率抽取,将N点DFT分解为两个N/2点DFT的组合(N=8)40 与时间抽取法的推导过程一样,由于 N=2L,N/2仍然是 一个偶数,因而可以将每个N/2点DFT的输出再分解为偶数组 与奇数组,这就将N/2点DFT进一步分解为两个N/4点DFT。N=841 5.4.2 频率抽取法与时间抽取法的异同?

频率抽取法输入是自然顺序,输出是倒位序 的;

时间抽取法正好相反。

频率抽取法的基本蝶形与时间抽取法的基本 蝶形有所不同。

频率抽取法运算量与时间抽取法相同。???

频率抽取法与时间抽取法的基本蝶形是互为转置的。42 5.5 快速傅里叶逆变换(IFFT)算法

IDFT公式 DFT公式1 x?n ? ? IDFT?X ?k ?? ? NN ?1 n ?0? nk X ( k ) W ? N k ?0N ?1

nk X (k ) ? DFT?x(n)? ? ? x(n)WN

nk 比较可以看出, WN? nk WN

IDFT多出N ?2M1 ?1? ?? ? N ? 2?MM个1/2可分解到M级蝶形运算中。43 例 频率抽取IFFT流图(N=8)44 快速傅里叶逆变换另一种算法1 IDFT [ X (k )] ? N? nk X ( k ) W ? N k ?0 N ?11? 1? ? nk ? ? ? nk ? ? ? ?? X (k )WN ? ? ?? X ( k ) W N ? N ? k ?0 N ? k ?0 ? ?N ?1???N ?1??? 1 ? IFFT [ X (k )] ? ?FFT [ X (k )]? N

X ( k ) ??? ? X ? (k ) ??? ?FFT [ X ? (k )]求共轭求FFT??? ? ?FFT [ X求共轭?

(k )]????? ? x(n)除以N45 5.8 Matlab实现? ?

用FFT进行谱分析的Matlab实现

用CZT进行谱分析的Matlab实现

在Matlab中使用的线性调频z变换函数为czt, 其调用格式为? ??

>>X= czt(x, M, W, A) 其中,x是待变换的时域信号x(n),其长度为N, M是变换的长度,W确定变换的步长,A确定变换 的起点。

若M= N,A= 1,则CZT变成DFT。46 5.8.1 用FFT进行谱分析的Matlab实现

例5.1 设模拟信号 x(t ) ? 2 sin(4? t ) ? 5 cos(8? t ) ,以 t= 0.01n (n=0: N-1) 进行取样,试用fft函数对其做频谱分析。N分别 为:(1) N=45;

(2) N=50;

(3) N=55;

(2) N=60。

程序清单如下 %计算N=45的FFT并绘出其幅频曲线 N=45;n=0:N-1;t=0.01*n; q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); figure(1) subplot(2,2,1) plot(q,abs(y)) title('FFT N=45')47 例5.1程序清单

%计算N=50的FFT并绘出其幅频曲线 N=50;n=0:N-1;t=0.01*n; q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); figure(1) subplot(2,2,2) plot(q,abs(y)) title('FFT N=50')48 %计算N=55的FFT并绘出其幅频曲线 N=55;n=0:N-1;t=0.01*n; q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); figure(1) subplot(2,2,3) plot(q,abs(y)) title('FFT N=55')49 %计算N=60的FFT并绘出其幅频曲线 N=60;n=0:N-1;t=0.01*n; q=n*2*pi/N; x=2*sin(4*pi*t)+5*cos(8*pi*t); y=fft(x,N); figure(1) subplot(2,2,4) plot(q,abs(y)) title('FFT N=60')50 例5.1程序运行结果

FFT N=45 150 150 FFT N=50 100 10050500024 FFT N=55680024 FFT N=60681501501001005050002468002468

从图中可以看出,这几种情况下均有较好的精度。51 例5.1程序运行结果分析

分析:由t=0.01n进行取样可得,采样频率fs=100Hz。

而连续信号的最高模拟角频率为Ω=8 π ,由Ω=2 πf可 得,最高频率为8 π /2 π=4Hz。

因此,满足采样定理的 要求。

采样序列为

x(n) ? 2cos(4? Tn) ? 5cos(8? Tn)即? 4? ? ? 8? ? x(n) ? 2cos ? n ? ? 5cos ? n? ? 100 ? ? 100 ?

为周期序列,周期N=50。

将程序中plot改为stem函数,则可以更清楚地看出频谱。52 例5.1修改程序运行结果

FFT N=45 150 150 FFT N=50 100 10050500024 FFT N=55680024 FFT N=6068150150100100505000246800246853

FFT(快速傅里叶变换)是DFT的一种特殊情况,就是当运算点的个数是2的整数次幂的时候进行的运算(不够用0补齐)。FFT计算原理及流程图:原理:FFT的计算要求点数必须为2的整数次幂,如果点数不够用0补齐。例如计算{2,3,5,8,4}的16点FFT,需要补11个0后进行计算。FFT计算运用蝶形运算,在蝶形运算中变化规律由W(N,p)推导,其中N为FFT计算点数,J为下角标的值。L=1时,W(N,p)=W(N,J)=W(2^L,J),其中J=0;L=2时,W(N,p)=W(N,J)=W(2^L,J),其中J=0,1;L=3时,W(N,p)=W(N,J)=W(2^L,J),其中J=0,1,2,3;所以,W(N,p)=W(2^L,J),其中J=0,1,.,2^(L-1)-1又因为2^L=2^M*2^(L-M)=N*2^(L-M),这里N为2的整数次幂,即N=2^M,W(N,p)=W(2^L,J)=W(N*2^(L-M),J)=W(N,J*2^(M-L))所以,p=J*2^(M-L),此处J=0,1,.,2^(L-1)-1,当J遍历结束但计算点数不够N时,J=J+2^L,后继续遍历,直到计算点数为N时不再循环。流程图:实现代码:123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161/*=*方法名:fft*方法功能:计算数组的FFT,运用蝶形运算*变量名称:*yVector-原始数据*length-原始数据长度*N-FFT计算点数*fftYreal-FFT后的实部*fftYImg-FFT后的虚部*返回值:是否成功的标志,若成功返回true,否则返回false*=*/(BOOL)fft:(floatfloat*)yVector andOriginalLength:(NSInteger)length andFFTCount:(NSInteger)N andFFTReal:(floatfloat*)fftYReal andFFTYImg:(floatfloat*)fftYImg {/确保计算时时2的整数幂点数计算 NSInteger N1=[self nextNumOfPow2:N];定义FFT运算是否成功的标志 BOOL isFFTOK=false;判断计算点数是否为2的整数次幂 if(N!N1){/不是2的整数次幂,直接计算DFT isFFTOK=[self dft:yVector andOriginalLength:length andFFTCount:N andFFTReal:fftYReal andFFTYImg:fftYImg];返回成功标志 return isFFTOK;}/如果计算点数位2的整数次幂,用FFT计算,如下/定义变量 float yVectorN[N1];N点运算的原始数据 NSInteger powOfN=log2(N1);N=2^powOfN,用于标记最大运算级数(公式中表示为:M)NSInteger level=1;运算级数(第几次运算),最大为powOfN,初值为第一级运算(公式中表示为:L)NSInteger lineNum;行号,倒序排列后的蝶形运算行号(公式中表示为:k)float inverseOrderY[N1];yVector倒序x NSInteger distanceLine=1;行间距,第level级运算每个蝶形的两个节点距离为distanceLine=2^(L-1)(公式中表示为:B)NSInteger p;旋转因子的阶数,旋转因子表示为 W(N,p),p=J*2^(M-L)NSInteger J;旋转因子的阶数,旋转因子表示为 W(2^L,J),J=0,1,2,.,2^(L-1)-1=distanceLine-1 float realTemp,imgTemp,twiddleReal,twiddleImg,twiddleTheta,twiddleTemp=PI_x_2/N1;NSInteger N_4=N1/4;判断点数是否够FFT运算点数 if(length){/如果N至少为length,先把yVector全部赋值 for(NSInteger i=0;i;i+){ yVectorN[i]=yVector[i];} if(length){/如果 N>length 后面补零 for(NSInteger i=length;i;i+){ yVectorN[i]=0.0;} } } else {/如果 N截取相应长度的数据进行运算 for(NSInteger i=0;i;i+){ yVectorN[i]=yVector[i];} }/调用倒序方法[self inverseOrder:yVectorN andN:N1 andInverseOrderVector:inverseOrderY];初始值 for(NSInteger i=0;i;i+){ fftYReal[i]=inverseOrderY[i];fftYImg[i]=0.0;}/三层循环/第三层(最里):完成相同旋转因子的蝶形运算/第二层(中间):完成旋转因子的变化,步进为2^level/第一层(最外):完成M次迭代过程,即计算出x(k)=A0(k),A1(k),.,Am(k)=X(k)/第一层循环 while(level){ distanceLine=powf(2,level-1);初始条件 distanceLine=2^(level-1)J=0;NSInteger pow2_Level=distanceLine*2;2^level NSInteger pow2_NSubL=N1/pow2_Level;2^(M-L)/第二层循环 while(J){ p=J*pow2_NSubL;lineNum=J;NSInteger stepCount=0;J运算的步进计数/求旋转因子 if(p=0){ twiddleReal=1.0;twiddleImg=0.0;} else if(p=N_4){ twiddleReal=0.0;twiddleImg=-1.0;} else {/计算尤拉公式中的θ twiddleTheta=twiddleTemp*p;计算复数的实部与虚部 twiddleReal=cos(twiddleTheta);twiddleImg=-11*sin(twiddleTheta);}/第三层循环 while(lineNum){/计算下角标 NSInteger footNum=lineNum+distanceLine;用复数运算计算每级中各行的蝶形运算结e68a84e8a2ade799bee5baa631333363383365果*X(k)=X(k)+X(k+B)*W(N,p)*X(k+B)=X(k)-X(k+B)*W(N,p)*-*/realTemp=fftYReal[footNum]*twiddleReal-fftYImg[footNum]*twiddleImg;imgTemp=fftYReal[footNum]*twiddleImg+fftYImg[footNum]*twiddleReal;将计算后的实部和虚部分别存放在返回数组中 fftYReal[footNum]=fftYReal[lineNum]-realTemp;fftYImg[footNum]=fftYImg[lineNum]-imgTemp;fftYReal[lineNum]=fftYReal[lineNum]+realTemp;fftYImg[lineNum]=fftYImg[lineNum]+imgTemp;stepCount+pow2_Level;行号改变 lineNum=J+stepCount;}/旋转因子的阶数变换,达到旋转因子改变的效果 J+;}/运算级数加一 level+;} isFFTOK=true;return isFFTOK;}内容来自www.egvchb.cn请勿采集。

www.egvchb.cn true http://www.egvchb.cn/wendangku/z0s/f01g/j37c29db94cv/kf7ec4afe04a1b0717fd5360cb20fl.html report 55479 因转码可能存在排版等问题,敬请谅解!以下文字仅供您参考:第五章 快速傅里叶变换 本章目录? ? ? ? ?直接计算DFT的问题及改进的途径 按时间抽取的基2-FFT算法 按频率抽取的基2-FFT算法 快速傅里叶逆变换(IFFT)算法 Matlab实现2 5.1 引言?DFT在实际应用中很重要: 可以计算信号的频谱、功率谱和线性卷积等。?直接按DFT变换进行计算,当序列长度N很大时,计算量非常大,所需时间会很长。 FFT并不是一种与DFT不同的变换
  • 猜你喜欢
马洪刚决战澳门 广西11选五投注网址 内蒙古彩票十一选五 专业股票配资平台有哪些 东方6 1开奖结果查询 中国十只最有价值股票 浙江体彩6十1走势图 融资融券股票一览表 直两码中特 福彩时时彩玩法中奖 上海时时乐走势图彩经网络 正规微信一元投资赚钱 湖北30选5开奖结果规定 今天股票行情查询 秒速时时彩走势图官方网 极速赛车app官方版 四川体育彩票网