第四章 快速傅里叶变换(FFT)

第四 章  快速傅里叶变换（ FFT）通信 1402班   肖太龙ContentsØ 从 DFT到 FFTØ FFT的基本概念Ø FFT算法的实现原理Ø FFT的物理意义  ØDFT与 FFT时间复杂度比较从 DFT到 FFT的缘由  依次 类 推， N/2个点的 时间 复 杂 度 为多少呢？ N/4呢？从哪个方面减少运算量呢？ 周期性表现为： 对称性表现为： 1965年库利和图基提出一种 DFT的快速算法，以此来解决 DFT变换的效率问题。 算法的提出的初衷主要是为了减少乘法次数，又因为旋转因子的对称性和周期性能达到这一目标，因此在理论上是可以实现的。 FFT算法就是不断地把长序列的 DFT分解为短序列的DFT的算法。最常用的是基2FFT算法。浅谈 FFT实现原理lFT算法基本上分为两大类 ：l时域抽取法 FFT(Decimation In Time FFT,简称 DIT-FFT)。l频域 抽取法 FFT(Decimation In Frequency FFT,简称DIFFFT)。l我们主要来分析 DIF-FFT算法。A 算法理论推导B 一个简单的算例 8点 DFT一次时域抽取分解C Matlab绘制图像DIF-FFT与 DIT-FFT算法有什么异同 ？算法理 论 推 导 ：设 序列 的 长 度 为 N，且 满 足 N=2M ， M为 自然数。按 n奇偶把 分解 为 两个 N/2点的子序列则 x(n)的 DFT为 因 为所以 其中 X1(k)和 X2(k)分 别为 x1(r)和 x2(r)的 N/2点 DFT， 即由于 X1(k)和 X2(k)均以 N/2为 周期， 且所以 X(k)又可表示 为至此，一次 时 域抽取法的理 论 推 导 就完成了。从上面的两个式子中，我 们 定 义 一种新的运算符 蝶形运算符。8点 DFT一 次 时 域抽取分解运算流 图完成一次蝶形运算需要一次复数乘法和两次复数加法运算。因此在 8点 DFT一 次 时 域抽取分解中，共需要 计算两个 N/2点 DFT运算和 N/2个蝶形运算。所以按照上 图 的 计 算 DFT时 ， 总 的复数乘法次数为复数加法次数 为类 似地，我 们 将 按奇偶分解成两个 N/4点子序列，其表达式分 别 如下：那么 ， 又 可表示 为式中同理， 由 的 周期性 和 的 对 称性最后 得到：用同 样 的方法可以 计 算出其中8点 DFT二次 时 域抽取分解运算8点 DIT-FFT运算流 图从上面的蝶形算法，当 时 ，其运算 应该 有M级 蝶形，每一 级 都由 N/2个蝶形运算构成。因此每一 级 运算都需要 N/2次复数乘法和 N次复数加法，所以 M级 蝶形运算的乘法次数 为 ：加法次数：一个 简单 的算例 计 算序列 的 8点 DFT。分 别用基本的 DFT算法和 FFT算法 实现 ，体会 计 算 过 程中的 时间 复 杂 度。 当然 针对计 算机来 说 ， 计 算乘法所需要的 时间远 大于加法。一般的 计 算机，差不多相差十倍左右。 用 计 算机 产 生随机的 1024个点构成序列，然后取N=1024.此 时计 算 时间 差距就会加大。 N=2048时 ， 时间 差距会更大。为 了更 进 一步的体会 FFT的物理意 义 ，引入一个算例：假 设对 某信号 经过 ADC之后，得到一个序列 ，此 时 我们 不知道其具体的函数表达式。但是我 们 可以 对 其做 FFT运算。 现 在我 们 做一个 验证 ： 假 设 我 们 有一个信号，它含有 2V的直流分量， 频 率 为 50Hz、相位 为 -30度、幅度 为 3V的交流信号，以及一个 频 率 为 75Hz、相位 为 90度、幅度 为 1.5V的交流信号。用数学表达式就是如下：S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180) 现 在 对 其 进 行 变换 ，取 样 点 为 1024，采 样频 率 为1024Hz 注意 这 里的取 样频 率只要 满 足原始 频 率的 2倍即可，且取样 点和取 样频 率根据 频 率分辨率来 选 取。从 图 中的 结 果可以得出，当 频 率 为 0、 50、 75时 ， 对应 的幅度 值 依次 为 2、 3、 1.5，相位依次 为 0、 90、 -30。当然， 这 看上去几乎是没有 误 差的，因 为 我 们 取 样频 率和所取点数比 较 大，当 N=Fs=256时 ，会存在 误 差，但是 这 个 误差完全不影响我 们对 信号函数的分析与判断。从而 进 一步的验证 了 DFT与 FFT算法的正确性。上 图 是从模 拟 信号到 频 域离散信号的完整的 过 程。 这 其中对应 有几个概念容易混淆。因此 对 此做出区分：数字 频 率 、 模 拟频 率与采 样频 率：模 拟频 率通常用 表示，数字 频 率用 表示。此 时 的数字 频 率主要是 针对 序列而言，因此没有采 样频 率 就不会有数字 频 率的概念，所以数字 频 率与采 样频 率和模 拟频 率一定 满 足某种关系。我们 知道有如下 过 程：因此其中 为 取 样间 隔。 对应 的采 样频 率 DFT取 样 点 N与采 样频 率和 频 率分辨率（步 进值 ）的关系 首先我 们 根据奈奎斯特定律得到： ， 为连续 信号的截至（最高） 频 率 （ 因 为 它可能有很多 频 率成分）。 现 在 问题 就是取 样 点 N的 选 取？ 现 在我 们 从序列的 FT说 起！ FT完成的任 务 就是 对时 域 连续 信号抽 样 之后通 过 序列的傅里叶 变换 得到的 频谱 关系，我 们 知道 这 个 谱 是 连续的。 此 时 的 谱 是一个周期 谱 ，那么周期是多少呢？是不是与 时 域采 样频 率有关系呢？答案是肯定的，此 时 的周期就是 （注意此 时 是 频 域， 频谱图 中的横坐 标 表示的 频 率）。对 于周期 连续 的 谱 ， 计 算机 还 是无法操作，因此我 们还 得 对 FT变换 之后的 谱进 行 频 域采 样 （ DFT的 实质 ） 。 类似于 时 域采 样 ， 频 域 FT谱线 的采 样频 率 为 多少合适呢？ 显然此 时 只需一个周期即可，即在一个周期中取 N个点， 满 足 ： 式 中 为频 率分辨率， 为 采 样频 率。因此关于 频 率分辨率和取 样 点 N的关系就得出来了。 首先必 须 得画出序列 对应 的 FT谱线图 ，横坐 标为 ， 纵坐 标为频 率。此 时 的 FT谱线应该 是 对 上 图 的周期延拓。 这就是我 们 一般只取其主 值 -pi,pi得原因。得到 FT谱线图 之后，接下来就是 DFT对 其 进 行取 样 ！由于 DFT的物理意 义 就是 对 0,2*pi的等 间 隔取 样 。即从 题 目 给 定的 频谱图 得出模 拟频 率 应该 是 连续分布的。但是离散化之后只能得到离散的模 拟频率，由 确定。但是不能超 过 200Hz。DFT算法与 FFT算法 耗 时 对 比 基于 matlab编 程 DFT函数 function Xk=dft(xn,N) n=0:1:N-1;k=0:1:N-1;WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.nk;Xk=xn*WNnk;endDFT测试 代 码 ：clcclearN=4096;xn=rand(1,4096); ticXk=dft(xn,4096);tocElapsed time is 10.993700 seconds.由于 Matlab工具箱自 带 FFT算法，因此 调 用 库 函数即可。但是可能存在一定的 问题 。 库 函数里面的 FFT算法性能可能会更加 优 化，耗 时 会更少。N=4096;xn=rand(1,4096);ticXk=fft(xn,4096);tocElapsed time is 0.002649 seconds.从耗 时 来看，同 样 是 4096个点， fft算法所需 时间远 小于 DFT算法 时间 。（做一次 实验 的 结 果）由于 FFT具体的算法不得而知，因此 该实验 只能 说 明， FFT快速算法大大减小了 时间 复 杂 度。由于自 带 的函数无法 显 示具体的 过 程，且与DFT可比性不 强 ，因此需要 编 写 FFT算法代 码显 示所需 时间 。通 过 运行代 码 ，得到耗 时为 0.020574秒。频谱 图 ：for L=1:MB=2(L-1);for R=0:B-1P=2(M-L)*R;for K=R:2L:N/2T=A(K+1)+A(K+B+1)*WN(P+1);A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);A(K+1)=T;endendendtock=0:4095;Xk=fft(xn,4096);subplot(2,1,1)plot(k,abs(A) xn=ones(1,4096);ticM=nextpow2(length(xn);N=2M;for m=0:N/2-1WN(m+1)=exp(-j*2*pi/N)m;endA=xn,zeros(1,N-length(xn);J=0;for I=0:N-1if I=KJ=J-K;K=K/2;endJ=J+K;end源代 码谢谢 大家！