双线性变换法设计巴特沃斯低通数字滤波器.

第一章 概述…………………………………………………………………………2

1.1课程设计的基本原理..........................................2 1.2课程设计目的................................................2

1.3课程设计的内容..............................................2

第二章 总体设计及关键技术分析…………………………………………………4

2.1卷积演示实验................................................4 2.2采样定理演示................................................5 2.3模拟滤波器设计演示..........................................7 2.4设计切比雪夫I型低通滤波器.................................10 2.5双线性变换法设计巴特沃斯低通数字滤波器.....................11 2.6用凯塞窗设计高通滤波器.....................................13 第三章 程序实现………………………………………………………………… 14

3.1卷积演示实验...............................................14 3.2采样定理演示...............................................19 3.3模拟滤波器设计演示.........................................23 3.4设计切比雪夫I型低通滤波器.................................26 3.5双线性变换法设计巴特沃斯低通数字滤波器.....................27 3.6用凯塞窗设计高通滤波器.....................................29 第四章 结束语.……………………………………………………………………32 参考文献……………………………………………………………………………33

第一章 概述

1.1 课程设计的基本原理

数字信号处理的直接对象是数字信号，处理的方式是数值运算的方式，它涉及到的内容非常丰富和广泛，它是应用最快、成效最显著的新学科之一，做为信息专业的专业基础课，需要掌握其基本理论和基本的分析方法，通过理论和实践、原理和应用结合让学生掌握知识。而matlab是集数学计算、图形处理和程序语言设计于一体的软件，通过matlab来实现数字信号处理的有关表达，最终实现课程设计的目的。 1.2课程设计目的

1、通过基于MATLAB的算法仿真实验及分析、基于DSP的算法综合实验等实践活动，进一步领会和深化课堂上学到的有关数字信号处理的基本概念、基本原理以及基本的信号处理操作及滤波器设计方法。

2、掌握线性卷积与圆周卷积软件实现的方法，并验证二者之间的关系；验证奈奎斯特取样定理，加深对时域取样后信号频谱变化的认识；掌握模拟滤波器的频率变换——模拟高通、带通、带阻滤波器的设计与对比演示；以及根据窗函数设计滤波器。 1.3课程设计内容

1. 设计一卷积演示程序。

（1）可输入任意2待卷积序列x(n)、x(n)，长度不做限定。测试数据为：

1

2

x1(n)={1,1,1,1,0,0,1,1,1,1,0,0}，x2(n)={0,1,2,1,0,0,0,1,2,1,0,0}； （2）分别动态演示2序列进行线性卷积x1(n)﹡x2(n)和圆周卷积x1(n) x2 (n)的过程；要求分别动态演示翻转、移位、乘积、求和的过程。

（3）圆周卷积默认使用2序列中的最大长度，但卷积前可以指定卷积长度N用以进行混叠分析。

（4）根据实验结果分析2类卷积的关系。

（5）假定时域序列x1(n)、x2(n)的长度不小于10000，序列内容自定义。利用 FFT实现快速卷积，验证时域卷积定理，并与直接卷积进行效率对比。 2. 演示采样定理，时域采样，频谱周期延拓，同时演示采样频率小于2fc时，产生的混叠效应。

（1）假设频域归一化三角波的频带宽度fc=100Hz，对应的时域信号为：

2

y(t)=fc[sinc(fct/4)]/2。

（2）要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示对应的信号频谱、采样脉冲及频域周期拓延。

注：sinc(x)=sin(πx)/(πx)。

3. 模拟滤波器设计演示—从模拟低通滤波器到模拟高通，带通，带阻的幅度特性对比演示。

（1）设计过程详见教材相关内容。

（2）使用巴特沃斯滤波器，其阶数N应该根据实际参数计算（计算公式和方法如教材所述），为方便作图，这里指定阶数为N=5，并假定通带截止频率f=1，

p

阻带截止频率fs=2。

（3）分别用不同颜色曲线绘制通带、过渡带和阻带。要求根据变换关系动态演示低通滤波器和目标滤波器的幅度特性。

4. 设计一切比雪夫I型低通滤波器，各参数要求：f=100Hz，α=2dB，f=120Hz，

p

p

s

αs=60dB。

给出所设计滤波器的幅度特性并分析是否满足设计需要。

5.使用双线性变换法设计巴特沃斯低通数字滤波器，各设计指标如下： ωs=0.4π，ωp=0.6π，αp=0.5 dB，αs=50 dB

