实验七：用 MATLAB 设计 FIR 数字滤波器

1、阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。
% % 实验示例，实验内容未要求这个，所以注释掉了
% % 各种窗函数特性的比较
% N = 64;
% beta = 7.865;
% n = 1: N;
% wbo = boxcar(N); % 矩形 1.8
% wtr = triang(N); % 三角形 6.1
% whn = hanning(N); % 汉宁 6.2
% whm = hamming(N); % 哈明 6.6
% wbl = blackman(N); % 布莱克曼 11
% wka = kaiser(N, beta); % 凯塞 10
% h = figure;
% plot(n, wbo, '-', n, wtr, '*', n, whm, '+', n, whm, '.', n, wbl, 'o', n, wka, 'd');
% axis([0, N, 0, 1.1]);
% legend('矩形', '三角形', '汉宁', '哈明', '布莱克曼', '凯塞');
% saveas(h, '各种窗函数特性的比较', 'svg');
% % 实验示例，实验内容未要求这个，所以注释掉了
% wp = 0.3 * pi;
% ws = 0.45 * pi;
% ap = 0.05;
% as = 50;
%
% % 运行时，将 FIR 函数里面改成 hamming ！！！默认 blackman
% % [N, wc, Rp, As] = FIR(wp, ws, ap, as, 0);
% [N, wc, Rp, As] = FIR(wp, ws, ap, as, 1);
% N, wc, Rp, As
N = 45
wc = 0.3750
Rp = 0.0428
As = 0
2、选择合适的窗函数设计 FIR 数字低通滤波器，要求：ωp=0.24 π ，Rp=0.1dB；ωs=0.3 π ，As=60dB。描绘该滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。
wp = 0.3 * pi;
ws = 0.5 * pi;
ap = 0.1;
as = 60;
%% 运行时，将 FIR 函数里面改成 blackman ！！！默认 blackman，不用改。
% [N, wc, Rp, As] = FIR(wp, ws, ap, as, 0);
[N, wc, Rp, As] = FIR(wp, ws, ap, as, 1);
N, wc, Rp, As
N = 55
wc = 0.4000
Rp = 0.0040
As = 0
3、调用信号产生函数 xtg 产生具有加性噪声的信号 x(t)，并显示信号及其频谱。
% h = figure;
[xt, t] = xtg(1000); % xt
4、采用实验内容步骤 2 中设计的 FIR 数字低通滤波器，调用 Matlab 快速卷积函数 fftfilt 实现对 x(t)的滤波，从高频噪声中提取 x(t)中的单频调幅信号。绘图显示滤波器的频率响应特性曲线、滤波器输出信号的幅频特性图和时域波形图。
wp = 0.24 * pi;
ws = 0.3 * pi;
deltaw = ws - wp;
N0 = ceil(11 * pi / deltaw); % 6.1 6.2 6.6 11
N = N0 + mod(N0 + 1, 2)
N = 185
use_fir1 = 0
use_fir1 = 0
if use_fir1 == 0
% 方法 1: 不用 fir1
windows = blackman(N)'; % triang hanning hamming blackman
wc = (ws + wp) / 2; % 截止频率取归一化通阻带频率的平均值
hd = ideal_lp(wc, N); % 不用 fir1 函数
b = hd .* windows;
else
% 方法 2: 用 fir1
windows = blackman(N); % triang hanning hamming blackman
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
b = fir1(N - 1, wc, windows); % 用 fir1 函数求系统函数系数，windows可省略
end
% % 以下不用求
% [H, w] = freqz(b, 1, 1000, 'whole');
% H = H(1:501)';
% w = w(1:501)';
%
% mag = abs(H);
% db = 20 * log10((mag + eps) / max(mag));
% pha = angle(H);
%
% n = 0:N - 1;
% dw = 2 * pi / 1000;
% Rp = -(min(db(1:floor(wp / dw) + 1))); % 检验通带波动
% As = -round(max(db(floor(wp / dw) + 1:501))); % 检验最小阻带衰减
yt = fftfilt(b, xt);
Hyk = abs(fft(yt));
h = figure;
subplot(2, 1, 1);
plot(t, yt);
% plot((0:999) / 1000, yt);
subplot(2, 1, 2);
stem(Hyk);
axis([80, 120, min(Hyk), max(Hyk)]);
saveas(h, '第4问-幅频特性图和时域波形图', 'svg');
5、选做题，读取音频信号 motherland.wav，得到 xn；
（1）对 xn 进行 I=2 的整数倍 0 值内插，得到音频信号 yn1；
（2）设计一个镜像低通滤波器（可在实验内容 2 的代码上进行修改）；
（3）对 yn1 进行滤波，得到音频信号 yn2。
%% ① 音频播放：依次播放原音频信号 xn、yn1 和 yn2，体验整数倍 0 值内插后的音质。
play = 0; % 是否播放音频
[xn, fs] = audioread('motherland.wav');
length(xn);
if play == 1
sound(xn, fs);
pause(length(xn) / fs);
end
I = 2;
yn1 = zeros(83200*I+1, 1);
for i = 1:length(xn)
yn1(I * i - 1) = xn(i);
% yn1(I * i) = 0;
end
if play == 1
sound(yn1, I * fs);
pause(length(xn) / (I * fs));
end
%% 设计镜像低通滤波器
wp = 0.45 * pi;
ws = 0.5 * pi;
deltaw = ws - wp;
N0 = ceil(6.6 * pi / deltaw);
N = N0 + mod(N0 + 1, 2);
windows = hamming(N);
wc = (ws + wp) / 2 / pi;
b = fir1(N - 1, wc, windows);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:floor(wp / dw) + 1)));
As = -round(max(db(floor(wp / dw) + 1:501)));
N, Rp, As
N = 133
Rp = 0
As = -54
%% ② 取原音频某段信号，如 n=8000~8199。画出该段信号模拟域幅度谱（横坐标为 f Hz）；
% 画出该段信号 I=2 内插后的模拟域幅度谱；画出该段信号内插后再经过镜像滤波后的模拟域幅度谱。
h = figure;
subplot(3, 1, 1);
N = 2048;
Xn = 1 / fs * fft(xn(8000:8199), N); % 从xn中取200点做谱分析，N可取2048
plot((0:N / 2 - 1) * fs / N, abs(Xn(1:N / 2))); % 模拟域幅度谱
xlabel('f / Hz');
ylabel('幅度');
title('原信号的模拟域幅度谱');
subplot(3, 1, 2);
Yn1 = 1 / (I * fs) * fft(yn1(16000:16399), N); % 内插后，200点长变成了400点长
plot((0:N / 2 - 1) * I * fs / N, abs(Yn1(1:N / 2)));
xlabel('f / Hz');
ylabel('幅度');
title('原信号内插后的模拟域幅度谱');
subplot(3, 1, 3);
yn2 = filter(b, 1, yn1); % 对yn1进行滤波，b为所设计的镜像滤波器
Yn2 = 1 / (I * fs) * fft(yn2(16000:16399), N); % 内插后，200点长变成了400点长
plot((0:N / 2 - 1) * I * fs / N, abs(Yn2(1:N / 2)));
xlabel('f / Hz');
ylabel('幅度');
title('原信号先滤波后内插的模拟域幅度谱');
if play == 1
sound(yn2, I*fs);
pause(length(xn) / (I * fs));
end
saveas(h, '第5问-音频信号滤波', 'svg');
思考题
根据下面指标要求设计四种不同类型的 FIR 线性相位数字滤波器，要求画出 h(n)，幅频特性曲线、幅频衰减特性曲线，相频特性曲线，标注相关信息，包括横坐标，纵坐标的单位，曲线名称。
（1） 设计数字低通滤波器，指标为：通带截止频率 𝜔𝑝 = 0.2𝜋，阻带截止频率 𝜔𝑠 = 0.3𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 24𝑑𝐵。（-25 三角 bartlett）
（2） 设计数字高通滤波器，指标为：阻带截止频率 𝜔𝑠 = 0.4𝜋，通带截止频率 𝜔𝑠 = 0.6𝜋，通带衰减 𝛼𝑝 = 0.2𝑑𝐵，阻带衰减 𝛼𝑠 = 43𝑑𝐵。（-44 汉宁 hanning）
（3） 设计数字带通滤波器，指标为：通带范围 0.2𝜋 ≤ 𝜔 ≤ 0.6𝜋，阻带范围 0 ≤ 𝜔 ≤ 0.15𝜋 和 0.65𝜋 ≤ 𝜔 ≤ 𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 50𝑑𝐵。（-53 哈明 hamming）
（4） 设计数字带阻滤波器，指标为：阻带范围 0.2𝜋 ≤ 𝜔 ≤ 0.6𝜋，通带范围 0 ≤ 𝜔 ≤ 0.15𝜋 和 0.65𝜋 ≤ 𝜔 ≤ 𝜋，通带衰减 𝛼𝑝 = 1𝑑𝐵，阻带衰减 𝛼𝑠 = 45𝑑𝐵。（-53 哈明 hamming）
%% 低通
wp = 0.2 * pi;
ws = 0.3 * pi;
ap = 1;
as = 24;
%% 默认 triang, 6.1 已经满足要求
[N, wc, Rp, As] = LPF(wp, ws, ap, as, 0);
% [N, wc, Rp, As] = LPF(wp, ws, ap, as, 1); % 效果不好
N, wc, Rp, As
N = 61
wc = 0.2500
Rp = 0.1105
As = 40
%% 高通
wp = 0.4 * pi;
ws = 0.6 * pi;
ap = 0.2;
as = 43;
%% 默认 hanning, 6.2 已经满足要求
[N, wc, Rp, As] = HPF(wp, ws, ap, as, 0);
% [N, wc, Rp, As] = HPF(wp, ws, ap, as, 1); % 低通
N, wc, Rp, As
N = 31
wc = 0.5000
Rp = 95.5778
As = 0
%% 带通
ws1 = 0.15 * pi;
wp1 = 0.2 * pi;
wp2 = 0.6 * pi;
ws2 = 0.65 * pi;
ap = 1;
as = 50;
%% 默认 hamming, 6.6 已经满足要求
[N, wcl, wch, Rp, As] = BPF([wp1, wp2], [ws1, ws2], ap, as);
N, wcl, wch, Rp, As
N = 133
wcl = 0.5498
wch = 1.9635
Rp = 148.8869
As = 0
%% 带阻
wp1 = 0.15 * pi;
ws1 = 0.2 * pi;
ws2 = 0.6 * pi;
wp2 = 0.65 * pi;
ap = 1;
as = 45;
%% 默认 hamming, 6.6 已经满足要求
[N, wcl, wch, Rp, As] = BSF([wp1, wp2], [ws1, ws2], ap, as);
N, wcl, wch, Rp, As
N = 133
wcl = 0.5498
wch = 1.9635
Rp = 0.0432
As = 0
function hd = ideal_lp(wc, N) % 点0到N-1之间的理想脉冲响应
% wc =截止频率(弧度)
% N = 理想滤波器的长度
% 通带在 ωc1～ωc2：
% Hd = ideal_lp(wc2,N) - ideal_lp(wc1,N)
% 阻带在 ωc1～ωc2：
% Hd = ideal_lp(pi,N) - ( ideal_lp(wc2,N) - ideal_lp(wc1,N) )
tao = (N - 1) / 2;
n = 0:(N - 1);
m = n - tao + eps; % 加一个小数以避免0作除数
hd = sin(wc * m) ./ (pi * m);
end
function [N, wc, Rp, As] = FIR(wp, ws, ap, as, use_fir1)
%% triang -> 6.1
%% hanning -> 6.2
%% hamming -> 6.6
%% blackman -> 11
%% kaiser -> 略，自己查表8.2.2
deltaw = ws - wp;
% N0 = ceil(6.6 * pi / deltaw); % hamming
% N0 = ceil(11 * pi / deltaw); % blackman
N0 = ceil(11 * pi / deltaw); % 修改1
N = N0 + mod(N0 + 1, 2); % 为实现FIR类型1偶对称滤波器，应确保N为奇数
if use_fir1 == 0
% 方法 1: 不用 fir1
% windows = hamming(N)'; % hamming
% windows = blackman(N)'; % blackman
windows = blackman(N)'; % 修改2
wc = (ws + wp) / 2; % 截止频率取归一化通阻带频率的平均值
% 使用 ideal_lp 除 2
hd = ideal_lp(wc, N); % 不用 fir1 函数
b = hd .* windows;
else
% 方法 2: 用 fir1
% windows = hamming(N); % hamming
% windows = blackman(N); % blackman
windows = blackman(N); % 修改3
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
% 使用 fir 除 2 pi
b = fir1(N - 1, wc, windows); % 用 fir1 函数求系统函数系数，windows可省略
end
[H, w] = freqz(b, 1, 1000, 'whole');
H = H(1:501)';
w = w(1:501)';
mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));
pha = angle(H);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:floor(wp / dw) + 1))); % 检验通带波动
As = -round(max(db(floor(wp / dw) + 1:501))); % 检验最小阻带衰减
h = figure;
subplot(2, 2, 1);
stem(n, b);
axis([0, N, 1.1 * min(b), 1.1 * max(b)]);
title('实际脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2, 2, 2);
stem(n, windows);
axis([0, N, 0, 1.1]);
title('窗函数特性');
xlabel('n');
ylabel('wd(n)');
subplot(2, 2, 3);
plot(w / pi, db);
axis([0, 1, -100, 10]);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-80, -as, -ap, 0]); grid;
subplot(2, 2, 4);
plot(w / pi, pha);
axis([0, 1, -4, 4]);
title('相位频率响应');
xlabel('频率（单位：\pi）');
ylabel('\phi(\omega)');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'XTick', [-3.1416, 0, 3.1416, 4]); grid;
saveas(h, '第2问-FIR数字低通滤波器-频率响应曲线', 'svg');
end
function [N, wc, Rp, As] = LPF(wp, ws, ap, as, use_win)
deltaw = ws - wp;
N0 = ceil(6.1 * pi / deltaw); % 修改1
N = N0 + mod(N0 + 1, 2); % 为实现FIR类型1偶对称滤波器，应确保N为奇数
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
if use_win == 0
b = fir1(N - 1, wc, 'low'); % 用fir1 函数求系统函数系数
else
windows = triang(N); % 修改2
b = fir1(N - 1, wc, windows); % 用fir1 函数求系统函数系数
end
[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';
mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));
pha = angle(H);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:wp / dw + 1))); % 检验通带波动
As = -round(max(db(ws / dw + 1:501))); % 检验最小阻带衰减
h = figure;
subplot(2, 2, 1);
stem(n, b);
axis([0, N, 1.1 * min(b), 1.1 * max(b)]);
title('实际脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2, 2, 2);
plot(w / pi, db);
title('幅频衰减特性曲线');
xlabel('频率（单位：\pi）');
ylabel('幅度（dB）');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 3);
plot(w / pi, mag);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 4);
plot(w / pi, pha);
axis([0, 1, -4, 4]);
title('相位频率响应');
xlabel('频率（单位：\pi）');
ylabel('\phi(\omega)');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]); grid;
saveas(h, 'FIR-LPF-频率响应曲线', 'svg');
end
function [N, wc, Rp, As] = HPF(wp, ws, ap, as, use_win)
deltaw = ws - wp;
N0 = ceil(6.2 * pi / deltaw); % 修改1
N = N0 + mod(N0 + 1, 2); % 为实现FIR类型1偶对称滤波器，应确保N为奇数
wc = (ws + wp) / 2 / pi; % 截止频率取归一化通阻带频率的平均值
if use_win == 0
b = fir1(N - 1, wc, 'high'); % 用fir1 函数求系统函数系数
else % 不行，这个还是低通！！！
windows = hanning(N); % 修改2
b = fir1(N - 1, wc, windows); % 用fir1 函数求系统函数系数
end
[H, w] = freqz(b, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';
mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));
pha = angle(H);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:wp / dw + 1))); % 检验通带波动
As = -round(max(db(ws / dw + 1:501))); % 检验最小阻带衰减
h = figure;
subplot(2, 2, 1);
stem(n, b);
axis([0, N, 1.1 * min(b), 1.1 * max(b)]);
title('实际脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2, 2, 2);
plot(w / pi, db);
title('幅频衰减特性曲线');
xlabel('频率（单位：\pi）');
ylabel('幅度（dB）');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 3);
plot(w / pi, mag);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 4);
plot(w / pi, pha);
axis([0, 1, -4, 4]);
title('相位频率响应');
xlabel('频率（单位：\pi）');
ylabel('\phi(\omega)');
set(gca, 'XTickMode', 'manual', 'XTick', [0, wp / pi, ws / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]); grid;
saveas(h, 'FIR-HPF-频率响应曲线', 'svg');
end
function [N, wcl, wch, Rp, As] = BPF(wp, ws, ap, as)
deltaw = min(abs(wp - ws));
N0 = ceil(6.6 * pi / deltaw); % 修改1
N = N0 + mod(N0 + 1, 2); % 为实现FIR类型1偶对称滤波器，应确保N为奇数
wcl = (ws(1) + wp(1)) / 2; % 截止频率取归一化通阻带频率的平均值
wch = (ws(2) + wp(2)) / 2; % 截止频率取归一化通阻带频率的平均值
hd = ideal_lp(wch, N) - ideal_lp(wcl, N);
b = hamming(N)'; % 修改2
h = hd .* b;
[H, w] = freqz(h, 1, 1000, 'whole');
H = H(1:501)';
w = w(1:501)';
mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));
pha = angle(H);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:wp / dw + 1))); % 检验通带波动
As = -round(max(db(ws / dw + 1:501))); % 检验最小阻带衰减
h = figure;
subplot(2, 2, 1);
stem(n, b);
axis([0, N, 1.1 * min(b), 1.1 * max(b)]);
title('实际脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2, 2, 2);
plot(w / pi, db);
title('幅频衰减特性曲线');
xlabel('频率（单位：\pi）');
ylabel('幅度（dB）');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 3);
plot(w / pi, mag);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 4);
plot(w / pi, pha);
axis([0, 1, -4, 4]);
title('相位频率响应');
xlabel('频率（单位：\pi）');
ylabel('\phi(\omega)');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]); grid;
saveas(h, 'FIR-BPF-频率响应曲线', 'svg');
end
function [N, wcl, wch, Rp, As] = BSF(wp, ws, ap, as)
deltaw = min(abs(wp - ws));
N0 = ceil(6.6 * pi / deltaw); % 修改1
N = N0 + mod(N0 + 1, 2); % 为实现FIR类型1偶对称滤波器，应确保N为奇数
wcl = (ws(1) + wp(1)) / 2; % 截止频率取归一化通阻带频率的平均值
wch = (ws(2) + wp(2)) / 2; % 截止频率取归一化通阻带频率的平均值
hd = ideal_lp(pi, N) - ideal_lp(wch, N) + ideal_lp(wcl, N);
b = hamming(N)'; % 修改2
h = hd .* b;
[H, w] = freqz(h, 1, 1000, 'whole');
H = (H(1:501))';
w = (w(1:501))';
mag = abs(H);
db = 20 * log10((mag + eps) / max(mag));
pha = angle(H);
n = 0:N - 1;
dw = 2 * pi / 1000;
Rp = -(min(db(1:wp / dw + 1))); % 检验通带波动
As = -round(max(db(ws / dw + 1:501))); % 检验最小阻带衰减
h = figure;
subplot(2, 2, 1);
stem(n, b);
axis([0, N, 1.1 * min(b), 1.1 * max(b)]);
title('实际脉冲响应');
xlabel('n');
ylabel('h(n)');
subplot(2, 2, 2);
plot(w / pi, db);
title('幅频衰减特性曲线');
xlabel('频率（单位：\pi）');
ylabel('幅度（dB）');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 3);
plot(w / pi, mag);
title('幅度频率响应');
xlabel('频率（单位：\pi）');
ylabel('H(e^{j\omega})');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-as-10, -as, -ap, 0]); grid;
subplot(2, 2, 4);
plot(w / pi, pha);
axis([0, 1, -4, 4]);
title('相位频率响应');
xlabel('频率（单位：\pi）');
ylabel('\phi(\omega)');
set(gca, 'XTickMode', 'manual', 'XTick', [0, sort([wp, ws]) / pi, 1]);
set(gca, 'YTickMode', 'manual', 'YTick', [-3.1416, 0, 3.1416, 4]); grid;
saveas(h, 'FIR-BSF-频率响应曲线', 'svg');
end