数字信号处理实验（吴镇扬版）matlab程序

（1）数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。

clear all; a=[1 2 3 4]; b=[3 4 5 6]; c=a+b; d=a-b; e=a.*b; f=a./b; g=a.^b; n=1:4; subplot(4,2,1);stem(n,a); xlabel('n');xlim([0 5]);ylabel('A'); subplot(4,2,2);stem(n,b); xlabel('n');xlim([0 5]);ylabel('B'); subplot(4,2,3);stem(n,c); xlabel('n');xlim([0 5]);ylabel('C'); subplot(4,2,4);stem(n,d); xlabel('n');xlim([0 5]);ylabel('D'); subplot(4,2,5);stem(n,e); xlabel('n');xlim([0 5]);ylabel('E'); subplot(4,2,6);stem(n,f); xlabel('n');xlim([0 5]);ylabel('F'); subplot(4,2,7);stem(n,g); xlabel('n');xlim([0 5]);ylabel('G');

（2）用MATLAB实现下列序列： a) x(n)=0.8n 0≤n≤15

b) x(n)=e(0.2+3j)n 0≤n≤15

c) x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π) 0≤n≤15 d) 将c)中的x(n)扩展为以16为周期的函数x16(n)=x(n+16),绘出四个周期。 e) 将c)中的x(n)扩展为以10为周期的函数x10(n)=x(n+10),绘出四个周期。 clear all; N=0:15;

xa=0.8.^N;

figure;subplot(2,1,1);stem(N,xa); xb=exp((0.2+3*j)*N);

xlabel('n');xlim([0 16]);ylabel('xa');

subplot(2,1,2);stem(N,xb);

xlabel('n');xlim([0 16]);ylabel('xb');figure;

xc=3*cos(0.125*pi*N+0.2*pi)+2*sin(0.25*pi*N+0.1*pi);

subplot(3,1,1);stem(N,xc);xlabel('n');xlim([0 16]);ylabel('xc'); k=0:3;m=0; for i=1:4 for j=1:16

m=m+1;

n(m)=N(j)+16*k(i);

x16(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.1*pi); end

end

subplot(3,1,2);stem(n,x16);xlabel('n');ylabel('x16'); for j=1:10 x10(j)=x16(j); end

for i=1:3

for m=1:10 x10(i*10+m)=x10(m); end end

n=1:40;

subplot(3,1,3);stem(n,x10); xlabel('n');ylabel('x10');

（3）x(n)=[1,-1,3,5],产生并绘出下列序列的样本: a) x1(n)=2x(n+2)-x(n-1)-2x(n)

5b) x2(n)?clear all n=1:4;

?nx(n?k)

k?1T=4;

x=[1 -1 3 5]; x(5:8)=x(1:4);

subplot(2,1,1);stem(1:8,x);grid;

for i=1:4 if i-1<0

x1(i)=2*x(i+2)-x(i-1)-2*x(i); else

x1(i)=2*x(i+2)-x(i-1+T)-2*x(i); end end

x1(5:8)=x1(1:4);

subplot(2,1,2);stem(1:8,x1);grid;

（4）绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注： a) x(t)=sin(2πt) 0≤t≤10s

b) x(t)=cos(100πt)sin(πt) 0≤t≤4s ta=0:0.05:10;

xa=sin(2*pi*ta);

subplot(2,1,1);plot(ta,xa); xlabel('t');ylabel('幅度'); tb=0:0.01:4;

xb=cos(100*pi*tb).*sin(pi*tb); subplot(2,1,2);plot(tb,xb); xlabel('t');ylabel('幅度');

（5）编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

n0=5;ns=1;nf=10;%ns为起点;nf为终点；在=n=n0处生成单位阶跃序列 n=[ns:nf];

x=[(n-n0)>=0]; stem(n,x);

（6）给一定因果系统应与相频响应。 clear all; b=[1,sqrt(2),1]; a=[1,-0.67,0.9];

[h,w]=freqz(b,a);

am=20*log10(abs(h));

subplot(2,1,1);plot(w,am); ph=angle(h);

subplot(2,1,2);plot(w,ph);

H(z)?(1?2z-1?1)/(1-0.67z-1?0.9z-2)求出并绘制H(z)的幅频响

