本文是模拟滤波器的设计。 如果你需要知道数字过滤器的内容,你可以按顺序浏览我写的另外两个博客。 如下所示。
2.MATLAB实现无限冲激响应数字滤波器(IIR )
3 .用3.MATLAB实现有限冲激响应数字滤波器(FIR )
目录1 .基础知识介绍2 .函数介绍2.1 buttord -求解滤波器阶数n和3dB截止频率wc2.2 butter -求解次滤波器的具体参数b和a,求解b和a后完成滤波器的设计。 2.3滤波器-滤波函数3 .代码实现: (1)低通滤波器;2 )高通滤波器;3 )带通滤波器;4 )带阻滤波器:
1 .基础知识介绍
让我们先澄清知识。 这非常重要。
某正弦波的信号,频率为50Hz
这意味着信号的模拟频率ff=50(Hz ),注意它的单位是Hz
信号的公式如下
y=sIn(2ft ) sin (250t ) sin (100t ) y=sIn(2ft )=sin (2) pi * 50t (=sin ) 100 ) pit (y=sin ) ) ) )
信号也可以以y=sIn(t ) y=sin ) omegat (y=sin )t )形式表示,所以在此
=2f=100 omega=2 pif=100pi=2f=100
这里的(omega是模拟角频率,它的单位是rad/s)。
注意模拟角频率 Ome
ga Ω 和模拟频率 f f f的关系 Ω = 2 π f Omega=2pi f Ω=2πf 2. 函数介绍首先介绍一些用到的MATLAB函数
2.1 buttord - 求解滤波器的阶数N和3dB截止频率wc [N,wc] = buttord(wp, ws, Rp, As, ‘s’)输入参数如下:
通带边界模拟频率wp、阻带边界模拟频率ws(模拟角频率,单位是rad/s)
通带最大衰减Rp、阻带最小衰减As(单位是dB)
‘s’指的就是模拟滤波器,设计数字滤波器时就没有’s’这个参数了。
2.2 butter - 求解N阶滤波器的具体参数B和A,求解完B和A后滤波器就设计完成了。 [B,A] = butter(N, wc, ‘ftype’, ‘s’) - 模拟滤波器设计输入参数如下:
N - 滤波器阶数
wc - 3dB截止模拟频率(单位rad/s,N和wc都是用buttord函数计算出来的)
ftype - 滤波器类型‘’:
(1)当输入wc为一维向量时:
默认情况下设计的低通滤波器,设计高通滤波器的话令ftype=high
(2)当输入wc为二维向量[wcl,wcu]时:
默认情况下设计的带通滤波器,设计带阻滤波器的话令ftype=stop
这个就是滤波函数了,
x是输入的有噪声的信号,
B,A就是设计好的滤波器参数
得到的输出y就是滤波后的信号了。
3. 代码实现: (1)低通滤波器:例: 设计通带截止频率5kHz,通带衰减2dB,阻带截止频率12kHz,阻带衰减30dB的文静的西牛低通滤波器
由题可知,设计的是模拟滤波器,所以用到下面三个函数:
[N,wc] = buttord(wp, ws, Rp, As, ‘s’)[B,A] = butter(N, wc, ‘ftype’, ‘s’)y = filter(B,A,x)划重点 ! ! !:
模拟滤波器的频率都是模拟角频率 Ω Omega Ω ,它和频率 f f f 的关系
Ω = 2 π f Omega = 2pi f Ω=2πf
所以这里
wp = 2 ∗ p i ∗ 5000 2*pi*5000 2∗pi∗5000,ws = 2 ∗ p i ∗ 12000 2*pi*12000 2∗pi∗12000,Rp = 2, As = 30
代码如下:
wp = 2 * pi * 5000;ws = 2 * pi * 12000;Rp = 2;As = 30;[N, wc] = buttord(wp, ws, Rp, As, 's');[B, A] = butter(N, wc, 's');上面这些代码就设计好了滤波器
如果有输入噪声信号x的话,调用y = filter(B,A,x),得到的y就是滤波后的信号了。
下面是绘图部分
为了让滤波器的结果得到更形象的表示,我们可以画出来它的幅频特性曲线,代码如下:
其中,我们使用了freqs这个函数,
它是用来计算当频率为wk时,对应的频率响应h的大小,主要是用来画图的。
绘图代码如下:
f = 0 : 10 : 14000;%取点,从0-14000,每隔10取一个点w = 2 * pi * f;%注意模拟滤波器用的频率都是模拟角频率,要乘上2pi的 Hk = freqs(B,A,w);%对于取的每个点,求该处的频率响应大小 %画图figureplot(f/1000, 20 * log10(abs(Hk)));%横坐标单位是kHz,纵坐标单位是dB,grid on;%设置横纵坐标标签xlabel('f/kHz');ylabel('-A(f)/dB');%设置横纵坐标轴范围axis([0, 14, -40, 5]);绘图结果如下:
(2)高通滤波器:高通滤波器与低通几乎完全一样,只要注意
[B,A] = butter(N, wc, ‘ftype’, ‘s’)中的 ftype=high
例: 设计通带截止频率4kHz,通带衰减0.1dB,阻带截止频率1kHz,阻带衰减40dB的文静的西牛高通滤波器
代码如下:
wp = 2 * pi * 4000;ws = 2 * pi * 1000;Rp = 0.1;As = 40; [N, wc] = buttord(wp, ws, Rp, As, 's');[B, A] = butter(N, wc,'high', 's');%注意这个'high'高通滤波器设计完成了
如果有输入噪声信号x的话,调用 y = filter(B,A,x),得到的y就是滤波后的信号了。
接着我们画出高通滤波器的幅频特性曲线
f = 0 : 10 : 14000;%取点,从0-14000,每隔10取一个点w = 2 * pi * f;%注意模拟滤波器用的频率都是模拟角频率,要乘上2pi的 Hk = freqs(B,A,w);%对于取的每个点,求该处的频率响应大小 %画图figureplot(f/1000, 20 * log10(abs(Hk)));%横坐标单位是kHz,纵坐标单位是dB,grid on;%设置横纵坐标标签xlabel('f/kHz');ylabel('-A(f)/dB');%设置横纵坐标轴范围axis([0, 14, -60, 5]);曲线图如下:
(3)带通滤波器:例: 设计文静的西牛带通滤波器,通带上下边界频率分别为4kHz和7kHz,通带衰减1dB,阻带上下边界频率2kHz和9kHz,阻带衰减20dB。
滤波器设计代码如下:
%带通wp = 2 * pi * [4000, 7000];ws = 2 * pi * [2000,9000];Rp = 1;As = 20; [N, wc] = buttord(wp, ws, Rp, As, 's');%此时输入wp和ws都是二维的,输出wc也是两维的[B, A] = butter(N, wc,'s');带通模拟滤波器设计完成了
如果有输入噪声信号x的话,调用y = filter(B,A,x),得到的y就是滤波后的信号了。
接着我们画出带通滤波器的幅频特性曲线,如下:
f = 0 : 10 : 14000;%取点,从0-14000,每隔10取一个点w = 2 * pi * f;%注意模拟滤波器用的频率都是模拟角频率,要乘上2pi的 Hk = freqs(B,A,w);%对于取的每个点,求该处的频率响应大小 %画图figureplot(f/1000, 20 * log10(abs(Hk)));%横坐标单位是kHz,纵坐标单位是dB,grid on;%设置横纵坐标标签xlabel('f/kHz');ylabel('-A(f)/dB');%设置横纵坐标轴范围axis([0, 14, -60, 5]);曲线图如下:
例: 设计文静的西牛带阻滤波器,通带上下边界频率分别为2kHz和9kHz,通带衰减1dB,阻带上下边界频率4kHz和7kHz,阻带衰减20dB。
%带阻wp = 2 * pi * [2000, 9000];ws = 2 * pi * [4000,7000];Rp = 1;As = 20; [N, wc] = buttord(wp, ws, Rp, As, 's');%此时输入wp和ws都是二维的,输出wc也是两维的[B, A] = butter(N, wc,'stop','s');带阻模拟滤波器设计完成了,如果有输入噪声信号x的话,调用y = filter(B,A,x),得到的y就是滤波后的信号了。接着我们画出带阻滤波器的幅频特性曲线,代码如下:
f = 0 : 10 : 14000;%取点,从0-14000,每隔10取一个点w = 2 * pi * f;%注意模拟滤波器用的频率都是模拟角频率,要乘上2pi的 Hk = freqs(B,A,w);%对于取的每个点,求该处的频率响应得下 %画图figureplot(f/1000, 20 * log10(abs(Hk)));%横坐标单位是kHz,纵坐标单位是dB,grid on;%设置横纵坐标标签xlabel('f/kHz');ylabel('-A(f)/dB');%设置横纵坐标轴范围axis([0, 14, -100, 5]);结果如下: