医学信号处理报告三

西安交通大学实验报告成绩实验日期：交报告日期：报告退发：教师审批签字：第页（共页）实验日期：交报告日期：报告退发：教师审批签字：第页（共页）2013年5月11日2013年5月11日（订正、重做）课程： 医学信号处理专业班号组别.姓名学号.同组者 实验名称： 陵机佶号分析一.实验要求1）产生1000点的白噪声信号，并计算它的均值、均方值、均方根值、方差。2）计算该1000点白噪声的自相关函数并作图显示白噪声与它的自相关函数。3）计算脑电信号的均值、均方值、均方根值、方差，计算脑电信号的自相关函数并作图显示脑电信号与它的自相关函数。（zuoye3_leeg.mat）4）计算含有噪声的心电信号的自相关函数并作图显示含有噪声的心电信号与它的自相关函数。（zuoye3_lzshecg.mat）二.实验原理自相关函数是信号在时域中特性的平均度量，它用来描述信号在一个时刻的取值与另一时刻取值的依赖关系，其定义式为：三.实验步骤与结果分析1）产生1000点的白噪声信号，并计算它的均值、均方值、均方根值、方差。根据题目的要求，利用randn产生均值为0,方差为1的随机白噪声。程序：clearcic%产生白噪声信号fs=1000;t=O:l/fs:l;x=randn(l,length(t));%计算均值meanx=mean(x)%计算均方值junfangx=sum(x(l,:).A2)/fs%计算均方根值junfanggx=junfangxA(l/2)%计算方差sum((x(l,:)-meanx).A2)/fs程序运行完成后，输出结果：均值：meanx=0.0305均方值：junfangx=0.9346均方根值：junfanggx=0.9667方差：fangchax=0.9337分析：和设计时均值为0,方差为1略有出入，但是从统计的角度看，是正确的。2)计算该1000点白噪声的自相关函数并作图显示白噪声与它的自相关函数根据题目的要求，求出白噪声的自相关函数，延迟为200。程序：%计算自相关函数lag=200;[c,lags]=xcorr(x,lag,'biased');%lag为M,lags为下标m%作图显示白噪声subplot(2,l,l)plot(t,x)xlabel('t');ylabel('x(t)');title。白噪声x(t)')%作图显示白噪声的自相关subplot(2,l,2)plot(lags/fs,c)xlabel('t');ylabel('Rxx(t)');title('白噪声的自相关')程序运行完成后，绘制出了白噪声的时域图像和他的自相关函数的图像，如图(1)：-E图(1)1000个点的白噪声和在maxlags=200时其自相关函数分析：从图中可以看出白噪声的自相关函数为一个在0时刻的冲击。3)计算脑电信号的均值、均方值、均方根值、方差，计算脑电信号的自相关函数并作图显示脑电信号与它的自相关函数加载eeg数据，loadzuoye3_leeg.mato其中ecgl变量存放的是原始脑电信号。根据要求计算各个参数的值。程序：clearcic%读入eeg数据load('E:\课程'医学信 乍业3\zuoye3_leeg.mat')fs=250;t=0:l/fs:10-l/fs;%计算均值meanx=mean(eegl)%计算均方值junfangx=sum(eegl(l,:).A2)/length(eegl)%计算均方根值junfanggx=junfangxA(l/2)%计算方差fangchax=sum((eegl(l,:)-meanx).A2)/length(eegl)%计算自相关函数[c,lags]=xcorr(eegl,'biased');%lag为M,lags为下标m%作图显示白噪声subplot(2,l,l)plot(t,eegl)xlabel('t');ylabel('eeg幅值')；title('脑电信号')%作图显示白噪声的自相关subplot(2,l,2)plot(lags/fs,c)xlabel('t');ylabel('eeg自相关')；title('白噪声的自相关')程序运行完成后，输出结果：均值：meanx=-0.0198均方值：junfangx=5.5061e-04均方根值：junfanggx=0.0235均方差：fangchax=1.5968e-04

