数字信号处理实验(吴镇扬版)matlab程序

(1)数组的加、减、乘、除和乘方运算。输入A=口234],B=[3456],求C=A+B,D=A-B,

E=A.*B,F=A./B,G=A/B并用stem语句画出A、B、C、D、E、F、Go

clearall;

a=[l234]：

b=[3456];

c=a+b;

d=a-b;

e=a.*b;

f=a./b;

g=a."b;

n=I:4;

subplot(4,2,1);stem(n,a);

xlabelCn');xlim(.[05])：ylabel("A');

subplot(4,2,2);stem(n,b);

xlabel('n');xlim([05]);ylabel('B');

subplot(4,2,3)：stem(n,o);

xlabelCn');xlim([05]);ylabel('C');

subplot(4,2,4);stem(n,d);

xlabelCn");xlim([05]);ylabel('D');

subplot(4,2,5);stem(n,e);

xlabelCn});xlim([05]);ylabelC;

subplot(4,2,6);stem(n,f);

xlabelCn');xlim(.[05]);ylabel('F');

subplot(4,2,7);stern(n,g);

xlabelCn');xlim([05]);ylabel('G');

(2)用MATLAB实现下列序列:

a)x(n)=0.8n0WnW15

b)x(n)=e(0.2+3j)n0WnW15

c)x(n)=3cos(0.125nn+0.2n)+2sin(0.25nn+0.1n)0Wn/15

d)将c)中的x(n)扩展为以16为周期的函数xl6(n)=x(n+16),绘出四个周期。

e)将c)中的x(n)扩展为以10为周期的函数xl0(n)=x(n+10),绘出四个周期。

clearall;

N=0:15;

xa=0.8."N;

figure;subplot(2,1,1);stem(N,xa);xlabel(*n*);xlim([016]);ylabel(*xa');

xb=exp((0.2+3*j)*N);

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

xlabel(*n*);xlim([016]);ylabel(*xb');figure:

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

subplot(3,1,1);stem(N,xc);xlabel(*n');xlim([016");ylabel(?xc');

k=0:3;m=0;

fori=l:4

forj=l:16

m=m+l;

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

xl6(m)=3*cos(0.125*pi*n(m)+0.2*pi)+2*sin(0.25*pi*n(m)+0.l*pi);

end

end

subplot(3,1,2);stem(n,x16);xlabel('n);ylabel('xl6');

forj=l:10

xl0(j)=xl6(j);

end

fori=l:3

form=1:10

xl0(i*10+m)=xl0(m);

end

end

n=l:40;

subplot(3,1,3);stem(n,xlO);xlabel(*n')jylabel(*xlO');

(3)x(n)=[l,-l,3,5],产生并绘出下列序列的样本:

a)xi(n)=2x(n+2)-x(n-l)-2x(n)

b)x2(n)=^nx(n-k)

k=l

clearall

n=l:4;

T=4;

x=[l-135];

x(5:8)=x(l:4);

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

fori=l:4

ifi-l<0

xl(i)=2*x(i+2)-x(i-l)-2*x(i);

else

xl(i)=2*x(i+2)-x(i-l+T)-2*x(i);

end

end

xl(5:8)=xl(l:4);

subplot(2,1,2);stem(l:8,xl);grid;

(4)绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注:

a)x(t)=sin(2JTt)OWtWIOs

b)x(t)=cos(100nt)sin(nt)OWtW4s

ta=0:0.05:10;

xa=sin(2*pi*ta);

subplot(2,1,1);plot(ta.xa);

xlabelCVJiylabelC幅度')；

tb=O:O.01:4;

xb=cos(100*pi*tb).*sin(pi*tb);

subplot(2,1,2);plot(tb.xb);

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

