资源描述：
《生物医学信号的参数建模及功率谱分析》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、《生物医学信号处理》实习报告学生姓名：学号：实验室名称：项目名称：生物医学信号的参数建模及功率谱分析项目内容：1）编写函数求解Y-W方程，运用函数x=filter(b,a,w)学习随机信号的建模；2）与Matlab的函数进行比较，检验正确性：[aE]＝aryule(x,p)；L-D算法——用不同的阶数进行比较。[aE]=arburg(x,p)；burg算法——用不同的阶数进行比较。3）根据提供的RR间期数据建立AR模型，阶数大于10不等；4）运用两种方法计算这两组RR间期数据的功率谱、LF/HF；原理（写出具体的计算公
2、式）一．AR参数模型参数模型谱估计作为现代功率谱估计的重要研究方向，主要包括AR模型、MA模型和ARMA模型。因AR模型参数的精确估计可通过解一组线性方程得到，计算工程量较小，最为常用。参数模型法的建模思路：1.假定所研究的过程x(n)是由一个输入序列u(n)激励一个线性系统H(z)的输出；2.有已知的x(n),或其自相关函数来估计H(z)的参数；3.由H(z)的参数来估计x(n)的功率谱。其中：1.1AR参数建模x(n)、u(n)间的关系可表述为下式的差分方程形式：(1)另外，(2)(1)式、(2)式两边z变换，并假
3、定，得：(3)其中：为保证H(z)为稳定且最小相位的系统，的零点均在单位圆内。假定u(n)为方差为的白噪声序列，由随机信号通过线性系统的理论可知，输出序列x(n)的功率谱：（4）(1)式中，当全为0，则(1)式、(3)式、(4)式分别变为：（5）（6）此三式给出的模型即为AR模型。AR模型是一个全极点模型，其当前输出时现在输入与过去p个输入的加权和。1.2AR模型的构建将1.1中(5)式两边同乘以x(n+m),并求均值，得：（7）因u(n)为白噪声，由(2)式，得：（8）如果是因果序列，即时，，则上式(8)可简化为：（
4、9）根据变换定义，，因此，综合(8)式、(9)式，有：所以式还可以继续简化为：（10）基于自相关函数的偶对称性()。上式可写成矩阵形式：（11）该式即为模型的方程。综上分析，一个p阶的AR模型共p+1个参数，即，仅需知道x(n)前p+1个自相关函数，由(9)、(10)、(11)式的线性方程即可求得这p+1个参数，将他们带入功率谱方程中，即可求出x(n)的功率谱。1.3AR模型阶数的选择阶次p未知，先选定稍大的值，在递推过程中确定。例如在使用Levinson递推时，给出由低阶到高阶的每一组参数，且模型的最小预测误差功率递
5、减。当达到所指定的希望值，或不再发生变化时，此时的阶数即最正确的阶次。最终预测误差准则：信息论准则：当FPE(k)和AIC(k)在某一k处取得极小值，此时k即为最合适的阶次p.二、随机信号功率谱分析方法2.1平均周期图法周期图法即把随机序列的个观测数据视为一能量有限的序列，直接计算的离散傅里叶变换，得，然后再取其幅值的平方，并除以，作为序列的真实功率谱的估计。平均周期图基本原理：观测一个随机变量，得到L组独立记录数据，用每一组数据求其均值，然后将L个均值加起来求均值，这样得到的均值，其方差是原始周期图法方差的1/L.平
6、均周期估计器定义：其中：为保证获得最大分辨率，L应尽可能选大一些，这里选定L=N-1。2.2自相关法自相关法的理论基础是W-H定理。先由序列估计出自相关函数，再对进行傅里叶变换：其中，为平稳随机过程的自相关函数。从而得到的功率谱估计。其中自相关函数的估计有两种方法：无偏自相关函数估计和有偏自相关函数估计，BT法功率谱估计采用的是有偏自相关函数估计法。编写的源程序：一．AR参数建模算法的实现代码1.直接运用观测值的自相关序列来求解方程，自定义的方程MATLAB实现代码如下：function[a,E]=YW(x,p)A=x
7、corr(x,x);%求x的自相关N=length(A);b=int16((N+1)/2);%b取R(0)在矩阵A中所对应的的下标fori=1:p+1%求解p阶的自相关矩阵k=b;%根据自相关序列左右对称求解forj=1:p+1v(i,j)=A(abs(k));k=k+1;endb=b-1;ends=zeros(p+1,1);%s初始化为零矩阵s(1)=1;%s(1)赋值为1a=(v^-1)*s;%解矩阵方程a=a/a(1);E=1/a(1);end2.主程序，其中包括与Matlab函数aryule(L-D算法)、ar
8、burg(burg算法)进行比较:clear;clcA=textread('NSR_rr1h.txt','%*s%*s%n%*s%*s','headerlines',1);%数据p=100;%阶数figure(1)subplot(2,2,1);plot(A);title('原始RR期间信号');gridon;xlim([0,500
