一、实验目的

1、加深对 IIR 数字滤波器设计方法和设计步骤的理解；

2、掌握用模拟滤波器原型设计 IIR 数字滤波器的方法；

3、能编写 MATLAB 函数，掌握设计 IIR 数字滤波器的函数调用方法；

4、根据不同的应用场景，确定不同的设计指标，设计出具有不同功能和性

能的滤波器。不同滤波器的设计方法具有不同的优缺点，因此要全面、客观看

待可能面对或出现的问题。

二、实验原理

2.1脉冲响应不变法的基本知识

用脉冲响应不变法设计 IIR 数字滤波器的步骤如下：

① 输入给定的数字滤波器的设计指标；

② 根据公式 Ω=ω/T 将数字滤波器设计指标转换为模拟滤波器设计指标；

③ 确定模拟滤波器的最小阶数和截止频率；

④ 计算模拟低通原型滤波器的系统传递函数；

⑤ 利用模拟域频率变换法求解实际模拟滤波器的系统传递函数；

⑥ 用脉冲响应不变法将模拟滤波器转换为数字滤波器。

 2.3双线性不变法的基本知识

双线性变换法克服了脉冲响应不变法从 s 平面到 z 平面的多值映射的缺点，

消除了频谱混叠现象。但其在变换过程中产生了非线性畸变，在设计 IIR 数字滤波器的过程中需要进行一定的修正。

用双线性变换法设计 IIR 数字滤波器的步骤如下：

① 输入给定的数字滤波器的设计指标；

② 根据公式 Ω=(2/T)tan(ω/2)进行预修正，将数字滤波器设计指标转换为模拟滤波器设计指标；

③ 确定模拟滤波器的最小阶数和截止频率；

④ 计算模拟低通原型滤波器的系统传递函数；

⑤ 利用模拟域频率变换法求解实际模拟滤波器的系统传递函数；

⑥ 用双线性变换法将模拟滤波器转换为数字滤波器。

2.6 信号的整数倍抽取

设 x(n)是连续信号 xa(t)的采样序列，其采样频率为 f1=1/T1(Hz)，T1 是采样

间隔。如果将其采样频率降低到原来的 1/D（D 为大于 1 的整数，称为抽取因

子），最简单的方法是对 x(n)每 D-1 个点抽取 1 点，组成一个新的序列 y(n)。由

于 y(n)的采样间隔 T2=DT1，除非抽取后仍能满足采样定理，否则会引起频谱混

叠现象。信号抽取前后的频谱关系见教材第 8 章的 8.2 节。为了避免抽取后的

频率混叠，在抽取前先采用一个抗混叠低通滤波器对信号滤波，把信号的频带

限制在某个频率以下。

抗混叠滤波器的系统函数为：

 三、实验程序及结果分析

第一题：

        1、用双线性变换法设计的巴特沃斯数字低通滤波器，要求：ωp=0.2π，Rp=1dB； 阻带：ωs=0.35π，As=15dB，滤波器采样频率 Fs=10Hz。

结果：

代码：

   %% 变量值
wp = 0.2*pi;
ws = 0.35*pi;
Rp = 1; As = 15;
ripple = 10^(-Rp/20);
Attn = 10^(-As/20);
Fs=10;T=1/Fs;
Omgp=(2/T)*tan(wp/2);%原型通带频率的预修正
Omgs=(2/T)*tan(ws/2);%原型阻带频率的预修正%% 模拟原型滤波器计算
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%计算阶数n和截止频率
[z0,p0,k0]=buttap(n);%设计归一化的巴特沃思模拟滤波器原型
bal=k0*real(poly(z0));%求原型滤波器的系数 b
aa1=real(poly(p0));%求原型滤波器的系数 a
[ba,aa]= lp2lp(bal,aa1,Omgc);%变换为模拟低通滤波器 %也可将以上4行替换为[bb,aa]=butter(n,Omgc,'s');直接求模拟滤波器系数
[bd,ad]= bilinear(ba,aa,Fs);%用双线性变换法计算数字滤波器系数
%% 绘制
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
figure(1);
subplot(2,2,1);plot(w/pi,abs(H));
ylabel("H");title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid;subplot(2,2,2);plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,3);plot(w/pi,dbH);title('幅度响应(dB)');
ylabel('dB');xlabel('频率(pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.35,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);gridsubplot(2,2,4);zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零极点图');%% 第二题
%序列
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...
0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...
4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
yn=filter(bd,ad,xn);
%% 绘制
figure(2);
% 绘制滤波前后的心电图波形
subplot(4, 1, 1);
stem(0:length(xn)-1, xn, 'filled');
title('滤波前心电图波形');
xlabel('时间 (采样点)');
ylabel('幅值');subplot(4, 1, 2);
stem(0:length(yn)-1, yn, 'filled');
title('滤波后心电图波形');
xlabel('时间 (采样点)');
ylabel('幅值');%% 绘制归一化频率幅频特性曲线
N = 1024; % FFT点数
Xf = fft(xn, N); % 原始信号的FFT
Yf = fft(yn, N); % 滤波后信号的FFT
k = 0:1023;
subplot(4, 1, 3);
plot(k*2/N, abs(Xf)); % 绘制幅频特性
title('滤波前幅频特性');
xlabel('\omega / \pi');
ylabel('幅值');
subplot(4, 1, 4);
plot(k*2/N, abs(Yf)); % 绘制幅频特性
title('滤波后幅频特性');
xlabel('\omega / \pi');
ylabel('幅值');