(5)编写函数stepshift(nO,nl,n2)实现u(n-nO),nl〈n(Kn2,绘出该函数的图形，起点为

nl,终点为n2。

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

n=[ns:nf];

x=[(n-nO)>=0];

stem(n,x);

⑹给一定因果系统H(z)=(1+缶"+1)«1.O67z」+0.9z=)求出并绘制H⑸的幅频响

应与相频响应。

clearall;

b=[l,sqrt(2),1];

a=[l,-0.67,0.9];

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

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

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

ph=angle(h);

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

（7）计算序列｛8-2-123｝和序列（23-1-3｝的离散卷积，并作图表示卷积结果。

clearall;

a=[8-2-123];

b=[23-1-3];

c=conv(a,b);为计算卷积

M=length(c)-1;

n=O:l:M;

stem(n,c);

xlabelCn');ylabel(R^lT);

(8)求以卜.差分方程所描述系统的单位脉冲响应h(n),0WnW50

y(n)+O.1y(n-l)-O.06y(n-2)=x(n)-2x(n-1)

clearall;

N=50;

a=[l-2];

b=[l0.1-0.06];

x=[lzeros(1,N1)];

k=0:1:N-1;

y=filter(a,b,x);

stem(k,y);

xlabel(Jn);ylabelC幅度');

(1)、观察高斯序列的时域和幅频特性，固定信号xa(n)中参数p=8,改变q的

值，使q分别等于2,4,8,观察它们的时域和幅频特性，了解当q取不同值时,

对信号序列的时域幅频特性的影响；固定q=8,改变p,使p分别等于8,13,

14,观察参数p变化而信号序列的时域及幅频特性的影响，观察p等于多少时，

会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出

相应的时域序列和幅频特性曲线。

解：程序见附录程序一：

n=0:l:15;

%p=8不变,q变化(2,4,8);

p=8;q=2;%p=8;q=2;

xal=exp(-((n-p).A2)/q);

subplot(5,2,l);

plot(n,xal,*-**);

xlabel('t/T');

ylabelCxatn),);

