请教，应用巴特沃斯低通滤波器进行滤波幅值变小问题

举个例子,希望有所帮助\x0d代码clc; clear all; close all;\x0dwp=50002pi;\x0dws=150002pi;ap=1;as=70;Fs=15000;Fp=9000;%选择滤波器的最小阶数\x0d[N,Wc]=buttord (wp,ws,ap,as,'s');\x0d%创建Butterworth低通滤波器原型\x0d[Z,P,K]=buttap(N);\x0d%零极点增益模型转换为状态空间模型\x0d[A,B,C,D]=zp2ss(Z,P,K);\x0d%实现低通向低通的转变\x0d[AT,BT,CT,DT]=lp2lp(A,B,C,D,Wc);\x0d%状态空间模型转换为传递函数模型\x0d[num1,den1]=ss2tf(AT,BT,CT,DT);\x0d%运用双线性变换法把模拟滤波器转换成数字滤波器\x0d[num2,den2]=bilinear(num1,den1,35000);\x0d%求频率响应\x0d[H,W]=freqz(num2,den2);\x0d%绘出频率响应曲线figure;plot(WFs/(2pi),abs(H));grid;\x0dxlabel('频率/Hz');ylabel('幅值');结果

课 程 设 计

课程设计名称： 数字信号处理

数字信号处理 专业课程设计任务书

题 目 用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器

主要内容

用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器，

要求通带边界频率为400Hz，500Hz，阻带边界频率分别为350Hz，550Hz，通带最大衰减1dB，阻带最小衰减40dB，抽样频率为2000Hz，用MATLAB画出幅频特性，画出并分析滤波器传输函数的零极点；

信号 经过该滤波器，其中 450Hz， 600Hz，滤波器的输出 是什么？用Matlab验证你的结论并给出 的图形。

任务要求

1、掌握用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器的原理和设计方法。

2、求出所设计滤波器的Z变换。

3、用MATLAB画出幅频特性图。

4、验证所设计的滤波器。

参考文献

1、程佩青著，《数字信号处理教程》，清华大学出版社，2001

2、Sanjit K Mitra著，孙洪，余翔宇译，《数字信号处理实验指导书（MATLAB版）》，电子工业出版社，2005年1月

3、郭仕剑等，《MATLAB 7x数字信号处理》，人民邮电出版社，2006年

4、胡广书，《数字信号处理 理论算法与实现》，清华大学出版社，2003年

审查意见 指导教师签字：

说明：本表由指导教师填写，由教研室主任审核后下达给选题学生，装订在设计（论文）首页

填 表 说 明

1．“课题性质”一栏：

A．工程设计；

B．工程技术研究；

C．软件工程（如CAI课题等）；

D．文献型综述；

E．其它。

2．“课题来源”一栏：

A．自然科学基金与部、省、市级以上科研课题；

B．企、事业单位委托课题；

C．校、院（系、部）级基金课题；

D．自拟课题。

1 需求分析

本课程设计要求用双线性变换法设计原型低通为巴特沃兹型的数字IIR带通滤波器，给定的设计参数为通带边界频率为400Hz，500Hz，阻带边界频率分别为350Hz，550Hz，通带最大衰减1dB，阻带最小衰减40dB，抽样频率为2000Hz。要求采用双线性变换法，使得s平面与z平面是单值的一一对应关系，不存在频谱混淆的问题。由于数字频域和模拟频域的频率不是线性关系，这种非线性关系使得通带截止频率、过渡带的边缘频率的相对位置都发生了非线性畸变。设计出巴特沃兹型的带通滤波器之后，让信号 经过该滤波器，其中 450Hz， 600Hz，分析滤波器的输出 是什么，用Matlab验证结论并给出 的图形。

2 概要设计

1>用双线性变换法设计IIR数字滤波器

脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用z=esT转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1-3所示。

图1-3双线性变换的映射关系

为了将S平面的整个虚轴jΩ压缩到S1平面jΩ1轴上的-π/T到π/T段上，可以通过以下的正切变换实现

（1-5）

式中,T仍是采样间隔。

当Ω1由-π/T经过0变化到π/T时，Ω由-∞经过0变化到+∞，也即映射了整个jΩ轴。将式（1-5）写成

将此关系解析延拓到整个S平面和S1平面，令jΩ=s，jΩ1=s1，则得

再将S1平面通过以下标准变换关系映射到Z平面

z=es1T

从而得到S平面和Z平面的单值映射关系为：

（1-6）

（1-7）

式（1-6）与式（1-7）是S平面与Z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换

式（1-5）与式（1-6）的双线性变换符合映射变换应满足的两点要求。

首先,把z=ejω，可得

（1-8）

即S平面的虚轴映射到Z平面的单位圆。

其次，将s=σ+jΩ代入式（1-8），得

因此

