特别说明：完整45期数字信号处理教程，原创高性能示波器代码全开源地址：链接
/ H) T; \2 T* N: H. w! ?

第39章 IIR滤波器的Matlab设计

% I5 k2 l7 _6 C# m8 q, ~

    本章节讲解IIR滤波器的Matlab设计，主要包括巴特沃斯滤波器，切比雪夫I型和II型滤波器以及椭圆滤波器。

    39.1 巴特沃斯滤波器的设计

    39.2 切比雪夫滤波器的设计

    39.3 椭圆滤波器的设计

    39.4 总结

+ F' Y. N! q6 t' p39.1 巴特沃斯滤波器的设计
39.1.1 butter函数

功能：用于设计Butterworth（巴特沃斯）滤波器

语法：[b, a] = butter(n, Wn);

说明：butter函数可以设计低通、带通、高通和带阻数字滤波器，其特性可以使通带内的幅度响应最大限度地平坦，但会损失截止频率处的下降斜度，使幅度响应衰减较慢。

l [b,a] = butter(n,Wn)可以设计截止频率为Wn的n阶低通butterworth滤波器，其中截止频率Wn应满足0Wn1，Wn=1相当于0.5fs（采样频率）。当Wn = [W1 W2]时，butter函数产生一个2n阶的数字带通滤波器，其通带为W1 < W < W2。

l [b,a] = butter(n,Wn,'ftype')可以设计高通或带阻滤波器。当ftype=high时，可设计截止频率为Wn的高通滤波器；当ftype=stop时，可设计带阻滤波器，此时Wn = [W1 W2]，阻带为W1 < W < W2。

    使用butter函数设计滤波器，可以使通带内的幅度响应最大地平坦，但会损失截止频率处的下降斜度。因此，butter函数主要用于设计通带平坦的数字滤波器。

39.1.2 buttord函数

功能：用来选择Butterworth滤波器的阶数。

语法：[n, Wn] = buttord（WP, WS, RP, Rs）；

说明：buttord函数可以在给定滤波器性能的情况下，选择Butterworth数字滤波器的最小阶数，其中WP和WS分别是通带和阻带的截止频率，其值为0Wp（或Ws）1，当该值为1时表示0.5fs（采样率）。RP和 Rs分别是通带和阻带区的波纹系数和衰减系数。

    [n, Wn] = buttord（WP, WS, RP, Rs）可以得到高通、带通和带阻滤波器的最小阶数n。

    当WP>WS时，为高通滤波器；当WP, WS为二元矢量时，若WP<WS，则为带通和带阻滤波器，此时Wn也为二元矢量。

    利用buttord函数可得到Butterworth数字滤波器的最小阶数n，并使通带（0，WP）内的波纹系数小于RP,阻带（WS, 1）内衰减系数大于Rs。buttord函数还可以得到截止频率Wn，再利用butter函数可产生满足指定性能的滤波器。

    使用butter函数设计数字滤波器，可以使通带内的幅度响应最大限度地平坦，但在截止频率附件幅度响应衰减慢。如果期望幅度响应下降斜度大，衰减快，可使用Elliptic(椭圆)或Chebyshev（切比雪夫）滤波器。

39.1.3 巴特沃斯低通滤波器设计

    下面我们通过一个实例来讲解巴特沃斯低通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组低通滤波器系数，其阶数是2，截止频率为0.25（也就是125Hz），采样率1Kbps。Matlab运行代码如下：

1. fs=1000;               %设置采样频率 1k
   
2. N=1024;               %采样点数
   . d* v/ D& W+ d3 o0 q/ \: `* A
3. n=0:N-1;
   , u3 q  Z# G( Q" }1 M6 s
4. t=0:1/fs:1-1/fs;         %时间序列         
   
5. f=n*fs/N;              %频率序列
   
6. 
   
7. x1=sin(2*pi*50*t);      %信号
   1 g& ^7 a  G7 R# U
8. x2=sin(2*pi*200*t);     %噪声
   
9. x=x1+x2;               %信号混合
   2 B+ Z+ Y% u8 d: I% L
10. 
    6 E% R- G* u' r/ }
11. subplot(221);
    5 ]" C# M1 b; h/ _! Y2 E
12. plot(t,x);               %绘制原始信号
    0 P2 i6 y/ u6 Z3 z
13. xlabel('时间');
    
14. ylabel('幅值');
    ( b4 B5 ?: H% |
15. title('原始信号');
    
16. grid on;
    
17. 
    
18. subplot(222);
    
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    
22. ylabel('振幅');
    
23. title('原始信号FFT');
    
24. grid on;
    + `! W" f* v8 z+ {4 `