第二题 

        2、用 1 设计的数字滤波器对实际心电图信号采样序列(实验原理中已给出)进行滤波处理，分别绘制出滤波前后的心电图波形图和其幅频特性曲线，观察总结滤波作用与效果。

结果：

        

结果分析：

        心电信号(频率一般在0.05HZ~100Hz范围)是一种基本的人体生理信号抗干扰性差[由于心电信号的微弱，在采集过程中，极易受人体内和体外环境的影响，采集到的心电信号常常都伴随着强烈的噪声。从给定的心电信号序列可以看出该心电信号存在高频干扰。可以看出经过IIR低通滤波器，噪声很好的被滤除掉了，原始信号也失真度很小。

代码：

        主体：

[bd,ad]=myfilter();
%% 第二题
%序列
xn=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,...
0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,...
4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];
yn=filter(bd,ad,xn);
%% 绘制
figure(2);
% 绘制滤波前后的心电图波形
subplot(4, 1, 1);
stem(0:length(xn)-1, xn, 'filled');
title('滤波前心电图波形');
xlabel('时间 (采样点)');
ylabel('幅值');subplot(4, 1, 2);
stem(0:length(yn)-1, yn, 'filled');
title('滤波后心电图波形');
xlabel('时间 (采样点)');
ylabel('幅值');%% 绘制归一化频率幅频特性曲线
N = 1024; % FFT点数
Xf = fft(xn, N); % 原始信号的FFT
Yf = fft(yn, N); % 滤波后信号的FFT
k = 0:1023;
subplot(4, 1, 3);
plot(k*2/N, abs(Xf)); % 绘制幅频特性
title('滤波前幅频特性');
xlabel('\omega / \pi');
ylabel('幅值');
subplot(4, 1, 4);
plot(k*2/N, abs(Yf)); % 绘制幅频特性
title('滤波后幅频特性');
xlabel('\omega / \pi');
ylabel('幅值');

        封装的函数:

function [bd,ad] = myfilter()
%% 变量值
wp = 0.2*pi;
ws = 0.35*pi;
Rp = 1; As = 15;
ripple = 10^(-Rp/20);
Attn = 10^(-As/20);
Fs=10;T=1/Fs;
Omgp=(2/T)*tan(wp/2);%原型通带频率的预修正
Omgs=(2/T)*tan(ws/2);%原型阻带频率的预修正%% 模拟原型滤波器计算
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%计算阶数n和截止频率
[z0,p0,k0]=buttap(n);%设计归一化的巴特沃思模拟滤波器原型
bal=k0*real(poly(z0));%求原型滤波器的系数 b
aa1=real(poly(p0));%求原型滤波器的系数 a
[ba,aa]= lp2lp(bal,aa1,Omgc);%变换为模拟低通滤波器 %也可将以上4行替换为[bb,aa]=butter(n,Omgc,'s');直接求模拟滤波器系数
[bd,ad]= bilinear(ba,aa,Fs);%用双线性变换法计算数字滤波器系数
end

 第三题：

3、设计一个抗混叠低通滤波器（可在实验内容 1 的代码上进行修改，截止频率的指标见 2.6 节，衰减指标与实验内容 1 一样）。（1）读取音频信号motherland.wav，得到 xn；（2）对 xn 进行 D=2 的整数倍抽取，得到整数倍抽取后的音频信号 yn1；（3）对 xn 先进行抗混叠滤波，再进行 D=2 的整数倍抽取，得到音频信号 yn2。