（7）计算序列{8 -2 -1 2 3}和序列{2 3 -1 -3}的离散卷积，并作图表示卷积结果。 clear all;

a=[8 -2 -1 2 3]; b=[2 3 -1 -3];

c=conv(a,b); %计算卷积 M=length(c)-1; n=0:1:M;

stem(n,c);

xlabel('n');ylabel('幅度');

（8）求以下差分方程所描述系统的单位脉冲响应h(n),0≤n≤50 y(n)+0.1y(n-1)-0.06y(n-2)=x(n)-2x(n-1) clear all; N=50;

a=[1 -2];

b=[1 0.1 -0.06]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y);

xlabel('n');ylabel('幅度 ');

(1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8，改变q的值，使q分别等于2，4，8，观察它们的时域和幅频特性，了解当q取不同值时，对信号序列的时域幅频特性的影响；固定q=8，改变p，使p分别等于8，13，14，观察参数p变化对信号序列的时域及幅频特性的影响，观察p等于多少时，会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。 解：程序见附录程序一： n=0:1:15;

%p=8不变，q变化（2,4,8）; p=8;q=2; %p=8;q=2; xa1=exp(-((n-p).^2)/q); subplot(5,2,1); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); title('p=8 q=2') xk1=abs(fft(xa1)); subplot(5,2,2); stem(n,xk1) xlabel('k'); ylabel('Xa(k)'); title('p=8 q=2')

p=8;q=4; %p=8;q=4; xa1=exp(-((n-p).^2)/q); subplot(5,2,3); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); title('p=8 q=4') xk1=abs(fft(xa1)); subplot(5,2,4); stem(n,xk1) xlabel('k'); ylabel('Xa(k)'); title('p=8 q=4')

p=8;q=8; %p=8;q=8; xa1=exp(-((n-p).^2)/q); subplot(5,2,5); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); xk1=abs(fft(xa1)); title('p=8 q=8') subplot(5,2,6);

stem(n,xk1) xlabel('k'); ylabel('Xa(k)'); title('p=8 q=8')

%q=8不变,p变化(8,13,14); p=8;q=8; %p=8;q=8; xa1=exp(-((n-p).^2)/q); subplot(5,2,5); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); xk1=abs(fft(xa1)); title('p=8 q=8') subplot(5,2,6); stem(n,xk1) xlabel('k'); ylabel('Xa(k)'); title('p=8 q=8')

p=13;q=8; %p=13;q=8; xa1=exp(-((n-p).^2)/q); subplot(5,2,7); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); xk1=abs(fft(xa1)); title('p=13 q=8') subplot(5,2,8); stem(n,xk1) xlabel('k'); ylabel('Xa(k)'); title('p=13 q=8')

p=14;q=8; %p=14;q=8; xa1=exp(-((n-p).^2)/q); subplot(5,2,9); plot(n,xa1,'-*'); xlabel('t/T'); ylabel('xa(n)'); title('p=14 q=8') xk1=abs(fft(xa1)); subplot(5,2,10); stem(n,xk1) xlabel('k'); ylabel('Xa(k)');

