快速傅里叶变换原理介绍及MATLAB实现

  • Post author:
  • Post category:其他


FFT算法主要有两种,按时间抽选的FFT的算法(DIT-FFT)和按频率抽选的FFT算法(DIF-FFT)。



时间抽取基-2 FFT算法





N

=

2

M

N={

{2}^{M}}






N




=











2












M














,有限长序列



x

(

n

)

x\left( n \right)






x





(


n


)






的离散傅里叶变换(DFT)为





X

(

k

)

=

n

=

0

N

1

x

(

n

)

e

j

n

k

2

π

N

,

0

k

N

1

X(k)=\sum\limits_{n=0}^{N-1}{x\left( n \right){

{e}^{-jnk\frac{2\pi }{N}}}},0\le k\le N-1






X


(


k


)




=

















n


=


0



















N





1





















x





(


n


)








e















j


n


k














N
















2


π































,




0













k













N













1







直接使用DFT,计算一个



X

(

k

)

X\left( k \right)






X





(


k


)






值,需计算:N次复数相乘和



(

N

1

)

\left( N-1 \right)







(


N









1


)






次复数相加。计算N点



X

(

k

)

X\left( k \right)






X





(


k


)






值,需计算:



N

2

{

{N}^{2}}









N












2














次复数相乘和



N

(

N

1

)

N\left( N-1 \right)






N





(


N









1


)






次复数相加。可以看出,在DFT计算中,不论是乘法还是加法,运算量均与



N

2

{

{N}^{2}}









N












2














成正比。这给处理器的运算性能带来了极大的挑战。令





W

N

=

e

j

2

π

N

{

{W}_{N}}={

{e}^{-j\frac{2\pi }{N}}}









W












N






















=











e















j














N
















2


π



































那么





X

(

k

)

=

n

=

0

N

1

x

(

n

)

W

N

n

k

,

0

k

N

1

X(k)=\sum\limits_{n=0}^{N-1}{x\left( n \right)W_{N}^{nk}},0\le k\le N-1






X


(


k


)




=

















n


=


0



















N





1





















x





(


n


)






W











N










n


k




















,




0













k













N













1







将序列



x

(

n

)

x\left( n \right)






x





(


n


)






按下标偶数项和奇数项分解为两组,偶数项为一组,奇数项为一组,得到两个



N

/

2

N/2






N


/


2





点的子序列,即