程序运行完成后，绘制出了脑电信号的时域图像和它的自相关函数的图像，如图(3)：图（3）脑电信号的原始波形和其自相关函数图象4）计算含有噪声的心电信号的自相关函数并作图显示含有噪声的心电信号与它的自相关函数加载eeg数据，loadzuoye3_leeg.mato其中ecgl变量存放的是原始脑电信号。程序：clearcic%读入eeg数据load('E:\课程、医学信号处理\作业3\zuoye3_lzshecg.mat')；fs=250;t=0:l/fs:10-l/fs;%计算自相关函数[c,lags]=xcorr(zshecgl,'biased')；%lag为M,lags为下标m%作图显示白噪声subplot(2,l,l)plot(t,zshecgl)xlabel('t')；ylabel('ecg幅值')；title。含有噪声的心电信号')%作图显示白噪声的自相关subplot(2,l,2)plot(lags/fs,c)xlabel('t');ylabel('eeg自相关')；title('含有噪声的心电信号的自相关')程序运行完成后，绘制出了有噪声的心电信号和其自相关函数的图像，如图(4)：含H吟声的心电偌：；8禽含H吟声的心电偌：；8禽W号声的4屯■:;力■相关图(4)含有噪声的心电信号和其自相关函数分析:可以从图中看出由于心电信号由于噪声的原因使得自相关函数在0处的冲击变得很大。四.实验总结通过本次试验，使我对随机信号的分析以及自相关函数等概念有了更加深入的认识。熟练地掌握了用Matlab求取随机信号的各个参数的方法。第页(共页)课程：_医学信号处理实验日期：2013年5月11日专业班号硕2102组别交报告日期：2013年5月11日姓名—陈星学号3112322010报告退发：(订正、重做)同组者教师审批签字:西安交通大学实验报告成绩西安交通大学实验报告成绩实验名称： 生理佶号的功率谱估计一.实验要求1)已知随机信号x(n)=sin(2*pi*fl*n)+sin(2*pi*f*2n)+w(n).n=l,,,,,N,fl=0.05,12=0.12,w(n)为白噪声。①N=1024时，当M=1024,M=256时，求自相关法功率谱密度并作图。②当N=256,N=1024时，求周期图法功率谱密度并作图。③用welch法计算功率谱密度并作图。2)信号为作业1获取的心电信号，取一个周期，求周期图法功率谱密度并作图。3)信号为作业1获取的心电信号，取10个周期，每段一个周期，用不重叠分段法求平均周期图法功率谱密度并作图。4)信号为作业1获取的脑电信号，当N=256,N=1024时，求周期图法功率谱密度并作图。5)信号为作业1获取的脑电信号，每段长256时，用welch法计算功率谱密度并作图。生理信号抽样频率250Hz。二.实验原理D直接法直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N个观察值直接进行傅里叶变换，得到，然后取其幅值的平方，再除以N,作为对功率谱的估计。2)改进的周期图法将N点的观察值分成L个数据段，每段的数据为M,然后计算L个数据段的周期图的平均，作为功率谱的估计，以此来改善用N点观察数据直接计算的周期图的方差特性。根据分段方法的不同，又可以分为Welch法和Bartlett法。3)Welch法所分的数据段可以互相重叠，选用的数据窗可以是任意窗。三.实验步骤与结果分析1)分别利用自相关法，周期图法，和welch法计算随机信号的功率谱。程序：clearde%信号产生N=1024;n=l:l:N;fl=0.05;f2=0.12；x=sin(2*pi*fl*n)+sin(2*pi*f2*n)+randn(l,N);%相关函数估计c=xcorr(x,'biased');%谱估计M=Nlenth=length(x);px2047=abs(fft(c,2047));%谱估计M=256c2=[c(lenth-255:lenth-l),c(lenth:lenth+255)];px511=abs(fft(c2));%绘图p=fftshift(px2047);p2=fftshift(px511);ml=-1023:1023;m2=-255:255;figure(l)subplot(2,l,l)plot(ml*N/2047,20*logl0(p))title('M=1024')subplot(2,l,2)plot(m2*N/511,20*logl0(p2))title('M=512')

