離散信號(hào)和系統(tǒng)分析的MATLAB實(shí)現(xiàn)

1、1.9 離散信號(hào)和系統(tǒng)分析的MATLAB實(shí)現(xiàn)1.9.1 利用MATLAB產(chǎn)生離散信號(hào) 用MATLAB表示一離散序列xk時(shí)，可用兩個(gè)向量來表示。其中一個(gè)向量表示自變量k的取值范圍，另一個(gè)向量表示序列xk的值。例如序列xk=2,1,1,-1,3,0,2可用MATLAB表示為 K=-2:4;x=2,1,1,-1,3,0,2可用stem(k,f)畫出序列波形。當(dāng)序列是從k=0開始時(shí)，可以只用一個(gè)向量x來表示序列。由于計(jì)算機(jī)內(nèi)寸的限制，MATLAB無法表示一個(gè)無窮長(zhǎng)的序列。例1-38 利用MATLAB計(jì)算單位脈沖序范圍內(nèi)各點(diǎn)的取值。解： %progran 1_1 產(chǎn)生單位脈沖序列 Ks=-4;ke=4;

2、n=2;K=ks:ke;X=(k-n)=0;Stem(k,x):xlabel(k);程序產(chǎn)生的序列波形如圖1-49所示。例1-39 利用MATLAB畫出信號(hào)Xk=10sin(0.02)+nk, 的波形。其中nk表示為均值為0方差為1的Gauss分布隨機(jī)信號(hào)。解：MALAB提供了兩個(gè)產(chǎn)生（偽）隨機(jī)序列的函數(shù)。Rand(1,N)產(chǎn)生1行N列的0，1均勻分布隨機(jī)數(shù)。Randn(1,N)產(chǎn)生1行N列均值為0方差為1的Gauss分布隨機(jī)數(shù)。%program 1_2 產(chǎn)生受噪聲干擾的正弦信號(hào)N=100;k=0:N;X=10*sin(0.02*pi*k)+randn(1,N+1);Plot(k,x);Xla

3、bel(k);Ylabel(xk);程序產(chǎn)生序列如圖1-50所示。 離散卷積是數(shù)字信號(hào)處理中的一個(gè)基本運(yùn)算，MTLAB提供的計(jì)算兩個(gè)離散序列卷積的函數(shù)是conv，其調(diào)用方式為 y=conv(x,h)其中調(diào)用參數(shù)x,h為卷積運(yùn)算所需的兩個(gè)序列，返回值y是卷積結(jié)果。MATLAB函數(shù)conv的返回值y中只有卷積的結(jié)果，沒有y的取值范圍。由離散序列卷積的性質(zhì)可知，當(dāng)序列x和h的起始點(diǎn)都為k=0時(shí)，y的取值范圍為k=0至length(x)+length(h)-2。當(dāng)序列x或（和）h的起始點(diǎn)不是k=0時(shí)，由例1-3知，y的非零值范圍可根據(jù)例1-4的結(jié)論進(jìn)行計(jì)算。例1-40試用MATLAB函數(shù)conv計(jì)算例

4、1-2中序列的卷積。解：program 1_3 計(jì)算離散卷積x=-0.5,0,0.5,1; %序列x的值kx=-1:2; %序列x的取值范圍h=1,1,1;kh=-2:0;y=conv(x,h); %計(jì)算卷積k=kx(1)+kh(1):kx(end)+kh(end); %計(jì)算y的取值范圍stem(k,y);xlabel(k);ylabel(y);程序的運(yùn)行結(jié)果如圖1-51所示。 離散LTI系統(tǒng)響應(yīng)MATLAB求解許多離散LTI都可用如下的線性常系數(shù)的差分方程描述 (1-151)其中xk、yk分別系統(tǒng)的輸入和輸出。在已知差分方程的N個(gè)初始狀態(tài)yk,和輸入xk,就可由下式迭代計(jì)算出系統(tǒng)的輸出。 y

5、k=-利用MATLAB提供的filter函數(shù)，可方便地計(jì)算出上述差分方程的零狀態(tài)響應(yīng)。filter函數(shù)調(diào)用形式為 y=filter(b,a,x)其中 a=a0,a1,aN, b=bo,b1,bM分別表示差分方程系數(shù)。X表示輸入序列，y表示輸出序列。輸出序列的長(zhǎng)度和序列相同。例1-40 試用M=9點(diǎn)滑動(dòng)平均系統(tǒng) yk= 處理例1-39中的受噪聲干擾的正弦信號(hào)。 解： 由式（1-151）可知，M點(diǎn)滑動(dòng)平均系統(tǒng)可看成是N=0的差分方程。在調(diào)用filter函數(shù)時(shí)，調(diào)用參數(shù)a=1,b為有M個(gè)元素的向量，b中每個(gè)元素的值均為1/M。 %program 1_4 滑動(dòng)平均去噪 M=9; N=100;k=0:N