由此看出，当σ<0时，|z|<1；当σ>0时，|z|>1。也就是说，S平面的左半平面映射到Z平面的单位圆内，S平面的右半平面映射到Z平面的单位圆外，S平面的虚轴映射到Z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。

2>双线性变换法优缺点

双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为S平面与Z平面是单值的一一对应关系。S平面整个jΩ轴单值地对应于Z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（1-8）所示，重写如下：

上式表明，S平面上Ω与Z平面的ω成非线性的正切关系，如图7-7所示。

由图7-7看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。

图1-4双线性变换法的频率变换关系

但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（1-8）及图1-4所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图1-5所示。

图1-5双线性变换法幅度和相位特性的非线性映射

对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。

根据以上IIR数字滤波器设计方法，下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器，通带截止频率fp1=500HZ,fp2=400HZ;阻带截止频率fs1=350HZ,fs2=550HZ；通带最大衰减αp=1dB;阻带最小衰减αs=40dB，抽样频率为2000HZ。

I设计步骤：

(1)根据任务,确定性能指标:

通带截止频率fp1=500HZ,fp2=400HZ；

阻带截止频率fs1=350HZ,fs2=550HZ；

通带最大衰减αp=1dB；

阻带最小衰减αs=40dB；

抽样频率为2000HZ。

(2)用Ω=(2/T)tan(w/2)对带通数字滤波器H(z)的数字边界频率预畸变，得到带通模拟滤波器H(s)的边界频率主要是通带截止频率fp1,fp2;阻带截止频率fs1，fs2的转换。

为了计算简便，对双线性变换法一般T=2s

通带截止频率fc1=(2/T)tan(2pifp1/2)

fc2=(2/T)tan(2pifp2/2)

阻带截止频率fr1=(2/T)tan(2pifs1/2)

fr2=(2/T)tan(2pifs2/2)

阻带最小衰减αs=1dB和通带最大衰减αp=40dB;

(3)运用低通到带通频率变换公式λ=(((f^2)-(f0^2))/ f)将模拟带通滤波器指标转换为模拟低通滤波器指标。

B=fc2-fc1

normwr1=(((fc1fc2)-(fc1^2))/fc1) normwr2=(((fc1fc2)-(fc2^2))/fc2)

normwc1=(((fc1fc2)-(fr1^2))/fr1)

normwc2=(((fc1fc2)-(fr1^2))/fr1)

模拟低通滤波器指标：normfc，normfr，αp，αs

(4)设计模拟低通原型滤波器。借助巴特沃斯(Butterworth)滤波器,用模拟低通

滤波器设计方法得到模拟低通滤波器的传输函数H(s);

(5)调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。

(6)利用双线性变换法将模拟带通滤波器H(s)转换成数字带通滤波器H(z)

II程序流程框图：

开始

↓

读入数字滤波器技术指标

↓

将指标转换成归一化模拟低通滤波器的指标

↓

设计归一化的模拟低通滤波器阶数N和1db截止频率

↓

模拟域频率变换，将G(P)变换成模拟带通滤波器H(s)

↓

用双线性变换法将H(s)转换成数字带通滤波器H(z)

↓

输入信号后显示相关结果

↓

结束

3 运行环境

PC机 MATLAB

4 开发工具和编程语言

MATLAB语言

5 详细设计

clc;clear all;

%设计滤波器

%所需设计的带通滤波器的参数设置

Fres=2000;

Ts=1/Fres;

Omgap1=500

Omgap2=400

Omgas1=350

Omgas2=550

%进行双线性变换，使得s平面与z平面是单值的一一对应关系

Omgap=(Omgap1-Omgap2)Ts

Omgap3=tan(Omgap1Ts2pi/2)

Omgap4=tan(Omgap2Ts2pi/2)

Omgas3=tan(Omgas1Ts2pi/2)

Omgas4=tan(Omgas2Ts2pi/2)

%对参数归一化

ap1=Omgap3/Omgap

ap2=Omgap4/Omgap

as1=Omgas3/Omgap

as2=Omgas4/Omgap

%将模拟带通滤波器指标转化成模拟低通滤波器指标

I1=(ap1ap2-as1as1)/as1

I2=(ap1ap2-as2as2)/as2

I3=(ap1ap2-ap1ap1)/ap1

I4=(ap1ap2-ap2ap2)/ap2

I5=min(I1,-I2)

alfpp=1;

alfps=40;

%设计巴特沃兹型低通滤波器HL（s）

[N,Wn]=buttord(I4,I5,alfpp,alfps,'s')

[Num,Den]=butter(N,Wn,'s')

%将设计的模拟低通滤波器传递函数HL（s）转化成模拟带通滤波器H（s）

[Num2,Den2]=lp2bp(Num,Den,sqrt(Omgap3Omgap4),Omgap);

%将设计的模拟带通滤波器H（s）转化成数字带通滤波器H（Z）

