滤波器设计及在心电信号滤波中的应用

[hs,ws]=freqz(s,1,4096); %求离散系统频响特性的函数freqz()
abshs=abs(hs);%hs取绝对值赋给abshs
figure(2);%新建图像2
subplot(2,1,1);
plot(s(1:4000));%绘出从1取到4000时s的图形
title('加入高频干扰信号后的时域图');xlabel('时间（s）');ylabel('幅值（A）');
plot(t);%绘出t的图形
title('原始波形图');xlabel('时间（s）');ylabel('幅值（A）');
y1=fft(a(:,1),4000);%行取全部
f1=100*(0:3999)/4000;%先生成一个0,1,2,...,3999的整数向量，然后对对这个向量的每一项乘以100除以4000.
对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz；由于电子设备采集到的信号经常会混有电源线干扰，电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。对于电源线干扰，可以让信号通过一个带阻滤波器。
某一个实际心电图信号的采样序列存于文件ecg.txt中，共4000个采样点。
⑤在 处平方幅度响应的各级导数均存在且等于0，因此 在该点上取得最大值，且具有最大平坦特性。
图1展示了2阶、4阶、8阶巴特沃斯低通滤波器的幅频特性。可见阶数n越高，其幅频特性越好，低频检测信号保真度越高，过渡带变窄，即衰减加剧，但半功率点不变。
图1 巴特沃斯低通滤波器的幅频特性
2.2 切比雪夫I型数字低通滤波器
3.3.2电源线干扰（50Hz）
用x表示加入电源线干扰后的信号，再对电源线信号的4000个频率点进行FFT运算，在进行相关的运算后，通过调用plot函数直接绘出加入电源线干扰后的时域图和频谱图。具体步骤如下：
x2=sin(2*pi*50*t);%x2表示正弦信号t=0:0.00025:0.00025*(4000-1);%从0开始间隔0.00025取值到 0.00025*(4000-1)
(3)按照模拟低通滤波器的技术指标设计模拟低通滤波器。
(4)利用双线性变换法将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。
（5）数字低通技术指标为：
=0.4πrad， =1dB； =0.5πrad， =40dB
（6）模拟低通的技术指标为：
，T=1
，
,
归一化截止角频率wp=2pi*Fs/Ft; ws=2pi*Fs/Ft
subplot(2,1,2);
plot(f1,abs(s3));
title('滤除高频信号后的频谱图');xlabel('频率(Hz)');ylabel('幅值(dB)');
3.4.2高通滤波器的设计
高通滤波器的设计可通过模拟低通滤波器再经频率变换而实现。其中模拟低通滤波器可根据已经存在的典型滤波器，如巴特沃斯滤波器等逼近实现，而由低通到高通转换理论依据在绪论部分已经进行了详细的论证，又本设计基于MATLAB仿真软件实现，可利用MATLAB信号处理工具箱提供的各种函数模型实现，可忽略其中的一些复杂的函数变换，从而简化理论设计和论证，具体如下：
说明：
1、使用MATLAB或其它编程软件编程实现。
2、如果滤波效果不明显，则需变动滤波器参数指标，重新设计滤波器，通过频谱分析，多次试验确定最合适的滤波器。
课程设计学生日志
时间
设计内容
2013.6.24-2013.6.30
查阅资料，确定初步方案
2013.6.31-2013.7.2
编写代码，调试改错
2013.7.3-2013.7.5
[B,A]=cheby1(n1,1,WN,'s');%计算出阶数为n1、截止频率为WN、通带波纹最大衰减为1的数字低通滤波器，它的返回值a、b分别表示数字低通滤波器的系统函数的分子和分母的多项式系数
[BT,AT]=lp2bp(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5);%实现模数的映射
subplot(2,1,2);%将2个图画到一个平面
plot(f1,abs(y1));%f1为横坐标，abs(y1)为纵坐标作图
title('原始频谱图');xlabel('频率（Hz）');ylabel('幅度(dB)');
3.3加入干扰
这里我加入的干扰是：白噪声高频干扰，50Hz的电源线干扰。
3.3.1 白噪声
二设计原理
2.1Butterworth低通数字滤波器的设计
巴特沃斯低通滤波器的平方幅度响应为
其中，n为滤波器的阶数， 为低通滤波器的截止频率。 该滤波器具有 一些特殊的性质：
① 对所有的n,都有当 时, ；
②对所有的n,都有当 时, ；
③ 是的单调递减函数,即不会出现幅度响应的起伏；
④当 时,巴特沃斯滤波器趋向于理想的低通滤波器；
freqz(num,den,64);
y=filter(num,den,s);
figure(5);
subplot(2,1,1);#39;);
xlabel('时间（s）');ylabel('幅值（A）');
s3=fft(y,4000);
f1=100*(0:3999)/4000;
subplot(2,1,2);
plot(ws/pi,abshs);%ws/pi为横坐标abshs为纵坐标作图
title('加入高频干扰信号后的频谱图');xlabel('Hz');ylabel('幅度');通过观察加入白噪声后的时域图和频域图，将它与未加入白噪声进行比较
可以发现频谱图在0Hz时的幅度增加的很大，而且又在没有谱线的频率上竟然出现了频谱，这是由于白噪声在所有频率上都有频率造成的。
通过用s来代表加入白噪声后的信号，并进行数字滤波器的频率响应，对s中4000个频率点调用plot函数画出加入白噪声后的时域图，再对ws/pi,abshs调用plot函数画出加入白噪声后的频谱图，具体操作如下：
q=50*rand(4000,1); %产生4000行1列的位于（0，1）区间的随机数
s=a(:,1)+q;%产生的随机数与原始信号叠加赋给s
figure(4);%新建图像4
fs=1000;%采样频率
f11=10;%通带频率
f12=25;%阻带频率
Wp1=(f11/fs)*2*pi;%通带角频率
Wp2=(f12/fs)*2*pi;%阻带角频率
Omegap1=2*fs*tan(Wp1/2);%数字转化为模拟
Omegap2=2*fs*tan(Wp2/2);%数字转化为模拟
要求根据以上指标进行如下设计：
1、设计一个低通滤波器滤除心电信号中的高频信号干扰，对滤波前后的心电信号的频谱进行分析比较。
2、设计一个高通滤波器滤除心电信号中的低频干扰，对滤波前后的心电信号的频谱进行分析比较。
3、设计一个带阻滤波器（50Hz陷波器）滤除心电信号中的电源线干扰，对滤波前后的心电信号的频谱进行分析比较。
BW=Omegap2-Omegap1;%频带宽度
W0=Omegap1*Omegap2;
W00=sqrt(W0);
WP=1;
[n1,WN]=buttord(WP,WS,1,50,'s');%参数WP和WS分别是通带边界频率和阻带边界频率，通带最大衰减为1阻带最小衰减50，返回的参数n1和WN分别为滤波器的阶数和3dB截止频率
subplot(2,1,1);
plot(t,x);%t为横坐标，x为纵坐标作图
title('加入电源线干扰后的时域图');xlabel('时间（s）');ylabel('幅值（A）');
subplot(2,1,2);
plot(f2,abs(y2));%f2为横坐标，abs(y2)为纵坐标作图
title('加入电源线干扰后的频谱图');xlabel('幅');ylabel('Hz');
3.4 滤波器的设计
3.4.1切比雪夫I型数字低通滤波器
用Wp1,Wp2,Ws1,Ws2表示分别用通带和阻带截止频率的角频率算出频带宽带，计算阶数n1和截止频率WN，再设计切比雪夫I型模拟滤波器，采用双线性法将模拟滤波器系数变为数字滤波器系数，画出切比雪夫I型数字滤波器的频率响应，调用filter实现对白噪声的滤波，再最后调用plot函数画出滤除白噪声后的时域图和频域图。具体过程如下：
x1=a(:,1);%将矩阵a用x1表示
x=x1+x2;%x1和x2叠加赋给x
y2=fft(x,4000);%对x信号做快速傅里叶变换
f2=100*(0:3999)/4000;%先生成一个0,1,2,...,3999的整数向量，然后对对这个向量的每一项乘以100除以4000.
figure(3);%新建图像3
（1）确定数字低通滤波器的技术指标：通带截止频率ωp、通带衰减p、阻带截止频率ωs、阻带衰减s
切比雪夫滤波器的振幅平方特性如图2所示：
图2 切比雪夫滤波器的振幅平方特性
(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。如果采用脉冲响不变法，边界频率的转换关系为:
如果采用双线性变换法，边界频率的转换关系为
设计要求：
一般正常人的心电信号频率在0.7～100HZ范围内，幅度为10μv（胎儿）~5mv（成人）。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。