令T=2，给出所设计滤波器的幅度衰减特性及其冲激响应。

6. 利用凯塞窗设计高通滤波器，设计指标分别是：ωs=0.4π，ωp=0.6π，α

p

=0.5dB，αs=60dB。给出窗函数及所设计滤波器的幅度特性。

第二章 总体设计及关键技术分析

2.1 卷积演示实验 2.1.1基本的原理 1、线性卷积：

线性时不变系统（Linear Time-Invariant System, or L. T. I系统）输入、输出间的关系为：当系统输入序列为x(n)，系统的单位脉冲响应为h(n)，输出序列为y(n)，则系统输出为：

y(n)y(n)mmx(m)h(nm)x(n)*h(n)h(m)x(nm)h(n)*x(n)

或

上式称为离散卷积或线性卷积。

2、圆周卷积

设两个有限长序列x1(n)和x2(n)，均为N点长

D F T X(k)x1(n) 1

x2(n) X2(k) 如果X3(k)X1(k)X2(k)

N1~x3(n)x1(m)~x2(nm)RN(n)m0则

x1(m)x2(nm)Nm0N1 x1(n)Nx2(n)

0nN1

上式称为循环卷积或圆周卷积

~~~x(n)R(n)x(n)x(n)1N1注：为序列的周期化序列；为x1(n)的主值序列。

编程计算时，x3(n)可表示如下：

1x3(n)x1(m)x2(nm)m0nmn1x(m)x1N12(Nnm)

3、两个有限长序列的线性卷积

序列x1(n)为L点长，序列x2(n)为P点长，x3(n)为这两个序列的线性卷积，则x3(n)为

x3(n)mx(m)x12(nm)

且线性卷积x3(n)的最大长LP1，也就是说当n1和nLP1时

x3(n)0。 2.1.2设计思想

首先建立一个基本的框架，制作一个菜单，其中包括主程序菜单和子程序的菜单，子程序菜单可以选择回到主程序菜单选择功能。菜单的框架完成后，实现可以任意输入两个序列，然后分别制作动态演示序列的线性卷积的程序、 动态演示序列的圆周卷积、以及验证时域卷机定理以及比较运行速率的程序。结合上面建立的框架完成菜单选择以及功能的调用，让整个设计完美。 2.1.3 设计流程图：

开始 输入x1（n）、x2（n）序列 调用程序主菜单、选择子菜单 动态演示2序列的线性卷积 动态演示序列的圆周卷积 验证时域卷机定理、比较速率 退出程序 结束

2.14 关键技术分析

本个设计主要要实现动态的演示，为实现动态演示，有很多不同的方法，而我采用的是for循环加上pause，在for循环中的一次执行中，实现一次绘图，表示某一时刻的状态，用pause暂停等待for循环中下一次执行，实现动态演示。

这是主要的一个设计，另外一个重要的技巧是如何将循环卷积表示出来，因为循环卷积的结果是周期的，其循环卷积的方法和线性卷积不同，并且循环卷积中要处理几种情况：如当x1序列和x2序列之间的最大长度比卷积N长度大时是一种画图方法，比它小时是另外一种方法，然而困难的是当循环卷积长度比序列最大长度小时还要分析x1和x2序列各自长度与N的关系而做不同的情况分析，具体的解决方法见程序实现。 2.2采样定理演示 2.2.1基本的原理

奈奎斯特取样定理指出：为了使实信号取样后能够不失真还原，取样频率必须大于信号最高频率的两倍。

若x(t)为有限带宽的连续信号，其频谱为X(j)，以T为取样间隔对x(t)ˆ(t)。xˆ(t)的频谱为： 理想取样，得到理想取样信号x

ˆX(j)1TX(jjm2)T

也就是说，一个连续信号经过理想取样后，它的频谱将沿着频率轴，从02sT重复出现一次，即频谱产生周期延拓。 开始，每个一个取样频率2.2.2设计思想

首先通过时域来表示出其波形的变化，从而观察采样前后的波形做比较，第二步通过频域检验其是否混叠，观察采样前后的频域，再通过改变不同的采样频率来分析采样定理。以此设计思想编写好相关程序之后，建立整体框架，实现完成一次采样后可以重新输入采样频率，以便进行分析。 2.2.3设计流程图：

开始

