1、Matlab 在数字信号处理中的应用（基础）一、数据类型：1、整数：Matlab 支持 8 位，16 位，32 位和 64 位的有符号和无符号整数数据类型。如：x=int8(50); %指定 x 的数据类型为 int8.x=502、浮点数：matlab 的默认数据类型是双精度类型（double） ，为了节省存蓄空间，matlab 也支持单精度数据类型的数组。Realmin(single)Ans=1.1755e-038Realmax(double)Ans=2.2251e-3083、复数：matlab 中虚数单位由 i 或者 j 表示。Z=6+7j另一种创建复数的方法可以通过 complex()函
2、数，complex()函数的调用格式：C=complex(a,b)，返回结果 c 为复数，实部是 a，虚部是 b。二、数组的创建1、一维数组的创建：创建一维行向量，只需要把所有数组元素用空格或者逗号分隔，并用方括号吧所有数组元素括起来即可。如用分号，即为列向量。创建等差的一维数组：格式 Var=start-val:step:stop-val。如果步长是 1，可以省略。2、二维数组的创建；在创建二维数组时，用逗号或者空格区分同一行的不同元素，用分号或者软回车区分不同的行。三、函数流程控制1、顺序结构。2、判断语句（if-else if-else-end）.3、循环语句（for-end）四、作图1
3、、二维图：plot(x,y,linespec)，linespec 参数，用于对图像外观属性的控制，包括线条的形状，颜色和点的形状，颜色。stem(x,y);绘制脉冲杆图图形。Stairs(x,y);绘制阶梯图图形。2、图像子窗口：subplot(m,n,p)，将图像分为 mn 个子区域，在第 p 个区域中绘制图像。3、坐标轴：axis(xmin,xmax,ymin,ymax).指定当前图像中 x 轴和 y 轴的范围。4、图形注释：1）标题：title(图形名字)。2）坐标轴名： xlabel(x 轴的名称)，ylabel(y 轴的名称 )。特殊符号的输入 : alpha 的输入，则自动转变成，
4、实验一、几种典型离散时间序列Matlab 中处理的数组，将下标放在变量后面的小括号内，且约定从 1 开始递增。例如：x=5,4,3,2,1,0，表示 x(1)=5, x(2)=4, x(3)=3, x(4)=2, x(5)=1, x(6)=0。要表示一个下标不由 1 开始的数组 x(n)，一般应采用两个矢量，如：n=-3:5;x=1,-1,3,2,0,-2,-1,2,1;这表示一个含有 9 个点的矢量，n 为一组时间矢量，对应x 有：x(-3)=1, x(-2)=-1.。 。 。 。 。 。连续信号作图使用 plot()函数，绘制线性图。离散信号作图使用 stem()函数，绘制脉冲杆图。一些常
5、用的函数：abs():求绝对值（幅值） 。调用格式：y=abs(x)。length():取某一变量的长度（采样点数） 。调用格式：N=length(n)，取 n 的点数，赋值给 N。real()：取一个复数的实部，调用格式：x=real(h)；取复数 h 的实部，赋值给变量 x。imag():取复数的虚部，调用格式：x=imag(h)；取复数 h 的实部，赋值给变量 yx=sawtooth(t);类似于 sin(t)，产生周期为 2pi，幅值从-1 到+1 的锯齿波。x=sowtooth(t,width);产生三角波，其中 width(0=0。Matlab 程序：n=0:49;x=ones(1
6、,50);close all;stem(n,x);title(单位阶跃信号序列) ；3、单位矩形序列：1）x=ones(1,N)，2）利用逻辑关系表达式产生：x=(n-n0=0)subplot(2,2,4), stem(na2,x2,filled); %作离散图形。title(实指数序列(a1);6、复指数序列：x=exp(sigma+jomega)*n);7、矩形波序列：y=rectpuls(t,width).该函数产生一个幅度为 1 宽度为 width，且以 t=0 为对称轴的矩形脉冲信号，width 的默认值为 1.y=square(t,DUTY),产生一个周期为 2*pi,幅值为+1（
7、-1）的周期性方波信号。其中DUTY 表示信号的占空比。默认值为 0.5.例：矩形脉冲信号的波形图：2 (01)MATLAB 程序如下：t=-0.5:0.01:3;t0=0.5;width=1;ft=2*rectpuls(t-t0,width);plot(t,ft);grid on;axis(-0.5,3,0.2,2.2);title(矩形脉冲信号) ；例：产生一个频率为 10HZ，占空比为 30%的周期方波信号。MATLAB 程序如下：t=0:0.001:3;y=square(2*pi*10*t,30);plot(t,y);aixs(0,0.3,-1.2,1.2);title(周期方波信号)
8、 ；实验二、序列的基本运算 1、序列的加法和乘法x=x1+x2; x=x1.*x2例：已知 x1(n)=u(n+2) (-4=0;n2=-5:8;n02=4;x2=(n2-n02)=0;%用 0 值来扩展它们的序列号，变成相同的起点和终点，原来的值不变。n=min(n1,n2):max(n1,n2);N=length(n);y1=zeros(1,N); y2=zeros(1,N);y1(find(n=min(n1) %单位阶跃信号gn=filter(b,a,x2,xi);subplot(1,2,2), stem(n,gn,k);title(系统的单位阶跃响应) ；ylabel(g(n);xla
9、bel(n);axis(0,N-1,-1.1*min(gn),1.1*max(gn);注：1、hold() 控制当前图形是否刷新的双向切换开关。hold on 使当前轴及图形保持而不被刷新，准备接受此后将绘制的新曲线。hold off 使当前轴及图形不再具备刷新的性质。3、 pause() 暂停执行文件，等待用户按任意键继续，pause(n) 在继续执行之前，暂停n 秒。实验四： 离散 LSI 系统的时域响应对于离散 LSI 系统的响应，MATLAB 为我们提供了多种求解方法：（1） 用 conv 子函数进行卷积积分，求任意输入的系统零状态响应。（2） 用 dlsim 子函数求任意输入的系统零
10、状态响应。（3） 用 filter 和 filtic 子函数求任意输入的系统完全响应 。1、 dlsim 子函数功能：求解离散系统的响应。调用格式：y=dlsim(b , a , x),求输入信号为 x 时系统的响应。其中，b 和 a 分别表示系统函数 H（z）中，由对应的分子项和分母项系数构成的数组，而且分母系数要归一化处理。例：（书本 P96 页 15 题）已知一个用以下差分方程表示的线性移不变因果系统为：)(2(cos)1(2) nxyrnry当激励 时，求系统的响应。(uax2、 用 filter 和 filtic 子函数求 LSI 系统对任意输入的响应Filter 子函数功能：对数字
11、系统的输入信号进行滤波处理。调用格式：y=filtet(b , a, x)，对于由 b 和 a 决定的数字系统（b 和 a 分别表示系统函数 H（z）中，由对应的分子项和分母项系数构成的数组，而且分母系数要归一化处理。 ）当输入信号为 x 时，对 x 中的数据进行滤波，结果存于 y 中，长度取max(na , nb).y , zf=filter(b , a, x);除得到结果矢量 y 外，还得到 x 的最终状态矢量 zf。y=filter(b , a, x, zi)；可在 zi 中指定 x 的初始状态。Filtic 子函数功能：为 filter 子函数选择初始条件。调用格式：zi=filtic
12、(b ,a ,y ,x)；求给定输入 x 和 y 时的初始状态。zi=filtic(b , a, y)；求 x=0，给定输入 y 时的初始状态。其中，x 和 y 分别是表示过去的输入和输出。例：已知一个因果系统的差分方程为：6y(n)-2y(n-4)=x(n)-3x(n-2)+3x(n-4)-x(n-6)满足初始条件 y(-1)=0,x(-1)=0，求系统的单位冲激响应和单位阶跃响应，时间轴上N 取 32 点作图。例：已知一个 IIR 数字低通滤波器的系统函数为：=)(zH321047.60439.349.0112zzz输入一个矩形信号序列 x=square(n/5) (-2=0+r(2).*
13、(n+1).*p(2).n.*n-1=0;subplot(1,2,2),stem(n,h);title(用部分分式法求反变换 h(n);%用冲击响应法求 h(n)h1=impz(b,a,N);subplot(1,2,2),stem(n,h1);title(用 impz()求反变换 h(n);例：已知一个离散系统的系统函数 ，输入序列5.0.1)(2zzH1)(zX求系统在变换域的响应 Y(z)及时间域的响应 y(n).MATLAB 程序：syms z;X=z./(z-1);H=z.2./(z.2-1.5*z+0.5);Y=X.*H;y=iztrans(Y);实验六： 离散系统的描述模型及其转换
14、系统传递函数（tf）模型 NMzazabbzXYH210)(系统零极点增益（zpk）模型 )()()(21NMpzzpqqkz极点留数（rpk）模型 01121)( kzrzrzrzHN二次分式（sos）模型lkkkzabbgz1210)(1、tf2zp功能：将系统的传递函数（tf）模型转换为系统函数的零极点增益（zpk）模型调用格式：z,p,k=tf2zp(num, den)；输入系统传递函数模型中分子（ num） 、分母（den） ；多项式的系数向量，求系统函数的零极点增益模型中的零点向量 z,极点向量 p 和增益系数 k.其中 z, p, k.为列量。2、zp2tf功能：将系统函数的零极
15、点增益（zpk）模型转换为系统的传递函数（tf）模型调用格式：num, den=zp2tf(z, p, k)；输入系统函数的零极点增益模型中的零点向量 z,极点向量 p 和增益系数 k.其中 z, p, k.为列量。求系统传递函数模型中分子（ num） 、分母（den） ，多项式的系数向量， 。3、tf2sos功能：将系统的传递函数（tf）模型转换为系统函数的二次分式（sos）模型调用格式：sos, g=tf2sos(num, den)；输入系统传递函数模型中分子（num） 、分母（den） ；多项式的系数向量，求系统函数的二次分式模型的系数矩阵 sos，增益系数 g.4、 sos2tf功能：
16、将系统函数的二次分式（sos）模型转换系统的传递函数（tf）模型调用格式：num, den=sos2tf(sos, g)；输入系统函数的二次分式模型的系数矩阵 sos，增益系数 g（默认值为 1） ，求系统传递函数模型中分子（num） 、分母（den） ，多项式的系数向量。5、sos2zp功能：将系统函数的二次分式（sos）模型转换系统函数的零极点增益（zpk）模型调用格式：z, p ,k=sos2zp(sos, g)；输入系统函数的二次分式模型的系数矩阵 sos，增益系数 g（默认值为 1） ，求系统函数的零极点增益模型中的零点向量 z,极点向量 p 和增益系数 k.其中 z, p, k.为
17、列量。6、zp2sos功能：将系统函数的零极点增益（zpk）模型转换为系统函数的二次分式（sos）模型调用格式：sos, g=zp2sos(z, p, k); 输入系统函数的零极点增益模型中的零点向量 z,极点向量 p 和增益系数 k.其中 z, p, k.为列量。 ，求系统函数的二次分式模型的系数矩阵 sos，增益系数 g.7、residuez 子函数功能：有理多项式的部分分式展开。调用格式：r p c=residuez(b, a)；把 b(z)/a(z)展开成部分分式的形式。b , a=residuez(r p c)；根据部分分式的 r p c 数组，返回有理多项式。其中：b ,a 为按降
18、幂排列的多项式的分子和分母的系数组；r 为余数数组；p 为极点数组；c 为无穷限多项式系数数组。例：已知离散时间系统的传递函数H(z)=10z-1/(1-3z-1+2z-2)求系统的零点向量 z，极点向量 p 和增益系数 k，并列出系统函数的零-极点增益模型。MATLAB 程序：num=0,10,0;den=1,-3,2;z,p,k=tf2zp(num,den);根据结果，零-极点增益模型的系统函数为：H(z)=10(z/(z-2)(z/(z-1).实验七： 离散系统的零极点分析1、zplane功能：显示离散系统的零极点分布图调用格式：zplane(z, p)；绘制由列向量 z 确定的零点，列
19、向量 p 确定的极点构成的零极点分布图。zplane(b, a)；绘制由行向量 b 和 a 构成的系统函数确定的零极点分布图。hz, hp, ht=zplane(z, p)；执行后可得到 3 个句柄向量，hz 为零点线句柄，hp 为极点线句柄，ht 为坐标轴，单位圆及文本对象的句柄。2、roots功能：求多项式的根调用格式：r=roots(a) ；由多项式的分子或分母系数向量求根向量。其中，多项式的分子或分母系数按降幂排列，得到的跟向量为列向量。例（p65）：已知系统的零极点增益模型分别为：， 85.0)(1zH1)(2zH5.1)(3zH求这些系统的零极点分布图以及系统的冲激响应，判断系统的
20、稳定性。结论：当极点处于单位圆内，系统的冲激响应曲线随着频率的增大而收敛；当极点处于单位圆上，系统的冲激响应曲线为等幅震荡；当极点处于单位圆外，系统的冲激响应曲线随着频率的增大而发散。实验八： 离散系统的频率响应一、离散系统频率响应的基本概念已知稳定系统传递函数的零极点增益（zpk）模型为：NnnMmdzcKzH1)()(则系统的频率响应函数为：)(11)()( wjjjwnMmjNnnjwMmjezjw eHeDCKdcKHejw 其中，系统的幅度响应特性为：NnMmjwDCKeH1)(系统的相位频响特性为：)()(11MwMmNn说明：系统函数与频率响应有着紧密的联系，适当地控制系统函数的
21、极点、零点的分布，就可以改变离散系统的频率响应特性。二、离散系统的频率响应子函数 freqz()功能：用于求离散时间系统的频率响应函数 。)(jweH调用格式：1）h, w=freqz(b ,a ,n)。可以得到数字滤波器的 n 点复频响应值，这 n个点均匀地分布之0, 上，并将这 n 个频点的频率记录在 w 中，相应的频响值记录在h 中，n 缺省时取 512 点。2）h f=freqz(b, a, n, Fs)；用于对 在0, Fs/2上等间隔采样 n 点，采样点频率)(jwe及相应频响值分别记录在 f 和 h 中，由用户指定 Fs（以 Hz 为单位）的值。3）h=freqz(b, a, w
22、)；用于对 在0, 上进行采样，采样频率点由矢量 w 指定。)(jH24）h=freqz(b, a, f, Fs)；用于对 在0, Fs上采样，采样频率点由矢量 f 指定。)(jwe2、angle( )功能：求相角。调用格式：p=angle(h)；用于求取复矢量或复矩阵 H 的相角（以弧度为单位） ，相角介于和 之间。3、grid功能：在指定的图形坐标上绘制分格线。调用格式：grid 紧跟在要绘制分格线的绘图指令后面。例如：plot(t, y)；grid。Grid on 绘制分格线。Grid off 不绘制分格线。4、hold功能：在当前轴或图形上多次叠绘多条曲线。调用格式：hold 使当前图
23、形具备刷新性质的双向开关。Hold on 使当前轴或图形保持而不被刷新，准备接受此后将绘制的新曲线。Hold off 使当前轴或图形不再具备不被刷新的性质。5、text功能：在图形上标注文字说明。调用格式：text(xt, yt, string)；在图面上（xt, yt）坐标处书写文字说明，其中文字说明字符串必须使用单引号标注。三、举例 1、已知离散时间系统的系统函数 4321.07.5 21020)( zzzH求系统在 频率范围内，归一化的绝对幅度频率响应，相对幅度频率响应，相位频率响应和零极点分部图。解：MATLAB 程序如下：b=0.2, 0.1, 0.3, 0.1, 0.2;a=1,
24、-1.1, 1.5, -0.7, 0.3;n=(0: 500)*pi/500; %在0, 的范围内取 501 个采样点h, w=freqz(b, a, n); %求系统的频率响应subplot(2, 2, 1), plot(n/pi, abs(h); grid %作系统的绝对幅度频响图axis(0, 1.1*min(abs(h), 1.1*max(abs(h);ylabel(幅度 )；title(幅频响应（V） );subplot(2,2, 2), plot(n/pi, angle(h); grid %作系统的相位频响图axis(0, 1.1*min(angle(h), 1.1*max(ang
25、le(h);ylabel(相位 )；xlabel(以 pi 为单位的频率) ；title(相频响应); db=20*log10(abs(h);subplot(2, 2, 3), plot(n/pi, abs(h); grid %作系统的相对幅度频响图title(幅频响应（db） );subplot(2, 2, 4), zplane(b, a); grid %作零极点分布图title(零极点分布); 四、系统零极点的位置对系统频率响应的影响系统的零极点的位置对系统响应有着非常明显的影响。观察系统极点位置对幅频响应的影响。已知一阶离散系统的传递函数为 ，假设系统的零点 在原点，极1)(pzqH1q
26、点 分别取 0.2,0.5,0.8。比较它们的幅频响应曲线，从中了解系统极点的位置对幅频响应1p有何影响。实验八、 时域抽样与信号的重建对连续信号进行采样：已知一个连续时间信号 f(t)=sin(2pif0)+1/3sin(6pif0), f0=1Hz,取最高有限带宽频率 fm=5f0,分别显示原连续时间信号和 Fs2f0, Fs=2f0, Fsws 时，为高通滤波器，当 wp, ws 为二元向量时，为带通或带阻滤波器，此时Wn 也为二元向量。2、cheb1ord功能：确定切比雪夫 I 型滤波器的阶数和通带截止频率调用格式：n, wn=cheblord(wp, ws, Rp, As )；计算切
27、比雪夫 I 型数字滤波器的阶数和通带截止频率。其中， ，其值为 1 时表示 0.5Fs，Rp 为通带最大衰减指标，As 为阻)(0wsp带最小衰减指标。n, wn= cheblord (wp, ws, Rp, As， s)；计算切比雪夫 I 型模拟滤波器的阶数和通带截止频率。wp, ws,可以是实际的频率值或角频率值， Wn 将取相同的量纲。Rp 为通带最大衰减指标，As 为阻带最小衰减指标。Wpws 时，为高通滤波器，当 wp, ws 为二元向量时，为带通或带阻滤波器，此时 Wn也为二元向量。3、cheb2ord功能：确定切比雪夫 II 型滤波器的阶数和阻带截止频率调用格式：n, wn=ch
28、eb2ord(wp, ws, Rp, As )；计算切比雪夫 II 型数字滤波器的阶数和通带截止频率。其中， ，其值为 1 时表示 0.5Fs，Rp 为通带最大衰减指标，As 为阻)(0wsp带最小衰减指标。n, wn= cheb2ord (wp, ws, Rp, As， s)；计算切比雪夫 II 型模拟滤波器的阶数和通带截止频率。wp, ws,可以是实际的频率值或角频率值， Wn 将取相同的量纲。Rp 为通带最大衰减指标，As 为阻带最小衰减指标。Wpws 时，为高通滤波器，当 wp, ws 为二元向量时，为带通或带阻滤波器，此时 Wn也为二元向量。4、ellipord功能：确定椭圆滤波器的
29、阶数和通带截止频率调用格式：n, wn=ellipord(wp, ws, Rp, As)；计算椭圆数字滤波器的阶数和通带截止频率。其中，其值为 1 时表示 0.5Fs，Rp 为通带最大衰减指标， As 为阻带最小衰减指)(0wsp标。n, wn=ellipord(wp, ws, Rp, As， s)；计算椭圆模拟滤波器的阶数和通带截止频率。wp, ws,可以是实际的频率值或角频率值，Wn 将取相同的量纲。 Rp 为通带最大衰减指标，As 为阻带最小衰减指标。Wpws 时，为高通滤波器，当 wp, ws 为二元向量时，为带通或带阻滤波器，此时 Wn也为二元向量。5、buttap功能：确定巴特沃思
30、（Butterworth）模拟低通滤波器原型调用格式：z, p, k=buttap(n)；设计巴特沃思模拟低通滤波器原型，其传递函数为： )()3()2()1(.10npssps mzzkHa 此时 z 为控阵，巴特沃思滤波器由通带内最平坦，总体上单调的幅度特性来表征。6、cheb1ap功能：切比雪夫 I 型模拟低通滤波器原型调用格式：z, p, k=cheb1ap(n, Rp)；设计切比雪夫 I 型模拟低通滤波器原型，其通带内的波纹系数为 Rp 分贝，传递函数为： )()3()2()1(.10npssps mzzkHa 此时 z 为控阵。切比雪夫 I 型滤波器是通带内等波纹，阻带内单调的滤波
31、器，其极点均为分布在左半平面的椭圆上。7、cheb2ap功能：切比雪夫 II 型模拟低通滤波器原型调用格式：z, p, k=cheb2ap(n, As)；设计切比雪夫 II 型模拟低通滤波器原型，其阻带内的波纹系数小于 As 分贝，传递函数为：)()2()1()( npsspzzksHa 切比雪夫 II 型滤波器是通带内单调，阻带内等波纹的滤波器，其极点位置为 cheb1ap 极点位置的倒数。8、ellipap功能：椭圆模拟低通滤波器原型调用格式：z, p, k=ellipap(n, Rp, As)；设计椭圆模拟低通滤波器原型，其通带内的波纹系数为 Rp分贝，其阻带内的波纹系数小于 As 分贝
32、。其传递函数为：)()2()1()( npsspzzksHa 椭圆滤波器是通带和阻带内均是等波纹的滤波器，它是有比巴特沃斯和切比雪夫更陡的下降斜率，但会损失通带和阻带的波纹指标。9、poly功能：求某根向量所对应的特征多项式调用格式：P=poly( )；求根向量 的特征多项式，产生多项式系数向量。例如：降幂多项式 P(x)= 其系数行向量表达式为：1121 nnaxxa,21naP若要表示( )( )-( )= ，可建立xx1121 nnxx，再利用指令： 。多项式 P 是一个特征多项式， 的元素,21n )(polyP被认为是多项式 P 的根。10、poly2str功能：以习惯方式显示多项式
33、。调用格式：Pa=poly2str(a,s)；以习惯方式显示 s 的多项式。例：输入程序：A=2 ; 6;7 PA=poly(A)PPA=poly2str(PA)得到：PA=1 -15 68 -84PPA=s3 - 15 s2 + 68 s 8411、pzmap功能：显示连续系统的零极点分布图。调用格式：pzmap(b, a)；绘制由行向量 b 和 a 构成的系统的系统函数确定的零极点分布图。例6-4：通过模拟滤波器原型设计一个巴特沃思模拟低通滤波器的系统函数，要求通带截止频率 ，通带最大衰减 ，阻带截止频率skradc/42dB31，阻带最小衰减st820解：程序如下： fp=4000;om
34、gp=2*pi*fp; fs=8000;omgs=2*pi*fs; Rp=3;As=20; n,omgc=buttord(omgp,omgs,Rp,As,s); %计算 n 阶模拟低通滤波器原型，得到左半平面零极点 z0,p0,k0=buttap(4); b0=k0*real(poly(z0); a0=real(poly(p0);%归一化 4 阶巴特沃斯低通滤波器，其分母多项式为 ppa=poly2str(a0,s);ppb=poly2str(b0,s);例6-6 给定模拟低通滤波器的性能指标为，在通带内，即在内，幅度函数的波纹（起伏） ,在阻带内，srad/10204 dB1时，幅度函数衰减
35、 。试求用切贝雪夫滤波器实现时，5. 52所需阶次 N 以及滤波器系统函数 的表达式。)(sHa解：程序如下： % 设计一个切贝雪夫 I 型低通滤波器原型 fp=1000;omgp=2*pi*fp; fs=1500;omgs=2*pi*fs; Rp=1;As=15; n,omgn=cheb1ord(omgp,omgs,Rp,As,s) %计算 n 阶模拟低通原型，得到左半平面零极点 z0,p0,k0=cheb1ap(n,Rp) %归一化 tf 模型的系数 b0=k0*real(poly(z0) a0=real(poly(p0) %习惯性的分母多项式 PPA=poly2str(a0,s);PPA
36、 =s4 + 0.95281 s3 + 1.4539 s2 + 0.74262 s + 0.27563 %求系统的幅频特性 H,omg=freqs(b0,a0); dbH=20*log10(abs(H)+eps)/max(abs(H); plot(omg*omgn/(2*pi),abs(H) %求系统的幅频特性 H,omg=freqs(b0,a0); dbH=20*log10(abs(H)+eps)/max(abs(H);/求分贝值，加 eps 以避开 0 点 plot(omg*omgn/(2*pi),abs(H)0 1000 2000 3000 4000 5000 6000 7000 800
37、0 9000 1000000.10.20.30.40.50.60.70.80.91 plot(omg*omgn/(2*pi),angle(H),grid0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000-4-3-2-101234 plot(omg*omgn/(2*pi),dbH),grid0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000-100-80-60-40-20020例：设计一个 切比雪夫 II 型滤波器原型。技术指标：通带 阻带 。,1,2dBRkHzfppdBAkHzfss4
38、0,5解：程序如下 fp=2000;Omgp=2*pi*fp; fs=5000;Omgs=2*pi*fs; Rp=1;As=40; %计算滤波器的阶数和阻带截止频率 n,Omgc=cheb2ord(Omgp,Omgs,Rp,As,s); %计算 n 阶模拟低通原型，得到左半平面的零极点 z0,p0,k0=cheb2ap(n,As); b0=k0*real(poly(z0); a0=real(poly(p0); %求滤波器分母系数 H,Omg=freqs(b0,a0); dBH=20*log10(abs(H)+eps)/max(abs(H);%化为分贝值 plot(Omg*Omgc/(2*pi)
39、,abs(H),grid n,Omgc=cheb2ord(Omgp,Omgs,Rp,As,s); %计算 n 阶模拟低通原型，得到左半平面的零极点 z0,p0,k0=cheb2ap(n,As); b0=k0*real(poly(z0); a0=real(poly(p0); %求滤波器分母系数 H,Omg=freqs(b0,a0); dBH=20*log10(abs(H)+eps)/max(abs(H);%化为分贝值 plot(Omg*Omgc/(2*pi),abs(H),grid0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 10400.10.20.30.40.50.60.7
40、0.80.91 plot(Omg*Omgc/(2*pi),angle(H),grid0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 104-4-3-2-101234 plot(Omg*Omgc/(2*pi),dBH),grid0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 104-100-80-60-40-20020 %程序运行结果如下： OmgpOmgp =1.2566e+004 OmgsOmgs =3.1416e+004 OmgcOmgc =2.9387e+004 z0z0 =0 + 1.0824i0 - 1.0824i0 + 2.6131i0 - 2.
41、6131i p0p0 =-0.1712 - 0.4761i-0.5045 - 0.2408i-0.5045 + 0.2408i-0.1712 + 0.4761i k0k0 =0.0100 b0b0 =0.0100 0 0.0800 0 0.0800 a0a0 =1.0000 1.3514 0.9139 0.3653 0.0800例：设计一个椭圆滤波器原型技术指标：通带 阻带 。,1,2dBRkHzfppdBAkHzfss40,5解：程序如下实验九（2） 模拟域频率变换法1、lp2lp功能：低通滤波器原型到低通滤波器。这种变换是使用 butter, chebyI, chebyII, ellip
42、函数设计数字低通滤波器的一个步骤。调用格式：bt, at =lp2lp(b, a, W0)。将传递函数表示的截止频率为 1rad/s 的模拟低通滤波器原型变换成截止频率为 W0 的低通滤波器。2、lp2hp功能：低通滤波器原型到高通滤波器。这种变换是使用 butter, chebyI, chebyII, ellip 函数设计数字高通滤波器的一个步骤。调用格式：bt, at =lp2hp(b, a, W0)。将传递函数表示的截止频率为 1rad/s 的模拟低通滤波器原型变换成截止频率为 W0 的高通滤波器。3、lp2bp功能：低通滤波器原型到带通滤波器。这种变换是使用 butter, cheby
43、I, chebyII, ellip 函数设计数字带通滤波器的一个步骤。调用格式：bt, at =lp2bp(b, a, W0,BW)。将传递函数表示的截止频率为 1rad/s 的模拟低通滤波器原型变换成中心频率为 W0，带宽为 BW 的带通滤波器。如果已知被设计的滤波器低端截止频率为 W1，高端截止频率为 W2，则可以计算出W0，BW 的值： W0=sqrt(W1*W2)；BW=W2-W1。4、lp2bs功能：低通滤波器原型到带阻滤波器。这种变换是使用 butter, chebyI, chebyII, ellip 函数设计数字带阻滤波器的一个步骤。调用格式：bt, at =lp2bs(b, a
44、, W0，BW) 。将传递函数表示的截止频率为 1rad/s 的模拟低通滤波器原型变换成中心频率为 W0，带宽为 BW 的带阻滤波器。如果已知被设计的滤波器低端截止频率为 W1，高端截止频率为 W2，则可以计算出W0，BW 的值： W0=sqrt(W1*W2)；BW=W2-W1。5、set功能：设置图形对象属性。 （本实验仅介绍设置坐标刻度的功能）调用格式：set(gca,Xtick,xs,Ytick,ys);二维坐标刻度设置。set(gca,Xtick,xs,Ytick,ys,Ztick,zs);三维坐标刻度设置。xs, ys , zs 可以是任何合法的实数向量，分别决定了 x,y,z 轴的
45、刻度。例：由模拟低通原型设计模拟低通滤波器设计要求：通带截止频率 ，通带最大衰减 ，阻带截止频率Hzfp20dBRp1，阻带最小衰减 。Hzfs50dBAs解：程序如下 fp=2000;Omgp=2*pi*fp; %输入实际滤波器的通带截止频率 fs=5000;Omgs=2*pi*fs; %输入实际滤波器的阻带截止频率 Rp=1;As=20; %输入实际滤波器的通带，阻带截止频率 %计算滤波器的阶数和 3dB 截止频率 n,Omgc=buttord(Omgp,Omgs,Rp,As,s); %计算 n 阶模拟低通原型，得到左半平面零极点 z0,p0,k0=buttap(n); %求归一化的低通滤
46、波器 b0=k0*real(poly(z0); %求归一化的滤波器分子系数 b0 a0=real(poly(p0); %求归一化的滤波器分母系数 a0 H, Omg0=freqs(b0,a0); %求归一化的滤波器频率特性 dbH=20*log10(abs(H)+eps)/max(abs(H); %幅度化为分贝值 %变换为实际模拟低通滤波器 ba,aa=lp2lp(b0,a0,Omgc); %从归一化低通变换到实际低通 Ha,Omga=freqs(ba,aa); %求实际系统的频率特性 dbHa=20*log10(abs(Ha)+eps)/max(abs(Ha); %求实际系统的幅度化为分贝值
47、 %为作图准备数据 Omg0p=fp/Omgc; %通带截止频率归一化 Omg0c=Omgc/2/pi/Omgc; %3 分贝截止频率归一化 Omg0s=fs/Omgc; %阻带截止频率归一化 fc=floor(Omgc/2/pi); %3dB 截止频率 %归一化模拟低通原型频率特性作图 subplot(2,2,1),plot(Omg0/2/pi,dbH); axis(0,1,-50,1);title(归一化模拟低通原型幅度); subplot(2,2,1),plot(Omg0/2/pi,dbH); axis(0,1,-50,1);title(归一化模拟低通原型幅度); set(gca,Xti
48、ck,0,Omg0p,Omg0c,Omg0s,1); set(gca,Ytick,-50,-20,-3,-1);grid subplot(2,2,2),plot(Omg0/2/pi,angle(H)/pi*180); axis(0,1,-200,200);title(归一化模拟低通原型相位); ylabel(phi); set(gca,Xtick,0,Omg0p,Omg0c,Omg0s,1); set(gca,Ytick,-180,-120,0,90,180);grid %实际模拟低通频率特性作图 subplot(2,2,3),plot(Omga/2/pi,dbHa); axis(0,2*fs,-5