6、; s=10*sin(0.02*pi*k); x=s+randn(1,N+1); b=ones(M,1)/M; a=1; y=filter(b,a,x); plot(k,y,s,:); xlabel(k); ylabel（幅度） legend(yk,sk);程序的運(yùn)行結(jié)果如圖1-52所示。圖中的虛線表示未受噪聲干擾的原信號(hào)sk,實(shí)線為9點(diǎn)滑動(dòng)的響應(yīng)yk。比較圖1-50的信號(hào)可以看出，系統(tǒng)濾出了受干擾信號(hào)中的大部分噪聲，輸出信號(hào)相對(duì)輸入信號(hào)有4個(gè)樣本的延遲。1.9.4 利用MATLAB計(jì)算DTET當(dāng)序列的DTET可寫成的有理多項(xiàng)式時(shí)，可用MATLAB信號(hào)處理工具箱提供的freqz函數(shù)計(jì)算DTFT

7、的抽樣值。另外，可用MATLAB提供的abs、angle、real、imag等基本函數(shù)計(jì)算 DTFT的幅度、相位、實(shí)部、虛部。若X（）可表示為 X（）=則freqz的調(diào)用形式為 X=freqz(b,a,w) (1-153)式（1-153）中的b和 a分別是表示式（1-152）中分子多項(xiàng)式和分母多項(xiàng)式系數(shù)的向量，即 b=b0,b1,bM a=a0,a1,aNw為抽樣的頻率點(diǎn)，在以式（1-153）形式調(diào)用freqz函數(shù)時(shí)，向量w的長(zhǎng)度至少為2。返回值X就是DTFT在抽樣點(diǎn)w上的值。注意一般情況下，函數(shù)freqz的返回值X是復(fù)數(shù)。例1-41 已知離散系統(tǒng)的H（z）為 H（z）= 試畫出該系統(tǒng)的幅度響

8、應(yīng)。解： %program 1_5 計(jì)算系統(tǒng)幅度響應(yīng) b1=0.5009 -1.0019 0.5009; b2=0.5320 1.0640 0.5320; a1=1.0000 -0.8519 0.4167; a2=1.0000 0.8519 0.4167; b=conv(b1,b2);%計(jì)算H（z）的分子多項(xiàng)式 a=conv(a1,a2);%計(jì)算H（z）的分母多項(xiàng)式 w=linspace(0,pi,512); H=freqz(b,a,w); plot(w/pi,abs(H); ylabel(幅度); xlabel(Normalized frequency);系統(tǒng)幅度響應(yīng)如圖1-531.9.5

9、部分分式法的MATLAB實(shí)現(xiàn)MATLAB的信號(hào)處理工具箱提供了一個(gè)計(jì)算X（z）的部分分式展開的函數(shù)，它的調(diào)用形式如下 r,p,k=resifuez(b,a)其中調(diào)用參數(shù)b,a分別表示用z表示X（z）的分子和分母多項(xiàng)式。如果X（z）的部分分式展開為 X(z)=則residuez的返回參數(shù)r,p,k分別為 r=r r r r p=p1 p2 p3 p3 k=k1 k2同一極點(diǎn)p3在向量p中出現(xiàn)了兩次，表示p3是一個(gè)二階的重極點(diǎn)。Residuez也用于由r,p,k計(jì)算z表示的X（z）的分子和分母多項(xiàng)式，其調(diào)用形式為 b,a=residuez(r,p,k)例1-43 試用MATLAB計(jì)算X（z）=的部

10、分分式展開。 解：計(jì)算部分分式展開的MATLAB程序如下 %program 1_6 部分分式展開 b=1.5,0.98,-2.608,1.2,-0.144; a=1,-1.4,0.6,-0.072; r,p,k=residuez(b,a); disp(留書)；disp(r) disp(極點(diǎn))；disp(p) disp(常數(shù))；disp(k);程序運(yùn)行結(jié)果為留數(shù)0.7000+0.0000i 0.5000-0.0000i 0.3000極點(diǎn)0.6000+0.0000i 0.6000-0.0000i 0.2000常數(shù)0 2部分分式展開的結(jié)果為 X(z)=z反變換的結(jié)果為+0.5x(k+1)0.6)uk

11、+2k-11.9.6 用MATLAB計(jì)算系統(tǒng)的零極點(diǎn)MATLAB信號(hào)處理工具箱提供的tf2zp和zp2tf函數(shù)可進(jìn)行系統(tǒng)函數(shù)的不同表示形式的轉(zhuǎn)換。z正冪有理多項(xiàng)式表示的系統(tǒng)函數(shù)為 (1-154)用零點(diǎn)、極點(diǎn)和常數(shù)表示的一階因子形式的系統(tǒng)函數(shù)為 （1-155）MATLAB函數(shù)z,p,k=tf2zp(b,a)將有理多項(xiàng)式表示的系統(tǒng)函數(shù)轉(zhuǎn)換為一階因子形式的系統(tǒng)函數(shù)。b,a=zp2tf(z,p,k)將一階因子形式的系統(tǒng)函數(shù)轉(zhuǎn)換為有理多項(xiàng)式表示的系統(tǒng)函數(shù)。 例1-44試將下面的系統(tǒng)函數(shù)表示為一階因子形式。 （1-156） 解：用z正冪有理多項(xiàng)式表示的系統(tǒng)函數(shù)為 %program 1_7 確定一階因子形式

