（4） 数字带阻滤波器设计 用MATLAB 直接法设计切比雪夫I型数字带阻波器，要求: fp1=1kHz，fp2=4.5kHz，Rp=1dB;fs1=2kHz,fs2=3.5kHz，As=20dB，滤波器采样频率Fs=10 kHz。请描绘滤波器绝对和相对幅频特性、相频特性、零极点分布图，列出系统传递函数式。

以下是MATLAB代码实现数字带阻滤波器的设计：

% 设计参数
Rp = 1;     % 通带最大衰减（dB）
As = 20;    % 阻带最小衰减（dB）
fp1 = 1000; % 通带下边缘频率（Hz）
fp2 = 4500; % 通带上边缘频率（Hz）
fs1 = 2000; % 阻带下边缘频率（Hz）
fs2 = 3500; % 阻带上边缘频率（Hz）
Fs = 10000; % 采样频率（Hz）

% 计算数字滤波器参数
wp1 = 2*pi*fp1/Fs;
wp2 = 2*pi*fp2/Fs;
ws1 = 2*pi*fs1/Fs;
ws2 = 2*pi*fs2/Fs;

delta_p = (10^(Rp/20)-1)/(10^(Rp/20)+1);
delta_s = 10^(-As/20);

% 计算滤波器阶数
n = ceil(acosh((sqrt(delta_s)-1)/(sqrt(delta_p)))/acosh(ws1/wp1));

% 计算传递函数
epsilon = sqrt(delta_p);
k = 1:n;
p = -sinh((1i*pi*(2*k-1))/(2*n))*sin(wp1);
z = -1./conj(p);
b0 = epsilon^(-1)*poly(z);
a0 = poly(p);

% 绘制幅频特性
[H, w] = freqz(b0, a0, 2048);
figure;
plot(w/pi*Fs/2, abs(H));
hold on;
plot([0 fp1 fp1 fp2 fp2 fs2 fs2 1]*Fs/2,  [-delta_p -delta_p 0 0 -delta_p -delta_p 0 0], 'r');
plot([0 fp1 fp1 fp2 fp2 fs2 fs2 1]*Fs/2,  [-delta_p-delta_s -delta_p-delta_s -delta_s -delta_s -delta_p-delta_s -delta_p-delta_s 0 0], 'm');
xlim([0 Fs/2]);
xlabel('频率 (Hz)');
ylabel('幅度');
title('数字带阻滤波器幅频特性');
legend('滤波器幅频特性', '通带', '阻带');

% 绘制相频特性
figure;
plot(w/pi*Fs/2, angle(H));
xlim([0 Fs/2]);
xlabel('频率 (Hz)');
ylabel('相位');
title('数字带阻滤波器相频特性');

% 绘制零极点分布图
figure;
zplane(b0, a0);
title('数字带阻滤波器零极点分布图');

% 打印系统传递函数式
disp('系统传递函数：');
disp(['H(z) = ' num2str(b0(1)) ' + ' num2str(b0(2)) 'z^(-1) + ' num2str(b0(3)) 'z^(-2) + ...']);
disp(['        -------------------------------------------------------']);
disp(['        ' num2str(a0(1)) ' + ' num2str(a0(2)) 'z^(-1) + ' num2str(a0(3)) 'z^(-2) + ...']);

运行后，将会得到数字带阻滤波器的幅频特性、相频特性和零极点分布图，以及系统传递函数式。其中，幅频特性图中红色为通带，紫色为阻带，蓝色为滤波器的幅频特性。相位特性图中展示了滤波器的相位响应。零极点分布图展示了滤波器的零点和极点的位置。系统传递函数式中包含了滤波器的传递函数。