傅立叶变换与频域分析.ppt

第九章傅立叶变换与频域分析,第一节傅立叶变换及其意义（FourierTransform）第二节快速傅立叶变换（FastFourierTransform）第三节傅立叶变换的性质（PropertiesoftheFourierTransform）第四节频域分析（FrequencyDomainanalysis）第五节频域分辨率和谱图表示（FrequencyResolutioninFrequencyDomain）第六节幅值平方相干函数（Magnitude-SquaredCoherentFunction）第七节频域滤波（FilteringinFrequencyDomain）,2,第一节傅立叶变换及其意义（FourierTransform）,9.1.1傅立叶变换的意义及各种变换对,如果一个LTI系统的输入可以表示为周期复指数的线性组合，则输出也一定能表示成这种形式，并且输出线性组合中的加权系数与输入中对应的系数有关。,图9.1,各种信号的傅立叶级数和傅立叶变换对：,傅立叶变换的意义,把一个无论多复杂的输入信号分解成复指数信号的线性组合，那么系统的输出也能通过图9.1的关系表达成相同复指数信号的线性组合，并且在输出中的每一个频率的复指数函数上乘以系统在那个频率的频率响应值。一个域离散必然另外一个域周期，相反的，如果一个域连续必然另外一个域是非周期的。,5,9.1.2离散傅立叶变换（DFT）,离散傅立叶变换的导出有多种方法，比较方便，同时物理意义也比较清晰，是从离散时间傅立叶变换（DTFT）和从离散傅立叶级数（DFS）入手。,【例9-1】试计算常用信号和的N点DFT。,解：,旋转因子具有下列性质：,周期性:,共轭对称性:,可约性:,第二节快速傅立叶变换（FastFourierTransform）,FFT是对计算DFT的快速算法的总称，FFT算法很多，最经典的一种就是库利图基算法，包括基于时间抽选和频率抽选的以二为基底的FFT算法；由以二为基底发展了任意基数的FFT算法。,设序列的长度N2m，其中m为正整数，如果不满足该条件，可以通过补零方法来达到该条件。既然点长为偶数，就先把序列分成两组，偶数项为一组，奇数项为一组，分别用两个序列来表示：(9-1)则N点DFT运算也相应分为两组：,根据的可约性，有，上式变成：,其中分别为的N/2点DFT：,(9-2),(9-3),利用的隐含周期性可以得到另外一半值.,从而得到N点DFT分解计算式：,(9-4),将式（9-4）用信号流图表示，如图9-2，左边表示输入，右边表示输出，支路上的箭头表示乘法运算，乘的因子只对有相位变换而没有幅度变换，所以被称为旋转因子，由于此图像蝴蝶，故称为蝶形运算。一个蝶形运算只包括一次复数乘法、两次复数加法。,图9-2蝶形流图,12,第三节傅立叶变换的性质（PropertiesoftheFourierTransform）,设序列和都是N点长，它们对应的N点DFT分别为和,来讨论傅立叶变换的一些性质。,1.线性,a，b为任意常数。如果两个序列的长度不同则短的序列补零使得两个序列长度相同即可。,13,2.时间翻转特性,证明：,这里需要补充,3.序列的循环移位,因而有,序列的循环移位在第六章详细介绍过，这里简单给出循环移位的定义：,14,即序列的循环移位相当于频域的相移。根据时域和频域的对偶性质，则频域的循环移位对应时域的调制：,上式表示的含义为，先将序列以N为周期进行周期性延拓，得到，然后再进行移位，得到，最后取主值序列，得到仍然是一个N点长的序列。,循环移位后的DFT为:,因此，序列循环移位后的DFT为：,15,4.循环卷积,第六章介绍了循环卷积的计算，这里考虑时域循环卷积结果和频域的关系。设,则有:,通常把式（9-5）称为循环卷积，它的结果仍然是N点长的序列,循环卷积交换序列的先后次序得到的结果都相同。时域和频域的对偶关系，可以得到频域循环卷积对应时域相乘：,(9-5),时域循环卷积对应于DFT的相乘，注意不要和线性卷积混淆，两个序列线性卷积对应于DTFT的相乘：,16,式中表示循环卷积运算符，式中表示线性卷积运算符。循环卷积和线性卷积存在一定关系，由第六章知道，循环卷积是N点循环卷积结果，序列长度为N，线性卷积序列长度为2N1。假设序列是两个序列的L点循环卷积，LN，就需要对补零，然后以L为周期进行周期延拓,则它们的L点循环卷积为：,(9-6),17,【例9-2】设有两序列分别为求它们的线性卷积和5点循环卷积。,式（96）表示循环卷积是线性卷积以L为周期进行周期延拓，然后取L点主值的结果。明显，如果线性卷积就等于循环卷积结果，如果,则循环卷积是线性卷积以L为周期延拓的混叠。,解：线性卷积,直接计算得到6点序列值：,循环卷积,，用表格法来计算，,如表9.2所示。,18,表9.2表格法求循环卷积,我们利用上述结果来验证式（96）是否正确.对线性卷积结果以5为周期进行周期延拓,则有,19,结果和5点循环卷积相同，比较这两个卷积结果，发现只有两点（n0，n5）发生了重叠，其它点结果都相同。,5.共轭对称性,我们知道任意一个信号可以表示成它的奇对称部分和偶对称部分之和，那里的对称是关于坐标原点或者纵坐标的对称性。DFT中的复序列和频域都是在0到N-1的范围内，因而它的对称是在主值范围内的对称，称为周期共轭对称和周期共轭反对称,它们的对称关系如下：,20,严格说上式当n0时有出现，已经超过主值范围，所以一般补充认为N点的值就等于在0点的值。,设任意有限长复序列可以分解成周期共轭对称分量和周期共轭反对称分量之和：,其中易证,当是实数序列时，共轭可以去掉,得：,21,同理，频域序列也可以分解成周期共轭对称分量和周期共额反对称分量之和：,周期共轭对称分量的含义是模数相等，幅角相反，周期共额反对称分量的含义是实部相反，虚部相等。,易证明DFT的共轭对称性可以用下式表示：,22,【例9-3】已知,若是实序列，并且，试证明也是实偶对称的。,证明：由于偶对称，则，由式（97）知，即为实序列。由于是实序列，则，因而，即为偶对称。证毕。,6.帕塞瓦尔（Parseval）定理,帕塞瓦尔（Parseval）定理是序列的能量定理,若则有,(9-7),23,证明：,计算序列的能量可以从时域或者频域入手。,【例2-4】已知，N=6，不求它的DFT结果，来计算下值：,24,（1）,（2）,（3）,解:（1）利用正变换公式，令k0，得,即所有序列值之和。,（2）利用反变换公式，令n0，得,（3）利用帕塞伐尔定理，得,表2.3列出了N点DFT的主要性质。,表2.3DFT的性质表,26,第四节频域分析（FrequencyDomainanalysis）,离散傅立叶变换作为傅立叶变换的一种近似而得到广泛应用，它的快速算法保证了DFT在实时信号处理中的应用。下面介绍频谱分析中常用的几种。,1.幅度谱,N点长序列的DFT结果,是离散的复序列，可以用下式表示：(9-8),离散傅立叶变换的模,表示信号的各复指数信号的频率分量（）的相对大小。例如，在k0附近小范围以外,那么所呈现的仅是相当低的频率。,27,如果序列是实序列，根据例题2-4，则偶对称，即模数相等，幅角相反.这时画出的幅度谱就是偶对称的，往往只需要画一半即可。,画幅度谱时，采用对数坐标也是很常用的，即幅度大小用来代替，这时纵坐标的单位就是分贝（dB），0分贝对应模等于1，20分贝就对应10倍的增益，20分贝对应于衰减0.1，等等。,【例2-6】已知信号,N取一个周期的大小，画出该信号的幅度谱并解释该图。,解：信号的第一个成分的周期为，,28,第二个成分的周期为，,因而的周期为N32。,幅度谱为：,29,因此信号的幅度谱如图2.6所示。,图2.6信号的幅度谱,由幅度谱可以看出信号只在k3，5，27，29有大小，它代表的含义是信号所包含的各个复指数频率分量的大小，即只有四个复指数频率分量存在：,30,k3，29的复指数分量大小是k5，27的复指数分量的一倍。这些和信号的幅度、频率信息相符合，但是没有给出该信号的相位信息。由于幅度谱的偶对称性，往往只画出一半的幅度谱即可。,2.相位谱,表示相位角，它的大小不会影响各复指数频率分量的大小，但能提供这些频率的初始相位信息。对信号的性质有着显著的影响,因此一般包含了信号的大量信息,用相同的幅度谱和不同的相位谱得到的信号完全不同。如果序列是实序列，即相位谱是奇对称。,31,【例2-6】画出例题2-6的相位谱并解释该图。,解：因为,相位谱为：,因此信号的相位谱如图2.7所示，纵坐标表示相位角除以的大小，由相位谱可以看出信号只在k3，5，27，29有值，它代表的含义是信号所包含的各个复指数频率分量的初相位，,32,例如k5表示信号的复指数频率分量为的初相位为.这些和信号的相位信息相符合，但是没有给出该信号的幅度信息。由于相位谱的奇对称性，往往只画出一半即可以得到另外一半的图形。,图2.7信号的相位谱,33,3.功率谱,信号的离散傅利叶变换一般是一个复数，与其共轭之积称为自功率谱，简称自谱或功率谱。其他文献多叫功率谱密度（函数）,其表示为：(2-33),系数1/N是为了满足式能量定理而进行的调整。反映的是信号的功率密度，在图形上与幅度谱类似，只是大小不同，功率谱不含相位信息，所以由功率谱不能恢复原始信号，因为存在多义性。,34,可以证明，线性自相关函数和功率谱是一对离散时间傅立叶变换对（DTFT），相应的循环自相关函数和功率谱是一对离散傅立叶变换对（DFT），考虑序列可能是复数，由于实际得到的是一段样本序列，因而它的自相关利用第三章的定义式（3-21），因而：,35,功率谱是正、实序列。当是实序列时，自相关序列也是实序列，则功率谱是偶对称的。从功率谱和自相关函数之间的关系，我们知道功率谱蕴涵着集合统计的实质，一个随机信号的自相关和功率谱都表达了随机信号的统计平均特性。,同样，可以给出两个信号之间的互功率谱的定义：,【例2-8】画出例题2-6的功率谱并解释该图。,36,解：由例题2-6知,因此功率谱如图2.8所示，该图表示信号各复指数频率分量的功率，波形与幅度谱类似，大小略有不同，也没有包含相位信息。,37,图2.8信号的功率谱,一般信号都是包含多种频率分量的，数据点也比较长，这时可以利用计算机的快速傅立叶变换来计算DFT，然后画出该信号的幅度谱、相位谱、功率谱进行分析。图2.9是心电信号及其功率谱。,38,图2.9（a）心电信号（b）功率谱,第五节频域分辨率和谱图表示,（FrequencyResolutioninFrequencyDomain）,利用离散傅立叶变换进行频谱分析时会引起误差，常见的就是由于离散傅立叶变换本身采样和截断过程所引起的混叠、泄漏、栅栏现象等。如果要分析的信号是连续信号，就必须先采样.当采样频谱比信号最高频率的两倍要小,发生混叠现象,可以提高采样率来避免混叠现象。如果要分析的信号是周期连续信号，就必须要对该信号截取一段来进行分析，即加了一个窗，便会发生泄漏现象，减少泄漏可以通过加不同的窗函数来截取信号.,40,离散傅立叶变换是对离散时间傅立叶变换的采样，它只给出频谱在离散点（）上的值，而无法反映这些点之间的频谱内容，这就是栅栏现象.改善栅栏效应的一种方法是信号后面补若干个零，通过补零来调整坐标上的位置。,在时域中，数字信号的时间分辨率由采样间隔t决定。在频域中，频域分辨率（两相邻谱线间的频率差值）由采样频率和采样点数N决定：,41,式中由采样定理决定，因此，要提高频域分辨率，就得增加采样点数N。如果数据量略有不足，传统方法是在数据尾部填0来解决，称为高密度频谱（thehighdensityspectrum）。,一般说来，功率谱、幅度谱、位相谱等皆为离散值而非连续值，把上述谱绘成连续谱线图是不恰当的，是误解，在由频率分辨率决定的两条谱线间的值是不能确定的或不存在的，如果将谱线绘成连续曲线则会误解为两条谱线间的值是确定的。,【例2-9】高密度频谱和高分辨率频谱的比较。设信号为：,利用有限长序列的FFT来分析下列情况的幅度谱,（a）采集数据长度N10，即，做10点的DFT，画出幅度谱。,（b）采集数据长度10点，但补90个零，做100点的DFT，画出幅度谱。,42,（c）采集数据长度100点，同样画出幅度谱。,解：利用FFT编程实现上述三种情况下的幅度谱，分别计算10点FFT、100点FFT，为了比较频率，把坐标k通过转化成频率的坐标，由于幅度谱的偶对称只画出一半的DFT结果即可，因此的范围是从0，把它归一化，则最后结果如图2.11所示。,图2.11的第一列是三种信号，第二列是对应信号的DFT幅度谱。第一行是采集10点，然后做DFT，由于取的点数太少，分辨率太低：，从幅度谱中无法确定该信号的频率分布情况。第二行是补90个零后的幅度谱，即高密度频谱，从图中看成最大成分是,这个结果也和原信号包含两个频率成分不相符合。前面两种情况都发生了泄漏现象。第三行就采集了100个数据，由于有足够的数据，幅度谱清楚的反映了原信号包含的两个频率成分,这就是高分辨率频率,能够分辩靠得很近的频率成分。我们可以计算它的分辨率达到：。,43,图2.11三种情况下的幅度谱图,44,通过上述讨论，利用FFT做谱分析时各参数的选择为：,1.采样频率应满足采样定理：,(2-36),或,(2-37),往往一段信号的频谱是无限大的，因而要先通过一个带通滤波器，来限制信号的最高频率。,2.离散傅立叶变换的N为：,一般，为了FFT计算快速，N都尽量取成2的幂次。,(2-38),45,3.采集信号的持续时间为：,(2-39),【例2-10】已知一连续信号为,试用DFT计算其幅度谱，并且与原信号进行比较。,解：由于是周期信号，首先要确定该信号的周期，如果截断信号时没有完整的一个周期的信号则会发生泄漏现象。第一部分的周期为0.8ms，第二部分的周期为1ms，第三部分周期为1.25ms，整个信号的周期为20ms，因此截取信号的持续时间最好为20ms的整数倍。信号的最高频率为，因此采样频率或采样间隔,46,取采样频率和采样间隔分别为：,频率分辨率为,很容易满足该条件，因此DFT的点长N从持续时间中算得：,令，代入原信号，利用FFT计算100点的DFT，幅度谱如图2.12所示。或者为了FFT计算，令N为2的幂次，但是要保证NT是20ms的整数倍,比如N取128，T取0.1563ms。只看前半图，从图2.12可以看出信号的频率成分包含k16、20、25，即包含三个数字频率成分,47,转换成模拟频率为,或者利用等式：,来求相应的Hz为单位的模拟频率。,幅度谱的大小和频率和原信号完全一致。如果N取得不合适就会发生泄漏现象，整个频率上都会有一定幅度大小。因此要注意取法，当然，正确的采样频率和N的组合是多种的。,图2.12信号的幅度谱,48,第六节幅值平方相干函数,（Magnitude-SquaredCoherentFunction）,频域相干函数，也称为幅值平方相干函数，设有两个信号，它们的幅值平方相干函数定义如下：,(2-40),其中表示两个信号的互功率谱,即这两个信号的循环互相关函数的离散傅立叶变换.为各自的功率谱,这里的k代表含义即频率。,的取值范围为01之间.,=1,说明两信号是完全相干的，即一个信号可以完全由另外一个信号决定；,=0,这两个信号不相干,即这两个信号是完全独立的；,在(0,1),说明这两个信号存在部分相干性,即非线性关系或者有外界的干扰存在。,49,可见，频域相干函数可以从频域上表示两个信号各频率成分互相关联的程度。,【例2-11】设有两个信号和,测量这两个信号时，假设含有了不相关的噪声和,即测量得到的两个信号为和,比较理想信号和测量信号的幅值相干函数。,解：设和的功率谱分别为,互功率谱为利用式(2-40)，理想的相干函数为,设噪声的功率谱为,由于信号和噪声不相关，则噪声和信号的互相关函数为零,因此测量信号的功率谱和互功率谱为：,50,因此测量信号的幅值相干函数：,(2-41),由于功率谱是非负的，因此上式满足,例2-11就是在有不相关噪声干扰下，频域相干函数将小于1。一般我们能得到的信号都是测量信号而非理想信号，因此要从测量信号的相干曲线来判断理想信号各频率成分互相关联的程度。,下面用实际算例来看相干曲线。,51,【例2-12】设信号,假设观测时引入的噪声均为白噪声，它们的功率谱密度都为1；观测记录N10点，比较理想信号和观测值的幅值相干函数。,解：先利用FFT求出和的功率谱，再利用式（2-40）、（2-41）计算相干函数，MATLAB程序如下：,N=10;n=0:N-1;x=N-1:-1:0;y=0:N-1;px=abs(fft(x).2/N;py=abs(fft(y).2/N;pxy=abs(conj(fft(x).*fft(y)/N).2;rxy=pxy./(px.*py);r=1./(1+1./px+1./py+1./(px.*py);figure(1);subplot(1,2,1);stem(n,px);xlabel(k);title(x(n)功率谱)subplot(1,2,2);stem(n,py);xlabel(k);title(y(n)功率谱)figure(2);subplot(1,2,1);stem(n,rxy);xlabel(k);title(理想相干函数)subplot(1,2,2);stem(n,r);xlabel(k);title(有白噪干扰的相干函数),52,图2.13x和y信号的功率谱图,图2.14理想的和有干扰的频域相干函数,53,结果如图2.13和2.14所示，2.13给出了和两个信号的功率谱图，两者完全相同，但要注意功率谱无法给出相位信息,事实上它们的相位谱是完全不同的。,第七节频域滤波,（FilteringinFrequencyDomain）,理想的数字滤波器幅度谱如图2.15所示，这些频率特性都是以2为周期的连续函数，当单位脉冲响应是实数序列时，幅度谱周期偶对称，相位谱周期奇对称，因而只需要给出一半的频谱图即可。,54,图2.15各类理想数字滤波器频率特性,55,数字频域滤波可以用硬件和软件的方法实现，滤波器的设计方法多种多样，大致分为IIR滤波器的设计和FIR滤波器设计，前者主要利用传统的模拟滤波器设计方法，后者多采用窗函数和频率取样设计法,频域滤波用软件方法实现,更灵活。,如不考虑信号因果性（非实时分析），则滤波特性十分陡直，实现方法也很简单，只需令欲滤波的频段对应的幅度为0，再作IFFT即可获得滤波后的时域波形，这对干净的滤去50Hz工频干扰是十分有效的。,如果要实现实时分析，就要设计出因果的滤波器，在Matlab中有滤波器设计工具箱，可以参考。,56,图2.16是房颤波的频域滤波例子。a)是房颤波，b)是滤去大于4Hz的成分（低通滤波器）而剩下相对的低频成分经IFFT后的时域图形，c)是保留412Hz的成分（带通滤波器）经IFFT后的时域图形，d)是滤去小于12Hz的成分（高通滤波器）而剩下相对的高频成分经IFFT后的时域图形。,图2.16房颤波的频域滤波,谢谢!,