实验六：用 MATLAB 设计 IIR 数字滤波器

实验内容

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

%% 双线性变换法

% wp = 0.25 * pi;

% Rp = 1; % dB

% ws = 0.4 * pi;

% As = 15; % dB

% Fs = 100; % Hz

wp = 0.2 * pi;

Rp = 1; % dB

ws = 0.35 * pi;

As = 15; % dB

Fs = 10; % Hz

[n, Omgc, B, A] = Bilinear_transformation(wp, Rp, ws, As, Fs);

n, Omgc, B, A

n = 4

Omgc = 7.9909

B = 1×5

    0.0092    0.0367    0.0550    0.0367    0.0092

A = 1×5

    1.0000   -2.0325    1.8204   -0.7706    0.1294

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

h = figure;

subplot(2, 2, 1);

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];

stem(0:length(xn) - 1, xn, '.');

title('滤波前');

subplot(2, 2, 3);

N = 1024;

k = 0:N - 1;

y1 = fft(xn, N);

plot(k / N * 2, abs(y1));

title('滤波前的幅频特性');

xlabel('\omega/\pi');

axis([0 2 0 500]);

subplot(2, 2, 2);

yn1 = filter(B, A, xn);

stem(0:length(yn1) - 1, yn1, '.');

title('滤波后');

subplot(2, 2, 4);

N = 1024;

k = 0:N - 1;

y2 = fft(yn1, N);

plot(k / N * 2, abs(y2));

title('滤波后的幅频特性');

xlabel('\omega/\pi');

axis([0 2 0 500]);

saveas(h, '双线性变换法-低通滤波器', 'svg');

%% 脉冲响应不变法

%% 1、

% wp = 0.25 * pi;

% Rp = 1; % dB

% ws = 0.4 * pi;

% As = 15; % dB

% Fs = 2000; % Hz

wp = 0.2 * pi;

Rp = 1; % dB

ws = 0.35 * pi;

As = 15; % dB

Fs = 10; % Hz

[n, Omgc, bd, ad] = Impulse_response_invariant_method(wp, Rp, ws, As, Fs);

n, Omgc, bd, ad

n = 5

Omgc = 7.8093

bd = 1×5

    0.0000    0.0071    0.0457    0.0277    0.0016

ad = 1×6

    1.0000   -2.5686    2.9696   -1.8328    0.5937   -0.0799

%% 2、

h = figure;

subplot(2, 2, 1);

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];

stem(0:length(xn) - 1, xn, '.');

title('滤波前');

subplot(2, 2, 2);

N = 1024;

k = 0:N - 1;

y1 = fft(xn, N);

plot(k, abs(y1));

title('滤波前的幅频特性');

% axis([0 1100 0 500]);

subplot(2, 2, 3);

yn1 = filter(bd, ad, xn);

stem(0:length(yn1) - 1, yn1, '.');

title('滤波后');

subplot(2, 2, 4);

N = 1024;

k = 0:N - 1;

y2 = fft(yn1, N);

plot(k, abs(y2));

title('滤波后的幅频特性');

% axis([0 1100 0 500]);

saveas(h, '脉冲响应不变法-低通滤波器', 'svg');

3、设计一个抗混叠低通滤波器（可在实验内容 1 的代码上进行修改，截止频率的指标见 2.6 节，衰减指标与实验内容 1 一样）。

（1）读取音频信号motherland.wav，得到 xn；

（2）对 xn 进行 D=2 的整数倍抽取，得到整数倍抽取后的音频信号 yn1；

（3）对 xn 先进行抗混叠滤波，再进行 D=2 的整数倍抽取，得到音频信号 yn2。

（1）音频播放：依次原始声音 xn、没有经抗混叠滤波进行整数倍抽取的音频 yn1、经过抗混叠滤波进行整数倍抽取的音频 yn2，体验音频有频域混叠时的音质。

参考代码：

[xn,fs]=audioread('motherland.wav');% 读取音频信号

sound(xn,fs);

pause(length(xn)/fs);% 播放音频信号，暂停执行程序 length(xn)/fs 秒