25. 
    4 i3 n# b, k8 M- Y
26. subplot(223);            
    
27. Wc=2*125/fs;           %设置截止频率125Hz                 
    & b6 ~8 T2 K: h5 [
28. [b,a]=butter(2,Wc);   %获取2阶IIR滤波系数
    
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    % N- y$ q* d+ Q$ i
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    3 k/ b, t4 n5 K
32. plot(f,abs(y3));
    2 j! A6 N( f) X5 e* p: ^  j
33. xlabel('频率/Hz');
    / N  U) D2 q; a
34. ylabel('振幅');
    . }; j' x1 C( w8 n
35. title('滤波后信号FFT');
    
36. grid on;
    9 d7 S+ R- w0 W4 \3 F
37. 
    $ h4 ?- ~7 G5 G3 \$ q8 U
38. [H,F]=freqz(b,a,512);
    " v7 d% K) H5 {' _0 L. W
39. subplot(224);
    
40. plot(F/pi,abs(H));
    8 \( s+ S3 a. {' s( b- h
41. xlabel('归一化频率');        %绘制绝对幅频响应
    
42. ylabel('幅度');
    / x" P" d3 A  r! E/ l5 D
43. Ps=sum(x1.^2);            %信号的总功率
    6 D8 S9 o# E6 Z/ [& d: b
44. Pu=sum((y2-x1).^2);       %剩余噪声的功率
    
45. SNR=10*log10(Ps/Pu);     %信噪比
    
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

运算Malab结果如下：

从滤波的效果来看，2阶的IIR滤波器能够达到将近20的信噪比，比使用FIR需要更少的阶数。

39.1.4 巴特沃斯高通滤波器设计

    下面我们通过一个实例来讲解巴特沃斯高通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组高通滤波器系数，其阶数是2，截止频率为0.25（也就是125Hz），采样率1Kbps。Matlab运行代码如下：

1. fs=1000;               %设置采样频率 1k
   7 z2 M7 m9 F. d% z7 U5 s2 q; N
2. N=1024;               %采样点数
   
3. n=0:N-1;
   - n9 R+ ~1 @/ H
4. t=0:1/fs:1-1/fs;         %时间序列         
   ; h6 ^, K: ^4 o8 g0 ?8 O$ w
5. f=n*fs/N;              %频率序列
   
6. 
   
7. x1=sin(2*pi*50*t);      %噪声
   7 P& N6 T: Y# ?& h& I7 V
8. x2=sin(2*pi*200*t);     %信号
   
9. x=x1+x2;              %信号混合
   , o# n5 F  ~9 U& i7 U+ Z  J5 h
10. 
    2 {8 n! d; h9 v3 K1 A. L. M
11. subplot(221);
    ) W0 g' G6 l. [/ r+ k
12. plot(t,x);               %绘制原始信号
    5 V) z* ^! k6 E4 ^: h9 q% e' J2 F& j
13. xlabel('时间');
    + v5 L! o  t' s0 B) q( m
14. ylabel('幅值');
    5 }5 ~% |5 ?; R  z" L9 e+ k
15. title('原始信号');
    ! X6 O! t  y1 E: s6 l0 c5 u. s4 u
16. grid on;
    
17. 
    7 a/ s/ {* J/ E( R$ I% ]9 i
18. subplot(222);
    
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    
22. ylabel('振幅');
    
23. title('原始信号FFT');
    
24. grid on;
    . }# M. O/ Q% ^7 A7 M0 r; w
25. 
    9 ^0 u) @6 }3 Y
26. subplot(223);            
    
27. Wc=2*125/fs;              %设置截止频率125Hz                 
    
28. [b,a]=butter(2,Wc, 'high');   %获取2阶IIR滤波系数
    
