实验目的：

实验类型：设计性实验。
掌握通过零、极点设计低阶数字滤波器的设计方法。
掌握IIR和FIR数字滤波器的设计方法，并利用所设计滤波器解决实际问题。
通过分析滤波前后信号频谱的变化，深刻理解滤波的概念。
对两种滤波器设计方法的进行比较。

实验主要仪器设备，软件

硬件准备：PC机
软件准备：Matlab语言环境

实验的基本原理与内容

实验原理

数字滤波器是具有某种特定频率特性的线性时不变系统。广义上，任何线性时不变离散系统都是一个数字滤波器。设计数字滤波器的任务就是需求一个因果稳定的线性是不变系统，并使系统函数H(z)具有指定的频率特性。

低阶数字滤波器一般指一阶或二阶滤波器。因其阶次较低，可用零极点的分布来调控其频率特性；其好处处理速度快，硬件简单。另外，高阶滤波器在许多情况下由多个低阶滤波器组合来实现。
LTI系统可在z域中用零极点图的形式来描述。这在设计简单的滤波器时很重要，利用其可进行谱分析，只要正确地配置零极点就可达到一定的设计要求。那么，建立零极点与频率特性的关系规律：

设计滤波器时，一定要保证结构的稳定，因此所有极点应该位于单位圆之内；
负相位越大引起系统的延时越大，为了减少系统时延，所设计的系统零点最好也在单位圆内(或上)；
极点在单位圆附近，对应频率幅度会出现波峰；零点在单位圆上，对应频率幅度会出现波谷。

基于零极点的低阶数据滤波器模型：

数字滤波器从实现的网络结构或者从单位冲激响应分类，可以分为无限长单位冲激响应（IIR）数字滤波器和有限长单位冲激响应（FIR）数字滤波器。
数字滤波器设计的基本问题：

根据实际要求确定数字滤波器的性能指标。
用一个稳定的系统函数去逼近这个指标。
用一个有限精度的运算去实现这个传输函数。

其中IIR数字滤波器主要是借助于模拟滤波器的设计方法。模拟滤波器已有一套成熟的设计方法，充分利用这些资源为IIR数字滤波器设计提供很大方便。但IIR数字滤波器的相位特征一般为非线性的。而在有些实际应用场合，如图像处理，对滤波器线性相位特性要求颇为严格，因此需要选择FIR数字滤波器，且FIR数字滤波器的幅度特征可任意设计，比较IIR数字滤波器的设计更灵活。

实验内容（一）：

设有一连续信号x(t)=sin(15t)+cos(40t)+cos(200t)，基于零、极点低阶数字滤波器模型，分别设计二阶数字波器和一阶数字滤波器，分别提取三个频率分量。设采样周期T=0.005秒。
基本流程：
a）构造连续信号x(t)并显示波形，信号的作用时间设定为：0~10s;

1
2
3

>> t = 0:0.005:10;
>> x=sin(15*t)+cos(40*t)+cos(200*t);
>> plot(x);   %% 绘出原始信号

b）需求分析并选定滤波的类型【低通，高通，带通，带阻】；基于二阶数字波器和一阶数字滤波器的模型传递函数设计其参数，并画出零、极点分布图及滤波器的频响；
c）用设计出的滤波器分别提取x(t)信号各频率成分；并与参考信号相比。

b=[1 1];
a=[20 -18];
y0=filter(b,a,x);  %%用低通滤除cos（200t）
a1 = poly([0.9,0.9]);
b1 = poly([exp(0.2j),exp(-0.2j)]);
w0 = 0.2;
alpha1 = sin(w0);
k1 = (1+alpha1)/2 ;
b1 = k1.*b1;
y1=filter(b1,a1,y0); %%用带阻滤除cos（40t）
figure,subplot(3,1,1),zplane(b1,a1); %%绘制零极点
[H,w] = freqz(b1,a1);
subplot(3,1,2), plot(w,abs(H)); %%绘制滤波器   
subplot(3,1,3),plot(y1); %%绘制提取的sin（15t）
hold on ;
plot(sin(15*t),'g'); %%绘出sin（15t）

alpha2 = 0.98;
beta2 = cos(0.2);
a2 = [1 -beta2*(1+alpha2) alpha2];
k2 = (1-alpha2)/2;
b2 = k2.*[1 0 -1];
y3 =filter(b2,a2,x);
figure,subplot(3,1,1),zplane(b2,a2);%%绘制零极点
[H,w] = freqz(b2,a2);
subplot(3,1,2), plot(w,abs(H)); %%绘制滤波器
subplot(3,1,3),plot(y3);%%绘制提取的cos（40t）
hold on;
plot(cos(40*t),'g');   %%绘出cos（40t）

alpha3 = 0.98;
beta3 = cos(1);
a3 = [1 -beta3*(1+alpha3) alpha3];
k3 = (1-alpha2)/2;
b3 = k3.*[1 0 -1];
y4 =filter(b3,a3,x);
figure,subplot(3,1,1),zplane(b3,a3);%%绘制零极点
[H,w] = freqz(b3,a3);
subplot(3,1,2),plot(w,abs(H));%%绘制滤波器
subplot(3,1,3),plot(y4);  %%绘制提取的cos（200t）
hold on;
plot(cos(200*t),'g');     %%绘制cos（200t）
axis([0,300,-1.1,1.1]);%%设置显示范围
hold off;

