西电计科院【嵌入式】数字信号处理实验

ZhiHL

已于 2024-08-15 20:26:05 修改

阅读量1k

点赞数 18

分类专栏：课程学习文章标签：信号处理学习

于 2023-12-12 23:09:36 首次发布

本文链接：https://blog.csdn.net/weixin_63554791/article/details/134960854

版权

课程学习专栏收录该内容

17 篇文章 0 订阅

订阅专栏

博主是21级计科院的，数字信号处理实验 97 分，笔者写这篇文章的动因在于当时自己写实验代码时，难以找到比较优秀的参考博客，很是疑惑之前的学长怎么没有流传下来。在此将自己当时的实验项目分享如下，希望各位学弟学妹以此为基础争取更优秀地完成实验！
复习资源以及详细内容见博客 DigitalSignalProcessing

文章目录

一、实验要求

利用傅立叶级数展开的方法，自由生成所需的 $x (t)$ ；
通过选择不同的采样间隔 $T$ （分别选 $T>\frac{1}{2f_c}$ 或 $T<\frac{1}{2fc}$ ），从 $x (t)$ 获得相应的 $x (n)$ （作出 $x (n)$ 图形）；
对获得的不同 $x (n)$ 分别作傅立叶变换，分析其频率响应特性（给出幅频与相频特性曲线）；
利用巴特沃思、切比雪夫或椭圆滤波器设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同 $x (n)$ 分别进行滤波（画出滤波结果），然后加以讨论；
利用窗函数设计法或频率采样法设计数字滤波器（滤波特性自定），要求通过改变滤波器参数或特性（低通、高通、带通或带阻等）设计至少两种数字滤波器，分析所设计滤波器（画出频率特性曲线），并对上述给出的不同 $x (n)$ 分别进行滤波（画出滤波结果），然后加以讨论。

二、源代码

I. 项目列表

项目列表

II. 各部分源代码

（0）项目主函数配置

%% Project Introduction:
%   This project is developed to design some signal filters based on digital
% signal processing.
clear, close all;
%% Produce and sample digital signal
f1 = 10;
f2 = 20;
f3 = 30;          % so the fc = f3 = 30Hz
Np = 200;          % number of periods for time-domain window
%% Experiment 1 (Choosing samling frequency fs = 3fc (fs > 2fs))
fs1 = 90;          % sampling frequency
xn1 = ProduceSamplingSignal(f1, f2, f3, fs1, Np, 'Sampling Analog Signal(fs = 3fc)');
DFTAnalysis(xn1, fs1, 'Frequency Response Characteristics(fs = 3fc)');

%% Experiment 2 (Choosing samling frequency fs = 2fc)
fs2 = 60;          % sampling frequency
xn2 = ProduceSamplingSignal(f1, f2, f3, fs2, Np, 'Sampling Analog Signal(fs = 2fc)');
DFTAnalysis(xn2, fs2, 'Frequency Response Characteristics(fs = 2fc)');

%% Experiment 3 (Choosing samling frequency fs < 2fc)
fs3 = 40;          % sampling frequency
xn3 = ProduceSamplingSignal(f1, f2, f3, fs3, Np, 'Sampling Analog Signal(fs < 2fc)');
DFTAnalysis(xn3, fs3, 'Frequency Response Characteristics(fs < 2fc)');

%% Experiment Description
% Experiment 4-7: Design a digital filter respectively with band pass, high
% pass, low pass, band stop based on ellipord.
%% Experiment 4: Design a digital filter with band pass using ellipord
fpl = 15; fpu=25; fsl=13; fsu=28;
rp = 1;             
rs = 40;
ellipBandPass(fpl, fpu, fsl, fsu, rp, rs, xn1, fs1, f1, Np, 'Digital Filter With Band Pass Using Ellipord(fs = 3fc)');
ellipBandPass(fpl, fpu, fsl, fsu, rp, rs, xn2, fs2, f1, Np, 'Digital Filter With Band Pass Using Ellipord(fs = 2fc)');
ellipBandPass(8, 10, 6, 12, rp, rs, xn3, fs3, f1, Np, 'Digital Filter With Band Pass Using Ellipord(fs < 2fc)');