29. % y2=filter(b,a,x);
    ! _' d' I9 z1 [0 d
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    ' q* e! }% g& ]& |3 a; y( B9 X/ q
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    ( H  W' c& e5 |! k- |0 ^
32. plot(f,abs(y3));
    
33. xlabel('频率/Hz');
    
34. ylabel('振幅');
    
35. title('滤波后信号FFT');
    
36. grid on;
    ' A3 @. S! }) u
37. 
    $ r7 I4 u* j; J; Q
38. [H,F]=freqz(b,a,512);
    
39. subplot(224);
    
40. plot(F/pi,abs(H));
    + r& d5 E! s( H& m! P& _0 u
41. xlabel('归一化频率');       %绘制绝对幅频响应
    ( I2 p$ v  z# h  ]) Z0 |  E$ x
42. ylabel('幅度');
    
43. Ps=sum(x2.^2);           %信号的总功率
    , P; z: |0 o/ H/ @/ l6 V: Y; _
44. Pu=sum((y2-x2).^2);       %剩余噪声的功率
    
45. SNR=10*log10(Ps/Pu);     %信噪比
    
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    & n1 u1 k, ~) ~4 R" {1 W- V
47. grid on;

复制代码

Matlab运行结果如下：

从滤波的效果来看，2阶的IIR滤波器效果还是比较好的。

+ H. w: k4 E* A+ `39.1.5 巴特沃斯带通滤波器设计

    下面我们通过一个实例来讲解巴特沃斯带通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组带通滤波器系数，其阶数是2，通带为125Hz到300Hz，采样率1Kbps。Matlab运行代码如下：

1. fs=1000;               %设置采样频率 1k
   
2. N=1024;               %采样点数
   ( h  S/ [& f. I; E$ N' H
3. n=0:N-1;
   
4. t=0:1/fs:1-1/fs;         %时间序列         
   
5. f=n*fs/N;              %频率序列
   
6. 
   - w. v" Z* X4 x3 \# m5 q
7. x1=sin(2*pi*50*t);      %噪声
   & y) L" T: Y+ T, [) Q
8. x2=sin(2*pi*200*t);     %信号
   2 @5 ~2 z0 ~% B1 i. z0 `
9. x=x1+x2;              %信号混合
   
10. 
    : Y. Y" C% t2 {4 ^1 q
11. subplot(221);
    
12. plot(t,x);               %绘制原始信号
    
13. xlabel('时间');
    
14. ylabel('幅值');
    
15. title('原始信号');
    
16. grid on;
    7 b% l+ X$ j) B1 n- B7 h
17. 
    ( L, F$ b1 a5 ], F1 |
18. subplot(222);
    
19. y=fft(x,N);             %绘制原始信号的幅频响应
    3 l, Z" `/ Q0 B* y, o
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    
22. ylabel('振幅');
    * Z2 I% [3 D0 u6 d, E# M1 B
23. title('原始信号FFT');
    
24. grid on;
    . F2 s4 L. `/ E+ f2 `# E
25. 
    ) X4 M3 @2 X' N$ _3 Y* s
26. subplot(223);            
    
27. Wn=[125*2 300*2]/fs;  %设置通带125Hz到300Hz                 
    + @7 ]0 C% L( ?
28. [b,a]=butter(1,Wn);     %注意第一个参数虽然是1，但生成的却是2阶IIR滤波器系数
    : C: w$ D5 E1 G8 L3 i
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    ( w4 y. w, i. ?: ?
32. plot(f,abs(y3));
    
33. xlabel('频率/Hz');
    $ K) D5 l1 j- C* I' E  w. l% p8 {
34. ylabel('振幅');
    ' ]' i3 U, U3 x5 \
35. title('滤波后信号FFT');
    - ?8 c/ L- s$ I& w9 {# Z
36. grid on;
    
37. 
    # @: Y" Y. e" U6 I# t2 q2 {% J$ i
38. [H,F]=freqz(b,a,512);
    $ H: @! Q( {8 f) d) F% H& y. p6 x
39. subplot(224);
    7 K. Y5 K: z: l9 }
40. plot(F/pi,abs(H));
    
41. xlabel('归一化频率');        %绘制绝对幅频响应
    
42. ylabel('幅度');
    
43. Ps=sum(x2.^2);            %信号的总功率
    3 \( c$ j5 |9 y( K# ?
44. Pu=sum((y2-x2).^2);        %剩余噪声的功率
    - v# G' z0 }( ~9 K5 A
45. SNR=10*log10(Ps/Pu);   %信噪比
    
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

Matlab的计算结果如下：

从滤波的效果来看，2阶的IIR滤波器效果还是比较好的。

9 l. ~2 D" ~2 P& |3 G' S4 P39.1.6 巴特沃斯带阻滤波器设计

    下面我们通过一个实例来讲解巴特沃斯带阻滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数butter设计一组带阻滤波器系数，其阶数是2，阻带为125Hz到300Hz，采样率1Kbps。Matlab运行代码如下：

1. fs=1000;              %设置采样频率 1k
   
2. N=1024;              %采样点数
   
3. n=0:N-1;
   ) X# W0 l7 M" B5 l
4. t=0:1/fs:1-1/fs;        %时间序列         
   
5. f=n*fs/N;             %频率序列
   : ~# r1 B1 J$ m! F$ f2 `) J: k7 l
6. 
   
7. x1=sin(2*pi*50*t);      %信号
   
8. x2=sin(2*pi*200*t);     %噪声
   2 m( K5 |+ @" y1 ~( @  k
9. x=x1+x2;             %信号混合
   0 h8 i; e/ `' I0 x- h+ b
10. 
    
11. subplot(221);
    9 v+ J- K8 V( q+ o" _
12. plot(t,x);               %绘制原始信号
    ' ~5 _7 j4 V# ^% f
13. xlabel('时间');
    
14. ylabel('幅值');
      _& u/ c9 j  }8 j, t5 w
15. title('原始信号');
    5 k( i6 s! p0 p7 X+ E! v
16. grid on;
    - n; I' `9 }' i& `
17. 
    - \+ x2 ?# k7 Z2 F& ]: v% V& F
18. subplot(222);
    
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    0 m0 I/ ?: \) h
22. ylabel('振幅');
    
23. title('原始信号FFT');
    % ?: j5 V" L/ L
24. grid on;
    : \; n. |  L0 w0 M
25. 
    " |) ]" \; ^! ~
26. subplot(223);            
    + F. W  i1 X3 H" f
27. Wn=[125*2 300*2]/fs;       %设置阻带125Hz到300Hz                 
    1 X  V6 L6 C% J4 L* [! W9 r) m0 B
28. [b,a]=butter(1,Wn, 'stop');   %注意第一个参数虽然是1，但生成的却是2阶IIR滤波器系数
    . ]! e2 G* y( W0 K2 }& @0 v
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    / [7 ~" k: B/ x( F' X
32. plot(f,abs(y3));
    
33. xlabel('频率/Hz');
    ! \5 K; a' Y  T$ L0 j
34. ylabel('振幅');
    
35. title('滤波后信号FFT');
    8 A7 f6 R3 s' J+ W; v8 s7 v
36. grid on;
    
37. 
    3 A! `+ }! y) D$ Y2 Y! ^/ Y
38. [H,F]=freqz(b,a,512);
    , Z- j2 n: R& \* S  X6 N+ v& ~. a6 w
39. subplot(224);
    ) d* o2 r" ^, h# h' T5 m0 z0 ?7 o
40. plot(F/pi,abs(H));
    
41. xlabel('归一化频率');        %绘制绝对幅频响应
    . Z6 {8 j7 [6 {# _6 h6 m6 D
42. ylabel('幅度');
    # \4 E8 m" k( b9 r7 \
43. Ps=sum(x1.^2);         %信号的总功率
    5 K, b/ e, ^3 e  |$ \
44. Pu=sum((y2-x1).^2);     %剩余噪声的功率
    3 q3 }$ t$ ]  A; g" x, O5 `; M1 k
45. SNR=10*log10(Ps/Pu);   %信噪比
    6 o9 w" @9 V& p" r6 v0 ]# q5 C. s
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    ! Z: p5 L% I+ i5 M  c- @2 q0 w
47. grid on;

复制代码

Maltab运行结果如下：

从滤波的效果来看，2阶的IIR滤波器效果还是比较好的。

39.2 切比雪夫滤波器的设计

    切比雪夫（Chebyshev）滤波器分为Chebyshev I型和Chebyshev II型，分别具有通带等纹波和阻带等纹波性能。

39.2.1 cheby1函数

功能：用来设计Chebyshev（切比雪夫）I型滤波器（通带等纹波）。

语法：[b,a] = cheby1(n, Rp, Wn);

            [b,a] = cheby1(n, Rp, Wn, 'ftype');

说明：cheby1函数可以设计低通，带通，高通和带阻Chebyshev I型数字滤波器，其通带内为等纹波，阻带内为单调。Chebyshev I型滤波器的下降斜度比Chebyshev II型大，但其代价是在通带内纹波较大。

    [b,a] = cheby1(n, Rp, Wn);可以设计n阶低通Chebyshev I型数字滤波器，其中RP用来确定通带内的纹波，Wn为改滤波器的截止频率。

    当Wn=[W1, W2]时，cheby1函数可产生一个2n的数字带通滤波器，其通带为W1<W<W2。

    [b,a] = cheby1(n, Rp, Wn, 'ftype');可用来设计n阶高通或带阻滤波器，其中Rp和Wn同上，ftype的定义与butter相同。

( X9 X3 h+ ?4 M' i3 i5 f' [39.2.2 cheby1ord函数

功能：用来选择Chebyshev I型滤波器的阶数。

语法：[n, Wn] = cheb1ord（WP, WS, RP, Rs）;

说明：cheb1ord函数可以在给定滤波器性能的情况下，选择Chebyshev I型数字滤波器的最小阶数，其中Wp和Ws分别是通带和阻带的截止频率，其值为0Wp（或Ws）1。Rp和Rs分别是通带和阻带区的波纹系数。

    [n, Wn] = cheb1ord（WP, WS, RP, Rs）;可以得到低通、高通、带通和带阻滤波器的最小阶数。

    利用cheblord函数，除了可以得到Chebyshev I型数字滤波器的最小阶数n外，还可以得到截止频率Wn，再利用cheby1函数可产生满足指定性能的滤波器，使滤波器通带(0, Wp)内的纹波系数小于Rp，阻带(WS, 1)内衰减系数大于RS

39.2.3 cheby2函数

功能：用来设计Chebyshev（切比雪夫）I型滤波器（通带等纹波）。

语法：[b,a] = cheby1(n, Rp, Wn);

            [b,a] = cheby1(n, Rp, Wn, 'ftype');

说明：cheby2函数与cheby1函数基本相同，只是用cheby2函数所设计的滤波器，其通带内为单调的，阻带内为等波纹，由RS指定阻带内的波纹。

    cheby2函数可以设计低通，带通，高通和带阻Chebyshev II型数字滤波器。

: Z9 s$ w- l7 l. u( c39.2.4 cheby2ord函数

功能：用来选择Chebyshev II型滤波器的阶数。

语法：[n, Wn] = cheb2ord（WP, WS, RP, Rs）;

说明：cheb2ord函数与cheb2函数类似，可以利用该函数确定Chebyshev II型数字滤波器的最小阶数n和截止频率Wn。

     cheb2ord函数和cheb2函数配合使用，可设计出最低阶数的Chebyshev II型数字滤波器。

39.2.5 切比雪夫I型低通滤波器设计

    下面我们通过一个实例来讲解切比雪夫I型低通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数cheby1设计一组低通滤波器系数，其阶数是4，截止频率为0.25（也就是125Hz），采样率1Kbps，通带波纹1db。Matlab运行代码如下：

1. fs=1000;               %设置采样频率 1k
   * t7 L% `. l- a, \
2. N=1024;               %采样点数
   
3. n=0:N-1;
   
4. t=0:1/fs:1-1/fs;         %时间序列         
   
5. f=n*fs/N;              %频率序列
   
6. 
   + c/ m: E4 w5 m6 T1 k  D
7. x1=sin(2*pi*50*t);      %信号
   
8. x2=sin(2*pi*200*t);     %噪声
   # I+ Z, n3 v) F5 \6 r
9. x=x1+x2;              %信号混合
   ) x3 X% n( d+ M; I' S. e
10. 
    9 j3 C7 ?/ |; ]' M: E$ t2 O- t
11. subplot(221);
    % q5 s& A, V: f* w
12. plot(t,x);                %绘制原始信号
    
13. xlabel('时间');
    % V6 c8 O; g1 P( D$ m
14. ylabel('幅值');
    
15. title('原始信号');
    1 R$ h7 j  f( q" m. D9 y
16. grid on;
    8 K7 v: J8 t. H5 |* }& {
17. 
    - Q5 G" ~+ C" c7 R5 X
18. subplot(222);
    ' E' E1 f  t2 V; N% i- ~1 x. J9 [
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    1 @" Z) [. ^" v4 m
22. ylabel('振幅');
    ( D8 R8 h- Z& h5 T; a: p2 |
23. title('原始信号FFT');
    
24. grid on;
      x; m2 R" L4 C0 d! c; O) r
25. 
    % J& t5 K7 \, V& P* G' J
26. subplot(223);            
    
27. Wc=2*125/fs;            %设置截止频率125Hz                 
    
28. [b,a]=cheby1(4, 3, Wc);   %获取2阶IIR滤波系数
    
29. % y2=filter(b,a,x);
    . P+ y! I: ^+ `3 e4 X4 @4 ]! C
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    % o. E- F/ B% e' N+ O3 m
31. y3=fft(y2,N);           %滤波后波形的幅频响应
    
32. plot(f,abs(y3));
    
33. xlabel('频率/Hz');
    
34. ylabel('振幅');
    " @5 W4 a3 b, L3 i# {, p
35. title('滤波后信号FFT');
    ; X1 ?8 I+ y$ V1 O! }/ F
36. grid on;
    
37. 
    * F3 h1 l7 I3 x: `7 ?( E7 o
38. [H,F]=freqz(b,a,512);
    " D6 j" g9 ^. E4 m
39. subplot(224);
    
40. plot(F/pi,abs(H));
    
41. xlabel('归一化频率');        %绘制绝对幅频响应
    ! P- v% D% M+ u8 F* N( q. x7 [7 Y
42. ylabel('幅度');
    / W: N" L( p1 H! p! o# ]  K
43. Ps=sum(x1.^2);            %信号的总功率
    3 a  `2 v# p$ X/ a, T0 j% y
44. Pu=sum((y2-x1).^2);       %剩余噪声的功率
    
45. SNR=10*log10(Ps/Pu);     %信噪比
    
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

Matlab计算结果如下：

从滤波的效果来看，4阶的切比雪夫I型滤波效果还是比较好的。

( V6 h* A& M. F* z39.2.6 切比雪夫I型高通滤波器设计

    下面我们通过一个实例来讲解切比雪夫I型高通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数cheby1设计一组高通滤波器系数，其阶数是2，截止频率为0.25（也就是125Hz），采样率1Kbps，通带波纹1db。Matlab运行代码如下：

1. fs=1000;           %设置采样频率 1k
   5 H7 o% [$ n* P+ L* q2 m+ F* x, t  |
2. N=1024;           %采样点数
   ! f7 c8 E# [; B# F  }# I' \
3. n=0:N-1;
   
4. t=0:1/fs:1-1/fs;     %时间序列         
   : c& g" v6 x+ t& i1 K% \
5. f=n*fs/N;          %频率序列
   
6. 
   
7. x1=sin(2*pi*50*t);      %噪声
   4 c: H; f' _1 S' L
8. x2=sin(2*pi*200*t);     %信号
   % |; z7 P: [( {) _9 q
9. x=x1+x2;              %信号混合
   % H: Y' q% [) s& W; K$ E7 g# }
10. 
    
11. subplot(221);
    & B& n5 C8 v1 L5 J* X* f: \9 Q
12. plot(t,x);               %绘制原始信号
      N! x4 N+ Z, p* S7 y1 c
13. xlabel('时间');
    
14. ylabel('幅值');
    " Q  \* ~. }1 R( \  X+ d
15. title('原始信号');
    & W0 o0 I% w8 @5 Q# c
16. grid on;
    
17. 
    
18. subplot(222);
    + H+ v3 u" b1 V, L4 I" M
19. y=fft(x,N);             %绘制原始信号的幅频响应
    0 R/ Q6 S' p. f0 _! P( C
20. plot(f,abs(y));
    6 S; ~6 S9 L0 S# m" }5 H0 q# X
21. xlabel('频率/Hz');
    ; b; b, q3 _- ?4 E
22. ylabel('振幅');
    4 J  Z9 q2 _5 W8 K# Z
23. title('原始信号FFT');
    
24. grid on;
    
25. 
    
26. subplot(223);            
    
27. Wc=2*125/fs;              %设置截止频率125Hz                 
    
28. [b,a]=cheby1(2, 1, Wc, 'high');  %获取2阶IIR滤波系数
    
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    + t. b: K: Z1 h0 I
32. plot(f,abs(y3));
    1 U1 A' ~! F/ ^' o+ `
33. xlabel('频率/Hz');
    
34. ylabel('振幅');
    
35. title('滤波后信号FFT');
    , Z/ @) ~# t# @. D5 n
36. grid on;
    * |1 B# Q" E& ~3 j* P8 c
37. 
    
38. [H,F]=freqz(b,a,512);
    : |8 i* c+ c& m) z, \/ v% {
39. subplot(224);
    
40. plot(F/pi,abs(H));
    
41. xlabel('归一化频率');       %绘制绝对幅频响应
    
42. ylabel('幅度');
    
43. Ps=sum(x2.^2);           %信号的总功率
    4 p- t8 ~3 X7 q" e; I
44. Pu=sum((y2-x2).^2);       %剩余噪声的功率
    
45. SNR=10*log10(Ps/Pu);     %信噪比
    
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

Matlab运行结果如下：

39.2.7 切比雪夫I型带通滤波器设计

    下面我们通过一个实例来讲解切比雪夫I型带通滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将50Hz的正弦波当做噪声滤掉，下面通过函数cheby1设计一组带通滤波器系数，其阶数是2，通带为125Hz到300Hz，采样率1Kbps，通带纹波1db。Matlab运行代码如下：

1. fs=1000;               %设置采样频率 1k
   + `4 h. N) Q8 U4 z9 |9 M* {
2. N=1024;               %采样点数
   ' M+ s* P- d# Q9 w: M
3. n=0:N-1;
   
4. t=0:1/fs:1-1/fs;         %时间序列         
   4 b+ y4 ]& \+ X  b/ O
5. f=n*fs/N;              %频率序列
   
6. 
   
7. x1=sin(2*pi*50*t);      %噪声
   
8. x2=sin(2*pi*200*t);     %信号
   ( i$ c2 f+ b' r( D/ ?
9. x=x1+x2;              %信号混合
   $ B) a( \, q' a2 V& W3 s. e
10. 
    
11. subplot(221);
      X, M) M& [+ A5 p- _0 ~
12. plot(t,x);               %绘制原始信号
    0 x! ^2 x: a7 v1 a* f) N, R
13. xlabel('时间');
    + o0 U+ g8 i8 `) Y9 t0 ^5 o2 @1 }
14. ylabel('幅值');
    
15. title('原始信号');
    
16. grid on;
    
17. 
    ) e$ j! J+ ], ]& K
18. subplot(222);
      ]+ H0 p' M5 F- m- I
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    4 A& F8 T9 Q5 c; z, y: v, a. ^
21. xlabel('频率/Hz');
    
22. ylabel('振幅');
    
23. title('原始信号FFT');
    
24. grid on;
    
25. 
    
26. subplot(223);            
    6 _% p5 o- |3 H- x9 k
27. Wn=[125*2 300*2]/fs;  %设置通带125Hz到300Hz                 
    
28. [b,a]=cheby1(1,1,Wn);  %注意第一个参数虽然是1，但生成的却是2阶IIR滤波器系数
    * @5 ]6 b4 P3 B) V
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    : c8 s- Z+ e9 N
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    
32. plot(f,abs(y3));
    % Q" b. _' y. I
33. xlabel('频率/Hz');
    
34. ylabel('振幅');
      R4 n- \( E* G5 K( h
35. title('滤波后信号FFT');
    4 X# O3 {2 [  E8 Y1 B' _
36. grid on;
    % v3 l! @, r( m! R$ Z1 V! C
37. 
    
38. [H,F]=freqz(b,a,512);
    
39. subplot(224);
    
40. plot(F/pi,abs(H));
    1 `5 A# O* F( ^
41. xlabel('归一化频率');      %绘制绝对幅频响应
    
42. ylabel('幅度');
    # Q  u; X- e9 X( Y1 b8 T3 L9 |
43. Ps=sum(x2.^2);          %信号的总功率
      ]  B3 z: K) Y2 L
44. Pu=sum((y2-x2).^2);     %剩余噪声的功率
    " L0 i' E* s2 U7 t
45. SNR=10*log10(Ps/Pu);   %信噪比
    - @% H* D- m$ W* h
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

Matlab的计算结果如下：

从滤波的效果来看，2阶的带通滤波器效果不够好，出现这种情况的时候，需要大家去重新的调节截止频率，滤波器阶数和通带波纹。

39.2.8 切比雪夫I型带阻滤波器设计

    下面我们通过一个实例来讲解切比雪夫I型带阻滤波器的设计。原始信号是由50Hz正弦波和200Hz的正弦波组成，将200Hz的正弦波当做噪声滤掉，下面通过函数cheby1设计一组带阻滤波器系数，其阶数是2，阻带为125Hz到300Hz，采样率1Kbps，通带波纹1db。Matlab运行代码如下：

1. fs=1000;              %设置采样频率 1k
   & n! r3 f( T6 s* c/ y
2. N=1024;              %采样点数
   
3. n=0:N-1;
   ; O4 @* \. K' x& V5 Y" g- A
4. t=0:1/fs:1-1/fs;        %时间序列         
   1 l6 }; t# U; Q- g+ ~( E
5. f=n*fs/N;             %频率序列
   % S" q4 S. M/ ?
6. 
   
7. x1=sin(2*pi*50*t);      %信号
   ; z. e0 i+ U$ I; l
8. x2=sin(2*pi*200*t);     %噪声
   % K* e0 q; E7 |# g
9. x=x1+x2;             %信号混合
   - c. d1 j& N: z) I* n
10. 
    
11. subplot(221);
    
12. plot(t,x);               %绘制原始信号
    % V+ g# ^% O3 m  O+ n  ^# n
13. xlabel('时间');
    
14. ylabel('幅值');
    
15. title('原始信号');
    9 Z* P4 }# ]1 `: W
16. grid on;
    ! X; Z& I6 F/ D
17. 
    
18. subplot(222);
    6 a- j  S1 @& L9 r. [
19. y=fft(x,N);             %绘制原始信号的幅频响应
    
20. plot(f,abs(y));
    
21. xlabel('频率/Hz');
    - y" f, p7 k/ L9 j( c! d0 [
22. ylabel('振幅');
    
23. title('原始信号FFT');
    
24. grid on;
    
25. 
    , K+ `4 q9 L9 A9 Y
26. subplot(223);            
    ' J& T/ X: R; ?) S6 @
27. Wn=[125*2 300*2]/fs;       %设置阻带125Hz到300Hz                 
    , m& l5 H) `6 c$ W, l
28. [b,a]=cheby1(1,1,Wn, 'stop'); %注意第一个参数虽然是1，但生成的却是2阶IIR滤波器系数
    . z5 b- `7 s0 G1 h, S
29. % y2=filter(b,a,x);
    
30. y2=filtfilt(b,a,x);      %计算滤波后的波形y2
    . Z3 v5 M& D2 _
31. y3=fft(y2,N);         %滤波后波形的幅频响应
    ; x" A% G3 A6 }- X) c' r
32. plot(f,abs(y3));
    * {* Y, r/ V* h/ q4 B+ J0 |
33. xlabel('频率/Hz');
    6 H: n  c0 f2 [4 D: s* S1 k# _8 S
34. ylabel('振幅');
    
35. title('滤波后信号FFT');
    0 \/ Y$ T; K0 E
36. grid on;
    * E. ^$ K: Q4 ~1 e* n  p
37. 
    
38. [H,F]=freqz(b,a,512);
    ! [5 |8 R  o8 S* C
39. subplot(224);
    
40. plot(F/pi,abs(H));
    / [# z7 i' z$ {8 g% x+ }2 S+ t, N$ u, A
41. xlabel('归一化频率');        %绘制绝对幅频响应
      _6 A! y+ T: E% ^; U6 N( G* P
42. ylabel('幅度');
    
43. Ps=sum(x1.^2);         %信号的总功率
    % {  G& `6 s6 \' Q9 W7 ]0 J; D* W
44. Pu=sum((y2-x1).^2);     %剩余噪声的功率
    ; P8 p, c5 D& u. B, `; n3 e3 @2 O# X
45. SNR=10*log10(Ps/Pu);   %信噪比
    ' F0 b0 g% ^: j/ p0 a
46. title(['Order=',int2str(2), ' SNR=',num2str(SNR)]);
    
47. grid on;

复制代码

Matlab计算结果如下：

从滤波的效果来看，2阶带阻滤波的效果较好。

39.3 椭圆滤波器的设计
7 _8 W3 h/ h- s) {9 G# s$ f39.3.1 ellip函数

功能：用来设计Elliptic(椭圆)型滤波器

语法：[b,a] = ellip(n, Rp, RS, Wn);

              [b,a] = ellip(n, Rp, RS, Wn, 'ftype');

说明：ellip函数与cheby1、cheby2函数类似，可以设计低通、高通、带通和带阻数字滤波器。参数RP和RS分别用来指定通带波纹和阻带波纹，Wn指定滤波器的截止频率，n为滤波器的阶数。

    与Butterworth和Chebyshev滤波器相比，ellip函数可以得到下降斜度更大、衰减更快的滤波器，但通带和阻带内均为等纹波。通常情况下，椭圆滤波器能以较低的阶数来实现指定的性能。

    [b,a] = ellip(n, Rp, RS, Wn);可设计n阶低通或带通滤波器。当Wn=[W1 W2]时，可设计带通滤波器。

    [b,a] = ellip(n, Rp, RS, Wn, 'ftype');可设计n阶高通或带阻滤波器。

    当ftype=high时，可设计截止频率为Wn的高通滤波器。

    当ftype=stop时，且Wn=[W1 W2]时，可设计带阻滤波器，阻带为W1<W<W2。

  r2 Y8 A+ X' D39.3.2 ellipord函数

功能：用来选择椭圆滤波器的阶数。

语法：[n, Wn] = cheb2ord（WP, WS, RP, Rs）；

说明：ellipord函数与cheb1ord函数类似，用于选择指定性能时的椭圆滤波器的最小阶数n和截止频率Wn，并与ellip函数配合可设计出最低阶数的椭圆滤波器。

39.3.3 椭圆滤波器设计

    关于椭圆滤波器的使用，大家参考前面的切比雪夫滤波器设计即可，使用方法基本是类似的。但是由于椭圆滤波器要同时给我通带和阻带的纹波，所以要得到满足要求的滤波器系数要花些时间去做测试。

' e/ E* D/ l- G8 s

39.4 总结

    本章节主要讲解了巴特沃斯，切比雪夫和椭圆滤波器的设计，如果想用好还需要大家多多做测试，并深入了解相关理论知识。

3 @9 `& T+ A2 O1 R0 p* z$ }

支持原创。。。。

收藏学习