结束 是否继续？1/0 否/0 是/1 继续执行程序，并且在figure2中绘画出时俞信号、采样信号以及采样后的信号的频谱图，如果混叠则画出混叠部分 调用程序并且绘制出时域信号、采样信号以及采样之后的信号的图形 输入采样频率fs

2.2.4 关键技术分析

本次的重点在于如何实现时域原信号的频谱，就是说如何画出三脚波，由于在实验中通过fft不能很好的表示出其频谱，即傅立叶变换，如果按照傅立叶变换求出其频谱也很难实现，所以在实验中采用画波形，得到其原信号的波形之后的关键就在于如何动态地表示出其采样过程以及混叠的过程。动态过程同样使用for循环以及hold on和pause，通过判断fc和fs的关系看是否混叠，混叠则绘制混叠部分，通过数学计算，求出其关系式，在程序中直接应用就可以实现绘制混叠部分的说明，具体见第三章的程序实现部分。 2.3模拟滤波器设计演示 2.3.1基本的原理 巴特沃斯滤波器：

特点：具有通带内最大平坦的振幅特性，且随f↗单调↘ ，其幅度平方函

数具有如下形式：

A(2)Ha(j)21j1jc2N式中，N为整数，称为滤波器的阶数，N越大，通带和阻带的近似性越好，

过渡带也越陡。如图。

图 巴特沃兹filter 振幅平方函数

1) 在通带，分母Ω/Ωc<1，随着N增加，( Ω/Ωc)2N→0，A(Ω2)→1。

2) 在过渡带和阻带，Ω/Ωc>1，随着N增加，Ω/Ωc>>1，A(Ω)快速下降。 3) Ω=Ωc时， 2.3.2 设计思想

依据低通滤波器和高通滤波器、带通滤波器和带阻滤波器的关系绘制动态图形对比其幅度特性。首先从所给的参数求出所需要的绘制低通滤波器的条件参数，然后通过各滤波器之间的关系表达式画图。由于他们之间的关系是通过横做标的关系来表达的，即通过转换横坐标就可以动态的画图。例如低通与高通之间，低通归一化频率w和高通归一化频率W之间的关系是w=1/W，所以，低通中w点对应于高通中1/W点的值是相同的，由此画图。 2.3.3 设计流程图

，幅度衰减

，相当于3db衰减点。

2

结束 1、低通到高通的变2、低通到带通的变换 3、低通到带阻的变换 4、退出程序 输入数字选择菜单，进行滤波器的转换动态表示 开始 2.3.4 关键技术设计

本实验关键的设计在于如何实现低通于其他的转换，在本设计中是通过改变横坐标的关系来设计的，在具体的设计中，采用到的设计如下： 1．模拟高通滤波器设计：

① 确定高通滤波器的技术指标：通带下限频率Ώp，阻带上限频率Ώs，

通带最大衰减αp，阻带最小衰减αs。

② 确定相应低通滤波器的设计指标：按λ＝1/η，将高通滤波器的边界频

率转换成低通滤波器的边界频率，各相设计指标为： 低通滤波器通带截止频率Ώp＝1/Ώp； 低通滤波器阻带截止频率Ώs＝1/Ώs； 通带最大衰减仍为αp，阻带最小衰减仍为α③ 设计归一化低通滤波器G(p)。

④ 求模拟高通的H(s)。将G(p)按λ＝1/η转换成归一化高通H(q)，为

去归一化，将q=s/Ώc代入H(q)中，得H(s)＝G(p)|p=Ώc/s

2．模拟带通滤波器设计思路：

s

① 确定模拟带通滤波器的设计指标，即：带通上限频率Ώu，带通下限频

率Ώl；下阻带上限频率Ώs1，上阻带下限频率Ώs2

通带中心频率Ώ02=ΏlΏu，通带宽度B=Ώu-Ώl

s1s1/B,s2s2/B与以上边界频率对应的归一化边界频率如下：

ll/B,uu/B2

0lu还需要确定的技术指标有：通带最大衰减αp，阻带最小衰减αs。

② 确定归一化低通技术要求； ③ 设计归一化低通G(p)； ④ 将G(p)转换成带通H(s)。 3．模拟带通阻滤波器设计思路：

① 确定模拟带通滤波器的设计指标，即： 下通带上限频率Ώl，上通带下限频率Ώu

阻带下限频率Ώs1，阻带上限频率Ώs2 阻带中心频率Ώ02=ΏlΏu，通带宽度B=Ώu-Ώl 它们相应的归一化边界频率为： s1s1,s2s2,llBBB u2,0luu ② 归一化模拟低通技术要求； B