title('p=14 q=8’)

1时域特性p=8 q=2幅频特性4p=8 q=20.500510t/Tp=8 q=410.500510t/Tp=8 q=810.5005t/T10151515Xa(k)xa(n)2005kp=8 q=441015Xa(k)xa(n)2005kp=8 q=8101015Xa(k)xa(n)5005k1015

时域特性p=8 q=810.5005101510幅频特性p=8 q=8Xa(k)xa(n)50051015t/Tp=13 q=810.500510155kp=13 q=8Xa(k)00xa(n)51015t/Tp=14 q=810.5005t/T10154kp=14 q=8Xa(k)xa(n)2005k1015

分析：

由高斯序列表达式知n=p为期对称轴；

当p取固定值时，时域图都关于n=8对称截取长度为周期的整数倍，没有发

生明显的泄漏现象；但存在混叠，当q由2增加至8过程中，时域图形变化越来越平缓，中间包络越来越大，可能函数周期开始增加，频率降低，渐渐小于fs/2,混叠减弱；

当q值固定不变，p变化时，时域对称中轴右移，截取的时域长度渐渐地不再是周期的整数倍，开始无法代表一个周期，泄漏现象也来越明显，因而图形越来越偏离真实值，

p=14时的泄漏现象最为明显，混叠可能也随之出现；

(2)、观察衰减正弦序列xb(n)的时域和幅频特性，a=0.1，f=0.0625，检查谱峰出现位置是否正确，注意频谱的形状，绘出幅频特性曲线，改变f，使f分别等于0.4375和0.5625，观察这两种情况下，频谱的形状和谱峰出现位置，有无混叠和泄漏现象？说明产生现象的原因。

?e?ansin?2?fn?,0?n?15xb(n)??0其他?

n1=0:1:15;

xb1=exp(-0.1*n1).*sin(2*pi*0.0625*n1);

subplot(3,2,1); plot(n1,xb1,'-*'); xlabel('n'); ylabel('x(n)'); title('f=0.0625'); xk1=abs(fft(xb1)); subplot(3,2,2); stem(n1,xk1) xlabel('k'); ylabel('X(k)'); title('f=0.0625'); n2=0:1:15;

xb2=exp(-0.1*n2).*sin(2*pi*0.4375*n2); subplot(3,2,3); plot(n2,xb2,'-*'); xlabel('n'); ylabel('x(n)'); title('f=0.4375'); xk2=abs(fft(xb2)); subplot(3,2,4); stem(n2,xk2) xlabel('k'); ylabel('X(k)'); title('f=0.4375'); n3=0:1:15;

xb3=exp(-0.1*n3).*sin(2*pi*0.5625*n3); subplot(3,2,5); plot(n3,xb3,'-*'); xlabel('n'); ylabel('x(n)'); title('f=0.5625'); xk3=abs(fft(xb3)); subplot(3,2,6); stem(n3,xk3) xlabel('k'); ylabel('X(k)'); title('f=0.5625');

时域特性f=0.062510-10510nf=0.437510-10510nf=0.562510-105n10154154154幅频特性f=0.0625X(k)x(n)200510kf=0.437515X(k)x(n)200510kf=0.562515X(k)x(n)2005k1015分析：

当f=f1=0.0625时，谱峰位置出现正确，存在在混叠现象，时域采样为一周期，不满足采样定理。

当f=0.4375和0.5625时，时域图像关于Y轴对称，频域完全相同。这是因为频域图是取绝对值的结果，所以完全相同。另外由于时域采样为6个半周期，满足采样定理，无混叠；但由于截取长度不是周期整数倍，出现泄漏。

(3)、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列xc(n)和xd(n)的幅频特性，观察两者的序列形状和频谱曲线有什么异同？绘出两序列及其幅频特性曲线。在xc(n)和xd(n)末尾补零，用N=16点FFT分析这两个信号的幅频特性，观察幅频特性发生了什么变化？两情况的FFT频谱还有相同之处吗？这些变化说明了什么？

三角波序列:

,0?n?3?n?xc(n)??8?n,4?n?7

?0,其他?反三角波序列：

?4-n?xc(n)??n-4?0?,0?n?3,4?n?7 ,其他clear all;

n=[0:3];k=[1:8]; %定义三角波序列 Xc(n+1) = n;Xc(n+5) =4-n; %定义反三角波序列

Xd(n+1) = 4-n;Xd(n+5) =n;

%%%%%%%%%%% 三角波特性 %%%%%%%%%%%%% subplot(2,2,1);plot(k-1,Xc);

xlabel('n');ylabel('时域特性');text(1,3,'三角波'); subplot(2,2,2);plot(k-1,abs(fft(Xc)));

xlabel('k');ylabel('幅频特性');text(4,10,'三角波');

%%%%%%%%%%% 反三角波特性 %%%%%%%%%% subplot(2,2,3);plot(k-1,Xd);

xlabel('n');ylabel('时域特性');text(3,3,'反三角波'); subplot(2,2,4);plot(k-1,abs(fft(Xd)));

xlabel('k');ylabel('幅频特性');text(4,10,'反三角波'); %末尾补0，计算32点FFT

Xc(9:32)=0;Xd(9:32)=0;k=1:32;figure;

%%%%%%%%%%% 三角波特性 %%%%%%%%%%%%% subplot(2,2,1);plot(k-1,Xc);

xlabel('n');ylabel('时域特性');text(1,3,'三角波');

subplot(2,2,2);plot(k-1,abs(fft(Xc)));

xlabel('k');ylabel('幅频特性');text(4,10,'三角波');

%%%%%%%%%%% fan三角波特性 %%%%%%%%%% subplot(2,2,3);plot(k-1,Xd);

xlabel('n');ylabel('时域特性');text(3,3,'反三角波'); subplot(2,2,4);plot(k-1,abs(fft(Xd)));

xlabel('k');ylabel('幅频特性');text(4,10,'反三角波');

N=8时域和幅度频谱图：

时域特性43210024n时域特性682015频域特性Xc(k)xc(n)1050024k频域特性6843210024n682015Xd(k)xd(n)1050024k68

分析：

由图知，三角波序列和反三角波序列的时域图像成镜像关系，但频域图像完全一样，只是因为幅频图是对x（k）的值取绝对值。 N=32时域和幅度频谱图：

13

时域特性432100102030n时域特性402015频域特性Xc(k)xc(n)10500102030k频域特性404321001020n30402010Xd(k)0-100xd(n)1020k3040

分析：

由实验所得的图形知，N=32点时xc(n)和xd(n)的幅频特性都更加密集，更多离散点的幅值显示，“栅栏效应”减小，分辨率提高，而对于xd(n)来说变化更加明显。在原序列的末端填补零值，变动了DFT的点数，人为的改变了对真实频谱采样的点数和位置，相当于搬动了“尖桩栅栏”的位置，从而使得频谱的峰点和谷点暴露出来。N=32时，xc(n)和xd(n)的频谱差别较大，但总体趋势仍然都是中间最小，两侧呈对称。

(4)、一个连续信号含两个频率分量，经采样得

x(n)=sin2π*0.125n+cos2π*(0.125+Δf)n n=0,1……,N-1

已知N=16，Δf分别为1/16和1/64，观察其频谱；当N=128时，Δf不变，其结果有何不同，为什么？

clear all;

%%%%%%% N = 16 %%%%%%%%%%%%

N=16;detf=1/16;n=[0:N-1];

x1(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);

detf = 1/64;x2(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n); %%%%%%% N = 16,detf = 1/16 %%%%%%%%%%%% subplot(2,2,1);stem(n,x1);hold;plot(n,x1);

xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/16'); subplot(2,2,2);stem(n,abs(fft(x1)));

14

xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/16'); %%%%%%% N = 16,detf = 1/64 %%%%%%%%%%%% subplot(2,2,3);stem(n,x2);

xlabel('n');ylabel('时域特性');text(6,1,'N=16,detf=1/64'); subplot(2,2,4);stem(n,abs(fft(x2)));

xlabel('n');ylabel('幅值特性');text(6,4,'N=16,detf=1/64'); %%%%%%% N = 128 %%%%%%%%%%%% N=128;detf=1/16;n=[0:N-1];

x3(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n); detf = 1/64;

x4(n+1)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n); %%%%%%% N = 128,detf = 1/16 %%%%%%%%%%%%