[Num3,Den3]=bilinear(Num2,Den2,05);

%画出H（Z）幅频特性曲线

w=0:pi/255:pi;

h=freqz(Num3,Den3,w);

H=20log10(abs(h));

plot(w,H);

%设计信号函数

f1=450;f2=600;

t=0:1/2000:40/f2;

f=sin(2pif1t)+sin(2pif2t);

plot(f);

%信号通过带通滤波器

g=invfreqz(h,w,40,50)

g1=conv(g,f)

plot(g1)

6 调试分析

在设计滤波器的时候，计算的是幅频特性。而相频特性没有进行频谱分析。而信号和信号通过滤波器的图形分析都为时与分析，没有进行频域分析，是不知道该用哪个函数对其进行傅里叶变换。试过freqz,fft等函数，但是都没有出来结果。说明调用格式不正确。

本次试验用的是巴特沃兹型低通滤波器设计的方法，还可以利用切比雪夫型低通滤波器设计，或者椭圆型等。

在本次试验中，我们只用了两种频率简单的叠加作为输入信号，但实际信号是多种频率的组合，可以再增加一些频率。

现实中还应该有噪声的影响，本实验中没有考虑。可以再加上噪声信号。

7 测试结果

Fres =

2000

Ts =

50000e-004

Omgap1 =

500

Omgap2 =

400

Omgas1 =

350

Omgas2 =

550

Omgap =

00500

Omgap3 =

10000

Omgap4 =

07265

Omgas3 =

06128

Omgas4 =

11708

ap1 =

200000

ap2 =

145309

as1 =

122560

as2 =

234170

I1 =

114562

I2 =

-110065

I3 =

-54691

I4 =

54691

I5 =

110065

alfpp =

1

alfps =

40

N =

8

Wn =

61894

Num =

10e+006

Columns 1 through 8

0 0 0 0 0 0 0 0

Column 9

21538

Den =

10e+006

Columns 1 through 8

00000 00000 00005 00052 00377 01984 07386 17837

Column 9

21538

Num2 =

10e-004

Columns 1 through 8

08413 00000 -00000 00000 -00000 00000 -00000 00000

Column 9

-00000

Den2 =

Columns 1 through 8

10000 15863 70705 87151 205005 199985 321353 248474

Columns 9 through 16

299185 180527 169631 76698 57123 17643 10400 01695

Column 17

00776

Num3 =

10e-004

Columns 1 through 8

00043 00000 -00341 00000 01194 00000 -02389 00000

Columns 9 through 16

02986 00000 -02389 00000 01194 00000 -00341 00000

Column 17

00043

Den3 =

Columns 1 through 8

10000 -22463 83946 -134519 276744 -336694 483386 -457295

Columns 9 through 16

495687 -364140 306576 -169904 111207 -42944 21332 -04524

Column 17

01603

h =

Columns 1 through 5

-00000 -00000 + 00000i -00000 + 00000i -00000 + 00000i -00000 + 00000i

……

Columns 251 through 256

-00000 - 00000i -00000 - 00000i -00000 - 00000i -00000 - 00000i -00000 - 00000i -00000 - 00000i

H =

Columns 1 through 8

-2792737 -2792732 -2792723 -2792813 -2793855 -2800020 -2830382 -2924507

……

Columns 249 through 256

-2810032 -2803352 -2801410 -2800979 -2800943 -2800973 -2800996 -2801004

f1 =

450

f2 =

600

t =

Columns 1 through 8

0 00005 00010 00015 00020 00025 00030 00035

……

Columns 129 through 134

00640 00645 00650 00655 00660 00665

f =

Columns 1 through 8

0 19387 -02788 -14788 03633 07071 -01420 01338

……

Columns 129 through 134

-03633 -07946 10000 11075 -15388 -10418

g =

Columns 1 through 8

00000 00000 -00000 -00000 00000 00000 -00001 -00002

……

Columns 33 through 41

-00002 00000 00001 -00000 -00000 -00000 00000 00000 00000

g1 =

Columns 1 through 8

0 00000 00000 -00000 -00000 00001 00001 -00003

……

Columns 169 through 174

00000 00000 -00000 -00000 -00000 -00000

综合上述分析，数据基本正确。

上图为带通滤波器幅频响应

、

上图为传输函数零极点分析

上图为输入信号时域

上图为输出波形的时域分析

参考文献

[1] 吴大正 《信号与线性系统分析》第四版 高等教育出版社

[2] 郑君里 《信号与系统》第二版 高等教育出版社

[3] Sanjit K Mitra 《数字信号处理—基于计算机的方法》第三版 清华大学出版社

[4] 余成波 《数字信号处理及MATLAB实现》清华大学出版社

[5] 周利清 《数字信号处理基础》 北京邮电大学出版社

[6] 美国莱昂 《数字信号处理》英文第二版 机械工业出版社

心得体会

00