%% Experiment 5: Design a digital filter with high pass using ellipord
fpz = 16; fsz = 13;
rp = 1;             
rs = 40;
ellipHighPass(fpz, fsz, rp, rs, xn1, fs1, f1, Np, 'Digital Filter With High Pass Using Ellipord(fs = 3fc)');
ellipHighPass(fpz, fsz, rp, rs, xn2, fs2, f1, Np, 'Digital Filter With High Pass Using Ellipord(fs = 2fc)');
ellipHighPass(15, 12, rp, rs, xn3, fs3, f1, Np, 'Digital Filter With High Pass Using Ellipord(fs < 2fc)');

%% Experiment 6: Design a digital filter with low pass using ellipord
fpz = 23; fsz=28; 
rp = 1;             
rs = 40;
ellipLowPass(fpz, fsz, rp, rs, xn1, fs1, f1, Np, 'Digital Filter With Low Pass Using Ellipord(fs = 3fc)');
ellipLowPass(fpz, fsz, rp, rs, xn2, fs2, f1, Np, 'Digital Filter With Low Pass Using Ellipord(fs = 2fc)');
ellipLowPass(12, 15, rp, rs, xn3, fs3, f1, Np, 'Digital Filter With Low Pass Using Ellipord(fs < 2fc)');

%% Experiment 7: Design a digital filter with band stop using ellipord
fpl = 15; fpu=25; fsl=17; fsu=22;
rp = 1;             
rs = 40;
ellipBandStop(fpl, fpu, fsl, fsu, rp, rs, xn1, fs1, f1, Np, 'Digital Filter With Band Stop Using Ellipord(fs = 3fc)');
ellipBandStop(fpl, fpu, fsl, fsu, rp, rs, xn2, fs2, f1, Np, 'Digital Filter With Band Stop Using Ellipord(fs = 2fc)');
ellipBandStop(5, 17, 8, 12, rp, rs, xn3, fs3, f1, Np, 'Digital Filter With Band Stop Using Ellipord(fs < 2fc)');

%% Experiment Description
% Experiment 8-11: Design a digital filter respectively with high pass, low
% pass, band pass, band stop based on hamming window.
%% Experiment 8: Design a digital filter with high pass using hamming window
fpz = 16; fsz = 13;
firlHighPass(fpz, fsz, xn1, fs1, f1, Np, 'Digital Filter With High Pass Using Hamming Window(fs = 3fc)');
firlHighPass(fpz, fsz, xn2, fs2, f1, Np, 'Digital Filter With High Pass Using Hamming Window(fs = 2fc)');
firlHighPass(15, 12, xn3, fs3, f1, Np, 'Digital Filter With High Pass Using Hamming Window(fs < 2fc)');
%% Experiment 9: Design a digital filter with low pass using hamming window
fpz = 23; fsz = 28;
firlLowPass(fpz, fsz, xn1, fs1, f1, Np, 'Digital Filter With Low Pass Using Hamming Window(fs = 3fc)');
firlLowPass(fpz, fsz, xn2, fs2, f1, Np, 'Digital Filter With Low Pass Using Hamming Window(fs = 2fc)');
firlLowPass(13, 17, xn3, fs3, f1, Np, 'Digital Filter With Low Pass Using Hamming Window(fs < 2fc)');
%% Experiment 10: Design a digital filter with band pass using hamming window
fpl = 15; fpu = 25;
firlBandPass(fpl, fpu, xn1, fs1, f1, Np, 'Digital Filter With Band Pass Using Hamming Window(fs = 3fc)');
firlBandPass(fpl, fpu, xn2, fs2, f1, Np, 'Digital Filter With Band Pass Using Hamming Window(fs = 2fc)');
firlBandPass(7, 15, xn3, fs3, f1, Np, 'Digital Filter With Band Pass Using Hamming Window(fs < 2fc)');
%% Experiment 11: Design a digital filter with band stop using hamming window
fsl = 15; fsu = 25;
firlBandStop(fsl, fsu, xn1, fs1, f1, Np, 'Digital Filter With Band Stop Using Hamming Window(fs = 3fc)');
firlBandStop(fsl, fsu, xn2, fs2, f1, Np, 'Digital Filter With Band Stop Using Hamming Window(fs = 2fc)');
firlBandStop(7, 15, xn3, fs3, f1, Np, 'Digital Filter With Band Band Stop Hamming Window(fs < 2fc)');