12、的系統(tǒng)函數(shù)b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;z,p,k=tf2zp(b,a);disp(零點(diǎn));disp(z);disp(極點(diǎn))；disp(p);disp(常數(shù))；disp(k);程序運(yùn)行結(jié)果為零點(diǎn) 0 0+0.2000i 0-0.2000i極點(diǎn) 0.8000 0.0000+0.4000i 0.0000-0.4000i常數(shù) 1系統(tǒng)函數(shù)的一階因子形式為 （1-157）為了在H（z）的表達(dá)式中不出現(xiàn)復(fù)數(shù)，也可將式（1-157）等價(jià)的寫成二階因子的形式 （1-158）當(dāng)H（z）是表示為式（1-154）的形式時(shí)，可用z,p,k=tf2zp(b,a)求出系統(tǒng)的零極點(diǎn)，從

13、而可根據(jù)極點(diǎn)的位置判斷系統(tǒng)的穩(wěn)定性，也可用函數(shù)zplane(b,a)畫出z平面的零極點(diǎn)分布來判斷系統(tǒng)的穩(wěn)定性。例-45已知一因果的系統(tǒng)的函數(shù)為試用函數(shù)zplane畫出系統(tǒng)的零極點(diǎn)分布，并判斷系統(tǒng)的穩(wěn)定性。解：%program 1_8 系統(tǒng)零極點(diǎn)分布 b=1 2 1; a=1 -0.5 -0.005 0.3; zplane(b,a);圖畫出了系統(tǒng)零極點(diǎn)分布。圖中符號(hào)表示零點(diǎn)。符號(hào)旁的數(shù)字表示零點(diǎn)的階數(shù)，符號(hào)x表示極點(diǎn)，圖中的虛線畫的是單位圓。由圖可知，該系統(tǒng)的極點(diǎn)全在單位圓內(nèi)，故系統(tǒng)是穩(wěn)定的。例1-46已知信號(hào)xk中含有角頻率分別為和的離散正弦信號(hào)。試設(shè)計(jì)一高通濾波器，濾出信號(hào)xk中低頻分量，保

14、留信號(hào)xk高頻分量。解：為了簡(jiǎn)單起見，假設(shè)數(shù)字濾波器是一個(gè)具有如下形式的長(zhǎng)度為的系統(tǒng)h0=h2=,h1=系統(tǒng)的查分方程+系統(tǒng)的頻率響應(yīng)為由上式可知系統(tǒng)的群延遲為為了濾除低頻分量，系統(tǒng)需滿足可用命令求解上面的的方程組得：下面的MATLAB程序?qū)崿F(xiàn)了濾波器的設(shè)計(jì)及信號(hào)的濾波。progam 1_9 簡(jiǎn)單數(shù)字濾波器的設(shè)計(jì)確定濾波器系數(shù)w1=0.015*pi;W2=0.4*pi;N=50;A=2*cos(w1) 1;2*cos(w1) 1;B=0;1x=A/B;a=1;b=x(1) x(2) x(1);%產(chǎn)生兩個(gè)余弦信號(hào)k=0:N;x1=cos(w1*k);x2=cos(w2*k);信號(hào)濾波y=filt

15、er(b,a,x1+x2);plot(k,y,r,k,x2,b-,x1+x2,k:);ylabel(幅度);xlabel(k)legend(yk,x2k,x1k+x2k);圖畫出了程序運(yùn)行結(jié)果。由圖可以看出在k=0，時(shí)，輸出響應(yīng)有波動(dòng)。這是因?yàn)橄到y(tǒng)響應(yīng)中存在瞬態(tài)響應(yīng)。當(dāng)時(shí)，系統(tǒng)進(jìn)入穩(wěn)態(tài)響應(yīng)。濾波后的信號(hào)相對(duì)與原信號(hào)有一個(gè)樣本的延遲，著是因?yàn)樵跁r(shí)，例1-47已知一段語音信號(hào)中混入了一頻率f=1200Hz正弦型干擾信號(hào)。語音信號(hào)的抽樣頻率f=22050Hz。試用二階帶阻濾波器濾除語音信號(hào)中的噪音信號(hào)。解：干擾信號(hào)的數(shù)字頻率為由式（1-110）可得取，由式（1-111）可得解上述方程得。當(dāng)?shù)闹禐?.

