常用FFT算法
常用FFT 算法总结
一、DFT 及IDFT 的定义
对于N 点有限长序列)(n x ，其DFT 及IDFT 定义如下： DFT ：∑-==1
0)()(N n nk N W n x k X
IDFT ：∑-=-=1
)(1)(N k nk N
W
k X N
n x
其中，N
p j p N
e
W
/2π-=称为旋转因子，p 称为旋转因子的指数。
二、基2-FFT 算法 设序列)(n x 的长度N 满足M
N
2
=，M 为自然数。
1、时间域抽取FFT （DIT-FFT ） （1）算法原理
按n 的奇偶把)(n x 分解为两个N/2点的子序列：
)
12()()2()(21+==r x r x r x r x
则)(n x 的DFT 为
)()()12()2()(211
2/0
)12(1
2/0
2k X W k X W
r x W
r x k X k
N N r k r N
N r rk N
+=++
=
∑∑-=+-=
其中)(1k X 和)(2k X 分别为)(1r x 和)(2r x 的N/2点DFT 。 利用)(1k X 和)(2k X 的周期性，)(k X 可以表示为
??
?
??-=++=)()()2()()()(2121k X W k X N k X k X W k X k X k
N k N 上式表明一个N 点的DFT 可以用两个N/2点的DFT 来表示。 （2）运算量 当M
N
2
=时，共有M 级蝶形，每级N/2个蝶形，每个蝶形有1次复数乘法2次复数加法。
复数乘法：N
N M N m
F
2
log
2
2
=
=
复数加法：N
N NM a F
2
log
==
N
N N
N N
FFT m DFT m F F 2
2
2
log
2log 2
)
()(=
=
（3）蝶形运算
?????-=+=----)
()()()
()()(1111j X W k X j X j X W k X k X m p
N m m m p
N m m m
表示第m 级迭代，k 和j 表示数据所在的行数，12-+=m k j ，m M k p -?=2。
2、频率域抽取FFT （DIF-FFT ） （1）算法原理
将)(n x 前后对半分开，得到两个子序列，其DFT 可表示为如下形式
nk N
N n kN N N n k N n N
N n nk N
N N n nk N
N n nk
N
W N n x W n x W
N n x W n x W
n x W n x k X ∑
∑
∑
∑
∑
-=-=+-=-=-=??
????++=
+
+
=
+
=
1
2/0
2
/1
2/0
)2/(1
2/0
1
2
/1
2/0
)2()()2
()()()()(
将)(k X 分解成偶数组与奇数组，可以得到
rn
N n N N n rn N N n W W N n x n x r X W N n x n x r X 2/1
2/0
2/1
2/0
)2()()12()2()()2(∑
∑
-=-=??
???
?
+-=
+??
????
++=
令
???
???
?
?????
?
+-=+
+=n N W N n x n x n x N n x n x n x )2()()()
2
()()(21
则)2(r X 和)12(+r X 分别是)(1n x 和)(2n x 的N/2点DFT ，这表明一个N 点的DFT 可以分解成两个N/2点的DFT 。 （2）运算量
运算量与时间域抽取FFT 相同。 （3）蝶形运算
[]???-=+=----p N
m m m m m m W j X k X j X j X k X k X )()()()
()()(1111
m
表示第m 级迭代，k 和j 表示数据所在的行数，m M k j -+=2，1
2
-?=m k p 。
三、基4-FFT 算法 设序列)(n x 的长度N 满足M
N
4
=，M 为自然数。
1、时间域抽取FFT （DIT-FFT ） （1）算法原理
将)(n x 分解为四个N/4点的子序列：
)
34()()24()()14()()4()(4321+=+=+==r x r x r x r x r x r x r x r x
则)(n x 的DFT 为
)
()()()()34()24()14()4()(4332211
4/0
)34(1
4/0
)24(1
4/0
)14(1
4/0
4k X W k X W k X W k X W
r x W
r x W
r x W
r x k X k
N k N k N N r k r N
N r k r N
N r k r N
N r rk N
+++=++
++
++
=
∑∑∑∑-=+-=+-=+-=
其中)(1k X 、)(2k X 、)(3k X 和)(4k X 分别为)(1r x 、)(2r x 、)(3r x 和)(4r x 的N/4点DFT 。 利用)(1k X 、)(2k X 、)(3k X 和)(4k X 的周期性，)(k X 可以分成4段：
?????
??????--+=+
-+-=++--=++++=)
()()()()4
3()
()()()()2()()()()()4()()()()()(433221433221433221433221k X jW k X W k X jW k X N
k X k X W
k X W
k X W k X N k X k X jW k X W k X jW k X N k X k X W k X W k X W k X k X k
N k
N k N k N
k N
k
N
k
N k
N k
N k
N k
N k
N
上式表明一个N 点的DFT 可以用四个N/4点的DFT 来表示。
2、频率域抽取FFT （DIF-FFT ）
将)(n x 按顺序平均分成四组，得到四个子序列，其DFT 可表示为如下形式
nk
N N n kN N kN N kN N N N n nk N
N N n nk N
N N n nk N
N n nk N
W N n x W N n x W N n x W n x W
n x W n x W n x W n x k X ∑
∑
∑
∑
∑
-=-=-=-=-=??
????++++++=
+
+
+
=
1
4/0
4
/32/4/1
4
/31
4/32/1
2/4/1
4/0
)43()2()4()()()()()()(
按r k
4=、14+=r k 、2
4+=r k 和34+=r k 展开上式，可以得到
??????
?????????
???
?
+-+-++=
+??
???
?+-+++-=+??
???
?
+++-+-=+??
????
++++++=
∑
∑
∑
∑
-=-=-=-=n N nr N N n n
N nr N N n n N nr N N n nr
N N n W W N n jx N n x N n jx n x r X W W N n x N n x N n x n x r X W W N n jx N n x N n jx n x r X W N n x N n x N n x n x r X 34/1
4/0
24/1
4/04/1
4/0
4/1
4/0
)43()2()4()()34()43()2()4()()24()43()2()4()()14()43()2()4()()4(
令
???
???
???
???
?
??
???
?
+-+-+==?????
?+-+++-=??
???
?
+++-+-=+
++
++
+=n N n
N n
N W N n jx N n x N n jx n x n x W N n x N n x N n x n x n x W N n jx N n x N n jx n x n x N n x N n x N n x n x n x 342321)43()2()4()()()43()2()4()()()43()2()4()()()
4
3()2
()4
()()(
则)4(r X 、)14(+r X 、)24(+r X 和)34(+r X 分别是)(1n x 、)(2n x 、)(3n x 和)(4n x 的N/4点DFT ，这表明一个N 点的DFT 可以分解成四个N/4点的DFT 。 四、分裂基算法
对偶序号使用基-2FFT 算法，对奇序号使用基-4FFT 算法，称分裂基FFT 算法 设序列)(n x 的长度N 满足M
N
2
=，M 为自然数。
1、时间域抽取FFT （DIT-FFT ） 将)(n x 分为三个序列：
)
34()()14()()2()(321+=+==l x l x l x l x r x r x ，则
)()()()34()14()2()()(33211
4/0
)34(1
4/0
)14(1
2/0
210
k X W k X W k X W
l x W
l x W
r x W
n x k X k
N k
N N l k l N
N l k l N
N r rk N
N n nk N
++=++
++
=
=
∑∑∑∑-=+-=+-=-=
其中)(1k X 为)(1r x 的N/2点DFT ，)(2k X 和)(3k X 为)(2l x 和)(3l x 的N/4点DFT 。 利用)(1k X 、)(2k X 和)(3k X 的周期性，)(k X 可分成4段：
?????
??????-++
=+
--=++-+
=+++=)
()()4
()4
3()
()()()2()
()()4
()4()()()()(3321332133213321k X jW k X jW N k X N
k X k X W k X W k X N k X k X jW k X jW N k X N k X k X W k X W k X k X k
N k
N k N
k
N
k
N k
N k
N k
N
上式表明一个N 点的DFT 可以用一个N/2点的DFT 和两个N/4点的DFT 来表示。 2、频率域抽取FFT （DIF-FFT ）
)(n x 按基
2和基4抽取，其DFT 可表示为如下形式
nk
N N n kN N W N n x W n x k X ∑
-=??
???
?++=
1
2/02
/)2()()(
nk N N n kN N kN N kN N W N n x W N n x W N n x W n x k X ∑
-=??
???
?++++++=
1
4/0
4
/32/4/)43()2()4()()(
按r k 2=和14+=r k 及34+=r k 分别展开上两式，可以得到
???
?
??
?
????
??
??
+
-+-++=+??
??
??+++-+-=+??????++=∑∑∑-=-=-=n
N nr N N n n
N nr N N n rn N N n W W N n jx N n x N n jx n x r X W W N n jx N n x N n jx n x r X W N n x n x r X 34/14/04/1
4/02/1
2/0)43()2
()4()()34()43()2()4()()14()2()()2(
令
??
??
?
????
??
???
?+-+-+==?????
?
+++-+-=+
+=n N n N W N n jx N n x N n jx n x n x W N n jx N n x N n jx n x n x N n x n x n x 3321)43()2()4()()()43()2()4()()()
2()()(
则)2(r X 是)(1n x 的N/2点DFT ，)14(+r X 和)34(+r X 分别是)(2n x 和)(3n x 的N/4点DFT ，这表明一个N 点的DFT 可以分解成1个N/2点的DFT 和2个N/4点的DFT 。
FFT 算法详解
2.3 快速傅立叶变换问题 1) 问题背景 在数值电路的传输中，为了避免信号干扰，需要把一个连续信号 x(t)先通过取样离散化为一列数值脉冲信号x(0), x(1), …… ,然后再通过编码送到传输电路中。如果取样间隔很小，而连续信号的时间段又很长，则所得到的数值脉冲序列将非常庞大。因此，传输这个编码信号就需要长时间的占用传输电路，相应地也需要付出昂贵的电路费用。 那么能否经过适当处理是使上述的数值脉冲序列变短，而同时又不会丧失有用的信息？的经过研究，人们发现，如果对上述数值脉冲序列作如下的变换处理： ∑-=--=-== 1 /21 ,1,...,1,0,)()(N k N nki i N n e k x n X π (1) 则所得到的新序列X(0), X(1) , ……将非常有序，其值比较大的点往往集中在某一很狭窄的序列段内，这将非常有利于编码和存储，从而达到压缩信息的目的。 公式（1）就是所谓的离散傅立叶变换，简称DFT 。现在我们来分析一下计算DFT 所需要的工作量。如果我们不考虑公式（7.1）中指数项的运算，那么计算其每一个点X (n) 需要N 次复数乘法和N-1次的复数加法。显然当N 很大时，这个工作量也非常巨大。正是由于这个原因，使得DFT 的应用范围在过去很长的时间里受到了严格的限制。注意到公式（1）是非常有规律性的，那么能否利用这种规律性来降低DFT 的计算时间？ 1965年，凯莱和塔柯的提出了一种用于计算DFT 的数学方法，大大减少了DFT 的计算时间，同时又特别适用于硬件处理，这就是所谓的快速傅里叶变换，简称FFT 。鉴于DFT 的数据结构可以通过傅立叶变换的离散化获得，亦可通过三角插值得到，而本质上又同连续傅里叶分析有着极为密切的关系。下面我们从傅立叶级数级数和傅立叶积分入手，导出DFT 结构的来源和FFT 的工作原理。 2） 傅立叶变换 如果x(t)是定义在整个实轴上的实值或复值函数，则其傅立叶变换可由下式给出： ? ∞ ∞ ---== 1 ,)()(/2i dt e t x f X T nift （2）
FFT算法的意义
FFT是离散傅立叶变换的快速算法，可以将一个信号变换 到频域。有些信号在时域上是很难看出什么特征的，但是如 果变换到频域之后，就很容易看出特征了。这就是很多信号 分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱 提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去 做，但是却不知道FFT之后的结果是什意思、如何决定要使用 多少点来做FFT。 现在圈圈就根据实际经验来说说FFT结果的具体物理意义。 一个模拟信号，经过ADC采样之后，就变成了数字信号。采样 定理告诉我们，采样频率要大于信号频率的两倍，这些我就 不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点， 经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT 运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT 之后结果就是一个为N点的复数。每一个点就对应着一个频率 点。这个点的模值，就是该频率值下的幅度特性。具体跟原始 信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT 的结果的每个点（除了第一个点直流分量之外）的模值就是A 的N/2倍。而第一个点就是直流分量，它的模值就是直流分量 的N倍。而每个点的相位呢，就是在该频率下的信号的相位。 第一个点表示直流分量（即0Hz），而最后一个点N的再下一个 点（实际上这个点是不存在的，这里是假设的第N+1个点，也 可以看做是将第一个点分做两半分，另一半移到最后）则表示 采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率 依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。 由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果 采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时 间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率
按时间抽取的基2FFT算法分析
第四章 快速傅里叶变换 有限长序列可以通过离散傅里叶变换(DFT)将其频域也离散化成有限长序列.但其计算量太大,很难实时地处理问题,因此引出了快速傅里叶变换(FFT). 1965年，Cooley 和Tukey 提出了计算离散傅里叶变换（DFT ）的快速算法，将DFT 的运算量减少了几个数量级。从此，对快速傅里叶变换（FFT ）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT 的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT 的多种算法，基本算法是基２DIT 和基２DIF 。FFT 在离散傅里叶反变换、线性卷积和线性相关等方面也有重要应用。 快速傅里叶变换（FFT ）是计算离散傅里叶变换（DFT ）的快速算法。 DFT 的定义式为 )(k X =)()(1 0k R W n x N N n kn N ∑-= 在所有复指数值kn N W 的值全部已算好的情况下，要计算一个)(k X 需要N 次复数乘法和N －1次复数加法。算出全部N 点)(k X 共需2 N 次复数乘法和)1(-N N 次复数加法。即计算量是与2 N 成正比的。 FFT 的基本思想：将大点数的DFT 分解为若干个小点数DFT 的组合，从而减少运算量。 N W 因子具有以下两个特性，可使DFT 运算量尽量分解为小点数的DFT 运算： （1） 周期性：k N n N kn N n N k N W W W )()(++== （2） 对称性：k N N k N W W -=+) 2/(
利用这两个性质，可以使DFT 运算中有些项合并，以减少乘法次数。例子：求当N ＝4时，X(2)的值 ) ()]3()1([)]2()0([)()]3()1([)]2()0([)3()2()1()0()()2(0 424046 44424043 24对称性－＝周期性W x x x x W x x W x x W x W x W x W x W n x X n n +++++=+++==∑= 通过合并，使乘法次数由4次减少到1次，运算量减少。 FFT 的算法形式有很多种，但基本上可以分为两大类：按时间抽取（DIT ）和按频率抽取（DIF ）。 4.1 按时间抽取（DIT ）的FTT 为了将大点数的DFT 分解为小点数的DFT 运算，要求序列的长度N 为复合数，最常用的是M N 2=的情况（M 为正整数）。该情况下的变换称为基2FFT 。下面讨论基2情况的算法。 先将序列x(n)按奇偶项分解为两组 ?? ?=+=) ()12()()2(21r x r x r x r x 12,,1,0-=N r Λ 将DFT 运算也相应分为两组 = =)]([)(n x DFT k X ∑-=1 )(N n kn N W n x ∑∑-=-== 1 n 0 1 0)()(N n kn N N n n kn N W n x W n x 为奇数为偶数＋ ∑∑-=+-=++ 1 2/0 )12(1 2/0 2)12()2(N r k r N N r rk N W r x W r x ＝
fft的发展、现状、典型算法
数字信号处理期末大作业FFT的发展史、现状及典型算法 班级 学号： 姓名：
FFT的发展史、现状及典型算法 傅里叶分析已有200多年的历史，目前FFT及其校正算法在工程实际中仍在广泛应用，展现了其不竭的生命力。本次作业我们论述FFT的现状，发展史以及一些算法，去详细了解、扩展这一算法，巩固所学知识。 一．FFT的简介 傅里叶变换是一种将信号从时域变换到频域的变换形式，然而当N很大的时候，求一个N点的DFT要完成N*N次复数乘法和N*（N-1）次复数加法，计算量非常大，所以人们开始探索一种简便的算法对于一个较大的N进行傅里叶变换。在20世纪60年代由Cooley和Tukey提出了快速傅里叶变换算法，它是快速计算DFT的一种简单高效的方法。 关于何为FFT，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。举个例子，设x(n)为N项的复数序列，由DFT 变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2*（N/2)^2=N+（N^2）/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。 所以使用FFT算法，可以大大提高傅里叶变换的运算速度，运算时间缩短一到两个数量级，从而使DFT变换应用迅速普及，不仅在频谱分析，而且在线性卷积、线性相关等方面得到广泛应用。
用matlab实现fft算法
A1=str2double(get(handles.edit8,'String')); A2=str2double(get(handles.edit9,'String')); F1=str2double(get(handles.edit10,'String')); F2=str2double(get(handles.edit11,'String')); Fs=str2double(get(handles.edit12,'String')); N=str2double(get(handles.edit13,'String')); t=[0:1/Fs:(N-1)/Fs]; x=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t); %信号x的离散值 axes(handles.axes1) %在axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x) ARM或单片机用的FFT算法，用于信号处理。经过愚昧的本人优化，提高了计算效率 （在aduc7026 @41MHz FFT256点环境下运算时间为0.06s左右） PS：以下有两部分（fft.h和fft.c） 【复制以下内容改名为fft.h】 #ifndef __FFT_H__ #define __FFT_H__ #include 按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。 按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 对N=M 2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p W N ，p 称为旋转因子的指数。例如N ＝8 ＝3 2 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：p W N =J /4W N =J 2L W J=0 第二级：L=2，有2个旋转因子：p W N =J /2W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N ＝M 2的一般情况，第L 级共有1 -L 2 个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1 -L 2－1 L 2=M 2×M -L 2 = N ·M -L 2 故： 按照上面两式可以确定第L 级运算的旋转因子 基2的DIT蝶形算法源代码及注释如下： /************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间 #include 利用FFT 计算卷积 一．线卷积的作用及定义 线卷积包括卷积积分和卷积和。 1．线卷积的作用 求解线性系统对任意激励信号的零态响应。 2．卷积积分 ) (*)(d )()()(t h t x t h x t y =-= ? ∞∞ -τττ 3．卷积和 离散系统的时域分析是，已知离散系统的初始状态和输入信号（激励），求离散系统的输出（响应），两种方法：递推解法和离散卷积法。 卷积和：)()()()()(n h n x m n h m x n y m *=-= ∑ ∞ -∞ = 二．圆周卷积的定义 圆周移位：一周期为N 的周期序列, 可视为一主值序列在圆周上的循环移位。周期序列在时间轴上左移 右移m 反时针 转称为圆周移位。 时域圆周卷积（循环卷积） )()()(n h n x n y ?=()()()∑ -=-= 1 )(N m N N n R m n h m x 条件：两序列实现圆卷积的条件是：长度相等，如果不相等, 可通过增补零值来使之相等。 特点：卷积求和范围只在10-≤≤N m 有限区间进行；卷积时不作反褶平移, 而是反褶圆移 步骤：量置换→反褶→圆移→相乘→求和。 三．两者的关系 有限长序列的圆卷积和线卷积的关系 在一般情况下，两序列的圆卷积和线卷积是不相等的，这是因为：线卷积是 平移, 结果长度为121-+=N N L ；而圆卷积是圆移，结果长度为2 1 N N L ==。只有 在两卷积的结果长度相时，二者才有相同的结果。解决方法是：在作圆卷积时，通过加零的方法，使两序列的长度都增加到121-+=N N L ，此时，圆卷积的结果和线卷积同。 四．利用FFT 计算卷积 工程实际需要解决的卷积：)()()(n h n x n y *=，但其计算量很大。 而圆卷积为：)()()(n h n x n y ?=，便于采用FFT 算法, 故计算速度快。若将线卷积的两个序列用增补零的方法将长度取为一致，此时两序列的离散线卷积和圆周卷积结果是相等的，这样就则可以通过圆卷积来快速计算线卷积。 1、 利用FFT 计算卷积的步骤 （1）设两序列原长度分别为：N 和M ，将长度增加到1-+≥M N L （L 为2的整数次幂）； （2）用FFT 法求加长序列的DFT 频谱； （3）计算两序列DFT 频谱的乘积； （4）用IFFT 求DFT 频谱乘积的逆变换，便得两序列的离散线卷积。 2、分段快速卷积 设)(n x 为长序列，)(n h 为短序列，长度为M ，则两序列的离散线卷积可以写成如 下 形 式 ， ∑∑∑-=-+=-=+-+ +-+ -= *=1 1 )1(1 2)()()()()()()()()(N m n K kN m N N m m N h m x m N h m x m N h m x n h n x n y 上述每个子段长度为N 。为便于圆卷积计算，将长度通过补零加长为：1-+=M N L x (n 0 n h (n 根据各子段()n x k 增补零的部位不一样而分两种算法。 /****************************************************************************/ /* */ /* 快速傅里叶变换/ 快速傅里叶逆变换测试*/ /* */ /* 2014年04月20日*/ /* */ /****************************************************************************/ #include 常用FFT 算法总结 一、DFT 及IDFT 的定义 对于N 点有限长序列)(n x ，其DFT 及IDFT 定义如下： DFT ：∑-==1 0)()(N n nk N W n x k X IDFT ：∑-=-=1 )(1)(N k nk N W k X N n x 其中，N p j p N e W /2π-=称为旋转因子，p 称为旋转因子的指数。 二、基2-FFT 算法 设序列)(n x 的长度N 满足M N 2 =，M 为自然数。 1、时间域抽取FFT （DIT-FFT ） （1）算法原理 按n 的奇偶把)(n x 分解为两个N/2点的子序列： ) 12()()2()(21+==r x r x r x r x 则)(n x 的DFT 为 )()()12()2()(211 2/0 )12(1 2/0 2k X W k X W r x W r x k X k N N r k r N N r rk N +=++ = ∑∑-=+-= 其中)(1k X 和)(2k X 分别为)(1r x 和)(2r x 的N/2点DFT 。 利用)(1k X 和)(2k X 的周期性，)(k X 可以表示为 ?? ? ??-=++=)()()2()()()(2121k X W k X N k X k X W k X k X k N k N 上式表明一个N 点的DFT 可以用两个N/2点的DFT 来表示。 （2）运算量 当M N 2 =时，共有M 级蝶形，每级N/2个蝶形，每个蝶形有1次复数乘法2次复数加法。 复数乘法：N N M N m F 2 log 2 2 = = 复数加法：N N NM a F 2 log == N N N N N FFT m DFT m F F 2 2 2 log 2log 2 ) ()(= = DFT 算法原理 快速傅氏变换（FFT ）是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设x(n)为N 项的复数序列，由DFT 变换，任一X （m ）的计算都需要N 次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的X （m ）,即N 点DFT 变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT 中，利用WN 的周期性和对称性，把一个N 项序列（设N=2k,k 为正整数），分为两个N/2项的子序列，每个N/2点DFT 变换需要（N/2）2次运算，再用N 次运算把两个N/2点的DFT 变换组合成一个N 点的DFT 变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要Nlog2N 次的运算，N 在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT 的优越性。 、FFT 的算法原理 FFT 算法的输出X(K)为自然顺序，但为了适应原位计算，其输入序列不是按x(n)的自然顺序排序，这种经过M-1次奇偶抽选后的排序为序列的倒序。因此，在运算之前应先对序列x(n)进行倒序。倒序的规律就是把顺序数的二进制位倒置，即可得到倒序值。倒序数是在M 位二进制数最高位加一，逢2向右进位。对于 ，M 位二进制数最高位的权值为N/2，且从左到右二进制位的权值依次为 你N/4，N/8，···，2，1。因此，最高位加一相当于十进制运算J+N/2。（J 表示当前倒序数的十进制数值） 实验原理与方法FIR 滤波器 FIR 滤波器的设计问题在于寻求一系统函数)(z H ，使其频率响应)(ω j e H 逼近滤波器 要求的理想频率响应)(ω j d e H ，其对应的单位脉冲响应)(n h d 。 1．用窗函数设计FIR 滤波器的基本方法 设计思想：从时域从发，设计)(n h 逼近理想)(n h d 。设理想滤波器)(ω j d e H 的单位脉 冲响应为)(n h d 。以低通线性相位FIR 数字滤波器为例。 一．FFT 算法思想： 1．DFT 的定义： 对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{x[k]}可以由离散付氏变换（DFT ）求得。DFT 的定义为： 21 [][]N j nk N n X k x n e π--==∑，k=0,1,…N-1 通常令2j N N e W π-=，称为旋转因子。 2．直接计算DFT 的问题及FFT 的基本思想： 由DFT 的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N 点DFT 需要（N-1）2次复数乘法和N （N-1）次加法。因此，对于一些相当大的N 值（如1024）来说，直接计算它的DFT 所作的计算量是很大的。 FFT 的基本思想在于，将原有的N 点序列分成两个较短的序列，这些序列的DFT 可以很简单的组合起来得到原序列的DFT 。例如，若N 为偶数，将原有的N 点序列分成两个（N/2）点序列，那么计算N 点DFT 将只需要约[(N/2)2 ·2]=N 2/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2代表必须完成两个DFT 。上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点的DFT （假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT 运算的情况。 3．基2按时间抽取（DIT ）的FFT 算法思想： 设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。 将长度为2L N =的序列[](0,1,...,1)x n n N =-，先按n 的奇偶分成两组： 12[2][][21][] x r x r x r x r =+=，r=0,1,…,N/2-1 DFT 化为： 1/21 /21 2(21)0 /21 /21 22120 0/21 /21 1/2 2/2 []{[]}[][2][21][][][][]N N N nk rk r k N N N n r r N N rk k rk N N N r r N N rk k rk N N N r r X k DFT x n x n W x r W x r W x r W W x r W x r W W x r W ---+===--==--===== + +=+= +∑∑ ∑ ∑ ∑ ∑ ∑ 上式中利用了旋转因子的可约性，即：2/2rk rk N N W W =。又令 /21 /21 11/2 22/20 [][],[][]N N rk rk N N r r X k x r W X k x r W --=== = ∑ ∑ ，则上式可以写成： 按时间抽取的基2FFT算法分析及MATLAB实现 二、 三、DIT-FFT算法的运算规律及编程思想 1.原位计算 对N=M2点的FFT共进行M级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M级运算后，原来存放输入序列数据的N个存储单元中可依次存放X(k)的N个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N点DIT―FFT运算流图中，每个蝶形都要乘以旋转因子p W N ，p称为旋转因子的指数。例如N＝8 ＝32时各级的旋转因子： 第一级：L=1，有1个旋转因子：p W N =J /4 W N =J 2L W J=0 第二级：L=2，有2个旋转因子：p W N =J /2 W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N＝M2的一般情况，第L级共有1-L2个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1-L2－1 L 2=M2×M-L2= N·M-L2 故：按照上面两式可以确定第L级运算的旋转因子 3、同一级中，同一旋转因子对应蝶形数目 第L级FFT运算中，同一旋转因子用在L-M2个蝶形中； 4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L级中，蝶距：D=L2； 5、同一蝶形运算两输入数据的距离 在输入倒序，输出原序的FFT变换中，第L 级的每一个蝶形的2个输入数据相距：B=1-L2。 数字信号处理实验报告 实验二 FFT 算法的MATLAB 实现 （一）实验目的：理解离散傅立叶变换时信号分析与处理的一种重要变换，特别是FFT 在数字信号处理中的高效率应用。 （二）实验原理： 1、有限长序列x(n)的DFT 的概念和公式： ??? ????-≤≤=-≤≤=∑∑-=--=1 01 01 0)(1)(10)()(N k kn N N n kn N N n W k x N n x N k W n x k x ) /2(N j N e W π-= 2、FFT 算法 调用格式是 X= fft(x) 或 X=fft(x,N) 对前者，若x 的长度是2的整数次幂，则按该长度实现x 的快速变换，否则，实现的是慢速的非2的整数次幂的变换；对后者，N 应为2的整数次幂，若x 的长度小于N ，则补零，若超过N ，则舍弃N 以后的数据。Ifft 的调用格式与之相同。 （三）实验内容 1、题一：若x(n)=cos(n*pi/6)是一个N=12的有限序列，利用MATLAB 计算它的DFT 并画出图形。 源程序： clc; N=12; n=0:N-1; k=0:N-1; xn=cos(n*pi/6); W=exp(-j*2*pi/N); kn=n'*k Xk=xn*(W.^kn) stem(n,Xk); xlabel('k'); ylabel('Xk'); grid on; 也可用FFT 算法直接得出结果，程序如下： clc; N=12; n=0:N-1; xn=cos(n*pi/6); Xk=fft(xn,N); stem(n,Xk); xlabel('k'); ylabel('Xk'); grid on; 实验结果： 24 681012 k X k 分析实验结果： 用DFT 和用FFT 对序列进行运算，最后得到的结果相同。但用快速傅立叶变换的运算速度可以快很多。 2、题二：一被噪声污染的信号，很难看出它所包含的频率分量，如一个由50Hz 和120Hz 正弦信号构成的信号，受均值随机噪声的干扰，数据采样率为1000Hz ，通过FFT 来分析其信号频率成分，用MATLAB 实现。 源程序： clc; fs=1000; N=1024; n=0:N-1; t=n/fs; x=sin(2*pi*50*t)+sin(2*pi*120*t)+rand(1,N); y=fft(x,N); mag=abs(y); f=n*fs/N; subplot(1,2,1),plot(f,mag); xlabel('频率/Hz'); 第五章 离散傅里叶变换及其快速算法 1 离散傅里叶变换(DFT)的推导 (1) 时域抽样： 目的：解决信号的离散化问题。 效果：连续信号离散化使得信号的频谱被周期延拓。 (2) 时域截断： 原因：工程上无法处理时间无限信号。 方法：通过窗函数(一般用矩形窗)对信号进行逐段截取。 结果：时域乘以矩形脉冲信号，频域相当于和抽样函数卷积。 (3) 时域周期延拓： 目的：要使频率离散，就要使时域变成周期信号。 方法：周期延拓中的搬移通过与)(s nT t -δ的卷积来实现。 表示：延拓后的波形在数学上可表示为原始波形与冲激串序列的卷积。 结果：周期延拓后的周期函数具有离散谱。 (4) 1。 图1 DFT 推导过程示意图 (5) 处理后信号的连续时间傅里叶变换：∑∑ ∞ -∞=-=π--δ???? ? ????= k N n N kn j s kf f e nT h f H )()()(~ 010/2 (i) )(~f H 是离散函数，仅在离散频率点S NT k T k kf f = ==00处存在冲激，强度为k a ，其余各点为0。 (ii) )(~ f H 是周期函数，周期为s s T NT N T N Nf 1 00= == ，每个周期内有N 个不同的幅值。 (iii) 时域的离散时间间隔(或周期)与频域的周期(或离散间隔)互为倒数。 2 DFT 及IDFT 的定义 (1) DFT 定义：设()s nT h 是连续函数)(t h 的N 个抽样值1,,1,0-=N n ，这N 个点的宽度为 N 的DFT 为：[])1,...,1,0(,)()(1 0/2-=??? ? ? ?==? -=π-∑N k NT k H e nT h nT h DFT s N n N nk j s s N (2) IDFT 定义：设??? ? ??s NT k H 是连续频率函数)(f H 的N 个抽样值1,,1,0-=N k ， 这N 个点的宽度为N 的IDFT 为： ())1,...,1,0(,11 0/21 -==??? ? ? ?=???????????? ???-=π--∑ N k nT h e NT k H N NT k H DFT s N k N nk j s s N (3) N nk j e /2π-称为N 点DFT 的变换核函数，N nk j e /2π称为N 点IDFT 的变换核函数。它们 互为共轭。 (4) 同样的信号，宽度不同的DFT 会有不同的结果。DFT 正逆变换的对应关系是唯一的， 或者说它们是互逆的。 (5) 引入N j N e W /2π-= (i) 用途： (a) 正逆变换的核函数分别可以表示为nk N W 和nk N W -。 (b) 核函数的正交性可以表示为：() )(* 1 0r n N W W kr N N k kn N -δ=∑-= (c) DFT 可以表示为：)1,,1,0(,)(10 -==? ??? ??∑ -=N k W nT h NT k H N n nk N s s (d) IDFT 可以表示为：)1,,1,0(,1 )(1 0-=??? ? ? ?= ∑ -=-N n W NT k H N nT h N k nk N s s (ii) 性质：周期性和对称性： (a) 12==π-j N N e W (b) 12 /-==π-j N N e W (c) r N r N N N r N N W W W W ==+ (d) r N r N N N r N N W W W W -=-=+2/2/ (e) )(1Z m W m N ∈?= (f) ),(/2/2Z n m W e e W n N N n j m N m n j m n m N ∈?===π-π- 3 离散谱的性质 (1) 离散谱定义：称)(Z k NT k H H S k ∈???? ? ?=? 为离散序列)0)((N n nTs h <≤的DFT 离散谱，简称离散谱。 (2) 性质： (i) 周期性：序列的N 点的DFT 离散谱是周期为N 的序列。 (ii) 共扼对称性：如果)0)((N n nTs x <≤为实序列，则其N 点的DFT 关于原点和N /2都 实验6 FFT 算法的应用 实验目的：加深对离散信号的DFT 的理解及其FFT 算法的运用。 实验原理：N 点序列的DFT 和IDFT 变换定义式如下： 1 [][]N kn N n X k x n W -== ∑ ， 1 1[][]N kn N k x n X k W N --== ∑ 利用旋转因子2j nk kn N N W e π-=具有周期性，可以得到快速算法（FFT ）。 在MATLAB 中，可以用函数X=fft （x ，N ）和x=ifft （X ，N ）计算N 点序列 的DFT 正、反变换。 例1 对连续的单一频率周期信号 按采样频率 采样，截取长度N 分别选 N =20和N =16，观察其DFT 结果的幅度谱。 解 此时离散序列 ，即k=8。用MATLAB 计 算并作图，函数fft 用于计算离散傅里叶变换DFT ，程序如下： k=8; n1=[0:1:19]; xa1=sin(2*pi*n1/k); subplot(2,2,1) plot(n1,xa1) xlabel('t/T');ylabel('x(n)'); xk1=fft(xa1);xk1=abs(xk1); subplot(2,2,2) stem(n1,xk1) xlabel('k');ylabel('X(k)'); n2=[0:1:15]; xa2=sin(2*pi*n2/k); subplot(2,2,3) plot(n2,xa2) xlabel('t/T');ylabel('x(n)'); xk2=fft(xa2);xk2=abs(xk2); subplot(2,2,4) stem(n2,xk2) xlabel('k');ylabel('X(k)');FFT算法(查表法)
按时间抽取的基2FFT算法分析及MATLAB实现
fft算法代码注释及流程框图
利用FFT计算卷积
FFT算法
常用FFT算法
DFT算法原理、FFT的算法原理
FFT算法思想
按时间抽取的基2FFT算法分析及MATLAB实现
FFT算法(用matlab实现)
FFT算法
实验6 FFT算法的应用