《生物医学信号处理》实验指导书.doc

生物医学信号处理实验指导书2009 年 8 月目录实验一 随机信号的数字特征分析 ............................................................................1实验二 数字相关和数字卷积 ....................................................................................4实验三 维纳－霍夫方程 ............................................................................................7实验四 Yule-Walker 方程 ........................................................................................11实验一 随机信号的数字特征分析（一）实验目的了解随机信号的特征掌握随机信号的数字特征分析算法二）实验原理对于平稳各态遍历随机过程，可以用单一样本函数的时间平均代替集总平均，即通过测量过程的单一样本来估计信号的统计特征量。

样本均值： nixm1ˆ样本均方值： nixE122样本方差： nixixm122ˆˆ（三）实验内容和步骤用 matlab 编制程序，分析信号的数字特征，包括均值、方差、均方值、协方差可以使用 matlab 自带函数观察信号的直方图，粗略估计其概率分布信号 1：利用 MATLAB 中的伪随机序列产生函数 randn()产生的长 1000 点的序列；信号 2：实际采集的生物医学信号（脑电，心电等）四）思考题（1）改变每段数据长度，观察各段数字特征的分布情况数据长度对于数字特征估计值有什么样的影响？（2）观察伪随机序列，心电信号和脑电信号的直方图，它们之间是否相似？（3）通过同一数据分段估计数字特征，大致判断该数据是否可以看作广义平稳五）实验报告要求简述实验原理及目的；按实验要求编程分析所给信号的数字特征，记录运行结果；简要回答思考题附：参考程序% 选择信号类型并设定参数，产生信号 x(n)clear; clc;disp('请选择信号');disp('1 ---- 伪随机序列 randn()');disp('2 ---- 实际测量的心电信号');disp('3 ---- 实际测量的脑电信号');b = input('信号：');% 输入序号，产生相应信号switch bcase 1L = input('每段数据长度 L \n');N = input('数据共多少段 N \n');x = randn(1, L*N);case 2load ecgdata;display(['数据总长度',num2str(length(x)),' 点 ']);L = input('每段数据长度 L \n');N = input('数据共多少段 N \n');x = x(1:(N*L));case 3load eegdata; display(['数据总长度',num2str(length(x)),' 点 ']);L = input('每段数据长度 L \n');N = input('数据共多少段 N \n');x = x(1:(N*L));end% 估计信号的统计特征量Xmean = zeros(1,N); % 每段数据均值Xms = zeros(1,N); % 每段数据均方值Xvar = zeros(1,N); % 每段数据方差for k = 1:Nxs = x(((k-1)*L+1):(k*L));Xmean(k) = mean(xs);Xms(k) = std(xs).^2;Xvar(k) = var(xs);end% 显示n = 1:N;figure;subplot(2,2,1); stem(n,Xmean,'.'); title('mean');subplot(2,2,2); stem(n,Xms,'.'); title('mean square');subplot(2,2,3); stem(n,Xvar,'.'); title('variance');xlabel(['L=',num2str(L),' ','N=',num2str(N)]);subplot(2,2,4); hist(x,100); % 绘制数据直方图，观察信号大致的概率分布实验二 数字相关（一）实验目的熟悉数字相关的运算，初步在信号处理中应用相关技术。

二）实验原理相关可以从时域角度表现信号间的相似（关联）程度，是统计信号处理最基本的手段之一设有离散信号 x(n)和 y(n)，线性相关函数定义为nxymyr实际采集的信号总是有限长度，用有限的样本估计相关（自相关）函数 L,21,01ˆ0xNRmnx求和项总数不是 N 而是 N-|m|，因为当 n=N-|m|-1 时，n+| m|=N-1此时 xn+m已经到了数据边沿这种估计是渐进无偏估计和一致估计计算中，只要将其中一个序列反转，就可以用计算线性卷积的程序计算线性相关 nyxmrxy因此可以用 FFT 来加速相关运算，即对序列补零后，用循环相关计算线形相关，然后用循环卷积的快速算法计算循环相关，得到最终结果三）实验内容和步骤已知发射波形，利用相关技术，在有强背景噪声的情况下检测回波的延时和强度首先使用已知信号模版及其若干次衰减延迟生成仿真回波波形，然后与白噪声背景叠加，构造仿真信号然后计算模版与仿真信号的相关函数，判断回波位置及相对强度四）思考题尝试修改程序，包括改变仿真信号中模版的形状，噪声的强弱，噪声的类型（对白噪声滤波可以获得各种有色噪声） ，哪些因素会影响相关函数的结果？（五）实验报告要求简述实验原理及目的；按实验要求编程进行相关分析并调整程序参数，记录运行结果；简要回答思考题。