figure;subplot(2,2,1);stem(n,x3);

xlabel('n');ylabel('时域特性');axis([0 128 -2 2]);text(6,1.5,'N=128,detf=1/16'); subplot(2,2,2);stem(n,abs(fft(x3)));

xlabel('n');ylabel('幅值特性');axis([0 128 -10 70]);text(40,60,'N=128,detf=1/16'); %%%%%%% N = 128,detf = 1/64 %%%%%%%%%%%% subplot(2,2,3);stem(n,x3);

xlabel('n');ylabel('时域特性');axis([0 128 -2 2]);text(6,1.5,'N=128,detf=1/16'); subplot(2,2,4);stem(n,abs(fft(x4)));

xlabel('n');ylabel('幅值特性');axis([0 128 -10 70]);text(40,60,'N=128,detf=1/16');

15

(5)、用FFT分别实现xa(n)（p＝8，q＝2）和 xb(n)（a＝0.1，f＝0.0625）的16点圆周卷积和线性卷积。 clear all; N=16; n=0:N-1;

p=8;q=2;

Xa(n+1)=exp(-(n-p).^2./q); a=0.1;f=0.0625;

16

Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); 点循环卷积

Fa=fft(Xa); Fb=fft(Xb); Fx=Fa.*Fb; X51=ifft(Fx); stem(n,X51); 点线性卷积 Xa(N+1:2*N-1)=0; Xb(N+1:2*N-1)=0; Fa=fft(Xa); Fb=fft(Xb); Fc=Fa.*Fb; X52=ifft(Fc);