③ 设计归一化模拟低通G(p)； ④ 将G(p)转换成带阻滤波器H(s) 2.4 设计切比雪夫I型低通滤波器 2.4.1 基本的原理

误差值在规定的频段上等波纹变化。 巴特沃兹滤波器在通带内幅度特性是单调下降的，如果阶次一定，则在靠近截止

处，幅度下降很多，或者说，

为了使通带内的衰减足够小，需要的阶次N很高，为了克服这一缺点，采用切比雪夫多项式来逼近所希望的

。切比雪夫I型低通滤波器的

在通带范围内是等幅起伏的，所以在同样的通常内衰减要求下，其阶数较巴特沃兹滤波器要小。

2.4.2 设计思想

可以通过直接调用有关切比雪夫I型低通滤波器的库函数，matlab本身就包含有这些函数供用户使用，另一种方法是按照公式一步一步地往下做，但是由于所学知识有限以及时间问题在本设计中采用了第一种方法，调用函数。然后通过定点绘制说明线来表达其特性。 2.4.3 设计流程图

绘制幅度特性图 不继续 结束 选择继续进行 调用库函数，求出N以及过度频率wc，并且求出等差频率所对应的幅度值 开始 输入有关参数fp、fs、ap、as

2.4.4 关键技术设计

由于采用了调用库函数的方法，所以关键的部分就不在于求出其转换的表达式了，重点在于如何更好的表达切比雪夫I型低通滤波器的特性，所以在本设计中田加了说明线，例如在其幅度的增益图象中，说明其通带的部分、过渡带以及阻带部分。

2.5 双线性变换法设计巴特沃斯低通数字滤波器 2.5.1 基本的原理

1.确定技术指标p，s，p，s

2.求滤波器阶数N

Nlgksplgsp

/103．求归一化极点pk

10p1,sps/p 其中 ksp10s/101pke12k1j()22N,k0,1,2,...,N1

再将pk代入式，得到归一化传输函数Ha(p)。

Ha(p)1(ppk0N1

k)

4．将Ha(p)去归一化，得到实际的Ha(s)，即Ha(s)Ha(p)2.5.2 设计思想

ps/p本设计中，由于我是选择先做这个双线性，所以采用了利用公式一步一步地做出来。首先利用双线性变换的公式w=2*tan(0.5*W)/T,其中w表示模拟角频率，而W表示数字角频率，再通过s=j*和双线性变换中s和z的关系，可以得到s=(2*(1-exp(-j.*W)))./(T*(1+exp(-j.*W)))，其中的s是数字滤波器的s，s=jw得到结果。这样利用原理可以求出它的几个极点s，带进公式，通过选取一定范围内的频率就可以求得对应的幅度值，结果就出来了。 2.5.3 设计流程图

绘制图形 结束 根据N和极点已经使用按照公式求出N以及各极点 开始 s=(2(1-exp(-jW)))/(T(1+exp(-jW)))带入，进行频率采样

2.5.4 关键技术设计

由于这个程序要按照公式一步一步地执行，所以关键在于如何实现极点与频率W的转换，后来看到了s=j*W，（这个s是数字滤波器的）这就是关键所在，从参考资料上可以知道，s和z-1的关系，而z-1=e-sT,所以便可以将W和模拟滤波器中的极点s的关系求出来了，这就是本实验的关键所在。有了这个关系之后使用for循环便可以将幅度函数和W的关系求出来了。 2.6用凯塞窗设计高通滤波器 2.6.1 基本的原理

凯塞窗函数的时域形式可表示为

22kI011N1 0kN1 w(k)I0()

其中， 是第1类变形零阶贝塞尔函数， 是窗函数的形状参数，由下式确定：

(8.7),500.1102 (21)0.40.07886(21),2150 0.54820,21

其中， 为凯塞窗函数的主瓣值和旁瓣值之间的差值(dB)。改变β的取值，可以对主瓣宽度和旁瓣衰减进行自由选择。β的值越大，窗函数频谱的旁瓣值就越小，而其主瓣宽度就越宽。 2.6.2 设计思想

按照公式首先要求出a和β，然后求出w（k）,这是前面必须要求出来的，然后通过相关的转换求得结果，其思想与2.5.2基本相同。 2.6.3 关键技术设计