mm图（1）M=1024和M=256的自相关法功率谱程序：clearcic%信号产生N=1024;n=l:l:N;fl=0.05;f2=0.12;x=sin(2*pi*fl*n)+2*sin(2*pi*f2*n)+randn(l,N);%谱估计N=1024pxl024=abs(fft(x,1024)),A2/1024;%N=256px256=abs(fft(x,256)).A2/256;%绘图p=fftshift(pxl024);p2=fftshift(px256);ml=-511:512;m2=-127:128;figure(l)subplot(2,l/l)plot(ml,20*logl0(p))title('N=1024')subplot。,1,2)plot(m2,20*logl0(p2))title('N=256')程序运行结果如图(2)：图(2)N=1024和N=256的周期图法功率谱程序：clearde%信号产生N=1024;n=l:l:N;fl=0.05;f2=0.12;xn=sin(2*pi*fl*n)+sin(2*pi*f2*n)+randn(l,N);%welch法海宁窗256个数据为一段重合128个数据L=7w=hanning(256),;pxx=(abs(fft(w.*xn(l:256))).A2+abs(fft(w.*xn(129:384))).A2+abs(fft(w.*xn(257:512))).A2+abs(fft(w.*xn(385:640))).A2+abs(fft(w.*xn(513:768))).A2+abs(fft(w.*xn(641:896))).A2+abs(fft(w.*xn(769:1024))).A2)/(256*7);p=fftshift(pxx);plot(-127:128,10*logl0(p));title('welch法')；程序运行结果如图(3)：法法图(3)welch法功率谱2)对心电信号取…个周期，用周期图法求功率谱密度并作图程序：clearcic%读入ecg数据load('E:\课程'医学信号处理、作业3\zuoye3_lzshecg.mat')；N=2500fs=250;m=-1249:1250;px2500=abs(fft(zshecgl,2500)).A2/N;p=fftshift(px2500);plot(m*fs/N,10*logl0(p));title('心电信号周期图法功率谱(一个周期)')；程序运行结果如图(4)：心电俱；青明阴法正章谓＜7WW!）图（4）心电信号一个周期周期图法功率谱3)对心电信号取10个周期，每段…个周期，用不重叠分段法求平均周期图法功率谱密度并作图。程序：clearcic%读入ecg数据load('E:\课程'医学信号处理、作业3\zuoye3_lzshecg.mat')；N=2500;fs=250;m=-124:125%t=0:l/fs:10-l/fs;figure(2)pxx=(abs(fft(zshecgl(l:250))).A2+abs(fft(zshecgl(251:500))).A2+abs(fft(zshecgl(501:750))).A2+abs(fft(zshecgl(751:1000))).A2+abs(fft(zshecgl(1001:1250))).A2+abs(fft(zshecgl(1251:1500))).A2+abs(fft(zshecgl(1501:1750))).A2+abs(fft(zshecgl(1751:2000))).A2+abs(fft(zshecgl(2001:2250))).A2+abs(fft(zshecgl(2251:2500))).A2)/250*10;p=fftshift(pxx);plot(m*fs/250,10*logl0(p));title('心电信号周期图法功率谱(10个周期)')；

