对心电图（ECG）的谱估计.doc

数字信号处课程小论文题目：功率谱估计方法与实现的研究对心电信号（ECG）谱估计的研究摘要：心血管疾病是威胁人类生命的最主要疾病之一, 而心电图(ECG)是诊断心血管疾病的主要依据。对其的特征分析一直是医学信号处理的热点，本文针对心电信号的谱估计做了一些分析讨论，首先是对来自MIT-BIH数据库的心电信号进行了预处理，然后分析了其AR模型的阶次问题，最后是在MATLAB中，用Burg算法实现了ECG信号的谱估计。实验结果显示，心电信号的谱能够反应隐藏在心电信号中的疾病问题。关键词：心电信号 谱估计 频谱心电图 Burg算法目 录一、课题研究背景与意义31 心电图（ECG）及其谱估计简介32 功率谱估计简介43 功率谱估计国内外的研究历史和现状53.1 基于二阶统计量的功率谱估计的方法53.1.1 经典功率谱估计方法的原理和算法63.1.2 现代功率谱估计方法的原理和算法73.2 基于高阶统计量(HOS)的谱估计方法93.2.1 非参数估计法103.2.2 参数模型估计法103.3 基于分数低阶统计量(FLOS)的谱估计方法114 总结125 参考文献12二、心电图谱估计问题的基本方法和技术141 心电图谱估计研究的现状与意义142 MIT-BIH 心电图数据库153 AR模型功率谱估计的有关方法153.1自相关法173.2 Burg算法183.3 改进的协方差方法193.4 总结概述214 本文主要的研究内容21三、MATLAB实验与讨论221 MIT/BIH 心电图数据的读取222 心电信号的简单预处理233 AR模型阶次的选取244 Burg算法的实现305 心电图谱估计的实现326 实验结果与分析34四、结束语36参考文献：36附件：38一、课题研究背景与意义1 心电图（ECG）及其谱估计简介心脏是人体循环系统中的重要器官。由于心脏不断地进行有节奏的收缩和舒张活动，血液才能在闭锁的循环系统中不停地流动。心脏在机械性收缩之前，首先产生电激动，心肌激动所产生的微小电流可经过身体组织传导到体表，使体表不同部位产生不同的电位。如果在体表放置两个电极，分别用导线联接到心电图机（即精密的电流计）的两端，它会按照心脏激动的时间顺序，将体表两点间的电位差记录下来，形成一条连续的曲线，这就是心电图（ECG）1。它是反映心脏兴奋的电活动过程，对心脏基本功能及其病理研究方面具有重要的参考价值。心电图可以分析与鉴别各种心律失常，也可以反映心肌受损的程度和发展过程以及心房、心室的功能结构情况3。图1.1典型心电信号的一个节拍2心电信号的功率谱估计图，是一种心脏检查的新技术4567。它是利用信息处理方法，将普通心电图信号的时域信息变换为多域信息，进行多参数、多指标和相关动态综合分析，显示心电信号频域特征的图形，是一无创性检测技术，对心脏病变（尤其冠心病)诊断可提供一个较好的方法8。同时，功率谱上的峰特征（峰与峰间关系）载有反映心功能的重要信息。通过对峰特征的分析可以判定心功能的情况9。在体育运动领域，这个被用来作为评价运动员心脏功能的重要指标。在对运动员进行训练时，可以依据此来寻求更加良好的训练计划。所以对其的研究，无论是医学理论上还是临床实践上，都具有非常重要的意义和价值。2 功率谱估计简介功率谱估计问题就是根据随机过程的一组有限观测值来确定该过程谱的内容20。对于平稳随机过程来说，功率谱理论上的数值不可能实现，只能用有限的观测数据来估计或者逼近真实值，从而比较真实地反映问题的本质。当然估计结果的好坏，与拟合的数学模型及所采用的处理方法等有很大关系20。功率谱估计在实际工程中有重要应用价值，如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域，发挥了重要作用。所以开展这方面的研究具有重要的理论和实践意义20。特别是随着现代社会科学技术的发展，对随机现象研究的深入，功率谱估计作为一种基本的方法得到了越来越多的挖掘应用。各项科学技术的应用中遇到的问题也迫切需要其来解决，如近年来灾害天气频发，如何根据已有的观测数据来研究预测天气问题，对我国的经济发展和社会稳定都很重要。在现代通信领域，随着人们对日常通信质量越来越高的要求，必须对复杂情况下的通信信号进行更彻底的研究，因此对通信电子领域的工作者来说，功率谱估计便是一个很好的研究载体。所以，如何改善功率谱估计的质量，得到更加稳定、质量更好的功率谱是一个很有意义的研究领域。3 功率谱估计国内外的研究历史和现状功率谱估计自诞生以来，就受到了人们的广泛重视。国内外对其研究的新理论及方法也是层出不穷。最早的谱估计方法是建立在基于二阶统计量，即自相关函数的功率谱估计的方法上。功率谱估计的方法经历了经典谱估计法和现代谱估计法两个研究历程， 在过去及现在相当长一段时间里，功率谱估计一直占据着谱估计理论里的核心位置。随着人们对随机信号和噪声特性地不断深入研究，又陆续提出了基于高阶统计量和基于分数低阶统计量的谱估计方法，以逐步完善谱估计理论体系。由谱估计的发展历程可看出,不管是在深度上还是在广度上谱估计的研究方法都在不断的前进。下面将就谱估计的各种研究方法展开进行简单介绍，从而分析下其在国内外的研究现状。3.1 基于二阶统计量的功率谱估计的方法功率谱(密度)是由随机信号的二阶统计量,即自相关函数作傅里叶变换(FT)得到的,其估计方法分两大类：一是以FT作为理论基础的线性方法,称之为经典功率谱估计法,它又包含周期图法和相关图法以及它们后来的改进方法;二是基于非线性运算的现代谱估计法, 它又包含参数模型法、MUDR谱估计法及特征分解法等。3.1.1 经典功率谱估计方法的原理和算法1.周期图法(又称直接法)19世纪末的Schuster直接利用FT对观察数据,首先提出了周期图的概念,定义信号x(n)的周期图为:（式中为有限长序列的FT),即直接将的模的平方除以N求得功率谱估计的方法,称为周期图法。该方法是利用数据的FT直接求得的,而没有计算自相关函数, 所以又称为直接法。Schuster鉴于周期图的剧烈起伏现象，提出了平均周期图的概念，并指出了在对有限长数据计算傅立叶系数时存在的“边瓣”问题，即后来的窗函数的影响。Schuster将周期图成功应用于了太阳黑子活动周期的计算中，与天文学中文献的观测十分接近111314。2.相关图法或BT法(又称间接法)1958年由Blackman和Tukey运用平稳随机过程的维纳辛钦定理提出了相关图法，其基本思想是：通过改善对自相关函数的估计方法，来对周期图进行平滑处理以改善周期图谱估计的方差性能。其先由x(n)估计出自相关函数，然后对其求傅立叶变换得到x(n)的功率谱，将其作为对原信号功率谱的估计，即：但是，周期图法和相关图法本质一样，都是将有限长的数据段作为无限长的抽样序列给予开窗截断的结果。它们原理简单，便于实现，并有可采用FFT等技术而使计算量大为减小等优点，适用于长数据记录。但它们的主要缺点是分辨率低、估计方差大,尽管可以通过平均、加窗平滑等方法来改进，可经典法始终无法解决频率分辨率和谱估计稳定性之间的矛盾，特别是在记录短数据情况下，这一矛盾显得更加突出。为了克服以上缺点的过程，即是现代谱估计方法的产生和发展1516。3.1.2 现代功率谱估计方法的原理和算法现代谱估计技术的研究和应用主要起始于60年代。1967年Burg在地震研究中受到线性预测滤波的启发，提出了最大嫡谱估计方法。1968年Parzne正式提出了自回归谱估计方法。1971年,VnaderBos证明了一维最大嫡谱估计与自回归谱估计等效。1972年出现的谱估计的Prnoy方法在数学上与自回归方法有着某些类似。 目前以自回归滑动平均模型为基础的谱估计己经比自回归模型具有更高的频率分辨率和更好的性能。1972年Piasernko提出的谐波分解方法提供了可靠的频率估计方法，实际上它是以自回归滑动平均为基础的谱估计的特例。在此期间，罗宾逊(于1978年)等许多人都在探讨适用性更广的一类有理谱的自回归滑动平均谱分析方法。其中除可采用各种ARMA模型的拟合算法外，卡佐夫C(adzow)、凯(Kya)等人同于1980年提出一种只估计ARMA模型的AR参数，而免于估计其撇参数的简化算法。格雷等(于1978年)也提出G谱估计法有利于这种ARMA谱分析。1981年Scridt提出了谱估计的多信号分类算法，它对于正弦分量频率估计很有效，塔夫(TuftS)和库玛勒森(Kualllersna)于1982年应用数据矩阵奇异值分解法求特征值，改进了前、后向线性预测的谱估计方法，又于1983年在相同原理下改进为最小模算法，可分别简记为MFBLP算法和min一norm算法。该方法有利于在较高的信噪比下分清线谱，且避免定阶的困难。然而其计算量较大。 除此之外，还有后来发展起来的基于高阶统计量和低阶统计量的谱估计，下面分别进行简单的介绍。1.参数模型估计法Yule于1927年提出了使用线性回归方程来模拟一个时间序列，从而发现隐含在该时间序列中的周期。Yule的这一工作实际上就成了现代谱估计中最重要的方法参数模型估计法的基础。Yule提出的自回归方程和线性预测有着密切的关系，Khintchine，Slutsky，Wold等人于1938年给出了线性预测的理论框架，并首次建立了自回归模型参数与自相关函数关系的Yule-Walker方程。Barlett于1948年首次提出了用自回归模型系数来计算功率谱。自回归模型和线性预测都用到了1911年提出的Toeplitz矩阵结构，Levinson根据该矩阵的特点于1947年提出了计算Yule-Walker方程的快速计算方法。现代谱估计的提出主要是为了解决经典谱估计的分辨率低和方差性能不好的问题。1967年，Burg提出的最大熵谱估计，即是朝着高分辨率谱估计所做的最有意义的努力20。参数模型估计法的基本思想：先根据随机信号的一系列观测数据建立一个模型，然后再根据已知随机信号的观测值来估计模型参数，使得模型产生的过程最佳近似于波研究的过程; 最后利用估计得到的模型参数计算功率谱。由该思路可知，这种估计方法的关键在于模型参数的估计。参数模型法又包括AR模型法、MA模型法、ARMA模型法等。其中由于AR 模型的参数计算是线性方程计算，相对于MA模型或ARMA模型参数的非线性运算而言要简单，另外任意的MA模型或ARMA模型都可以用无限阶或阶数足够大的AR模型来表示，所以AR模型法是用得最广泛的一种模型谱估计法。下面简单介绍下AR参数模型估计法。随机信号x(n)可看作由一个均值为零、方差为2的高斯白噪声u(n)通过具有最小相位的线性时不变系统H(Z)得到的输出信号,其模型如图1所示20:图1.2 AR参数模型法示意图设p阶AR模型的差分方程和系统函数分别为和，则信号的估计频谱为：，由此可以知道，只要求得了模型参数a1,a2,.,ap及就可得到信号的功率谱估计了。对其模型参数及的计算可采用Yule-Walker法、协方差法、Burg法、Levinson-Durbin法及Itakura法等等。目前公认的最好的算法还是Burg方法，本文下面的实验中就会采用这种方法来进行功率谱的估算。在实践中，在对AR模型的参数估计中，不仅不同的参数估计方法对功率谱估计的精度有所不同，而且谱估计的质量受模型阶次的影响，阶次过低，谱太平滑，反映不出谱峰; 阶次过大，容易产生虚假峰值。阶次的选择需要在实际中对所得结果作多次比较后, 予以确定1718。AR谱估计在实践中存在的一些主要问题有下面几个。在AR谱估计中可能会出现“谱线分裂”现象，这就要求适当增加数据的记录长度，以减小谱峰位置对相位的依赖性。AR谱估计的效果受到制约的另一个很重要的因素就是它对观测噪声的敏感度。为了减小噪声对谱估计的恶化影响，可以对数据先进行滤波处理。AR 谱估计不足的地方是当信噪比(SNR)较低时其性能并不好，此时就要采用超分辨方法, 即多种基于矩阵奇异值分解或特征值分解的改进谱估计方法。2.特征分解法谱估计特征分解法谱估计主要是针对白噪声的复正弦信号的功率谱估计问题, 其主要思想是：把自相关矩阵中信息空间分成两个子空间，即信号子空间和噪声子空间，这两个子空间中的某些矢量函数在正弦波频率上有尖锐的峰，据此估计正弦波的频率。但这些函数并不是功率谱，因为它们并不反映随机过程的功率。特征分解法谱估计又包含Pisarenko谱估计法、MUSIC 谱估计法、MN 谱估计法等, 这些方法可以得到比AR模型法更高的分辨率和更准确的谱估计，尤其是在信噪比较低时此优点更为明显。除此之外，在基于特征分解技术和最大似然技术上，人们又提出了很多新的谱估计方法，使随机信号的谱及其淹没在噪声中的正弦波频率估计得到更好的估计效果1221。3.2 基于高阶统计量(HOS)的谱估计方法在前面叙述的基于二阶统计量的功率谱估计的理论和方法虽已相当成熟, 而且应用也极为广泛，可是在功率谱估计过程中，人们通常假设是观测数据由高斯白噪声序列激励线性最小相位系统产生的。另外自相关函数及相应的功率谱中只包含信号的幅度信息，不包含相位信息。但是在许多实际应用中，如地质勘探、水声信号处理与识别、远距离通信等场合下，信号或噪声可能不服从高斯分布; 系统可能不是线性最小相位系统；需要的可能不仅仅为信号的幅度信息，还需相位信息等，此时，功率谱估计法便显得无能为力。正因为如此, 基于高阶统计量()和分数低阶统计量()的谱估计方法越来越受到重视和应用。基于HOS的谱包含高阶矩谱和高阶累积量谱, 它们都是对相应的高阶统计量作多重FT得到的，高阶累积量谱又简称为高阶谱或多谱。由于累积量中含有相位信息，所以人们热衷于研究多谱估计方法, 现有的多谱估计法中又重点围绕三阶谱(双谱)和四阶谱(三谱) 研究，其中双谱是阶数最低的高阶谱()，处理方法简单，同时它包含了相位信息，其应用相对广泛，理论成熟，所以这里简单介绍下，其他不在叙述。双谱估计法与功率谱估计法大致相同, 都可分为非参数法和参数模型法。其中非参数法又包含直接法和间接法；参数模型法包含AR模型法、MA模型法、ARMA 模型法等。3.2.1 非参数估计法1.直接法直接使用FFT将有限长数据段x(n)变换至频域，得到x()，然后取代H()便得到多谱估计了。2.间接法先计算累量的估计，再求其傅里叶变换便得。非参数估计法的性能特点是，估计的方差较大，仅对较长的观测数据段比较适用，若用平均和平滑方法来减小估计方差，则其副作用是分辨率下降。3.2.2 参数模型估计法图1.3 基于高阶谱的随机信号模型如图所示，假设u(n)是平稳非高斯白噪声序列，且均值为零；系统H(Z)不是最小相位的，其为有理传输函数；测量噪声v(n)为高斯有色噪声，且均值为零。该参数模型估计法与基于自相关函数的参数模型估计法类似，都重在模型参数的辨识估计上, 在估计模型参数时要依据一定的准则，通常采用最小均方误差准则。现将几种参数模型估计法的方法简单介绍如下:1）估计AR 模型参数可采用Yule-Walker法。当AR 模型不稳定时, 即求非因果系统的AR 模型参数可采用穷搜索法、最优法、转换法等。2）估计MA 模型参数可采用线性代数法、闭合解法、非线性优化法, 其中闭合解法抗噪能力差, 难以实际使用：非线性优化法涉及到非线性问题，也难以实现，所以最实用的是线性代数法。3）估计ARMA 模型参数可采用剩余时间序列法、非线性优化法及相位恢复法等。3.3 基于分数低阶统计量(FLOS)的谱估计方法自20世纪80年代中期以来，基于HOS的谱估计法能够有效地解决以前的功率谱估计方法所无法解决的问题，如当信号为非高斯信号时，能充分获取隐含在数据中的信息；当系统为非最小相位时，能反映出原信号非最小相位的特点；当测量噪声为较大的有色噪声时，具有很好的抗有色噪声的能力等。但是基于HOS 谱估计应用的前提是所分析处理的信号和噪声要具有有限的二阶矩和高阶矩。实际上, 有许多信号和噪声不仅不存在高阶矩，甚至连二阶矩也不存在，所以就需要研究适用于稳定分布随机过程的新的谱，那么在其基础上基于FLOS的谱估计理论和方法就应运而生了。分数低阶统计量(FLOS)，是指低于二阶的，即(0,2)范围内的分数阶统计量。在实际应用中，我们所遇到的大量非高斯信号或噪声具有显著的尖峰脉冲特性，其概率密度函数具有较厚的拖尾。为了分析这类信号，必须选择与之匹配的稳定分布的参数模型，在此基础上定义的谱统称为分数低阶谱，包括谱、共变谱、分数低阶协方差谱、几何功率谱、分数低阶循环谱、分数低阶多谱等。对上述这些分数低阶谱的估计，发展出了许多的方法，但是从根本上而言，这些方法都是在功率谱的经典谱估计法的基础上改造的。基于各种模型参数的分数低阶协方差谱估计法还有待进一步研究。除了上面介绍的这些分数低阶谱估计外, 还需要寻找更多其它的基于FLOS的新谱及其估计方法，全方位广义化非高斯分数阶白噪声激励线性或非线性过程的参数化的谱估计理论和方法，构造FLOS- AR、FLOS- MA 及FLOS- ARMA 模型的基本参数化方法。这也是当前研究的重点问题222324。4 总结随着人们对随机信号的特性研究和谱的概念的建立, 新的谱估计方法不断产生和更新，以逐步形成一个完整的理论体系。谱估计方法来源于信号特性，最后又服务于信号的研究, 尤其是在分析信号的频域特性、时频特性方面起了不可估量的作用。由上面的叙述可知，信号的谱估计方法大致可分为两大类：一是非参数法；二是模型参数法。而非参数法谱估计性能不及模型参数法，因为它总存在估计方差大，分辨率低等缺点。所以，最近几十年里，人们总是热衷于模型参数的估计方法研究。但是不管是哪种谱估计方法,，基本上都是源于功率谱估计的理论和思想。毕竟功率谱估计理论已相当成熟，现在我们需要做的就是建立新的谱概念，找出新的谱估计方法，尤其是对基于FLOS或稳定分布的谱估计探索。随着现今人们对非高斯稳定分布研究的深入, 有望在不久的将来在这一领域有更大的进展和突破。5 参考文献1段新显，ECG监视过程中的消噪问题J，国外医学生物医学工程分册，1992, 15(3)。2来自网络，3王顺宝. 频域心电图的临床应用价值A. 第四届全国心功能学术研讨会论文摘要集C, 1994 .4沈庆乐，李培英，张存琪，原玉兰.频谱心电图对冠心病的诊断价值探讨J. 福建医药杂志, 1996,(01) .5王颖,杨建昌. 用频谱心电图评价心脏功能的研究进展J. 内江科技, 2009,(02) . 6黄佐贵. 频域心电图的临床意义及其诊断标准的研究J. 中国基层医药, 2002,(01)7陈范元,顾琴芳. 慢性阻塞性肺疾病患者频谱心电图测定的临床意义J. 上海医学, 2000,(05) .8付小娜,何勇. 频谱心电图诊断冠心病的临床价值的研究J. 河南医药信息, 2001,(02) .9徐志明,申行运,刘国印,唐志忠. 体位改变对飞行员频域心电图的影响J. 航天医学与医学工程, 2003,(03) .10张文泉,李彦斌. 最大熵谱估计及应用研究. 现代电力,2001,18(1):41- 45.11陈海英. AR模型功率谱估计常用算法的性能比较J. 漳州师范学院学报(自然科学版), 2009,(01) . 12储彬彬,王琛,漆德宁. AR模型功率谱估计的典型算法比较及MATLAB实现J. 中国新通信, 2008,(17) . 13黄志宇,刘保华,陈高平,张先军,冯松立. 随机信号的功率谱估计及Matlab的实现J. 现代电子技术, 2002,(03) . 14王福杰,潘宏侠. MATLAB中几种功率谱估计函数的比较分析与选择J. 电子产品可靠性与环境试验, 2009,(06) . 15姚文俊. 自相关法和Burg法在AR模型功率谱估计中的仿真研究J. 计算机与数字工程, 2007,(10) . 16伊鑫,曲爱华. 基于Welch算法的经典功率谱估计的Matlab分析J. 现代电子技术, 2010,(03) .17宋宁,关华. 经典功率谱估计及其仿真J. 现代电子技术, 2008,(11) .18王凤瑛,张丽丽. 功率谱估计及其MATLAB仿真J. 微计算机信息, 2006,(31) . 19傅广操,樊明捷. Matlab在现代功率谱估计中的应用J. 电脑学习, 2003,(06) .20胡光书.数字信号处理理论算法与实现.清华大学出版社,2003.21冯西安, 黄建国, 谢一清. 白噪声中高分辨频率估计特征分解法的统计性能研究J. 信号处理, 1995,(02)22罗静, 查代奉, 武良丹. 基于分数低阶统计量的双谱及其估计方法J. 信号处理, 2010,(03)23罗静,游芳. 基于分数低阶的双谱及其估计方法J.九江学院学报, 2009,(06) .24孙永梅. 稳定分布参数估计与谱分析理论及应用研究D. 大连理工大学, 2006二、心电图谱估计问题的基本方法和技术1 心电图谱估计研究的现状与意义心电图的研究和应用已有五十多年的历史。早在1960年，Schey和Ycung曾用扫描仪进行ECG频率分析。1961年Cady应用傅立叶级数表示ECG图形，并求出一个判断函数式，用来判定正常与左室肥厚的心电图形。1965年Cooley和Tukey创建了快速傅立叶转换的高效算法。1978年我国学者封根泉等从人体工程学观点出发，应用生物控制论概念，通过微机对心电信号的分析，提出了频谱心电图分析的新方法。臧益民、赖世忠等学者也分别对频谱心电图与心脏功能的关系进行了研究。目前，ECG的分析技术已能够对1000Hz以内的心电信息作全面的分析。常规心电信号周期性反映心电活动从心房到心室的除极和复极过程, 通过对P、Q、R、S、T、U 等波的分析, 可以诊断心律失常、心肌缺血、心肌梗死等多种疾病。因此常规心电图( ECG) 在临床上有着不可替代的作用。从信号处理角度看, 心电图波形又是一种正弦波组合而成的时间序列,而心电信号的谱图是常规ECG的频域转化图形，采集分析时间域、空间域和频率域的整个信息为一体, 尤其突出频率域的特点。能检出常规ECG 所不能反映的动态变化, 突破了传统ECG系统时间分析概念，是我国首创的一种无创性电检测技术。FCG从1985年问世以来，即作为心脏病诊断新的手段引起国内外医学界的重视。二十多年来经过一定数量的临床验证和研究，多数人认为其比常规ECG的信息量大，敏感性高，有助于冠心病的早期诊断。2 MIT/BIH 心电图数据库MIT/BIH心电数据库是由美国麻省理工学院和Beth Israel医院合作建立的数据库。MIT/BIH 数据库共有48个病例，每个病例数据长30分钟，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。并附有医学专家对每个心拍作出的识别和标注。在美国，其常被用作心电仪器的检验标准。MIT/BIH数据库中每个病例的ECG数据，是按照每通道360Hz的采样率、12位采样精度进行量化的，每个采样值2个字节。数据库中每个病例的ECG数据，都有两个导联的值，即MLII导联和V5导联。网站上可以直接下载得到。3 AR模型功率谱估计的有关方法设p阶AR模型的差分方程和系统函数分别为和。从上文可以知道，当知道系统的模型参数及和噪声功率就可得到信号的功率谱估计了。对其系数的求取是通过求解方程的正则方程，也就是Yule-Walker方程得到的。其如下所示： (2.1)对于某随机过程来说，如果已知，就可由上式求解得到了。但实际上，对于一个所要估计的随机过程来说，通常己知的仅仅是有限的观察数据，如x(0)，x(1)，x(N-1)，而其自相关函数值通常是未知的，因此AR模型参数的求得，通常是首先要求得自相关函数的估计值，然后再通过某种算法，得到AR模型参数的估计值。用线性方程组的常用算法(例如高斯消元法)求解(2.1)式，需要的运算量的数量级为。但若利用系数矩阵的对称性和Toeplitz矩阵的性质，则可得到一些高效算法，Lveinsno-Durbni算法是其中最著名、应用最广泛的一种。这种算法的运算量数量级为。这是一种按阶次进行递推的算法，即首先以AR(0)和AR(l)模型参数作为初始条件，计算AR(2)模型参数;然后根据这些参数计算AR(3)模型的参数，等等，一直到计算出AR(P)模型参数为止。这样，当整个迭代计算结束后，不仅求得了所需要的阶AR模型的参数，而且还得到了所有各低阶模型的参数。定义为阶AR模型在阶次为m时的第k个系数，k=1,2,m。m=1,2,p。为m阶时的前向预测的最小误差功率(且有)。由(2.1)式，当m=1时，有：求解得到：，。定义初始条件为：那么： 再定义m阶时的第m个系数，即为。称为反射系数，那么，由Toeplitz矩阵的性质，可得到如下Levinson-Durbin递推算法：Levinson-Durbin算法从低阶开始递推，直到阶次p，给出了在每一个阶次时的所有参数。这一特点特别有利于选择AR模型的合适的阶次。但由于进行自相关序列估计时，仍对0-（N-1)以外的观测数据为零的假设，必然引入较大误差，尤其对短序列功率谱估计将导致谱线分裂与谱峰频率偏移现象。以上接介绍的L-D算法，后来更发展除了更多的算法。下面介绍三种常用的求取模型参数的方法。3.1自相关法自相关法是AR模型参数求解中最简单的一种方法。其递推算法是在满足前向预测均方误差最小的前提下，先求得观测数据的自相关函数,然后利用方程（3.1）的递推性质求得模型参数，进而求得功率谱的估值。它是模型阶次逐次加大的一种算法，即先计算阶次m=1时的预测系数和和，再计算m=2时的，和，按此依次计算下去，直到m=p时，满足精度要求时即可停止递推。递推公式为： 其中，。3.2 Burg算法为了克服L-D算法中因估计相关函数给功率谱带来的影响，Burg提出了一种新的算法，其基本思想是直接从观测的数据中利用线性预测器得前向和后向预测的总均方误差之和为最小的准则来估计预测系数，进而通过L-D算法的递推公式求出AR模型优化的参数。前、后向预测误差的定义式分别为： （2.2）Burg算法与自相关法不同，它是使序列x(n)的前后向预测误差功率之和： 最小。在上式中，当阶次由1至p时，和有如下的递推关系： （2.3）其中，初始条件为，。式子中，表示m阶预测器结构模型用格型实现时的反射系数等于直接实现形式时m次的预测系数。在阶次为m时，令总的功率相对于为最小，则可以估计出反射系数。由此可以得到在前后向误差功率之和为最小时的。令，则求得反射系数等于：，其中，。 （2.4）按照上式估计出后，在阶次为m时的AR模型系数仍然由Levinson算法递推求出，即有： （2.5）其中，上面式子是假设在第（m-1）阶时的AR参数已经求出。Burg算法的优点在于频率分辨率高;所得的AR模型稳定;计算效率高。然而，Burg方法也有公认的不足。首先，在高信噪比时，它的谱线会呈现出分裂;其次对于高阶模型来讲，该方法也可能引入假谱峰;再次，对于噪声中的正弦信号，Burg方法对正弦信号的初始相位呈现出敏感性。特别是对短的数据记录更是如此。不管如何，一般情况下而言，Burg算法的分辨率是好于自相关法的。Burg算法的递推步骤如下：由初始条件，在由式（3.4）求出。由得到m=1时的参数，；由和式（3.3）求出和，再有式（3.4）求出。依照式（3.5）的Levinson递推关系，求出m=2时的，及。重复上述过程，直到m=p时，求出了所有阶次时的AR参数。3.3 改进的协方差方法确定AR模型参数的Burg方法基本上是具有预测系数满足Levinson递推的附加约束的最小二乘格型算法。由于这一约束的结果，增加AR模型的阶，仅需要在各级最优化一个参数。与这种方法相比我们可以利用无约束改进协方差法(最小二乘算法)确定AR参数。该算法的主要思路是为了摆脱因采用递推运算对确定预测系数的约束，让每一预测系数(模型参数)的确定直接让每一预测系数(模型参数)的确定直接与前、后向预测的总的平方误差最小(最小二乘法)联系起来。总的平方误差为： （2.6）它与Burg方法有相同的性能指标。但是，我们并不为了求AR参数而把Lveinsno递推加入上式。令，可以得到一组线性方程： （2.7）其中，自相关的定义为： （2.8）所得到的剩余最小二乘误差为： （2.9）所以，可得改进协方差功率谱估计为: (2.10)元素为的相关矩阵不是Toeplitz阵，因此，不能利用Levinson算法求解。但是，相关矩阵的结构足以设计复杂性与成正比的有效计算方法。Marple(1980)设计了这样的算法，他的算法具有格型结构，利用了Levinson类型的阶递推和另外的时间递推。这种形式的修正协方差(无约束最小二乘法)方法也一直被称为数据未加窗最小二乘方法。修正协方差方法基本上克服了谱线分裂、谱峰偏移和出现伪峰等缺点，在此意义上优于Burg方法。但此法不能保证AR模型稳定，并且所需的运算量也偏大，但是在谱估计中，并不认为这是一个大问题。3.4 总结概述上述的三种算法是求取模型参数的常用方法，都具有各自的特点。自相关法算法虽然简单，但是分辨率较差；Burg 算法容易出现谱线分裂和伪峰，受初相的影响产生频率偏移，但是它能满足大多数的应用要求，且计算不复杂；Marple 算法运算量较大，但性能最好，分辨率最高。因此，实际应用时，应根据实际的需求选择合适的算法。本文主要研究的是心电信号的谱估计，选择Burg算法已经能够满足其要求，故下文中使用Burg算法。所以对Burg算法做了详细介绍，并给出了算法的执行步骤。4 本文主要的研究内容本文的研究目标在于利用MIT/BIH数据库中的ECG数据，研究其频谱心电图的估计，在正确读取ECG信号的基础上，对其进行谱估计，以供医学临床上参考。本文将首先将讨论下AR模型阶次的选取，然后在MATLAB中实现burg算法，最后求取几个心电信号的FCG，并进行比较评价，与以前研究者的结果进行对比，以证明FCG的重要性。本文拟解决的关键问题如下：1. 读取MIT/BIH数据库中的ECG数据，将正常情况和各种病例区别开来。2. 对ECG信号进行简单预处理，去除其中的工频干扰。3. 针对ECG信号讨论其AR模型的阶次。4. 在MATLAB中使用基本语言实现Burg算法。5. 求取FCG。6. 对获得的FCG进行分析讨论。三、MATLAB实验与讨论1 MIT/BIH 心电图数据的读取18MIT/BIH 是由美国麻省理工学院提供的研究心律失常的数据库。从下面这个官方链接页面可以下载到所有48组MIT-BIH心电数据：/physiobank/database/mitdb/MIT为了节省文件长度和存储空间，使用了自定义的格式，所以没有通用的读取方式。其所记录的一个心电记录由三个部分组成： （1）头文件.hea，存储方式ASCII码字符。（2）数据文件.dat，按二进制存储，每三个字节存储两个数，一个数12bit。 （3）注释文件.art，按二进制存储，格式定义比较复杂，但是只要知道第（2）个如何读入matlab就可以类推出如何读第（3）个。 从网络上下载了专门的读取程序后，读取了数据库中100数据的信息，显示前1.6秒数据如下，其中的数字标注具有自己的含义，可以从网上查询得知。1表示正常，图中的28表示心率改变，心电记录起始的标志，其他数字分别表示各种心脏的问题。图3.1 MIT/BIH 心电图数据库中100段前1.6秒的心电图使用网页上自带的工具，将数据保存为txt文本格式后存盘，用作下面的实验。2 心电信号的简单预处理心电信号有随机性和背景噪声非常强的特点，是一种非线性、非平稳的微弱信号。常规心电信号是mV 级信号，频带范围是0.05100 Hz。这种微弱信号在一般的测量环境中，其检测不可避免的受到电网工频电磁辐射的干扰，工频干扰是在心电信号自身频率内的强干扰，会淹没原始心电波形中的细小特征，使心电波形和工频干扰相互混杂，因此，工频干扰抑制能力是决定心电信号检测质量的重要参数之一12。本文使用的MIT/BIH心电数据库中的心电信号，必然也含有心电信号以及其他的噪声，所以在使用前有必要进行预处理。由于本文的重点是心电信号的谱估计问题，所以对于滤波处理不做过多的要求，只进行简单的处理（心电信号的预处理本身是一个较复杂的问题）。MIT/BIH心电数据库采样于美国，所以其工频干扰是60Hz，本文使用简单的FIR陷波器对其进行预处理345。图所示是陷波器的幅频响应，滤波器的系数使用MATLAB数字信号处理工具箱得到，将其作为M代码后，见附件。将系数保存为mat文件后读入，与原始信号做卷积即得到预处理过的心电信号。图3.3，3.4所示是处理前后的心电信号对比。图3.2 心电信号预处理的陷波器，中心频率为60Hz图3.3 处理前的心电信号图3.4 处理后的心电信号3 AR模型阶次的选取AR模型的阶次事先是不知道的，而且其对AR模型的影响比较重大，阶次太小时，谱的分辨率不够，阶次过大，则会带来虚假谱峰的现象。在阶次的选择上，会有一些参考的准则，但是具体的值还是需要实验来对比研究得到。本文采用Burg算法，对下载得到的ECG信号进行实验验证，比较不同阶次时，FCG信号质量的好坏，评价标准是其的频率分辨特性和方差特性。频率分辨率越高，则表示信（a） 源信号（b）阶次p=1时估计得到的功率谱（c）阶次p=10时估计得到的功率谱（d）阶次p=50时估计得到的功率谱（d）阶次p=200时估计得到的功率谱图3.5 信号使用不同阶次的AR模型得到的谱估计对比号的谱越好，方差越大，则谱的杂乱程度越大，表示可能出现虚假谱峰810。为了验证上述结论，图3.2为正弦信号中加入白噪声所做的实验。原信号为：x=cos(2*pi*f1*n*Ts)+cos(2*pi*f2*n*Ts)+r(n); r(n)表示白噪声。从图中可以看出，当Burg算法阶次选择较小时，谱线是直的，但是随着阶次的增加，就会出现上述虚假谱线的现象，如图3.5（b）（c）（d）（e）所示。由此，我们可以得到如下结论：当阶数比较小时，Burg功率谱比较平坦，信号被噪声掩盖；随着阶数的增加，分辨率增加；当阶数进一步增加时，虚假峰出现，并逐渐增多，功率谱质量退化。给我们的启示是：随着阶数的增加，功率谱从平坦变尖锐，再从尖锐被越来越多的虚假峰淹没，其中一定存在某个阶数使得功率谱的质量比较好1217。因此，寻找一个较好的阶次来构造AR模型是十分有必要的。由文献12可知道，对心电谱，较为合适的阶次应该在20左右。本文分别选择10、20、26、30、60、100阶AR模型对进行上述信号进行频谱分析，图为其使用不同阶次AR模型谱估计的效果，从图中可以看出，使用20阶AR模型已经可以获得比较光滑的功率谱估计，但是频率分辨率不高，谱峰不容易分辨，随着模型阶次的增大，频率分辨率逐渐增加，到60阶AR模型时，虽然频率分辨率较大了，但是随之出现的是大量的虚假谱峰，因此综合权衡谱图的虚假谱峰以及频率分辨率的要求，本文选择26阶AR模型对ECG信号进行谱估计。（a）AR模型阶次为10得到的估计谱（b）AR模型阶次为20得到的估计谱（c）AR模型阶次为26得到的估计谱（d）AR模型阶次为30得到的估计谱（e）AR模型阶次为60得到的估计谱（f）AR模型阶次为100得到的估计谱图3.6 ECG信号使用不同阶次的AR模型得到的谱的对比4 Burg算法的实现确定了AR模型的阶次后，下面在MATLAB中用基本的代码实现Burg算法，而不是像上文中那样使用工具箱中的函数。下面依照上文给出的算法步骤，给出下MATLAB代码并进行分析解释。K=26; %AR模型的阶次N=max(size(x); %输入数据的长度r=;%存储格型结构的反射系数ef(1,:)=x; %存储格型结构的前向系数eb(1,:)=x; %存储格型结构的后向系数sigma(1)=sum(x.2)/N; %计算r（0）a=zeros(K+1,K+1); %存储每个阶次的系数a(:,1)=1;for p=2:K+1; r(p-1)=2*(ef(p-1,p:N)*eb(p-1,p-1:N-1)/(ef(p-1,p:N)*ef(p-1,p:N)+eb(p-1,p-1:N-1)*eb(p-1,p-1:N-1); %Burg算法中反射系数的计算，也就是公式（3.4）sigma(p)=(1-r(p-1).2)*sigma(p-1); %计算每个阶次的功率误差 for i=2:p-1 a(p,i)=a(p-1,i)-r(p-1)*a(p-1,p-i+1);%按照公式计算每个阶次的模型参数 end %也就是公式（3.5） a(p,p)=-r(p-1); for n=2:N ef(p,n)=ef(p-1,n)-r(p-1)*eb(p-1,n-1); %按照递推求解下个参数，也就是 eb(p,n)=-r(p-1)*ef(p-1,n)+eb(p-1,n-1);%公式（3.3） end endan=a(2:K+1,2:K+1); %取得所有的阶次的模型参数 这样就在MATLAB中，用较为简单的函数实现了的Burg算法，上述代码符合书上的介绍，下面我们使用MATLAB中自带的函数来对比，以验证上述算法的正确性。从图中可以看出，上述代码得到的估计谱图与MATLAB中自带的pburg函数得到估计谱是一致的，这就从侧面证明了上述代码的正确性。从而使用基本语言实现了Burg算法，这有利于采用其他语言进行算法的移植。图（a）使用自己编写的Burg算法得到的功率谱图（b）使用MATLAB中自带的pburg函数得到的功率谱图3.7 自编函数与库函数效果对比5 心电图谱估计的实现在上述Burg算法实现后，下面使用MIT/BIH数据库中的数据来求解其功率图3.8 100段数据前7秒的心电图图3.9 100段数据2-4秒正常心电信号的功率谱图3.10 100段数据5-7秒含有房性早搏（APB）心电信号的功率谱谱。以对比验证频谱心电图的功效，为医学上对心脏的研究提供一些参考。读取数据库中编号为100的数据，前5秒是标注为“1”的心电信号，表示正常心跳，提取两个周期的正常心电信号（720个），绘制为正常的心电功率谱，第五秒与第六秒之间的标注数字为“8”，表示房性早搏（APB），提取第五秒到第七秒之间的两个周期的数据，绘制为心电功率谱图。6 实验结果与分析从图3.9和图3.10上来看，两者之间还是存在一定的差别的，只是这种差别表现的不是很明显。在低频部分两者的区别稍微明显，从图中可以看出，病变的功率谱的第二个谱峰要比正常的谱峰突出。实际上，FCG作为一种处理结果，在医学上有较大的作用，是医生进行诊治的参考因素13141516，很多医学的参考资料都对其进行了讲解，由于非医学专业，这里不在进行医学方面的分析，下面给出101段的几组心电信号的功率谱图，作为对Burg算法的应用。从读取的数据上来看，这几组数据都是正常心跳。图3.11 101段数据10-12秒正常心电信号的功率谱图3.12 101段数据20-22秒正常心电信号的功率谱图3.11 101段数据22-24秒正常心电信号的功率谱 从图上来看，正常心电信号的功率谱基本上是相似的，特别是在低频部分，基本一致，而在高频部分则存在少量的差异，原因可能是存在高频的噪声干扰，而预处理没有将其滤除的缘故。四、结束语本文中涉及的实验全部都在MATLAB中实现，通过编写并运行相应的实验代码，实现了预定的实验目标。首先是分析了AR模型阶次的选择，其次是在MATLAB中实现了Burg算法，最后使用Burg算法对ECG信号进行了谱估计，并进行了分析。心电信号的谱估计作为一种无创性心电检测技术，其操作简便、快速、费用低、重复性好、病人容易接受。其在将来一定具有广阔的研究和应用价值，值得深入研究。参