实验内容（二）

实验背景:
心电( Electrocardiograph, ECG) 信号是心脏电活动在人体体表的表现，它一般比较微弱，正常的心电信号在0.01Hz-100Hz频带范围内。而90%的ECG频谱能量又集中在0.25-35Hz之间。因此, 心电信号在提取的过程中极易受到内、外环境的干扰, 而其中最而其中最主要的干扰是近场 50 Hz工频干扰源。因此, 如何消除 50 Hz工频干扰, 成为处理心电信号的首要任务。
实验数据：
实验数据是从美国麻省理工学院提供的MIT-BIH数据库下载的212文件格式带工频干扰的心电信号。该212文件格式包括头、数据和注释三个文件，以ASCll码方式存储。其采样率为360Hz，共有30000个样本点，数据文件xindiantu.xls。
实验任务：
分别设计IIR数字带阻滤波器（陷波器）和FIR数字带阻滤波器（陷波器），陷波器的模拟指标如下图所示：阻带从49-51Hz间，其中心频率为f 0=50 Hz，过渡带宽B=2Hz，通带的波动控制在1dB内，而阻带的衰减至少到45dB，采样频率为fs=360Hz；并将所设计滤波器对实际被50Hz干扰后的心电图信号进行滤波处理仿真实验，滤除其中的干扰成分。观察总结滤波作用与效果。

实验步骤：
（1）读取给定心电图信号的数据文件xidiantu.xls，并绘出心电图信号波形图；

将xindiantu.xls导入matlab，并转换成行向量data。

运行以下程序绘出原始心电图

1
2
3

T=1/360;
t=0:T:29999*T;
plot(t,data)；

（2）对心电图信号的频谱进行分析，得出滤波器的设计指标参数；并绘制其模拟域频谱图。
根据实验要求，所设计的滤波器阻带为49-51Hz，其中心频率为f 0=50 Hz，过渡带宽B=2Hz，通带的波动控制在1dB内，而阻带的衰减至少到45dB，采样频率为fs=360Hz；所以设计滤波器的阻带为60dB。
（3）根据所分析的指标参数，设计50Hz IIR数字陷波器（带阻滤波器）；

将数字域频率转化成模拟域频率
1
2
3
Rp=1;
As=60;
fs=360;

设定滤波器各项参数

wpl=fs*2*tan((47*pi/360));
wph=fs*2*tan((53*pi/360));
wsl=fs*2*tan((49*pi/360));
wsh=fs*2*tan((51*pi/360));
wp=[wpl,wph];
ws=[wsl,wsh];

设置巴特沃思(Butterworth)低通滤波器原型

[N,wc]=buttord(wp,ws,rp,as,'s');
[z,p,k]=buttap(N);
[b0,a0]=zp2tf(z,p,k);
[H,w] = freqs(b0,a0);
plot(w,abs(H));

低通滤波器转换成带阻滤波器

1
2
3

B=wc(2)-wc(1);
wz=(wc(2).*wc(1)).^(1/2);
[bt,at]=lp2bs(b0,a0,wz,B);

绘制滤波器形状

[b,a]=bilinear(bt,at,fs); 
subplot(1,2,1) ,zplane(b,a);
[H,w]=freqz(b,a);
subplot(1,2,2),plot(w*360/(2*pi),20*log10(abs(H)));

开始滤波

1
2
3

y=filter(b,a,data);
plot(t,data);
figure,plot(t,y);

(4)根据所分析的指标参数，设计50Hz FIR数字陷波器（带阻滤波器）；用给定指标进行验证设计的合理性。

将信号转到换数字域
1
2
3
rp=1;
As=60;
fs=360;

设置滤波器参数

wp1=47*2*pi/360;
wph=53*2*pi/360;
ws1=49*2*pi/360;
wsh=51*2*pi/360;

生成滤波器并绘制出滤波器形状

deltaw=2*2*pi/360;
N0=ceil(6.6*pi/deltaw);
N=N0+mod(N0+1,2);
wd=(hamming(N));
wc1=(ws1+wp1)/2; wc2=(wph+wsh)/2;
hd=[zeros(1,(N-1)/2),1,zeros(1,(N-1)/2)]-ideallp(wc2,N)+ideallp(wc1,N)
h=hd.*wd';
a=1;
H=freqz(h,a,0:0.01:pi);
subplot(1,2,1), zplane(h,a);
subplot(1,2,2), plot([0:0.01:pi]*360/(2*pi),abs(H)) ;

进行滤波

T=1/360;
t=0:T:29999*T;
y=filter(h,a,data);
plot(t,data);
figure,plot(t,y);

（5）分别利用两种所设计的数字陷波器对给定的心电图信号进行滤波，滤除50Hz的工频干扰；
（6）绘制滤波后信号的频谱图，并原始心电信号的频谱进行比对，分析结果；

（7）两种滤波器设计结果的比对，分析其原因。
IIR数字滤波器的设计方法简单方便，但其相位特征一般为非线性的。FIR数字滤波器，具有系统稳定、其幅度特征可任意设计等优点。