程序运行完成后，结果如图（5）：图（5）心电信号一个周期周期图法功率谱分析：可以看出10个周期的功率谱密度图像要比一个周期的平滑许多。4）求脑电信号当N=256,N=1024时的周期图法功率谱密度并作图。程序：clearcic%读入eeg数据load('E:\课程'医学信号处理、作业3\zuoye3_leeg.mat')fs=250%谱估计N=1024pxl024=abs(fft(eegl,1024)).A2/1024;%N=256px256=abs什ft(eegl,256)).△2/256;%绘图p=fftshift(pxl024);p2=fftshift(px256);ml=-5U:512;m2=-127:128；figure(l)subplot(2,l,l)plot(ml*fs/1024,20*logl0(p))title('N=1024周期图法eeg功率谱')subplot(2,l,2)plot(m2*fs/256,20*logl0(p2))title('N=256周期图法eeg功率谱')程序运行完成后，结果如图(6)：图（6）N=1024和N=256时的脑电信号周期图法功率谱5）求脑电信号每段长256时welch法的功率谱密度并作图。程序:clearcic%读入eeg数据load('E:\课程'医学信号处理\作业3\zuoye3_leeg.mat')fs=250;xn=eegl;%welch法海宁窗256个数据为一段重合128个数据L=18w=hanning(256),;pxx=(abs(fft(w.*xn(l:256))).A2+abs(fft(w.*xn(129:384))),A2+abs(fft(w.*xn(257:512))).A2+abs(fft(w.*xn(385:640))).A2+abs(fft(w.*xn(513:768))).A2+abs(fft(w.*xn(641:896))).A2+abs(fft(w.*xn(769:1024))).A2+abs(fft(w.*xn(897:1152))).A2+abs(fft(w.*xn(l025:1280))).A2+abs(fft(w.*xn(1153:1408))).A2+abs(fft(w.*xn(1281:1536))).A2+abs(fft(w.*xn(1409:1664))).A2+abs(fft(w.*xn(1537:1792))).A2+abs(fft(w.*xn(1665:1920))).A2+abs(fft(w.*xn(1793:2048))).A2+abs(fft(w.*xn(1921:2176))).A2+abs(fft(w.*xn(2049:2304))).A2+abs(fft(w.*xn(2177:2432))).A2)/(256*7);p=fftshift(pxx);plot((-127:128)*fs/256,10*logl0(p));title('welch法eeg功率谱')；程序运行完成后，结果如图(7)：图（7）N=1024和N=256时的脑电信号周期图法功率谱四.实验总结通过本次的实验，使我对自相关法、周期图法和welch法求功率谱有了较为深刻的理解，同时使我对利用matlab进行功率谱估计方法有了比较深入地认识，对以后的学习具有非常大的帮助。实验日期：交报告日期：报告退发：教师审批签字：第页（共页）2013年实验日期：交报告日期：报告退发：教师审批签字：第页（共页）2013年5月11日2013年5月H日（订正、重做）实验名称:脑电佶号的AR参数模型谱佶计西安交通大学实验报告课程：医学信号处理 专业班号硕2102组别姓名陈星学号3112322010同组者 一.实验要求1）已知32点观测值x（n）=[0.4282 1.14541.55971.89941.6854 2.30752.4679 1.97901.60631.2804 -0.20830.0577 0.02060.3572 1.65720.74881.6666 1.98302.6914 1.25211.8691 1.68550.62420.1763 1.34900.6955 1.2941-0.6177]A①计算自相关序列氏-团）1.0475 0.4319并作图显示。0.03120.5802②用求出的自相关序列Rx（⑼来估计3阶AR模型参数a3⑴，a3（2）,a3（3）及A2方差0③求该观测序列的3阶AR模型功率谱估计并作图显示。④若要求Ep<=0.3783,求AR模型的阶数并计算功率谱，要求作图显示。2）信号为作业1获取的脑电信号，求3阶AR模型功率谱估计并作图显示。二.实验原理AR模型谱估计原理：该模型除bo=l外，勺（'=1八1%）的取值均为零。其系统函数为“（z）=-—k=\x(n)=w(n)-2\akx(n-k)其差分方程为 *=i ,可知其功率谱密度为' A⑵A—)所以可推出其功率谱为2P(°)=嘘——AR模型又称为p阶自回归模型，简称AR模型。其系统函数只有极点，没有零点,因此也称为全极点模型。三.实验步骤与结果分析A1)计算32点观测值的自相关序列，然后利用公式计算3阶AR模型参数a3⑴，A A Aa3(2),a3(3)及方差继而可以求出3阶AR模型功率谱估计并作图显示，最后求解在限制条件下Ep<=0.3783,求出AR模型的阶数并计算功率谱，作图显示。将AR模型计算功率谱写成函数puguji.m。程序：functionP=puguji(a,p/E);%a为AR模型参数，p为阶数,E为方差w="pi:0.01:pi;forn=l:l:length(w)sum=l;fork=l:l:psum=sum+a(k+l)*exp(-j*w(n)*k);endP(n)=E/abs(sum)A2;endplot(w/2/pi,P)clearcicxn=[0.42821.14541.55971.89941.68542.30752.46791.97901.60631.2804-0.20830.05770.02060.35721.65720.74881.66661.98302.69141.25211.86911.68550.62420.17631.34900.69551.29411.04750.43190.03120.5802-0.6177];%求自相关函数并作图[rx,lags]=xcorr(xn,'biased');subplot(3,l,l)plot(rx)title('自相关函数')；%m=0e0=rx(32);a0=l;%m=lall=-rx(33)/e0;el=eO*(l-allA2);%m=2a22=-(rx(34)+all*rx(33))/el;a21=all*(l+a22);e2=el*(l-a22A2);%m=3a33=-(rx(35)+a21*rx(34)+a22*rx(33))/e2;a31=a21+a33*a22;a32=a22+a33*a21;e3=e2*(l-a33A2);%显示所要求的参数a31=a31a32=a32a33=a33e3=e3%求3阶功率谱并作图A=[l,a31,a32,a33];subplot(3,l,2)puguji(A,3,e3);title('3阶AR模型谱估计')；