这个设计的关键在于1、如何求得w（k）；2、如何求得实际的冲击响应。如何求得w（k）在设计思想已经大概讨论过了，如何求得冲击响应，可以通过

hn=fir1(N,wc/pi,'high',kaiser(N+1,b))来求得

第三章 程序实现

3.1卷积原理演示 3.1.1 程序的实现

本设计中线性卷积的主要实现主要由如下片段实现：

p=length(x1);q=length(x2);n=p+q-1; a=0:q-1;

y2(a+1)=x2(q-a); for n=1:p+q-1

k=-q+n:1:-1+n;subplot(3,1,2) stem(k,y2)

title('x2(n-m)');axis([-16,16,0,24]);

这一部分是实现翻转并移位，在设计中是将序列x2进行翻转和移位 --------------------------------------------------------------

y=conv(x1,x2); t=1:1:n h(t)=y(t);

subplot(3,1,3) t=0:n-1;

stem(t,h);

title('线性卷积y(n)') axis([-16,16,0,24]); pause(1) end

以上整个部分就是实现线性卷积的过程

--------------------------------------------------------------- 对于循环卷积，要求我们进行判断并根据情况做不同的分析：

p=length(x1);q=length(x2);k=max(p,q);%pN

if px11=[x1,zeros(1,N-p)]; n=0:1:N-1; x22(n+1)=x2(n+1); elseif p==q|p>N

n=0:1:N-1; x11(n+1)=x1(n+1)

x22(n+1)=x2(n+1);

else disp('错误，x1的长度要比x2短') end

由于原程序比较长，在执行分情况讨论做循环卷积时程序很长所以不显示出来。而第三部分，验证时域卷积定理和效率对比，通过观察通过fft实现出来的线性卷积和直接进行线性卷积运行的时间长短来进行效率的对比，程序长度可以通过rand来实现足够大的长度的随机数来进行验证。 3.1.2 结果及分析

1、开始运行程序，会进入主菜单，按照提示进行选择： 请输入x1:[1,1,1,1,0,0,1,1,1,1,0,0] 请输入x2:[0,1,2,1,0,0,0,1,2,1,0,0] 1、动态演示2序列的线性卷积 2、动态演示2序列的圆周卷积 3、验证时域卷机定理 4、退出

请选择菜单项目:1

选择1之后可以见到其动态的过程，其最后的结果如下图1：

x1(m)10.50210x2(m)051015051015x2(n-m)20100-15-10-50线性卷积y(n)20100-15-10-505101551015

图1 线性卷积的结果

2、运行完成之后可以看到子菜单 1、重新演示

2、返回主采单修改2序列的值 3、返回主菜单 请选择，输入：3

1、动态演示2序列的线性卷积 2、动态演示2序列的圆周卷积 3、验证时域卷机定理 4、退出

请选择菜单项目:2

请输入执行循环卷积的数N:10 其最后的结果如下图2所表示：

x1(m)10.50210x2(m)05100510x2(n-m)20100-15-10-50y（n）20100-15-10-505101551015 图2 循环卷积的结果

3、选择返回主菜单并修改x1和x2进行第三步实验，如下运行： 请输入x1:rand(1,10000) 请输入x2:rand(1,10000) 1、动态演示2序列的线性卷积

2、动态演示2序列的圆周卷积 3、验证时域卷机定理 4、退出

请选择菜单项目:3 其结果如下图所示：

由于是采用随机10000个进行设计，所以得到的图在粘贴时有点不稳定，并且在每个图都有说明标题，所以图下面就不标面图号了。通过分析，第一个图出来比第二个图慢，说明通过fft求线性卷积会比直接卷积快，通过后面2个图可以验证时域卷积定理。

直接线性卷积结果y125002000150010005000020004000600080001000012000140001600018000

频域相乘结果的时域y225002000150010005000020004000600080001000012000140001600018000

2.5x 107频域相乘得到的频域幅度特性|Y2|21.510.50020004000600080001000012000140001600018000

2.5x 107直接线性卷积y1的频域幅度特性|Y1|21.510.50020004000600080001000012000140001600018000

3.2 采样定理演示 3.2.1 程序的实现

程序中主要包括时域部分的分析和频域部分的分析，在时域部分主要将时间定位在t=-0.05:0.001:0.05;

将采样时间T=1/fs的应用，就可以实现时域的分析 n=0.1/T+1 h=ones(1,n); stem(t1,h,'^')

