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

生物医学信号处理

实验指导书

2009年8月

目录

实验一 随机信号的数字特征分析............................................................................ 1 实验二 数字相关和数字卷积.................................................................................... 4 实验三 维纳－霍夫方程............................................................................................ 7 实验四

Yule-Walker方程 ........................................................................................ 11 实验一 随机信号的数字特征分析

（一）实验目的

了解随机信号的特征。

掌握随机信号的数字特征分析算法。

（二）实验原理

对于平稳各态遍历随机过程，可以用单一样本函数的时间平均代替集总平均，即通过测量过程的单一样本来估计信号的统计特征量。

1n?x??xi 样本均值：mni?11n2样本均方值：Ex??xi

ni?1??2n1n????xi?m?x?2 样本方差：?ni?12x

（三）实验内容和步骤

用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 b case 1

L = input('每段数据长度 L \\n'); N = input('数据共多少段 N \\n'); x = randn(1, L*N); case 2

load ecgdata;

display(['数据总长度',num2str(length(x)),'点']); L = input('每段数据长度 L \\n'); N = input('数据共多少段 N \\n'); x = x(1:(N*L)); case 3

load 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:N

xs = 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)，线性相关函数定义为

rxy?m??n????x?n?y?n?m?

?实际采集的信号总是有限长度，用有限的样本估计相关（自相关）函数

??m??1RxNN?m?1?xxn?0nn?mm?0,?1,?2,?

求和项总数不是N而是N-|m|，因为当n=N-|m|-1时，n+|m|=N-1。此时xn+m已经到了数据边沿。这种估计是渐进无偏估计和一致估计。

计算中，只要将其中一个序列反转，就可以用计算线性卷积的程序计算线性相关

rxy?m??x?n??y??n?

因此可以用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)');

实验三 维纳－霍夫方程

（一）实验目的

学习求解维纳-霍夫方程，寻找最小均方误差意义下的最优滤波器。

（二）实验原理

根据正交原理可以推导出维纳-霍夫方程，满足该方程的滤波器输出信号的估计值与信号在最小均方误差意义下最接近。

Rxs?j??m????h?m?R?j?m?optxx??j???,?,0,?,??

根据滤波器的形式，维纳滤波器可以分为三种情况：非因果IIR型，因果IIR型，FIR型，对于实时性有要求的情况下用后两种形式。

x?n??s?n??w?n? h?n???n?y?n??s 图1 维纳滤波器

对于FIR型维纳滤波器，维纳-霍夫方程的形式为

Rxs?j???hopt?m?Rxx?j?m?m?0N?1j?0,1,?,N?1

或者写成矩阵形式

RxxHopt?Rxs其中

Hopt??h?0?h?1??h?N?1??T

TRxs??Rxs?0?Rxs?1??Rxs?N?1??

Rxx?1??Rxx?0??R?1?Rxx?0?xxRxx???????Rxx?N?1?Rxx?N?2??Rxx?N?1???Rxx?N?2????????Rxx?0??

这样，如果信号和噪声的二阶统计特性已知，则易求解

1H?R?xxRxs

维纳滤波的均方误差是

Ee?n??Rss?0???hopt?m?Rxs?m?

2m?0??N?1

（三）实验内容和步骤

已知信号的自相关函数和噪声的能量，编写程序求解维纳-霍夫方程，寻找最优滤波器。

编写程序仿真信号，噪声和观察波形，然后把观察信号通过滤波器得到的信号估计与原始信号比较，观察是否达到了去噪的目的。

选择不同信号（仿真信号，实际采集的心电，脑电信号），人工添加噪声，调整噪声的相对强度，观察滤波效果。

（四）思考题

观察实验结果，对于几种不同的信号，维纳滤波是否都取得了较好的效果？如果效果不好，试分析原因。

（五）实验报告要求

简述实验原理及目的；

按要求编程求解维纳-霍夫方程，并对信号去噪，记录运行结果； 简要回答思考题。

附：参考程序

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:N for n = 1:N

R_xx(j,n) = Rxx(abs(j-n)+1); end end

h = inv(R_xx)*Rxs'; e1 = Rss(M)-h'*Rxs'; end

N % 显示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 = xcorr(s,s); %估计信号自相关函数

Rww = xcorr(w,w); % 估计噪声自相关函数 [h,e] = WH(Rss,Rww,M);

ss = filter(h,1,x); %用维纳滤波器滤波 figure; subplot(2,2,1); plot(n,s); title('信号'); subplot(2,2,2); plot(n,w); title('噪声'); subplot(2,2,3); plot(n,x); title('观测值'); subplot(2,2,4); plot(n,ss); title('信号估计'); figure; plot(n,ss-s); title('估计误差');

实验四 Yule-Walker方程

（一）实验目的

学习求解Yule-Walker方程，建立随机信号的AR模型。

（二）实验原理

随机信号可以看作是由当前激励白噪声w(n)以及若干次以往信号x(n-k)的线性组合产生，即所谓自回归模型（AR模型）

x?n??w?n???akx?n?k?k?1pH?z??11??akz?kk?1p

模型参数满足Yule-Walker方程

Rxx?m????akRxx?m?k?m?0

k?1p矩阵形式

R??1??R?0??R?1?R?0???????R?p?1?R?p?2??R?1?p???a1??R?1???a??R?2???R?2?p??2??????? ??????????????a???R?0???Rp?p????求解Yule-Walker方程，就可以得到AR模型系数

当模型阶次较大时，直接用矩阵运算求解的计算量大，不利于实时运算。利用系数矩阵的特性，人们提出了如L-D算法等快速算法。

（三）实验内容和步骤

编写求解Yule-Walker方程的程序，并对实际生理信号（例如脑电）建立AR模型。

对同一数据，使用matlab信号处理工具箱自带函数aryule计算相同阶数AR模型系数，检验程序是否正确。

用伪随机序列（白噪声）驱动AR模型，观察输出是否与真实脑电信号相似，对比真实信号与仿真信号的功率谱。

（四）思考题

对EEG建模后，该模型产生的信号是否能反映EEG信号的特征？

（五）实验报告要求

简述实验原理及目的；

按要求编程求解Yule-Walker方程，并对脑电信号建立AR模型； 使用白噪声驱动生成仿真脑电； 简要回答思考题。

附：参考程序

% 实验四 Yule-Walker方程 clear; clc; load eegdata;

x = x(1:1024); % 长度可以任意选择，但信号越长计算量越大

p = 12; % 尝试改变模型阶数，观察效果 Rxx = xcorr(x,'biased'); Rtemp = zeros(1,p); Rl = zeros(p,1); for k = 1:length(Rtemp)

Rtemp(k) = Rxx(length(x)-1+k); Rl(k) = Rxx(length(x)+k); end

Rs = toeplitz(Rtemp); % 生成自相关系数矩阵(Toeplitz型) A = -inv(Rs)*Rl; % AR模型系数估计

Sw = [Rtemp(1),Rl']*[1;A]; % 白噪声方差估计

% 采用malab自带函数估计模型系数

[a,E] = aryule(x,p); % a--系数，E--预测误差，k--反射系数 da = a(2:end)-A' % 自编程序求解是否正确？

w = randn(size(x));

x2 = filter(1,a,w); % 仿真数据 figure; subplot(3,1,1); plot(x); title('真实数据'); subplot(3,1,2); plot(x2); title('仿真数据');

本文来源：https://www.bwwdw.com/article/fi5g.html