生物医学信号心电信号预处理

学生姓名： 学号： 实验室名称： 项目名称：心电信号的预处理 项目内容： 1) 阅读相关文献资料，理解噪声的相关知识，重点是噪声有哪些分类，据此综述抑制信号中干扰噪声信号的基本理论和基本方法； 2) 设计两种分别用于抑制不同噪声的滤波器； 3) 基于W-H方程，设计最优滤波器； 4) 运用前面设计的几种滤波器，对加有噪声的模拟ECG信号进行去噪； 5) 总结不同滤波器的去噪效果。 原理（写出具体的计算公式） 说明： 模拟时加入的噪声：高斯白噪声 信噪比求解： 1.低通滤波器（巴特沃斯） 特点：通频带：频率响应曲线最大限度平坦，无起伏，理想的单位响应； 阻频带：逐渐下降为零，下降速度正比于滤波器的阶数，截断频率处有3dB衰减（ ）。 幅频特性：随滤波器阶次N的增加变得越来越好，在通频带内外有平稳的幅频特性，但过渡带较长，易造成失真。 对理想低通滤波的逼近：以原点附近的最大平坦响应来逼近理想低通滤波器，即巴特沃斯函数来近似滤波器的系统函数。 设计思路，如下图1-1所示： 图1-1 Step1:Fs:根据MIT-BIH数据库的心电信号，采样频率为360Hz; fp、fs:心电信号频率范围为0.05~100Hz,因此对于设计的巴特沃斯低通滤波器的通带截止频率和阻带截止频率定为100Hz、105Hz,滤去频率大于100Hz的心电信号； Rp、Rs:对于低通滤波器：通带波纹小于1dB,阻带衰减大于40dB（ ）； Step2: Step3：滤波器阶次: 幅频响应为： Step4:利用双线性变换法将模拟滤波器转换为数字滤波器； 2.带阻滤波器(巴特沃斯) 截止频率包括通带上下限截止频率fp1、fp2,下阻带截止频率fs1、上阻带截止频率fs2。 两边z变换，得数字滤波器H(z): 再根据零极点反向设计带馅数字滤波器，在频率w处出现凹陷（即滤波），把极点设置在零的径向上距圆点的距离为1-u处，得陷波器H(z): 其中：u越小，极点越靠近单位圆，滤掉的波形越接近w,根据采集到的MIT-BIH数据库，工频干扰的是60Hz,令1-u=0.9,经计算得： 从而确定fp1、fp2、fs1、fs2分别为：58.5，61.5,59.5,60.5. 3.小波去噪（小波阈值消噪） 特点：时间-频率分析方法，时间窗、频率窗可变，因此本质为函数逼近问题，即寻找含噪声信号对原信号的最佳逼近。 如下图1-2所示，实际是特征提取和低通滤波功能的综合，按照一定的阈值压缩信号的小波变换系数，然后用被压缩后的系数重构信号达到消噪目的，小波去噪效果好，可消除心电信号的基线漂移、工频干扰和肌电干扰。 图1-2 小波阈值去噪原理： 图1-3 Step1:对含噪声信号Signal_N做小波变换，将信号分解到不同程度频率的线性变换，即分层处理，并得到一组小波系数。 分层处理：进行小波分解选取的层数过高会导致边缘信息丢失，过度滤波，通过规定不同层数计算比较其信噪比大小，我们设定分解层数为3层达到最佳滤波效果； 小波函数：不同的小波函数导致分层结果不同，滤波效果相差迥异，db5具有很好的正规性，在本滤波器中，我们选取db5作为心电信号的小波进行处理，其效果最佳。 Step2:阈值处理即通过对第一步得到的小波系数用阈值准则计算处理，得到各层的阈值。 软阈值确定： 即把信号绝对值与指定的阈值进行比较，大于阈值的点变位该点与阈值的差，否则变0. 阈值准则：本滤波器设计时采用Sureshrink阈值。即基于Stein的无偏似然估计原理的自适应阈值选择，给定一个阈值 ，得到他的似然估，再将非似然 最小化，就得到所选的阈值，具体计算： 1).根据分层求取信号长度N; 2).将第i层的小波系数的平方由小到大排列，得到新向量 ； 3).计算风险向量 ，其中 4).以R中最小元素 为风险值，由 的相应位置B求出相应的 ，则阈值： Step3:利用以上两步得到的阈值和小波分解结构重建信号，达到阈值去噪效果。 4.最佳滤波器（基于W-H方程的维纳滤波器） 特点：利用信号和干扰的统计特征（自相关函数、功率谱等），从噪声中提取信号的滤波，以估计的结果与信号真值之间的均方误差最小作为最佳准则。 设计维纳滤波器实际上即选择h(n)，使其输出信号y(n)与期望信号d(n)误差的均方值为最小。 维纳滤波器原理： 以均方误差最小（LMS）为准则的，根据过去信号和当前信号来估计信号的当前值，其解形式是系统的传递函数或单位脉冲响应，其实质就是解维纳-霍夫（Wiener-Hopf）方程。 图1-4 根据时域求解的方法，有W-H方程： 简化形式为： ，对上式求逆，即 上式中，H为滤波器的系数向量， 为含有噪声的混合信号的自相关矩阵， 为混合信号和原始信号的互相关向量。 因此，先求出原始信号与含噪声信号的互相关函数及含噪声的自相关函数时（这两个函数可由样本得到），再通过矩阵求逆运算，得到维纳滤波器的最佳解。 当 为 的白噪声时， 可得： 即： 其中 表示单位阶跃函数。 相应的 变换为 。 预白化的系统函数为： 最优滤波器 最小均方误差为： 编写的源程序： 信噪比求解： function snr=SNR(I,In) % I: original signal % In: noisy signal(ie. original signal + noise signal) % snr=10*log(sigma2(I2)/sigma2(I2-I1)) snr=0; Ps=sum(abs(I).^2);%signal power Pn=sum(abs(In-I).^2);%noise power snr=10*log10(Ps/Pn); 一．巴特沃斯低通滤波器 1.巴特沃斯低通滤波器设计（去肌电干扰和高于100Hz频率的信号）： clear,clc fp = 100; fs = 105; %通带截止频率 阻带截止频率 Fs=360; Rp=1; Rs=40; %通带波纹 阻带衰减 wp=fp/(Fs/2); ws=fs/(Fs/2); %计算归一化角频率 wp =2*pi*tan(pi*wp/2);ws = 2*pi*tan(pi*ws/2); %数字转模拟 [n, Wn]=buttord(wp,ws,Rp,Rs,'s'); %计算阶数和截止频率 [B,A]=butter(N,wn,'s'); %计算模拟低通H(s) [b, a] = bilinear(B, A, 0.5); %双线性变换模拟转数字 [h, w] = freqz(b, a, 256, Fs); figure(1) pha=angle(h); subplot(2,1,1); plot(w,20*log10(pha)); grid on; title('巴特沃斯低通滤波器相频特性'); xlabel('频率（Hz）');ylabel('相位（dB）'); subplot(2,1,2); plot(w, abs(h)); title('巴特沃斯低通滤波器幅频特性'); xlabel(t');ylabel('幅值（dB）'); grid on; 2.巴特沃斯带阻滤波器（去60Hz工频干扰） clear,clc; Fs=360;T=1/Fs; Rp=1;Rs=40; fp1=58.5;fp2=61.5;fs1=59.5;fs2=60.5; wp1=fp1/(Fs/2); wp2=fp2/(Fs/2); %归一化 ws1=fs1/(Fs/2); ws2=fs2/(Fs/2); wc1 =(2/T)*tan(pi*wp1/2); wc2 =(2/T)*tan(pi*wp2/2); %频率预畸变 wr1 =(2/T)*tan(pi*ws1/2); wr2 =(2/T)*tan(pi*ws2/2); w0=sqrt(wc1*wc2); B=wc2-wc1; wp=1; %归一化通带截止频率 ws=wp*(wr1*B)/(w0^2-wr1^2); %归一化阻带截止频率 [N,wc]=buttord(wp,ws,Rp,Rs,'s'); %求滤波器阶数和3dB截止频率 [Z,P,K]=buttap(N); [Md,Nd]=zp2tf(Z,P,K); %将零极点形式转换为传输形式 [M,N]=lp2bs(Md,Nd,w0,B);%对低通滤波器进行频率变换，转换为带阻滤波器 [h,w]=freqs(M,N); %模拟带阻滤波器的幅频响应 figure(1); plot(w/(2*pi),abs(h)); xlabel('频率/Hz');ylabel('幅度');title('模拟带阻滤波器'); grid on; figure(2); [b,a]=bilinear(M,N,360); %对模拟滤波器双线性变换 [H,W]=freqz(b,a); %绘出频率响应 plot(W*Fs/(2*pi),abs(H)); title('数字滤波器幅频响应|H(ejOmega)| '); xlabel('频率/Hz');ylabel('幅值'); grid on; 3.小波去噪（小波阈值消噪） clear;clc %读数据 load('118e00m.mat'); %噪声信号数据 Signal=val(1,:); Fs=2*max(Signal); %采样频率确定 Signal=Signal/(Fs/2); %归一化 Signal_N=awgn(Signal,8); % 加入高斯噪声，信噪比为8dB. snr1=SNR(Signal,Signal_N) %加噪时SNR n=size(Signal_N); %小波去噪 [C L]=wavedec(Signal_N,3,'db5'); %分解层数3 多层次分解 （db5、db1、symlet8三种小波） cA3=appcoef(C,L,'db5',3); %提取尺度3下的近似小波系数 [cD1,cD2,cD3]=detcoef(C,L,[1,2,3]); %对信号做层数3的多尺度分解 cD1,cD2,cD3为各层小波系数 %使用stein的无偏似然估计原理进行选择各层的阈值 rigrsure为无偏似然估计阈值类型 thr1=thselect(cD1,'rigrsure'); thr2=thselect(cD2,'rigrsure'); thr3=thselect(cD3,'rigrsure'); TR=[thr1,thr2,thr3]; %各层的阈值 SORH='s'; %'s'为软阈值;'h'硬阈值。 [Signal_N_p,CXC,LXC,PERF0,PERF2]=wdencmp('lvd',Signal_N,'db5',3,TR,SORH); %去噪 小波分解结构、恢复和压缩的范数% lvd允许设置各层的阈 snr2=SNR(Signal,Signal_N_p) MSE=M/N hold on figure(1) subplot(3,1,1); plot(Signal); axis([0 2000 -3 2]) ; %xlim([0, 2000]); title('原ECG信号');ylabel('幅值（dB）'); subplot(3,1,2); plot(Signal_N); axis([0 2000 -3 2]) ; %xlim([0, 2000]); title('加入高斯白噪声的ECG信号'); ylabel('幅值（dB）'); subplot(3,1,3); plot(Signal_N_p); axis([0 2000 -3 2]) ; %xlim([0, 2000]); title('小波去噪后的ECG信号');ylabel('幅值（dB）'); 我们以118e00的数据为例，对比小波去噪前后信号的幅频特性图和信噪比进行说明： 小波去噪前SNR1=9.7706,小波去噪后SNR2=17.2769，信噪比增大，去噪效果增强。 4.最佳滤波器（基于W-H方程的维纳滤波器） M文件（W-H方程）： function [H,Emin]=WH(Rs,Rw,M) L1=(length(Rs)+1)/2; Rss=zeros(1,L1); %生成1*L1矩阵 for k=1:L1 Rss(k)=Rs(k+L1-1); %原信号的自相关 end L2=(length(Rw)+1)/2; Rww=zeros(1,L2); for k=1:L2 Rww(k)=Rw(k+L2-1); %白噪声的自相关 end Rx=zeros(1,M); for k=1:M Rx(k)=Rss(k)+Rww(k); end Rxx=zeros(M,M); for i=1:M for j=1:M if(i<=j) Rxx(i,j)=Rx(j-i+1); else Rxx(i,j)=Rx(i-j+1); end end end H=pinv(Rxx)*Rss'; %求矩阵的伪逆，转置 Emin=Rss(1)-sum(H*Rss); 主程序： clear,clc; n = input('信号的长度 n = '); %读数据 load('118e00m.mat'); %噪声信号数据 Signal=val(1,:); Fs=2*max(Signal); %采样频率确定 Signal=Signal/(Fs/2); %归一化 Signal_N=awgn(Signal,8); % 加入高斯噪声，信噪比为8dB. snr1=SNR(Signal,Signal_N) %加噪时SNR %维纳滤波 Mlag=100; %相关函数长度变量 Rxn=xcorr(Signal_N,Mlag,'biased'); %计算含噪声信号的自相关函数，有偏的互相关函数估计； Rxnx=xcorr(Signal_N,Signal,Mlag,'biased'); %产生含噪声信号与原始信号的互相关函数 %rxnx矩阵 rxnx=zeros(n,1); %产生n*1阶矩阵 rxnx(:)=Rxnx(101:101+n-1); %Rxx矩阵 Rxx=zeros(n,n); %产生含噪声信号的自相关矩阵 Rxx=diag(Rxn(101)*ones(1,n));%矩阵对角元素提取和创建对角阵 for i=2:n c=Rxn(101+i)*ones(1,n+1-i); Rxx=Rxx+diag(c,i-1)+diag(c,-i+1); end Rxx; h=zeros(n,1); h=pinv(Rxx)*rxnx; %计算维纳滤波器的h(n) Signal_N_p=filter(h,1,Signal_N); %将含噪声信号通过维纳滤波器滤波 snr2=SNR(Signal,Signal_N_p) figure(1) subplot(2,2,1) plot((-Mlag:Mlag),Rxn) %绘制自相关函数图像 title('118e00加噪后自相关函数图像') [f,p]=ksdensity(Signal_N); %计算输入信号的概率密度，f为样本点xi处的概率密度 subplot(2,2,2) plot(p,f) %绘制概率密度图像 title('118e00加噪后概率密度图像') X=fft(Signal_N); %离散傅里叶变换 Px=X.*conj(X)/600; %计算信号频谱 subplot(2,2,3) semilogy(Px) %绘制在半对数坐标系下频谱图像 title('加噪后频谱图像') %xlim([0, 2000]); axis(**************]);xlabel('x轴单位：w/rad','color','b') ylabel('y轴单位：w/HZ','color','b') pxx=periodogram(Signal_N); %计算输入信号的功率谱密度 subplot(2,2,4) semilogy(pxx) %绘制在半对数坐标系下功率谱密度图像 title('加噪后功率谱密度图像') axis(**************]);xlabel('x轴单位：w/rad','color','b') ylabel('y轴单位：w/HZ','color','b') figure(2) subplot(2,1,1) plot(Signal_N) %绘制输入信号图像 axis([0 2000 -2 1.5]); title('加噪后ECG信号图像'); xlabel('t');ylabel('f/HZ'); Rxn=xcorr(Signal_N,Mlag,'biased');%计算输入信号自相关函数 subplot(2,1,2) plot(Signal_N_p); axis([0 4000 -1 0]); %xlim([0, 2000]); title('维纳滤波器去噪后的ECG信号'); xlabel('频率（Hz）');ylabel('幅值（dB）'); 对含噪声信号的分析： 小波去噪前SNR1=9.8316,小波去噪后SNR2=15.7252，信噪比明显增大，去噪效果增强。 结论（画出要求的图形） 1.分别运用巴特沃斯低通滤波器、巴特沃斯带阻滤波器、小波去噪和最优滤波器对加有两种不同噪声的信号（无噪声心律失常ECG数据心电信号118e00.dat）进行去噪，并运用信噪比指标定量分析滤波器的去噪能力； 1).巴特沃斯低通滤波器结果（代码如文件夹中呈现的ECG-butt118.m、ECG-butt118e00.m）： 对心电信号118处理(加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =4.8706；（消噪前） snr2 =5.3493；（消噪后） 对心电信号118e00处理(加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =9.7975；（消噪前） snr2 =10.2573；（消噪后） 消噪后信噪比稍有增加，说明该巴特沃斯低通滤波器有一定的滤波作用，由于该滤波器仅对大于100Hz的心电信号进行去噪，不能完全滤除所有噪声，所以滤波效果不明显。 2).巴特沃斯带阻滤波器结果（代码如文件夹中呈现的ECG-butt（daizu118).m、ECG-butt(daizu118e00).m）： 对心电信号118处理(加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =4.8153；（消噪前） snr2 =4.8867；（消噪后） 对心电信号118e00处理(加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =9.9232；（消噪前） snr2 =9.9778；（消噪后） 消噪后信噪比稍有增加，说明该带阻滤波器有一定的滤波作用，由于该滤波器仅对60Hz工频干扰进行消噪，不能完全滤除所有噪声，所以滤波效果不明显。 3).小波去噪结果（代码如文件夹中呈现的ECG-wnoise110.m、ECG-wnoise118e00.m）： 对心电信号118e00的验证在上一部分已作说明处理，因此在这里呈现心电信号118的处理结果(同样加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =4.7810；（消噪前） snr2 =13.9620；（消噪后） 对比发现，消噪后信噪比增加较大，说明小波去噪相比于巴特沃斯低通滤波器和带阻滤波器，效果很好。 4).最佳滤波器去噪结果（代码如文件夹中呈现的Wiener118.m、Wiener118e00.m）： 对心电信号118e00的验证在上一部分已作说明处理，因此在这里呈现心电信号118的处理结果(同样加入高斯白噪声后消噪)： 消噪前后信噪比变化： snr1 =4.8350；（消噪前） snr2 =19.1876；（消噪后） 对比发现，采用维纳滤波器对含噪声信号进行滤波处理后，信噪比增加较大，说明维纳滤波器的消噪能力强。 总结 在此次心电信号预处理中，对于高于100Hz信号，我分别采用了巴特沃斯和切比雪夫滤波器进行滤波比较，对于60Hz工频干扰，我采用了巴特沃斯低通滤波器进行滤波处理，另外，我还加入了小波去噪对含噪声的心电信号进行处理，最后，结合W-H方程，给出了维纳滤波器。 以下是5种滤波器的信噪比呈现结果的汇总： 118 去噪后 5.3493 去噪前 9.7975 118e00 去噪后 10.2573 9.9778 10.1435 17.2769 15.7252 4.8867 7.3273 9.9232 9.7608 13.962 9.7706 19.1876 9.8316 巴特沃斯 低通 带通 切比雪夫 小波去噪 低通 4.781 4.835 4.8153 4.7776 维纳滤波器 去噪前 4.8706 综上所述，对于心电信号的滤波处理，维纳滤波器的去噪效果是最佳的，小波去噪的效果同样优秀，它们不仅可以去掉工频干扰，还可以去掉基线漂移和肌电信号等，另外，还对QRS波进行“细节”修正。但是小波去噪与维纳滤波器相比存在的不4.足之处在于易造成过度去噪，因此在小波变换部分的分层中，为增强其滤波能力，需要依次设定为2~7层左右，比较每一分解层数下滤波后的信噪比，从而确定最优分别层数，增加了操作的复杂性。对于另外三个滤波器，因为都只是对其中一部分噪声处理，其滤波效果不强，不宜采纳。 实习报告分数： 指导教师：