實(shí)驗(yàn)4Z變換和系統(tǒng)頻域特性的MATLAB實(shí)現(xiàn)實(shí)驗(yàn)?zāi)康膶W(xué)習(xí)通過Z變換來分析離散系統(tǒng)的頻率響應(yīng)，并用MATLAB實(shí)現(xiàn)。加深對(duì)系統(tǒng)的零、極點(diǎn)分布概念的理解。實(shí)例分析2.1通過Z變換分析求解系統(tǒng)的沖激響應(yīng)h(n)已知用線性常系數(shù)差分方程：y(n)=^Lay(n-k)+￡bx(n-r)表示的線性時(shí)k=1r=0不變系統(tǒng)其系統(tǒng)函數(shù)為:Ebz-rEbz-rY(z)rH(z)==r=0x(z)iE,1-"az-kkk=1上式為兩個(gè)關(guān)于z-i的多項(xiàng)式之比(3-1)即H(z)為有理分式。同時(shí)，式可以表示成部分分式的形式：H(z)=￡Rk+區(qū)Cz-kk=i1-Pk"k=0k則可以通過所熟悉的常見序列的Z變換形式求得H(z)的Z反變換，從而求得系統(tǒng)的沖激響應(yīng)函數(shù)h(n)。MATLAB提供了一個(gè)內(nèi)部函數(shù)residuez()，來計(jì)算有理多項(xiàng)式的留數(shù)和直接項(xiàng)，residuez()函數(shù)有幾種調(diào)用方式：[R,p,C]=residuez(b,a)，在已知以分子行向量b和分母行向量a下，得到列向量R含有留數(shù)，列向量p是極點(diǎn)位置，行向量C包含直接項(xiàng)；[b,a]=residuez(R,p,C)，將部分分式展開式轉(zhuǎn)換到分子行向量b和分母行向量a.MATLAB還提供了一個(gè)內(nèi)部函數(shù)impz(b,a,N)，在已知分子行向量b和分母行

向量a下，計(jì)算N點(diǎn)的單位沖激響應(yīng)h(n)。例3.1求系統(tǒng)：、0.001836+0.007344z-1+0.011016z-2+0.007374z-3+0.001836z-4H(z)=1-3.0544z-1+3.8291z-2-2.2925z-3+0.55075z-4的單位沖激響應(yīng)h(n)。解：由上述基本原理和MATLAB提供的函數(shù)，將下列指令編輯到“exe3impz.m〃文件中，可以得到所求系統(tǒng)的沖激響應(yīng)h(n)。%exe3impz.mh(n)求解b=[0.001836,0.007344,0.011016,0.007374,0.001836]；a=[1,-3.0544,3.8291,-2.2925,0.55075];[h,n]=impz(b,a,40);stem(n,h,'.');ylabel('h[n]');gridon_|-|1IIIIIII'■0510152025_|-|1IIIIIII'■0510152025303540圖4-1系統(tǒng)沖激響應(yīng)運(yùn)行“exe3impz.m〃文件將產(chǎn)生如圖3-1所示的序列。U.25IIIIIII例3.2求X(z)=(1V|z|V8)的Z反變換。解:X(z)=3z2-4z+1z-13-4z-1+z-2%exe3zinver.mz反變換求解

