本文是IIR数字滤波器设计，如果需要了解模拟滤波器或者FIR的内容，可以看我写的另外两篇博客，如下：

1.巴特沃斯模拟滤波器（低通，高通，带通，带阻）设计-MATLAB实现

3.MATLAB实现有限脉冲响应数字滤波器（FIR）

1. 基础知识介绍
我们首先明确一个知识（这个非常重要）：

某正弦信号，频率为50Hz
这意味着 信号的模拟频率 f = 50 (Hz)，注意它的单位是Hz

信号的表达式为

注意！！！
模拟滤波器设计中用的频率是指模拟角频率 Ω 
数字滤波器设计中用的频率是指归一化数字角频率 w 

模拟滤波器设计中用的频率是指模拟角频率 Ω 
数字滤波器设计中用的频率是指归一化数字角频率 w 

模拟滤波器设计中用的频率是指模拟角频率 Ω 
数字滤波器设计中用的频率是指归一化数字角频率 w 

2. 数字滤波器简介
数字滤波器特点：

输入信号中有用的频率成分和希望滤除的频率成分占有不同的频带，我们需要通过一个合适的选频滤波器滤除干扰，得到纯净信号。

3.IIR与FIR介绍
3.1 如何选择
数字滤波器可以分为：无限脉冲响应(IIR)数字滤波器、有限脉冲响应(FIR)数字滤波器两种
我们一般使用IIR即可，因为它设计起来较为方便，当某个滤波器要求线性相位时，我们才使用FIR。

3.2 设计方法
对IIR滤波器，我们一般采用间接法来设计，即先按照参数设计好一个模拟滤波器，然后将其变换为数字滤波器，因为模拟滤波器的理论和设计方法比较成熟。

对FIR滤波器，我们一般采用窗函数设计法。

4.无限脉冲响应(IIR)数字滤波器设计
4.1 函数介绍
（1）buttord - 求解滤波器的阶数N和3dB截止频率wc

[N，wc] = buttord(wp, ws, Rp, As)

该函数用于求解巴特沃斯数字滤波器的阶数N和3dB截止频率wc

输入参数如下：

通带边界数字频率wp、阻带边界模拟频率ws（对pi归一化数字角频率，单位是rad）

通带最大衰减Rp、阻带最小衰减As（单位是dB）

注：因为设计的是数字滤波器时，所以没有’s’这个参数了。

（2）butter - 求解N阶滤波器的具体参数B和A，求解完B和A后滤波器就设计完成了。

[B, A] = butter(N, wc, ‘ftype’)

该函数用于求解N阶巴特沃斯数字滤波器的具体参数B和A，求解完B和A后滤波器就设计完成了。

输入参数如下：
N - 滤波器阶数
wc - 3dB截止数字频率wc （单位rad，N和wc都是用buttord函数计算出来的）

ftype - 滤波器类型‘’：

当输入wc为一维向量时：
设计低通滤波器时ftype不填，设计高通滤波器的话令ftype=high

当输入wc为二维向量[wcl,wcu]时：
设计带通滤波器时ftype不填，设计带阻滤波器的话令ftype=stop

（3）filter - 滤波函数

y = filter(B,A,x)

这个就是滤波函数了，x是输入的有噪声的信号，B，A就是设计好的滤波器参数，得到的输入y就是滤波后的信号了。

4.2 直接调用Matlab函数求解
（1）低通滤波器
例： 设计数字低通滤波器，wp = 0.2 π rad,Rp = 1 dB, ws = 0.3 π  rad, As = 15 dB，采样频率1Hz。

求解这个问题我们需要用到以下两个函数：

[N，wc] = buttord(wp, ws, Rp, As) 
[B, A] = butter(N, wc, ‘ftype’)

注意！！！
函数涉及到的频率都是对 pi 归一化后的数字角频率，单位是rad

所以 wp = 0.2 pi / pi = 0.2 ws = 0.3 pi / pi = 0.3

代码如下

wp = 0.2;
ws = 0.3;
Rp = 1;
As = 15
 
[N, wc] = buttord(wp, ws, Rp, As);
[B, A] = butter(N, wc);

至此，数字低通滤波器设计完成了，如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

接着我们画出数字低通滤波器的幅频特性曲线，代码如下：

%画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);

得到的结果如下所示：

低通数字滤波器

（2）高通滤波器
例： 设计数字低通滤波器，wp = 0.3 π rad,Rp = 1 dB, ws = 0.2 π rad, As = 15 dB.采样频率1Hz。

求解这个问题我们需要用到以下两个函数：

[N，wc] = buttord(wp, ws, Rp, As) 
[B, A] = butter(N, wc, ‘ftype’)

函数涉及到的频率都是对 π 归一化后的数字角频率，单位是rad

所以wp = 0.3 pi / pi = 0.3 ws = 0.2 pi / pi = 0.2

代码如下：

wp = 0.3;
ws = 0.2;
Rp = 1;
As = 15
 
[N, wc] = buttord(wp, ws, Rp, As);
[B, A] = butter(N, wc,'high');