title('采样序列h(t)') %这部分是实现采样序列的时域画图 而频谱的分析则首先画出原序列的频谱 x=[0 1 0];t=[-0.1 0 0.1];

plot(t,x)；title('时域信号的频谱Y(jw)') 以下的这一部分就是实现频域采样的程序部分 wc=0.1;ws=(fs/fc)*wc;N=fix(0.5/ws) for i=0:1:N subplot(3,1,2)

stem(-ws*i,1,'^');axis([-0.5,0.5,0,1]);hold on; stem(ws*i,1,'^');axis([-0.5,0.5,0,1]) subplot(3,1,3)

t1=[-ws*i-wc,-ws*i,-ws*i+wc];t2=[ws*i-wc,ws*i,ws*i+wc] x=[0,1,0];

plot(t1,x,t2,x);hold on;axis([-0.5,0.5,0,1]);pause(1) end

3.2.2 结果及分析

1、运行程序，程序进如采样频率输入，如下 请输入采样频率fs:150 其运行结果如下图3和图4

时域信号ya(t)500-0.0510.50-0.0550-0.04-0.03-0.02-0.0100.010.020.030.040.05采样序列h(t)-0.04-0.03-0.02-0.0100.010.020.030.040.05采样后序列ya(t)0-0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05

图3 采样fc=150HZ时的时域部分

时域信号的频谱Y(jw)1Y(jw)0.50-0.510.50-0.510.50-0.5-0.4-0.3-0.2-0.10w0.10.20.30.40.5-0.4-0.3-0.2-0.100.10.20.30.40.5-0.4-0.3-0.2-0.100.10.20.30.40.5 图4 采样fc=150HZ时的频域部分

重新选择输入采样频率

1、选择1改变采样频率fs，继续 2、按其他输入退出程序 请选择：1

请输入采样频率fs:300 其结果运行如图5和图6

时域信号ya(t)500-0.0510.50-0.0550-0.04-0.03-0.02-0.0100.010.020.030.040.05采样序列h(t)-0.04-0.03-0.02-0.0100.010.020.030.040.05采样后序列ya(t)0-0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05 图5 采样fc=300HZ时的时域部分

时域信号的频谱Y(jw)1Y(jw)

0.50-0.510.50-0.510.50-0.5-0.4-0.3-0.2-0.10w0.10.20.30.40.5-0.4-0.3-0.2-0.100.10.20.30.40.5-0.4-0.3-0.2-0.100.10.20.30.40.5

图6 采样fc=300HZ时的频域部分

分析：从图3到图6对比可以知道，当采样频率大于等于2倍的模拟频率时才

能保证采样不混碟，从而正常的恢复到原信号。 3.3模拟滤波器设计演示 3.3.1 程序的实现

在本设计中对每个转换都使用3个情况，即通带、过渡带和阻带，将它们分别讨论分别绘制图形，如，低通到高通m=1时为通带、m=2时过渡带、m=3时阻带部分，现将它们的部分程序列出来：

if m==1

for i=0:1:10 w1=0:0.1:i*0.1; if i==0

w2=3.2; elseif i==1 w2(2)=3.1 elseif i==2 w2(3)=3.0 elseif i==3 w2(4)=2.8 else

w2(i+1)=1/w1(i+1); end

G=1./((j.*w1).^5+3.236.*(j.*w1).^4+5.2361.*(j.*w1).^3+5.2361.*(j.*w1).^2+3.2361.*(j.*w1)+1); G1=abs(G);

subplot(2,1,1)

plot(w1,G1);hold on axis([0,3.5,0,1.2]) subplot(2,1,2) plot(w2,G1); hold on axis([0,3.5,0,1.2]) pause(0.4) end

end

这一部分就是实现完整的低通到高通中的通带部分，为了使其图形更接近资料中的图形，所以在开始的时候故意取了一些固定的值。 -----------------------------------------------------

for i=10:1:20

w1=1:0.1:0.1*i; w2(i-10+1)=1/w1(i-10+1);

这部分就是实现m=2，过度带的主要函数。

-----------------------------------------------------

for i=20:1:30

w1=2:0.1:0.1*i;

w2(i-20+1)=1/w1(i-20+1);

这部分就是实现阻带部分的主要函数。 3.3.2 结果以及分析 运行程序可以看到主菜单： ******主菜单******

