生物医学信号处理实习报告模版.docx
《生物医学信号处理实习报告模版.docx》由会员分享，可在线阅读，更多相关《生物医学信号处理实习报告模版.docx（24页珍藏版）》请在冰豆网上搜索。
生物医学信号处理实习报告模版
《生物医学信号处理》实习报告一
姓名：
学号：
实验室名称：
实习名称：
心电信号的预处理
内容：
1）根据相关文献资料，综述信号消噪的基本理论和基本方法；
2）设计两种分别用于抑制不同噪声的滤波器；
3）基于W-H方程，设计最优滤波器；
4）运用前面设计的几种滤波器，对加有噪声的模拟ECG信号进行去噪；
5）总结不同滤波器的去噪效果；
滤波器设计原理
一：
两种能消噪的滤波器的设计原理
1.巴特沃斯滤波器的设计原理
其特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零（对理想低通滤波的逼近：
巴特沃思滤波器是以原点附近的最大平坦响应来逼近理想低通滤波器）。
如下图一表示，可以看出滤波器的幅频特性是随着滤波器的阶次N的增加而变得越来越好，在截止频率
有：
（1）
衰减具有
不变性。
通带、阻带均具有单调下降的特性。
图一
低通巴特沃斯滤波器的幅频响应为：
（2）
其中，
为滤波器的阶数，取正整数；
为滤波器的截止频率，为其
的频率；
为滤波器的通带截止频率；
2.切比雪夫滤波器的设计原理
切比雪夫滤波器有两类，一类是在同带内幅度频率响应呈现等波纹特性，而阻带内是单调的全极点滤波器；另一类是在通带内幅度频率响应呈现单调特性，而阻带内是等波纹特性，同时有零点和极点的滤波器，这类滤波器的零点位于s平面的虚轴上。
如图二，为第Ⅰ类切比雪夫滤波器的幅频特性曲线。
图二
第Ⅰ类切比雪夫滤波器的幅度平方函数定义为
（3）
其中，
为切比雪夫多项式（
），
为滤波器阶数；
为通带截止角频率，此处是指被通带波纹所限制的最高角频率，；
为小于1的正数，表示通带内幅度波动的程度，
越小，幅度波动越小；
其特点为：
1）当
时，
在
之间呈等波纹变化，
越小，波动幅度越小；
2）无论
为何值，所有的曲线都通过
点，该点被定义为
点；
3）当
时，若
为奇数，
；若
为偶数，
（4）；
4）当
时，曲线呈单调下降；
5）通带内的起伏使对应的相频特性呈现非线性；
3.维纳滤波器的设计原理
维纳滤波器属于现代滤波器，传统的滤波器只能滤除信号和干扰频带没有重叠的情况，当信号和干扰频带有重叠的时候传统滤波器将无能为力，这时就需要用到现代滤波器，现代滤波器利用信号和干扰的统计特征（如自相关函数、功率谱等）导出一套最佳估计算法，然后用硬件或软件予以实现。
下面是对维纳滤波器的设计仿真：
维纳滤波器是以均方误差最小（LMS）为准则的，它根据过去观测值和当前观测值来估计信号的当前值，因此它的解形式是系统的传递函数或单位脉冲响应
图三
设计维纳滤波器的过程就是寻求在最小均方误差下滤波器的单位脉冲响应或传递函数的表达式，其实质就是解维纳-霍夫（Wiener-Hopf）方程。
均方误差为：
（5）
维纳—霍夫方程最小均方误差是下的解为：
（6）
其中，
为滤波器的系数向量，
为含有噪声的混合信号的自相关矩阵，
为混合信号和原始信号的互相关向量，为此先求出混合信号的自相关函数以及混合信号和原始信号的互相关函数，这两个函数，我们需要有样本得到。
滤波器实现的代码
1.巴特沃斯滤波器（Butter.m）
[B,A]=butter（3,2*pi*1000,'s'）;%脉冲响应不变法低通滤波器
[num1,den1]=impinvar（B,A,4000）;
[h1,w]=freqz（num1,den1）;
[B,A]=butter（3,2/0.00025,'s'）;%双线性变换法低通滤波器
[num2,den2]=bilinear（B,A,4000）;
[h2,w]=freqz（num2,den2）;
f=w/pi*2000;
plot（f,abs（h1）,'GREEN',f,abs（h2）,'RED'）;
legend（'脉冲响应不变法','双线性变换法'）;
grid;
xlabel（'频率/Hz'）
ylabel（'幅值/dB'）
title（'巴特沃斯低通滤波器'）;
图四
2.切比雪夫滤波器（qiebixuefu1.m）
clc,clear
%确定数字滤波器指标
Fs=10000;%采样频率
fp=280;fs=450;
wp=2*pi*fp/Fs;%通带截止频率
Ap=0.1;%通带内的最大衰减
ws=2*pi*fs/Fs;%阻带截止频率
As=40;%阻带内的最小衰减
%数字角频率转换成模拟角频率
%双线性变换法
wps=2*Fs*tan（wp/2）;
wss=2*Fs*tan（ws/2）;
%设计模拟滤波器
k1=sqrt（（10^（0.1*As）-1）/（10^（0.1*Ap）-1））;
qsp=wss/wps;
N=ceil（acosh（k1）/acosh（qsp））;
e=sqrt（10^（0.1*Ap）-1）;
Ap=10*log（1+e^2）;
%计算得到模拟低通H（s）分子和分母的系数
[z,p,k]=cheb1ap（N,Ap）;
[bp,ap]=zp2tf（z,p,k）;
[bs,as]=lp2lp（bp,ap,Fs）;
%复变量映射s-->z
[b,a]=bilinear（bs,as,Fs）;
%画出切比雪夫1型低通滤波器的幅频响应曲线图
[H,w]=freqs（bp,ap,2000）;%按n指定的频率点给出频率响应
magH2=（abs（H））.^2;%给出传递函数的幅度平方
%输出波形
plot（w,magH2）;title（'切比雪夫1型低通滤波器的幅频响应特性'）;
xlabel（'w/wc'）;
ylabel（'Chebyshev1|H（jw）|^2'）;
图五
3.维纳滤波器（Winner.m）
function[H,Emin]=Winner（Rs,Rw,N）
L1=（length（Rs）+1）/2;
Rss=zeros（1,L1）;
fork=1:
L1
Rss（k）=Rs（k+L1-1）;
end
L2=（length（Rw）+1）/2;
Rww=zeros（1,L2）;
fork=1:
L2
Rww（k）=Rw（k+L2-1）;
end
Rx=zeros（1,N）;
fork=1:
N
Rx（k）=Rss（k）+Rww（k）;
end
Rxx=zeros（N,N）;
fori=1:
N
forj=1:
N
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）;
测试程序：
（WinnerTest.m）
clc;
t=-50:
0.1:
50;
w=0.1*pi;
s=cos（w.*t）;
w=randn（1,length（t））;
x=s+w;
Rs=xcorr（s）;%信号的自相关函数
Rw=xcorr（w）;%白噪声的自相关函数
[h,E]=Winner（Rs,Rw,length（t））;%维纳滤波
se=filter（h,1,x）;
plot（t,s）;
holdon
plot（t,x,'g'）;
holdon
plot（t,se,'r'）;
e=s-se;
plot（t,e,'BLACK'）;
legend（'原始信号','加噪后信号','维纳滤波后信号','误差'）;
图六
4.分别运用上面设计的滤波器和最优滤波器分别对含噪声的模拟ECG信号进行去噪，并运用SNR指标定量分析滤波器的去噪能力：
1）读入原始ECG模拟信号（代码见ecgsyn.m和derivsecgsyn.m），并存入ECG.txt中；
2）提取ECG.txt中的（500-2500个信号），保存至ECG（2000）.txt中，并将生成的心电信号中加入高斯噪声，使其生成含噪的心电信号，保存至ECGWithGauss（2000）.txt中（代码见ECG（2000）.m），如图七；
图七
3）设计SNR（信噪比）的计算方法
functionsnr=SNR（I,In）
%计算信噪比函数
%I:
originalsignal
%In:
noisysignal（ie.originalsignal+noisesignal）
snr=0;
Ps=sum（sum（（I-mean（mean（I）））.^2））;%signalpower
Pn=sum（sum（（I-In）.^2））;%noisepower
snr=10*log10（Ps/Pn）
%其中I是纯信号，In是带噪声的信号，snr是信噪比
4）利用所设计的三种滤波器对所生成的ECG信号除噪；
a.巴特沃斯滤波器消噪（SNRButter.m）
clear,clc
X=importdata（'ECG（2000）.txt'）;%加噪心电数据读入
X1=X（:
）;
X2=importdata（'ECGWithGauss（2000）.txt'）;
X3=X2（:
）;
[B,A]=butter（3,2*pi*1000,'s'）;
[num1,den1]=impinvar（B,A,4000）;
[h1,w]=freqz（num1,den1）;
f=w/pi*2000;
y=filter（num1,den1,X3）;
figure;
subplot（3,1,1）;
plot（X1）;%原ECG信号
title（'原始信号'）;
subplot（3,1,2）;
plot（X3,'g'）
title（'加入高斯噪声后的ECG信号'）;%加入高斯噪声之后的噪声
subplot（3,1,3）;
plot（y,'r'）;
title（'除噪后的ECG信号'）;
snr1=SNR（X1,X3）;
sne2=SNR（X1,y）;
b.切比雪夫1型滤波器（SNRQie.m）
clc,clear
%确定数字滤波器指标
Fs=1000;%采样频率
fp=5;%通带带截止频率
fs=200;%阻带截止频率
wp=2*pi*fp/Fs;%通带截止频率
Ap=0.1;%通带内的最大衰减
ws=2*pi*fs/Fs;%阻带截止频率
As=40;%阻带内的最小衰减
%数字角频率转换成模拟角频率
%双线性变换法
wps=2*Fs*tan（wp/2）;
wss=2*Fs*tan（ws/2）;
%设计模拟滤波器
k1=sqrt（（10^（0.1*As）-1）/（10^（0.1*Ap）-1））;
qsp=wss/wps;
N=ceil（acosh（k1）/acosh（qsp））;
e=sqrt（10^（0.1*Ap）-1）;
Ap=10*log（1+e^2）;
%计算得到模拟低通H（s）分子和分母的系数
[z,p,k]=cheb1ap（N,Ap）;
[bp,ap]=zp2tf（z,p,k）;
[bs,as]=lp2lp（bp,ap,Fs）;
%复变量映射s-->z
[b,a]=bilinear（bs,as,Fs）;
X=importdata（'ECG（2000）.txt'）;%心电信号数据导入
X1=X（:
）;
X2=importdata（'ECGWithGauss（2000）.txt'）;%含噪心电信号的读入
X3=X2（:
）;
%调用滤波器实现函数filter进行滤波
y=filter（b,a,X3）;
%画图比较
subplot（3,1,1）;
plot（X1）;
title（'原始心电信号'）;
subplot（3,1,2）;
plot（X3,'g'）%加入高斯噪声后的信号
title（'含噪心电信号'）;
subplot（3,1,3）;
plot（y,'r'）;
title（'经切比雪夫滤波后的心电信号'）;
%调用SNR
snr1=SNR（X1,X3）;
snr2=SNR（X1,y）;
c.维纳滤波器（SNRWinner.m）
clear,clc
Date=load（'ECG（2000）.txt'）;
mk=zeros（1,2000）;
mk（1,:
）=Date（1:
2000）;
mk=mk（1,:
）;
b=max（mk）;
a=min（mk）;
n=0:
1999;
mk=（mk-a）/（b-a）;
Date1=load（'ECGWithGauss（2000）.txt'）;
k=zeros（1,2000）;
k（1,:
）=Date1（1:
2000）;
k=k（1,:
）;
b1=max（k）;
a1=min（k）;
n1=0:
1999;
k=（k-a1）/（b1-a1）;
Rs=xcorr（mk）;%自相关
Rw=xcorr（k-mk）;
[h,E]=Winner（Rs,Rw,length（n1））;
y=filter（h,1,k）;%维纳滤波器
axis（[0200001]）;
subplot（3,1,1）;
plot（n,mk）;
title（'原始信号'）;
subplot（3,1,2）;
plot（n,k,'g'）%加入高斯噪声后的信号
%plot（n,x）;
axis（[0200001]）;
title（'含噪信号'）;
subplot（3,1,3）;
plot（n,y,'r'）;
axis（[0200001]）;
title（'经维纳滤波后的信号'）;
snr1=SNR（mk,k）%滤波前信噪比
snr2=SNR（mk,y）%滤波后信噪比
5.分别运用上面设计的滤波器和最优滤波器对含高斯噪声心律失常ECG数据消噪，并运用SNR指标定量分析滤波器的去噪能力
1）数据的获取：
运行PhysioNet_ECG_Exporter_modified.m文件，读入118.dat,118.hea,118.atr，获取120秒的无噪声心律失常信号数据，保存在118_ECG_0_120.mat文件中。
实验中用的是它的第一张表格数据；
2）利用三种滤波器滤波：
clear,clc;
realecg=load（'118_ECG_0_120.mat'）;
mk=zeros（5,1000）;
mk（1,:
）=realecg.ECG_1（1:
1000）;%获取第一张表格数据
mk=mk（1,:
）;
N=length（mk）;
b=max（mk）;
a=min（mk）;
fs=N/（b-a）;
n=0:
999;
mk=（mk-a）/（b-a）;
x=awgn（mk,8,'measured'）;%加入高斯噪声，信噪比为8dB.
Rs=xcorr（mk）;
Rw=xcorr（x-mk）;
[h,E]=Winner（Rs,Rw,length（n））;
se1=filter（h,1,x）;
e=mk-se1;
holdon
%切比雪夫滤波
%确定数字滤波器指标
Fs=1000;%采样频率
fp=5;%通带带截止频率
fs=200;%阻带截止频率
wp=2*pi*fp/Fs;%通带截止频率
Ap=0.1;%通带内的最大衰减
ws=2*pi*fs/Fs;%阻带截止频率
As=40;%阻带内的最小衰减
%数字角频率转换成模拟角频率
%双线性变换法
wps=2*Fs*tan（wp/2）;
wss=2*Fs*tan（ws/2）;
%设计模拟滤波器
k1=sqrt（（10^（0.1*As）-1）/（10^（0.1*Ap）-1））;
qsp=wss/wps;
N=ceil（acosh（k1）/acosh（qsp））;
e=sqrt（10^（0.1*Ap）-1）;
Ap=10*log（1+e^2）;
%计算得到模拟低通H（s）分子和分母的系数
[z,p,k]=cheb1ap（N,Ap）;
[bp,ap]=zp2tf（z,p,k）;
[bs,as]=lp2lp（bp,ap,Fs）;
[b,a]=bilinear（bs,as,Fs）;%实现S域到Z域的变换
%调用滤波器实现函数filter进行滤波
se2=filter（b,a,x）;
e3=mk-se2;
e0=mk-x;
holdon
[B,A]=butter（1,1/0.00025,'s'）;
[num1,den1]=bilinear（B,A,4000）;%双线性变换法
[h1,w]=freqz（num1,den1）;
se3=filter（num1,den1,x）;
subplot（511）;
plot（n,mk）;
title（'原序列'）;
subplot（512）;
plot（n,x）;
axis（[0100001]）;
title（'混入噪声的信号'）;
subplot（513）;
plot（n,se1）;
axis（[0100001]）;
title（'经维纳滤波后的信号'）;
subplot（514）;
plot（n,se3）;
axis（[050001]）;
e3=mk-se3;
title（'巴特沃斯低通滤波后的信号'）;
subplot（515）;
plot（n,se2）;
axis（[050001]）;
title（'切比雪夫低通滤波后的信号'）;
snr1=SNR（mk,x）
snr2=SNR（mk,se1）
snr3=SNR（mk,se3）
snr4=SNR（mk,se2）
6.分别运用上面设计的滤波器和最优滤波器对含噪声心律失常ECG数据消噪，并运用SNR指标定量分析滤波器的去噪能力：
1）数据的获取：
运行PhysioNet_ECG_Exporter_modified.m文件，读入118.dat,118.hea,118.atr，获取120秒的无噪声心律失常信号数据，保存在118_ECG_0_120.mat文件中。
以同样方法读入118e06，获取含噪声心律失常信号数据，在这里直接用118e06.dat文件读取心律失常含噪信号数据，用于滤波实验。
实验中用的均为第一张表格的数据。
取1000个采样点数据；
2）代码实现：
clear,clc;
realecg=load（'118_ECG_0_120.mat'）;
mk=zeros（5,1000）;
mk（1,:
）=realecg.ECG_1（1:
1000）;%获取第一张表格数据
mk=mk（1,:
）;
N=length（mk）;
b=max（mk）;
a=min（mk）;
fs=N/（b-a）;
n=0:
999;
mk=（mk-a）/（b-a）;
x=awgn（mk,8,'measured'）;%加入高斯噪声，信噪比为8dB.
Rs=xcorr（mk）;
Rw=xcorr（x-mk）;
[h,E]=Winner（Rs,Rw,length（n））;
se1=filter（h,1,x）;
e=mk-se1;
holdon
%切比雪夫滤波
%确定数字滤波器指标
Fs=1000;%采样频率
fp=5;%通带带截止频率
fs=200;%阻带截止频率
wp=2*pi*fp/Fs;%通带截止频率
Ap=0.1;%通带内的最大衰减
ws=2*pi*fs/Fs;%阻带截止频率
As=40;%阻带内的最小衰减
%数字角频率转换成模拟角频率
%双线性变换法
wps=2*Fs*tan（wp/2）;
wss=2*Fs*tan（ws/2）;
%设计模拟滤波器
k1=sqrt（（10^（0.1*As）-1）/（10^（0.1*Ap）-1））;
qsp=wss/wps;
N=ceil（acosh（k1）/acosh（qsp））;
e=sqrt（10^（0.1*Ap）-1）;
Ap=10*log（1+e^2）;
%计算得到模拟低通H（s）分子和分母的系数
[z,p,k]=cheb1ap（N,Ap）;
[bp,ap]=zp2tf（z,p,k）;
[bs,as]=lp2lp（bp,ap,Fs）;
[b,a]=bilinear（bs,as,Fs）;%实现S域到Z域的变换
%调用滤波器实现函数filter进行滤波
se2=filter（b,a,x）;
e3=mk-se2;
e0=mk-x;
holdon
[B,A]=butter（1,1/0.00025,'s'）;
[num1,den1]=bilinear（B,A,4000）;%双线性变换法
[h1,w]=freqz（num1,den1）;
se3=filter（num1,den1,x）;
subplot（511）;
plot（n,mk）;
title（'原序列'）;
subplot（512）;
plot（n,x）;
axis（[0100001]）;
title（'混入噪声的信号'）;
subplot（513）;
plot（n,se1）;
axis（[0100001]）;
title（'经维纳滤波后的信号'）;
subplot（514）;
plot（n,se3）;
axis（[050001]）;
e3=mk-se3;
title（'巴特沃斯低通滤波后的信号'）;
subplot（515）;
plot（n,se2）;
axis（[050001]）;
title（'切比雪夫低通滤波后的信号'）;
%计算其信噪比
snr1=SNR（mk,k）
snr2=SNR（mk,se1）
snr3=SNR（mk,se3）
snr4=SNR（mk,se2）
实验数据分析及结果
1.运用上面设计的滤波器和最优滤波器分别对含噪声的模拟ECG信号进行去噪，并运用SNR指标定量分析滤波器的去噪能力：
1）图八即为运用巴特沃斯滤波器对模拟ECG信号进行去噪的效果图：
其中：
snr1=9.0667（除噪前）
snr2=9.0937（除噪后）
图八
2）图九即为运用切比雪夫滤波器对模拟ECG信号进行去噪的效果图：
其中：
snr1=9.0667（除噪前）
snr2=10.8936（除噪后）
图九
3）图九即为运用维纳滤波器对模拟ECG信号进行去噪的效果图：
其中：
snr1=8.5580（除噪前）
snr2=10.6860（除噪后）
升级会员 