至此，数字高通滤波器设计完成了，
如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

接着我们画出数字高通滤波器的幅频特性曲线，代码如下：

%画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);

结果如下：

在这里插入图片描述

（3）数字带通滤波器

例： 设计数字带通滤波器，采样频率Fs = 8kHz，保留2025 - 2225 Hz的部分，幅度失真小于1dB，滤除1500Hz以下和2700Hz以上部分的噪声，衰减大于40dB。

数字频率

代码如下：

Fs = 8000 %采样频率
fpl = 2025; %通带较小频率
fpu = 2225; %通带较大频率
fsl = 1500; %阻带较小频率
fsu = 2700; %阻带较大频率
 
ws = [2 * fpl / Fs, 2 * fpu / Fs];
wp = [2 * fsl / Fs, 2 * fsu / Fs];
Rp = 1;
As = 40;
 
[N, wc] = buttord(wp, ws, Rp, As);%此时输入wp和ws都是二维的，输出wc也是两维的
[B, A] = butter(N, wc);

至此，数字带通滤波器设计完成了，如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

接着我们画出数字带通滤波器的幅频特性曲线，代码如下：

%画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);

画图结果如下：

在这里插入图片描述

（4）数字带阻滤波器

例： 设计数字带阻滤波器，采样频率Fs = 8kHz，保留1500Hz以下和2700Hz以上部分的部分，幅度失真小于1dB，滤除2025 - 2225 Hz的噪声，衰减大于40dB。

代码如下：

Fs = 8000; %采样频率
fpl = 1500; %通带较小频率
fpu = 2700; %通带较大频率
fsl = 2025; %阻带较小频率
fsu = 2225; %阻带较大频率
 
wp = [2 * fpl / Fs, 2 * fpu / Fs];
ws = [2 * fsl / Fs, 2 * fsu / Fs];
Rp = 1;
As = 40;
 
[N, wc] = buttord(wp, ws, Rp, As);%此时输入wp和ws都是二维的，输出wc也是两维的
[B, A] = butter(N, wc,'stop');

至此，数字带阻滤波器设计完成了，如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

接着我们画出数字带阻滤波器的幅频特性曲线，代码如下：

%画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);

画图结果如下：

在这里插入图片描述

4.3 脉冲响应不变法设计

  • 脉冲响应不变法就是先计算模拟滤波器对应的参数，设计好模拟滤波器，接着使用impinvar函数转为数字滤波器。
  • 脉冲响应不变法一般只能用来设计低通和带通滤波器，对于高通和带阻滤波器容易出现频谱混叠现象。
  • 按照这两个参数设计模拟低通滤波器
  • % 脉冲响应不变法低通
Fs = 1; %采样频率
%此处的wp和ws都是模拟角频率
wp = 0.2 * pi * Fs; 
ws = 0.35 * pi * Fs;
Rp = 1;
As = 10;
 
%低通模拟滤波器设计，注意加上's'
[N, wc] = buttord(wp, ws, Rp, As,'s');
[Bs, As] = butter(N, wc, 's');
 
%把模拟滤波器转换为数字滤波器
[B,A] = impinvar(Bs, As, Fs);

  • 至此，脉冲响应不变法的数字低通滤波器设计完成了

    如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

    接着我们画出数字低通滤波器的幅频特性曲线，代码如下：

  • %画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);
  • 结果如下：
  • 在这里插入图片描述

  • 4.4 双线性变换法设计
    与脉冲响应不变法不同，双线性变化法利用下式

  • 然后设计模拟滤波器，再转化为数字滤波器滤波器设计代码如下：
  • %双线性变化法低通
T = 1;
wpz = 0.2 * pi;
wsz = 0.3 * pi;
Rp = 1;
As = 15;
 
wp = (2 / T) * tan(wpz / 2); %双线性变换
ws = (2 / T) * tan(wsz / 2);%双线性变换
 
[N, wc] = buttord(wp, ws, Rp, As,'s');
[Bs, As] = butter(N, wc, 's');
 
[B,A] = bilinear(B,A, 1/T);%将模拟滤波器转化为数字滤波器

  • 至此，双线性变换法的数字低通滤波器设计完成了

    如果有输入噪声信号x的话，调用y = filter(B,A,x)，得到的y就是滤波后的信号了。

    接着我们画出数字低通滤波器的幅频特性曲线，代码如下：

  • %画图
w = 0 : 0.01 * pi : pi;%取点，从0-pi，每隔 0.01 pi 取一个点
Hk = freqz(B,A,w);%对于取的每个点，求该处的频率响应得下
 
figure
plot(w / pi, 20 * log10(abs(Hk)));%横坐标单位是数字频率对pi的归一化值,纵坐标单位是dB,
grid on;
%设置横纵坐标标签
xlabel('\omega / \pi');
ylabel('幅度/dB');
%设置横纵坐标轴范围
axis([0, 1, -100, 5]);
  • 绘图结果如下：
  • 在这里插入图片描述