yn1=xn(1:D:length(xn)); % 每个 D-1 个点抽取 1 点，这里 D=2

sound(yn1,fs/D); % 采用频谱降低到 fs/D

（2）取原音频某段信号，如 n=8000~8199。画出该段信号模拟域幅度谱（横坐标为 f Hz）；画出该段信号 D=2 抽取后的模拟域幅度谱；画出该段信号先经过抗混叠滤波再进行 D=2 抽取的模拟域幅度谱。

参考代码：

Xn=1/fs*fft(xn(8000:8199),N); % 从xn中取200点，N可取2018

plot((0:N/2-1)*fs/N,abs(Xn(1:N/2))); % 模拟域幅度谱

Yn1=D/fs*fft(yn1(8000:8099),N); % 2点取1点后，200点长变成了100点长

plot((0:N/2-1)*fs/(N*D),abs(Yn1(1:N/2)));% 模拟域幅度谱

%% 1、读取音频信号

D = 2;

play = 0;

[xn, fs] = audioread('motherland.wav'); % 读取音频信号

if play == 1

sound(xn, fs); % 播放音频信号

pause(length(xn) / fs); % 暂停执行程序 length(xn)/fs 秒

end

%% 2、整数倍抽取

yn1 = xn(1: D: length(xn)); % 每个 D-1 个点抽取 1 点，这里 D=2

if play == 1

sound(yn1, fs / D); % 播放音频信号，采用频谱降低到 fs/D

pause(length(yn1) * D / fs); % 暂停执行程序 length(xn)*D/fs 秒

end

%% 3.1、抗混叠滤波；3.2、进行 D=2 的整数倍抽取

yn1_filter = filter(B, A, xn); % 低通滤波

yn2 = yn1_filter(1: D: length(yn1_filter)); % 抽取

if play == 1

sound(yn2, fs / D); % 播放音频信号通过抗混叠低通滤波器后的音乐

pause(length(yn2) * D / fs); % 暂停执行程序 length(xn)*D/fs 秒

end

h = figure;

N = 2048;

begin = 8000;

ending = 8199;

L = ending - begin + 1;

subplot(3, 1, 1);

Xn = 1 / fs * fft(xn(begin: begin + L - 1), N);

plot((0:N / 2 - 1) * fs / N, abs(Xn(1: N / 2))); % 1. 4kHz

xlabel('f/Hz');

% plot((0:N / 2 - 1) * 2 / N, abs(Xn(1: N / 2))); % 2. pi

% xlabel('\omega/\pi');

% plot((0:N - 1) * fs / N, abs(Xn(1: N))); % 3. 8kHz

% xlabel('f/Hz');

% plot((0:N - 1) * 2 / N, abs(Xn(1: N))); % 4. 2pi

% xlabel('\omega/\pi');

title('原信号-频率响应');

subplot(3, 1, 2);

Yn1 = D / fs * fft(yn1(begin: begin + L / 2 - 1), N);

plot((0:N / 2 - 1) * fs / (N * D), abs(Yn1(1: N / 2))); % 1. 2kHz

xlabel('f/Hz');

% plot((0:N / 2 - 1) * 2 / (N * D), abs(Yn1(1: N / 2))); % 2. pi

% xlabel('\omega/\pi');

% plot((0:N - 1) * fs / (N * D), abs(Yn1(1: N))); % 3. 4kHz

% xlabel('f/Hz');

% plot((0:N - 1) * 2 / (N * D), abs(Yn1(1: N))); % 4. 2pi

% xlabel('\omega/\pi');

title('原信号2点采样-频率响应');

subplot(3, 1, 3);

Xn2 = D / fs * fft(yn2(begin: begin + L / 2 - 1), N);

plot((0:N / 2 - 1) * fs / (N * D), abs(Xn2(1: N / 2))); % 1. 2kHz

xlabel('f/Hz');

% plot((0:N / 2 - 1) * 2 / (N * D), abs(Xn2(1: N / 2))); % 2. pi

% xlabel('\omega/\pi');

