滤波器设计及在心电信号滤波中的应用.docx
《滤波器设计及在心电信号滤波中的应用.docx》由会员分享，可在线阅读，更多相关《滤波器设计及在心电信号滤波中的应用.docx（16页珍藏版）》请在冰豆网上搜索。
滤波器设计及在心电信号滤波中的应用
目录
一课题概述.......................................................1
1.课题名称........................................................1
2.题目及要求......................................................1
二设计的详细步骤.................................................1
1.心电信号采集....................................................1
2.心电信号分析....................................................1
3.含噪心电信号合成................................................1
4.数字滤波器设计及滤波............................................2
5.心电信号波形观察、频谱观察......................................2
三具体模块过程及结果（图形）.....................................3
1.设计的模块图....................................................3
2.模块语句及实现结果..............................................3
四设计总结.......................................................15
一课题概述
1.课题名称
数字滤波器设计及在心电信号滤波中的应用
2.题目及要求
采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪心电信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。
。
二设计的详细步骤
1.心电信号采集
心电信号作为心脏电活动在人体体表的表现，信号一般比较微弱，幅度在10μV~5mV，频率为0.05~100Hz。
在心电信号的采集、放大、检测及记录过程中，有来自外界的各种干扰。
记录一段时间内的人体心电信号波形，要求长度不小于10秒，并对记录的信号进行数字化，保存为数据文件；这里，我们使用美国的MIT/BIH心电原始数据，由实验老师给出一定长度的心电原始数据，数据保存在文件“a01.txt～a12.txt”中，在MATLAB中通过如下语句读取：
%从当前路径下的a01.txt文件读取心电原始数据到变量a01中，a01为二维数据，第一列为心电信号时间，第二列为心电信号幅度。
a=load（'a15.txt'）;
2.心电信号分析
使用MATLAB绘出数字化后的心电信号的时域波形和频谱图。
根据频谱图求出其带宽，并说明心电信号的基本特征。
3.含噪心电信号合成
在MATLAB软件平台下，给原始的心电信号叠加上噪声或干扰，干扰类型分为如下几种：
1）白噪声；2）工频干扰（50Hz）；3）谐波干扰（二次、三次谐波为主，分别为100Hz、150Hz）；4）其它干扰，可设置为低频、高频、带限噪声，或冲激干扰。
绘出叠加噪声后的心电信号时域和频谱图，在视觉上与原始心电信号图形对比，绘出其时域波形差，分析频域基本特征变化。
4.数字滤波器设计及滤波,完成以下题目中的一个
给定滤波器的规一化性能指标（参考指标,实际中依据每个同学所叠加噪声情况而定）例如：
通带截止频率wp=0.25*pi,阻通带截止频率ws=0.3*pi;通带最大衰减Rp=1dB;阻带最小衰减Rs=15dB，每个题目至少设计出3个用不同方法的不同类型滤波器。
题目
（1）：
采用窗函数法与等波纹法分别设计各型FIR滤波器（低通、高通、带通、带阻中的至少3种类型）来对叠加干扰前后的心电信号进行滤波处理，绘出滤波器的频域响应，绘出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比，分析信号的变化；在相同的性能指标下比较各方法的滤波效果，并从理论上进行分析（或解释）。
题目
（2）：
采用双线性变换法与脉冲响应不变法，分别利用不同的原型低通滤波器（Butterworth型与切比雪夫I型）来设计各型IIR滤波器（低通、高通、带通、带阻中的至少3种类型），绘出滤波器的频域响应；并用这些数字滤波器对含噪心电信号分别进行滤波处理，比较不同方法下设计出来的数字滤波器的滤波效果，并从理论上进行分析（或解释）。
5.心电信号波形观察、频谱观察
对滤波后的心电信号观察其时域、频域特征变化。
绘出滤波后、滤波前、加噪后三个心电信号的差值波形，观察相互间的差异性；同时，分析频谱变化。
学生也可选用持续时间更长的心电原始数据，加上干扰后按上述要求设计滤波器。
三具体模块过程及结果（图形）
1.设计的模块图
2.模块语句及实现结果
（1）源数据的导入
a=load（'a15.txt'）;
这里我选择第十五组数据
（2）绘出源信号的时域波形图和频域波形图
%绘出时域波形图
t=0:
0.001:
0.001*（1000-1）;
figure
（1）;
subplot（2,1,1）;
plot（a（:
1）,a（:
2））;
title（'时域波形图'）;xlabel（'时间t/s'）;ylabel（'幅值/A'）;
%绘频域波形图
y1=fft（a（:
2）,1000）;
x1=a（:
2）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（y1））;%画频谱图
title（'频域波形图'）;xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
图一
（３）加入噪声干扰
这里我加入了三个噪声，分别为白噪声、单频正弦干扰、chirp噪声
1）白噪声
q=0.8*rand（1000,1）;%干扰信号为白噪声
s=a（:
2）+q;
[hs,ws]=freqz（s,1,1024）;
abshs=abs（hs）;
figure
（2）;
subplot（2,1,1）,plot（s（1:
1000））;
title（'加白噪音后的时域图'）;xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（2,1,2）;,plot（ws/pi,abshs）;
title（'加白噪音后的频谱图'）;
ylabel（'幅值'）;xlabel（'Hz'）;
图二
说明：
观察加白噪声后的时域和频域图，与未加噪声前的比较，会发现频谱图中在0Hz处的幅值增加很大,同时在没有谱线的频率上也出现了频谱，这是由于白噪声产生的，因为白噪声在所有的频率上都有频谱.
2）单频正弦干扰
Y=[sin（2*pi*80*t）]';%干扰信号为单频函数（正弦）
x2=x1+Y;
y2=fft（x2,1000）;
f2=100*（0:
999）/1000;%频率为100Hz
figure（3）;
subplot（2,1,1）;
plot（t,x2）;
title（'加单频正弦干扰后的时域波形图'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（2,1,2）;
plot（f2,abs（y2））;
title（'加单频正弦干扰后的信号频谱'）;
ylabel（'幅值'）;
xlabel（'Hz'）;
图三
说明：
观察图一和图三，加噪后的时域图在整个时间轴上的幅值大多集中在0附近，而且比较密集，对于频谱图则可以很明显的看出在大约8Hz和92Hz（由于对称性，后半个周期不用分析）附近幅值增长很大，这个是由80Hz的正弦噪声信号产生的，在其他的频率范围内，较不易观察。
3）chirp噪声
p=0.5*chirp（a（:
1）,0,a（1000,1）,200）;%干扰信号为chirp
k=a（:
2）+p;
[hc,wc]=freqz（k,1,1024）;
abshc=abs（hc）;
figure（4）;
subplot（2,1,1）;plot（k（1:
1000））;
title（'加chirp噪音后的时域图'）;xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（2,1,2）;plot（wc/pi,abshc）;
title（'加chirp噪音后的频谱图'）;ylabel（'幅值'）;xlabel（'Hz'）;
图四
说明：
加chirp噪音干扰，加噪前时域波形图中的幅值对各个时间点成均匀分布，加了噪音后，则前疏后密，就加噪前后的频谱图而言，与加白噪声后有些许类似，即在个别点处幅值增加很多（加噪前，信号频谱出现在0到100Hz范围内，加噪后，在0到200Hz都会出现谱线，因为所加的chirp噪声的频率是时变的，变化范围是0到200Hz，频谱幅度与所加信号的幅度有关），在其他点则不明显，原因也与横坐标取值有关。
（4）滤波器的设计（三种）
1）双线性变换法实现巴特沃斯低通滤波器
figure（5）;
fs=100;
f1=5;f2=10;%通带、阻带截止频率
wp=（f1/fs）*2*pi;%临界频率采用角频率表示
ws=（f2/fs）*2*pi;%临界频率采用角频率表示
Omegap=2*fs*tan（wp/2）;
Omegas=2*fs*tan（ws/2）;
[n1,Wn]=buttord（Omegap,Omegas,1,50,'s'）;%计算阶数和截止频率
[b,a]=butter（n1,Wn,'s'）;%设计低通巴特沃斯模拟滤波器
[bz,az]=bilinear（b,a,fs）;%将模拟滤波器系数变为近似等价的数字滤波器系数
freqz（bz,az,512,fs）;%滤波器频谱图
y=filter（bz,az,x2）;%利用低通巴特沃斯数字滤波器实现对x2的滤波
figure（6）;
subplot（2,1,1）;
plot（y）;
title（'滤单频正弦干扰波后的信号时域波形图'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
y3=fft（y,1000）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（y3））;
title（'滤白噪声后的信号频域波形图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
图五
说明：
此图为巴特沃斯低通滤波器频率响应图
图六
说明：
此图为用低通滤波器滤单频正弦干扰信号的时域和频域波形图，由于信号主要处于频段的低频部分，而单频正弦干扰信号的频谱在整个频段上呈对称分布，即与源信号的频段分布是类似的，对于与源信号共存的低频部分用选频滤波器无法滤除，因此用低通滤波器将噪声信号的高频部分滤掉。
观察图三和图六我们会发现对于加入了正弦干扰波的源信号，滤波后的时域波形图的幅值显然比滤波前的小，对频域波形图则更加明显的可以看出频率在0~10和90~100之间的幅值变化比较大，而对于10~90间的区域，滤波后的频域波形图中更趋于直线。
2）双线型变换法实现切比雪夫1数字滤波器
figure（7）;
fs=100;
f11=10;f12=25;f21=5;f22=30;%通带、阻带截止频率
Wp1=（f11/fs）*2*pi;%临界频率采用角频率表示
Ws1=（f21/fs）*2*pi;%临界频率采用角频率表示
Wp2=（f12/fs）*2*pi;%临界频率采用角频率表示
Ws2=（f22/fs）*2*pi;%临界频率采用角频率表示
Omegap1=2*fs*tan（Wp1/2）;Omegap2=2*fs*tan（Wp2/2）;
Omegas1=2*fs*tan（Ws1/2）;Omegas2=2*fs*tan（Ws2/2）;
BW=Omegap2-Omegap1;%频带宽度
W0=Omegap2*Omegap1;W00=sqrt（W0）;
WP=1;WS=WP*（W0^2-Ws1^2）/（Ws1*BW）;
[n1,WN]=buttord（WP,WS,1,50,'s'）;%计算阶数和截止频率
[B,A]=cheby1（n1,1,WN,'s'）;%设计低通切比雪夫1模拟滤波器
[BT,AT]=lp2bp（B,A,W00,BW）;
[num,den]=bilinear（BT,AT,0.5）;%将模拟滤波器系数变为近似等价的数字滤波器系数
freqz（num,den,64）;%低通切比雪夫1数字滤波器的频率响应
y=filter（num,den,s）;%用低通切比雪夫1数字滤波器实现对白噪声的滤波
figure（8）;
subplot（2,1,1）;plot（y）;
title（'滤白噪声后的信号时域波形图'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
s3=fft（y,1000）;
f1=100*（0:
999）/1000;
subplot（2,1,2）;
plot（f1,abs（s3））;
title（'滤白噪声后的信号频域波形图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
图七
说明:
此图为切比雪夫1型滤波器的频率响应图.
图八
说明：
此图为用低通滤波器滤白噪声后的时域和频域波形图，要选择低通是由于信号主要处于频段的低频部分，而白噪声的频谱在整个频段上都存在，对于与源信号共存的低频部分用选频滤波器无法滤除，因此用低通滤波器将噪声信号的高频部分滤掉。
观察图二和图七，二者在时间域波形图上的幅值大致相当，但在频域波形图中图二在0Hz处的幅值是很高的，但滤波后就明显变小了，在其他频率值范围内由于二图的横坐标的取值范围不同而不太容易比较。
3）带阻滤波器
figure（9）;
d=fir1（1000,[0.1570.17],'stop'）;%设计带阻滤波器
freqz（d,512）;%带阻滤波器的频率响应图
y4=filter（d,1,x2）;%滤波器实现滤波
figure（10）;
subplot（2,1,1）;plot（y4）;
title（'滤单频正弦干扰波后的信号时域波形图'）;
xlabel（'时间t/s'）;ylabel（'幅值/A'）;
subplot（2,1,2）;
y5=fft（y4,1000）;
f1=100*（0:
999）/1000;
plot（f1,abs（y5））;
title（'滤单频正弦干扰波后的信号频域波形图'）;
xlabel（'频率f/Hz'）;ylabel（'幅值/db'）;
图九
说明:
此图为带阻滤波器的频率响应图.
图十
说明：
选择带阻滤波器滤单频正弦干扰波后，效果最显著的地方在于8Hz和92Hz附近原来的很大的振幅被滤掉了。
之所以选择带阻滤波器是因为单频政宣干扰是对称的，且与原信号的高幅值部分是在相同的频率段上，加噪声后在高幅值部分的频率段上幅值就会更高，所以我们选择带阻滤波器以滤除噪声的干扰作用。
升级会员 