(實驗三)連續(xù)時間LTI系統(tǒng)的頻域分析

1、實驗三 連續(xù)時間LTI系統(tǒng)的頻域分析一、實驗?zāi)康?、掌握系統(tǒng)頻率響應(yīng)特性的概念及其物理意義；2、掌握系統(tǒng)頻率響應(yīng)特性的計算方法和特性曲線的繪制方法，理解具有不同頻率響應(yīng)特性的濾波器對信號的濾波作用；3、學(xué)習(xí)和掌握幅度特性、相位特性以及群延時的物理意義；4、掌握用MATLAB語言進(jìn)行系統(tǒng)頻響特性分析的方法。基本要求：掌握LTI連續(xù)和離散時間系統(tǒng)的頻域數(shù)學(xué)模型和頻域數(shù)學(xué)模型的MATLAB描述方法，深刻理解LTI系統(tǒng)的頻率響應(yīng)特性的物理意義，理解濾波和濾波器的概念，掌握利用MATLAB計算和繪制LTI系統(tǒng)頻率響應(yīng)特性曲線中的編程。二、實驗原理及方法1 連續(xù)時間LTI系統(tǒng)的頻率響應(yīng)所謂頻率特性，也稱為

2、頻率響應(yīng)特性，簡稱頻率響應(yīng)（Frequency response），是指系統(tǒng)在正弦信號激勵下的穩(wěn)態(tài)響應(yīng)隨頻率變化的情況，包括響應(yīng)的幅度隨頻率的變化情況和響應(yīng)的相位隨頻率的變化情況兩個方面。連續(xù)時間LTI系統(tǒng)的時域及頻域分析圖 上圖中x(t)、y(t)分別為系統(tǒng)的時域激勵信號和響應(yīng)信號，h(t)是系統(tǒng)的單位沖激響應(yīng)，它們?nèi)咧g的關(guān)系為：，由傅里葉變換的時域卷積定理可得到： 3.1或者： 3.2為系統(tǒng)的頻域數(shù)學(xué)模型，它實際上就是系統(tǒng)的單位沖激響應(yīng)h(t)的傅里葉變換。即 3.3由于H(jw)實際上是系統(tǒng)單位沖激響應(yīng)h(t)的傅里葉變換，如果h(t)是收斂的，或者說是絕對可積（Absolutly

3、integrabel）的話，那么H(jw)一定存在，而且H(jw)通常是復(fù)數(shù)，因此，也可以表示成復(fù)數(shù)的不同表達(dá)形式。在研究系統(tǒng)的頻率響應(yīng)時，更多的是把它表示成極坐標(biāo)形式： 3.4上式中，稱為幅度頻率相應(yīng)（Magnitude response），反映信號經(jīng)過系統(tǒng)之后，信號各頻率分量的幅度發(fā)生變化的情況，稱為相位特性（Phase response），反映信號經(jīng)過系統(tǒng)后，信號各頻率分量在相位上發(fā)生變換的情況。和都是頻率w的函數(shù)。對于一個系統(tǒng)，其頻率響應(yīng)為H(jw)，其幅度響應(yīng)和相位響應(yīng)分別為和，如果作用于系統(tǒng)的信號為，則其響應(yīng)信號為 3.5若輸入信號為正弦信號，即x(t) = sin(w0t)，則系

4、統(tǒng)響應(yīng)為 3.6可見，系統(tǒng)對某一頻率分量的影響表現(xiàn)為兩個方面，一是信號的幅度要被加權(quán)，二是信號的相位要被移相。由于和都是頻率w的函數(shù)，所以，系統(tǒng)對不同頻率的頻率分量造成的幅度和相位上的影響是不同的。2 LTI系統(tǒng)的群延時從信號頻譜的觀點看，信號是由無窮多個不同頻率的正弦信號的加權(quán)和（Weighted sum）所組成。正如剛才所述，信號經(jīng)過LTI系統(tǒng)傳輸與處理時，系統(tǒng)將會對信號中的所有頻率分量造成幅度和相位上的不同影響。從相位上來看，系統(tǒng)對各個頻率分量造成一定的相位移（Phase shifting），相位移實際上就是延時（Time delay）。群延時（Group delay）的概念能夠較好地反