% plot((0:N - 1) * fs / (N * D), abs(Xn2(1: N))); % 3. 4kHz

% xlabel('f/Hz');

% plot((0:N - 1) * 2 / (N * D), abs(Xn2(1: N))); % 4. 2pi

% xlabel('\omega/\pi');

title('原信号先经抗频谱混叠滤波再2点采样-频率响应');

saveas(h, '频率响应', 'svg');

% 滤波前，pi附近还有频率分量，滤波后，pi 附近没有频率分量，

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

% 实验：理解理论

思考题

按照如下指标要求设计四种选频数字滤波器，要求画出滤波器的幅频特性、相频特性和幅度衰减曲线，标注相关信息，如横坐标，纵坐标的单位，曲线名称等。（设计方法自己查阅资料完成）

（1）设计数字低通滤波器，指标为：通带截止频率 𝜔𝑝 = 0.2𝜋，阻带截止频率 𝜔𝑠 = 0.3𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 20𝑑𝐵。

（2）设计数字高通滤波器，指标为：阻带截止频率 𝜔𝑠 = 0.4𝜋，通带截止频率 𝜔𝑠 = 0.6𝜋，通带衰减 𝛼𝑝 = 2𝑑𝐵，阻带衰减 𝛼𝑠 = 30𝑑𝐵。

（3）设计数字带通滤波器，指标为：通带范围 0.2𝜋 ≤ 𝜔 ≤ 0.6𝜋，阻带范围 0 ≤ 𝜔 ≤ 0.15𝜋 和 0.65𝜋 ≤ 𝜔 ≤ 𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 45𝑑𝐵。

（4）设计数字带阻滤波器，指标为：阻带范围 0.2𝜋 ≤ 𝜔 ≤ 0.6𝜋，通带范围 0 ≤ 𝜔 ≤ 0.15𝜋 和 0.65𝜋 ≤ 𝜔 ≤ 𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 45𝑑𝐵。

%% （1）LPF

wp = 0.2 * pi;

ws = 0.3 * pi;

Rp = 1; % dB

As = 20; % dB

Fs = 10; % Hz

[n1, Omgc1, B1, A1] = LPF(wp, Rp, ws, As, Fs);

n1, Omgc1, B1, A1

n1 = 7

Omgc1 = 7.3392

B1 = 1×8

    0.0002    0.0012    0.0037    0.0062    0.0062    0.0037    0.0012    0.0002

A1 = 1×8

    1.0000   -3.8471    6.7860   -6.9629    4.4452   -1.7542    0.3944   -0.0388

%% （2）HPF

wp = 0.6 * pi;

ws = 0.4 * pi;

Rp = 2; % dB

As = 30; % dB

Fs = 10; % Hz

[n, Omgc, B2, A2] = HPF(wp, Rp, ws, As, Fs);

n, Omgc, B2, A2

n = 6

Omgc = 25.8378

B2 = 1×7

    0.0129   -0.0774    0.1935   -0.2581    0.1935   -0.0774    0.0129

A2 = 1×7

    1.0000    0.9575    1.1203    0.5115    0.2083    0.0376    0.0038

%% （3）BPF

wp1 = 0.15 * pi;

ws1 = 0.20 * pi;

ws2 = 0.60 * pi;

wp2 = 0.65 * pi;

Rp = 1; % dB

As = 45; % dB

Fs = 10; % Hz

% [n, Wc, B3, A3] = BPF([ws1, ws2], Rp, [wp1, wp2], As, Fs);

% n, Wc, B3, A3

[n, Wc, B3, A3] = BPF([wp1, wp2], Rp, [ws1, ws2], As, Fs);

n, Wc, B3, A3

n = 23

Wc = 1×2

    0.1738    0.6439

B3 = 1×47

    0.0000         0   -0.0000         0    0.0001         0   -0.0009         0    0.0043         0   -0.0164         0    0.0493         0   -0.1197         0    0.2394         0   -0.3990         0    0.5586         0   -0.6602         0    0.6602         0   -0.5586         0    0.3990         0   -0.2394         0    0.1197         0   -0.0493         0    0.0164         0   -0.0043         0    0.0009         0   -0.0001         0    0.0000         0   -0.0000