1、观察低通与高通的频域特性 2、观察低通与带通的频域特性 3、观察低通与阻带的频域特性 0、退出程序 请选择功能：1

其运行结果如图7所示：

10.80.60.40.2000.511.522.533.510.80.60.40.2000.511.522.533.5

图7 低通到高通的转变结果

图7第一图表示低通、第二个表示高通，蓝色的线表示通带部分，红色表示过渡带，绿色表示的是阻带部分。

10.80.60.40.20-3-2-1012310.80.60.40.20012345678

图8 低通到带通的图形转变

图8第一图表示低通、第二个表示带通，蓝色的线表示阻带部分，红色表示过渡带，绿色表示的是通带部分。

10.80.60.40.20-3-2-1012310.80.60.40.20012345678

图9 低通到带阻的图形转变

图9中第一个图表示的是低通滤波器的幅度特性，第个图是带阻的幅度特性，

其中绿色表示的是通带部分，红色表示是过度代部分，蓝色的是阻带部分，其实阻带在负半轴还有一个跟它一样的波形，只是没有表示出来。 3.4 切比雪夫I型低通滤波器的幅频特性设计 3.4.1 程序的实现

程序中主要的部分就是求Cn，程序如下： if abs(w/wp)<=1 cn=cos(N*acos(w/wp)); else

cn=N*log(w/wp+sqrt((w/wp).^2-1)); cn=(exp(cn)+exp(-cn))/2; end

所用到的库函数有

[N,Wc]=cheb1ord(wp,ws,ap,as,'s'); [b,a]=cheby1(N,ap,Wc,'s'); w=linspace(1,400,100)*2*pi; H=freqs(b,a,w);

这一部分是实现求出阶数N以及截止频率Wc，[b,a]=cheby1(N,ap,Wc,'s')是求得系统函数H(z)的分子和分母的多项式的系数，H=freqs(b,a,w)用于取得数字滤波器的频率所对应的幅度。 3.4.2 结果及分析

运行程序后可以看到菜蛋提示，按照要求输入 请输入有关的参数 fp的值100 fs的值120 ap的值2 as的值50

其结果运行如下图10：

切比雪夫I型低通滤波器增益20log|H(esp(jw))| (dB)0-50-100-15002040100120140160频率(Hz)切比雪夫I型低通滤波器的幅频特性60801802001幅度|H(jw)|1/sqrt(1+e2)0.500200400wp600w(rad)ws80010001200

图10 切比雪夫I型低通滤波器的幅频特性

图10中第一个图是带增益即分贝db的形式求其特性图，图2是幅度即频谱的绝对值的特性。

3.5 双线性变换法设计巴特沃斯低通数字滤波器 3.5.1 程序的实现

这个程序的实现主要是求N和数字系统函数，其各自的主要程序部分如下： 首先求得阶数N

Us=2*tan(0.5*Ws)/T;Up=2*tan(0.5*Wp)/T;Ksp=Us/Up; ksp=sqrt((10^(0.1*ap)-1)/(10^(0.1*as)-1)); N=-log10(ksp)/log10(Ksp);

N=fix(N)+1;

-------------------------------------

这一部分是实现将数字滤波器的频率W与模拟滤波器中的系统函数相互转换得到数字滤波器的系统函数 Uc=Up*(10^(0.1*ap)-1)^(-1/(2*N));

W=linspace(0,1,100)*pi;

s=(2*(1-exp(-j.*W)))./(T*(1+exp(-j.*W))); Ha=Uc^N; for i=0:1:N-1;

sk(i+1)=Uc*exp(j*pi*(0.5+(2*i+1)./(2*N))); Ha=Ha./(s-sk(i+1)); end H=abs(Ha);

3.5.2 结果以及分析

1幅度|H(jw)|0.5000.511.5w/rad22.53增益20log|H(esp(jw))| (dB)0-20-40-60-8011.11.21.31.41.5w/rad1.61.71.81.92

图11 双线性变换法实现的数字低通滤波器的幅度特性

y在平面上的冲击响应0.60.40.20-0.20204060100120ty的绝对值的图示80140160180200y绝对值|y|0.40.30.20.10020406080100t120140160180200

图12 双线性变换法实现的数字低通滤波器的冲击响应

从图11可以看出，这个通过双线性变换法求得的低通滤波器是符合要求的数字滤波器，图12是它的冲击响应，由于冲击响应是虚数的，所以分别做了在平面上和绝对值这2个图示。 3.6利用凯塞窗设计高通滤波器 3.6.1 程序的实现

