完结篇—数字信号时频分析研究与分析软件设计.doc

67附录学号代号0713321602 10701密级分类号 TP391 公 开数字信号时频分析研究与分析软件设计题（中、英文）目Research on Time-Frequency Analysis of Digital Signalsand Design of the Software Applied to Analyze Signals宣宗强 高工指导教师姓名、职称 苏 鹏 博作者姓名测试计量技术及仪器学科、专业 工学学科门类二一一年一月提交论文日期西安电子科技大学学位论文独创性（或创新性）声明秉承学校严谨的学风和优良的科学道德，本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加以标注和致谢中所罗列的内容以外，论文中不包含其他人已经发表或撰写过的研究成果；也不包含为获得西安电子科技大学或其它教育机构的学位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中做了明确的说明并表示了谢意。申请学位论文与资料若有不实之处，本人承担一切的法律责任。本人签名： 日期 西安电子科技大学关于论文使用授权的说明本人完全了解西安电子科技大学有关保留和使用学位论文的规定，即：研究生在校攻读学位期间论文工作的知识产权单位属西安电子科技大学。学校有权保留送交论文的复印件，允许查阅和借阅论文；学校可以公布论文的全部或部分内容，可以允许采用影印、缩印或其它复制手段保存论文。同时本人保证，毕业后结合学位论文研究课题再攥写的文章一律署名单位为西安电子科技大学。（保密的论文在解密后遵守此规定）本学位论文属于保密，在 年解密后适用本授权书。本人签名： 日期 导师签名： 日期 摘 要时频分析作为分析非平稳时变信号的有力工具，成为现代信号处理研究的一个热点，近年来受到越来越多的重视。时频分析方法提供了时间域与频率域的联合分布信息，清楚地描述了信号频率随时间变化的关系。本文对时频分析的三种算法进行了原理分析和研究，通过Matlab得出三种算法的特点和他们之间的区别，并用C语言开发了算法的实现代码。基于时频分析算法，本文以Visual Studio 2005为工具，开发设计了数字信号的频谱分析软件。分析软件对信号进行数字下变频(DDC)、抽取和信号在时域、频域和调制域的分析，并对分析的结果进行图形显示，以实现软件在Windows操作系统下对数字信号的测量和分析。与扫描频谱分析仪相比，该软件可以对信号进行块模式的无缝采集，它无缝采集每个块，并存储在内存中，然后进行信号的后期处理。这就克服了在信号变化迅速时，扫描频谱可能会漏掉信号的高频部分的缺点。该软件在数据采集系统中的测试表明，它测量精度和数据分析能力都达到了立项的设计要求，具有较高的实用价值。关键词：时频算法 软件设计 频谱分析仪 AbstractAs a powerful tool for analyzing non-stationary signals, Time-Frequency Analysis（TFA） has been a hot spot of modern signal research in recent years and got more and more attentions. Time-frequency analysis method provides the joint distribution of a time domain and frequency domain. It clearly describes that how the frequency of signal varies with time.In this paper, the theories of three time-frequency algorithms have been analyzed and researched. The features of three algorithms and differences between them have been obtained by Matlab software, and C code of three algorithms developed.Based on time-frequency algorithms, we have designed software of analyzing digital signal spectrum with Visual Studio 2005 as a development tool.The software processes the signal through the mode of Digital Data Converter and of extracting and then analyzes the signal in time domain, frequency domain and modulation domain; the results are displayed by graphical style.The software realizes the spectrum measurement and analysis of digital signal.Compared with the Scanning Spectrum Analyzer（SPA）, the software can get signals seamlessly in block mode. Signals are got and stored in memory, and then processed. It overcomes the problem that the high-frequency signal may be missed in SPA when the signals have a section of rapid changes.The tests in the data acquisition system show that precision and data analysis capabilities of the software have reached the projects design requirements, with a high practical value.Keyword: Time-frequency algorithms Design software Spectrum analyzer 目录第一章绪论11.1课题的背景和意义11.2数字信号时频分析的现状及发展11.3本文内容和结构41.3.1 数字下变频（DDC）模块61.3.2 信号时频分析算法的选择71.3.3 数字信号时频分析软件的开发8第二章数字下变频92.1正交数字下变频原理92.2滤波抽取原理简述102.2.1 FIR低通滤波器原理102.2.2 FIR滤波器特点112.2.3 FIR数字滤波器的窗函数设计112.3滤波器系数的计算12第三章时频分析方法的理论与实现153.1 STFT算法理论与实现153.1.1 傅里叶变换153.1.2 短时傅里叶变换163.1.3 窗函数以及窗函数宽度对时间-频率分辨率的影响193.1.4 离散信号的短时傅立叶变换及其实现223.2小波变换算法223.2.1 连续小波的定义233.2.2 Morlet复值小波243.2.3 连续Morlet复值小波变换的计算253.2.4 连续小波变换的性质263.3 Wigner分布293.3.1 Wigner分布的定义293.3.2 WVD的性质303.3.3 WVD的缺点323.3.4 Wigner 分布的实现323.3.5 Wigner分布中交叉项的的行为333.4 MATLAB对三种算法的验证333.4.1 基于Matlab的算法验证333.4.2 三种时频分析算法的区别以及时频分析算法的优点37第四章时频分析软件的结构和设计394.1时频分析软件整体结构设计图：394.2信号处理过程过涉及的概念：样点、帧和块404.3软件的关键模块介绍：414.4程序运行的主要介绍和操作功能介绍434.4.1 菜单功能介绍434.4.2 测试软件中参数设置方法和频谱显示图44第五章总结与展望515.1论文工作的总结515.2后续工作的展望525.3时频分析技术的应用前景52致谢55参考文献57研究成果59附录61第一章 绪论1.1课题的背景和意义时间和频率是描述信号的两个重要的物理量，并且信号的时域和频域之间有紧密的联系。信号处理的任务之一是认识客观世界中存在的信号的本质特征，并找出规律。从不同的角度去认识、分析信号有助于了解信号的本质特征。信号最初是以时间（空间）的形式来表达的。除了时间以外，频率是一种表示信号特征最重要的方式。频率的表示方法是建立在傅里叶分析（Fourier Analysis）基础之上的。时频分析（JTFA）即时频联合域分析（Joint Time-Frequency Analysis）的简称。时频分析方法提供了时间域与频率域的联合分布信息，清楚地描述了信号频率随时间变化的关系。时频分析的基本思想是：设计时间和频率的联合函数，用它同时描述信号在不同时间和频率的能量密度或强度。时间和频率的这种联合函数简称为时频分布。利用时频分布来分析信号，能给出各个时刻的瞬时频率及其幅值，并且能够进行时频滤波和时变信号研究。由于傅里叶分析是一种全局的变换，要么完全在时间域，要么完全在频率域，因此无法表述信号的时频局部性质，而时频局部性质恰好是非平稳信号最基本和最关键的性质。为了分析和处理非平稳信号，在傅里叶分析理论基础上，提出并发展了一系列新的信号分析理论：短时傅里叶变换（Short Time Fourier Transform）、连续小波变换、Wigner-Ville分布等。随着现代测量技术的发展，对数字信号进行硬件采集、软件分析的现代测量仪器也飞速革新。将时频分析算法应用于频谱分析仪就可以实现信号的时频联合分析。1.2数字信号时频分析的现状及发展频谱分析仪的主要功能是在频域里显示输入信号的频谱特性，是对无线电信号进行测量的必备手段，是从事电子产品研发、生产、检验的常用工具。因此，频谱分析仪的应用领域相当广泛，诸如卫星接收系统、无线电通信系统、行动电话系统基地台辐射场强的测量、电磁干扰等高频信号的侦测与分析，同时也是研究信号成份、信号失真度、信号衰减量、电子组件增益等特性的主要仪器。频谱分析仪经历了扫频频谱分析仪、矢量信号分析仪和实时频谱分析仪，由于矢量信号分析仪与实时频谱分析仪在结构上有很多相似之处，我们一下只介绍第一种和第三种。扫频频谱分析仪：传统频域分析扫描调试的超外差频谱仪是几十年前第一个使工程师能够进行频域测量的传统结构。扫频频谱分析仪最初是使用纯模拟器件构建的，之后一直随着其服务的应用不断发展。当前一代扫频频谱分析仪包括各种数字单元，如ADC、DSP 和微处理器。但是，基本扫描方法在很大程度上保持不变，其最适合观测受控的静态信号。扫频频谱分析仪通过把感兴趣的信号向下变频，并扫描通过分辨率带宽(RBW)滤波器的传输频带，来测量功率随频率变化。RBW 滤波器后面跟有一个检测器，检测器计算选择的跨度内每个频率点的幅度。尽管这种方法可以提供很高的动态范围，但其缺点在于，它一次只能计算一个频率点的幅度数据。分析仪在频率跨度内扫描需要很长的时间，在某些情况下要达到几十秒。这种方法基于这样一个假设，即分析仪能够完成多次扫描，而被测信号没有明显变化。结果，这种方法要求输入信号相对稳定及不变。图1.1：扫频分析仪漏掉扫描波段之外发生的瞬时事件如果信号迅速变化，那么在统计上可能会漏掉变化。如图1.1 所示，扫描查看频段Fa，而在Fb(左图)上发生了一个瞬时频谱事件(如图中用黄色突出显示的部分)。在扫描到达频段Fb时，事件已经消失，没有检测到事件(右图)。频谱分析仪没有提供触发这个瞬时信号的方式，也不能存储全面的信号在不同时间上的行为记录。图1.2：典型的扫频分析仪结构图1.2说明了典型的现代扫描频谱分析仪结构。它采用数字技术替换较窄的滤波器，而不用继承下来的模拟分辨率带宽(RBW)滤波器。对BW1、BW2或BW3范围内的带宽，ADC之前的滤波、混频和放大都是模拟处理。在需要窄于“BW3”的滤波器时，将在模数转换之后的步骤中通过数字信号处理(DSP)进行滤波。实时频谱分析仪（RTSA）：实时频谱分析的基本概念是能够触发RF信号，把信号无缝地捕获到内存中，并在多个域中分析信号。这可望可靠地检测和检定随时间变化的RF信号。图1.3：典型的实时频谱分析仪结构图1.3 说明了RTSA 结构简化的方框图。可以在仪器的整个频率范围内调谐RF 前端，它把输入信号下变频到固定IF，固定IF与RTSA的最大实时带宽有关。然后ADC 对信号进行滤波、数字化，然后传送到DSP引擎上，DSP引擎负责管理仪器的触发、内存和分析功能。RTSA提供了实时触发、无缝信号捕获和时间相关多域分析功能。此外，ADC的技术进步可以实现高动态范围和低噪声转换，使RTSA的性能相当于或超过许多扫频频谱分析仪的基本RF性能。对小于等于实时带宽的测量跨度，通过数字化RF信号，并把时间连续的样点存储在内存中，RTSA结构能够无缝地捕获没有时间间隔的输入信号。这较扫频频谱分析仪的采集过程有多种优势，因为扫频频谱分析仪通过连续扫描频率跨度，构建频域图像。泰克公司的频谱分析仪RSA6000 系列为这种类型的代表。实时频谱分析系统在硬件上要具备信号采集的能力，采集到的中频信号首先进行调整、A/D采样，然后对采样得到的跨度内的时域数据进行快速傅立叶变换（FFT），得到频域数据，或者进行时频分析，得到信号的时频特性数据。为达到实时测量的目的，实时频谱分析借鉴软件无线电的设计思想，最大限度地利用软件来分析各种信号，实现频谱分析的各项功能。实时频谱分析系统硬件结构将宽带A/D变换尽可能地靠近天线，即尽可能早地将接收到的模拟信号数字化。硬件增加ADC模块，计算机控制部分使用微软的操作系统。信号分析如频谱分析、参数测量、时间测量、分析等功能由软件实现。实时频谱分析系统的软件是基于实时操作系统的控制和分析软件，是整个实时频谱分析的最顶层软件。系统软件除对整个实时频谱分析系统进行控制调度以外，还要提供全面的时域、频域和调制域的信号分析功能。如：频谱分析、参数测量、时间测量、脉冲信号测试等。1.3本文内容和结构信号采集系统的硬件结构与实时频谱分析仪相类似，我们将对采集的数据在分析软件中处理，得出信号的基本时域、频域特点，重点将应用时频分析算法分析信号的时频联合分布。本系统中的输入信号是有限带宽信号,其频率范围为（275MHz475MHz），在采集卡中以500MHz的频率对其进行采样。数据采集及分析系统设计的框架如图1.4所示。图1.4 数据采集及分析系统总体结构本文主要讨论对采集信号进行处理中使用的主要数字信号处理模块，包括：PCI9054接口卡的驱动设计、数字下变频（DDC）、信号时频分析算法的选择和实现以及开发信号频谱测试应用软件。其中，500MHz是指每一路的采样频率。硬件中总共有四路500MHz的采样频率分时对信号进行采集，所以就相当于采集频率为2GHz。1.3.1 数字下变频（DDC）模块数字变频技术是软件无线电的核心技术之一。与模拟变频器相比，数字变频不存在模拟变频器中混频器的非线性和模拟本地振荡器的频率稳定度、边带、相位噪声、温度漂移、转换速率等人们关心但是难以彻底解决的问题，而且数字变频中频率步进和频率间隔也具有理想的性能，并且数字变频器的控制和修改比较容易，实现比较简单。数字变频器由数字混频器、数字控制振荡器（NCO）和低通滤波器三部分组成，如图1.5所示。图1.5 数字下变频结构方框图数字下变频的基本功能是从输入的宽带高数据流的数字信号中提取所需的窄带信号，将其下变频为数字基带信号，并转换成较低的数据流。图1.5是数字下变频的基本模型。图中高速ADC转换器的输出信号送入数字下变频器；经两个相乘器所构成的混频器，将输入来的数字信号和复正弦信号产生器产生的正交正弦信号相乘，相乘结果为I 、Q 两路信号；再分别经高抽取滤波器和有限长冲击响应(FIR) 滤波器的处理，由这两个滤波器构成的复合滤波器的功能是低通滤波和抽取，其输出是数据流降低了的数字基带信号。本文中，对采样数据进行数字下变频，输出的混频频率为375MHz，这时信号的频谱搬移到基频（-100MHz100MHz）。有利于数据的进一步处理，抽取后可以进一步降低要处理的数据量，提高运算速度。1.3.2 信号时频分析算法的选择对DDC模块得出的I、Q数据，选择不同的算法进行分析，得出采样数据的时频分布特性。短时傅里叶变换是一种单一分辨率的信号分析方法，它的思想是：选择一个时频局部化的窗函数，假定分析窗函数g(t)在一个短时间间隔内是平稳（伪平稳）的，移动窗函数，使f(t)g(t)在不同的有限时间宽度内是平稳信号，从而计算出各个不同时刻的功率谱。短时傅里叶变换使用一个固定的窗函数，窗函数一旦确定了以后，其形状就不再发生改变，短时傅里叶变换的分辨率也就确定了。如果要改变分辨率，则需要重新选择窗函数。短时傅里叶变换用来分析分段平稳信号或者近似平稳信号犹可，但是对于非平稳信号，当信号变化剧烈时，要求窗函数有较高的时间分辨率；而波形变化比较平缓的时刻，主要是低频信号，则要求窗函数有较高的频率分辨率。短时傅里叶变换不能兼顾频率与时间分辨率的需求。短时傅里叶变换窗函数受到W.Heisenberg不确定准则的限制，时频窗的面积不小于2。这也就从另一个侧面说明了短时傅里叶变换窗函数的时间与频率分辨率不能同时达到最优。小波变换使用一个窗函数（小波函数），时频窗面积不变，但形状可改变。小波函数根据需要调整时间与频率分辨率，具有多分辨分析（Multi-resolution Analysis）的特点，克服了短时傅里叶变换分析非平稳信号单一分辨率的困难。小波变换是一种时间-尺度分析方法，而且在时间、尺度（频率）两域都具有表征信号局部特征的能力，在低频部分具有较高的频率分辨率和较低的时间分辨率，在高频部分具有较高的时间分辨率和较低的频率分辨率，很适合于探测正常信号中夹带的瞬间反常现象并展示其成分。所以，小波变换被称为分析信号的显微镜。小波变换不会“一叶障目，不见泰山”，又可以做到“管中窥豹，略见一斑”。但是小波分析不能完全取代傅里叶分析，小波分析是傅里叶分析的发展，而时频分析是一种非线性二次变换，与线性的小波变换相去甚远。Wigner分布是具有双线性形式的时频分布。所谓双线性形式，是指所研究的信号在时频分布的数学表达式中以相乘的形式出现两次，又称为非线性时频分布。Wigner-Ville分布有一系列好的性质，这是它得到广泛应用的主要原因，但是也存在严重的交叉项。本论文旨在在实际的项目实践中，根据信号的不同特征，采用不同的算法，对信号进行分析，以期得到最精确的测量结果。同时理解各种不同算法的区别和特性。1.3.3 数字信号时频分析软件的开发在基于windows操作系统下，本文以Visual Studio 2005为工具，开发设计了数字信号的频谱分析软件。使用Visual Studio 2005在MFC框架下生成单文档程序，添加数据接收和保存模块，数据处理模块，图形显示模块，菜单功能模块等。使得软件具有接收、保存数据、分析继而显示信号在时域、频域和调制域的特性。通过反复调试，我们分析获得了信号的幅度特性、瞬时频率特性、瞬时频谱特性和时频特性。可以选择不同的时频分析算法来分析信号的时频特性。论文将在下面的四章对这四部分详细介绍和分析。第二章 数字下变频2.1正交数字下变频原理数字下变频的基本功能是从输入的高宽带数据流的数字信号中提取所需的窄带信号，将其下变频为数字基带信号，并转换成较低的数据流。高速转换器的输出信号送入数字下变频器；经两个相乘器所构成的混频器，将输入来的数字信号和复正弦信号产生器产生的正交正弦信号相乘，相乘结果为I、Q两路信号；再分别经高抽取滤波器和有限长冲击响应(FIR)滤波器的处理，由这两个滤波器构成的复合滤波器的功能是低通滤波和抽取，其输出是数据流降低了的数字基带信号。本文中，源信号以取样速率Fs=500MHz 进行数字转换。然后把数字转换后的源信号发送到DDC。DDC中的数字振荡器在感兴趣的波段中心频率上生成正弦和余弦。正弦和余弦以数字方式乘以数字化源信号，生成I 和Q 基带样点流，其中包含原始信号中包含的所有信息。然后I和Q流经过可变带宽低通滤波器。正交数字下变频包括两个部分：一是乘法器，二是数控振荡器（NCO）。正交数字下变频是将数字化后的实信号分为两路，一路乘以，下变频到0中频，形成与原始信号相位相同的信号；另一路乘以，下变频到0中频，形成与原信号正交的信号。其数学表达式为： 式（2.1）其中为中频信号的载频（中心频率）为375MHz，是采样间隔为1/500MHz，表示ADC后输出实信号，为数字下变频后输出。正交数字下变频中的正弦波和余弦波由NCO产生。NCO主要有三部分组成：相位累加器、相位加法器和sin/cos表只读存储器。相位累加器的作用是将数字本振频率转换成相位，相位加法器的功能是设置一定初始相位，相位的正余弦值由查表求得。此时、两路数据可以表达为式（2.2）和式（2.3）： 式（2.2） 式（2.3）根据公式式（2.2）和式（2.3）即可实现信号的数字下变频中的混频部分。2.2滤波抽取原理简述经过正交数字下变频之后，得到了零中频的基带复信号，此时信号的采样率仍然是ADC采样频率（即中频信号的采样频率500MHz），数据率很高，信号的带宽远小于采样频率，所以滤波抽取（采样率降低）成为了数字下变频的一个必不可少的组成部分。抽取就是从信LPFx(n)y(m)号中除去高频信息，降低抽样频率而不导致频谱混叠的过程。图2.1 抽取电路方框图抽取的基本过程如图2.1所示。先进行低通滤波，再作抽取。FIR数字滤波器也是在DDC中经常用到抽取滤波器，它的优点是FIR滤波器的参数可设置，达到比较好的阻带衰减特性和通带特性。2.2.1 FIR低通滤波器原理所谓数字滤波器是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分的相对比例，或者滤除某些频率成份的器件模型。 FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。2.2.2 FIR滤波器特点有限长单位冲激响应（FIR）滤波器有以下特点：1、既具有严格的线性相位，又具有任意的幅度；2、FIR滤波器的单位抽样响应是有限长的，因而滤波器性能稳定；3、只要经过一定的延时，任何非因果有限长序列都能变成因果的有限长序列，因而能用因果系统来实现；4、FIR滤波器由于单位冲击响应是有限长的，因而可用快速傅里叶变换(FFT)算法来实现过滤信号，可大大提高运算效率。5、FIR也有利于对数字信号的处理，便于编程，用于计算的时延也小，这对实时的信号处理很重要。6、FIR滤波器比较大的缺点就是阶次相对于IIR滤波器来说要大很多。FIR数字滤波器是一个线性时不变系统（LTI），N阶因果有限冲激响应滤波器可以用传输函数H（z）来描述， 式（2.4）在时域中，上述有限冲激响应滤波器的输入输出关系如下： 式（2.5）其中，xn和yn分别是输入和输出序列。2.2.3 FIR数字滤波器的窗函数设计FIR数字滤波器的设计是选择有限长度(长度为H)的单位脉冲响应，使其传输函数满足技术要求。FIR数字滤波器的设计问题就是要求所设计的FIR数字滤波器的频率响应去逼近所要求的理想滤波器的响应。从单位取样序列来看，就是使所设计的滤波器的逼近单位取样响应序列。理想低通滤波器的单位取样响应是无限长的，n从-+，且是非因果序列。的波形如图2.2。因此为了构造一个有限长度为N的线性滤波器，只有将截取一段，使之关于a=（N1）/2对称。对的截取实际上是对其做加窗处理，即让其与一个窗函数相乘，窗函数的形状有矩形窗、海明窗和布莱克曼窗等。这里以矩形窗为例，设矩形窗的长度为N，如图2.3示。对 加窗处理后，得，的波形如图2.4示，这样就用一个有限长序列去代替。图2.2 矩形函数窗 图2.3 理想低通滤波器冲击响应图图2.4 低通滤波器冲击响应图2.3滤波器系数的计算基于MATLAB信号处理工具箱函数可以方便、快捷地设计数字滤波器，在设计过程中，根据实际滤波器特性，可以随时更改参数，以达到设计要求9。本文基于Matlab语言，对采用窗函数法设计的FIR 低通数字滤波器进行了仿真，得出滤波器的系数。滤波器设计相关的Matlab函数fir1()、freqz()介绍：fir1 函数实现线性相位FIR数字滤波器的窗函数法设计。利用这一函数可以设计出标准的低通、高通、带通和带阻线性相位FIR滤波器。调用格式：b=fir1(n,wn, type, window)wn是截止频率，以“弧度”为单位。fir1函数返回FIR滤波器的系数向量b。如果wn是一个标量，则返回的是一个n阶的低通FIR数字滤波器。type为滤波器的类型，当type=high 时，设计高通FIR滤波器；当type=stop 时，设计带阻FIR滤波器。向量window 用来指定窗函数类型，其长度为n+1，缺省时为海明窗。为了检验所设计的数字滤波器是否正确，可以画出其幅频特性。MATLAB提供了求解数字滤波器频率响应的函数freqz，其调用格式如下：h,w=freqz(b,a,n)返回数字滤波器的n点频率响应h和频率向量w。若n缺省，则自动取n为默认值512。设计一个线性相位 FIR 低通滤波器，技术指标如下：通带截止频率fp=100MHz，阻带起始频率fst=140MHz。通带允许的最大衰减为Rp=0.35dB，阻带应达到的最小衰减为As=50dB. 滤波器的采样频率为fs=500MHz。首先根据阻带衰减As=50dB 来选择窗形状，海明窗和布拉克曼窗等窗函数均可提供大于50dB 的衰减。由于海明窗可提供较小的过渡带宽，所以选择海明窗。然后，编写程序，求出技术指标中fp、fst相对应的数字频率，求出过渡带宽，由过渡带宽确定窗口长度N，求出低通滤波器的截止频率，由fir1函数求出滤波器的系数。用Matlab计算滤波器系数，求得滤波器的Rp=0.0324dB，As=54.5162dB，满足设计要求。并验证Rp 和As 是否满足设计要求。画出幅频特性及归一化分贝幅特性如图2.5所示。可以看出，设计是符合要求的。图2.5 幅频特性 归一化分贝幅频特性根据MATLAB工具可以求出此FIR低通滤波器的系数，分别为：H= 0.00211382494510535 -0.000517781080573372-0.004279317569740170.00322741368479888 0.0103944368631471 -0.0121778889206461 -0.01957329553390050.03387884259277840.0294627902974621-0.0851974534051769-0.0370698757123113 0.310255075044648 0.5389664575888170.310255075044648 -0.0370698757123113 -0.08519745340517690.02946279029746210.0338788425927784-0.0195732955339005-0.01217788892064610.01039443686314710.00322741368479888 -0.00427931756974017 -0.000517781080573372 0.00211382494510535在算法中，我们将使用这组FIR低通滤波器的系数来设置滤波器。第三章 时频分析方法的理论与实现3.1 STFT算法理论与实现3.1.1 傅里叶变换傅里叶（Fourier）变换从本质上看无非是研究如何利用简单、初等的函数近似表达复杂函数（信号）的方法和手段5。1807年，法国科学家傅里叶提出周期为2的函数f(x)可以表示为系列三角函数之和，即 式（3.1.1）其中，表达式（3.1.1）可以理解为信号f(x)是由正弦波（含余弦与正弦函数）叠加而成，其中ak，bk为叠加的权值，表示信号在不同频率时刻的谱幅值大小。显然，当信号具有对称性（偶）特征时， 式（3.1.2）而当信号具有反对称性（奇）特征时， 式（3.1.3）在研究热传导方程的过程中，为了简化原问题，傅里叶建议将热导方程从时间域变换到频率域，为此他提出了著名的傅里叶变换的概念。信号f(x)的傅里叶变换定义为： 式（3.1.4）傅里叶变换建立了信号时域与频域之间的关系，频率是信号的物理本质之一。随着计算机技术的发展与完善，科学与工程中的所有计算问题跟计算机已经密不可分，计算机计算的一个典型特征是离散化。而式（3.1.4）定义的傅里叶变换本质上是一个积分计算，体现为连续化特征，同时在实际应用中信号都是通过离散化采样得到的。为了通过离散化来采样信息以及有效地利用计算机实现傅里叶变换的计算，需要对式（3.1.4）实现高效、高精度的离散化。为此，需要导出离散傅里叶变换（DFT）的概念。为简单计，设f(x)为-,上的有限信号，则f(x)的傅里叶变换可简化为： 式（3.1.5）再假设采用等间距采样，其采样点数为N，输入时域信号为fk，要求输出频率信号为。为了利用采样点fk得到尽可能符合式（3.1.4）的输出值，DFT的思想是根据fk拟合出f(x)的最佳逼近多项式S(x)，然后在式（3.1.4）中利用S(x)代替f(x)，从而得到。下面简要讨论与的求法。给定一组正交基：，。直接验证向量满足内积关系：，其中IN为N阶单位矩阵。 式（3.1.6）设，利用正交基求解最小二乘问题：= 式（3.1.7）求解式（3.1.4）得到： 式（3.1.8）现在利用S(x)的定义，以及由式（3.1.8）得到的系数值ck来近似计算。将式（3.1.8）中的系数值代入多项式函数S(x)中，并利用S(x)作为f(x)的近似，则有： 式（3.1.9）除开常数2外，式（3.1.9）即为通常意义的离散傅里叶变换（DFT），其中输入fn与输出分别为信号的时域与频域信息。3.1.2 短时傅里叶变换尽管傅里叶变换及其离散形式DFT已经成为信号处理，尤其是时频分析中最常用的工具，但是，傅里叶变换存在信号的时域与频域信息不能同时局部化的问题。例如，从定义式（3.1.4）我们看到，对于任一给定频率，根据傅里叶变换不能看出该频率发生的时间与信号的周期（如果有的话），即傅里叶变换在频率上不能局部化。同时，在傅里叶变换将信号从时域上变换到频域上时，实质上是将信息在整个时间轴上的叠加，其中起到频限的作用，因此，傅里叶变换不能够观察信号在某一时间段内的频域信息。而另一方面，在信号处理，尤其是非平稳信号处理过程中，如音乐、地震信号等，人们经常需要对信号的局部频率以及该频率发生的时间段有所了解。由于标准傅里叶变换只在频域有局部分析的能力，而在时域内不存在局部分析的能力，故Dennis Gabor于1946年引入短时傅里叶变换（Short-Time Fourier Transform）。短时傅里叶变换的基本思想是：把信号划分成许多小的时间间隔，用傅里叶变换分析每个时间间隔，以便确定该时间间隔存在的频率。信号分析和处理的一个重要任务，一方面是要了解信号所包含的频谱信息，另一方面还希望知道不同频率所出现的时间。给定一信号，其STFT定义为 式（3.1.10）式中 式 (3.1.11）及，并且窗函数应取对称函数。STFT的含义可解释如下：在时域用窗函数去截（注：将，的时间变量换成），对截下来的局部信号作傅立叶变换，即得在时刻得该段信号得傅立叶变换。不断地移动，也即不断地移动窗函数的中心位置，即可得到不同时刻的傅立叶变换。这些傅立叶变换的集合，即是，如图3.1所示。显然，是变量的二维函数。图3.1 STFT算法示意图由于是窗函数，因此它在时域应是有限支撑的，又由于在频域是线谱，所以STFT的基函数在时域和频域都应是有限支撑的。这样，（3.1.10）式内积的结果即可实现对进行时频定位的功能。当然，我们自然要关心这一变换时域及频域的分辨率。对（3.1.11）式两边作傅立叶变换，有式（3.1.12）式中是和等效的频率变量。由于 式(3.1.13)所以 式（3.1.14)该式指出，对在时域加窗，引导出在频域对加窗。由图3.1可以看出，基函数的时间中心（注意，是移位变量），其时宽 式(3.1.15)即的时间中心由决定，但时宽和无关。同理，的频率中心，而带宽 式(3.1.16)也和中心频率无关。这样，STFT的基函数具有时频平面上的一个如下的分辨“细胞”：其中心在处，其大小为，不管取何值（即移到何处），该“细胞”的面积始终保持不变。该面积的大小即是STFT的时频分辨率。如图3. 2所示。图3.2 STFT的时频分辨率当我们对信号作时频分析时，一般，对快变的信号，我们希望它有好的时间分辨率以观察其快变部分（如尖脉冲等），即观察的时间宽度要小，受时宽带宽积的影响，这样，对该信号频域的分辨率必定要下降。由于快变信号对应的是高频信号，因此对这一类信号，我们希望有好的时间分辨率，但同时就要降低高频的分辨率。反之，对慢变信号，由于它对应的是低频信号，所以我们希望在低频处有好的频率分辨率，但不可避免的要降低时域的分辨率。3.1.3 窗函数以及窗函数宽度对时间-频率分辨率的影响假设对信号f (x)在时间x=附近内的频率感兴趣，显然一个最简洁的方法是仅取式（3.1.2）中定义的傅里叶变换在某个时间段I内的值，即定义 式(3.1.17)其中|I|表示区域I的长度。如果定义方波函数g(x)为 式(3.1.18)则式(3.14)又可以表示为: 式(3.1.19)其中R表示整个实轴。从式（3.1.4）、式（3.1.18）与式（3.1.19）很容易看到，为了分析信号f (x)在时刻的局部频域信息，式（3.1.17）实质上是对函数f(x)加上窗口函数g(x)。显然，窗口的长度|I|越小，则越能够反映出信号的局部频域信息。容易得到下面的简单性质： 式(3.1.20) ， 式(3.1.21)将函数g(x)与著名的“函数”及其性质 式(3.1.22)以及 式(3.1.23)比较不难发现，“函数”(x)实际上可以视为函数g(x)的极限函数。从另外一个角度来看，窗口函数可以看作对于原信号在区域上的加权，而利用方波函数g(x)作为窗口函数时存在的一个明显缺陷就是在区域I上平均使用权值，不符合权值应该重点位于时刻且距离该时刻越远和权值越小的特点。也就是权函数主值位于时刻，在该时刻的两端函数图像迅速衰减的特点。在满足上述特性并保持函数的光滑性质的前提下，Dennis Gabor于1946年提出了利用具有无穷次可微的高斯函数 式(3.1.24)作为窗口函数。图3.3（b）给出了取几种不同的值时高斯函数的图像，显然高斯函数具有窗口函数所需要的性质。下面讨论高斯函数与函数的关系。 （a）窗口函数g(x)的图形 （b）a取值不同时高斯函数的图形图3.3 窗口函数与高斯函数的图形例如设由两个高斯幅度调制的chirp 信号迭加而成，这两个信号一个时间中心在处，时宽，另一个时间中心在处，时宽也是32，调制信号的归一化频率都是0.25。在时频分布中，往往都称为一个“时频原子（atom）”，在该例的中，包含了两个时频原子信号。选择为Hanning窗，取窗的宽度为55，其STFT如图3.4a所示，这时频率定位是准确的，而在时间上分不出这两个“原子”信号的时间中心，我们将窗函数的宽度减为13，所得STFT如图3.4b所示，这时，在时间上也实现了两个中心的定位。以上例子说明了窗函数宽度的选择对时间频率分辨率的影响。总之，由于受不定原理的制约，我们对时间分辨率和频率分辨率只能取一个折中，一个提高了，另一个就必然要降低，反之亦然。（a）窗函数宽度为55（b）窗函数宽度为13图3.4窗函数宽度对时频分辨率的影响3.1.4 离散信号的短时傅立叶变换及其实现当我们要在计算机上实现一个信号的短时傅立叶变换时，该信号必须是离散的，且为有限长。设给定的信号为，对应（3.1.10）式，有 式（3.1.25）式中是在时间轴上窗函数移动的步长，是圆周频率，为由得到的抽样间隔。该式对应傅立叶变换中的DTFT，即时间是离散的，频率是连续的。为了在计算机上实现，应将频率离散化，令式（3.1.26）则 式（3.1.27）上式将频域的一个周期分成了M点，显然，上式是一个标准的M点DFT，若窗函数的宽度正好也是M点，那么上式可写成， 式（3.1.28）若的宽度小于M，那么可将其补零，使之变成M，若的宽度大于M，则应增大M使之等于窗函数的宽度。总之，（3.1.28）式为一标准DFT，时域、频域的长度都是M点。式中N的大小决定了窗函数沿时间轴移动的间距，N越小，上面各式中m的取值越多，得到的时频曲线越密。若，即窗函数在的时间方向上每隔一个点移动一次，这样按（3.1.28）式，共应做个M点DFT。当然，这时前和后个DFT所截的数据不完全，得到的效果不够好。3.2小波变换算法在短时傅立叶变换中我们使用了固定的时间窗函数，这就引出了时间分辨率和频率分辨率的概念，时间分辨率和频率分辨率是一对矛盾。根据海森堡的测不准原理，即时间窗函数的长度越长，频率分辨率就越高，而对于时间分辨率则越差。为了平衡时间分辨率和频率分辨率这个矛盾，可以采取对存在高频分量的部分采用高的时间分辨率和低的频率分辨率，而对于低频分量则采用高的频率分辨率和低的时间分辨率的方法，这就是多分辨分析的思想。 小波变换是一种在时间-尺度平面内，利用多分辨率分析思想分析非平稳号的方法.3.2.1 连续小波的定义 前面讨论的短时傅里叶变换（STFT）其窗口函数通过函数时间轴的平移与频率限制得到，由此得到的时频分析窗口具有固定的大小。对于非平稳信号而言，需要时频窗口具有可调的性质，即要求在高频部分具有较好的时间分辨率特性，而在低频部分具有较好的频率分辨率特性。为此特引入窗口函数，并定义变换5： 式（3.2.1）其中，aR且a0。式（3.2.1）定义了连续小波变换，a为尺度因子，表示与频率相关的伸缩，b为时间平移因子。当满足允许条件： 式（3.2.2）时，称为一个基本小波或母小波。母小波经过伸缩和平移后得到窗口函数。 式（3.2.3）式中：a为尺度因子；b为平移因子。若定义小波母函数的窗口宽度为，窗口的中心为，则相应可求得连续小波的窗口中心为，窗口宽度为。的频域窗口中心为，其频域窗口中心，对于任意函数的连续小波变