A3 = 1×47

10³ ×

    0.0010   -0.0093    0.0433   -0.1360    0.3304   -0.6710    1.1923   -1.8999    2.7552   -3.6772    4.5569   -5.2767    5.7359   -5.8748    5.6870   -5.2158    4.5400   -3.7560    2.9568   -2.2168    1.5835   -1.0780    0.6996   -0.4328    0.2550   -0.1431    0.0764   -0.0388    0.0187   -0.0086    0.0037   -0.0015    0.0006   -0.0002    0.0001   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000

%% （4）BSF

ws1 = 0.15 * pi;

wp1 = 0.20 * pi;

wp2 = 0.60 * pi;

ws2 = 0.65 * pi;

Rp = 1; % dB

As = 45; % dB

Fs = 10; % Hz

% [n, Wc, B4, A4] = BSF([ws1, ws2], Rp, [wp1, wp2], As, Fs);

% n, Wc, B4, A4

[n, Wc, B4, A4] = BSF([wp1, wp2], Rp, [ws1, ws2], As, Fs);

n, Wc, B4, A4

n = 23

Wc = 1×2

    0.1965    0.6057

B4 = 1×47

10⁴ ×

    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0001    0.0003   -0.0011    0.0034   -0.0093    0.0229   -0.0518    0.1083   -0.2095    0.3772   -0.6339    0.9966   -1.4693    2.0351   -2.6518    3.2543   -3.7646    4.1075   -4.2284    4.1075   -3.7646    3.2543   -2.6518    2.0351   -1.4693    0.9966   -0.6339    0.3772   -0.2095    0.1083   -0.0518    0.0229   -0.0093    0.0034   -0.0011    0.0003   -0.0001    0.0000   -0.0000    0.0000   -0.0000    0.0000

A4 = 1×47

10⁴ ×

    0.0001   -0.0010    0.0056   -0.0210    0.0615   -0.1497    0.3145   -0.5854    0.9821   -1.5034    2.1199   -2.7727    3.3827   -3.8663    4.1546   -4.2090    4.0293   -3.6514    3.1367   -2.5572    1.9799   -1.4567    1.0187   -0.6772    0.4279   -0.2568    0.1464   -0.0791    0.0406   -0.0197    0.0090   -0.0039    0.0016   -0.0006    0.0002   -0.0001    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000   -0.0000    0.0000

观察信号经过滤波器前后

%% 原信号

N = 2048;

begin = 8000;

ending = 8799;

L = ending - begin + 1;

[xn, fs] = audioread('motherland.wav'); % 读取音频信号

xn = xn(begin: begin + L - 1);

h = figure;

subplot(211);

plot(xn);

title('原信号');

subplot(212);

Xn = 1 / fs * fft(xn, N);

plot((0:N - 1) * 2 / N, abs(Xn(1: N))); % 4. 2pi

xlabel('\omega/\pi');

title('原信号-频率响应');

saveas(h, '原信号-频率响应', 'svg');

%% LPF信号

h = figure;

subplot(221);

xn_LPF = filter(B1, A1, xn); % 低通滤波

Xn_LPF = 1 / fs * fft(xn_LPF, N);

plot(xn_LPF);

title('经过LPF的信号');

subplot(222);

plot((0:N - 1) * 2 / N, abs(Xn_LPF(1: N)));

xlabel('\omega/\pi');

title('经过LPF的信号-频率响应');

%% HPF信号

subplot(223);

xn_HPF = filter(B2, A2, xn); % 高通滤波

Xn_HPF = 1 / fs * fft(xn_HPF, N);

plot(xn_HPF);

title('经过HPF的信号');

subplot(224);

plot((0:N - 1) * 2 / N, abs(Xn_HPF(1: N)));

xlabel('\omega/\pi');

title('经过HPF的信号-频率响应');

saveas(h, '经过LPF,HPF的信号-频率响应', 'svg');