figure;stem(1:2*N-1,X52);

（7）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的16点循环相关和线性相关，问一共有多少种结果，他们之间有何异同点。 clear all; N=16; n=0:N-1; p=8;q=2;

Xa(n+1)=exp(-(n-p).^2./q); a=0.1;f=0.0625;

Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); N=length(Xa);

点循环相关

Fa=fft(Xa,2*N); Fb=fft(Xb,2*N); Fx=conj(Fa).*Fb; X71=real(ifft(Fx));

X71=[X71(N+2:2*N) X71(1:N)]; n=(-N+1):(N-1); stem(n,X71); 点线性相关 Xa(N+1:2*N-1)=0; Xb(N+1:2*N-1)=0; Fa=fft(Xa); Fb=fft(Xb); Fc=conj(Fa).*Fb; X72=real(ifft(Fc)); figure;stem(1:2*N-1,X72);

（8）用FFT分别计算xa(n)(p=8,q=2)和xb(n)(a=0.1,f=0.0625)的自相关函数。

17

clear all;

N=16; n=0:N-1; p=8;q=2;

Xa(n+1)=exp(-(n-p).^2./q);

a=0.1;f=0.0625;

Xb(n+1)=exp(-a.*n).*sin(2*pi*f.*n); %自然对数的底:e=:2.71828 18284 59045 23536 N=length(Xa);

% Xa(n) 16点自相关

Fa=fft(Xa,2*N); Fb=fft(Xb,2*N); F1=conj(Fa).*Fa; X81=real(ifft(F1));

X81=[X81(N+2:2*N) X81(1:N)]; n=(-N+1):(N-1);

subplot(2,1,1);stem(n,X81); xlabel('n'); ylabel('幅度'); % Xb(n) 16点自相关 Fb=fft(Xb,2*N); F2=conj(Fb).*Fb;

X82=real(ifft(F2));

X82=[X82(N+2:2*N) X82(1:N)]; % n=(-N+1):(N-1);

subplot(2,1,2);stem(n,X82); xlabel('n'); ylabel('幅度');

18

(1)fc?0.3kHz，??0.8dB，fr?0.2kHz，At?20dB,T?1ms;设计一切比雪夫高通滤波器，观察其通带损耗和阻带衰减是否满足要求。 解： 程序： clear;

fc=300;fr=200;fs=1000;rp=0.8;rs=20; wc=2*fs*tan(2*pi*fc/(2*fs)); wt=2*fs*tan(2*pi*fr/(2*fs)); [N,wn]=cheb1ord(wc,wt,rp,rs,'s'); [B,A]=cheby1(N,rp,wn,'high','s'); [bz,az]=bilinear(B,A,fs); [h,w]=freqz(bz,az);

100-10-20幅度/dBf=w*fs/(2*pi);

plot(f,20*log10(abs(h))); axis([0,fs/2,-80,10]); grid;

xlabel('频率/Hz'); ylabel('幅度/dB');

-30-40-50-60-70-80050100150200250300频率/Hz350400450500 19

分析：f=200Hz时阻带衰减大于30dB，通过修改axis([0,fs/2,-80,10])为axis([200,fs/2,-1,1]) 发现通带波动rs满足<0.8。

bz =[0.0262 -0.1047 0.1570 -0.1047 0.0262] az =[1.0000 1.5289 1.6537 0.9452 0.2796] 系统函数为：

H(z)?0.0262-0.1047z1?1.5289z?1?1?0.1570z?2?2?0.1047z?3?3?0.0262z?4?4?1.6537z?0.9452z?0.2796z

