連續(xù)周期信號時域采樣頻譜分析程序設(shè)計

連續(xù)周期信號時域采樣頻譜分析程序設(shè)計[摘要]信號的頻譜分析是信號與系統(tǒng)分析的基礎(chǔ)，本文分析了用數(shù)值計算以及使用MATLAB軟件與mathematica軟件的方法實現(xiàn)連續(xù)周期信號的時域采樣頻譜分析。數(shù)值計算即使用傅里葉級數(shù)展開方法進行連續(xù)周期信號的頻譜分析，使用軟件計算即設(shè)計實用的程序進行周期信號的時域采樣頻譜分析。[關(guān)鍵詞]頻譜分析；傅里葉級數(shù)；MATLAB軟件；mathematica軟件ContinuousperiodicsignalspectrumanalysisintimedomainsamplingprogramdesignAbstract:RepairingSpectralanalysisofthesignalisthebasisofsignalsandsystemsanalysis,thispaperanalyzesthespectralanalysisoftime-domainsamplestoachievecontinuousperiodicsignalswithnumericalcalculationsandtheuseofMATLABsoftwareandmathematicasoftwaremethod.NumericalcalculationsthattheuseofFourierseriesexpansionmethodforspectralanalysisofcontinuousperiodicsignal,calculatedusingthesoftwareprogramthatisdesignedandpracticalperiodicsignalinthetimedomainsamplingspectrumanalysis.Keywords:SpectrumAnalysis;Fourierseries;MATLABsoftware;mathematicasoftware目錄1引言 12周期信號頻譜分析的傳里葉級數(shù)展開方法 12.1傅里葉級數(shù) 12.1.1周期信號的分解 22.1.2指數(shù)形式傅里葉級數(shù) 32.2連續(xù)時間傅里葉級數(shù)的基本性質(zhì) 52.2.1線性 52.2.2卷積特性 52.2.3微分特性 52.3技巧性方法 52.3.1頻移方法 52.3.2求導(dǎo)方法 73周期信號的頻譜分析方法 73.1已知變化規(guī)律的周期連續(xù)信號頻譜識別方法 73.2未知變化規(guī)律的周期連續(xù)信號頻譜識別方法 83.3未知周期信號的頻譜的理論方法和計算技巧 84Mathematica軟件計算方法 105結(jié)束語 146致謝 15參考文獻 151引言隨著計算機的發(fā)展，數(shù)字信號處理的理論與技術(shù)得到飛速發(fā)展，20世紀(jì)60年代以來，我國形成了一系列的數(shù)字信號處理的理論與算法，比如，數(shù)字濾波器，快速傅立葉變換（FFT），這些都是數(shù)字信號處理的技術(shù)基礎(chǔ)，隨著信息科技的飛速發(fā)展，信號處理取得了重大的飛躍。信號的頻譜分析是數(shù)字信號處理中的一個很重要的研究課題，對信號進行頻譜分析，是對其進行傅里葉變換，得到其振幅譜與相位譜。對于信號來說，分模擬信號與數(shù)字信號。進行頻譜分析時，對于模擬信號來說，首先對其進行抽樣，使其離散化，然后利用離散傅里葉變換（DFT）或者快速傅里葉變換（FFT），然后對其幅度（ABS）和相位（ANGLE）的圖像進行分析，而對于數(shù)字信號來說，則可直接進行離散傅里葉變換或快速傅里葉變換。在現(xiàn)實生活中對于信號進行頻譜分析具有重要的意義。通過對信號頻譜的分析，可以得到信號的頻率結(jié)構(gòu)，了解信號的頻率成分或系統(tǒng)的特征。在此基礎(chǔ)之上，可實現(xiàn)對信號的跟蹤控制，從而實現(xiàn)對系統(tǒng)狀態(tài)的早期預(yù)測，發(fā)現(xiàn)潛在的危險并診斷可能發(fā)生故障的原因，對系統(tǒng)參數(shù)進行識別及校正。因此，頻譜分析是揭示信號特征的重要方法，也是處理信號的重要手段。頻譜分析工作在許多科學(xué)技術(shù)領(lǐng)域都會遇到。首先，在無線電技術(shù)的許多方面，例如通訊、地震測量、電信、導(dǎo)航、雷達、電子對抗、空間技術(shù)等；此外，由于光象波、機械震動、沖擊、響聲等各種非電量都可以通過所謂換能器轉(zhuǎn)換成電流或電壓的變化來方便地分析，所以，頻譜分析在各種振動、噪聲、電聲、發(fā)動機、建筑、生物、醫(yī)學(xué)等領(lǐng)域也起重要作用。許多工程設(shè)計都要對振動頻率進行認(rèn)真分析。如果外界強迫振動的頻率與某物體的固有頻率相近會發(fā)生共振。如電動機轉(zhuǎn)動時會引起臺座的受迫振動，作用力的頻率如果接近機器某部分的固有頻率，就可能因共振而引起損壞；輪船和碼頭會經(jīng)常收到海浪的周期性沖擊，沖擊頻率如果與船體或碼頭結(jié)構(gòu)的固有頻率相近，也會引起船頭或碼頭的損壞。這就需要研究清楚各海域海浪的頻率范圍，避免共振的發(fā)生；在地震過程中，建筑物的固有頻率如果與地面振動的固有頻率相近，及時地震烈度較低，也可能激發(fā)建筑物產(chǎn)生劇烈震動而導(dǎo)致?lián)p壞。所以，建筑工程師對于一棟建筑物的動力學(xué)特征，如固有頻率，振型，阻尼等是需要認(rèn)真研究清楚的，特別要避開經(jīng)常出現(xiàn)起主導(dǎo)作用的地震頻率；研制發(fā)動機的科學(xué)家們?yōu)榱伺宄l(fā)動機活塞在不同運行速度下的性能，他們在活塞邊緣事先嵌入一個微型無線電發(fā)射機，這個發(fā)射機可以在實際運行條件下把活塞所收到的各種應(yīng)力和振動的詳細信息通過無線電信號發(fā)射過來，科學(xué)家對這些信號進行頻譜分析和其他分析就可以對活塞性能進行正確直觀的評價；研制飛機的科學(xué)家們需經(jīng)常做一中所謂疲勞實驗。因為在飛機飛行過程中，他的機翼一方面收到發(fā)動機的振動，一方面受到不均勻氣流的不斷沖擊，引起機翼金屬應(yīng)力不斷變化。雖然這成千上萬次的小振動幅度和能量均很小，但是累計下來也是很大的能量，為了避免發(fā)生意外，科學(xué)家們也要對這些振動進行頻譜分析，以保證飛機的正常運行。在工程領(lǐng)域中，MATLAB軟件和mathematica軟件是一種倍受程序開發(fā)人員青睞的語言，對于一些需要做大量數(shù)據(jù)運算處理的復(fù)雜應(yīng)用以及某些復(fù)雜的頻譜分析算法MATLAB顯得游刃有余NOTEREF_Ref388997070\f\h[1]。2周期信號頻譜分析的傳里葉級數(shù)展開方法2.1傅里葉級數(shù)周期信號是定義在（）區(qū)間，每隔一定時間T，按相同規(guī)律重復(fù)變化的信號，如圖2.1所示，它可以表示為圖2.1周期信號（2.1）式中m為任意整數(shù)。時間T被稱為該信號的重復(fù)周期，簡稱周期。周期的倒數(shù)稱為該信號的頻率NOTEREF_Ref388997070\f\h[3]。由式：（2.2）可知，周期信號f(t)在區(qū)間可以展開成在完備正交信號空間中的無窮級數(shù)。如果完備的正交函數(shù)集是三角函數(shù)集或指數(shù)函數(shù)集，那么，周期信號所展開的無窮級數(shù)就分別稱為“三角型傅里葉級數(shù)”或“指數(shù)型傅里葉技術(shù)”，統(tǒng)稱為傅里葉級數(shù)。需要指出，只有當(dāng)周期信號滿足狄里赫利條件時，才能展開成傅里葉級數(shù)。通常遇到的周期信號都滿足該條件。狄里赫利條件是：（1）函數(shù)在任意有限區(qū)間內(nèi)連續(xù)，或只有有限個第一類間斷點（當(dāng)t從左或右趨于這個間斷點時，函數(shù)有有限的左極限或右極限）；（2）在一周期內(nèi)，函數(shù)有有限個極大值或極小值。2.1.1周期信號的分解設(shè)有周期信號，他的周期是T，角頻率,它可分解為(2.3)式(2.3)中的系數(shù)a，b稱為傅里葉系數(shù)。它可由式(2.4)求得。為簡便，式(2.2)的積分區(qū)間取為或。由式(2.3)可得傅里葉系數(shù)n=0，1，2，…(2.5)n=0，1，2，…(2.6)式中T為函數(shù)的周期，為角頻率。由式(2.5)和式(2.6)可見，傅里葉系數(shù)和都是n的函數(shù)，其中是n的偶函數(shù)，即；而是n的奇函數(shù)，即有.將式（2.3）中同頻率項合并，可寫成如下形式(2.7)式中,n=1，2，…(2.8)如將式（2.8）的形式轉(zhuǎn)化為式（2.1-1）的形式，他們系數(shù)之間的關(guān)系為n=1，2，…(2.9)由式（2.8）可見，是n的偶函數(shù)，即有；而是n的奇函數(shù)，即有。傅里葉系數(shù)的這些重要性質(zhì)是很有用的。式（2.9）表明，任何滿足狄里赫利條件的周期函數(shù)可分解為直流和許多余弦或正弦分量。其中第一項是常數(shù)項，他是周期信號中所包含的直流分量；式中第二項稱為基波或一次諧波，它的角頻率與元周期信號相同，是基波振幅，是基波初相角；式中第三項稱為二次諧波，他的頻率是基波頻率的二倍，是二次諧波的振幅，是其初相角。以此類推，還有三次，四次、…諧波。一般而言，稱為n次諧波，是n次諧波的振幅，是其初相角。式（2.7）表明，周期信號可以分解為各次諧波分量。2.1.2指數(shù)形式傅里葉級數(shù)三角函數(shù)形式的傅里葉級數(shù)含義比較明確，但運算常感不便，因而經(jīng)常采用指數(shù)形式的傅里葉級數(shù)。由于（2.10）所以式（2.9）可以寫為（2.11）將上式第三項中的n用-n代換，并考慮到到n是偶函數(shù)，即；是n的奇函數(shù)，即，則上式可寫為如將上式中的寫成（其中=0），則上式可以寫為（2.12）令復(fù)數(shù)量，稱其為復(fù)傅里葉系數(shù)，簡稱傅里葉系數(shù)，其模為，相角為，則傅里葉級數(shù)的指數(shù)形式為（2.13）根據(jù)式（2.9），傅里葉系數(shù)（2.14）將式（2.5）和式（2.6）代入上式，得n=o,1，2，…（2.15）這就是求指數(shù)形式傅里葉級數(shù)的復(fù)系數(shù)的公式。式（2.13）表明，任意周期信號可分解為許多不同頻率的虛指數(shù)信號（）之和，其各分量的復(fù)數(shù)幅度（或相量）為。2.2連續(xù)時間傅里葉級數(shù)的基本性質(zhì)2.2.1線性若，，則對于任意常數(shù)，，有（2.16）2.2.2時移特性若，則（2.17）式中，為常數(shù)。2.2.2卷積特性若，，則（2.18）圖2.2時域卷積運算圖2.3頻域相乘運算2.2.3微分特性設(shè)是以為周期的周期信號，其傅里葉系數(shù)為（2.19）則導(dǎo)數(shù)的傅里葉系數(shù)是（2.20）若的傅里葉系數(shù)是，則的傅里葉系數(shù)是(2.21)2.3技巧性方法2.3.1頻移方法有些周期信號的周期規(guī)律不明顯或周期變化計算過于復(fù)雜，則可通過頻移和增加（減小）直流增量來改變信號位置，由此可使信號由非奇非偶信號變成奇信號或者偶信號，然后進行傅里葉級數(shù)展開時就會比較快NOTEREF_Ref388997070\f\h[5]。某信號在時域中乘以，相當(dāng)于頻域中頻譜向右移；某信號在時域中乘以，相當(dāng)于頻域中頻譜向左移；也就是：時域頻域表2.1頻移特性前面兩個結(jié)論可以這么直觀的理解：信號由無數(shù)不同旋轉(zhuǎn)角速度（）的旋轉(zhuǎn)向量疊加而成。圖2.4信號頻移當(dāng)與相乘之后，原來的所有旋轉(zhuǎn)向量的旋轉(zhuǎn)角速度都將增加，達到，對應(yīng)的頻譜也將由原來的位置向右搬移。圖2.5信號向右頻移當(dāng)與相乘之后，原來的所有旋轉(zhuǎn)向量的旋轉(zhuǎn)角速度都將減少，達到，對應(yīng)的頻譜也將由原來的位置向左搬移。圖2.6信號向右頻移2.3.2求導(dǎo)方法有些周期信號的周期規(guī)律不明顯或周期變化計算過于復(fù)雜，則可對其進行求導(dǎo)（求導(dǎo)只是對其一個周期求導(dǎo)，因為信號的所有信息在一個周期均有完全的反應(yīng)），這樣就是將復(fù)雜的信號分析改變成一個周期內(nèi)的函數(shù)進行分析，可以減少計算過程，若經(jīng)過一次求導(dǎo)以后函數(shù)依舊復(fù)雜，可再次對其進行求導(dǎo)。給出一個周期的表達式，，周期信號表達式為，對其求一次導(dǎo)則有轉(zhuǎn)換可得：。對其求二次導(dǎo)則有轉(zhuǎn)換可得：。3周期信號的頻譜分析方法3.1已知變化規(guī)律的周期連續(xù)信號頻譜識別方法通過上面分析我們已經(jīng)知道，對于這類信號要獲得其頻譜，我們只需將此周期連續(xù)信號做傅里葉級數(shù)展開，即Ｆ(s)．設(shè)f(t)為一個周期連續(xù)信號，T為其周期，f1=1/T是信號的基波頻率，為信號的基波角頻率，則f(t)可以表示為，（n=0,1,2………）。對f(t)進行傅里葉級數(shù)展開．(式2.22)(式2.23)(式2.24)(式2.25)(式2.26)通過上面推導(dǎo)即可計算已知變化規(guī)律的周期連續(xù)信號的頻譜和振，從而畫出相應(yīng)信號f(t)的幅頻特性桿狀圖。對于周期連續(xù)信號f(t)，若要計算其連續(xù)頻譜，可對f(t)取傅里葉變換FT[f(t)]，其定義為：(式2.27)其中密度函數(shù)F(jw)是一個復(fù)函數(shù)，它可寫為(式2.28)式中和分別是頻譜函數(shù)的模和相位。發(fā)現(xiàn)是由頻率位置各不相同的一系列狄拉克(-k1)帶權(quán)疊加，對于這樣的頻譜函數(shù)是不能直接畫幅頻特性圖的，可在各頻率坐標(biāo)k1處以狄拉克函數(shù)前面的權(quán)系數(shù)為高度畫幅頻桿狀圖來圖示周期連續(xù)信號所含頻率成份。3.2未知變化規(guī)律的周期連續(xù)信號頻譜識別方法未知變化規(guī)律的周期連續(xù)信號通常是通過對信號采樣，然后分析采樣信號的頻譜來獲得連續(xù)信號的頻譜，這種做法首先要滿足采樣定理。為了使實際信號在采樣后能夠不失真的再現(xiàn)，采樣頻率必須大于信號最高頻率的兩倍。其次，為獲得連續(xù)信號的頻譜，要對采樣數(shù)據(jù)進行離散傅里葉變換（DFT）來獲得離散信號的頻譜。我們假設(shè)周期連續(xù)信號采樣得到的序列為x(n)，設(shè)x(n)為有限長序列，長度為N，但它取自周期信號，我們可以認(rèn)為在采樣參數(shù)選擇恰當(dāng)時，對它進行周期延拓所得周期序列的最小周期與被采樣信號的周期是一致的。通過采樣定理的學(xué)習(xí)我們得知采樣信號的頻譜是被采樣信號頻譜的周期延拓，所以我們只需要對一個采樣周期進行研究便可得到信號的全部信息。現(xiàn)在我們?nèi)《ú蓸有蛄衳(n)，長度為N（稱為一個主值區(qū)間）對它進行頻譜分析即可。對x(n)取DFT，計算過程如下：具體計算步驟如下：(式2.29)(式2.30)式中是待采樣得周期連續(xù)信號，T是采樣周期，是采樣后的有限長離散序列，與是等同的。由上式可知有限長離散信號頻譜的采樣是由其連續(xù)頻譜經(jīng)過采樣得到的。由此得DTFT是計算離散序列的頻譜函數(shù)的，而DFT則是計算離散序列頻譜的采樣的。3.3未知周期信號的頻譜的理論方法和計算技巧例試畫出的振幅譜和相位譜解：的周期信號，題中所給的表達式可視為的傅里葉級數(shù)展開式。據(jù)可知，其基波頻率=（rad/s），其基本周期T=2s，=2，3，6分別為二、三、六次諧波頻率。且有：其余=0圖3.1信號的振幅譜圖3.2信號的相位譜圖3.3信號的雙邊頻譜振幅譜圖3.4信號的雙邊頻譜相位譜4Mathematica軟件計算方法例1利用mathematica軟件快速對周期信號進行頻譜選擇NOTEREF_Ref388997070\f\h[2]x[t_]=5Cos[100t+/4]+2Cos[300t-/8]+0.5Cos[500t];Ts=0.001;fs=1/Ts;n=1000;ts=Range[0,n-1]*Ts;xn=x[ts];ListPlot[Transpose[{ts,xn}][[20;;42]],Joined->True]n=20;T=n*Ts;yn=xn[[1;;20]];Y=2Fourier[yn,FourierParameters->{-1,-1}];Y[[1]]/=2;AY=Abs[Y[[1;;n/2]]];Am=Max[AY];th=0.08;F=Range[0,Floor[n/2-1]]*fs/n;data1=Transpose[{F,AY}];data2=Transpose[{AY,F,Arg[Y[[1;;n/2]]]}];ListPlot[data1,Filling->Axis]sp=Select[data2,#[[1]]>th*Am&]//Chop運行程序可得到一下結(jié)論：周期信號波形圖：圖4.1周期信號波形圖圖4.2信號幅頻桿狀圖頻點、振幅、相位分別為：例1利用mathematica軟件對圖4.2周期信號進行頻譜分析。圖4.3信號圖形Clear["f*","t*",u,s,T,Ts,L];u[t_]=UnitStep[t];T=0.001;q=5;Ts=0.00001;L=Round[T/Ts];s[t_]=(4q/T)((t(u[t]-u[t-T/8])+(t/2+T/16)(u[t-T/8]-u[t-3T/8])+(T/4)(u[t-3T/8]-u[t-5T/8])+(9T/16-t/2)(u[t-5T/8]-u[t-7T/8])+(T-t)(u[t-7T/8]-u[t-T])));ts=Range[0,L-1]*Ts;sn=s[ts];data=Transpose[{ts,sn}];p1=ListPlot[data,AspectRatio1/3,AxesLabel{"t","ladder(t)"}]S=2*Fourier[sn,FourierParameters{-1,-1}];S[[1]]/=2;A=Abs[S[[1;;(L/2)]]]//Chop;Print["{f0,f1,f2,f3}=",{0,1,2,3}/T];Print["{A0,A1,A2,A3}=",A[[1;;4]]];Print["{0,1,2,3}=",[[1;;4]]];=Arg[S[[1;;(L/2)]]];ListPlot[Table[{k/T,A[[k+1]]},{k,0,L/2-1}],AspectRatio1/3,AxesLabel{"f/Hz","A(f)"},FillingAxis,PlotRangeAll]ListPlot[Table[{k/T,[[k+1]]},{k,0,L/2-1}],AspectRatio1/3,AxesLabel{"f/Hz","(f)"},FillingAxis]sc[t_]=Sum[A[[k+1]]Cos[2*(k/T)*t+[[k+1]]],{k,0,40}];p2=Plot[sc[t],{t,-3T,5T},AxesLabel{"t","sc(t)"},AspectRatio1/3,PlotRangeAll,PlotStyleRed]Show[p2,p1]運行程序可得出以下結(jié)論：信號采樣圖：圖4.4信號采樣圖頻點、振幅、相位分別為振幅頻譜圖：圖4.5振幅頻譜圖相頻圖：圖4.6相頻圖信號波形圖與采樣圖疊加：圖4.7信號波形圖和采樣圖疊加5結(jié)束語本次畢業(yè)設(shè)計至此已經(jīng)接近尾聲，在這幾個月的時間里，我通過利用MATLAB和mathematica強大的數(shù)據(jù)運算功能以及圖像處理功能對于連續(xù)時間的復(fù)頻域分析進行深入的研究。在整個設(shè)計過程中，我首先對于所學(xué)的基礎(chǔ)信號知識進行溫習(xí)鞏固，比如傅立葉級數(shù)、時域采樣、信號頻譜分析等；其次，整個實現(xiàn)過程是通過MATLAB軟件與mathematica軟件完成的，MATLAB的圖形功能十分強大，具有良好的設(shè)計平臺，在此次設(shè)計過程中，我熟練了MATLAB的編程方法，掌握了很多函數(shù)的表示含義及使用方法；mathematica軟件數(shù)學(xué)實現(xiàn)方法特別的強，并且內(nèi)容豐富，操作簡便。最后，通過此次畢業(yè)設(shè)計，我對設(shè)計所用到的軟件有了更加深刻的認(rèn)識，MATLAB與mathematica不僅在數(shù)值計算方面的功能十分強大，而且其圖形仿真功能能夠滿足各個領(lǐng)域的需要，因此MATLAB與mathematica已經(jīng)成為我們工作學(xué)習(xí)中不可或缺的軟件。由于MATLAB與mathematica軟件是專業(yè)性較強的軟件，所以剛開始使用是比較

困難的，通過老師的指導(dǎo)、翻閱了大量的相關(guān)資料，不能說已經(jīng)通悟了MATLAB與mathematica應(yīng)用，但是對一些基本知識還是有所了解的，

在學(xué)習(xí)的過程中學(xué)習(xí)的過程中我們進一步對MATLAB與mathematica編程中的常用語句、過程已經(jīng)初步掌握。這次畢業(yè)設(shè)計，使我懂得了理論與實際相結(jié)合是很重要的，只有理論知識是遠遠不夠的，只有把所學(xué)的理論知識與實踐相結(jié)合起來，從理論中得出結(jié)論，才能真正為社會服務(wù)，從而提高自己的實際動手能力和獨立思考的能力。這次課程設(shè)計終于順利完成，在設(shè)計中遇到的運行和調(diào)試問題，最后在老師的耐心指導(dǎo)下，終于游逆而解。在以后的學(xué)習(xí)過程中我要不段的學(xué)習(xí)，不斷豐富自己的知識。6致謝在此，對畢業(yè)設(shè)計中給