常见滤波器简要介绍及Matlab实现
常见滤波器简要介绍及Matlab实现
简要收集关于“椭圆、切比雪夫、巴特沃斯、贝塞尔”滤波器的使用与特点，不完善之处还请提出意见与谅解！
1.椭圆滤波器Elliptic Filter
①椭圆滤波器又称考尔滤波器（Cauer filter），是在通带和阻带等波纹特性的一种滤波器，因此通带，阻带逼近特性良好。
②相较其他类型的滤波器，椭圆滤波器在阶数相同的条件下有着最小的通带和阻带波动。
③相比其他滤波器其阻带下降最快，过渡带更为陡峭、更窄的特性（通带和阻带的起伏为代价来换取）。相同阶数时，椭圆滤波 器的幅频曲线下降最陡，其次为切比雪夫滤波器，再次为巴特沃斯滤波器，下降最平缓的为贝塞尔滤波器。
Matlab实现方法：
[n,Wn]=ellipord(wp,ws,Rp,Rs)；
椭圆滤波器最小阶数n和截止频率Wn的确定函数，并使滤波器在通带内(0 , Wp)的波纹系数小于通带最大衰减 R p ,阻带内(W s , 1)的波纹系数大于阻带最小衰减 R s 。Rp、Rs分别为通带最大波纹和阻带最小衰减；wp、ws分别为为通带边界频率和阻带边界频率，单位为rad/s。这四个参数为滤波器的基本性能指标.
[ b, a] = ellip(n , R p , R s , Wn)；
函数的功能是设计滤波器 ，参n表示最小阶数，参数R p表示通带最大波纹度（单位dB），参数R s表示阻带最小衰减（单位dB），参数Wn表示归一化的低通滤波器截止频率（如果阶数n与截止频率不由ellipord计算，可以假设截止频率为F，则有Wn =F/(Fs/2)=F*2/Fs）。
如 : 设计一个带通椭圆数字滤波器 , 通带为 100 ～200 H z,过渡带均为 50 H z,通带波纹小于 3 db ,阻带衰减为 30 db ,采样频 率 f s =1 000 H z 。其程序为 :
-
fs
=
1000
;
-
-
Rp
=3;
Rs
=30
;
-
-
Wp
=2
*[
100
200
]
/fs;
-
-
W
s
=
2
*[
80
220
]
/f
s;
-
-
[
n ,
W
n]
=
ellipo
rd
(Wp
,
W
s
,
Rp
,
Rs);
-
-
[
b ,
a]
=
ellip(n
,
Rp,
Rs
Wn);
-
-
dataOut
=
filter(b,a,dataIn);
2.巴特沃斯滤波器
① 巴特沃斯滤波器是一种在通带和阻带都平坦的滤波器，通带率响应曲线最平坦，没有起伏，阻带频带则逐渐下降为零，下降慢，在过渡带上很容易造成失真。
②在振幅对数与角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少，趋向负无穷大。
③一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝，三阶巴特沃斯滤波器的衰减率为每倍频18分贝，如此类推。
④巴特沃斯滤波器的频率特性曲线，无论在通带内还是阻带内都是频率的单调函数。因此，当通带的边界处满足指标要求时，通带内肯定会有裕量。所以，更有效的设计方法应该是将精确度均匀的分布在整个通带或阻带内，或者同时分布在两者之内。这样就可用较低阶数的系统满足要求。
下图是巴特沃斯滤波器（左上）和同阶第一类切比雪夫滤波器（右上）、第二类切比雪夫滤波器（左下）、椭圆函数滤波器（右下）的频率响应图。
巴特沃斯低通滤波器可用如下振幅的平方对频率的公式表示：
巴特沃斯低通滤波器可用如下振幅的平方对频率的公式表示：
其中, = 滤波器的阶数
= 截止频率 = 振幅下降为 -3分贝时的频率
wp= 通频带边缘频率
在通频带边缘的数值
说明：buttord函数使用阻带指标计算3dB截止频率，这样阻带会刚好满足要求，而通带会有富余。
-
% [B,A] = butter(N,Wn,
'high')
---用来设计高通滤波器
-
% [B,A] = butter(N,Wn,
'low')
--低通滤波器
-
[B,A] = butter(N,Wn)
--带通滤波器
-
-
out_sign=filter(B,A,x_sign); %输入信号最好是转化为double型
为了确定滤波器的最小阶数N，通过[N,Wn]=buttord(Wp,Ws,Rp,Rs)函数获得 。
-
完整示例
8-
30Hz带通滤波：
-
-
fs =
1000 ;
%信号的采样频率
-
wp=[
8
30]*
2/fs;
%通带边界频率 ，单位为rad/s
-
ws=[
7
32]*
2/fs;
%阻带边界频率 ，单位为rad/s
-
Rp=
1;
%通带最大波纹度 ，单位dB (不要太小)
-
Rs=
30;
%表示阻带最小衰减，单位dB
-
[N,Wn]=buttord(Wp,Ws,Rp,Rs);
%巴特沃斯数字滤波器的阶数n和-3dB归一化截止频率Wn
-
-
[B,A]=butter(N,Wn);
%得到n阶巴特沃斯滤波的分子分母
-
dataOut = filter(B,A,dataIn);
3.切比雪夫滤波器
①切比雪夫滤波器是在通带或阻带上频率响应幅度等波纹波动（通带平坦、阻带等波纹或是阻带平坦、通带等波纹）的滤波器，振幅特性在通带内是等波纹。
②切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。
③切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。
④在通带（或称“通频带”）上频率响应幅度等波纹波动的滤波器称为“I型切比雪夫滤波器”，在阻带（或称“阻频带”）上频率响应幅度等波纹波动的滤波器称为“II型切比雪夫滤波器”
-
fs =
1000; % Hz采样频率
-
wp =
55/(fs/
2); %通带截止频率,取
50~
100中间的值,并对其归一化
-
ws =
90/(fs/
2); %阻带截止频率,取
50~
100中间的值,并对其归一化
-
Rp =
3; %通带允许最大衰减为 db
-
Rs =
40; %阻带允许最小衰减为 db
-
[n,Wn]=cheb1ord(Wp,Ws,Rp,Rs); % 获取阶数和截止频率
-
[b,a]=cheby1(n,Rp,Wn, 'low'); %获得转移函数系数
-
dataOut = filter(b,a,dataIn); %信号滤波运算
-
4.贝塞尔滤波器Bessel filter
贝赛尔(l)滤波器是具有最大平坦的群延迟（线性相位响应）的线性过滤器，即最平坦的幅度和相位响应，常用在音频天桥系统中。具有最平坦的幅度和相位响应。带通（通常为用户关注区域）的相位响应近乎呈线性。
由于具有向其截止频率以下的所有频率提供等量延时的特性，才被用于音频设备中，在音频设备中，必须在不损害频带内多信号的相位关系前提下，消除带外噪声。另外，贝塞尔滤波器的阶跃响应很快，并且没有过冲或振铃，这使它在作为音频DAC输出端的平滑滤波器，或音频ADC输入端的抗混叠滤波器方面，是一种出色的选择。贝塞尔滤波器还可用于分析D类放大器的输出，以及消除其它应用中的开关噪声，来提高失真测量和示波器波形测量的精确度。
[b,a] = besself(n,Wn) %低通滤波器（Wn为二元向量时，为带通）
% Wn是一个实际频率值，不是相对量
[b,a] = besself(n,Wn,’ftype’); % ftype类型low，hight，stop
dataOut = filter(b,a,dataIn); %信号滤波运算
四种滤波器的区别对比
①贝塞尔滤波器具有最平坦的幅度和相位响应。带通（通常为用户关注区域）的相位响应近乎呈线性。
②巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。
③切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。
④相同阶数时：椭圆滤波器的幅频曲线下降最陡，其次为切比雪夫滤波器，再次为巴特沃斯滤波器，下降最平缓的为贝塞尔滤波器。
⑤相同阶数时：巴特沃斯滤波器通带最平坦，阻带下降慢；切比雪夫滤波器通带等纹波，阻带下降较快；贝塞尔滤波器通带等纹波，阻带下降慢。也就是说幅频特性的选频特性最差。但是，贝塞尔滤波器具有最佳的线性相位特性；椭圆滤波器在通带等纹波（阻带平坦或等纹波），阻带下降最快。
转插两篇他人博客供大家参考：
博客地址 http://blog.sina.com.cn/s/blog_49c02a8c0100yszh.html （matlab提供了简便的产生各种滤波器的方法）
https://blog.csdn.net/LYduring/article/details/80443573 （几种常见空间滤波器MATLAB实现）
https://blog.csdn.net/colapin/article/details/52840075 （七种滤波方法的matlab实现和测试）
https://www.cnblogs.com/alimy/p/9140695.html （切比雪夫Ⅰ型滤波器设计：低通、高通、带通和带阻）
-
% Cheby1Filter.m
-
% 切比雪夫Ⅰ型滤波器的设计
-
%
-
-
clear;
-
close all;
-
clc;
-
-
fs =
1000;
%Hz 采样频率
-
Ts =
1/fs;
-
N =
1000;
%序列长度
-
t = (
0:N-
1)*Ts;
-
delta_f =
1*fs/N;
-
f1 =
50;
-
f2 =
100;
-
f3 =
200;
-
f4 =
400;
-
x1 =
2*
0.5*sin(
2*pi*f1*t);
-
x2 =
2*
0.5*sin(
2*pi*f2*t);
-
x3 =
2*
0.5*sin(
2*pi*f3*t);
-
x4 =
2*
0.5*sin(
2*pi*f4*t);
-
x = x1 + x2 + x3 + x4;
%待处理信号由四个分量组成
-
-
X = fftshift(abs(fft(x)))/N;
-
X_angle = fftshift(angle(fft(x)));
-
f = (-N/
2:N/
2-
1)*delta_f;
-
-
figure(
1);
-
subplot(
3,
1,
1);
-
plot(t,x);
-
title(
'原信号');
-
subplot(
3,
1,
2);
-
plot(f,X);
-
grid on;
-
title(
'原信号频谱幅度特性');
-
subplot(
3,
1,
3);
-
plot(f,X_angle);
-
title(
'原信号频谱相位特性');
-
grid on;
-
-
%设计一个切比雪夫低通滤波器,要求把50Hz的频率分量保留,其他分量滤掉
-
wp =
55/(fs/
2);
%通带截止频率,取50~100中间的值,并对其归一化
-
ws =
90/(fs/
2);
%阻带截止频率,取50~100中间的值,并对其归一化
-
alpha_p =
3;
%通带允许最大衰减为 db
-
alpha_s =
40;
%阻带允许最小衰减为 db
-
%获取阶数和截止频率
-
[ N1 wc1 ] = cheb1ord( wp , ws , alpha_p , alpha_s);
-
%获得转移函数系数
-
[ b a ] = cheby1(N1,alpha_p,wc1,
'low');
-
%滤波
-
filter_lp_s = filter(b,a,x);
-
X_lp_s = fftshift(abs(fft(filter_lp_s)))/N;
-
X_lp_s_angle = fftshift(angle(fft(filter_lp_s)));
-
figure(
2);
-
freqz(b,a);
%滤波器频谱特性
-
figure(
3);
-
subplot(
3,
1,
1);
-
plot(t,filter_lp_s);
-
grid on;
-
title(
'低通滤波后时域图形');
-
subplot(
3,
1,
2);
-
plot(f,X_lp_s);
-
title(
'低通滤波后频域幅度特性');
-
subplot(
3,
1,
3);
-
plot(f,X_lp_s_angle);
-
title(
'低通滤波后频域相位特性');
-
-
-
-
%设计一个高通滤波器,要求把400Hz的频率分量保留,其他分量滤掉
-
wp =
350/(fs/
2);
%通带截止频率,取200~400中间的值,并对其归一化
-
ws =
380/(fs/
2);
%阻带截止频率,取200~400中间的值,并对其归一化
-
alpha_p =
3;
%通带允许最大衰减为 db
-
alpha_s =
20;
%阻带允许最小衰减为 db
-
%获取阶数和截止频率
-
[ N2 wc2 ] = cheb1ord( wp , ws , alpha_p , alpha_s);
-
%获得转移函数系数
-
[ b a ] = cheby1(N2,alpha_p,wc2,
'high');
-
%滤波
-
filter_hp_s = filter(b,a,x);
-
X_hp_s = fftshift(abs(fft(filter_hp_s)))/N;
-
X_hp_s_angle = fftshift(angle(fft(filter_hp_s)));
-
figure(
4);
-
freqz(b,a);
%滤波器频谱特性
-
figure(
5);
-
subplot(
3,
1,
1);
-
plot(t,filter_hp_s);
-
grid on;
-
title(
'高通滤波后时域图形');
-
subplot(
3,
1,
2);
-
plot(f,X_hp_s);
-
title(
'高通滤波后频域幅度特性');
-
subplot(
3,
1,
3);
-
plot(f,X_hp_s_angle);
-
title(
'高通滤波后频域相位特性');
-
-
-
%设计一个带通滤波器,要求把50Hz和400Hz的频率分量滤掉,其他分量保留
-
wp = [
65
385 ] / (fs/
2);
%通带截止频率,50~100、200~400中间各取一个值,并对其归一化
-
ws = [
75
375 ] / (fs/
2);
%阻带截止频率,50~100、200~400中间各取一个值,并对其归一化
-
alpha_p =
3;
%通带允许最大衰减为 db
-
alpha_s =
20;
%阻带允许最小衰减为 db
-
%获取阶数和截止频率
-
[ N3 wn ] = cheb1ord( wp , ws , alpha_p , alpha_s);
-
%获得转移函数系数
-
[ b a ] = cheby1(N3,alpha_p,wn,
'bandpass');
-
%滤波
-
filter_bp_s = filter(b,a,x);
-
X_bp_s = fftshift(abs(fft(filter_bp_s)))/N;
-
X_bp_s_angle = fftshift(angle(fft(filter_bp_s)));
-
figure(
6);
-
freqz(b,a);
%滤波器频谱特性
-
figure(
7);
-
subplot(
3,
1,
1);
-
plot(t,filter_bp_s);
-
grid on;
-
title(
'带通滤波后时域图形');
-
subplot(
3,
1,
2);
-
plot(f,X_bp_s);
-
title(
'带通滤波后频域幅度特性');
-
subplot(
3,
1,
3);
-
plot(f,X_bp_s_angle);
-
title(
'带通滤波后频域相位特性');
-
-
-
%设计一个带阻滤波器,要求把50Hz和400Hz的频率分量保留,其他分量滤掉
-
wp = [
65
385 ] / (fs/
2);
%通带截止频率?,50~100、200~400中间各取一个值,并对其归一化
-
ws = [
75
375 ] / (fs/
2);
%阻带截止频率?,50~100、200~400中间各取一个值,并对其归一化
-
alpha_p =
3;
%通带允许最大衰减为 db
-
alpha_s =
20;
%阻带允许最小衰减为 db
-
%获取阶数和截止频率
-
[ N4 wn ] = cheb1ord( wp , ws , alpha_p , alpha_s);
-
%获得转移函数系数
-
[ b a ] = cheby1(N4,alpha_p,wn,
'stop');
-
%滤波
-
filter_bs_s = filter(b,a,x);
-
X_bs_s = fftshift(abs(fft(filter_bs_s)))/N;
-
X_bs_s_angle = fftshift(angle(fft(filter_bs_s)));
-
figure(
8);
-
freqz(b,a);
%滤波器频谱特性
-
figure(
9);
-
subplot(
3,
1,
1);
-
plot(t,filter_bs_s);
-
grid on;
-
title(
'带阻滤波后时域图形');
-
subplot(
3,
1,
2);
-
plot(f,X_bs_s);
-
title(
'带阻滤波后频域幅度特性');
-
subplot(
3,
1,
3);
-
plot(f,X_bs_s_angle);
-
title(
'带阻滤波后频域相位特性');
效果：
原始信号：
生成的低通滤波器和滤波后的效果：
生成的高通滤波器和滤波后的结果：
生成的带通滤波器和滤波后的结果：
生成的带阻滤波器和滤波后的结果：
上一篇： 一步一步 的理解设计模式
下一篇： JS 异步（上）