（2）fc?0.2kHz，??1dB，fr?0.3kHz，At?25dB,T?1ms;分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。 解： 程序：

clear;

fs=1000;fc=200;fr=300;rp=1;rs=25; %脉冲响应不变法 wp=2*pi*fc; ws=2*pi*fr;

[N, wn] = buttord(wp, ws, rp, rs, 's'); [b1 a1]=butter(N,wn,'s'); [bz1,az1]=impinvar(b1,a1,fs); [h1,w]=freqz(bz1,az1); %双线性变换法

wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2); [N, wn] = buttord(wp, ws, rp, rs, 's'); [b2 a2]=butter(N,wn,'s'); [bz2,az2]=bilinear(b2,a2,fs); [h2,w]=freqz(bz2,az2); f=w/(2*pi)*fs;

20

figure; plot(f,abs(h1),'-.r',f,abs(h2),'-b'); grid; xlabel('频率/Hz'); ylabel('幅度');

legend('脉冲响应不变法','双线性变换法'); title('巴特沃思低通滤波器，线性幅度谱');

巴特沃思低通滤波器，线性幅度谱1.4脉冲响应不变法双线性变换法 1.210.8幅度0.60.40.20 050100150200250300频率/Hz350400450500

bz1 =[0.0000 0.0002 0.0153 0.0995 0.1444 0.0611 0.0075 0.0002 0.0000 0] az1 =[1.0000 -1.9199 2.5324 -2.2053 1.3868 -0.6309 0.2045 -0.0450 0.0060 -0.0004] 因此脉冲响应不变法的系统函数为:

Himp(z)?0.0002z1?1.9199z?1?1?0.0153z?2?0.0995z?3?3?0.1444z?4?0.0611z?5?0.0075z?6?0.0002z?7?8

?0.0004z?9?2.5324z?2?2.2053z?1.3868z?4?0.6309z?5?0.2045z?6?0.0450z?7?0.0060zbz2 =[0.0179 0.1072 0.2681 0.3575 0.2681 0.1072 0.0179] az2 =[1.0000 -0.6019 0.9130 -0.2989 0.1501 -0.0208 0.0025] 因此双线性变换法的系统函数为:

Hbil(z)?0.0179?0.1072z1?0.6019z?1?1?0.2681z?2?2?0.3575z?3?3?0.2681z?4?4?0.1072z?5?5?0.0179z?6?6?0.9130z?0.2989z?0.1501z?0.0208z?0.0025z

分析：

脉冲响应不变法的N=9，双线性变换法的N=6，由图知它们都满足要求，但脉冲响应的衰减较快，双线性变换的过渡带窄一些，且阶数比脉冲小，容易实现。

（3）利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型和椭圆型数字低通滤波器，并作图验证设计结果：fc?1.2kHz，??0.5dB，fr?2kHz，

At?40dB,fs?8kHz。

解： 程序： clear;

fs=8000;fc=1200;fr=2000;rp=0.5;rs=40; %巴特沃思低通滤波器 wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2);

21

[N, wn] = buttord(wp, ws, rp, rs, 's'); [b1 a1]=butter(N,wn,'s'); [bz1,az1]=bilinear(b1,a1,fs); [h1,w]=freqz(bz1,az1); H1=20*log10(abs(h1));

f=w/(2*pi)*fs;

figure; plot(f,H1);%对数幅度谱

axis([0,fs/2,-100,10]);

grid; xlabel('频率/Hz'); ylabel('幅度');

title('巴特沃思低通滤波器，对数幅度谱'); %切比雪夫低通滤波器 wc=2*fs*tan(2*pi*fc/(2*fs)); wt=2*fs*tan(2*pi*fr/(2*fs));

[N,wn]=cheb1ord(wc,wt,rp,rs,'s'); [b2,a2]=cheby1(N,rp,wn,'low','s'); [bz2,az2]=bilinear(b2,a2,fs); [h2,w]=freqz(bz2,az2); H2=20*log10(abs(h2)); f=w*fs/(2*pi); figure; plot(f,H2);

axis([0,fs/2,-100,10]); grid;

xlabel('频率/Hz');

ylabel('幅度/dB');