title(*p=8q=2')

xkl=abs(fft(xal));

subplot(5,2,2);

stcm(n,xkl)

xlabel('k');

ylabeI('Xa(k),);

title('p=8q=2')

p=8;q=4;%p=8;q』

xa1=exp(-((n-p).A2)/q);

subplot(5,2,3)；

plot(n,xal/-*');

xlabelCt/T,);

ylabel(,xa(n),);

title('p=8q=4')

xkl=abs(fft(xal));

subplot(5,2,4);

stem(n,xkl)

xlabcl('k');

ylabel('Xa(k),);

title('p=8q=4')

p=8;q=8;%p=8;q=8;

xal=exp(-((n-p).A2)/q);

subplot(5,2,5);

plot(n,xal,

xlabel('t/T');

ylabel('xa(n)');

xkl=abs(fft(xal));

title('p=8q=8')

subplot(5,2,6);

stem(n,xkl)

xlabel('k');

ylabel('Xa(k),);

title('p=8q=8')

%q=8不变,p变化(8,13,14);

p=8;q=8;%p=8;q=8;

xa1=exp(-((n-p).A2)/q);

subplot(5,2,5);

plot(n,xal,*-*')；

xlabel('t/T');

ylabcl('xa(n)');

xkl=abs(fft(xal));

title('p=8q=8')

subplot(5,2,6);

stem(n,xkl)

xlabel('k');

ylabel('Xa(k),);

title(*p=8q=8')

p=13;q=8;%p=13;q=8;

xa1=exp(-((n-p).A2)/q);

subplot(5,2,7);

plot(n,xal,

xlabel('t/T');

ylabel(,xa(n)');

xkl=abs(fft(xal));

title('p=13q=8')

subplot(5,2,8);

stem(n,xkl)

xlabel('k');

ylabel('Xa(k)');

title(*p=13q=8')

p=14;q=8;%p=14;q=8;

xa1=exp(-((n-p).A2)/q);

subplot(5,2,9);

plot(n,xal/-*');

xlabelCt/T');

ylabel('xa(n)');

title('p=14q=8f)

xkl=abs(fft(xal));

subplot(5,2,10);

stcm(n,xkl)

xlabel('k');

ylabelCXa(k),);

title(*p=14q=8,)

时域特性、幅频特性

时域特性

t/Tk

k

分析：

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

当P取固定值时，时域图都关于产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,观察这两种情况下，频谱的形状和谱峰出现位置，有无

混叠和泄漏现象？说明产生现象的原因。

“叫。其他

nl=0:l:15;

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

subplot(3,2,l);

plol(nl,xbl,

xlabel('n');

ylabeiCxCn)1);

title('f=0.0625');

xkl=abs(fft(xbl));

subplot。,2,2);

stem(nl,xkl)

xlabel('k');

ylabelCXCk)1);

title(T=0.0625');

n2=0:1:15;

xb2=exp(-0.1*n2).*sin(2*pi*0.4375*n2);

subplot。,2,3);

plot(n2,xb2,*-*');

xlabei('n');

ylabel('x(n)");

titleCf^.4375');

xk2=abs(fft(xb2));

subplot(3,2,4)；

stem(n2,xk2)

xlabel('k');

ylabel('X(k)');

title('f=0.4375');

n3=0:l:15;

xb3=exp(-0.1*n3).*sin(2*pi*0.5625*n3);

subplot(3,2,5);

plot(n3,xb3,

xlabel('n');

ylabelCxCn)1);

title('f=().5625')；

xk3=abs(fft(xb3));

subplot。,2,6);

stem(n3,xk3)

xlabcl('k');

ylabelCXCk)1);

title('仁0.5625);

时域特性:

f=C.O625

X

X

f=C.5625

1

^2v

X

051015

k

分析：

当f二门=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频谱还有

相同之处吗？这些变化说明了什么？

三角波序列：

n,0<//<3

儿.(〃)=，8-z/,4</?<7

0，其他

反三角波序列：

4-/?,0<n<3

儿(〃)=1〃-4A<n<7

0，其他

clearall;

n=fO:3J;k=[l: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);

工1处川心》1处理时域特性);妙〃1,3:三角波);

subplot(2,2,2);plot(k-l,abs(ffi(Xc)));

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

%%%%%%%%%%%反三角波特性%%%%%%%%%%

subplot(2,2,3);piot(k-l,Xd);

HabelC心;ylabelC时域特性);text(3,3：反三角波);

subplot(2,2,4);plot(k-1,abs(ffi(Xd)));

xlabe/kXylabelC幅频特性*exl(4,10,反三角波);

%末尾补0,计算32点FFT

Xc(9:32)=0;Xd(9:32)=0;k=l:32;fiSi(re;

%%%%%%%%%%%三角波特性%%%%%%%%%%%%%

subplot(2,2,1);plot(k-l,Xc);

xlcibeic心;ylabeic时域特性);三角波);

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

xlabeK'kXylabelC幅频特性);text(4,10,'三角波上

%%%%%%%%%%%必〃三留波潜性%%%%%%%%%%

subplot(2,2,3);plot(k-l,Xd);

xlabelCnXylabeU：时域特性);心|(3,3；反三角波);

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

HabelCkXylabelC幅频特性);text(4J0；反三角波上

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

时域特性频域特性

nk

分析：

由图知，三角波序列和反三角波序列的时域图像成镜像关系，但频域

图像完全一样，只是因为幅频图是对x(k)的值取绝对值。

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

13

时域特性频域特性

010203040010203040

nk

分析：

由实验所得的图形如，N=32点时x,(〃)和4(〃)的幅频特性都更加密集，更多

离散点的幅值显示，“栅栏效应”减小，分辨率提高，而对于猫(〃)来说变化更加

明显。在原序列的末端填补零值，变动了DFT的点数，人为的改变了对真实频

谱采样的点数和位置，相当于搬动了“尖桩栅栏”的位置，从而使得频谱的峰点

和谷点暴露出来。N=32时，&(〃)和与5)的频谱差别较大，但总体趋势仍然都

是中间最小，两侧呈对称。

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

x(n)=sin2n*0.125n+cos27i*(0.125+Af)nn=0,1........,N-1

已知N=16,Af分别为1/16和1/64,观察其频谱；当N=128时，Af不变，其结

果有何不同，为什么？

clearall;

%%%%%%%N=/6%%%%%%%%%%%%

N=16;detf=l/16;n=[0:N-lJ：

xl(n+l)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);

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

%%%%%%%N=I6,detf=1/16%%%%%%%%%%%%

subplot(2,2,1);stem(n,xl);hold;plot(n,xl);

xlabeK'nXylabelC时域特性,'N=16,detf=I/16V

subplot(2,2,2);stem(n,abs(fft(xl)));

14

xlabelCn');ylabel('幅值特性);lext(6,4，N=16,detf=1/16');

%%%%%%%N=16,detf=1/64%%%%%%%%%%%%

subplot(2,2,3);stem(ntx2);

xlabelCnXylabe*时域特性);text⑹l.'N=16,detf=J/64');