%% BPF信号

h = figure;

subplot(221);

xn_BPF = filter(B3, A3, xn); % 带通滤波

Xn_BPF = 1 / fs * fft(xn_BPF, N);

plot(xn_BPF);

title('经过BPF的信号');

subplot(222);

plot((0:N - 1) * 2 / N, abs(Xn_BPF(1: N)));

xlabel('\omega/\pi');

title('经过BPF的信号-频率响应');

%% BSF信号

subplot(223);

xn_BSF = filter(B4, A4, xn); % 带阻滤波

Xn_BSF = 1 / fs * fft(xn_BSF, N);

plot(xn_BSF);

title('经过BSF的信号');

subplot(224);

plot((0:N - 1) * 2 / N, abs(Xn_BSF(1: N)));

xlabel('\omega/\pi');

title('经过BSF的信号-频率响应');

saveas(h, '经过PBF,BSF的信号-频率响应', 'svg');

function [n, Omgc, bd, ad] = Impulse_response_invariant_method(wp, Rp, ws, As, Fs)
%     wp = 0.25 * pi;                                         % 滤波器的通带截止频率
%     ws = 0.4 * pi;                                          % 滤波器的阻带截止频率
%     Rp = 1; As = 15;                                        % 滤波器的通阻带衰减指标
    ripple = 10^(-Rp / 20);                                 % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);                                   % 滤波器的阻带衰减对应的幅度值
                                                            % 转换为模拟滤波器的技术指标
%     Fs = 2000;
%     T = 1 / Fs; 
    Omgp = wp * Fs;
    Omgs = ws * Fs;                                         % 模拟原型滤波器计算
    [n, Omgc] = buttord(Omgp, Omgs, Rp, As, 's');           % 计算阶数n和截止频率
    [z0, p0, k0] = buttap(n);                               % 设计归一化的巴特沃思模逆滤波器原型，零点，极点和增益
    ba1 = k0 * real(poly(z0));                              % 求原型滤波器的系数b
    aa1 = real(poly(p0));                                   % 求原型滤波器的系数a
    [ba, aa] = lp2lp(ba1, aa1, Omgc);                       % 变换为模拟低通滤波器
    
    % 用脉冲响应不变法计算数字滤波器;系数
    [bd, ad] = impinvar(ba, aa, Fs);                        % 求数字系统的频率特性
    [H, w] = freqz(bd, ad);                                 % 求数字系统的频率特性
    dbH = 20 * log10((abs(H) + eps) / max(abs(H)));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, abs(H));
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 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, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
    ylabel('dB'); xlabel('频率(\pi)'); axis([0, 1, -40, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    
    subplot(2, 2, 4); zplane(bd, ad); axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, '脉冲响应不变法', 'svg');
end
function [n, Omgc, bd, ad] = Bilinear_transformation(wp, Rp, ws, As, Fs)
%     wp = 0.25 * pi;                                % 滤波器的通带截止频率
%     ws = 0.4 * pi;                                 % 滤波器的阻带截止频率
%     Rp = 1; As = 15;                               % 滤波器的通阻带衰减指标
    ripple = 10^(-Rp / 20);                        % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);                          % 滤波器的阻带衰减对应的幅度值
                                                   % 转换为模拟滤波器的技术指标