（1）利用傅里叶级数展开法生成 $x (t)$ ，采样生成 $x (n)$

function xn = ProduceSamplingSignal(f1, f2, f3, fs, Np, Alltitle)
% Function Description: 
%        We want to make a digital signal composed of three frequency
%        components and sample the produced signal.
% Inputs: 
%        f1, f2, f3: means our selected frequency components, fs
%            represents the sampling frequency.
%        Np: means the number of periods.
% Outputs:
%        xn: represents the sampled signal.

    period = 1/f1;         % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;           % sampling timestep
    t = 0: Ts : T;         % samping sequence of discrete sampling points
    t0 = 0: 0.0001: T;     % analog time sequence

    % Step I: Produce digital signal
    xt = cos(2*pi*f1*t0) + cos(2*pi*f2*t0) + cos(2*pi*f3*t0);
    % Step II: Sample produced signal
    xn = cos(2*pi*f1*t) + cos(2*pi*f2*t) + cos(2*pi*f3*t);

    % Step III: Visualize produced signal and sampled signal
    figure('Position', [210, 80, 950, 750]);
    subplot(2, 1, 1);
    plot(t0, xt);
    title('Time-domain signal $x(t)$', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('$t/s$', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-2.5, 3.5]);
    grid on

    subplot(2, 1, 2);
    stem(t, xn);
    title('Time-domain sampled signal $x(n)$', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('$t/s$', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-2.5, 3.5]);
    grid on
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

（2）对采样序列进行傅里叶变换并进行频谱分析

