(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,