结果：

实验结果分析：
% 滤波前，pi附近还有频率分量，滤波后，pi 附近没有频率分量，

% 滤波前，突变，高频成分多，滤波后，光滑，高频成分减少，低频成分多

总结：滤波后的信号更容易恢复，失真低。

代码：

%% 参数初始化
D = 2; % 抽取因子
[xn, fs] = audioread("motherland.wav");
sound(xn, fs);%% 第一问：直接抽取
yn1 = xn(1:D:end); % 每 D 点抽取 1 点
sound(yn1, fs / D);%% 第二问：频谱对比
N = 2018;
figure(1);
Xn = 1/fs * fft(xn(8000:8199), N);
subplot(3, 1, 1);
plot((0:N/2-1)*fs/N, abs(Xn(1:N/2))); % 原始信号频谱
title('原始信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度');Yn1 = D/fs * fft(yn1(4000:4099), N);
subplot(3, 1, 2);
plot((0:N/2-1)*fs/(N*D), abs(Yn1(1:N/2))); % 抽取信号频谱
title('直接抽取信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度');%% 第三问：抗混叠滤波器设计和应用
wc = pi / D; % 截止频率
wp = wc * 0.9999; % 通带截止频率
ws = wc * 1.04; % 阻带截止频率，确保阻带超出 Nyquist
Rp = 1; % 通带衰减 (dB)
As = 15; % 阻带衰减 (dB)
Fs = fs; % 原始采样频率
T = 1 / Fs;% 模拟频率转换
Omgp = (2 / T) * tan(wp / 2);
Omgs = (2 / T) * tan(ws / 2);% Butterworth 滤波器设计
[n, Omgc] = buttord(Omgp, Omgs, Rp, As, 's');
[z0, p0, k0] = buttap(n); % 原型低通滤波器
bal = k0 * real(poly(z0));
aa1 = real(poly(p0));
[ba, aa] = lp2lp(bal, aa1, Omgc); % 转换为目标低通滤波器
[bd, ad] = bilinear(ba, aa, Fs); % 双线性变换离散化% 抗混叠滤波器频率响应
figure(2);
subplot(2, 1, 1);
[h, w] = freqz(bd, ad, 1024, Fs);
plot(w, abs(h), 'LineWidth', 1.5);
title('抗混叠滤波器幅频特性');
xlabel('频率 (Hz)');
ylabel('幅度');
grid on;subplot(2, 1, 2);
plot(w, angle(h)*180/pi, 'LineWidth', 1.5);
title('抗混叠滤波器相频特性');
xlabel('频率 (Hz)');
ylabel('相位 (度)');
grid on;% 抗混叠滤波
yn2 = filter(bd, ad, xn); % 滤波
yn2_downsampled = yn2(1:D:end); % 滤波后抽取
sound(yn2_downsampled, fs / D);% 滤波后信号频谱
Yn2 = D/fs * fft(yn2_downsampled(4000:4099), N);
figure(1);
subplot(3, 1, 3);
plot((0:N/2-1)*fs/(N*D), abs(Yn2(1:N/2)));
title('抗混叠滤波后信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度');

思考题：

结果：

代码：

