生物医学信号的参数建模及功率谱分析.doc

生物医学信号处理实习报告学生姓名：学号：实验室名称：项目名称：生物医学信号的参数建模及功率谱分析项目内容：1） 编写函数求解Y-W方程，运用函数x=filter(b,a,w)学习随机信号的建模；2） 与Matlab的函数进行比较，检验正确性：a Earyule(x,p)；L-D算法用不同的阶数进行比较。a E=arburg(x,p)；burg算法用不同的阶数进行比较。3） 根据提供的RR间期数据建立AR模型，阶数大于10不等；4） 运用两种方法计算这两组RR间期数据的功率谱、LF/HF；原理（写出具体的计算公式）1 AR参数模型参数模型谱估计作为现代功率谱估计的重要研究方向，主要包括AR模型、MA模型和ARMA模型。因AR模型参数的精确估计可通过解一组线性方程得到，计算工程量较小，最为常用。参数模型法的建模思路：1. 假定所研究的过程x(n)是由一个输入序列u(n)激励一个线性系统H(z)的输出；2. 有已知的x(n),或其自相关函数来估计H(z)的参数；3. 由H(z)的参数来估计x(n)的功率谱。其中：1.1 AR参数建模x(n)、u(n)间的关系可表述为下式的差分方程形式： (1)另外， (2)(1) 式、(2)式两边z变换，并假定，得： (3) 其中： 为保证H(z)为稳定且最小相位的系统，的零点均在单位圆内。 假定u(n)为方差为的白噪声序列，由随机信号通过线性系统的理论可知，输出序列x(n)的功率谱： （4）(1)式中，当全为0，则(1)式、(3)式、(4)式分别变为： （5） （6） 此三式给出的模型即为AR模型。AR模型是一个全极点模型，其当前输出时现在输入与过去p个输入的加权和。1.2 AR模型的构建将1.1中(5)式两边同乘以x(n+m),并求均值，得： （7）因u(n)为白噪声，由(2)式，得： （8）如果是因果序列，即时，则上式(8)可简化为： （9）根据变换定义，因此，综合(8)式、(9)式，有：所以式还可以继续简化为： （10）基于自相关函数的偶对称性()。上式可写成矩阵形式： （11）该式即为模型的方程。综上分析，一个p阶的AR模型共p+1个参数，即，仅需知道x(n)前p+1个自相关函数，由(9)、(10)、(11)式的线性方程即可求得这p+1个参数，将他们带入功率谱方程中，即可求出x(n)的功率谱。1.3 AR模型阶数的选择阶次p未知，先选定稍大的值，在递推过程中确定。例如在使用Levinson递推时，给出由低阶到高阶的每一组参数，且模型的最小预测误差功率递减。当达到所指定的希望值，或不再发生变化时，此时的阶数即最正确的阶次。最终预测误差准则：信息论准则：当FPE(k)和AIC(k)在某一k处取得极小值，此时k即为最合适的阶次p.2、 随机信号功率谱分析方法2.1平均周期图法周期图法即把随机序列的个观测数据视为一能量有限的序列，直接计算的离散傅里叶变换，得，然后再取其幅值的平方，并除以，作为序列的真实功率谱 的估计。平均周期图基本原理：观测一个随机变量，得到L组独立记录数据，用每一组数据求其均值，然后将L个均值加起来求均值，这样得到的均值，其方差是原始周期图法方差的1/L.平均周期估计器定义：其中： 为保证获得最大分辨率，L应尽可能选大一些，这里选定L=N-1。2.2自相关法自相关法的理论基础是W-H定理。先由序列估计出自相关函数，再对进行傅里叶变换：其中，为平稳随机过程的自相关函数。从而得到的功率谱估计。其中自相关函数的估计有两种方法：无偏自相关函数估计和有偏自相关函数估计，BT法功率谱估计采用的是有偏自相关函数估计法。编写的源程序：1 AR参数建模算法的实现代码 1.直接运用观测值的自相关序列来求解方程，自定义的方程MATLAB实现代码如下：function a,E=YW(x,p) A=xcorr(x,x); %求x的自相关N=length(A);b=int16(N+1)/2); %b取R(0)在矩阵A中所对应的的下标for i=1:p+1 %求解p阶的自相关矩阵 k=b; %根据自相关序列左右对称求解 for j=1:p+1 v(i,j)=A(abs(k); k=k+1; end b=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);end 2.主程序，其中包括与Matlab函数aryule(L-D算法)、arburg(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期间信号);grid on;xlim(0,500);subplot(2,2,2);a1 E1=aryule(A,p);B=filter(-1*a1(2:end),1,A);plot(B);title(L-D算法建模模拟);grid on;xlim(0,500);subplot(2,2,3);a2 E2=arburg(A,p);C=filter(-1*a1(2:end),1,A);plot(C);title(burg算法建模模拟);grid on;xlim(0,500);subplot(2,2,4);a,E=YW(A,p);D=filter(-1*a(2:end),1,A);plot(D);title(Y-W方程建模模拟);grid on;xlim(0,500); figure(2)N=length(A);subplot(3,1,1);plot(1:N),A,r,(1:N),B,b);title(原信号与L-D算法建模信号比较);axis(0,2000,0.4,0.9);grid on;subplot(3,1,2);plot(1:N),A,r,(1:N),C,b);title(原信号与burg算法建模信号比较);axis(0,2000,0.4,0.9);grid on;subplot(3,1,3);plot(1:N),A,r,(1:N),D,b);title(原信号与Y-W算法建模信号比较);axis(0,2000,0.4,0.9);grid on;2 随机信号功率谱的实现代码1. 平均周期图法functionk=Bartlett(Fs,A,B) xn1=A;xn2=B;n=0:1/Fs:1;nfft=1024;window=boxcar(length(n); %矩形窗noverlap=0; %数据无重叠p=0.9; %置信概率Pxx1,Pxxc1=psd(xn1,nfft,Fs,window,noverlap,p);Pxx2,Pxxc2=psd(xn2,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx1=10*log10(Pxx1(index+1);plot_Pxx2=10*log10(Pxx2(index+1);figure(1);subplot(2,1,1);plot(k,plot_Pxx1);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(正常心律Bartlett法);grid on;subplot(2,1,2);plot(k,plot_Pxx2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(心律不齐Bartlett法);grid on;2. 自相关法functionk=zixiangguan(Fs,A,B) xn1=A;xn2=B;n=0:1/Fs:1;nfft=1024;index=0:round(nfft/2-1);k=index*Fs/nfft;%正常心律序列cxn1=xcorr(xn1,biased); %计算正常心律序列的自相关函数CXk1=fft(cxn1,nfft);Pxx1=abs(CXk1);plot_Pxx1=10*log10(Pxx1(index+1);%心律不齐序列cxn2=xcorr(xn2,biased); %计算心律不齐序列的自相关函数CXk2=fft(cxn2,nfft);Pxx2=abs(CXk2);plot_Pxx2=10*log10(Pxx2(index+1); subplot(2,1,1);plot(k,plot_Pxx1);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(正常心律自相关法);grid on;subplot(2,1,2);plot(k,plot_Pxx2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(心律不齐自相关法);grid on;结论（画出要求的图形）1. AR参数建模图1-1中4图分别表示当阶次p=100时，采用L-D算法、burg算法和Y-W方程算法建模模拟出的RR间期信号：图1-1 为反映三种算法的预测效果，分别给出三种算法和原信号的比较图，如图1-2所示，可知当阶次为100时，三种算法的建模预测均较好：图1-22. 随机信号功率谱计算： 图2-3和2-4分别显示了采用平均周期图法和自相关法的功率谱估计结果：图2-3图2-4总结1. 关于AR建模中的阶次选择从MATLAB运行结果可以看出：.L-D算法的AR模型功率谱曲线基本可以区分间隔比较远的频率点，不会出现假峰现象，而且L-D递推算法估计出的时频分布相比之下性能优越于其他谱估计方法。但是对于间距比较近的频率点的分辨率和检出能力不是很好，这时就需要提高AR模型的阶次；.Burg算法中的功率谱曲线的分辨率要明显好于L-D算法，在AR模型阶次比较低的情况下也能分辨出间距比较近的频率，得到的谱线也比较光滑。且在谱估计分辨率高阶情况下没有频率分量的区域比较平伏。但是若模型阶次选的太低，得到的谱太平滑，可能反映不出谱峰；若选的过大，可能会出现虚假的峰值，也就是说谱估计的质量是受到模型阶次的影响。.此外，经实验验证，根据Y-W方程计算估计的功率谱曲线与L-D算法结果基本相同。综上，阶次决定了模型功率谱估计的优劣，当阶数小时，三种算法结果大致相同，当阶数比较大时，结果有差异，更能体现出burg算法相对L-D算法的优越性。2. 关于两种功率谱估计方法的比较说明平均周期图法：基于周期图法估计存在的功率不平滑、估计协方差大的确定，我们改进的平均周期图法估计得出的功