subplot(2,2,4);stem(n,abs(ffi(x2)));

xlabelCnXylabelC幅值特性);text(64；N=16,detf=1/64'};

%%%%%%%N=128%%%%%%%%%%%%

N=128;detf=///6;n=[():N-l];

x3(n+])=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detf).*n);

detf=1/64;

x4(n+l)=sin(2*pi*0.125.*n)+cos(2*pi*(0.125+detfi.*n);

%%%%%%%N=128,detf=1/16%%%%%%%%%%%%

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

xlabdCnXykibeK'时域特性);axis([O128-22J);text(6,1.5,'N=128,detf=1/16');

subplot(2,2,2);stem(n,abs(fft(x3)));

xlabel(H);ylabelC幅值特性);axis([O128-1070J);text(40,60,'N=128,detf=l/16");

%%%%%%%N=128,detf=1/64%%%%%%%%%%%%

subplot(2,2,3);stem(n,x3);

xlabelCn');ykibel('时域特性);axis([O128-22]);text(6,1.5,'N=128,detf=1/161);

subplot(2,2,4);stem(n,abs(ffi(x4)));

xlabdCnXylabelC幅值特性');axis([0128-1070]);text(40,60,'N=128,detf=l/16');

挈

爨

督

fc

当

婢

出

理

Q5

270

QQQQ

nN=128.detf=1/16

2040GO80100120

N=128,det^1/16

70r

60l-

50l-

40l-

30l-

20l-

10h

0*«

■0

20406080100120

⑸、用FFT分别实现xa(n)(p=8,q=2)和xb(nj(a=0.1,f=0.0625)的16点圆

周卷积和线性卷积。

clearal1;

N=16;

n=0:N-l;

p=8;q=2;

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

a=0.1;f=0.0625;

Xb(n+l)=exp(-a.*n).*sin(2*pi*f.*n);

耨6点循环卷积

Fa=fft(Xa);Fb=fft(Xb);

Fx=Fa.*Fb;

X51=ifft(Fx);

stcm(n,X51);

%16点线性卷积

Xa(N+l:2*N-l)=0;

Xb(N+l:2*N-l)=0;

Fa=fft(Xa);Fb=fft(Xb);

Fc=Fa.*Fb;

X52=ifft(Fc);

figure;stcm(l:2*NT,X52);

(7)用FFT分别计算Xa(n)(p=8,q=2)和xb(n)(a=0.l,f=O.0625)的16点循环相关和线性

相关，问一共有多少种结果，他们之间有何异同点。

c1earal1;

N=16;

n=0:N-l;

p=8;q=2;

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

a=0.1;f=0.0625;

Xb(n+l)=exp(-a.*n).*sin(2*pi*f.*n);

N=length(Xa);

%16点循环相关

Fa=fft(Xa,2*N);Fb=:ft(Xb,2*N);

Fx=conj(Fei).*Fb;

X71=real(ifft(Fx));

X71=[X71(N+2:2*N)X71(1:N)];

n=(-N+l):(N-l);stem(n,X71);

知6点线性相关

Xa(N+l:2*N-l)=0;

Xb(N+l:2*N-l)=0;

Fa=fft(Xa);Fb=fft(Xb);

Fc=conj(Fa).*Fb;

X72=real(ifft(Fc));

figure;stcmd:2*N-1,X72);

(8)用FFT分别计算x.(n)(p二8,q=2)和xMn)(a=0.1,f=0.0625)的自相关函数。

17

clearall;

N=16;

n=O:N-l;

P=8;q=2;

Xa(n+l)=exp(-(n~p).2./q);

a=0.1;f=0.0625;

Xb(n+l)=exp(-a.*n).*sin(2*pi*f.*n);%自然对数的底:e=:2.71828182845904523536

N=length(Xa);

%Xa(n)16点自相关

Fa=fft(Xa,2*N);Fb=fft(Xb,2*N);

Fl二conj(Fa).*Fa;

X81=real(ifft(Fl));

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

n=(-N+l):(N-l);

subplot(2,1,1);stcm(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+l):(N-l);

subplot(2,1,2);stem(n,X82);

18

(1)fc=0.3kHzfb=0.&/8,fr=0.2kHz,A，=2(H8.T=1〃时;设计一切比雪夫高通

滤波器，观察其通带损耗和阻带衰减是否满足要求。

解：

程序：

clear;

fc=3OO;fr=2OO;fs=1000;rp=0.8;rs=20;f=w*fs/(2*pi);

wc=2*fs*tan(2*pi*fc/(2*fs));plot(f,20*log10(abs(h)));

wl=2*fs*tan(2*pi*fr/(2*fs)];axis(r0.fs/2,-80J0]);

lN,wn]=chcb1ord(wc,wt,rp.rs,M);grid;