%     Fs = 100; 
    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
    aal = real(poly(p0));                          % 求原型滤波器的系数a
    [ba, aa] = lp2lp(bal, aal, Omgc);              % 变换为模拟低通滤波器
    %也可将以上4行替换为[ba,aa]=butter(n,Omgc,');直接求模拟滤波器系数
                                                   % 用双线性变换法计算数字滤波器系数
    [bd, ad] = bilinear(ba, aa, Fs);               % 求数字系统的频率特性
    [H, w] = freqz(bd, ad);
    dbH = 20 * log10((abs(H) + eps) / max(abs(H)));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, abs(H));
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 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, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
    ylabel('dB'); xlabel('频率(pi)'); axis([0, 1, -40, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    
    subplot(2, 2, 4); zplane(bd, ad);
    axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, '双线性变换法', 'svg');
end
function [n, Omgc, bd, ad] = LPF(wp, Rp, ws, As, Fs)
%     wp = 0.25 * pi;                                % 滤波器的通带截止频率
%     ws = 0.4 * pi;                                 % 滤波器的阻带截止频率
%     Rp = 1; As = 15;                               % 滤波器的通阻带衰减指标
    ripple = 10^(-Rp / 20);                        % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);                          % 滤波器的阻带衰减对应的幅度值
                                                   % 转换为模拟滤波器的技术指标
    % 使用修正值
    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
    aal = real(poly(p0));                          % 求原型滤波器的系数a
    [ba, aa] = lp2lp(bal, aal, Omgc);              % 变换为模拟低通滤波器
    %也可将以上4行替换为[ba,aa]=butter(n,Omgc,');直接求模拟滤波器系数
                                                   % 用双线性变换法计算数字滤波器系数
    [bd, ad] = bilinear(ba, aa, Fs);               % 求数字系统的频率特性
    [H, w] = freqz(bd, ad);
    dbH = 20 * log10((abs(H) + eps) / max(abs(H)));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, abs(H));
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 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, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
%     p = find(abs(w - wp / pi) < eps);
%     plot(w(p), dbH(p), 'o', 'color', 'g');
%     text(w(p), dbH(p), ['(',num2str(w(p)),',',num2str(dbH(p)),')'],'color','b');
%     p = find(abs(w - ws / pi) < eps);
%     plot(w(p), dbH(p), 'o', 'color', 'g');
%     text(w(p), dbH(p), ['(',num2str(w(p)),',',num2str(dbH(p)),')'],'color','b');
    ylabel('dB'); xlabel('频率(pi)'); axis([0, 1, -50, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    
    subplot(2, 2, 4); zplane(bd, ad);
    axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, 'LPF', 'svg');
end
function [n, Omgc, bd, ad] = HPF(wp, Rp, ws, As, Fs)
%     wp = 0.25 * pi;                                % 滤波器的通带截止频率
%     ws = 0.4 * pi;                                 % 滤波器的阻带截止频率
%     Rp = 1; As = 15;                               % 滤波器的通阻带衰减指标
    ripple = 10^(-Rp / 20);                        % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);                          % 滤波器的阻带衰减对应的幅度值
                                                   % 转换为模拟滤波器的技术指标
    % 使用修正值
    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
    aal = real(poly(p0));                          % 求原型滤波器的系数a
    [ba, aa] = lp2hp(bal, aal, Omgc);              % 变换为模拟低通滤波器
    %也可将以上4行替换为[ba,aa]=butter(n,Omgc,');直接求模拟滤波器系数
                                                   % 用双线性变换法计算数字滤波器系数
    [bd, ad] = bilinear(ba, aa, Fs);               % 求数字系统的频率特性
    [H, w] = freqz(bd, ad);
    dbH = 20 * log10((abs(H) + eps) / max(abs(H)));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, abs(H));
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, ws / pi, wp / pi, 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, ws / pi, wp / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
    ylabel('dB'); xlabel('频率(pi)'); axis([0, 1, -50, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, ws / pi, wp / pi, 1]);
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    subplot(2, 2, 4); zplane(bd, ad);
    axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, 'HPF', 'svg');
end
function [N, Wc, bd, ad] = BPF(wp, Rp, ws, As, Fs)
    ripple = 10^(-Rp / 20);            % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);              % 滤波器的阻带衰减对应的幅度值
                                       % 转换为模拟滤波器的技术指标
    % 模拟原型滤波器计算
    % 模拟带通性能指标转换成模拟低通性能指标
    WP = wp / pi;
    WS = ws / pi;
    
    % 不使用修正值，使用上面归一化值
%     T = 1 / Fs;
%     WP = (2 / T) * tan(wp / 2);                  % 原型通带频率的预修正
%     WS = (2 / T) * tan(ws / 2);                  % 原型阻带频率的预修正
    % 模拟低通滤波器的构造
    [N, Wc] = buttord(WP, WS, Rp, As);   % 求阶数和边缘频率
    [bd, ad] = butter(N, Wc);            % 求极点、零点和增益
    
