信號處理綜合設(shè)計性實驗

1、課 程 設(shè) 計第1章 信號的時域分析1.1 連續(xù)信號的時域分析連續(xù)時間信號的時域分析就是將不同形式的信號波形用不同的時間函數(shù)來描述，例如：連續(xù)周期信號通過CTFS分解為不同諧波成分的三角函數(shù)或指數(shù)函數(shù)之和，連續(xù)非周期信號分解為頻率無限密集的虛指數(shù)函數(shù)的線性組合，即ICTFT。在時域還可以把實信號分解為奇信號與偶信號之和等等。其中常用的連續(xù)時間信號是信號時域分析的基礎(chǔ)，在此用Matlab提供的函數(shù)可以產(chǎn)生此類常用的連續(xù)時間信號。例1. 產(chǎn)生幅度為2，頻率為4Hz，初相為/6的正弦信號解：A=2;f0=4;phi=pi/6;w0=2*pi*f0;t=0:0.001:1; x=A*sin(w0*t+

2、phi);plot(t,x);ylabel('x(t)');xlabel('t');例2：產(chǎn)生一階躍信號x(t)=u(t).解：t=-1:0.01:5;x=(t>=0);plot(t,x);axis(-2,6,-0.1,1.1);例3：產(chǎn)生一語音信號解：t=0:0.01:1;plot(t,randn(1,length(t); 設(shè)計題目:用Matlab產(chǎn)生下列信號并與人工分析結(jié)果進行比較：(1) r (t)= t u (t) -1<t<10(2) x (t)=1+cos10 t -1<t<1(3) x (t)=(5e - t - 5e

3、 -3 t ) u (t) -1<t<5(4) x (t)=cos(2t) cos(20t) 0<t<5(5) x (t)=sin(t)/t -10<t<101.2 離散時間序列的時域分析及信號的運算常見基本離散序列有正弦序列、指數(shù)序列、矩形脈沖序列、單位脈沖序列與單位階躍序列等。這些序列可以由Matlab產(chǎn)生： 例1：產(chǎn)生一正弦序列 解：A=2;N=15;phi=pi/6;omega=pi/6;n=-10:10;y=A*sin(omega*n+phi);stem(n,y);ylabel('y(n)');xlabel('n')

4、;在信號的時域分析中，重要的一項內(nèi)容為信號的運算，包括信號的相加、相乘、乘方、卷積、相關(guān)等運算，以及求解信號的功率與能量。連續(xù)時間信號能量及功率的求解公式為：信號的能量：E= |x(t)|2dt=x(t) x(t)dt信號的功率：p=|x(t)|2dt=x(t)x(t)dt 因果離散序列能量及功率的求解公式為：序列的能量：E=|x(n)|2 =x(n)x(n)序列的功率：p=|x(n)|2 =x(n)x(n) 例2：求解x(t)=e -t u(t)-u(t-1)在0，1時間內(nèi)的能量。解：首先建立一內(nèi)部函數(shù)計算信號的瞬時能量：function f=powert(t)f=(abs(exp(-1.*

5、t).2; 計算信號在0，1時間內(nèi)的能量：power-t=quad(powert,0,1) power-t=04323例3：已知序列x(n)=0.8 n u(n),計算前十點的能量。解：N=10;n=0:N-1;x=(0.8).n;e=sum(abs(x).2)e =2.7458通過計算，序列前十點的能量占總能量的98.85%。離散序列的線性卷積和運算：例4：若x(n)=1,1,1,1,0,0, 計算離散序列的線性卷積和y(n)=x(n)*x(n)。解：x=1,1,1,1,0,0,;y=conv(x,x);subplot(2,1,1);stem(0:length(x)-1,x);ylabel(

6、'x(n)');xlabel('Time index n');subplot(2,1,2);stem(0:length(y)-1,y);ylabel('y(n)=x(n)*x(n)');xlabel('Time index n');常見的Matlab在信號處理應(yīng)用中產(chǎn)生信號的函數(shù)函 數(shù) 名功 能 sawtooth產(chǎn)生周期鋸齒波或三角波 square產(chǎn)生周期方波sinc產(chǎn)生sin(t)/t函數(shù)chirp產(chǎn)生掃頻信號pulstran產(chǎn)生脈沖串rectpuls產(chǎn)生矩形波tripuls產(chǎn)生三角形波diric產(chǎn)生dirichlet函數(shù)或周期

7、sinc函數(shù)函 數(shù) 名功 能 zeros(1,N)產(chǎn)生單位脈沖序列(n)ones(1,N)產(chǎn)生單位階躍序列u(n)rand(1,N)產(chǎn)生在0，1區(qū)間均勻分布的隨機信號randn(1,N)產(chǎn)生均值為0，方差為1的白噪聲信號運算的常用Matlab的實現(xiàn)運算名稱數(shù)學表達式MATLAB實現(xiàn)信號幅度變化y(n)=x(n)Y=alpha*x信號時移y(n)=x(n-k)Y=zeros(1,k),x信號翻轉(zhuǎn)y(n)=x(-n)Y=fliplr(x)信號累加y(n)=x(k)Y=cumsum(x)信號差分（或近似微分）y(n)=x(n+1)-x(n)Y=diff(x)信號求和y= x(n)Y=sum(x(n1

8、:n2)兩個信號相加y(n)= x1(n)+x2(n)Y=x1+x2兩個信號相乘y(n)= x1(n)x2(n)Y=x1.*x2兩個信號卷積y(n)=x1(k)x2(n-k)Y=conv(x1,x2)兩個信號相關(guān)Rxx(m)=x1(n)x2(n+m)Y=xcorr(x1,x2)設(shè)計題目：1.使用Matlab產(chǎn)生下列序列、作圖并與理論值進行比較：（1）x(n)=2(n+n0) (2) x(n)=(0.9) n sin(0.25n)+cos(0.25n) (3) 已知LTI離散系統(tǒng)，x(n)=1 1 1,h(n)=0 1 2 3,求y(n)(4) 已知x(t)=e 2 tu (t), y(t)=e

9、 - tu (t) , 求：x(t) * y(t)(5) 已知信號x(t)=(1+t/2)u(t+2)-u(t-2), 求x(t+2),x(t-2),x(-t),x(2t),-x(t)第2章 信號的頻域分析任一信號可以在時域?qū)ζ溥M行分析和描述，利用傅立葉變換理論也可以對其進行頻域分析，以便更好地對信號進行存儲、傳輸和處理，達到提取有用信號的目的。信號可分為四大類，與之對應(yīng)存在四種類型的傅立葉變換，成為信號頻譜分析的基礎(chǔ)。歸納如下表：信號類型變換名稱變換公式頻譜特性連續(xù)周期CTFS非周期離散連續(xù)非周期CTFT非周期連續(xù)離散周期DFS周期離散離散非周期DTFT周期連續(xù)（0=2）四種信號的變化規(guī)律為

10、：周期信號的頻譜是離散的、互為諧波關(guān)系的；非周期信號的頻譜是連續(xù)的；離散信號的頻譜是為周期的；連續(xù)信號的頻譜是非周期的。所謂信號的頻譜分析就是利用傅立葉變換的分析方法，找出與信號時域波形相對應(yīng)的頻譜函數(shù)的幅度、相位以及能量或功率的分布規(guī)律等，以便在頻域提取信號的特征。實際工程中，通過積分公式求取復雜信號的頻譜函數(shù)本身就比較困難，何況在許多情況下只是記錄了實際信號的一段波形或數(shù)據(jù)，而沒有對應(yīng)的解析表達式。若對這些信號進行頻譜分析，就必須利用離散傅里葉變換（DFT）。DFT表征一個在時域為N點有限長的序列x(n) 經(jīng)過傅里葉變換到頻域成為另一個N點有限長序列X (k)，即 ：=離散傅里葉反變換（I

11、DFT）定義為 可見，由于DFT變換對在時域、頻域都是離散的，可以通過計算機實現(xiàn)數(shù)值計算。而且DFT存在快速算法FFT，可以高速、高效地完成DFT運算。Matlab中提供了相應(yīng)函數(shù)以實現(xiàn)DFT變換對的計算，調(diào)用格式為： X=fft(x) 其按照基2時間抽取快速算法計算序列x（n）的傅里葉變換，當x (n) 的長度為2的整數(shù)次冪或者x(n)為實序列時，計算的時間會大大縮短。 X=fft(x,n) 其是補零或截短的n點傅里葉變換，當x(n)的長度小于n時，在x(n)的尾部補零使x(n)的長度達到n點；當x(n)的長度大于n時，將x(n)截短使x(n)的長度成n點；然后對補零或截短的數(shù)據(jù)進行快速傅里

12、葉變換。 x=ifft(X)和x=ifft(X，n)為相應(yīng)的反變換。 fftshift(x)將fft計算輸出的零頻移到輸出的中心。 DFT解決了用計算機對各類信號進行頻譜分析的問題，因此可以使用DFT對連續(xù)信號進行頻譜分析。但必須要對連續(xù)信號進行離散化，并且當信號長度為無限長時需要作截短處理。恰當?shù)卮_定取樣時間間隔T和相應(yīng)的時間長度L，是決定DFT結(jié)果是否符合實際的關(guān)鍵因素。如果不滿足取樣定理的約束條件，在時域欠取樣的情況下，會出現(xiàn)頻譜混疊而無法恢復原信號頻譜，因而不能從時域取樣點準確地重建原連續(xù)時間信號；如果截斷和選取的長度不合適，造成頻譜擴散使信號的能量和功率產(chǎn)生泄漏。利用DFT分析各類信

13、號的頻譜時，參數(shù)的選擇主要滿足時域抽樣定理、頻域采樣定理，依據(jù)信號時域、頻域能量分布情況，恰當?shù)剡x擇時域取樣間隔T和取樣長度L。當信號的頻帶寬度無限或很寬，則取頻域能量集中在95%-98%的頻帶寬度。即要在頻域計算信號的能量，選取fm （信號的最高截頻），即選擇時域取樣間隔T；在時域計算信號的能量，則可選取窗函數(shù)，對時間信號進行截短，即選擇時間長度L。2.1 利用DFT分析連續(xù)信號頻譜首先， 調(diào)用Matlab的fourier（x,t,w）函數(shù)可以在理論上求連續(xù)時間信號的頻譜。例1：求出下列信號的傅立葉幅度譜(1). 門函數(shù)脈沖信號x1(t)=(t+0.5)-(t-0.5)(2). 高斯信號x2

14、(t)=e解：(1)利用符號函數(shù)求門函數(shù)的頻譜syms t w %符號函數(shù)ut=sym('Heaviside(t+0.5)-Heaviside(t-0.5)');Fw=fourier(ut,t,w);FFw=maple('convert',Fw,'piecewise');FFP=abs(FFw);ezplot(FFP,-10*pi 10*pi)axis(-10*pi 10*pi 0 1) (2) 利用符號函數(shù)求高斯函數(shù)的頻譜syms w tft=exp(-t.2);fw=fourier(ft,t,w);ezplot(fw)2.1.1 連續(xù)周期信號

15、的頻譜分析應(yīng)用離散傅里葉變換DFT，分析連續(xù)周期信號xT(t)的頻譜。 連續(xù)周期信號在滿足一定條件下，可以通過傅立葉級數(shù)（CTFS）展開為一系列正弦信號的線性疊加，其頻譜函數(shù)X（k）是離散頻率的復函數(shù)，因而周期信號的頻譜結(jié)構(gòu)具有離散性和諧波性。對x(t)以T 為間隔進行取樣，長度為一個周期T0 ，dtT, , T0=NT,得到=連續(xù)周期信號的頻譜求解步驟： (1)根據(jù)取樣定理，確定時域取樣間隔T； (2)計算一個周期內(nèi)的取樣點數(shù)N。 (3)使用FFT命令作N點FFT計算，求得X (k)； (4)最后求得連續(xù)周期信號的頻譜為X （k）= X (k) 例2： 已知一個連續(xù)周期信號x (t) = 1

16、 cos (t) + 2sin (2t) + cos (3t), 用DFT計算其頻譜。解： 信號基頻 ，周期T02；最高次諧頻為33，所以N7。可取N8，16，32，64，進行比較。T02；N16；TT0N；周期為T0，F(xiàn)FT的點數(shù)為N，時域取樣間隔為T：X = 1 / N * fft (x , N); 用FFT計算其頻譜，其結(jié)果為：T0=2; N=16;T=T0/N;t=0:T:T0;x=1-cos(pi*t)+2*sin(2*pi*t)+cos(3*pi*t);X=1/N*fft(x,N);f=1/T/N*(-N/2:(N/2-1);%N,f=1/T/N*(-(N-1)/2:(N-1)/2

17、);subplot(2,1,1);stem(f,abs(fftshift(X);xlabel('Frequency(Hz)');ylabel('magnitude');subplot(2,1,2);stem(f,angle(fftshift(X);xlabel('Frequency(Hz)');ylabel('phase');理論分析如下：按x (t) = 1 cos (t) + 2sin (2t) + cos (3t)ej30t + e-j30t 2ej0t + e-j0t 2ej20t + e-j20t 2j= 1 - 2 故

18、得： 1, n = 0 -1/2, n = ±1 X ( k0) = -+j, n = ± 21/2, n = ± 30, 其他可見，MATLAB作出的幅頻特性與理論分析很吻合。例3：XT（t）為幅度A10，周期T01秒，脈沖寬度0.2秒的周期脈沖串。1） 計算基波和多少次諧波的功率和可以占信號總功率的90以上。2） 作N256點的fft，計算其頻譜。查看直流及各次諧波k = 0, 1,2的幅度。為便于觀察，在30Hz,30Hz內(nèi)繪出其頻譜。解： 1）周期信號是功率信號，信號總功率為：理論分析該信號的頻譜為： X(k0) = (A / T0)Sa(k0)/2 =

19、2Sa(0.2k) 下面計算前2k + 1項的總功率：k=input('k=');n=-k:k;x=2*sinc(0.2*n);p0=20;p=sum(x.2)/p0基波和前5次諧波的功率之和占信號90以上的能量。XT（t）的傅里葉級數(shù)：XT（t）2 + 1.87 ej0t + 1.51ej20t + 1.01ej30t + + 1.87 e - j0t + 1.51 e - j20t + 1.01e - j30t 2)計算N=256點的FFT, 基波和前127次諧波所包含的能量占總能量的99.6%。 N=256;T0=1; T=T0/N; t=0:T:T0; %t = -T0

20、 / 2: T: T0 / 2; x=5*square(2*pi*(t+0.1),20)+1; X=1/N*fft(x,N); f=1/T/N*(-30:30) ; %在 -30Hz, 30Hz內(nèi)繪出其頻譜。 subplot(2,1,1); plot(t,x); xlabel('Time (s)'); ylabel('Amplitude'); title('square wave and its spectrum'); y=abs(fftshift (X); subplot(2,1,2); stem (f,y(129-30:129+30); xl

21、abel('Frequency(Hz)'); ylabel('Magnitude');2.1.2 連續(xù)非周期信號的頻譜分析應(yīng)用離散傅里葉變換DFT，通過MATLAB實現(xiàn)對連續(xù)信號頻譜分析，以滿足理論分析和工程實際的需要。 連續(xù)非周期信號通過連續(xù)時間傅里葉反變換(ICTFT)，可以表達為頻率為無限密集的虛指數(shù)信號的線性組合，即其中，頻率函數(shù) X(j)可按下式求得：對x(t)以T 為間隔進行取樣，x(t)的時間長度為T 0 ，由dtT, , T 0= N T,= k F0得到= = 連續(xù)非周期信號頻譜的求解步驟： (1)根據(jù)取樣定理，確定時域取樣間隔T； (2)確定

22、時域截取的長度L(或窗函數(shù)的點數(shù)M) (3)確定頻域取樣點數(shù)N，要求NM (4)使用FFT命令做N點FFT計算X (k) 。 (5)X ()| TX(k)例4: 已知信號x (t) = e - t u (t) u (t 1), 利用FFT計算其幅度頻譜。解 : 1）此信號時域為無限長，因此要確定數(shù)值計算選取的長度T0。從能量的觀點來看，總能量E 2dt = e-2tdt=使用MATLAB計算0，T0時間內(nèi)信號的能量：首先做一個函數(shù)powert.m從時域計算信號的能量：function f = powert (t) f = (abs (exp (-1.*t) ) ). 2;計算0，1時間內(nèi)的能量

23、： >>power_t = quad (powert, 0, 1) power_t = 0.4323計算0，4時間內(nèi)的能量： >>power_t = quad (powert, 0, 4) power_t = 0.4998計算0，5時間內(nèi)的能量： >>power_t = quad (powert, 0, 5) power_t = 0.5000 >>power_t = quad (powert, 0, 14) power_t = 0.5000接著試下去可知：在0，T05內(nèi)，信號所包含的能量非常接近于原無限長信號的總能量，因此可以選取T0 5。2）計

24、算時域取樣間隔T此信號的頻譜為無限寬，在此選取占信號總能量95的頻帶寬度（fm, fm）, fm為其最高頻率。作一個函數(shù)powerw.m，從頻域計算信號的能量： function f = powerw (w) f = (abs (1./ (j.*w + 1) ) ) . 2;使用下面名為p.m的函數(shù)計算（m, m）內(nèi)信號包含的能量占總能量的比例： w_max = input (w_max=); power_t = quad ( powert, 0, 5) power_w = 1 / (2 * pi) * quad ( powerw , -1 * w_max, w_max) propotion

25、= power_w / power_ t信號頻譜X（j） 1（j + 1）信號總能量E =在0，m內(nèi)，信號所包含的能量原信號總能量的95： Em = (1/)tan-1()|0m = (1/2) ×0.95 tan-1 (m) = 0.95/2 所以m 12.706，T(m )12.706 0.2473。 T pi / 12.706 T = 0.24733)確定FFT點數(shù)N由于 N M T0T 5T N 5 T N 20.2222所以N應(yīng)大于20 ,取N=22。N=22; T=5/N t=0:T:5; x=exp(-1*t); X=T*fft(x,N); x0=1/T.*ifft(X

26、,N); subplot(2,1,1); plot(t,x,t(1:N),real(x0),'r-'); axis(0,5,0,1.2); xlabel( 'Time (s)' ); ylabel('x(t)'); w=2*pi/T/N.*(-N/2:N/2-1); y=1./(j*w+1); subplot(2,1,2); plot(w,abs(fftshift(X),'r-'); legend('|TXm| ','|X(w)|',0); xlabel('Frequency(rad/s)&

27、#39;); ylabel('|X(w)|and|TXm|');設(shè)計內(nèi)容：1用fourier函數(shù),理論上求下列連續(xù)時間信號的頻譜。(1). 三角脈沖信號x1(t)= (2). 單邊指數(shù)信號x2(t)=eu(t)2用DFT計算下列信號的頻譜：（1） （2） (3) 3 利用FFT計算信號x(t)=e- 2 tu(t)- u(t-2)的離散頻譜X(m)；利用FFT計算信號x(t)=e- 2 t u(t)的離散頻譜X(m)；要求： （1） 確定DFT計算的各參數(shù)（取樣間隔T，時域長度Tp，頻譜分辨率F等）； （2） 進行理論值與計算值比較，分析各信號頻譜分析的計算精度； 4產(chǎn)生一個淹

28、沒在噪聲中的信號x(t)，例如由50Hz和120Hz的正弦信號以及一個零均值的隨機噪聲疊加而成。確定分析長度和取樣速度，計算信號的頻譜；計算其功率譜密度并作圖，指出50Hz和120Hz的正弦成分以及噪聲；詳細列出檢測信號的步驟和原理。)(tx 5假設(shè)一實際測得的一段信號的長度為0.4秒，其表達式為 x(t)=cos(2f1t)+0.75 cos(2f2t) 其中f1=100Hz，f2=110Hz。為了用FFT近似計算信號的頻譜，必須對信號進行取樣。試確定一合適取樣頻率fs，用MATLAB的fft函數(shù)近似計算信號x(t)的頻譜。 y(t)=cos(2f1t)+0.2 cos(2f2t)，計算信號

29、y(t)的頻譜。 若使用Hamming窗時，由實驗確定計算出的頻譜能分辨最小譜峰間隔f和信號長度TL的關(guān)系。對不同參數(shù)的Kaiser窗，重新確定頻譜分辨最小譜峰間隔f和信號長度TL的關(guān)系。 思考題： (1)既然可以直接計算CTFT，為什么利用DFT分析連續(xù)信號譜？ (2)若信號持續(xù)時間無限，且無解析表達式，如何利用DFT分析其頻譜？ (3)在利用DFT分析連續(xù)信號頻譜時，會出現(xiàn)哪些誤差？如何克服或減弱？ (4)在利用DFT分析連續(xù)信號頻譜時，如何選擇窗函數(shù)？ (5)如何選擇取樣頻率？ (6)討論補零對計算結(jié)果的影響。 (7) 窗函數(shù)對頻譜分辨率有何影響嗎？怎樣提高頻譜分辨率？2.2 利用DFT

30、分析離散序列頻譜應(yīng)用傅里葉變換DFT，分析各種離散序列x(n)的頻譜。2.2.1 離散周期信號的分析離散周期信號可以展開成離散傅里葉級數(shù)，其中傅里葉系數(shù)如下所示：式中：N是序列的周期，n為時間離散變量， k為數(shù)字頻率離散變量，是k次諧波的數(shù)字頻率。所以離散周期信號的頻譜X (k)是一個以N為周期的周期性離散頻譜，各譜線之間的間隔為0=2/ N，而且存在著諧波的關(guān)系。 離散周期序列在時域與頻域都是離散的、以N為周期的序列。正確地在一個周期內(nèi)選取N，即可以準確地求取周期序列的頻譜。離散周期序列頻譜的求解步驟： （1）確定離散周期序列的基本周期N； （2）使用fft命令作N點FFT計算。頻率分辨率F

31、0=2/N（3）令X(k0 )=1/NX（k）。 例1：已知一個周期序列 x(n)=cos(n+)用FFT計算其頻譜。解 ：序列基頻0 216，周期N16，用程序求解頻譜如下：N=16; n=0:N-1;x=cos(pi/8*n+pi/3);X=1/N*fft(x,N);omega=2*pi/N*(n-N/2);subplot(2,1,1); stem(omega,abs(fftshift(X); axis(-pi,pi,0,1);ylabel('Magnitude'); xlabel('Frequency (rad)' );subplot(2,1,2); st

32、em(omega,angle(fftshift(X); axis(-pi,pi,-4,4);ylabel('Phase'); xlabel('Frequency (rad) ' );2.2.2 離散非周期序列的頻譜分析通過離散時間傅里葉變換(DTFT)可求得非周期序列的頻譜密度函數(shù)，是數(shù)字頻率的周期連續(xù)函數(shù)。離散非周期信號的頻譜結(jié)構(gòu)是連續(xù)的且具有以2為周期的周期性。類似于對連續(xù)信號的譜分析，可以使用MATLAB提供的fft函數(shù)計算離散非周期信號的頻譜。對于離散非周期信號，當序列長度有限時，可以求得準確的頻譜樣值。若序列很長或無限長，則由于截短必然產(chǎn)生泄漏誤差以及

33、混疊誤差，使計算的結(jié)果只能是頻譜樣值的近似值。 離散非周期信號頻譜的求解步驟： （1）確定序列的長度L。根據(jù)能量分布，當序列為無限長時需要進行截短。 （2）確定作FFT的點數(shù)N；根據(jù)頻域取樣定理，為使時域波形不產(chǎn)生混疊必須?。?NL 。（3）使用fft命令作N點FFT計算X（k）。 1 -4n4,例2：已知一個有限長脈沖序列x（n）= 0 其它用FFT計算其頻譜。解： 該序列的持續(xù)時間長度為L2M19，所以要求FFT的點數(shù)N 9，取N51。 M=4; N=51; n=-(N-1)/2:(N-1)/2; %若N為偶，n = -N / 2: (N / 2 - 1); x=ones(1,M+1),z

34、eros(1,N-2*M-1),ones(1,M); X=fft(x,N); subplot(2,1,1); stem(n,fftshift(x); ylabel('x n'); xlabel('Time n'); omega = 2 * pi / N * n;由于X（）為實數(shù)，可直接利用下列程序畫出圖形： subplot (2, 1, 2 ); stem (omega, real (fftshift(X); ylabel('X k'); xlabel('Frequency (rad)'); hold on;利用下列程序在同一張圖

35、上作出理論分析得到的頻譜， N1=51; for n=-(N1-1)/2:(N1-1)/2 omega=2*pi/N1*n; if omega=0 X0(n+(N1-1)/2+1)=2*M+1; else X0(n+(N1-1)/2+1)=sin(2*M+1)/2*omega)./sin(1/2*omega); end end n=-(N1-1)/2:(N1-1)/2; omega=2*pi/N1*n; c=plot(omega,real(X0),'r-'); legend(c,'X (omega)'); hold off ; 設(shè)計內(nèi)容： 1 利用DFT計算序列

36、的頻譜；2 利用DFT計算序列的頻譜；3 有限長序列， 0n31，分別用N=32，N=60，N=120點DFT計算其頻譜。要求：(1)確定DFT計算的各參數(shù)； (2)進行理論值與計算值比較，分析各信號頻譜分析的計算精度； (3)詳細列出利用DFT分析離散信號頻譜的步驟； (4)寫出實驗原理。 4.離散序列，0n63，利用FFT分析其頻譜。 (1) 對x(n)做64點FFT，繪出信號的譜N；能分辨出兩個頻率嗎？ (2) 對x(n)補零到256點計算FFT，能分辨出兩個頻率嗎？(3) 改用非矩形窗計算FFT，能分辨出兩個頻率嗎？ 思考題： (1)既然可以直接計算DTFT，為什么利用DFT分析離散信

37、號譜？ (2)若信號持續(xù)時間無限，且無解析表達式，如何利用DFT分析其頻譜？ (3)在利用DFT分析離散信號頻譜時，會出現(xiàn)哪些誤差？如何克服或減弱？ (4)在利用DFT分析離散信號頻譜時，如何選擇窗函數(shù)？ (5)補零和增加信號長度都可以提高頻譜分辨率，兩者有何本質(zhì)區(qū)別？ 第3章 系統(tǒng)分析與設(shè)計信息系統(tǒng)就是為了有效的傳輸和交換信息而對載荷信息的信號進行加工、處理和變換的設(shè)備的總稱。系統(tǒng)根據(jù)輸入、輸出信號的不同分為連續(xù)系統(tǒng)和離散系統(tǒng)。其分析方法有卷積法，零輸入、零狀態(tài)法，變換域法等。3.1 連續(xù)系統(tǒng)分析3.1.1 卷積積分法求解系統(tǒng)的響應(yīng) 對LTI連續(xù)系統(tǒng)，輸入信號為x(t)，輸出信號為y(t)，

38、系統(tǒng)的沖激響應(yīng)為h(t)，則輸入、輸出之間的關(guān)系為：時域：y(t)=x(t)*h(t)頻域：Y(j)=X(j) H(j)其中線性卷積的計算引入了DFT分析連續(xù)信號的頻譜。例1：圖（A）和（B）分別表示一個LTI連續(xù)系統(tǒng)的輸入x(t)和單位沖激響應(yīng)h(t)，試利用DFT的快速算法計算該系統(tǒng)的零狀態(tài)響應(yīng)。 x(t) h(t)（A） （B） 1 1 3 4 解：T=1/8; Tx=3 ; Mx=Tx/T;Th=4; Mh =Th/T;L=Mh+Mx-1;N=64;x=ones(1,Mx),zeros(1,N -Mx);h=ones(1,Mh),zeros(1,N -Mh);X=T.*fft(x,N)

39、;H=T.*fft(h,N);Y=X.*H;y=1/T.*ifft(Y,N);figure(1);t=0:N-1*T;subplot(3,1,1); plot(t,x); ylabel('x(t)'); title('Time domain');subplot(3,1,2); plot(t,h); ylabel('h(t)');subplot(3,1,3); plot(t(2:N), real(y(1:N-1);ylabel('y(t)');xlabel('Time (s)');figure(2);Fs= 1/T;

40、w=0:N-1*Fs/N;subplot(3,1,1); plot(w-Fs/2,fftshift(abs(x); ylabel('|X(w)|');title('Frequency domain');subplot(3,1,2); plot(w-Fs/2,fftshift(abs(H); ylabel('|H (w)|');subplot(3,1,3); plot(w-Fs/2,fftshift(abs(Y);ylabel('|Y(w)|');xlabel ('Frequency (Hz)');3.1.2 零輸入

41、、零狀態(tài)法求解系統(tǒng)的響應(yīng)MATLAB的控制工具箱（control toolbox）里包含了許多可用于分析線性時不變（LTI）系統(tǒng)的函數(shù)，使用命令help control可以查看控制工具箱里的這些函數(shù)。 對于線性時不變的連續(xù)系統(tǒng)，在時域中其數(shù)學模型用常系數(shù)線性微分方程來描述。如下：在拉普拉斯變換中，可用系統(tǒng)函數(shù)H(s)來描述，且根據(jù)其零，極點的分布情況，可以決定系統(tǒng)的結(jié)構(gòu)及系統(tǒng)的穩(wěn)定性。 另外還可以用狀態(tài)變量分析法進行系統(tǒng)描述。此設(shè)計完成用atlab實現(xiàn)幾種分析方法的相互轉(zhuǎn)換。在Matlab中，描述系統(tǒng)的傳遞函數(shù)（系統(tǒng)函數(shù)）型tf(transfer function),零極點型zp(zero p

42、ole)以及狀態(tài)變量型ss(state space)三種方式可以方便的轉(zhuǎn)換。Matlab相應(yīng)的語句為：tf2zp系統(tǒng)函數(shù)型轉(zhuǎn)換到零極點型tf2ss系統(tǒng)函數(shù)型轉(zhuǎn)換到狀態(tài)變量型zp2tf零極點型轉(zhuǎn)換到系統(tǒng)函數(shù)型zp2ss零極點型轉(zhuǎn)換到狀態(tài)變量型ss2tf狀態(tài)變量型轉(zhuǎn)換到系統(tǒng)函數(shù)型ss2zp狀態(tài)變量型轉(zhuǎn)換到零極點型例2：已知描述LTI系統(tǒng)的數(shù)學模型為常系數(shù)線性微分方程：+8+19+12y(t)=2+10 f ( t )將其轉(zhuǎn)換為零,極點型。狀態(tài)變量型。解：1 .num=2,10;den=1 8 19 12;z,p,k=tf2zp(num,den) 2num=2,10；den=1 8 19 12a b

43、 c d=tf2ss(num,den)MATLAB提供了用于分析LTI系統(tǒng)的函數(shù)，可對微分方程描述的系統(tǒng)的響應(yīng)進行仿真。使用MATLAB命令可以計算系統(tǒng)的各種時間響應(yīng)：如單位沖激響應(yīng)、單位階躍響應(yīng)、零輸入響應(yīng)、零狀態(tài)響應(yīng)和完全響應(yīng)。 在調(diào)用這些函數(shù)時，需要用系數(shù)向量表示系統(tǒng)。在MATLAB中可使用向量b=bm,bm-1,b1,b0和向量a=an,an-1,a1,a0，分別為系統(tǒng)函數(shù)分子多項式和分母多項式的系數(shù)，注意均按s的降冪排列其系數(shù)直至s0。 （1） 用命令impulse可以求解系統(tǒng)的單位沖激響應(yīng)： impulse(sys)：計算并畫出系統(tǒng)的沖激響應(yīng)，sys可以是用命令tf、zpk或ss建

44、立的系統(tǒng)函數(shù)。 impulse(sys, t)：計算并畫出系統(tǒng)在向量t定義的時間內(nèi)的沖激響應(yīng)。 Y=impulse(sys, t)：向量Y保存對應(yīng)時間的系統(tǒng)輸出值。例3： 已知系統(tǒng)的微分方程為,計算該系統(tǒng)的單位沖激響應(yīng)。解：b=2,8;a=1,5,6; t=0:0.1:10; y=impulse(b,a,t); （2） 用命令step求解系統(tǒng)的單位階躍響應(yīng) step (sys)：計算并畫出系統(tǒng)的沖激響應(yīng)，sys可以是用命令tf、zpk或ss建立的系統(tǒng)函數(shù)。 step (sys, t)：計算并畫出系統(tǒng)在向量t定義的時間內(nèi)的沖激響應(yīng)。 Y=step (sys, t)：向量Y保存對應(yīng)時間的系統(tǒng)輸出值

45、。 （3） 使用命令lsim求解系統(tǒng)在任意輸入下的響應(yīng) lsim(sys, x, t)計算并畫出任意輸入下系統(tǒng)的零狀態(tài)響應(yīng)，sys可以是用命令tf、zpk或ss建立的系統(tǒng)函數(shù)，t為時間向量定義時間范圍，x為系統(tǒng)的輸入。 lsim(sys, x, t, zi) 計算并畫出系統(tǒng)的完全響應(yīng)，但sys必須是狀態(tài)空間形式的系統(tǒng)函數(shù)，zi為系統(tǒng)的初始狀態(tài)。 Y=lsim(sys, x, t, zi)：向量Y保存對應(yīng)時間的系統(tǒng)輸出值。 b=2,8;a=1,5,6; sys=tf(b,a); t=0:10/300:10; x=exp(-t); y=lsim(sys,x,t); plot(t,y); 當系統(tǒng)具有

46、初始狀態(tài)時，使用lsim(sys, x, t, zi)命令格式，其中sys必須以狀態(tài)方程形式求解，而且系統(tǒng)的初始條件必須轉(zhuǎn)化為狀態(tài)向量在0時刻的初始值。 例如求上述系統(tǒng)在初始條件為y(0-)=-3,y(0-)=0，輸入x(t)=e t u(t)時的完全響應(yīng)。使用下面命令可得到狀態(tài)方程形式的系統(tǒng)函數(shù)。 b=2,8;a=1,5,6; A,B,C,D=tf2ss(b,a); sys=ss(A,B,C,D) 輸出為各矩陣的值，依次為：ABCD 計算零輸入響應(yīng)的程序如下： b=2,8;a=1,5,6; A,B,C,D=tf2ss(b,a); sys=ss(A,B,C,D); t=0:10/300:10;

47、 x=zeros(1,length(t); zi=-9/2,3/4; y=lsim(sys,x,t,zi); 3.1.3 系統(tǒng)函數(shù)的零、極點對系統(tǒng)的頻響特性H(j)的影響：連續(xù)因果LTI連續(xù)系統(tǒng)的系統(tǒng)函數(shù)H(s)的極點全部位于左半平面時，系統(tǒng)的頻響特性H (j)可由H(s)求出。使用MATLAB命令freqs(b,a)可計算系統(tǒng)的頻率特性H (j)。 例4：已知某系統(tǒng)的系統(tǒng)函數(shù)為：H(S) = ，求系統(tǒng)的幅頻率曲線|H(j) |。 解： b=1;a=conv(1,1,1,1,1); H,w=freqs(b,a); plot(w,abs(H); xlabel(Frequency(rad/s);y

48、label(Amplitude); title(Magnitude reponse);由圖可見，該系統(tǒng)具有低通濾波特性。 設(shè)計內(nèi)容：1已知某LIT系統(tǒng)的微分方程為 ，用Matlab求該系統(tǒng)的單位沖激響應(yīng)和階躍響應(yīng),并與理論值進行比較。2研究具有以下零極點連續(xù)系統(tǒng)： (1)1個極點s=-0.1，增益 k =1；(2)1個極點s =0，增益k =1； (3)2個共軛極點s=±j5，增益 k =1 ；(4)2個共軛極點s=-0 .1± j5，增益k =1； (5)零點在s =0. 5，極點在s=-0 .1±j5，增益k =1； (6)零點在s =0. 5，極點在s=0

49、.1±j5，增益k =1； 試完成下列任務(wù)： （1） zpk和tf命令建立系統(tǒng)的系統(tǒng)函數(shù)；作出系統(tǒng)的零極點圖； （2） 分析系統(tǒng)是否穩(wěn)定？若穩(wěn)定，作出系統(tǒng)的幅頻特性曲線； （3） 作出系統(tǒng)的沖激響應(yīng)波形； （4） 詳細列出根據(jù)零極點分析系統(tǒng)特性的過程。 3已知系統(tǒng)函數(shù)為 ,y(0-)=1，輸入x (t) = u(t)，計算該系統(tǒng)的響應(yīng)，并與理論分析結(jié)果做比較；詳細列出系統(tǒng)響應(yīng)分析的步驟。 1. 已知一個連續(xù)系統(tǒng)的s域模型為：H(S)=,求x(t)輸入分別為u(t),sintu(t),e-tu(t)時，求解系統(tǒng)的輸出，并與理論分析結(jié)果做比較. 2. .尋找一個實際系統(tǒng)（如馬達、受控過程

50、等）的數(shù)學模型，利用MATLAB求解系統(tǒng)的響應(yīng)； 思考題: (1) 系統(tǒng)函數(shù)零極點對系統(tǒng)頻響的影響？ (2) 設(shè)計因果穩(wěn)定、實系數(shù)的低通、高通、帶通、帶阻濾波器，如何布置零極點？ (3)零極點對沖激響應(yīng)波形有何影響？ (4)因果穩(wěn)定系統(tǒng)對系統(tǒng)函數(shù)零極點的要求是什么？ (5)零極點知識在實際技術(shù)中的應(yīng)用？(6) 連續(xù)系統(tǒng)響應(yīng)的計算機求解可以分為哪些方法？各是什么原理？3.2 離散系統(tǒng)分析3.2.1 卷積和法求解系統(tǒng)的響應(yīng)對LTI離散系統(tǒng)，若輸入序列為x(n)，輸出序列為 y(n)，系統(tǒng)的單位響應(yīng)為h(n)，則輸入、輸出之間的關(guān)系為：時域：y(n)=x(n)h(n)頻域：Y(k) =X(k) H(

51、k)利用循環(huán)卷積與線性卷積的關(guān)系，采用DFT計算系統(tǒng)的響應(yīng)。其中線性卷積的計算引入了DFT分析離散序列的頻譜。 3.2.2 零輸入、零狀態(tài)法求解系統(tǒng)的響應(yīng)LTI離散系統(tǒng)可用系統(tǒng)函數(shù)H（z）的分子分母多項式形式、零極點-增益形式或狀態(tài)空間形式來描述。 （1） 系統(tǒng)函數(shù)的分子分母多項式形式H (z) ：離散系統(tǒng)的差分方程為： 對上式兩邊同取z變換，可得：在MATLAB中可使用向量b=b0,b1,bM和向量a=1,a1,a2,aN分別表示分子多項式和分母多項式的系數(shù)，注意從z按z的降冪排列其系數(shù)。 離散系統(tǒng)響應(yīng)的求解 除可以使用MATLAB命令lsim求解外，還可以使用命令filter來求解系統(tǒng)響應(yīng)

52、。 例1：已知系統(tǒng)函數(shù)為H(z)=, 求 :(1)系統(tǒng)的脈沖響應(yīng)h(n)； (2)輸入x(n)=u(n)，求系統(tǒng)的零狀態(tài)響應(yīng)y(n): (3) 輸入x(n)=u(n)，初始條件y(-1)=1,y(-2)=2,求系統(tǒng)的完全響應(yīng)y(n)解: (1)計算前11個時刻的h(n);N=11;b=1,2;a=1,0.4,-0.12; x=1,zeros(1,N-1); y=filter(b,a,x) Columns 1 through 5 1.0000 1.6000 -0.5200 0.4000 -0.2224 Columns 6 through 10 0.1370 -0.0815 0.0490 -0.0294 0.0176 Columns 11 through 15 -0.0106 0.0063 -0.0038 0.0023 -0.0014 Column 16 0.0008(2)計算前11個時刻的零狀態(tài)響應(yīng) N=11;b=1,2;a=1,0.4,-0.12; x=ones(1,N); y=filter(b,a,x) Columns 1 through 5 1.0000 2.6000 2.0800 2.4800 2.2576 Columns 6 through 10 2.3946 2.3131 2.3621 2.