16、0619時(shí)，系統(tǒng)有一個(gè)極點(diǎn)在單位圓外，當(dāng)?shù)闹禐?.941時(shí)，系統(tǒng)的二個(gè)極點(diǎn)在單位圓內(nèi)。為了保證系統(tǒng)的穩(wěn)定，取0.9417。由和的值可得系統(tǒng)函數(shù)為%program 1_10 語音信號(hào)濾波Fn=1200;w=0.06;%讀入語音信號(hào)x,Fs=wavread(sample);%播放語音信號(hào)sound(x,Fs);pause;%設(shè)計(jì)濾波器w0=2*pi*Fn/Fs;beta=cos(w0);alpha=min(roots(1-2/cos(w) 1);a=1,-beta*(1+alpha),alpha;b=1 -2 *beta 1*(1+alpha)/2 ;%信號(hào)濾波y=filter(b,a,x);%播

17、放濾波后語音信號(hào)sound(y,Fs);圖畫出了濾波前后語音信號(hào)的頻譜。濾波前的語音信號(hào)的頻譜在z有一高峰，這是干擾所造成的。由濾波后語音信號(hào)的頻譜可看出，濾波器成功地制造了語音信號(hào)中的干擾。2.6 利用MATLAB實(shí)現(xiàn)信號(hào)DFT的計(jì)算在MATLAB信號(hào)處理工具箱中，函數(shù)dftmtx(N)可用來產(chǎn)生NN的DFT矩正D。NN的IDFT矩正D可用函數(shù)conj(dfmtx(N)/N來確定。此外，MATLAB提供了4個(gè)內(nèi)部函數(shù)用于計(jì)算DFT和IDFT，它們分別是： fft(x), fft(x,N), ifft(X), ifft(X,N)fft(x) 計(jì)算M點(diǎn)的DFT。M是序列x的長(zhǎng)度，即M=lengt

18、h(x)。fft(x,N) 計(jì)算N點(diǎn)的DFT。若M>N，則將原序列截短為N點(diǎn)序列，再計(jì)算其N點(diǎn)DFT；若M<N，則將原序列補(bǔ)零至N點(diǎn)，然后計(jì)算其N點(diǎn)DFT。ifft(X) 計(jì)算M點(diǎn)的IDFT。M是序列X的長(zhǎng)度。ifft(X,N) 計(jì)算N點(diǎn)IDFT。若M>N，則將原序列截短為N點(diǎn)序列，再計(jì)算其N點(diǎn)IDFT；若M<N，則將原序列補(bǔ)零至N點(diǎn)，然后計(jì)算其N點(diǎn)IDFT。為了提高fft與ifft的運(yùn)算效率，應(yīng)盡量使序列長(zhǎng)度M=2，或?qū)π蛄醒a(bǔ)零使N=2。實(shí)現(xiàn)例2-7的程序如下：%program 2_1 計(jì)算有限長(zhǎng)余弦信號(hào)頻譜N=30；%數(shù)據(jù)長(zhǎng)度L=512；%DFT的點(diǎn)數(shù)f1=100;

19、 f2=120; fs=600;T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;f=cos(2*pi*f1*t)+cos(2*pi*f2*t);F=fftshift(fft(f,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(F);ylabel(幅度譜)實(shí)現(xiàn)例2-8的程序?yàn)椋?%program 2_2 利用Hamming計(jì)算有限長(zhǎng)余弦信號(hào)頻譜 N=50；%數(shù)據(jù)長(zhǎng)度 L=512；%DFT的點(diǎn)數(shù) f1=100; f2=120; fs=600; %N=25；T=1/fs;ws=2*pi*fs;t=(0:N-1)*T;f=cos(2*pi*f1*t)

20、+cos(2*pi*f2*t); wh=(hamming(N); f=f.*wh; F=fftshift(fft(f,L);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);plot(w,abs(F);ylabel(幅度譜)例2-11 已知一長(zhǎng)度為16點(diǎn)的有限序列 試用MATLAB計(jì)算序列xk的16點(diǎn)和512點(diǎn)DFT。 解： %progam 2_3 Numerical Computation of DTFT Using FFT k=0:15; f=cos(2*pi*k*4./16); F_16=fft(f); F_512=fft(f,512); L=0:511; plot(L/512

21、,abs(F_512); hold on; plot(k/16,abs(F_16),0); set(gca,xtick,0,0.25,0.5,0.75,1); set(gca,ytick,0,2,4,6,8); grid on; xlabel(Normalized frequency); ylabel(Magnirude); hold off 從圖2-19可見，對(duì)長(zhǎng)度為16的序列xk進(jìn)行512點(diǎn)DFT，可以獲得頻譜函數(shù)X（）更多的細(xì)節(jié)，因?yàn)樾蛄蠳點(diǎn)的離散傅立葉變換就是序列X（）一個(gè)周期的N個(gè)等間隔抽樣，盡管從信息表達(dá)的角度，由序列xk16點(diǎn)DFT Xm足以恢復(fù)原序列xk。 利用MATLAB實(shí)現(xiàn)

22、由DFT計(jì)算線性卷積例2-12 試?yán)肕ATLAB實(shí)現(xiàn)由DFT計(jì)算有限序列線性卷積，并與線性卷積直接計(jì)算的結(jié)果進(jìn)行比較。 解：在MALAB中，序列線性卷積的直接結(jié)果是通過函數(shù)conv(x,h)來實(shí)現(xiàn)的。 %program 2_4 linear Convolution through DFT x=1 2 0 1; h=2 2 1 1; %Determine the length of the result of convolution L=length(x)+length(h)-1; %Compute the DFT by zero-padding XE=fft(x,L); HE=fft(h,L

23、); %Determine the IDFT of the product y1=ifft(XE.*HE); %plot the sequence generated by circular convolution %and the error from direct linear convolution k=0:L-1; subplot(2,1,1); setm(k,real(y1);axis(0 6 0 7); title(Result of Linear Convolution); xlabel(Time index k);ylabel(Amplitude); y2=conv(x,h);

24、 error=y1-y2;subplot(2,1,2); stem(k,abs(error); xlabel(Time index k);ylabel(Amplitude); title(Error Magnitude);由圖2-20（b）可見，由DFT計(jì)算的線性卷積的誤差非常小，這些誤差主要是計(jì)算機(jī)在計(jì)算DFT時(shí)，由計(jì)算機(jī)有限字長(zhǎng)而產(chǎn)生的舍入誤差與計(jì)算誤差。此外，MATLAB信號(hào)處理工具箱中提供了fftfilt函數(shù)，該函數(shù)是基于重疊相加法的原理，可以實(shí)現(xiàn)長(zhǎng)序列與短序列之間的線性卷積，其調(diào)用形式為 y=fftfilt(h,x,n)其中h:短序列；x：長(zhǎng)序列；n：DFT的點(diǎn)數(shù)。一般選擇n為2正整

25、次冪，以提高DFT運(yùn)算的效率。3.6 線性調(diào)頻z變換算法例3-3 已知序列xk=,(),試計(jì)算序列xk在單位圓上的CZT，其中。解：由于是計(jì)算序列xk在單位圓上的CZT，故A0=1，W0=1。計(jì)算序列czt的MATLAB程序如下：%program 3_1% Compute CZT of sequence xkN=13; %length of xkn=0:N-1;xn=(0.8).n;sita=pi/4; %phase of start pointfai=pi/6; %phase intervalA=exp(j*sita); %complex starting pointW=exp(-j*fai

26、); %complex ratio between pointsM=8; %length of czty=czt(xn,M,W,A); %call czt functionsubplot(2,1,1)stem(n,xn);xlabel(k);title(sequence xk);m=0”M-1;subplot(2,1,2);stem(m,abs(y);xlabel(m);title(CZT of xk);運(yùn)行結(jié)果如圖3-22所示。 4.5 用MATLAB實(shí)現(xiàn)濾波器設(shè)計(jì)MATLAB信號(hào)處理工具箱提供了常用的設(shè)計(jì)模擬低通濾波器的MATLAB函數(shù)。在實(shí)現(xiàn)應(yīng)用中，可方便地調(diào)用這些函數(shù)完成模擬濾波器的設(shè)

27、計(jì)。關(guān)于這些函數(shù)現(xiàn)分別介紹如下：Butterworth濾波器 N，wc=buttord(wp,ws,Ap,As,s) num,den=butter(N,wc,s)根據(jù)BW型濾波器的設(shè)計(jì)指標(biāo)，利用MATLAB函數(shù)buttord獲得BW型濾波器參數(shù)N和 wc。函數(shù)buttord的輸入?yún)?shù)wp和ws(rad/s)分別表示濾波器的通帶和阻帶截頻，Ap和As(dB)表示濾波器的通帶和阻帶衰減。s表示所設(shè)計(jì)的是模擬濾波器。函數(shù)buttord返回參數(shù)N為BW濾波器的階數(shù)，wc(rad/s)等于BW濾波器3 dB截頻。由于wc由阻帶方程確定，故由參數(shù)N、wc得出的濾波器在阻帶剛好滿足設(shè)計(jì)指標(biāo)，在通帶將超過設(shè)計(jì)

28、指標(biāo)。當(dāng)BW型濾波器的參數(shù)N、wc確定后，可用MATLAB函數(shù)butter獲得BW型濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式（num）和分母多項(xiàng)式（den）。Chebyshev I型濾波器 N，wc=cheblord(wp,ws,Ap,As,s) num,den=cheby1(N,Ap,wc,s)MATLAB函數(shù)cheblord返回參數(shù)N表示CB I型濾波器的階數(shù)，wc(rad/s)等于wp。參數(shù)N、wc的取值可使cheby1設(shè)計(jì)出的CB I型濾波器在通帶剛好滿足設(shè)計(jì)指標(biāo)。cheby1函數(shù)利用參數(shù)N、wc和Ap確定CB I型濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式（num）和分母多項(xiàng)式（den）。Chebyshev II

29、型濾波器 N，wc=cheb2ord(wp,ws,Ap,As,s) num,den=cheby2(N,As,wc,s)MATLAB函數(shù)cheb2ord返回參數(shù)N表示CBII型濾波器的階數(shù)，wc(rad/s)的取值可使cheby2設(shè)計(jì)出的濾波器在通帶剛好滿足設(shè)計(jì)指標(biāo)。cheby2函數(shù)利用參數(shù)N、wc和As確定CBII濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式（num）和分母多項(xiàng)式（den）。橢圓濾波器 N，wc=ellipord(wp,ws,Ap,As,s) num,den=ellip(N,Ap,wc,s)MATLAB函數(shù)ellipord返回參數(shù)N表示橢圓濾波器的階數(shù)，wc=wp。MATLAB函數(shù)ellip確定

30、階數(shù)為N，通帶衰減為Ap(dB),阻帶衰減為As(dB),通帶截頻為wc的橢圓濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式和分母多項(xiàng)式。在濾波器的階數(shù)確定后，采用了節(jié)中方案A重新計(jì)算濾波器的參數(shù)k1。故設(shè)計(jì)出的橢圓濾波器阻帶衰減超過設(shè)計(jì)指標(biāo)，而其他指標(biāo)剛好滿足設(shè)計(jì)要求。例 4-15 設(shè)計(jì)一個(gè)滿足下列指標(biāo)的模擬BW低頻濾波器。 解：%program 4_1 設(shè)計(jì)模擬Butterworth低頻濾波器%濾波器指標(biāo)wp=2*pi*1000;ws=2*pi*3000;Ap=1;As=50;%設(shè)計(jì)濾波器N，wc=buttord(wp,ws,Ap,As,s)num,den=butter(N,wc,s)fprintf(Orde

31、r of the fiter=%.Ofn,N)disp(Numerator polynomial);fprintf(%.4en,num);disp(Denominator polynomial);fprintf(%.4en,den);%計(jì)算Ap,As及增益響應(yīng)omega1=linspace(0,wp,500);omega2=linspace(wp,ws,200);omega3=linspace(ws,5*1000*pi*2,500);h1=20*log10(abs(freqs(num,den,omega1);h2=20*log10(abs(freqs(num,den,omega2);h3=20

32、*log10(abs(freqs(num,den,omega3);fprintf(Ap=%.4fn,max(-h1);fprintf(As=%.4fn,min(-h3);plot(omega1 omega2 omega3/(2*pi),h1 h2 h3);grid;xlabel(Frequency in Hz);ylabel(Gain in dB);程序運(yùn)行后的結(jié)果如下：order of the filter=6 Numerator polynomial0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+0000.0000e+000 0.0000e+000

33、1.4184e+023Denominator polynomial1.000e+000 2.7902e+004 3.8927e+008 3.4429e+0122.0301e+016 7.5889e+019 1.4184e+023Ap=0.7488 As=50.0000圖4-25畫出了該濾波器的增益響應(yīng)。例4-16 利用模擬橢圓低通濾波器重新設(shè)計(jì)例4-15。解：只需將例4-15程序的第6行和第7行改為N，wc=ellipord(wp,ws,Ap,As,s) num,den=ellip(N,Ap,wc,s)即可。程序的運(yùn)行結(jié)果如下： Order of the filter=4 Numerator

34、polynomial3.1615e-003 1.8190e-012 3.2366e+006 1.5259e-004 4.4736e+014Denominator polynomial1.0000e+000 5.9359e+003 5.8689e+007 1.9320e+011 5.0195e+014Ap=1.0000 As=51.3315圖4-26畫出了該濾波器的增益響應(yīng)。比較例4-15和例4-16可知，實(shí)現(xiàn)同樣設(shè)計(jì)指標(biāo)的濾波器，橢圓濾波器所需階數(shù)較低。 MATLAB信號(hào)處理工具箱提供了實(shí)現(xiàn)四種模擬域頻率變換的函數(shù)，它們分為：低通到低通的變換 numt,dent=1p21p(num,den,w

35、0)低通到高通的變換 numt,dent=1p2hp(num,den,w0)低通到帶通的變換 numt,dent=1p2bp(num,den,w0,B)低通到帶阻的變換 numt,dent=1p2bs(num,den,w0,B)其中num、den分別表示變換前模擬濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式和分母多項(xiàng)式。numt、dent分別表示變換后模擬濾波器系統(tǒng)函數(shù)的分子多項(xiàng)式和分母多項(xiàng)式。w0和B為變換中的參數(shù)。 例4-47 試設(shè)計(jì)一個(gè)滿足下列指標(biāo)的模擬BW型帶阻濾波器 解： （1）由式（4-77）確定低通到帶阻變換中的參數(shù) B=2， （2）由式（4-79） 由式（4-80）得原型低通濾波器的通帶截頻wp

36、_bar為 wp_bar=0.3714 （3）由 N，wc=buttord(wp_bar,1,Ap,As,s); num,den=butter(N,wc,s);得原型低通濾波器為 （4）由 numt,dent=1p2bp(num,den,w0,B)；得帶阻低通濾波器為 圖4-27畫出了該濾波器的增益響應(yīng)曲線。該濾波器的通帶和阻帶衰假分別為Ap1=0.05dB, Ap2=0.68dB, As1=As2=10dB。 MATLAB信號(hào)處理提供的impinvar(num,den,Fs)函數(shù)，可實(shí)現(xiàn)用脈沖響應(yīng)不變法將模擬濾波器轉(zhuǎn)化為數(shù)字濾波器，起調(diào)用形式為： numd,dend=impinvar(num

37、,den,Fs)式中num和den分別表示模擬濾波器系數(shù)函數(shù)H（s）的分子多項(xiàng)式和分母多項(xiàng)式，F(xiàn)s是脈沖響應(yīng)不變法中的抽樣頻率，單位是Hz。輸出變量numd和dend分別表示數(shù)字濾波器的系統(tǒng)函數(shù)H(z)分子多項(xiàng)式和分母多項(xiàng)式。若某個(gè)因果模擬濾波器的系統(tǒng)函數(shù)為 取T=1，則由 numd,dend=impinvar(2 14,1 2 5,1)得 numd=2.0000 2.3133 dend=1.0000 0.3062 0.1353所以由脈沖響應(yīng)不變法獲得的數(shù)字濾波器為 例4-18 利用Butterworth低通濾波器及脈沖響應(yīng)不變法設(shè)計(jì)滿足下列指標(biāo)的數(shù)字濾波器 解：設(shè)計(jì)滿足上述指標(biāo)數(shù)字濾波器的M

38、ATLAB程序如下： %program 4_2 脈沖響應(yīng)不變法設(shè)計(jì)BW型LP DF %DF BW LP 指標(biāo)Wp=0.1*pi;Ws=0.4*pi;Ap=1;As=25; Fs=1;%抽樣頻率（Hz）%確定模擬BW 指標(biāo)Wp=Wp*Fs;Ws=Ws*Fs;%確定AF階數(shù)N=buttord(Wp,Ws,Ap,As,s);%由通帶指標(biāo)確定3-dB截頻Wc=Wp/(10(0.1*Ap)-1)(1/2/N);%確定BW AFnuma,dena=butter(N,Wc,s);%確定 DFnumd,dend=impinvar(numa,dena,Fs);w=linspace(0,pi,512);h=fre

39、qz(numd,dend,w);%幅度歸一化DF的幅度響應(yīng)norm=max(abs(h);numd=numd/norm;plot(w/pi,20*log10(abs(h)/norm);grid;xlabel(Normalized frequency);ylabel(Gain,dB);disp(Numerator polynomial);fprintf(%.4en,numd);disp(Denominator polynomial);fprintf(%.4en,dend);%計(jì)算Ap和Asw=Wp Ws;h=freqz(numd,dend,w);fprintf(Ap=%.4fn,-20*log1

40、0(abs(h(1);fprintf(As=%.4fn,-20*log10(abs(h(2);運(yùn)行結(jié)果為Numerator polynomial-5.5515e-0172.3231e-0021.7880e-0020.0000e+000Denominator polynomial1.0000e+000-2.2230e+0001.7193e+000-4.5520e-001Ap=0.9985As=30.3240Ap=0.9985As=30.3240圖4-28畫出了該濾波器的增益響應(yīng) 雙線性變換法的MATLAB實(shí)現(xiàn) MATLAB信號(hào)處理工具箱提供的（num,den,Fs）函數(shù)可以用來實(shí)現(xiàn)雙線性變換，其

41、調(diào)用形式為： numd,dend=bibinear(num,den,Fs)式中num和den分別表示模擬濾波器系統(tǒng)函數(shù)H（s）的分子多項(xiàng)式和分母多項(xiàng)式，F(xiàn)s=1/T.輸出變量numd和dend分別數(shù)字濾波器系統(tǒng)函數(shù)H（s）的分子多項(xiàng)式和分母多項(xiàng)式。例4-19 三階歸一化Butterworth濾波器的系統(tǒng)函數(shù)為 解：取T=2，則由 den=conv(1 1,1 1 1); numd,dend=bilinear(1,den,0.5)得經(jīng)雙線性變換后的數(shù)字濾波器的分子多項(xiàng)式和分母多項(xiàng)式為 numd=0.1667 0.5000 0.5000 0.1667 dend=1.0000 0.0000 0.33

42、33 0.0000注意0.1667是1/6的近似值，0.3333是1/3的近似值，所以雙線性變換后的數(shù)字濾波器的系統(tǒng)函數(shù)H（z）為 例4-20 利用CB I型模擬濾波器和雙線性變換法，設(shè)計(jì)一滿足下列指標(biāo)的數(shù)字低通濾波器。 解：取T=2，由得模擬濾波器的頻率指標(biāo) 由 N,wc=cheb1ord(0.1854,0.7265,1,10,s) num,den=cheby1(N,1,wc,s)得模擬濾波器的分子多項(xiàng)式和分母多項(xiàng)式為 num= 0 0 0.0246 den=1.0000 0.1739 0.0277即 再由 numd,dend=bilinear(num,den,0.5)得雙線性變換后的數(shù)字濾

43、波器的分子多項(xiàng)式和分母多項(xiàng)式 numd=0.0205 0.0410 0.0205 dend=1.0000 -1.6185 0.7106即 上式表示的數(shù)字濾波器的增益響應(yīng)如圖4-29所示。該濾波器的通帶和阻帶衰減分別為。 用MATLAB實(shí)現(xiàn)數(shù)字濾波器設(shè)計(jì)MATLAB也提供了直接設(shè)計(jì)IIR數(shù)字濾波器的函數(shù)。設(shè)計(jì)的基本思想是用式（4-103）將數(shù)字濾波器的頻率指標(biāo)轉(zhuǎn)化為模擬濾波器的指標(biāo)，然后設(shè)計(jì)模擬濾波器，最后用雙線性變換把模擬濾波器轉(zhuǎn)換為數(shù)字濾波器。由于在設(shè)計(jì)中用的模擬濾波器的類型不同，所以給出了以下四組不同的函數(shù)BW型數(shù)字低通濾波器 N,Wc=buttord(Wp,Ws,Ap,As) num,d

44、en=butter(N,Wc)在函數(shù)buttord函數(shù)中，調(diào)用參數(shù)Wp、Ws是數(shù)字低通濾波器的歸一化的通帶和阻帶截頻。例如要求數(shù)字濾波器的，則Wp=0.1,Ws=0.4。Ap和As為濾波器的通帶和阻帶的衰減（dB）。返回的參數(shù)為和Wc，其中N為濾波器階數(shù)。函數(shù)butter中，調(diào)用參數(shù)N和Wc由函數(shù)buttord確定。返回參數(shù)num和den分別是數(shù)字濾波器的分子多項(xiàng)式和分母多項(xiàng)式的系數(shù)。CB I型數(shù)字低通濾波器 N,Wc=cheb1ord(Wp,Ws,Ap,As) num,den=cheby1(N,Ap,Wc)CB II 型數(shù)字低通濾波器 N,Wc=cheb1ord(Wp,Ws,Ap,As) n

45、um,den=cheby1(N,Ap,Wc)橢圓型數(shù)字低通濾波器 N,Wc=ellipord(Wp,Ws,Ap,As) num,den=ellip(N,Ap,As,Wc)利用上述四組MATLAB函數(shù)，也可設(shè)計(jì)數(shù)字高通、帶通和帶阻濾波器。調(diào)用方法可參見MATLAB手冊(cè)。 例4-21 設(shè)計(jì)滿足下列指標(biāo)的BW型和橢圓型數(shù)字低通濾波器。 解： %program 4_3 設(shè)計(jì)數(shù)字低通濾波器 %DF 指標(biāo) Wp=0.2;Ws=0.4; Ap=0.5;As=20; %設(shè)計(jì)BW型DF LP N,Wc=buttord(Wp,Ws,Ap,As); b1,a1=butter(N,Wc); %設(shè)計(jì)C型DF LP N,

46、Wc=ellipord(Wp,Ws,Ap,As); b2,a2=ellip(N,Ap,As,Wc); disp(BW型分子多項(xiàng)式)； fprintf(%.4en,b1); disp(BW型分母多項(xiàng)式)； fprintf(%.4en,a1); disp(C型分子多項(xiàng)式)； fprintf(%.4en,b2); disp(C型分母多項(xiàng)式)； fprintf(%.4en,a2); w=linspace(0,0.6*pi,200); h1=freqz(b1,a1,w); h2=freqz(b2,a2,w); plot(w/pi,20*log10(abs(h1),w/pi,20*log10(abs(h2

47、); legend(BW型，C型)； xlabel(Normalized frequency); ylabel(Gain,dB); axis(0 0.6 -50 1);程序運(yùn)行的結(jié)果為 BW型分子多項(xiàng)式 4.7792e-003 2.3896e-002 4.7792e-002 4.7792e-002 2.3896e-002 4.7792e-003 BW型分母多項(xiàng)式 1.0000e+000 -2.2360e+000 2.4061e+000 -1.3906e+000 4.2865e-001 -5.5151e-002 C型分子多項(xiàng)式 9.4379e-002 -1.5627e-002 -1.5627e-

48、002 9.4379e-002 C型分母多項(xiàng)式 1.0000e+000 -1.9853e+000 1.6032e+000 -4.6040e-001 圖4-30畫出了所設(shè)計(jì)的濾波器的增益響應(yīng)。由運(yùn)行結(jié)果知對(duì)同樣的設(shè)計(jì)指標(biāo)，用BW型模擬濾波器設(shè)計(jì)出的數(shù)字濾波器是五階的，而用C型模擬濾波器設(shè)計(jì)出的數(shù)字濾波器只有三階。 5.5 利用MATLAB實(shí)現(xiàn)FIR濾波器設(shè)計(jì)窗函數(shù)的計(jì)算 MATLAB提供了許多常用的窗函數(shù)，其中部分窗函數(shù)的調(diào)用形式為 w=hanning(N) w=hamming(N) w=blackman(N) w=Kaiser(N,beta)其中N是窗函數(shù)的長(zhǎng)度,beta是控制kaiser窗形狀的參數(shù)。返回的變量w是一個(gè)長(zhǎng)度為N的列向量，給出窗函數(shù)在N點(diǎn)的取值。 窗函數(shù)法設(shè)計(jì)FIR濾波器窗函數(shù)法設(shè)計(jì)FIR濾波器一般分為3個(gè)步驟：第1步估計(jì)FIR濾波器階數(shù)M（或長(zhǎng)度N）。如果用Kaiser窗時(shí)可用式（5-41）估計(jì)FIR濾波器階數(shù)；第2步確定所用的窗函數(shù)并計(jì)算出窗函數(shù)的值；第3步計(jì)算理想濾波器的