基于MATLAB的谱估计实现毕业设计

西 安 邮 电 学 院毕 业 设 计(论 文)题 目: 基于MATLAB的谱估计实现 系 别: 通信工程 专 业: 通信工程 班 级:通工0509班学生姓名: 钮金鑫导师姓名: 邢务强 职称: 讲师起止时间: 2009 年 3 月 10 日至2009年6 月11日 西 安 邮 电 学 院毕业设计论文任务书 学生姓名钮金鑫指导教师 邢务强职称 讲师系别通信工程系专业通信工程题目基于Matlab的谱估计实现 任务与要求 掌握随机信号功率谱的相关理论及其性质,在此基础上讨论功率谱估计的各种原理方法,比较分析功率谱估计的各自特点。利用MATLAB对谱估计的方法进行仿真实现,比较各自的性能指标及参数选取,最后通过实例加以说明。开始日期完成日期系主任签字年月日西 安 邮 电 学 院毕 业 设 计 论文 工 作 计 划 学生姓名_钮金鑫_指导教师_ 邢务强 _职称 讲师_ 系别_通信工程系_ _专业_通信工程_ _ 题目_基于MATLAB的谱估计实现_ _ 工作进程 第3周熟悉谱估计的基本理论知识第4周熟悉谱估计的基本理论知识,提交开题报告第5周熟悉MATLAB的相关基本知识。第6周阅读相关参考文献。第7周研究经典谱估计的方法。第8周进行经典谱估计系统仿真。第9周比较经典谱估计的各种方法,并进行中期检查。第10周研究现代经典谱估计的方法。第11周进行现代谱估计系统仿真。第12周比较现代谱估计的各种方法,通过实例加以说明。第13周总结以上讨论结果,写出论文初稿。第14周对毕业论文进行修改,并进行后期检查。第15周总结、完成毕业设计论文。第16周完成毕业设计答辩。 主要参考书目 离散随机信号处理 清华大学出版社 信号检测与估计 化学工业出版社 数字信号处理?时域离散信号处理 西安电子科技大学出版社 随机信号的功率谱估计及MATLAB实现 现代电子技术2002年第3期 经典谱估计方法的MATLAB分析 华中理工大学学报 主要仪器设备及材料 带有MATLAB软件的计算机一台 论文设计过程中教师的指导安排 每周指导一次,帮助理解较为繁杂的公式和每种方法的关键步骤, 并指出其中的重点难点,检查上次布置的任务。 对计划的说明 本计划为开题之初所定,后续会根据具体情况随时调整。最终一定按照规定结束日期完成毕业设计 西安邮电学院毕业设计论文开题报告 通信工程 系 通信工程 专业 05 级 09 班 课题名称: 基于MATLAB的谱估计实现 学生姓名: 钮金鑫 学号:05052112 指导教师:邢务强报告日期: 2009年3月16日 功率谱估计PSD 是用有限长的数据来估计信号的功率谱, 它对于认识一个随机信号或其他应用方面来讲都是重要的, 是数字信号处理的重要研究内容之一,在军事,生物医学,通信等领域得到了较为广泛的应用。 功率谱估计可以分为经典谱估计非参数估计 和现代谱估计参数估计。前者的主要方法有BT PSD估计法和周期图法; 后者的主要方法有最大熵谱分析法AR 模型法、Pisarenko 谐波分解法、Prony 提取极点法、Prony 谱线分解法以及Capon 最大似然法。其中周期图法和AR 模型法是用得较多且最具代表性的方法。 对于本课题来讲,应首先了解随机信号的基本知识,掌握信号与系统中傅里叶变换的相关知识,掌握MATLAB中一些基本语法,全面了解MATLAB中提供用于功率谱估计的函数。2.本课题需要重点研究的关键问题、解决的思路及实现预期目标的可行性分析 本课题需要研究的关键问题是: (1),研究掌握随机信号的基本理论; (2),研究每种功率谱估计的方法。 (3),学习MATLAB的相关知识。 解决的思路: (1)查看信号处理方面的书籍,复习关于随机信号的基本理论(2)上网查找关于一些功率谱估计的资料并参考其中的一些方法,重点研究有关理论并进行分析。(3)学习MATLAB的基本知识及用MATLAB进行功率谱估计的方法,了解每种方法的优缺点,并针对被估计信号的特点选择合适的方法。 预期目标的可行性分析:本课题要求对随机信号、信号与系统、MATLAB等学科的知识有较全面的了解,在大学期间,对随机信号分析,MATLAB这两门课程已有初步了解,对信号与系统了解较深,在导师的指导下,完成本课题是完全有可能的。 完成本课题要认真学习信号处理、随机过程的相关知识,掌握随机信号功率谱估计基本原理和计算机仿真的基本方法,并对其进行谱估计有一个清晰的认识,故应复习随机信号基本理论和MATLAB语言的使用。为此制定完成本课题的工作方案如下: 第 3 周 熟悉随机信号的基本理论知识 第 4 周复习信号与系统中的相关知识,提交开题报告 第56周学习谱估计的基本方法。 第78周阅读相关参考文献。 第910周 学习MATLAB的相关知识。 第1112周 对谱估计的每种方法进行仿真。 第 13 周总结仿真结果,写出论文初稿。 第 14 周 对毕业论文进行修改,并进行后期检查。 第 15 周 总结、完成毕业设计论文。 第 16 周 完成毕业设计答辩。指导教师签字: 2009 年 3 月 20 日说明: 本报告必须由承担毕业论文设计课题任务的学生在毕业论文设计 正式开始的第1周周五之前独立撰写完成,并交指导教师审阅。西安邮电学院毕业设计 论文成绩评定表学生姓名钮金鑫性别男学号05052112专 业班 级通工0509班课题名称 基于MATLAB的谱估计实现 课题类型理论分析难度一般毕业设计(论文)时间2009 年 3 月 10 日 6 月 11 日指导教师邢务强职称讲师 课题任务完成情况论文 20 千字; 设计、计算说明书 千字; 图纸 张;其它含附件:指导教师意见分项得分:开题调研论证 分; 课题质量(论文内容) 分; 创新分;论文撰写(规范)分; 学习态度分; 外文翻译 分 指导教师审阅成绩: 指导教师签字: 年 月 日评阅教师意见分项得分:选题分; 开题调研论证 分; 课题质量(论文内容) 分; 创新分;论文撰写(规范)分; 外文翻译 分评阅成绩: 评阅教师签字: 年 月 日验收小组意见分项得分:准备情况 分; 毕业设计(论文)质量 分; (操作)回答问题 分验收成绩: 验收教师组长签字: 年 月 日答辩小组意见分项得分:准备情况 分; 陈述情况 分; 回答问题 分; 仪表分 答辩成绩: 答辩小组组长签字: 年 月 日成绩计算方法填写本系实用比例指导教师成绩 20 % 评阅成绩 30 % 验收成绩 30 % 答辩成绩 20 %学生实得成绩百分制指导教师成绩 评阅成绩 验收成绩答辩成绩总评答辩委员会意见 毕业论文设计总评成绩等级: 系答辩委员会主任签字: 系签章 年月日备 注 西安邮电学院毕业论文设计成绩评定表续表目录摘要IABSTRACTII1 引言12 随机过程及估计理论基本知识22.1 随机过程的基本知识22.1.1 随机过程的分类22.1.2 离散随机过程的数字特征22.2 离散随机过程的功率谱密度32.2.1 离散随机过程功率谱密度定义32.2.2 随机过程功率谱密度的性质42.2.3 功率谱密度与自相关函数的关系42.3 估计量的评价标准42.3.1 无偏性42.3.2 有效性52.3.3 相合性53 经典谱估计及其仿真53.1 周期图53.1.1 周期图的定义53.1.2 周期图的性能73.1.3 平均周期图93.1.4 加窗周期图103.1.5 周期图的仿真123.1.6 平均周期图的仿真133.1.7 加窗周期图的仿真143.2 BT法163.2.1 BT法163.2.2 BT法的仿真194 AR模型法214.1 AR模型法的基本理论214.1.1 AR模型的Yule-Walker方法224.1.2 AR模型的协方差方法与修正协方差方法224.1.3 burg递推法254.2 AR模型法的仿真254.2.1 Yule-Walker方法254.2.2 协方差方法与修正协方差方法264.2.3 burg递推法274.2.4 AR模型阶的选择285 结论30摘要 功率谱估计是信号处理领域的重要问题之一。在实际应用中往往不能获得具体信号的表达式,需要根据有限的数据样本来获得较好的谱估计效果。 本论文研究了功率谱估计的几种常用的方法,包括经典法和参数法中的AR模型法,对每种方法的估计质量做了数学推导,给出仿真程序及仿真图,在仿真图的基础上对每种方法的性能进行了讨论。 经典法主要包括周期图法和BT法,但这两种方法都存在明显缺陷,即认为观测数据之外的数据都为零,所以对经典法中的周期图法进行了加窗,平均等修正。 从仿真图中可以看出,经典法无论是在谱分辨率还是在主瓣宽度,旁瓣幅度等方面都存在不足。AR模型法通过应用随机过程可用一白噪声通过一物理网络表示的理论,避免了经典法中除了观测数据之外的数据都为零的假设,因此能得到较好的谱估计效果。 关键字:功率谱估计;周期图;BT法;AR模型法。ABSTRACT Power Spectral Estimation is a fundamental topic in the signal processing field. Actually, we cant get the precise expression of a specific signal, so we need to estimate the power spectral of a signal according to some sample data sequences. In this thesis, some common methods of Power Spectral Estimation, such as classical spectral estimation、AR model of modern spectral estimation, are studied. The quality of each estimation method is derived, and the program and some simulated figures of these methods are also provided. In addition, the performance of each method is discussed according to the simulated figures. Classical methods of Power Spectral Estimation mainly include The Periodogram and the BT method. But both of them have a common drawback: the data sequences, beyond the area of the observed sequences, are all presumed to zero. So the Windows and the average method are introduced to improve the quality of The Periodogram From the simulated figures, we see that the classical methods have some conspicuous disadvantages, no matter in the spectral resolution or the main lobe width、the side lobe amplitude. AR model, which supposed that any random processes can be expressed as white noise crossing a physical union, avoids the drawback of the classical methods. So, better results of the estimation of the power spectrum can be obtained. Keywords: Power Spectrum Estimation;The Periodogram;the BT methods;AR model.1 引言 所谓功率谱估计,就是用已观测到的一定数量的样本数据估计一个平稳随机信号的功率谱,它在随机信号的分析中起着类似于频谱在确定性信号分析中所起的作用。功率谱是随机信号的一种重要的表征形式,在雷达信号处理中,由回波信号功率谱密度、谱峰的密度、高度和位置,可以确定运动目标的位置、辐射强度和运动速度等。在电子战中,谱分析可用来对目标进行分类识别等。功率谱估计在各种随机信号处理中得到了十分广泛的应用。例如,根据信号、干扰与噪声的功率谱,可以设计适当的滤波器,以尽量不失真地重现信号,而最大幅度地抑制干扰与噪声。 功率谱估计分为两大类。一类是非参数化方法,另一类是参数化方法。非参数化方法又叫经典谱估计法,它实质上仍依赖于传统的傅里叶变换法。经典谱估计法通常又分为两种,一种是间接法,它是由布莱克曼(/.ckman)和图基(/.ey)提出,又称为BT法。它首先对信号的自相关函数进行估计,然后再估计信号的功率谱。另一种是直接法,通过对观测到的数据样本进行傅里叶变换,然后将所得到结果的幅值平方后得到功率谱估计,这种方法又称为周期图法。 经典法的缺点表现为,除了得到的个数据以外,序列的其他值均被认为是零(或者等效,序列的自相关函数值除了能估计出的有限个值之外,其他的值被当做零),但序列或其自相关函数的那些未能观测到或未估计出来的值,实际上并不全是零。在现代谱估计的方法中,不再认为在观测到的个数据以外的数据全为零,因此克服了经典法的缺点,并且提高的谱估计的分辨率,能给出较好的估计。 本篇文章将首先讨论各种方法的理论依据,然后针对每种方法进行MATLAB仿真,给出详细的仿真程序,并在仿真图的基础上分析它的性能。2 随机过程及估计理论基本知识2.1 随机过程的基本知识2.1.1 随机过程的分类 按照时间和状态(一般称随机过程在的可能取值为它的状态)是连续的还是离散的来分类,可以分成以下四类:连续型随机过程:对于任意的,都是连续型随机变量,也就是时间和状态都是连续的情况。离散型随机过程:对于任意的,都是离散型随机变量,也就是时间连续,状态离散的情况。连续随机序列,随机过程在任一离散时刻的状态是连续型随机变量,也就是时间离散,状态连续的情况。例如在时间域上对接收机输出噪声电压过程进行采样,就可得到一个连续序列,其中。离散型随机序列:相应于时间和状态都是离散的情况。例如,对连续型随机序列再进行量化,即得到离散随机序列。 本篇论文主要涉及的信号类型属于离散随机信号,它一般是对一个连续信号进行抽样,对抽样所得的序列进行功率谱估计,故它属于连续随机序列。2.1.2 离散随机过程的数字特征 离散时间随机过程的均值或数学期望定义为 若是单值函数,则构成一个新的离散时间随机过程,其均值可定义为 均值有下列性质: (1) ,即和的均值等于均值的和。 (2) ,即乘以一个常数的均值等于的均值乘以此常数。 (3) 若,则称与是统计独立的。 离散时间随机过程的均方值定义为 离散时间随机过程的方差定义为 由于和的均值等于均值的和,所以容易证明上式可写成 为非负函数,其平方根称作离散时间随机过程的标准差或均方差,即 一般来说,均值,均方值,方差都是的函数,但对平稳离散随机过程来说,它们与无关,都是常数,即 离散时间随机过程的自相关函数定义为 2.2 离散随机过程的功率谱密度2.2.1 离散随机过程功率谱密度定义 设为宽平稳离散时间随机过程,或简称为宽平稳随机序列,具有零均值,其自相关函数为 或简写为 2-1当满足的条件时,我们定义的功率谱密度为的离散傅里叶变换,并记为式中,是随机序列相邻各值的时间间隔。2.2.2 随机过程功率谱密度的性质功率谱密度为非负的。 功率谱密度是的实函数。对于实随机过程来说,功率谱密度是的偶函数。 4,功率谱密度可积。2.2.3 功率谱密度与自相关函数的关系 可以证明,平稳随机过程的自相关函数与功率谱密度之间构成傅里叶变换对,即 以及 这一关系就是著名的维纳?辛钦定理。它给出了平稳随机过程的时域特性与频域特性之间的联系,它是分析随机信号的一个最重要,最基本的公式。2.3 估计量的评价标准2.3.1 无偏性 设是总体的一个样本。是包含在总体的分布中的待估参数,这里是的取值范围。 无偏性的描述如下 若估计量的数学期望存在,且对于任意有 则称是的无偏估计量。 可以证明,样本均值(也就是总体的一阶矩)是总体均值的无偏估计量,样本方差是总体方差的无偏估计量。2.3.2 有效性 设与都是的无偏估计量,若对于任意,有 且至少对于某一个上式中的不等号成立,则称较有效。2.3.3 相合性 无偏性和有效性都是在样本容量固定的前提下提出的。我们自然希望随着样本容量的增大,一个估计量的值稳定于待估参数的真值。这样,对估计量又有下述相合性的要求。 设为参数的估计量,若对于任意,当时依概率收敛于,则称是的相合估计量。 即,若对于任意都满足:对于任意,有 则称是的相合估计量。 可以证明,样本阶矩是总体的阶矩的相合估计量。 相合性是对一个估计量的基本要求,若估计量不具有相合性,那么不论将样本容量取得多么大,都不能将估计得足够准确,这样的估计量是不可取的。3 经典谱估计及其仿真 经典谱估计方法分为两种:BT法和周期图法。下面逐一阐述。3.1 周期图3.1.1 周期图的定义 对于一个离散随机过程,由公式2-1,以及样本的一阶矩就是对总体期望的无偏估计可得: 3-1其中是信号的采样点数,是对此离散时间随机过程自相关函数的估计量。 定义矩形窗函数如下 则式3-1可写成令,上式可简写为即可看成是与的卷积。而的傅里叶变换 是有限长序列的傅里叶变换,显然是周期性的。当为实序列时,也为实序列,根据傅里叶变换的性质,的傅里叶变换为。 根据功率谱密度与自相关函数的关系,功率谱估计为自相关函数估计的傅里叶变换,即 直接将模的平方除以求得的功率谱估计的方法称为周期图法,其结果用表示,即 3-2 周期图法是利用数据的傅里叶变换直接求得的,而不再计算自相关函数,所以又称直接法。由于序列的傅里叶变化可利用计算而提高运算效率,这是周期图法的主要优点。3.1.2 周期图的性能 为了了解周期图法的谱估计效果,我们来讨论它的估计均值和方差。 令代入上式得 (3-3)式中 由于为两个矩形窗函数的卷积,因此它是一个三角窗函数,其傅里叶变换为 式3-3为两个函数乘积的傅里叶变换,等于他们各自傅里叶变换的卷积,即 由上式可知,除非是一个冲激函数,否则将不等于,故周期图作为功率谱的估计是有偏差的。 当时 故此时 3-4 因此,周期图作为功率谱估计,当时是无偏的,即渐近无偏。 下面讨论谱估计的方差。 为分析简单起见,通常假设零均值,方差为的实高斯白噪声序列,即功率谱密度为常数。按方差的定义应有 由于假设零均值,方差为的实高斯白噪声序列,根据式3-4可知,当时,。 为了求,先求在两个频率和处的协方差,最后令。 根据式3-2可得 3-5其中为矩形窗函数。利用正态白噪声、多元正态随机变量的多阶矩公式,有 将上式代入式3-6中,得将代入上式,可得 由上式可见,当时,。而且,无论取何值,总与在同一数量级。显然,周期图不是功率谱估计的一致估计,所以这种估计不是最好的估计方法。3.1.3 平均周期图 由于周期图估计的方差总是不趋于零,所以需要对其进行改进,其中一种方法就是平均周期图法。平均周期图法就是将长度为的观测数据截断成段,每一段长度为,分别计算每一段的周期图,最后进行平均得到平均周期图的估计,其数学表达式如下: 由此可见,估计的方差很明显多了一项,这样方差便随着分段数目的增加而下降。 为了在平均周期图的基础上更进一步降低用周期图法进行功率谱估计的方差,可以采用的方法是使每一个分段的数据有一些重叠,使得观测的数据段的数目增多,以此来达到进一步减小方差的目的。3.1.4 加窗周期图 周期图法相当于只取了原信号的一段进行估计,也就是说对原信号加了一个矩形窗。时域上乘以一矩形窗相当于频域与一个函数进行卷积。在理想的情况下应该是原信号乘以幅值始终为1的信号,此信号在频域表现为一冲激函数,主瓣无限窄。而函数的主瓣不是无限窄,主瓣宽度反比于数据记录长度,而在实际中又不可能获得很长的数据长度。如果原来真是的功率谱是窄的,那么与主瓣卷积后会使功率向附近频域扩散,卷积以后使得加窗的信号的频谱变成函数频谱的形状,即有一定宽度,存在若干旁瓣的光滑脉冲那样。这就是说,若信号中包含两个以上正弦信号且频率非常接近时,两个分量的频谱叠加的结果,两个正弦信号不可能分辨,使信号模糊,降低了分辨率。 对其进行改进的办法是将矩形窗换为其它的窗函数。 几种常用窗函数的数学表达式如下: 矩形窗: 三角窗(bartlett): 汉宁窗: 哈明窗: 布莱克曼窗: 在MATLAB中还给出了几种常用的窗函数,他们的时域和频域特性分别如下: 图3-1 图3-2 从这几种窗函数的频域特性可以看出,矩形窗的主瓣宽度最小,因此运用矩形窗函数以后,谱分辨率是最高的,但旁瓣幅度较大,功率泄露情况也比较严重。旁瓣较低的窗是blackman窗和hamming窗,但这两种窗的主瓣宽度比矩形窗的主瓣宽度大,因此分辨率没有矩形窗的分辨率高。 采用加窗周期图方法的数学表达式为 其中,为窗函数。3.1.5 周期图的仿真 在MATLAB中用周期图方法进行功率谱估计的函数为periodogram。 仿真程序如下: Fs600; n0:1/Fs:1; xncos2*pi*40*n+5*cos2*pi*90*n+0.1*randnsizen; nfft1024; Pxx,fperiodogramxn,window,nfft,Fs; figure1 plotf,10*log10Pxx;图3-3 由前述讨论可知,在周期图方法中,无论怎样增大数据采样点数,功率谱估计的方差一直不为零。所以周期图不是功率谱的一致估计。如图4-1,真实的谱线应该是在频率为40Hz和90Hz处为两个冲激,而周期图法估计出的功率谱在这两个频率出的主瓣非常宽,失真较明显。3.1.6 平均周期图的仿真 仿真程序如下: Fs600; %采样频率 n0:1/Fs:1;%产生含有噪声的序列 xncos2*pi*40*n+ 5*cos2*pi*90*n+0.1*randnsizen; K10; Pxxbartxn,10; index0:round256/2- 1;kindex*Fs/256;%k相当于频率,256是PERIODOGRAM(在bart函数中调用的)函数中默认的nfft的值 plot_Pxx10*log10Pxxindex+1; figure1 plotk,plot_Pxx; function sx,fbartxn,K Lfloorlengthxn/K; sx0; f0; n11; for k1:Kp,fperiodogramxnn1:n1+L-1;pp/K;sxsx+p;n1n1+L; end; 图3-4图3-5从图3-4和图3-3的对比中很明显可以看出,图3-4频率为40Hz和90Hz处的两个主瓣的宽度较图3-3而言明显变窄,更接近于原信号。 程序如下: Fs600; n0:1/Fs:1; xncos2*pi*40*n+5*cos2*pi*90*n+0.1*randnsizen; nfft512; windowboxcar100; noverlap50; Pxx,fpwelchxn,window,noverlap,nfft; index0:roundnfft/2- 1; kindex*Fs/nfft; plot_Pxx10*log10Pxxindex+1; figure1 plotk,plot_Pxx 数据重叠方法的仿真图如图3-5所示。从此图与图3-4的对比中可以看出,采用数据重叠方法后,主瓣宽度变得更窄,更接近于原信号。3.1.7 加窗周期图的仿真 采用加窗周期图方法进行功率谱估计的仿真程序如下: Fs600; %采样频率 n0:1/Fs:1;%产生含有噪声的序列 xncos2*pi*40*n+cos2*pi*41*n+3*cos2*pi*90*n+0.1*randnsizen; nfft1024; window1boxcarlengthxn; window2bartlettlengthxn; window3blackmanlengthxn; window4hamminglengthxn; window5hannlengthxn; window6kaiserlengthxn,2.5; window7trianglengthxn; Pxx1,f1periodogramxn,window1,nfft,Fs;%直接法 Pxx2,f2periodogramxn,window2,nfft,Fs; Pxx3,f3periodogramxn,window3,nfft,Fs; Pxx4,f4periodogramxn,window4,nfft,Fs; Pxx5,f5periodogramxn,window5,nfft,Fs; Pxx6,f6periodogramxn,window6,nfft,Fs; Pxx7,f7periodogramxn,window7,nfft,Fs; figure1 subplot331;plotf1,10*log10Pxx1;titleboxcar窗; subplot332;plotf2,10*log10Pxx2;titlebartlett窗; subplot333;plotf3,10*log10Pxx3;titleblackman窗; subplot334;plotf4,10*log10Pxx4;titlehamming窗; subplot335;plotf5,10*log10Pxx5;titlehann窗; subplot336;plotf6,10*log10Pxx6;titlekaiser窗; subplot337;plotf7,10*log10Pxx7;titletriang窗; %plotf,10*log10Pxx; figure2 plotf1,10*log10Pxx1;titleboxcar窗; figure3 plotf2,10*log10Pxx2;titlebartlett窗; figure4 plotf3,10*log10Pxx3;titleblackman窗; figure5 plotf4,10*log10Pxx4;titlehamming窗; figure6 plotf5,10*log10Pxx5;titlehann窗; figure7 plotf6,10*log10Pxx6;titlekaiser窗; figure8 plotf7,10*log10Pxx7;titletriang窗; 仿真图如下:图3-6 从图3-6可以看出,采用blackman窗和hamming窗后,功率谱的主瓣宽度比较窄,频率泄露的情况相对不太明显(也就是除主瓣以外的谱线的幅度较主瓣幅度相比较小)。 在上述程序中包含两个频率相差较小的正弦信号,但是运用不同窗函数,估计的分辨率也不同。对图3-6的主瓣进行放大后如图3-7所示. 图3-7 从图3-7可以看出,只有运用矩形窗和凯撒窗才能分辨出这两个频率相差较小的信号。所以运用矩形窗和凯撒窗进行功率谱估计的频率分辨率较高。但是在实际应用中,幅度相差100倍甚至更高的正弦分量同时存在于一个信号中是可能的,这种情况下,由于矩形窗的旁瓣幅度较大,难以区分这样的两个正弦分量,故估计的效果较差,故还应考虑旁瓣幅度较小的窗函数。 所以,在实际应用中,具体选择哪种窗函数,应该从旁瓣幅度和分辨率两个角度综合考虑到底选用哪一个窗函数进行功率谱估计的效果较好。 在采用周期图法进行功率谱估计时,可以将平均周期图法和加窗周期图法这两种方法结合在一起,效果会更好。对于此方法,在进行计算机仿真时,只需将平均周期图法中,window参数改成相应的窗函数即可。3.2 BT法3.2.1 BT法 BT法首先计算随机信号的自相关函数,若,可由的傅里叶变换的得到它的功率谱,即: BT法的估计步骤如下: 1,先估计随机信号的自相关函数。自相关函数由下式估计: 3-6 式中,的长度为。 2,求的离散傅里叶变换,得到功率谱的估值: 由于式3-6仅利用了的个有限值得到自相关函数的估计,因此它与信号本身的自相关函数会有一定程度的差别。下面讨论接近的程度。 首先求的均值,有 即自相关函数的估计值的数学期望等于序列的自相关函数的真实值。因此是自相关函数的无偏估计。 其次求的方差,根据方差的定义 3-7由式(3-6)可知,上式等号右边第一项可表示为 3-8 当是零均值正态序列时,它的各高阶矩都可以用其一阶和二阶矩来表示。对于零均值正态序列,有 所以 3-9将式3-9代入式3-8得3-10将式3-10代入式3-7,并注意到,得令,显然的最小值为,最大值为,且(即)的情况将出现次,的情况将出现此。依此类推,对不同的的情况,出现的次数将为,于是上式可写成当时,上式以趋于零。即故满足一致估计的条件。 式3-6估计自相关函数的方法,虽然当较大且时能得到一致估计,但当接近时,的估计方差变大,因而不能得到有用的估计。因此可以按下式估计自相关函数 其均值这相当于真值用三角窗函数加权。上式表明:当时,由于,是的有偏估计;若取有限值,则当时,的偏差为零,即,是的渐近无偏估计,其估计方差 由于是的渐近无偏一致估计,且其估计方差小于的估计方差,所以一般用作为自相关函数的估计,而较少使用。所以可表示为:3.2.2 BT法的仿真 仿真程序如下: Fs500;%采样频率 n0:1/Fs:1; %产生含有噪声的序列 xncos2*pi*40*n+3*cos2*pi*90*n+randnsizen; nfft512; cxnxcorrxn; %计算序列的自相关函数 CXkfftcxn,nfft PxxabsCXk;%CXk估计出来以后是复数,此处是对它取模 index0:roundnfft/2- 1; kindex*Fs/nfft;%k相当于频率 plot_Pxx10*log10Pxxindex+1; figure1; plotk,plot_Pxx; 仿真图如图3-8所示 图3-8(采样点数为500) 图3-9采样点数为5000 从前述讨论可知,BT法可通过增加采样点数N来使估计的方差减小,从而使性能得到改进。如图3-9所示,当采样点数为5000时,主瓣宽度明显变窄,旁瓣的起伏程度也较小。 同周期图法一样,BT法也存在谱分辨率的问题。解决的方法也可以在时域采用不同的窗函数,从而使得分辨率有所提高。4 AR模型法4.1 AR模型法的基本理论 任何具有功率谱密度的随机信号都可以看成由一白噪声激励一物理网络所形成。如果能根据已观察到的数据估计出这一物理网络的模型参数,就不必认为个以外的数据全为零,这就有可能克服经典谱估计的缺点。 实际应用中所遇到的随机过程常常可以看成由白噪声经意线性系统形成的。这个线性系统的传递函数为 当输入白噪声的功率谱密度时,输出的功率谱密度为 如果能确定与,的值,通过上式就可得到所需信号的功率谱密度。 若,除外其余的都为0,此时系统函数就变为 4-1其差分方程为 4-2其中,为白噪声通过物理网络后的输出信号,也就是待估计的信号。这样的模型称为阶AR模型。其系统函数只有极点,没有零点。模型输出功率谱密度为4-3只要能求得及所有的值,就可以得到随机信号的功率谱。4.1.1 AR模型的Yule-Walker方法 阶AR模型的系统函数和差分方程分别如式4-1,4-2所示,模型输出的功率谱密度如式4-3所示。将AR模型的差分方程代入的自相关函数表达式,得 4-4 按式4-2所示,只与相关,而与无关,故将式4-2代入式4-4中的第二项可得 4-5将式4-5代入式4-4得4-6将代入式4-6,并将两式合并后写成矩阵形式,得上式就是AR模型的Yule-Walker方程。对于实序列,由于,因此只要已知或估计出个自相关函数的值,可由该方程解出个模型参数,根据这些参数即可得到随机信号的功率谱估计。4.1.2 AR模型的协方差方法与修正协方差方法 介绍两种方法之前首先介绍一下前、后向线性预测误差滤波器的基本知识。 前向预测就是指利用预测,其估计值和预测误差可用下式表示 4-7其中就是指预测系数与阶数有关。式4-7表明是由时刻以前的个数据得到的估计。用均方误差最小准则求前向预测误差滤波器的最佳系数从而求得估计系数。 如果利用数据预测,则成为后向预测,其估计值用表示 一般,前向,后向预测用同一数据进行,即利用进行预测,为此,将上式改为 这样,前向预测是由预测,后向预测是由预测。 设后向预测误差用表示(实际表示的是信号在时刻的预测误差),这样 4-8 同样,利用最小均方误差准则可求得。 下面介绍AR模型的协方差方法和改进协方差方法。 在Yule-Walker中,没有用到前、后向预测的概念,所以由观测数据求预测误差功率的公式为: 假设信号的数据区在范围,有个预测系数,个数据经过冲激相应为的滤波器,输出预测误差的长度为,应用下式计算:显然,的长度长于数据的长度,上式中数据的两端需补充零点,这相当于无穷长的信号经过加窗处理,得到长度为的数据。 而在协方差算法中运用了前向预测算法求模型参数。由于使用的个数据都是已知的,对以外的数据不需要做任何假设,相当于避免了加窗概念的引入,所以得到的预测误差只存在于之间,所以,预测误差功率的公式如下式: 运用最小均方误差准则可求得模型参数。 修正协方差法使用前向和后向预测误差平均值最小的方法,估计AR模型的参数,又称为前向?后向法。其预测误差的平均功率为式中,和分表表示前向和后向预测误差功率。 由数字滤波器的相关理论可知,前向线性预测误差滤波器的系数和后向线性预测误差滤波器的系数互为共轭。再根据最小均方误差准则可求得修正协方差方法中的模型参数。4.1.3 burg递推法 在Yule-Walker中需要先由信号的观测数据估计自相关函数,这是它的缺点。伯格递推法则由信号的观测数据直接计算AR模型的参数。 设信号观测数据区间为:,前向,后向预测误差功率分别用和表示,公式分别为: 4-9 4-10 我们已经知道,前向,后向预测误差公式分别为式4-7和式4-8。在这两个公式中,信号项的自变量最大的是,最小的是。未来保证计算范围不超出给定的数据范围,在4-9和4-10中,选择的求和范围为。 为了求出最小时的模型参数,可先求出参数,再根据解Yule-Walke方程中导出的递推公式,进一步解出其他参数,从而可根据观测数据直接计算AR模型的参数。4.2 AR模型法的仿真4.2.1 Yule-Walker方法 Yule-Walker方法的仿真程序如下: fs200; n0:1/fs:1; xn cos 2*pi*40*n+3*cos 2*pi*90*n+0.5*randnsizen; nfft512; pyulearxn,10,nfft,fs 如前所述,对于Yule-Walker方法,当数据很长,估计较准确,但是当估计数据较短时,估计误差较大,会出现谱峰频率偏移与虚假谱线的问题。 Yule-Walker方法的仿真图如下: 图4-1(采样点为50) 图4-2(采样点为200) 从图4-1可以看出,当采样点为50时,在频率为40Hz处的谱峰明显向左偏移,且在频率为60-70Hz之间出现了虚假谱峰。当采样点位200时,如图4-2所示,虽然频率为40Hz处的谱峰偏移不明显,但是依然存在虚假谱峰。4.2.2 协方差方法与修正协方差方法 协方差方法和修正协方差方法克服了Yule-Walker方法中存在的频率偏移的问题且在分辨率方面有良好的性能,但虚假谱峰的问题依然没有较好的解决。 仿真程序如下: fs200; n0:1/fs:1; xncos2*pi*40*n+cos2*pi*41*n+3*cos2*pi*90*n+0.5*randnsizen; nfft512; order50; subplot211; pcovxn,order,nfft,fs;%AR模型协方差算法 subplot212; pmcovxn,order,nfft,fs%AR模型改进协方差算法 仿真图如下: 图4-3 从上图可以看出,协方差方法和修正协方差方法在分辨率方面的性能甚至比加窗周期图的分辨率还要好。所以这两种方法适合于估计正弦波的频率,但产生的虚假谱峰也比较多。4.2.3 burg递推法 如前所述,burg递推法由于避免了由有限个数据估计自相关函数的计算,因此基于burg递推法的谱估计要比自相关函数法准确。 burg递推法的仿真程序如下:fs200; n0:1/fs:1; xncos2*pi*40*n+cos2*pi*41*n+3*cos2*pi*90*n+0.5*randnsizen; windowboxcarlengthxn; nfft512; order150;%order1就是AR模型数学公式中的P figure1 pburgxn,order1,nfft,fs,half%half表示单边谱 图4-4 从图4-4可以看出,burg递推法估计出的谱线在宽度和幅度上面的性能都要优于Yule-Walker方法。4.2.4 AR模型阶的选择 在AR模型中,模型阶次的选择是一个关键问题。一般模型的最好选择是先验未知的,实际中需预先选定模型的阶次。选择模型阶次时,将产生对谱的平滑作用,降低谱的分辨率。当选择的模型阶次较大时会出现虚假的频率成分。这一点可由下图说明。如图4-5,估计选择的是阶数为50,估计的信号接近真实信号,可以分辨出频率为40Hz和41Hz的两个谱峰。如图4-6,估计选择的阶数是30,此时,估计的结果对频率为40Hz和41Hz的两个谱峰产生了平滑,降低了估计的分辨率。如图4-7,估计选择的阶数是80,虽然频率为40Hz和41Hz的两个谱峰较为明显,但是其他频率谱峰的起伏性也明显增加,因此,增加了虚假的频率成分。 图4-5(阶数为50) 图