{

x

1

(

r

)

=

x

(

2

r

)

x

2

(

r

)

=

x

(

2

r

+

1

)

r

=

0

,

1

,


,

N

/

2

1

\left\{ \begin{aligned} & {

{x}_{1}}\left( r \right)=x\left( 2r \right) \\ & {

{x}_{2}}\left( r \right)=x\left( 2r+1 \right) \\ \end{aligned} \right.r=0,1,\cdots ,N/2-1








{
















































x












1























(


r


)





=




x





(


2


r


)














x












2























(


r


)





=




x





(


2


r




+




1


)

























r




=








0


,




1


,











,




N


/


2













1













X

(

k

)

=

r

=

0

N

/

2

1

x

(

2

r

)

W

N

2

r

k

+

r

=

0

N

/

2

1

x

(

2

r

+

1

)

W

N

(

2

r

+

1

)

k

X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{x\left( 2r \right)}W_{N}^{2rk}+\sum\limits_{r=0}^{N/2-1}{x\left( 2r+1 \right)}W_{N}^{(2r+1)k}






X





(


k


)





=

















r


=


0



















N


/


2





1





















x





(


2


r


)





W











N










2


r


k





















+

















r


=


0



















N


/


2





1





















x





(


2


r




+




1


)





W











N










(


2


r


+


1


)


k
























又因为





W

N

2

r

k

=

e

j

2

π

N

2

r

k

=

e

j

2

π

N

/

2

r

k

=

W

N

/

2

r

k

W_{N}^{2rk}={

{e}^{-j\frac{2\pi }{N}2rk}}={

{e}^{-j\frac{2\pi }{N/2}rk}}=W_{N/2}^{rk}







W











N










2


r


k





















=











e















j














N
















2


π





















2


r


k













=











e















j














N


/


2
















2


π





















r


k













=









W











N


/


2










r


k
























所以





X

(

k

)

=

r

=

0

N

/

2

1

x

(

2

r

)

W

N

/

2

r

k

+

r

=

0

N

/

2

1

x

(

2

r

+

1

)

W

N

k

W

N

/

2

r

k

X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{x\left( 2r \right)}W_{N/2}^{rk}+\sum\limits_{r=0}^{N/2-1}{x(2r+1)W_{N}^{k}W_{N/2}^{rk}}






X





(


k


)





=

















r


=


0



















N


/


2





1





















x





(


2


r


)





W











N


/


2










r


k





















+

















r


=


0



















N


/


2





1





















x


(


2


r




+




1


)



W











N










k




















W











N


/


2










r


k

























上式又可以写作





X

(

k

)

=

r

=

0

N

/

2

1

x

1

(

r

)

W

N

/

2

r

k

+

W

N

k

r

=

0

N

/

2

1

x

2

(

r

)

W

N

/

2

r

k

=

X

1

(

k

)

+

W

N

k

X

2

(

k

)

\begin{aligned} & X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{

{

{x}_{1}}}\left( r \right)W_{N/2}^{rk}+W_{N}^{k}\sum\limits_{r=0}^{N/2-1}{

{

{x}_{2}}}\left( r \right)W_{N/2}^{rk} \\ & \begin{matrix} \begin{matrix} {} & {} \\ \end{matrix} & ={

{X}_{1}} \\ \end{matrix}\left( k \right)+W_{N}^{k}{

{X}_{2}}\left( k \right) \\ \end{aligned}















































X





(


k


)





=













r


=


0



















N


/


2





1
























x












1
























(


r


)






W











N


/


2










r


k





















+





W











N










k






























r


=


0



















N


/


2





1
























x












2
























(


r


)






W











N


/


2










r


k


























































































































=







X












1








































(


k


)





+





W











N










k






















X












2























(


k


)

























同理可得





X

(

k

+

N

2

)

=

X

1

(

k

)

W

N

k

X

2

(

k

)

k

=

0

,

1

,

.

.

.

,

N

/

2

1

X\left( k+\frac{N}{2} \right)={

{X}_{1}}\left( k \right)-W_{N}^{k}{

{X}_{2}}\left( k \right)\begin{matrix} {} & k=0,1,…,N/2-1 \\ \end{matrix}






X






(



k




+















2














N





















)






=











X












1























(


k


)















W











N










k






















X












2























(


k


)












































k




=




0


,




1


,




.


.


.


,




N


/


2









1
























那么,傅里叶变换可以写作





{

X

(

k

)

=

X

1

(

k

)

+

W

N

k

X

2

(

k

)

X

(

k

+

N

2

)

=

X

1

(

k

)

W

N

k

X

2

(

k

)

k

=

0

,

1

,

.

.

.

,

N

/

2

1

\left\{ \begin{aligned} & X\left( k \right)={

{X}_{1}}\left( k \right)+W_{N}^{k}{

{X}_{2}}\left( k \right) \\ & X\left( k+\frac{N}{2} \right)={

{X}_{1}}\left( k \right)-W_{N}^{k}{

{X}_{2}}\left( k \right) \\ \end{aligned} \right.\begin{matrix} {} & k=0,1,…,N/2-1 \\ \end{matrix}
















































































































X





(


k


)





=







X












1























(


k


)





+





W











N










k






















X












2























(


k


)











X






(



k




+















2














N





















)






=







X












1























(


k


)











W











N










k






















X












2























(


k


)
































































k




=




0


,




1


,




.


.


.


,




N


/


2









1
























通过M次分解最后分解成2点的DFT运算。则,时间抽取基-2 FFT总共需要的运算量为:复乘:



N

/

2

M

=

N

/

2

log

2

N

N/2\cdot M=N/2\cdot {

{\log }_{2}}N






N


/


2













M




=








N


/


2
















lo

g













2




















N





,复加:



N

M

=

N

log

2

N

N\cdot M=N{

{\log }_{2}}N






N













M




=








N





lo

g













2




















N





。显然,要比直接DFT运算快得多。



频率抽取基-2 FFT算法

DIF-FFT算法是将输入序列



X

(

k

)

X\left( k \right)






X





(


k


)






分成前后两个部分





X

(

k

)

=

n

=

=

0

N

1

x

(

n

)

W

N

n

k

=

n

=

0

N

2

1

x

(

n

)

W

N

n

k

+

n

=

N

2

N

1

x

(

n

)

W

N

n

k

=

n

=

0

N

1

x

(

n

)

W

N

n

k

+

n

=

0

N

1

x

(

n

+

N

2

)

W

N

k

(

n

+

N

2

)

=

n

=

0

N

2

1

[

x

(

n

)

+

x

[

n

+

N

2

]

W

N

N

k

/

2

]

W

N

n

k

\begin{aligned} & X\left( k \right)=\sum\limits_{n==0}^{N-1}{x\left( n \right)W_{N}^{nk}}=\sum\limits_{n=0}^{\frac{N}{2}-1}{x\left( n \right)W_{N}^{nk}}+\sum\limits_{n=\frac{N}{2}}^{N-1}{x\left( n \right)W_{N}^{nk}} \\ & =\sum\limits_{n=0}^{N-1}{x\left( n \right)W_{N}^{nk}}+\sum\limits_{n=0}^{N-1}{x\left( n+\frac{N}{2} \right)W_{N}^{k\left( n+\frac{N}{2} \right)}} \\ & =\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+x\left[ n+\frac{N}{2} \right]W_{N}^{Nk/2} \right]W_{N}^{nk}} \end{aligned}





















































X





(


k


)





=













n


=


=


0



















N





1





















x





(


n


)






W











N










n


k






















=













n


=


0































2
















N
























1





















x





(


n


)






W











N










n


k






















+













n


=














2
















N






































N





1





















x





(


n


)






W











N










n


k






























=













n


=


0



















N





1





















x





(


n


)






W











N










n


k






















+













n


=


0



















N





1





















x






(



n




+















2














N





















)







W











N










k




(



n


+














2
















N






















)
































=













n


=


0































2
















N
























1























[



x





(


n


)





+




x






[



n




+















2














N





















]







W











N










N


k


/


2




















]







W











N










n


k








































因为



W

N

N

/

2

=

1

W_{N}^{N/2}=-1







W











N










N


/


2





















=











1





,则



W

N

N

k

/

2

=

(

1

)

k

=

{

1

k

1

k

W_{N}^{Nk/2}={

{\left( -1 \right)}^{k}}=\left\{ \begin{matrix} 1 & k为偶数 \\ -1 & k为奇数 \\ \end{matrix} \right.







W











N










N


k


/


2





















=












(





1


)













k













=










{














1











1





























k

















k


































,所以可以得到





X

(

k

)

=

n

=

0

N

2

1

[

x

(

n

)

+

(

1

)

k

x

(

n

+

N

2

)

]

W

N

n

k

X(k)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+{

{\left( -1 \right)}^{k}}x\left( n+\frac{N}{2} \right) \right]W_{N}^{nk}}






X


(


k


)




=

















n


=


0































2
















N
























1























[



x





(


n


)





+








(





1


)













k











x






(



n




+















2














N





















)





]







W











N










n


k























将k按照奇数和偶数来分,得到





{

k

=

2

r

k

k

=

2

r

+

1

k

r

=0,1,

N/2-1

\left\{ \begin{matrix} k=2r & k为偶数 \\ k=2r+1 & k为奇数 \\ \end{matrix} \right.r\text{=0,1,}\cdots \text{N/2-1}








{














k




=




2


r








k




=




2


r




+




1





























k

















k

































r



=0,1,











N/2-1








将X(k)分为两部分





{

X

(

2

r

+

1

)

=

n

=

0

N

2

1

[

x

(

n

)

x

(

n

+

N

2

)

]

W

N

n

W

N

/

2

r

n

X

(

2

r

)

=

n

=

0

N

2

1

[

x

(

n

)

+

x

(

n

+

N

2

)

]

W

N

/

2

r

n

\left\{ \begin{aligned} & X\left( 2r+1 \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)-x\left( n+\frac{N}{2} \right) \right]W_{N}^{n}W_{N/2}^{rn}} \\ & X\left( 2r \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+x\left( n+\frac{N}{2} \right) \right]W_{N/2}^{rn}} \\ \end{aligned} \right.










































































































































































































X





(


2


r




+




1


)





=













n


=


0































2
















N
























1























[



x





(


n


)










x






(



n




+















2














N





















)





]







W











N










n




















W











N


/


2










r


n




























X





(


2


r


)





=













n


=


0































2
















N
























1























[



x





(


n


)





+




x






(



n




+















2














N





















)





]







W











N


/


2










r


n

















































x

1

=

x

(

n

)

+

x

(

n

+

N

2

)

{

{x}_{1}}=x\left( n \right)+x\left( n+\frac{N}{2} \right)









x












1






















=








x





(


n


)





+








x






(



n




+
















2
















N






















)











x

2

(

n

)

=

[

x

(

n

)

x

(

n

+

N

2

)

]

W

N

n

{

{x}_{2}}\left( n \right)=\left[ x\left( n \right)-x\left( n+\frac{N}{2} \right) \right]W_{N}^{n}









x












2























(


n


)





=










[



x





(


n


)










x






(



n




+
















2
















N






















)





]







W











N










n






















,可以得到





{

X

(

2

r

+

1

)

=

n

=

0

N

2

1

x

2

(

n

)

W

N

/

2

r

n

X

(

2

r

)

=

n

=

0

N

2

1

x

1

(

n

)

W

N

/

2

r

n

\left\{ \begin{aligned} & X\left( 2r+1 \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{

{

{x}_{2}}\left( n \right)W_{N/2}^{rn}} \\ & X\left( 2r \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{

{

{x}_{1}}\left( n \right)W_{N/2}^{rn}} \\ \end{aligned} \right.










































































































































































































X





(


2


r




+




1


)





=













n


=


0































2
















N
























1
























x












2























(


n


)






W











N


/


2










r


n




























X





(


2


r


)





=













n


=


0































2
















N
























1
























x












1























(


n


)






W











N


/


2










r


n













































则经过M次分解,同样可以将DFT分解成为N/2个两点DFT。以8点DIF-FFT运算为例,DIF-FFT运算流图如下图所示

在这里插入图片描述

从上图可以看出,输出序列的排列规律不是从小到大按顺序的,而是按照倒叙规则排序的,即先将0-7转换为二进制数,然后将二进制数左右倒序,再转为十进制就可以得到新的数列,即:0,4,2,6,1,5,3,7。



利用MATLAB实现DFT

DFT的矩阵表示形式为





[

X

(

0

)

X

(

1

)

X

(

2

)

X

(

N

1

)

]

=

[

W

N

0

W

N

0

W

N

0

W

N

0

W

N

0

W

N

1

W

N

2

W

N

N

1

W

N

0

W

N

2

W

N

4

W

N

2

(

N

1

)

W

N

0

W

N

N

1

W

N

2

(

N

1

)

W

N

(

N

1

)

(

N

1

)

]

[

x

(

0

)

x

(

1

)

x

(

2

)

x

(

N

1

)

]

\left[ \begin{matrix} X\left( 0 \right) \\ X\left( 1 \right) \\ X\left( 2 \right) \\ \begin{matrix} \cdots \\ X\left( N-1 \right) \\ \end{matrix} \\ \end{matrix} \right]=\left[ \begin{matrix} W_{N}^{0} & W_{N}^{0} & W_{N}^{0} & \begin{matrix} \cdots & W_{N}^{0} \\ \end{matrix} \\ W_{N}^{0} & W_{N}^{1} & W_{N}^{2} & \begin{matrix} \cdots & W_{N}^{N-1} \\ \end{matrix} \\ W_{N}^{0} & W_{N}^{2} & W_{N}^{4} & \begin{matrix} \cdots & W_{N}^{2\left( N-1 \right)} \\ \end{matrix} \\ \begin{matrix} \cdots \\ W_{N}^{0} \\ \end{matrix} & \begin{matrix} \cdots \\ W_{N}^{N-1} \\ \end{matrix} & \begin{matrix} \cdots \\ W_{N}^{2\left( N-1 \right)} \\ \end{matrix} & \begin{matrix} \begin{matrix} \cdots & \cdots \\ \end{matrix} \\ \begin{matrix} \cdots & W_{N}^{\left( N-1 \right)\left( N-1 \right)} \\ \end{matrix} \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} x\left( 0 \right) \\ x\left( 1 \right) \\ x\left( 2 \right) \\ \begin{matrix} \cdots \\ x\left( N-1 \right) \\ \end{matrix} \\ \end{matrix} \right]


























































































X





(


0


)









X





(


1


)









X





(


2


)




























X





(


N









1


)

















































































































=





























































































W











N










0


























W











N










0


























W











N










0













































W











N










0
































































W











N










0


























W











N










1


























W











N










2













































W











N










N





1
































































W











N










0


























W











N










2


























W











N










4













































W











N










2




(



N





1



)










































































































W











N










0



















































































W











N










N





1



















































































W











N










2




(



N





1



)

































































































































































W











N












(



N





1



)






(



N





1



)








































































































































































































































x





(


0


)









x





(


1


)









x





(


2


)




























x





(


N









1


)




















































































































根据DFT公式和矩阵展开形式,利用MATLAB可以得到如下结果

在这里插入图片描述



利用MATLAB实现FFT(8点DIF-FFT为例)

由DIF-FFT的流程图可以得到,FFT运算的基本单元是蝶形运算单元。蝶形运算单元可以用简单的语句实现,然后可以用循环语句来进行各个蝶形运算单元的运算。

8点FFT的蝶形运算有3级,第1级有1组,每组4个蝶形运算单元,旋转因子是



W

8

0

W_{8}^{0}







W











8










0


























W

8

1

W_{8}^{1}







W











8










1


























W

8

2

W_{8}^{2}







W











8










2


























W

8

3

W_{8}^{3}







W











8










3






















;第2级有2组,每组2个蝶形运算单元,循环因子是



W

4

0

W_{4}^{0}







W











4










0


























W

4

1

W_{4}^{1}







W











4










1






















;第3级有4组,每组1个蝶形运算单元,旋转因子是



W

2

0

W_{2}^{0}







W











2










0






















。总结运算规律,来设定循环语句。

第一层循环在1到3级间循环,循环变量



m

m

=

1

,

2

,

3

mm=1,2,3






m


m




=








1


,




2


,




3





。旋转因子下标



N

m

=

2

4

m

m

=

8

,

4

,

2

N_{m}=2^{4-mm}=8,4,2







N











m





















=









2











4





m


m












=








8


,




4


,




2







第二层循环在该级的各组间循环,每级有



2

m

m

1

2^{mm-1}







2











m


m





1













组,每组第一行对应的



x

x






x





值为:第1级是



x

(

0

)

x(0)






x


(


0


)





,第2级是



x

(

0

)

x(0)






x


(


0


)









x

(

4

)

x(4)






x


(


4


)





,第3级是



x

(

0

)

x(0)






x


(


0


)









x

(

2

)

x(2)






x


(


2


)









x

(

4

)

x(4)






x


(


4


)









x

(

6

)

x(6)






x


(


6


)





。设第二层循环变量为



p

p






p





,则在MATLAB中,



p

=

0

:

N

m

:

7

p=0:N_{m}:7






p




=








0




:









N











m





















:








7







第三层循环在该组的各个蝶形运算单元间循环,每组有



N

m

2

\frac{

{

{N}_{m}}}{2}


















2



















N












m










































个蝶形运算单元,则循环变量



k

k






k





从1到



N

m

2

\frac{

{

{N}_{m}}}{2}


















2



















N












m










































,旋转因子是



e

2

π

(

k

1

)

j

N

m

{

{e}^{\frac{-2\pi (k-1)j}{Nm}}}









e
























N


m



















2


π


(


k





1


)


j

































,每次蝶形运算跨越的行数是



N

m

2

\frac{

{

{N}_{m}}}{2}


















2



















N












m










































,则参加蝶形运算的



x

x






x





值为



x

(

k

+

p

)

x(k+p)






x


(


k




+








p


)









x

(

k

+

p

+

N

m

2

)

x(k+p+\frac{

{

{N}_{m}}}{2})






x


(


k




+








p




+




















2



















N












m







































)







循环完成后则FFT运算完成,再将x序列按倒叙规律重新排列就可以得到



X

(

k

)

X(k)






X


(


k


)





序列。



涉及到的MATLAB函数



exp(x)计算e的x次方,abs(x)计算x绝对值。



倒序运算主要用到的函数有,dec2bin(x,m),是把十进制序列x转换为m位二进制数,bin2dec(x,m)也是类似功能。



fliplr函数是将一个矩阵左右颠倒,则程序中求倒序的语句就是:d=bin2dec(fliplr(dec2bin([0:N-1],m)))+1;y=x(d); 其中,x是N点序列,执行完语句后,y序列就是x序列的倒序排列。

根据上述流程,可以得到结果如下

在这里插入图片描述



调用MATLAB的FFT函数

MATLAB自带的FFT函数语法格式为:A=fft(X,N)。A返回信号X的FFT变换矩阵,输入信号X和输出信号A大小相同。N为FFT的点数,当X长度不足N,自动在数据末尾补零,当X长度超过N,则进行自动截取。结果如下
在这里插入图片描述

通过对比可以发现,和上述实现的FFT结果是一致的。

实验一代码如下:

clear;
close all;
clc;
N=256;                      %信号点数
n=0:N-1;                    %时间采样
xn=cos(pi*n/6);           %信号
k=0:N-1;                    %频域采样
WN=exp(-1j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;            %计算DFT
figure(1);
subplot(2,1,1);stem(xn);title('变换序列');xlim([0 50])
subplot(2,1,2);stem(abs(Xk));title('序列的DFT');

实验二代码如下:

clear;
close all;
clc;
N=8;                                                           %FFT点数
n=0:N-1;
x=[0 1 2 2 1 0 1 2];                                      %变换序列
figure(1);
subplot(2,1,1);stem(x);title('变换序列');
x1=x;
m=log2(N);                                                 %蝶形运算的级数
for mm=1:m                                               %循环蝶形运算
    Nm=2^(m-mm+1);                                 %旋转因子下标Nm=842
    for p=0:Nm:N-1                                       %循环该级1~2mm-1组蝶形运算
        for k=1:Nm/2
            kp=k+Nm/2+p;                               %蝶形运算的下标
            a=x(kp);
            x(kp)=(x(k+p)-a)*exp(-1j*2*pi*(k-1)/Nm);
            x(k+p)=x(k+p)+a;
        end
    end
end
d=bin2dec(fliplr(dec2bin([0:N-1],m)))+1;               %倒序排列
y=x(d);
mag=abs(y);
subplot(2,1,2);stem(mag);title('FFT结果');
x=[0 1 2 2 1 0 1 2];
yy=fft(x,N);
figure(2);
subplot(2,1,1);stem(x);title('变换序列');
subplot(2,1,2);stem(abs(yy));title('FFT结果');



版权声明:本文为yhcwjh原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。