5、映系統(tǒng)對不同頻率分量造成的延時。LTI系統(tǒng)的群延時定義為： 3.7群延時的物理意義：群延時描述的是信號中某一頻率分量經(jīng)過線性時不變系統(tǒng)傳輸處理后產(chǎn)生的響應(yīng)信號在時間上造成的延時的時間。如果系統(tǒng)的相位頻率響應(yīng)特性是線性的，則群延時為常數(shù)，也就是說，該系統(tǒng)對于所有的頻率分量造成的延時時間都是一樣的，因而，系統(tǒng)不會對信號產(chǎn)生相位失真（Phase distortion）。反之，若系統(tǒng)的相位頻率響應(yīng)特性不是線性的，則該系統(tǒng)對于不同頻率的頻率分量造成的延時時間是不同的，因此，當(dāng)信號經(jīng)過系統(tǒng)后，必將產(chǎn)生相位失真。3 用MATLAB計算系統(tǒng)頻率響應(yīng)在本實驗中，表示系統(tǒng)的方法仍然是用系統(tǒng)函數(shù)分子和分母多項式系數(shù)

6、行向量來表示。實驗中用到的MATLAB函數(shù)如下：H,w = freqs(b,a)：b,a分別為連續(xù)時間LTI系統(tǒng)的微分方程右邊的和左邊的系數(shù)向量（Coefficients vector），返回的頻率響應(yīng)在各頻率點的樣點值（復(fù)數(shù)）存放在H中，系統(tǒng)默認(rèn)的樣點數(shù)目為200點；Hm = abs(H)：求模數(shù)，即進(jìn)行運算，求得系統(tǒng)的幅度頻率響應(yīng)，返回值存于Hm之中。real(H)：求H的實部；imag(H)：求H的虛部；phi = atan(-imag(H)./(real(H)+eps)：求相位頻率相應(yīng)特性，atan()用來計算反正切值；或者phi = angle(H)：求相位頻率相應(yīng)特性；tao =

7、grpdelay(num,den,w)：計算系統(tǒng)的相位頻率響應(yīng)所對應(yīng)的群延時。計算頻率響應(yīng)的函數(shù)freqs()的另一種形式是：H = freqs(b,a,w)：在指定的頻率范圍內(nèi)計算系統(tǒng)的頻率響應(yīng)特性。在使用這種形式的freqs/freqz函數(shù)時，要在前面先指定頻率變量w的范圍。例如在語句H = freqs(b,a,w)之前加上語句：w = 0:2*pi/256:2*pi。 下面舉例說明如何利用上述函數(shù)計算并繪制系統(tǒng)頻率響應(yīng)特性曲線的編程方法。假設(shè)給定一個連續(xù)時間LTI系統(tǒng)，下面的微分方程描述其輸入輸出之間的關(guān)系 編寫的MATLAB范例程序，繪制系統(tǒng)的幅度響應(yīng)特性、相位響應(yīng)特性、頻率響應(yīng)的實部

8、和頻率響應(yīng)的虛部。程序如下：% Program3_1% This Program is used to compute and draw the plots of the frequency response % of a continuous-time systemb = 1; % The coefficient vector of the right side of the differential equationa = 1 3 2; % The coefficient vector of the left side of the differential equationH,w = f

9、reqs(b,a); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the real part of the frequency responseHi = imag(H); % Compute the imaginary part of the frequency responsesubplot(221)plot(w,Hm), gr

10、id on, title('Magnitude response'), xlabel('Frequency in rad/sec')subplot(223)plot(w,phai), grid on, title('Phase response'), xlabel('Frequency in rad/sec')subplot(222)plot(w,Hr), grid on, title('Real part of frequency response'), xlabel('Frequency in rad/

11、sec')subplot(224)plot(w,Hi), grid on, title('Imaginary part of frequency response'), xlabel('Frequency in rad/sec')三、實驗內(nèi)容及步驟實驗前，必須首先閱讀本實驗原理，了解所給的MATLAB相關(guān)函數(shù)，讀懂所給出的全部范例程序。實驗開始時，先在計算機上運行這些范例程序，觀察所得到的信號的波形圖。并結(jié)合范例程序所完成的工作，進(jìn)一步分析程序中各個語句的作用，從而真正理解這些程序。實驗前，一定要針對下面的實驗項目做好相應(yīng)的實驗準(zhǔn)備工作，包括事先編寫好相

12、應(yīng)的實驗程序等事項。給定三個連續(xù)時間LTI系統(tǒng)，它們的微分方程分別為系統(tǒng)1： Eq.3.1系統(tǒng)2： Eq.3.2系統(tǒng)3： Eq.3.3Q3-1 修改程序Program3_1，并以Q3_1存盤，使之能夠能夠接受鍵盤方式輸入的微分方程系數(shù)向量。并利用該程序計算并繪制由微分方程Eq.3.1、Eq.3.2和Eq.3.3描述的系統(tǒng)的幅度響應(yīng)特性、相位響應(yīng)特性、頻率響應(yīng)的實部和頻率響應(yīng)的虛部曲線圖。抄寫程序Q3_1如下：%Q3_1b = input('請輸入右邊向量系數(shù)'); % The coefficient vector of the right side of the differe

13、ntial equationa = input('請輸入左邊向量系數(shù)'); % The coefficient vector of the left side of the differential equationH,w = freqs(b,a); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the real p

14、art of the frequency responseHi = imag(H); % Compute the imaginary part of the frequency responsesubplot(221)plot(w,Hm), grid on, title('Magnitude response'), xlabel('Frequency in rad/sec')subplot(223)plot(w,phai), grid on, title('Phase response'), xlabel('Frequency in ra

15、d/sec')subplot(222)plot(w,Hr), grid on, title('Real part of frequency response'), xlabel('Frequency in rad/sec')subplot(224)plot(w,Hi), grid on, title('Imaginary part of frequency response'), xlabel('Frequency in rad/sec')執(zhí)行程序Q3_1，繪制的系統(tǒng)1的頻率響應(yīng)特性曲線如下：從系統(tǒng)1的幅度頻率響應(yīng)曲線看，

16、系統(tǒng)1是低通、高通、全通、帶通還是帶阻濾波器？答：執(zhí)行程序Q3_1，繪制的系統(tǒng)2的頻率響應(yīng)特性曲線如下：從系統(tǒng)2的幅度頻率響應(yīng)曲線看，系統(tǒng)2低通、高通、全通、帶通還是帶阻濾波器？答：執(zhí)行程序Q3_1，繪制的系統(tǒng)3的頻率響應(yīng)特性曲線如下：從系統(tǒng)3的幅度頻率響應(yīng)曲線看，系統(tǒng)3是低通、高通、全通、帶通還是帶阻濾波器？答：這三個系統(tǒng)的幅度頻率響應(yīng)、相位頻率相應(yīng)、頻率響應(yīng)的實部以及頻率響應(yīng)的虛部分別具有何種對稱關(guān)系？請根據(jù)傅里葉變換的性質(zhì)說明為什么會具有這些對稱關(guān)系？答：Q3-2 編寫程序Q3_2，使之能夠接受鍵盤方式輸入的輸入信號x(t)的數(shù)學(xué)表達(dá)式，系統(tǒng)微分方程的系數(shù)向量，計算輸入信號的幅度頻譜，系

17、統(tǒng)的幅度頻率響應(yīng)，系統(tǒng)輸出信號y(t)的幅度頻譜，系統(tǒng)的單位沖激響應(yīng)h(t)，并按照下面的圖Q3-2的布局，繪制出各個信號的時域和頻域圖形。圖Q3-2你編寫的程序Q3_2抄寫如下：% Q3_2b = input('Type in the right coefficient vector of differential equation：'); % The coefficient vector of the right side of the differential equationa = input('Type in the left coefficient vect

18、or of differential equation：'); % The coefficient vector of the left side of the differential equationw=-10:0.01:10;H= freqs(b,a,w); % Compute the frequency response HHm = abs(H); % Compute the magnitude response Hmphai = angle(H); % Compute the phase response phaiHr = real(H); % Compute the rea

19、l part of the frequency responseHi = imag(H); % Compute the imaginary part of the frequency responsesubplot(221);plot(w,Hm); grid on, title('Magnitude response'), xlabel('Frequency in rad/sec');subplot(223);plot(w,phai); grid on, title('Phase response'), xlabel('Frequency

20、 in rad/sec');subplot(222);plot(w,Hr); grid on, title('Real part of frequency response'); xlabel('Frequency in rad/sec');subplot(224);plot(w,Hi); grid on, title('Imaginary part of frequency response'); xlabel('Frequency in rad/sec');執(zhí)行程序Q3_2，輸入信號x(t) = sin(t) + si

21、n(8t)，輸入由Eq.3.3描述的系統(tǒng)。得到的圖形如下：此處粘帖執(zhí)行程序Q3_2所得到的圖形請手工繪制出信號x(t) = sin(t) + sin(8t) 的幅度頻譜圖如下：你手工繪制的信號x(t) = sin(t) + sin(8t) 的幅度頻譜圖與執(zhí)行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度頻譜圖是否相同？如不同，是何原因造成的？答：執(zhí)行程序Q3_2得到的x(t) = sin(t) + sin(8t) 的幅度頻譜圖實際上是另外一個信號x1(t)的幅度頻譜，這個信號的時域數(shù)學(xué)表達(dá)式為 x1(t) = 請利用傅里葉變換的相關(guān)性質(zhì)計算并繪制信號x1(t)的幅度頻

22、譜圖。 計算過程： 手工繪制的x1(t) 的幅度頻譜圖如下：結(jié)合所學(xué)的有關(guān)濾波的知識，根據(jù)上面所得到的信號的時域和頻域圖形，請從時域和頻域兩個方面解釋濾波的概念。答：Q3-3 編寫程序Q3_3，能夠接受從鍵盤輸入的系統(tǒng)微分方程系數(shù)向量，并分別繪制所給三個系統(tǒng)的群延時曲線圖。抄寫程序Q3_3如下：% Q3_3b1 = input('Type in the right coefficient vector of differential equation：'); % The coefficient vector of the right side of the differenti

23、al equationa1= input('Type in the left coefficient vector of differential equation：'); % The coefficient vector of the left side of the differential equationb2 = input('Type in the right coefficient vector of differential equation：'); % The coefficient vector of the right side of the

24、 differential equationa2= input('Type in the left coefficient vector of differential equation：'); % The coefficient vector of the left side of the differential equationb3 = input('Type in the right coefficient vector of differential equation：'); % The coefficient vector of the right

25、side of the differential equationa3= input('Type in the left coefficient vector of differential equation：'); % The coefficient vector of the left side of the differential equationw=-10:0.01:10;H1 = freqs(b1,a1,w); % Compute the frequency response Hphi1 = angle(H1);H2 = freqs(b2,a2,w); % Compute the frequency response Hphi2 = angle(H2);H3 = freqs(b3,a3,w); % Compute the frequency response Hphi3 = angle(H3);tao1= grpdelay(b1,a1,w);tao2= grpdelay(b2,a2,w);tao3= grpdelay(b3,a3,w);subplot(321);plot(w,phi1); grid on, title('Phase response of num1');subplot(323)