b=[0,1,0];a=[3,-4,1];[R,p,C]=residuez(b,a)[b,a]=residuez(R,p,C)運(yùn)行結(jié)果如下:R=0.5000-0.5000P=1.00000.3333因此得到因式分解后的X(因此得到因式分解后的X(z)=1/21/2所以Z反變—,3-4zt+z-21-zt1-1/3ztC0b-0.00000.33330a=1.0000-1.33330.3333換的結(jié)果為x(n)=1u(n)一1(；)nu(n)。JJJ2.2通過Z變換分析離散系統(tǒng)的頻率響應(yīng)和零、極點(diǎn)分布將式(3-1)表示的H(z)的分子、分母進(jìn)行因式分解，可采用根的形式表示多項(xiàng)式，即HnY(z)山(1-?-1)H(z)==X(z)「(1-dz-1)k=1其中，七為分子多項(xiàng)式的根，為系統(tǒng)函數(shù)的極點(diǎn)，A為比例稱為系統(tǒng)函數(shù)的零點(diǎn)，dk為分母多項(xiàng)式的根，稱。這樣系統(tǒng)函數(shù)就表示成了零極點(diǎn)的形式，零極點(diǎn)在2其中，七為分子多項(xiàng)式的根，為系統(tǒng)函數(shù)的極點(diǎn)，A為比例根據(jù)頻率響應(yīng)的定義，將z=ej^代入上式，系統(tǒng)頻率響應(yīng)為:H(eg)=H(乙)1.z=emMATLAB提供了一個(gè)內(nèi)部函數(shù)zplane(b,a),在已知系統(tǒng)函數(shù)的分子行向量b和分母行向量a下，畫出極點(diǎn)和零點(diǎn)。MATLAB還提供了一個(gè)內(nèi)部函數(shù)freqz()來求系統(tǒng)的頻率響應(yīng)，freqz()函數(shù)有幾種調(diào)用方式：[H,w]=freqz(b,a,N)，在已知分子行向量b和分母行向量a情況下得到N點(diǎn)的頻率向量w和N點(diǎn)的系統(tǒng)復(fù)頻率響應(yīng)向量H，頻率響應(yīng)在單位圓的上半圓的N個(gè)等分點(diǎn)上求值；[H,w]=freqz(b,a,N,’whole’)，在已知分子行向量b和分母行向量a情況下得到N點(diǎn)的頻率向量w和N點(diǎn)的系統(tǒng)復(fù)頻率響應(yīng)向量H，頻率響應(yīng)用環(huán)繞整個(gè)單位圓的N個(gè)等分點(diǎn)求值；H=freqz(b,a,w)，在已知分子行向量b和分母行向量a情況下得到在頻率向量w上的頻率響應(yīng)向量H。例3.3畫出因果系統(tǒng)j(n)=0.9j(n-1)+x(n)的零極點(diǎn)圖和頻率響應(yīng)曲線，并求出脈沖響應(yīng)h(n)。解：由上述的基本原理和MATLAB提供的函數(shù)，將下列指令編輯到“exe3freqz.m〃文件中，即可求得系統(tǒng)的零極點(diǎn)圖和頻率響應(yīng)。%exe3freqz.m系統(tǒng)求解b=[1,0];a=[1,-0.9];subplot(2,2,1);zplane(b,a);[H,w]=freqz(b,a,100);magH=abs(H);phaH=angle(H);subplot(2,2,3);plot(w/pi,magH);gridonxlabel('FrequencyinPiuints');ylabel('Magnitude');title('Magnituderesponse');subplot(2,2,4);plot(w/pi,phaH/pi);xlabel('FrequencyinPiuints');ylabel('Phase');

title('Phaseresponse');[h,n]=impz(b,a,40);subplot(2,2,2);stem(n,h,'.');ylabel('h[n]');gridon運(yùn)行“exe3freqz.m〃文件將產(chǎn)生如圖3-2所示的序列。-1-0.500.51RealPartMagnituderesponsePhaseresponse-0.2

Phaseresponse-□.4■00.5FrequercyinPiuints□51FrequencyinPiuints圖4-2系統(tǒng)零極點(diǎn)及頻率響應(yīng)例3.4定性畫出H（-□.4■00.5FrequercyinPiuints解：由上述基本原理和MATLAB提供的函數(shù)，將下列指令編輯到“exe3freqz2.m〃文件中，即可求得系統(tǒng)的零極點(diǎn)圖和頻率響應(yīng)。%exe3freqz2.m系統(tǒng)求解N=8;b=[1,0,0,0,0,0,0,0,-1];a=[1,0,0,0,0,0,0,0,0];subplot(2,1,1);zplane(b,a);title('零點(diǎn)-極點(diǎn)圖');[H,w]=freqz(b,a,100);magH=abs(H);phaH=angle(H);

subplot(2,2,3);plot(w/pi,magH);gridonxlabel('FrequencyinPiuints');ylabel('Magnltude');title('Magnituderesponse');subplot(2,2,4);plot(w/pi,phaH/pi);xlabel('FrequencyinPiuints');ylabel('Phase');title('Phaseresponse');運(yùn)行“exe3freqz2.m〃文件將產(chǎn)生如圖3-3所示的序列。極點(diǎn)-零點(diǎn)圖□00.51FrequencyinPiuints00.51FrequencyinPiuints圖4-3零極點(diǎn)及頻率響應(yīng)圖2.3由H(z)求解零極點(diǎn)及由零極點(diǎn)得到H(z)的表達(dá)式例3.5已知系統(tǒng)函數(shù)H(z)有5個(gè)零點(diǎn)，一個(gè)在z運(yùn)行“exe3freqz2.m〃文件將產(chǎn)生如圖3-3所示的序列。極點(diǎn)-零點(diǎn)圖□00.51FrequencyinPiuints00.51FrequencyinPiuints圖4-3零極點(diǎn)及頻率響應(yīng)圖解：MATLAB提供了一個(gè)內(nèi)部函數(shù)tf2zp()，來求H(z)=B(z)/A(z)的零極點(diǎn)和增益；提供了一個(gè)內(nèi)部函數(shù)zp2tf()，在零極點(diǎn)已知時(shí)求B(z)和A(z)的系數(shù)；提供了一個(gè)內(nèi)部函數(shù)roots(),來求一個(gè)多項(xiàng)式的根；提供了一個(gè)內(nèi)部函數(shù)poly(),由給定的根求出相應(yīng)多項(xiàng)式的系數(shù)。因此，將下列指令編輯到“exe2zphz.m〃文件中，可求解和驗(yàn)證。%exe3zphz.m零極點(diǎn)與系統(tǒng)函數(shù)的轉(zhuǎn)換B=poly([0?5*exp(j*pi/3),0?5*exp(-j*pi/3),2*exp(j*pi/3),2*exp(-j*pi/3),0?8]);L=length(B);A=zeros(1,L);A(1)=1;[Z,P,K]=tf2zp(B,A);sort(Z)[b,a]=zp2tf(Z,P,K)Z1=roots(B)運(yùn)行“exe3zphz.m〃文件結(jié)果如下，可自行比較驗(yàn)證。ans=0.2500-0.4330i0.2500+0.4330i0.80001.0000-1.7321i1.0000+1.7321ib=1.0000-3.30007.2500-6.70003.0000-0.8000a=100000Z1=1.0000+1.7321i1.0000-1.7321i0.80000.2500+0.4330i0.2500-0.4330i3.實(shí)驗(yàn)內(nèi)容3.1求系統(tǒng)〃/、0.0528+0.0797z-1+0.1295z-2+0.1295z-3+0.797z-4+0.0528z一5H(z)=1-1.8