快速傅里叶变换.ppt

数字信号处理,2/62,为了分析连续信号频谱，需要使用DFT来逼近连续时间信号的傅里叶变换。但由于DFT要对连续信号进行采样和截断等处理，因此，常产生误差分析，体现在三个方面。 混叠现象：由于不满足抽样定理而引起； 截断效应：由于序列截断引起； 栅栏效应：由于对频谱的有限点采样引起。,信号逼近所产生的问题,回顾,3/62,回顾 频谱混叠,1. 混叠现象,4/62,回顾 频谱混叠,原因：时域采样不满足奈奎斯特准则： 其中 为抽样频率， 为信号最高频率。 此条件只规定了 的下限，其上限要受频率分辨率（抽样间隔） 的约束。 也可表示为记录长度的倒数，即： 其中， 为抽样点数。,措施：采样前对信号进行预滤波，采样时满足采样定理，以避免信号在w=处附近的混迭，再滤去信号中频率高于 的频率分量。,5/62,回顾 频谱混叠,需要注意的是： 当 给定时, 若为了避免混叠现象而一味提高抽样频率 , 导致 增加，即导致频率分辨能力下降; 反之, 若打算要提高频率分辨率即减小 , 则导致减小 , 最终必须减小信号的高频容量。 显然，针对高频容量 与频率分辨率 , 想保持其中一个不变而使另一个性能得以提高的唯一办法是： 增加记录长度内的点数N, 即： 这是未采用任何特殊数据处理（例如加窗）情况下，为实现基本DFT算法所必须满足的基本条件。,6/62,2. 截断效应,回顾 截断效应,窗函数不可能取无限宽，即其频谱不可能为一冲激函数，信号的频谱与窗函数的卷积必然产生拖尾、变宽的现象，从而造成截断效应 。如下图所示。,7/62,“泄漏”是由矩形窗函数带来的，,求极限：,的连续频谱。,为,处的一根谱线，变成了以,原来在,为中心，形状,可见：,回顾 截断效应,8/62,小、主瓣窄的窗函数。,接近,为了尽量减少泄漏，需要寻找频域中窗函数,具有旁瓣，并且主瓣也,不过，泄漏的产生是由于,所以不能用无限加宽窗口来减少泄漏。,，意味着需要无限加宽窗函数，等于未截短。,若使泄漏至0，则,回顾 截断效应,占有一定宽度。,，即旁瓣,9/62,回顾 栅栏效应,N点DFT是在区间0, 2上的N点等间隔采样，采样点之间的频谱函数值是不知道的，就好像从N+1个栅栏缝隙中观看信号的频谱特性，得到的是N个缝隙中看到的频谱函数值，这种现象称为栅栏效应。 原因：对信号的频谱进行有限点采样； 后果：栅栏效应可能漏掉（挡住）大的频谱分量； 措施：对原序列补0，增大N，以增加采样点。,3. 栅栏效应,管中窥豹,10/62,减小栅栏效应的一个方法是，在所取数据的末端添加一些零值点，使一个周期内点数增加, 但是不改变原有的记录数据。 这种方法相当于加长了周期T0，T0增加, 使得频域抽样间隔变小, 从而能保持原来频谱形式不变的情况下使谱线变密, 也就使频谱抽样点数增加。 这样，原来看不到的频谱分量就有可能看到了。,回顾 栅栏效应,减小栅栏效应的方法 补零,11/62,回顾 栅栏效应,后面补4个零值，则有：,12/62,【解】,例 1： 栅栏效应,已知：,13/62,例 1： 栅栏效应,x=2,3,3,2; subplot(2,2,1); stem(x,fill);title(x(n); N=8*length(x); NFFT = 2nextpow2(N); Xk=fft(x,NFFT); subplot(2,2,2); f = 0:2*pi/(NFFT-1):2*pi; stem(f,abs(Xk); title(X(k) DFT点数 N= num2str(NFFT); x1=2,3,3,2,0,0,0,0; subplot(2,2,3);stem(x1,fill); title(x1(n); N1=8*length(x1); NFFT1 = 2nextpow2(N1); X1k=fft(x1,NFFT1); subplot(2,2,4); f1 = 0:2*pi/(NFFT1-1):2*pi; stem(f1,abs(X1k); title(X1(k) DFT点数 N= num2str(NFFT1);,14/62,注意： 补加零点以改变周期时, 所用窗函数宽度却不能变。亦即必须按数据记录原长来选择窗函数, 而不能按补了零值点后的长度来选择窗函数。 即：“先加窗, 再补零。”,回顾 栅栏效应,15/62,一般来说，信号长度 越长，抽样点数 越大，则 越小，即分辨力越强。但此处的 是指实际的信号长度， 是指这个长度上的抽样点数，而不是进行补零等辅助措施后的信号长度及抽样点数。,回顾 分辨率,频率分辨率,所谓频率分辨率，是指能够将信号中两个靠得很近的谱峰分开的能力：,16/62,由上两式可见，当采样频率由,变为,时，在数据长度,变为,，频率采样间隔不变。,不变情况下，采样点数由,仅增加采样率并不能改变频率分辨率。,如果数据长度,不变，,则有：,注意（1）：单纯增加采样率并不能改变频率分辨率 ！,回顾 分辨率,分析：,显然，只有“增加点数的同时导致数据有效长度的增加”才能使分辨率越好。,17/62,例 2： 分辨率,已知某信号如下： 假设最高频率 ，请用 对其抽样，并观察其频谱情况。,【解】 假设取信号长度为： ，经抽样所得的 的点数为： 对 做DFT变换，所得到的频率最大分辨率为：,18/62,例 2： 分辨率,clear all; L=256; % 信号长度 N=256; % DFT点数 n=0:1:L-1; f1=2;f2=2.02;f3=2.07 x=sin(2*pi*f1*0.1*n)+sin(2*pi*f2*0.1*n)+sin(2*pi*f3*0.1*n); subplot(2,1,1);stem(x); title(x(n) 信号长度 tp= num2str(L); X=fft(x,N);Xk=abs(X(1:N/2); k=0:1:N/2-1;w=2*pi/N*k; subplot(2,1,2);plot(5*w/pi,Xk/N); title(DFT 点数N= num2str(N); axis(1,3,0,1);,现象： 只显示了两个频率：f1、f3 原因： f2-f1=0.02F0，所以能分辨出来 改进方法： 增加采样点数 N，即增加数据长度。,19/62,例 2： 分辨率,clear all; L=1024; % 信号长度 N=1024; % DFT点数 ,20/62,注意（2）：并非补零后就能提高频率分辨率 ！ 设原数据长度T1，抽样点数N1，补零后的数据长度T2，抽样点数N2，则：,显然，N2 N1，故F2 F1。 但并不代表：“补零后，频率分辨率提高了！” 因为：补零不能增加数据的有效长度，上面实际数据的有效长度仍为T1（有效抽样点数仍为N1）。因此，补零是不能提高频率分辨率的。,回顾 分辨率,21/62,例 3： 分辨率,求出序列x(n)=cos(0.48n)+cos(0.52n)基于有限个样点n=10的频谱; 求n=100时，取x(n)的前10个，后90个设为零，得到x(n)的频谱; 增加x(n)有效的样点数，取100个样点得到x(n)的频谱。,22/62,例 3： 分辨率,x(n)基于10个样点的频谱 figure(1) n=0:1:99; x=cos(0.48*pi*n)+cos(0.52*pi*n); n1=0:1:9;y1=x(1:1:10); subplot(2,1,1);stem(n1,y1);title(signal x(n), 0 = n = 9);xlabel(n) axis(0,10,-2.5,2.5) Y1=fft(y1);magY1=abs(Y1(1:1:6); k1=0:1:5;w1=2*pi/10*k1; subplot(2,1,2);stem(w1/pi,magY1);title(10点DFT); xlabel(w/pi)，axis(0,1,0,10),23/62,例 3： 分辨率,在10个样点的基础上添90个零，得到密度高的频谱 figure(2) n3=0:1:99;y3=x(1:1:10) zeros(1,90); 添90个零。得到100个数据 subplot(2,1,1);stem(n3,y3);title(signal x(n), 0 = n = 9 + 90 zeros);xlabel(n) axis(0,100,-2.5,2.5) Y3=fft(y3);magY3=abs(Y3(1:1:51); k3=0:1:50;w3=2*pi/100*k3; subplot(2,1,2);stem(w3/pi,magY3); title(100点DFT);xlabel(w/pi) axis(0,1,0,10),24/62,例 3： 分辨率,增加x(n)有效的样点数，取100个样点 figure(3) n=0:1:99; x=cos(0.48*pi*n)+cos(0.52*pi*n); subplot(2,1,1);stem(n,x); title(signal x(n), 0 = n = 99); xlabel(n) axis(0,100,-2.5,2.5) X=fft(x);magX=abs(X(1:1:51); k=0:1:50;w=2*pi/100*k; subplot(2,1,2);stem(w/pi,magX);title(100点DFT);xlabel(w/pi) axis(0,1,0,60),25/62,4 快速傅里叶变换（FFT）,基本思路 DIT基-2FFT 编程思想 程序设计技术 实序列的DFT,26/62,图基 库利 FFT,引言,“机器计算傅里叶级数的一种算法”,傅里叶 FT,27/62,一、问题的提出,k=0, 1, , N-1,n=0, 1, , N-1,4.1 基本思路,计算量,复乘：,复加：,实乘：,实加：,DFT,都是复数,28/62,二、改善原理 系数Wnk具有以下特性： （1） 对称性,（2） 周期性,4.1 基本思路,（3） 可约性,（4）其它性质,29/62,合并法：合并DFT运算中的某些项。 分解法：将长序列DFT利用对称性和周期性，分解为短序列DFT。,三、改善DFT运算效率的基本途径,4.1 基本思路,快速傅里叶变换正是基于这样思想而发展起来的。 它的算法形式有很多种，但基本上可以分成两大类： 按时间抽取（Decimation in Time，DIT） 按频率抽取（Decimation in Frequency， DIF）,30/62,一、算法原理 设序列x(n)长度为N，且满足N=2M，M为正整数。按n的奇偶把x(n)分解为两个N/2点的子序列：,4.2 按时间抽取的基 -2 FFT,若不满足这个条件，可以人为地加上若干零值（加零补长）使其达到 N=2M,则DFT转化为下式：,31/62,由于,4.2 按时间抽取的基 -2 FFT,故上式可表示成：,可见，一个N点DFT分解成两个N/2点的DFT。 但由于X1(k)，X2(k)只有N/2个点，而X(k)却有N个点，故计算得到的只是X(k)的前一半的结果，要用X1(k)、X2(k)来表达全部的X(k)值，还必须应用系数的周期性。,32/62,得到：,同理：,说明：后半部分k值（N/2kN-1）所对应的X1(k)、2(k)分别等于前半部分k值（0kN/2-1）所对应的X1(k)、X2(k)。,4.2 按时间抽取的基 -2 FFT,即,33/62,我们已知：,因此得：,4.2 按时间抽取的基 -2 FFT,= -1,34/62,按时间抽取将一个N点DFT分解为两个N/2点DFT（N=8）,4.2 按时间抽取的基 -2 FFT,35/62,一个N点DFT分解为两个N/2点DFT，每一个N/2点DFT只需(N/2)2=N2/4次复数乘法，N/2(N/2-1)次复数加法。两个N/2点DFT共需2(N/2)2=N2/2次复数乘法和N(N/2-1)次复数加法。此外，把两个N/2点DFT合成为N点DFT时，有N/2个蝶形运算，还需要N/2次复数乘法及2N/2=N次复数加法。因而通过第一步分解后，总共需要 (N2/2)+(N/2)=N(N+1)/2N2/2 次复数乘法和N(N/2-1)+N=N2/2次复数加法。由此可见，通过这样分解后运算工作量差不多节省了一半。 既然这样分解是有效的，由于N=2M，因而N/2仍是偶数，可以进一步把每个N/2点子序列再按其奇偶部分分解为两个N/4点的子序列。 依次类推，得到下图所示的DFT图：,4.2 按时间抽取的基 -2 FFT,36/62,一个N=8点DFT分解为四个N/4点DFT,4.2 按时间抽取的基 -2 FFT,37/62,N=8 按时间抽取的蝶形运算流图,4.2 按时间抽取的基 -2 FFT,38/62,二、运算量比较 由按时间抽取法FFT的流图可见，当N=2M时，共有M级蝶形，每级都由N/2个蝶形运算组成，每个蝶形需要一次复乘、二次复加，因而每级运算都需N/2次复乘和N次复加，这样M级运算总共需要：,式中，数学符号：lb=log2。,4.2 按时间抽取的基 -2 FFT,39/62,由于计算机上乘法运算所需的时间比加法运算所需的时间多得多，故以乘法为例，直接DFT复数乘法次数是N2，FFT复数乘法次数是(N/2) lbN。 直接计算DFT与FFT算法的计算量之比为,当N=2048时，这一比值为372.4，即直接计算DFT的运算量是FFT运算量的372.4倍。当点数N越大时，FFT的优点更为明显。,4.2 按时间抽取的基 -2 FFT,40/62,如果某通用单片计算机的速度为平均每次复数乘需要4s，每次复数加需要1s，用来计算N=1024点DFT，问直接计算需要多少时间。用FFT计算呢？照这样计算，用FFT进行快速卷积对信号进行处理时，估计可实现实时处理的信号最高频率。,习题 4.1,【解】当N=1024=210时，直接计算DFT的复数乘法运算次数为： N 2 =10241024=1 048 576次 复数加法运算次数为： N（N1）=10241023=1 047 552次 所以，直接计算所用计算时间TD为： TD=41061 048 576+1 047 5521106=5.241 856 s 用FFT计算1024点DFT所需计算时间TF为：,41/62,快速卷积时，需要计算一次N点FFT、N次频域复数乘法和一次N点IFFT。 所以，计算1024点快速卷积的计算时间Tc约为：,习题 4.1,42/62,所以， 每秒钟处理的采样点数（即采样速率）,由采样定理知， 可实时处理的信号最高频率为,习题 4.1,43/62,研读教材 P115的4.2.4节。,一、原位运算（同址运算） 可以看出这种运算是很有规律的，其每级（每列）计算都是由N/2 个蝶形运算构成的，每一个蝶形结构完成下述基本迭代运算：,式中，m表示第m列迭代，k, j为数据所在行数。,4.3 编程思想,44/62,蝶形运算单元,4.3 编程思想,45/62,某一列的任何两个节点k和j的节点变量进行蝶形运算后，得到结果为下一列k, j两节点的节点变量，而和其他节点变量无关。因而可以采用原位运算，即某一列的N个数据送到存储器后，经蝶形运算，其结果为下一列数据，它们以蝶形为单位仍存储在这同一组存储器中，直到最后输出，中间无需其他存储器。也就是蝶形的两个输出值仍放回蝶形的两个输入所在的存储器中。 每列的N/2个蝶形运算全部完成后，再开始下一列的蝶形运算。 这样存储器数据只需N个存储单元。下一级的运算仍采用这种原位方式，只不过进入蝶形结的组合关系有所不同。 这种原位运算结构可以节省存储单元， 降低设备成本。,4.3 编程思想,46/62,二、倒位序规律 观察同址计算结构，发现当运算完成后，FFT的输出X(k)按正常顺序排列在存储单元中，即按X(0)，X(1)，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，而是按x(0)，x(4)， ， x(7)的顺序存入存储单元，看起来好像是“混乱无序”的，实际上是有规律的，我们称之为倒位序。,4.3 编程思想,47/62,造成倒位序的原因是输入x(n)按标号n的偶奇的不断分组。 如果n用二进制数表示为(n2n1n0)2（当N=8=23时，二进制为三位）， 第一次分组，n为偶数（相当于n的二进制数的最低位n0=0）在上半部分，n为奇数（相当于n的二进制数的最低位 n0=1）在下半部分。下一次则根据次最低位n1为“0”或是“1”来分偶奇（而不管原来的子序列是偶序列还是奇序列）。如此继续分下去，直到最后N个长度为1的子序列。下图的树状图描述了这种分成偶数子序列和奇数子序列的过程。,4.3 编程思想,48/62,倒位序的形成,4.3 编程思想,49/62,一般实际运算中，总是先按自然顺序将输入序列存入存储单元，为了得到倒位序的排列，我们通过变址运算来完成。如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示，则其倒位序J对应的二进制数就是（n0n1n2）。这样，在原来自然顺序时应该放x(I)的单元，现在倒位序后应放x(J)。 例如, N=8时，x(3)的标号是I=3，它的二进制数是011，倒位序的二进制数是110，即J=6，所以原来存放在x(011)单元的数据现在应该存放在x(110)内。下表列出了N=8时的自然顺序二进制数以及相应的倒位序二进制数。,4.3 编程思想,50/62,N=8时的自然顺序二进制数和相应的倒位序二进制数,4.3 编程思想,51/62,由表可见，自然顺序数I增加1，是在顺序数的二进制数最低位加1，向左进位。而倒序数J则是在二进制数最高位加1， 逢2向右进位。 例如，在（000）最高位加1，则得（100），再在（100）最高位加1，向右进位，则得（010）。因（100）最高位为1， 所以最高位加1要向次高位进位， 其实质是将最高位变为0，再在次高位加1。用这种算法，可以从当前任一倒序值求得下一个倒序值。,4.3 编程思想,52/62,对于N=2M，M为二进制数最高位，其权值为N/2，且从左向右二进制位的权值依次为N/4，N/8， 2，1。 因此，最高位加1相当于十进制运算J+N/2。如果最高位是0（JN/2）， 则直接由J+N/2得下一个倒序值；如果最高位是1（JN/2），则要将最高位变为0（J J-N/2），次高位加1（J+N/4）。但次高位加1时，同样要判断次高位0、1值，如果为0（JN/4），则直接加1（J J+N/4）; 否则， 将次高位变为0（J J-N/4），再判断下一位；以此类推，直到完成最高位加1，逢2向右进位的运算。,4.3 编程思想,53/62,把按自然顺序存放在存储单元中的数据，换成FFT原位运算流图所要求的倒位序的变址功能如下图所示，当I=J时，不必调换；当IJ时，必须将原来存放数据x(I)的存储单元内调入数据x(J)，而将存放x(J)的存储单元内调入x(I)。为了避免把已调换过的数据再次调换，保证只调换一次（否则又回到原状），我们只需看J是否比I小。若J比I小，则意味着此x(I)在前边已和x(J)互相调换过，不必再调换了; 只有当JI时，才将原存放x(I)及存放x(J)的存储单元内的内容互换。这样就得到输入所需的倒位序列的顺序。,4.3 编程思想,54/62,N=8 倒位序的变址处理,4.3 编程思想,55/62,3. 蝶形运算两节点的“距离” 输入是倒位序的，输出是自然顺序的。其第一级（第一列）每个蝶形的两节点间“距离”为1， 第二级每个蝶形的两节点“距离”为2，第三级每个蝶形的两节点“距离”为4。 由此类推得，对N=2M点FFT，当输入为倒位序，输出为正常顺序时，其第m级运算，每个蝶形的两节点“距离”为2m-1。,4.3 编程思想,56/62,四、WNr 的确定 由于对第m级运算，一个DFT蝶形运算的两节点“距离”为2m-1， 因而:,为了完成上两式运算，还必须知道系数Nr的变换规律： 把式中蝶形运算两节点中的第一个节点标号值，即k值，表示成M位（注意N=2M）二进制数；,4.3 编程思想,57/62, 把此二进制数乘上2M-m，即将此M位二进制数左移M-m位（注意m是第m级运算），把右边空出的位置补零， 此数即为所求r的二进制数。 从图看出，WNr因子最后一列有N/2种，顺序为WN0， WN1， ， 其余可类推。,4.3 编程思想,58/62,4.4 FFT程序设计,程序框图,function A=myFFT_1(x,m) N=2m; if length(x)N % 补0 x=x,zeros(1,N-length(x); end A=myBitRevorder(x,N); % 倒序函数 for L=1:m % L级循环 B=2(L-1); % B个旋转因子 Nmr=2L; u=1; % 旋转因子u初始化 WN=exp(-i*2*pi/Nmr); % 因子WN for j=1:B % B次蝶形运算 for k=j:Nmr:N % 本次碟形运算 t=A(k+B)*u; % 碟形运算 A(k+B)=A(k)-t; A(k)=A(k)+t; end u=u*WN; % 修改旋转因子 end end,59/62,function A=myBitRevorder(A,N) j=N/2; for i=2:1:N-1 if i=k j=j-k; k=k/2; end j=j+k; end,倒序程序,4.4 FFT程序设计,参数 A：源序列 N：点的数目,倒序程序,60/62,4.4 FFT程序设计,MatLab提供的倒序函数： nxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1; y=x(nxd); 其中， fliplr的含义如下： B = fliplr(A) % 前后倒序 例如： 若：A = 1 3 5 7 9 则：B = 9 7 5 3 1 bin2dec：二进制十进制 例如：bin2dec(10111) ，结果： 23 dec2bin：十进制二进制,61/62,4.4 FFT程序设计,测试主程序 clear all; M=8; nn=2M; t=linspace(0,3,nn); % 在0-3之间产生nn个等间隔的行向量 Xn=sin(2*pi*100*t)+sin(2*pi*200*t)+rand(size(t); Yk=myFFT_1(Xn,M); % 调用自编FFT函数 plot(t,Yk);,62/62,4.5 内置FFT函数,Y=fft(X) 如果X是向量，则采用傅立叶变换来求解X的离散傅立叶变换； 如果X是矩阵，则计算该矩阵每一列的离散傅立叶变换； Yfft(X,N) N是进行离散傅立叶变换的X的数据长度，可以通过对X进行补零或截取来实现。 Yfft(X,dim) 或Yfft(X,N,dim) 在参数dim指定的维上进行离散傅立叶变换； 当X为矩阵时，dim用来指定变换的实施方向：dim=1,表明变换按列进行；dim=2表明变换按行进行。 函数ifft的参数应用与函数fft完全相同。,63/62,4.5 内置FFT函数,例如： t=0:0.001:1; %采样周期为0.001s,即采样频率为1000Hz; %产生受噪声污染的正县正弦波信号； x=sin(2*pi*100*t)+sin(2*pi*200*t)+rand(size(t); subplot(2,1,1) plot(x(1:50); %画出时域内的信号； Y=fft(x,512); %对X进行512点的傅立叶变换； f=1000*(0:256)/512; %设置频率轴坐标，1000为采样频率； subplot(2,1,2) plot(f,Y(1:257); %画出