%     [N, Wc] = cheb2ord(WP, WS, Rp, As);   % 求阶数和边缘频率
%     [bd, ad] = cheby2(N, Rp, Wc);            % 求极点、零点和增益
    
%     [N, Wc] = cheb1ord(WP, WS, Rp, As);   % 求阶数和边缘频率
%     [bd, ad] = cheby1(N, Rp, Wc);            % 求极点、零点和增益
    
    [H, w] = freqz(bd, ad);
    m = abs(H);
    p = angle(H);
    dbH = 20 * log10((m + eps) / max(m));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, m);
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [0, Attn, ripple, 1]); grid;
    
    subplot(2, 2, 2); plot(w / pi, p / pi);
    ylabel('\phi'); title('相位响应'); axis([0, 1, -1, 1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
    ylabel('dB'); xlabel('频率(pi)'); axis([0, 1, -50, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    
    subplot(2, 2, 4); zplane(bd, ad);
    axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, 'BPF', 'svg');
end
function [N, Wc, bd, ad] = BSF(wp, Rp, ws, As, Fs)
    ripple = 10^(-Rp / 20);            % 滤波器的通带衰减对应的幅度值
    Attn = 10^(-As / 20);              % 滤波器的阻带衰减对应的幅度值
                                       % 转换为模拟滤波器的技术指标
    
    % 模拟原型滤波器计算
    % 模拟带阻性能指标转换成模拟低通性能指标
    WP = wp / pi;
    WS = ws / pi;
    
    % 不使用修正值，使用上面归一化值
%     T = 1 / Fs;
%     WP = (2 / T) * tan(wp / 2);                  % 原型通带频率的预修正
%     WS = (2 / T) * tan(ws / 2);                  % 原型阻带频率的预修正
    % 模拟低通滤波器的构造
    [N, Wc] = buttord(WP, WS, Rp, As);  % 求阶数和边缘频率
    [bd, ad] = butter(N, Wc, 'stop');   % 求极点、零点和增益
    
%     [N, Wc] = ellipord(WP, WS, Rp, As);   % 求阶数和边缘频率
%     [bd, ad] = butter(N, Rp, Wc, 'stop');            % 求极点、零点和增益
    
%     [N, Wc] = cheb2ord(WP, WS, Rp, As);   % 求阶数和边缘频率
%     [bd, ad] = cheby2(N, Rp, Wc, 'stop');            % 求极点、零点和增益
    
%     [N, Wc] = cheb1ord(WP, WS, Rp, As);  % 求阶数和边缘频率
%     [bd, ad] = cheby1(N, Rp, Wc, 'stop');   % 求极点、零点和增益
    
    [H, w] = freqz(bd, ad);
    m = abs(H);
    p = angle(H);
    dbH = 20 * log10((m + eps) / max(m));
    
    h = figure;
    subplot(2, 2, 1); plot(w / pi, m);
    ylabel('|H|'); title('幅度响应'); axis([0, 1, 0, 1.1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [0, Attn, ripple, 1]); grid;
    
    subplot(2, 2, 2); plot(w / pi, p / pi);
    ylabel('\phi'); title('相位响应'); axis([0, 1, -1, 1]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [-1, 0, 1]); grid;
    
    subplot(2, 2, 3); plot(w / pi, dbH); title('幅度响应(dB)');
    
    ylabel('dB'); xlabel('频率(pi)'); axis([0, 1, -50, 5]);
    set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([WP, WS]), 1]);  % [0, 0.15, 0.2, 0.6, 0.65, 1]
    set(gca, 'YTickMode', 'manual', 'YTick', [-50, -As, -Rp, 0]); grid;
    
    subplot(2, 2, 4); zplane(bd, ad);
    axis([-1.1, 1.1, -1.1, 1.1]); title('零极点图');
    
    saveas(h, 'BSF', 'svg');
end