这个设计可以采用库函数法也可以采用公式法，本设计采用利用公式的做法，其部分重要的程序如下：

以下这部分是求函数中幅度函数w(n) for n=0:63

b=a*sqrt(1-((2*n/63-1))^2);x=1; for k=1:20 s=1; for i=1:k s=s*i;end

x=((b/2)^(k)/s)^(2)+x;end

y=1;for k=1:20 s=1; for i=1:k s=s*i;end

y=((a/2)^(k)/s)^(2)+y;end w(n+1)=x/y;end

----------------------- 这一部分是设计实际单位取样响应函数 wc=(wp+ws)/2;

N=-20*log10(min(1-10^(-ap/20),10^(-as/20))); N=N+rem(N,2);

hn=fir1(N,wc/pi,'high',kaiser(N+1,b)); 3.6.2 结果以及分析

窗函数的幅度|Wn(jw)|特性10.50010203040506070滤波器的幅度|H(jw)|特性500-50-10000.10.20.30.40.50.60.70.80.91 图13 凯塞窗窗函数以及其高通滤波器的幅度特性

理想冲激响应h(n)0.10.050-0.05-0.1010203040506070实际冲激响应h(n)0.10.050-0.05-0.1010203040506070 图14 凯塞窗下的高通滤波器的理想和实际的冲击响应

分析：图13体现了凯塞窗的窗口函数的频域特性以及其对应的高通滤波器的幅度特性，图15展现了其理想和实际中冲击响应的差别，可以看出实际的冲击响应和理想的冲击响应还是有区别的。

第四章 结束语

通过2个星期紧张的课程设计，不用说，收获肯定是不少的。虽然自己在选修课上有学习matlab，但是所做的实验还与数字信号处理所做的实验不同，平时数字信号处理实验中让我感受到的提高已经很充裕了，这次的课程设计更让我受益非浅。

经过总结，在本次的课程设计中我的收获成果主要有：

1、巩固我所学过的数字信号处理课程的有关知识，同时也让我发现了平时掌握不牢固或错误的地方。例如平时很不注意的数字角频率和频率有可能混了，以前对双线性法的原理不大懂，通过这次的实验让我对它的公式进行了研究，从而掌握了这点知识，通过课程设计，让我们平时不大会的从低通到高通、带通、阻带的变换在这次设计中学会了如何去转变，并且掌握了其原理，印象更深刻。

2、让我对matlab这软件有了更深的了解已经它与数字信号处理这门课程之

间的紧密关系，matlab中是采用数组和距阵的方式处理数据，如何将数字信号处理有关的资料以数组和距阵进行编程是我们学习的一个方面，通过这次的课程设计，让我发现了数字信号处理在matlab中的应用，同时也激发了我利用这软件来实现数字信号处理有关问题的兴趣。

3、它提高我检查错误以及调试的能力。在这次课程设计中错误是难免的，有时候2个小时都在想同一个问题，但当自己能够将错误检测出来并且逐渐积累了调试的经验时，愉悦的心情是无法比拟的。

本次课程设计中，有机地结合了理论与实践，既考察了我们对理论知识的掌握情况，还反映出我们实际动手能力和编程能力，更主要的是它激起我们创新思维，提高了自己独立分析问题和解决问题的能力，这就是它的吸引力所在，我会为了能编写出程序好几个小时都在研究如何编程，这在无形中以及提高了我各方面的能力。

时间虽然短暂，但这次课程设计所给的绝不会是短暂的效果，无论是在知识上，还是在思想上都给我烙下了深刻的印象，我想，这次对于这次课程设计我获得的已经满足了，谢谢老师在这段时间的指导。

参考文献

[1] Sanjit K. K. Mitra，Digital Signal Processing: A Computer-Based Approach，

2000

[2] Richard G. Lyons，Understanding Digital Signal Processing，科学出版社 [3]赵树杰等，数字信号处理，西安电子科技大学出版社，1997.10

[4] 丁玉美等，数字信号处理—时域离散随机信号处理，西电出版社，2002.11 [5] 陈怀琛等，MATLAB及在电子信息课程中的应用，电子工业出版社出版，2002.4 [6]李丽 王振领，MATLAB工程计算及应用，人民邮电出版社，2001.9

[7] 陈怀琛等译，《数字信号处理及其MATLAB实现》，电子工业出版社；