function DFTAnalysis(xn, fs, Alltitle)
% Function Description:
%       This function calculates the DFT[x(n)] and do spectral analysis.
% Inputs:
%       xn: digital discrete signal
%       fs: sampling frequency
% Outputs:
%       No return

    N = length(xn);     % number of sampling points
    df = fs / N;        % spectral resolution
    f1 = (0:N-1)*df;    % tranverse to the frequncy sequence
    f2 =  2*(0:N-1)/N;
    fprintf("采样点数: %d\n", N);
    % DFT using FFT algorithm
    Xk = fft(xn, N);   
    % Tranverse to the real amplitude
    RM = 2*abs(Xk)/N;
    Angle = angle(Xk);

    figure('Position', [210, 80, 950, 750]);
    % Amplitude-Frequency Characteristics
    subplot(4,1,1);
    stem(f1, RM,'.');
    title('Amplitude-Frequency Characteristics', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('$f$/Hz', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    
    % Phase-Frequency Characteristics
    subplot(4,1,2);
    stem(f1, Angle,'.'); 
    line([(N-1)*df, 0],[0,0]);
    title('Phase-Frequency Characteristics', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('$f$/Hz', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    % Amplitude-Frequency Characteristics
    subplot(4,1,3);
    plot(f2, abs(Xk));
    title('Amplitude-Frequency Characteristics', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    
    % Phase-Frequency Characteristics
    subplot(4,1,4);
    plot(f2, Angle);
    title('Phase-Frequency Characteristics', 'Interpreter', 'latex', 'FontSize', 12);
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

（3）利用椭圆滤波器设计四种滤波特性的数字滤波器

function ellipBandPass(fpl, fpu, fsl, fsu, rp, rs, x, fs, f1, Np, Alltitle)
    wp = [2*fpl/fs, 2*fpu/fs];
    ws = [2*fsl/fs, 2*fsu/fs];
    [N, wn] = ellipord(wp, ws, rp, rs);    % 获取阶数和截止频率
    [B, A] = ellip(N, rp, rs, wn, 'bandpass');         % 获得转移函数系数

    filter_bp_s = filter(B, A, x);
    X_bp_s = abs(fft(filter_bp_s));
    X_bp_s_angle = angle(fft(filter_bp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    % 带通滤波器频谱特性
    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 512;
    wk = 0:pi/M:pi;
    Hz = freqz(B,A,wk);
    plot(wk/pi, 20*log10(abs(Hz)));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('带通滤波器频谱特性');
    axis([0.2,0.9,-80,20]);set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20);
    grid on;

    subplot(4,4,[3,4,7,8]);
    plot(t, filter_bp_s);
    xlabel('t/s', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_bp_s);
    title('带通滤波后频域幅度特性');
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_bp_s_angle);
    title('带通滤波后频域相位特性');    
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function ellipHighPass(fpz, fsz, rp, rs, x, fs, f1, Np, Alltitle)
    wpz = 2*fpz/fs;
    wsz = 2*fsz/fs;
    [N, wn] = ellipord(wpz, wsz, rp, rs);    % 获取阶数和截止频率
    [B, A] = ellip(N, rp, rs, wn, 'high');         % 获得转移函数系数

    filter_hp_s = filter(B, A, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    % 高通滤波器频谱特性
    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 512;
    wk = 0:pi/M:pi;
    Hz = freqz(B,A,wk);
    plot(wk/pi, 20*log10(abs(Hz)));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('高通滤波器频谱特性');
    axis([0.2,0.8,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20);
    grid on;

    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('高通滤波后时域图形');
    xlabel('t/s', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('高通滤波后频域幅度特性');
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('高通滤波后频域相位特性');    
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function ellipLowPass(fpz, fsz, rp, rs, x, fs, f1, Np, Alltitle)
    wpz = 2*fpz/fs;
    wsz = 2*fsz/fs;
    [N, wn] = ellipord(wpz, wsz, rp, rs);    % 获取阶数和截止频率
    [B, A] = ellip(N, rp, rs, wn, 'low');         % 获得转移函数系数

    filter_hp_s = filter(B, A, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    % 低通滤波器频谱特性
    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 512;
    wk = 0:pi/M:pi;
    Hz = freqz(B,A,wk);
    plot(wk/pi, 20 * log10(abs(Hz)));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('低通滤波器频谱特性');
    axis([0.2,0.9,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20)
    grid on;

    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('低通滤波后时域图形');
    xlabel('t/s', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('低通滤波后频域幅度特性');
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('低通滤波后频域相位特性');   
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function ellipBandStop(fpl, fpu, fsl, fsu, rp, rs, x, fs, f1, Np, Alltitle)
    wp = [2*fpl/fs, 2*fpu/fs];
    ws = [2*fsl/fs, 2*fsu/fs];
    [N, wn] = ellipord(wp, ws, rp, rs);    % 获取阶数和截止频率
    [B, A] = ellip(N, rp, rs, wn, 'stop');         % 获得转移函数系数

    filter_bp_s = filter(B, A, x);
    X_bp_s = abs(fft(filter_bp_s));
    X_bp_s_angle = angle(fft(filter_bp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    % 带阻滤波器频谱特性
    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 512;
    wk = 0:pi/M:pi;
    Hz = freqz(B,A,wk);
    plot(wk/pi, 20*log10(abs(Hz)));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('带阻滤波器频谱特性');
    axis([0.2,0.9,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20);
    grid on;
   
    subplot(4,4,[3,4,7,8]);
    plot(t, filter_bp_s);
    title('带阻滤波后时域图形');
    xlabel('t/s', 'Interpreter', 'latex', 'FontSize', 12);
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_bp_s);
    title('带阻滤波后频域幅度特性');
    ylabel('Amplitude', 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_bp_s_angle);
    title('带阻滤波后频域相位特性');    
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('Phase', 'Interpreter', 'latex', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

（4）利用窗函数设计法设计四种滤波特性的数字滤波器

function firlHighPass(fpz, fsz, x, fs, f1, Np, Alltitle)
    wpz = 2 * pi * fpz / fs;
    wsz = 2 * pi * fsz / fs;
    DB = wpz - wsz;              % 计算过渡带宽度
    N0 = ceil(6.2 * pi / DB);    % 计算所需h(n)长度N0
    N = N0 + mod(N0 + 1, 2);     % 确保h(n)长度N是奇数

    wc = (wpz + wsz) /2 / pi;    % 计算理想高通滤波器通带截止频率
    hn = fir1(N-1, wc, 'high', hamming(N));

    filter_hp_s = filter(hn, 1, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 1024;
    k = 1:M / 2;
    wk = 2*(0:M/2-1)/M;
    Hz = freqz(hn, 1);
    plot(wk, 20*log10(abs(Hz(k))));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('高通滤波器频谱特性')
    axis([0.2,0.8,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20)
    grid on;
   
    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('高通滤波后时域图形');
    txt = xlabel('t/s', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('高通滤波后频域幅度特性');
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('高通滤波后频域相位特性');    
    txt = ylabel('Phase', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    xlabel('\omega/\pi', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function firlLowPass(fpz, fsz, x, fs, f1, Np, Alltitle)
    wpz = 2 * pi * fpz / fs;
    wsz = 2 * pi * fsz / fs;
    DB = wsz - wpz;           
    N0 = ceil(6.2 * pi / DB);
    N = N0 + mod(N0 + 1, 2);     

    wc = (wpz + wsz) / 2 / pi;
    hn = fir1(N-1, wc, 'low', hamming(N));

    filter_hp_s = filter(hn, 1, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 1024;
    k = 1:M / 2;
    wk = 2*(0:M/2-1)/M;
    Hz = freqz(hn, 1);
    plot(wk, 20*log10(abs(Hz(k))));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('低通滤波器频谱特性')
    axis([0.2,0.9,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20)
    grid on;
   
    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('低通滤波后时域图形');
    txt = xlabel('t/s', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('低通滤波后频域幅度特性');
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('低通滤波后频域相位特性');    
    txt = ylabel('Phase', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    xlabel('\omega/\pi', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function firlBandPass(fpl, fpu, x, fs, f1, Np, Alltitle)
    wpl = 2 * fpl / fs;
    wpu = 2 * fpu / fs;
    fpass = [wpl, wpu];
    N = 111;
    hn = fir1(N-1, fpass, 'bandpass', hamming(N));

    filter_hp_s = filter(hn, 1, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 1024;
    k = 1:M / 2;
    wk = 2*(0:M/2-1)/M;
    Hz = freqz(hn, 1);
    plot(wk, 20*log10(abs(Hz(k))));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('带通滤波器频谱特性')
    axis([0.2,0.9,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20)
    grid on;
   
    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('带通滤波后时域图形');
    txt = xlabel('t/s', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('带通滤波后频域幅度特性');
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('带通滤波后频域相位特性');    
    txt = ylabel('Phase', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    xlabel('\omega/\pi', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end

function firlBandStop(fsl, fsu, x, fs, f1, Np, Alltitle)
    wsl = 2 * fsl / fs;
    wsu = 2 * fsu / fs;
    fstop = [wsl, wsu];
    N = 111;
    hn = fir1(N-1, fstop, 'stop', hamming(N));

    filter_hp_s = filter(hn, 1, x);
    X_hp_s = abs(fft(filter_hp_s));
    X_hp_s_angle = angle(fft(filter_hp_s));
    
    % plot the graphs    
    period = 1/f1;        % the period of analog signal(assuming f1 is the minimal)
    T = Np*period;         % sampling time-domain window(several periods)
    Ts = 1 / fs;          % sampling timestep
    t = 0: Ts : T;       % samping sequence of discrete sampling points
    N = length(x);    % number of sampling points
    f =  2*(0:N-1)/N;

    figure('Position', [210, 80, 950, 750]);
    subplot(4,4,[1,2,5,6]);
    M = 1024;
    k = 1:M / 2;
    wk = 2*(0:M/2-1)/M;
    Hz = freqz(hn, 1);
    plot(wk, 20*log10(abs(Hz(k))));
    xlabel('\omega/\pi', 'FontSize', 12);
    ylabel('$20lg|Hg(\omega)|$', 'Interpreter', 'latex', 'FontSize', 12);
    title('带阻滤波器频谱特性')
    axis([0.2,0.9,-80,20]);
    set(gca,'Xtick',0:0.1:1,'Ytick',-80:20:20)
    grid on;
   
    subplot(4,4,[3,4,7,8]);
    plot(t, filter_hp_s);
    title('带阻滤波后时域图形');
    txt = xlabel('t/s', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[9, 10, 11, 12]);
    plot(f, X_hp_s);
    title('带阻滤波后频域幅度特性');
    txt = ylabel('Amplitude', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;

    subplot(4,4,[13, 14, 15, 16]);
    plot(f, X_hp_s_angle);
    title('带阻滤波后频域相位特性');    
    txt = ylabel('Phase', 'FontSize', 12);
    ylim([-3.5, 3.5]);
    xlabel('\omega/\pi', 'FontSize', 12);
    set(txt, 'Interpreter', 'latex', 'FontSize', 12);
    grid on;
    sgtitle(Alltitle, 'FontName', 'Times New Roman', 'FontSize', 14);
end