%利用matlab的函数来求解阶数nE=l;i=4;while(E>03783)[aE]=aryule(xn,i);i=i+l;endi=i-l%求i阶Ar功率subplot(3,l,3)puguji(a,i,E);title('12阶AR模型谱估计')程序运行结果：Aa3(l)：a31=-0.6984Aa3(2)：a32=-0.2748Aa3(3)：a33=0.0915A方差6？：e3=0.4678在Ep〈=0.3783最小阶数为:i=12图(1)中做出了序列的自相关函数，3阶和12阶的AR模型功率谱:图(1)AR模型功率谱估计2）信号为脑电信号，求3阶AR模型功率谱估计并作图显示。程序：cleardeload（'E:\课程'医学信号处理'作业3\zuoye3_leeg.mat'）xn=eegl;[a,Ep]=aryule（xn,3）%求3阶功率谱并作图puguji（a,3,Ep）;title（'3阶AR模型谱估计，）；程序运行结果如图（2）：图（1）脑电信号的3阶AR模型功率谱估计四.实验讨论经典谱估计法可以利用FFT计算，因而有计算效率高的优点，在谱分辨力要求不是太高的地方常用这种方法。但频率分辨率地是经典谱估计的一个无法回避的缺点。如周期图法在计算中把观测到的有限长的N个数据以外的数据认为是零,而BT法仅利用N个有限的观测数据作自相关函数估计，实质上也就是假设除已知数据外的自相关函数全为零，这些显然都是与事实不符的。为了克服以上缺点,人们提出了平均，加窗平滑等方法，在一定程度上改善了经典谱估计的性能。但是，经典谱估计，始终无法解决，频率分辨率与谱估计稳定性之间的矛盾，特别是在数据记录长度比较短时，这一矛盾尤其突出。现代谱估计理论也就是在这种背景下产生的，以1967年Burg提出的最大燧谱分析法为代表的现代谱估计法，不认为在观察到的N个数据以外的数据全为零。因此克服了经典法的这个缺点，提高了谱估计的分辨率。后来发现线性预测自回归模型法（简称AR模型法）与Burg的最大烯谱分析法是等价的，它们都可归结为通过Yule-Walker方程求解自回归模型的系数问题。五.实验总结通过本次试验，使我对AR模型有了更加深入的认识和了解，让我对现代谱估计法的优势有了更深层次的认识，并掌握了利用matlab来计算AR模型谱估计的参数的方法对今后的学习大有裨益。西安交通大学实验报告第页(共页)课程：_医学信号处理实验日期：2013年5月11日专业班号硕2102组别交报告日期：2013年5月11日姓名—陈星学号3112322010报告退发：(订正、重做)同组者_教师审批签字：实验名称： 自适应滤波春设计与仿真一.实验要求1)噪声抵消器图⑴2)按照以上噪声抵消器框图设计自适应滤波器。其中s(k)为心电信号，如ecg0410,x(k)=sin(2n50nT)为50Hz干扰，T=l/250o3)使用以下迭代公式求滤波器系数h(n)及e(k)0ep)ad(n)~戈，巾He-H"+/xe(i)写成分量形式：hn+,(k)=h„(k)+ux(n-k)e(n) k=0,…M可取M=5其中步长入max,入max为矩阵R的最大特征值，口可取0.1左右。初始值可设h°(k)=Ok=0,-Mo当n足够大时，有e(n)=s(n)。要求显示观测数据d(k)与滤波器输出e(k)的波形。可参考matlab帮助―>Demos—ToolboxesfFilterDesignfAdaptiveFiIter中的有关例子。二.实验原理D维纳滤波从连续的(或离散的)输入数据中滤除噪声和干扰以提取有用信息的过程称为滤波，而相应的装置称为滤波器。根据滤波器的输出是否为输入的线性函数，可将它分为线性滤波器和非线性滤波器两种。滤波器研究的一个基本课题就是：如何设计和制造最佳的或最优的滤波器。所谓最佳滤波器是指能够根据某•最佳准则进行滤波的滤波器。20世纪40年代，维纳奠定了关于最佳滤波器研究的基础。即假定线性滤波器的输入为有用信号和噪声之和，两者均为广义平稳过程且知它们的二阶统计特性，维纳根据最小均方误差准则(滤波器的输出信号与需要信号之差的均方值最小)，求得了最佳线性滤波器的参数，这种滤波器被称为维纳滤波器。在维纳研究的基础上，人们还根据最大输出信噪比准则、统计检测准则以及其他最佳准则求得的最佳线性滤波器。实际上，在一定条件下，这些最佳滤波器与维纳滤波器是等价的。因而，讨论线性滤波器时，一般均以维纳滤波器作为参考。维纳滤波理论用于解决最小均方误差下的线性滤波问题。设接收到(或观测到)的信号为随机信号x@)■项)+4) ⑴其中s(t)是未知的实随机信号，n(t)是噪声。要设计的线性滤波器，其冲击响应为h(t,t),输入为x(t),输出为W),即喇・04小》 ⑵令为估计误差。冲击响应h(t,C按最小均方误差准则确定，即h(t,T)必须满足使

达到最小。根据最小均方误差估计的正交条件，有以下关系成立TOC\o"1-5"\h\z啡(0-匚坤加网则｝.0 -84-48 ⑷令凡(5)见1)・小仲崎 ⑹则有口k款&。< -8$- ⑦上述方程通常称为非平稳随机过程条件下的维纳-霍甫(Wiener-Kolmogorov)积分方程。特别当x(t),s(t)均为广义(或宽)平稳随机信号，而滤波器是线性时不变系统的情况下，x(t)与s(t)必为联合平稳，式(7)可写为入0口1咪(，-印・ (8)令”(■工,则有《(切■口(砒y°旬3・喇叫(动此处，“*”号表示卷积，对上式两边取Fourier变换，可得TOC\o"1-5"\h\zs上WWS) (io)3霸 (11)对于因果线性系统，有»(/)-fh(r-rX(r)dr (⑵采用完全相同的分析方法，推得因果平稳维纳-霍甫积分方程如下(jy)-rii(4)R.(y-4)(U (13)则郡 (14)其中$<•)・£(•)£(•