% 通用参数设置
alpha_p = 1;  % 通带衰减 (dB)
alpha_s = 45; % 阻带衰减 (dB)% ----------------- 低通滤波器设计 -----------------
wp_lp = 0.2 * pi; % 低通通带截止频率
ws_lp = 0.3 * pi; % 低通阻带截止频率
[n_lp, wn_lp] = buttord(wp_lp/pi, ws_lp/pi, 1, 20); % 计算阶数
[b_lp, a_lp] = butter(n_lp, wn_lp, 'low'); % 低通滤波器
[h_lp, w_lp] = freqz(b_lp, a_lp, 1024); % 频率响应% 绘图：低通滤波器
figure;
subplot(3,1,1);
plot(w_lp/pi, abs(h_lp), 'LineWidth', 1.5);
title('低通滤波器 - 幅度响应');
xlabel('归一化频率 (×π rad/sample)');
ylabel('|H|');
grid on;subplot(3,1,2);
plot(w_lp/pi, 20*log10(abs(h_lp)), 'LineWidth', 1.5);
title('低通滤波器 - 幅度响应 (dB)');
xlabel('归一化频率 (×π rad/sample)');
ylabel('幅度 (dB)');
grid on;subplot(3,1,3);
plot(w_lp/pi, angle(h_lp)*180/pi, 'LineWidth', 1.5);
title('低通滤波器 - 相频特性');
xlabel('归一化频率 (×π rad/sample)');
ylabel('相位 (度)');
grid on;% ----------------- 高通滤波器设计 -----------------
wp_hp = 0.6 * pi; % 高通通带截止频率
ws_hp = 0.4 * pi; % 高通阻带截止频率
[n_hp, wn_hp] = buttord(wp_hp/pi, ws_hp/pi, 2, 30); % 计算阶数
[b_hp, a_hp] = butter(n_hp, wn_hp, 'high'); % 高通滤波器
[h_hp, w_hp] = freqz(b_hp, a_hp, 1024); % 频率响应% 绘图：高通滤波器
figure;
subplot(3,1,1);
plot(w_hp/pi, abs(h_hp), 'LineWidth', 1.5);
title('高通滤波器 - 幅度响应');
xlabel('归一化频率 (×π rad/sample)');
ylabel('|H|');
grid on;subplot(3,1,2);
plot(w_hp/pi, 20*log10(abs(h_hp)), 'LineWidth', 1.5);
title('高通滤波器 - 幅度响应 (dB)');
xlabel('归一化频率 (×π rad/sample)');
ylabel('幅度 (dB)');
grid on;subplot(3,1,3);
plot(w_hp/pi, angle(h_hp)*180/pi, 'LineWidth', 1.5);
title('高通滤波器 - 相频特性');
xlabel('归一化频率 (×π rad/sample)');
ylabel('相位 (度)');
grid on;% ----------------- 带通滤波器设计 -----------------
wp_bp = [0.2*pi, 0.6*pi]; % 带通通带范围
ws_bp = [0.15*pi, 0.65*pi]; % 带通阻带范围
[n_bp, wn_bp] = buttord(wp_bp/pi, ws_bp/pi, alpha_p, alpha_s); % 计算阶数
[b_bp, a_bp] = butter(n_bp, wn_bp, 'bandpass'); % 带通滤波器
[h_bp, w_bp] = freqz(b_bp, a_bp, 1024); % 频率响应% 绘图：带通滤波器
figure;
subplot(3,1,1);
plot(w_bp/pi, abs(h_bp), 'LineWidth', 1.5);
title('带通滤波器 - 幅度响应');
xlabel('归一化频率 (×π rad/sample)');
ylabel('|H|');
grid on;subplot(3,1,2);
plot(w_bp/pi, 20*log10(abs(h_bp)), 'LineWidth', 1.5);
title('带通滤波器 - 幅度响应 (dB)');
xlabel('归一化频率 (×π rad/sample)');
ylabel('幅度 (dB)');
grid on;subplot(3,1,3);
plot(w_bp/pi, angle(h_bp)*180/pi, 'LineWidth', 1.5);
title('带通滤波器 - 相频特性');
xlabel('归一化频率 (×π rad/sample)');
ylabel('相位 (度)');
grid on;% ----------------- 带阻滤波器设计 -----------------
wp_bs = [0.15*pi, 0.65*pi]; % 带阻通带范围
ws_bs = [0.2*pi, 0.6*pi]; % 带阻阻带范围
[n_bs, wn_bs] = buttord(wp_bs/pi, ws_bs/pi, alpha_p, alpha_s); % 计算阶数
[b_bs, a_bs] = butter(n_bs, wn_bs, 'stop'); % 带阻滤波器
[h_bs, w_bs] = freqz(b_bs, a_bs, 1024); % 频率响应% 绘图：带阻滤波器
figure;
subplot(3,1,1);
plot(w_bs/pi, abs(h_bs), 'LineWidth', 1.5);
title('带阻滤波器 - 幅度响应');
xlabel('归一化频率 (×π rad/sample)');
ylabel('|H|');
grid on;subplot(3,1,2);
plot(w_bs/pi, 20*log10(abs(h_bs)), 'LineWidth', 1.5);
title('带阻滤波器 - 幅度响应 (dB)');
xlabel('归一化频率 (×π rad/sample)');
ylabel('幅度 (dB)');
grid on;subplot(3,1,3);
plot(w_bs/pi, angle(h_bs)*180/pi, 'LineWidth', 1.5);
title('带阻滤波器 - 相频特性');
xlabel('归一化频率 (×π rad/sample)');
ylabel('相位 (度)');
grid on;