[B,Al=cheby1(N,rp,wn,'high','s');xlabel(濒率/Hz');

[bz,az]=bilinear(B,A,fs);ylabel，幅度/dB');

[h,vv]=freqz(bz,az);

8

?

例-34

关

-60

50100150200250300350400450500

嫉率「Hz

19

分析:f=200Hz时阻带衰减大于30dB,通过修改axis([0,fs/2,-80,10D为axis([200,fs/2,-lJJ)

发现通带波动rs满足＜0.8。

bz=[0.0262-0.10470.1570-0.10470.0262]

az=[1.00001.52891.65370.94520.2796]

系统函数为：

0.0262-0.10474-0.1570z-2-0.1047z-34-0.0262z-4

H(z)=

I+1.5289z-1+1.6537z-2+0.9452z-3+0.2796z-4

(2)Af=25d&7=bm;分别用脉冲响应不变

f,=02kHz,5=ldB,fr=0.3kHz,

法及双线性变换法设计一巴特沃思数字低通滤波器，观察所设计数字滤波器的幅频特性曲

线，记录带宽和衰减量，检查是否满足要求。比较这两种方法的优缺点。

解：

程序：

clear;figure;plot(f,abs(hl),'-.r',f,abs(h2),'-b');

fs=1(X)();fc=2(X);fr=300;rp=l;rs=25;grid;xlabelC频率/Hz);ylabel。幅度);

%脉冲响应不变法legend(脉冲响应不变法:双线性变换法》

wp=2*pi*fc;title，巴特沃思低通滤波器，线性幅度谱)

ws=2*pi*fr;

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

[blal]=butter(N,wn,'s');

[bzLazl]=impinvar(bl,a1,fs);

[h1,w]=freqz(bzI,az1);

%双线性变换法

wp=2*fis*tan(2*pi*fc/fs/2);

ws=2*fs*tan(2*pi*fr/fs/2);

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

[b2a2]=butler(N.wn,'s');

[bz2,az2]=bilinear(b2,a2,fs);

[h2,w]=freqz(bz2,az2);

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

20

bzl=[0.00000.00020.01530.09950.14440.06110.00750.00020.00000|

azl=[1.0000-1.91992.5324-2.20531.3868-0.63090.2045-0.04500.0060-0.00041

因此脉冲响应不变法的系统函数为：

H______().(XX)2ZT+().OI53Z-2+().0995z7+0.1444z"4—().061-O.OO75ZY-().0002z'___________

,nv，'-1-1.9199r-1+2.5324z-2-2.2053z-?+1.3868r-*-0.6309z-5+0.2O45Z-6-0.0450z-7+0.0060Z-8-0.0C04z^

bz2=[0.01790.10720.26810.35750.26810.10720.0179]

az2=[1.0000-0.60190.9130-0.29890.1501-0.02080.00251

因此双线性变换法的系统函数为：

0.0179+09072z"+0.268lz“+0.3575ZT+0.2681z~410.10722。+0.0179z>

⑶一一]-0.6019Z-!+0.9130z-2-0.2989z-3+0.I501Z-4-0.0208z-5+0.0025z6-

分析：

脉冲响应不变法的N=9,双线性变换法的N=6,由图知它们都满足要求，但脉冲响应

的衰减较快，双线性变换的过渡带窄一些，且阶数比脉冲小，容易实现。

(3)利用双线性变换法分别设计满足下列指标的巴特沃思型、切比雪夫型和椭圆型数字低

通滤波器，并作图验证设计结果：Z=12kHz,6<().5dB,fr=2kHz,

At>40dBJs=8kHz.

解：

程序：

clear;[N.wn]=buttord(wp,ws.rp,rs,'s');

fs=8000;fc=1200;fr=2000;rp=0.5;rs=40;[blal]=butler(N,wn,'s,);

%巴特沃思低通滤波器[bzl,az1J=bilinear(bl,al,fs);

wp=2*fs*tan(2*pi*fc/fs/2);[h1,w]=freqz(bz1,azI);

ws=2*fs*tan(2*pi*fr/fs/2);Hl=20*logl0(abs(hl));

21

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

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

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

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]=chcbyl(N,rp,wn,'low','s');

(bz2,az2]=bilinear(b2,a2,fs);

[h2,w]=freqz(bz2.az2):

H2=20*logl0(abs(h2));

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

figure;

plot(f,H2);

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

grid;

xlabel(濒率/Hz');

ylabelC幅度/dB');

titl贝切比雪夫低通滤波器，对数幅度谱')；

%椭圆型数字低通滤波器

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

ws=2*fs*tan(2*pi*fr/fs/2);

IN,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*logl0(abs(h3));

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

figure;plot(f,H3);

axis([0.fs/2rl00,101);

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

【ill或椭圆型数字低通滤波器，对数幅度谱力

22

巴特沃思低通滤波器:

频率/Hz

bzl=[0.0004O.(X)320.0129().03020.04530.04530.03020.01290.00320.0004]

azl=[1.0000-2.79964.4582-4.54123.2404-1.63300.5780-0.13700.0197-0.0013]

系统函数为:

0.00044-O.OO32Z-1-hO.O129z~2-O.O3O2z~34-0.0453Z-44-0.04532-54-0.0302z-6-»-0.0129z~7+0.00322-84-Q.OJ04z^

-|2-J57-9

i-2.7996z+4,4582z--4.5412z+3.24O4z-*-1.6330z-+0.56802-*-0.l370z-+0.0197z^-0.001Jz

分析：

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

23

切比雪夫低通滤波器：

切比雪夫低通滤波器，对数幅度谱

10e------------1------:------:------:------:------l

cpo

«、

鹫

-70-

-80-

-90•

-100

05001000150020002500300035004000

频率/Hz

bz2=[0.00260.01320.02640.02640.01320.0026J

az2=[1.0000-2.97754.2932-3.51241.6145-0.3334]

系统函数为：

0.0026+0.C132Z-1+0.0264z~24-0.0264z-3+0.0132z-44-0.0026z~5

1-2.9775Z-1+4.2932z-2-3.5124z-3+1.6154z-4-0.3334z~5

分析：

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

24

椭圆型数字低通滤波器:

椭圆型数字低通滤波器，对数幅度谱

10e------------1------:-------:-------:-------:------1~~

-1001-rrrr•rr

05001000150020002500300035004000

频率/Hz

bz=[0.03890.03630.06650.03630.0389]

az=f1.0000-2.14442.3658-1.32500.33321

系统函数为：

“/、0.0389-0.0363Z-1+0.0665z-2-0.0363z-3+0.0389^-4

HU)=---------------------------：----------------------：---------------------；----------------------：—

31-2.1444z-1+2.3658z-2-1.3250^-3+0.3332z-4

分析：

N=4,为四阶椭圆型数字低通滤波器，从图中可以看出通带波动和阻带衰减都满足设计要求。

(4)分别用脉冲响应不变法及双线性变换法设计一巴特沃思数字带通滤波器，已知

fx=30kHz,其等效的膜拟滤波器指标为5<3d8,2kHz<f<3kHz,At>5dB,

f>6kHz,At>2(klB.f<\,5kHz.

解：程序：

clear;[h1,w]=freqz(bzI,azI);

fs=3OOOO:fc=[2OOO.3OOOl;%双线性变换法

fr=[l5OO,6O(M)];rp=3;rs=2();wp=2*fs*tan(2*pi*fc/fs/2);

%脉冲响应不变法ws=2*fs*tan(2*pi*fr/fs/2);

wp=2*pi*fc;[N,wn]=buttord(wp,ws,rp,r

ws=2*pi*fr;[b2,a2]=butler(N,wn,,s');

[N,wn]=buttord(wp,