title('切比雪夫低通滤波器，对数幅度谱'); %椭圆型数字低通滤波器

wp=2*fs*tan(2*pi*fc/fs/2); %双线性变换法 ws=2*fs*tan(2*pi*fr/fs/2);

[N,wp]=ellipord(wp,ws,rp,rs,'s'); [b3,a3]=ellip(N,rp,rs,wp,'low','s'); [bz3,az3]=bilinear(b3,a3,fs); [h3,w]=freqz(bz3,az3); H3=20*log10(abs(h3)); f=w/(2*pi)*fs;

figure;plot(f,H3); axis([0,fs/2,-100,10]);

grid; xlabel('频率/Hz'); ylabel('幅度/dB');

title('椭圆型数字低通滤波器，对数幅度谱');

22

巴特沃思低通滤波器:

巴特沃思低通滤波器，对数幅度谱100-10-20-30幅度-40-50-60-70-80-90-1000500100015002000频率/Hz2500300035004000

bz1 =[0.0004 0.0032 0.0129 0.0302 0.0453 0.0453 0.0302 0.0129 0.0032 0.0004] az1 =[1.0000 -2.7996 4.4582 -4.5412 3.2404 -1.6330 0.5780 -0.1370 0.0197 -0.0013] 系统函数为:

H1(z)?0.0004?0.0032z1?2.7996z?1?1?0.0129z?2?2?0.0302z?3?3?0.0453z?4?4?0.0453z?5?5?0.0302z?6?6?0.0129z?7?7?0.0032z?8?8?0.0004z?9?9?4.4582z?4.5412z?3.2404z?1.6330z?0.5680z?0.1370z?0.0197z?0.0013z

分析：

N=9，为九阶巴特沃思低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。

23

切比雪夫低通滤波器:

切比雪夫低通滤波器，对数幅度谱100-10-20-30幅度/dB-40-50-60-70-80-90-1000500100015002000频率/Hz2500300035004000

bz2=[0.0026 0.0132 0.0264 0.0264 0.0132 0.0026] az2= [1.0000 -2.9775 4.2932 -3.5124 1.6145 -0.3334] 系统函数为： H2(z)?0.0026?0.0132z1?2.9775z?1?1?0.0264z?2?2?0.0264z?3?3?0.0132z?4?4?0.0026z?5?5?4.2932z?3.5124z?1.6154z?0.3334z

分析：

N=5，为五阶切比雪夫低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。

24

椭圆型数字低通滤波器:

椭圆型数字低通滤波器，对数幅度谱100-10-20-30幅度/dB-40-50-60-70-80-90-1000500100015002000频率/Hz2500300035004000

bz =[0.0389 0.0363 0.0665 0.0363 0.0389] az =[1.0000 -2.1444 2.3658 -1.3250 0.3332] 系统函数为：

H3(z)?0.0389-0.0363z1-2.1444z?1?1?0.0665z?2?2?0.0363z?3?3?0.0389z?4?4?2.3658z-1.3250z?0.3332z

分析： N=4，为四阶椭圆型数字低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。 （4）分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字带通滤波器，已知

fs?30kHz，其等效的模拟滤波器指标为??3dB，2kHz?f?3kHz，At?5dB，

f?6kHz，At?20dB，f?1.5kHz。

解：程序： clear;

fs=30000;fc=[2000,3000]; fr=[1500,6000];rp=3;rs=20; %脉冲响应不变法 wp=2*pi*fc; ws=2*pi*fr;

[N,wn]=buttord(wp, ws, rp, rs, 's'); [b1,a1]=butter(N,wn,'s');

[bz1,az1]=impinvar(b1,a1,fs);

25

[h1,w]=freqz(bz1,az1); %双线性变换法

wp=2*fs*tan(2*pi*fc/fs/2); ws=2*fs*tan(2*pi*fr/fs/2); [N,wn]=buttord(wp, ws, rp, rs, 's'); [b2,a2]=butter(N,wn,'s'); [bz2,az2]=bilinear(b2,a2,fs); [h2,w]=freqz(bz2,az2); f=w/(2*pi)*fs;

本文来源：https://www.bwwdw.com/article/6rs8.html