附：参考程序% 估计两个相似信号间的时间延迟np = 0:99;% p = sin(pi/5*np); % 正弦% p = exp(-0.06*np); % 指数衰减% p = sin(pi/5*np).*exp(-0.06*np); % 指数衰减正弦p = ones(size(np)); % 方波figure;subplot(2,2,1); plot(np,p);n = 0:1000;w = randn(size(n));s = zeros(size(n));A = 3; % 衰减系数s(100:199) = s(100:199)+A*p;s(500:599) = s(500:599)+A/3*p;s(800:899) = s(800:899)+A/3/3*p;x = s+w;figure;subplot(3,1,1); plot(n,w); title('Noise');subplot(3,1,2); plot(n,s); title('Signal');subplot(3,1,3); plot(n,x); title('Signal with Noise');p = [p,zeros(1,length(x)-length(p))]; % 如果要求归一化相关系数（相干系数） ，两个序列要同样长Rps = xcorr(s,p,'coeff');Rpw = xcorr(w,p,'coeff');Rpx = xcorr(x,p,'coeff');n2 = (n(1)-np(end)):(np(end)-n(1));figure;subplot(3,1,1); plot(Rpw); axis([0,2001,-1,1]); title('Rxy of p(n) and w(n)');subplot(3,1,2); plot(Rps); axis([0,2001,-1,1]); title('Rxy of p(n) and s(n)');subplot(3,1,3); plot(Rpx); axis([0,2001,-1,1]); title('Rxy of p(n) and x(n)');% subplot(3,1,1); plot(Rpw); title('Rxy of p(n) and w(n)');% subplot(3,1,2); plot(Rps); title('Rxy of p(n) and s(n)');% subplot(3,1,3); plot(Rpx); title('Rxy of p(n) and x(n)');实验三 维纳－霍夫方程（一）实验目的学习求解维纳-霍夫方程，寻找最小均方误差意义下的最优滤波器。

二）实验原理根据正交原理可以推导出维纳-霍夫方程，满足该方程的滤波器输出信号的估计值与信号在最小均方误差意义下最接近  ,0,opt LjmjRhjRmxxs根据滤波器的形式，维纳滤波器可以分为三种情况：非因果 IIR 型，因果IIR 型， FIR 型，对于实时性有要求的情况下用后两种形式 nwsxhnsyˆ图 1 维纳滤波器对于 FIR 型维纳滤波器，维纳-霍夫方程的形式为,010opt  NjmjRhjRNmxxs L或者写成矩阵形式 xsoptxH其中 Topt Nhh110Lxsxsxss RR021201xxx xxx RNRLMOM这样，如果信号和噪声的二阶统计特性已知，则易求解xsRH1维纳滤波的均方误差是 102NmxsoptsRhneE（三）实验内容和步骤已知信号的自相关函数和噪声的能量，编写程序求解维纳-霍夫方程，寻找最优滤波器编写程序仿真信号，噪声和观察波形，然后把观察信号通过滤波器得到的信号估计与原始信号比较，观察是否达到了去噪的目的。

选择不同信号（仿真信号，实际采集的心电，脑电信号） ，人工添加噪声，调整噪声的相对强度，观察滤波效果四）思考题观察实验结果，对于几种不同的信号，维纳滤波是否都取得了较好的效果？如果效果不好，试分析原因五）实验报告要求简述实验原理及目的；按要求编程求解维纳-霍夫方程，并对信号去噪，记录运行结果；简要回答思考题附：参考程序function [h,e] = WH(Rss,Rww,M)% 求解维纳 -霍夫方程的函数，其中 M 为信号的长度e1 = 10; e0 = 0;N = 0; % 给 e1,e0,N 赋初值% 以下循环的目的是找出 FIR 滤波器合适的阶数% 判据是当阶数增加而均方误差没有明显下降时，则认为阶数足够while abs(e0-e1)>1e-6 % e1 和 e0 不够接近则循环N = N+1;e0 = e1;Rxs = Rss(M:(M+N-1));Rxx = Rww(M:(M+N-1))+Rss(M:(M+N-1));R_xx = zeros(N);for j = 1:Nfor n = 1:NR_xx(j,n) = Rxx(abs(j-n)+1);endendh = inv(R_xx)*Rxs';e1 = Rss(M)-h'*Rxs';endN % 显示 N 的最终值e = e1;% 主程序clear; clc;M = input('信号的长度 M = ');n = 1:M;s = exp(-0.002*n).*sin(pi*n/50); % 仿真信号，可以自己生成，任意形式% load ecgdata; % 实际心电信号% s = x(1:M)';w = 0.4*randn(1,M); % 白噪声，系数代表噪声相对强度x = s+w; % 仿真信号Rss = xcor。