畢業設計說明書基于Matlab的實時聲信號采集與譜分析設計學生姓名：學號：學院：系名：專業：指導教師：2011年6月GatherandanalyzethedesignwiththetableinacousticalsignalbasedonMATLABAbstractInthispaperwilldiscusstheprinciplesoffrequencyanalysis,analysisoftheclassicalandthemodernspectralandtheadvantageanddisadvantageofthem,specificallycomparePeriodogram,BartlettwithWelchalgorithm.AndthenusingMATLABtorealizationofthefrequencyofacousticsignalanalysis,contrastadvangeanddisadvangevisually.ThedatawhichwillbecollectedthroughserialhasreadbackandspectralanalysisundertheenvironmentofMATLAB.Usingthismethodtheacousticsignalacquisitionandspectralanalysissystemhasthestrongandreal-timeprocessingspeedadvantage.Keywords:DSP,spectralanalysis,spectralestimation,andMATLAB基于Matlab的實時聲信號采集與譜分析設計摘要本文介紹了基于MATLAB的聲信號采集與譜分析的設計過程，在闡述頻率分析及功率譜估計原理的基礎上，分析了經典功率譜估計和現代功率譜估計的兩大類算法，主要采用了經典功率譜估計的周期法，間接法和改進的周期法及現代功率譜估計的AR模型參數法對聲信號進行譜分析，并通過大量試驗對各種算法進行驗證對比。該設計利用串口將采集到的聲信號上傳，并在MATLAB環境下進行譜分析，用此方法實現的聲信號采集與譜分析系統，具有實時性強、處理速度快點優點。關鍵詞：DSP，頻譜分析，功率譜估計，MATLAB目錄1緒論 11.1數字信號處理的簡介及發展狀況 11.2頻譜分析的發展現狀及趨勢 11.3本設計的主要任務 32聲信號的采集及傳輸 42.1TMS320F2812的概述 42.2串口傳輸模塊 53譜分析在MATLAB中的實現 73.1MATLAB軟件簡介 73.2譜分析的幾種算法 73.2.1經典功率譜估計的幾種典型算法 73.2.2現代功率譜估計分析方法 163.3用戶界面的設計 203.31窗口界面的實現過程 203.32用戶界面 224總結 25附錄A基于MATLAB用戶界面的譜分析的程序 26參考文獻 33致謝 351緒論1.1數字信號處理的簡介及發展狀況數字信號處理（DigitalSignalProcessing，簡稱DSP）是將信號以數字方式表示并處理的理論和技術。它是一門涉及眾多學科而又廣泛應用于多個領域的新興學科。20世紀60年代以來，隨著計算機和信息技術的飛速發展，數字信號處理技術應運而生并得到迅速的發展。DSP數字信號處理是一種通過使用數學技巧執行轉換或提取信息來處理現實信號的方法。它的核心算法是離散傅立葉變換（DFT），DFT使信號在數字域和頻域都實現了離散化，從而可以用通用計算機處理離散信號。而使數字信號處理從理論走向實用的是快速傅立葉變換（FFT），FFT的出現大大減少了DFT的運算量，使實時的數字信號處理成為可能，極大促進了該學科的發展[1]。DSP（數字信號處理器）是一種高速專用的微處理器，DSP芯片建立在數字信號處理的各種理論和算法基礎上，專門完成各種實時數字信息處理，DSP系統所選用的算法是經過各種實踐檢驗的通用算法的組合和改進而來，它的運算功能強大，專門處理以運算為主，不允許延遲的實時信號；它有特殊的尋址方式，可高效的進行FFY運算；靈活的輸入輸出接口和片內輸入輸出管理；有高速的并行數據處理算法的優化指令集，修改，升級都很方便；靈活的使用C語言或匯編語言編程；集成化程度高，成本低，可靠性好，硬件簡化，有完整的開發和調試工具，開發周期短[2]。其應用領域也由最初的軍用的尖端產品擴展到計算機、通信、家電、辦公自動化、儀器儀表、汽車電子等各個領域，發展速度快，應用實效之大，是目前任何一種器件不能與之相提并論的。在當今的數字化時代快速發展的背景下，DSP已成為通信、計算機、消費類電子產品等領域的基礎器件，被譽為信息社會革命的旗手。DSP將是未來集成電路中發展最快的電子產品，并成為電子產品更新換代的決定因素，它將徹底變革人們的工作、學習和生活方式[2]。1.2頻譜分析的發展現狀及趨勢信號分析主要包括時域分析和頻域分析。譜分析就是頻域分析，具體是指將信號源發出的信號強度按頻率順序展開，使其成為頻率的函數，并考察變化規律。頻譜分析在生產實踐和科學研究中獲得日益廣泛的應用。例如，對汽車、飛機、輪船、汽輪機等各類旋轉機械、電機、機床等機器的主體或部件進行實際運行狀態下的譜分析，可以提供設計數據和檢驗設計效果，或者尋找振源和診斷故障，保證設備的安全運行等；在聲納系統中，為了尋找海洋水面船只或潛艇，需要對噪聲信號進行譜分析，以提供有用信息，判斷艦艇運動速度、方向、位置、大小等。因此對譜分析方法的研究，受到人們的普遍注意和重視，是當前信號處理技術中一個十分活躍的課題。1965年庫利一圖基在《計算數學》雜志上發表快速傅里葉變換（FFT）算法，FFT和頻譜分析很快發展成為機械設備故障診斷、振動分析、無線電通信、信息圖像處理和自動控制等多種學科重要的理論基礎。然而長期的應用和近年來的理論分析表明：經快速傅立葉變換得到的離散頻譜，頻率、幅值和相位均可能產生較大誤差，單諧波加矩形窗時最大誤差從理論上分析可達36.4%；即使加其他窗時，也不能完全消除此影響，在加Hanning窗時，只進行幅值恢復時的最大幅值誤差仍高達15.3%，相位誤差高達90度。因此，頻譜分析的結果在許多領域只能定性而不能精確的定量分析和解決問題，大大限制了該技術的工程應用，特別是在機械振動和故障診斷中的應用受到極大限制。從上世紀70年代中期，有關學者開始致力于頻譜校正理論的研究以期解決離散頻譜誤差較大的問題,通過加窗、局部細化、多點卷積幅值修正等方法來提高頻率識別精度，解決了離散高次諧波參數的精確測量等問題[3]。從目前國內外學者所進行的大量研究工作來看，主要是對單頻率信號（或頻率間隔較大的多頻率信號）離散頻譜的自動識別和校正方法進行探討，密集頻率的校正也只限于兩個鄰近頻率成分的校正，未能深入到連續頻率成分頻譜的誤差和校正方法的研究。而實際工程中的很多信號是密集頻率成分或連續頻率成分的信號，比如“拍振”信號是最簡單的密集頻率信號。在旋轉機械、故障診斷和非線性動力系統分析中，常常也會出現密集頻譜現象。在有限樣本長度下，僅僅由此類信號的FFT頻譜很難識別其頻率構成，也不能確定其各頻率的參數。而且由于旁瓣泄露或主瓣干涉的影響，基于單頻率信號頻譜校正的比值法不再適合此類多頻信號。對此類信號在進行離散頻譜分析時所產生誤差的分析方法與頻率間隔較大的信號誤差分析方法存在巨大差異，校正方法也不相同，校正的難度極大。因此，只有對這類信號在進行離散傅里葉變換時所產生的誤差進行深入系統的分析與研究，并找到一種較完善的頻率、幅值和相位的校正方法，才能使離散傅里葉變換和頻譜分析在機械工程中得到更廣泛應用，同時也擴大其在無線電通信、信息圖像處理、自動控制、多媒體、機械設備故障診斷等技術的應用范圍[4]。當前，具有密集頻譜的頻譜校正問題是目前頻譜校正技術最難解決的問題之一，成為工程界和研究離散頻譜校正的學者們關注的焦點。1.3本設計的主要任務DSP聲DSP聲傳感器譜分析A/DSPI串口通信譜分析本設計主要任務圖1.1本設計的原理框圖本設計的原理框圖如圖1.1所示，首先通過傳感器采集目標聲信號，聲信號進入DSP芯片TMS320F2812內部的AD轉換器實現模數轉換，之后數據通過串口實時上傳至計算機運用Matlab軟件進行頻譜和功率譜分析。本設計的主要任務是通過串口把實時采集的數據傳至計算機后，在MATLAB環境下將數據讀回并實現信號的頻譜分析。本設計主要研究了各種頻譜分析算法實現后的效果圖，直觀比較不同的算法的不同效果。其中包括經典功率譜分析和現代功率譜分析。經典譜估計是將數據工作區外的未知數據假設為零，相當于數據加窗，主要方法有直接法，間接法和改進的直接法。現代譜估計是通過觀測數據估計參數模型再按照求參數模型輸出功率的方法估計信號功率譜，主要是針對經典譜估計的分辨率低和方差性能不好等問題提出的，主要的參數模式是自回歸（AR）模型、移動平均（MA）模型和自回歸/移動平均（ARMA）模型，其中AR模型由線性方程描述，而MA和ARMA模型則由非線性方程描述。由于MA和ARMA模型均可用高階的AR模型來近似，本文使用的是AR參數模型。2聲信號的采集及傳輸2.1TMS320F2812的概述TMS320F2812是美國TI公司推出的C2000平臺上的定點32位DSP芯片，主頻150MHZ、處理性能可達150MIPS，每條指令周期6.67ns。TMS320F2812采用哈佛總線結構，具有統一的存儲模式，包括4M可尋址程序空間和4M可尋址數據空間。同時片內具有128×16位的FLASH存儲器和18K×16位的SRAM，以及4K×16位的引導ROM。最大支持外擴512K×16位的SRAM和512K×16位的FLASH。具有兩個事件管理器（EVA、EVB）以及外設中斷模塊（PIE），最大支持96個外部中斷。TMS320F2812的外部存儲器接口（XINTF）被映射到5個獨立的存儲空間[5]。圖2.1F2812功能組成框圖TMS320F28lx系列DSP的處理器的器件上集成了多種先進的外設，為電機及其他運動控制領域應用的實現提供了良好的平臺。同時代碼和指令與F24x系列數字信號處理器完全兼容，從而保證了項目或產品設計的可延續性[6]。與F24x系列數字信號處理器相比，F2812系列數字信號處理器提高了運算的精度（32位）和系統的處理能力（達到15OMIPS）。該數字信號處理器還集成了128KB的Flash存儲器，4KB的引導ROM，數學運算表以及ZKB的OTPROM，從而大大改善了應用的靈活性。128位的密碼保護機制有效地保護了產品的知識產權。兩個事件管理器模塊為電機及功率變換控制提供了良好的控制功能。16通道高性能12位ADC單元提供了兩個采樣保持電路，可以實現雙通道信號同步采樣。2.2串口傳輸模塊本設計通過TMS320F2812內部的A/D轉換器實現信號的采集，理解串口通信的原理并編寫串口讀數的程序，將采樣結果通過串口上傳到計算機。數據的各位逐位按順序傳送稱為串行通信[7]。串行通信可分為異步傳送和同步傳送兩種方式。異步傳送方式方式采用每個字符按照一個獨立的整體進行發送，字符的間隔時間可以任意變化，即每個字符作為獨立的信息單位（幀），可以隨機地出現在數據流中。所謂“異步”，就是指通信時兩字符之間的間隔事先不能確定，也沒有嚴格的定時要求。在異步通信中，CPU與外部之間的通信遵循以下規定：（1）字符格式。字符格式是指字符的編碼形式及規定。例如，規定每個串行字符由4個部分組成：1個起始位、5～8個數據位、1個奇偶校驗位以及1～2個停止位。（2）傳輸速率。傳輸速率是指每秒鐘傳送的二位進制數，通常稱為波特率。國際上規定了標準波特率系列，最常用的標準波特率是：110、300、600、1200、1800、2400、4800、9600、115200、19200波特等。（3）字符速率。字符速率是指每秒鐘傳送的字符數，它與波特率是兩個相關但表達的意義不相同的概念。例如，若異步通信的數據格式由1位起始位、8位數據位、1位奇偶位、2位停止位組成，波特率為9600b/s，則每秒鐘能夠最多傳送9600/（1+8+1+2）=800個字符[8]。SCI（SerialCommunicationInterface），即串行通信接口，是一個雙線的異步串口，即具有接收和發送兩根信號線的異步串口，一般可以看作是UART（通用異步接收/發送裝置）。SCI發送數據的速度是由波特率來決定的。TMS320F2812的每個SCI都具有兩個8位的波特率寄存器，SCIHBAUD和SCILBAUD，通過編程，可以實現達到64K不同的速率。在進行通信的時候，雙方都必須以相同的數據格式和波特率進行通信，否則通信會失敗[9]。TMS320F2812和PC機上的串口調試軟件進行通信時，TMS320F2812采用了什么樣的數據格式和波特率，那么串口調試軟件也需要設定成相同的數據格式和波特率，反之也一樣。TMS320F2812的SCI模塊支持CPU與采用NRZ（non-return-to-zero不歸零）標準格式的異步外圍設備之間進行數字通信。設計時我們的SCI使用的是RS232串行接口，TMS320F2812就能和其他使用RS232接口的設備進行通信。TMS320F2812內部的兩個SCI之間，或者TMS320F2812的SCI和其他DSP的SCI之間均能實現通信。TMS320F2812內部具有兩個相同的SCI模塊，SCIA和SCIB，每一個SCI模塊都各有一個接收器和發送器。SCI的接收器和發送器各具有一個16級深度的FIFO（Firstinfistout先入先出）隊列，它們還都有自己獨立的使能位和中斷位，可以在半雙工通信中進行獨立的操作，或者在全雙工通信中同時進行操作[10]。根據信息的傳送方向，串行通信可以分為單工、半雙工和全雙工三種[9]。如果在通信過程的任意時刻，信息只能由一方A傳到另一方B，則稱為單工。如果在任意時刻，信息即可由A傳到B，又能由B傳到A，但只能由一個方向上的傳輸存在，稱為半雙工傳輸。如果在任意時刻，線路上存在A到B和B到A的雙信號傳輸，則稱為全雙工。本次設計采用的波特率是115200bps，無校驗，8個數據位，1個停止位。圖2.4顯示了串口調試的界面。圖2.4串口調試窗口3譜分析在MATLAB中的實現在本次設計中要求對從串口接收的數據進行譜分析，由于MATLAB軟件相對其他軟件具有語言簡潔緊湊，庫函數豐富，運算符豐富，功能強勁等優點，所以本設計選擇MATLAB軟件進行譜分析。3.1MATLAB軟件簡介MATLAB的名稱來源自MatrixLaboratory，它是一種科學計算軟件，專門以矩陣的形式處理數據。MATLAB將高性能的數值計算和可視化集成在一起，并提供了大量的內置函數，從而被廣泛應用于科學計算、控件系統、信息處理等領域的分析、仿真和設計工作，而且利用MATLAB產品的開放式結構，可以非常容易地對MATLAB的功能進行擴充，從而在不斷深化對問題認識的同時，不斷完善MATLAB產品以提高產品自身的競爭能力[10]。MATLAB語言是一種基于C語言內核的工程計算語言，集數值分析、矩陣運算、信號處理和圖形顯示于一體，構成了一個使用方便、界面友好的用戶環境。其優點是：可擴展性強,允許用戶自行建立指定功能的M文件、開發自己的工具箱或利用現有的數百種商用Toolbox，易學易用[11]。3.2譜分析的幾種算法信號的頻譜分析是研究信號特性的重要手段之一，對于聲信號，由于它一般是非平穩隨機信號，通常是求其功率譜來進行頻譜分析。功率譜估計（PSD）是用有限長的數據來估計信號的功率譜，它對于認識一個隨機信號或其他應用方面來講都是非常重要的，是數字信號處理的重要研究內容之一。功率譜估計分為經典譜估計和現代譜估計。經典譜估計是將數據工作區外的未知數據假設為零，相當于數據加窗，主要方法有直接法和間接法；現代譜估計是通過觀測數據估計參數模型再按照求參數模型輸出功率的方法估計信號功率譜，主要是針對經典譜估計的分辨率低和方差性能不好等問題提出的，應用最廣的是AR參數模型[12]。本章將分別介紹經典功率譜估計中的直接法、間接法、改進算法和現代功率譜估計中的基于AR模型的幾種相關算法。3.2.1經典功率譜估計的幾種典型算法經典譜估計具有物理概念明確、算法簡單的特點，是目前經常使用的譜估計方法。在經典譜估計中，主要方法有周期圖法、間接法，和直接法的改進算法Bartlett法及Welch法。（1）周期圖法周期圖法又稱直接法，利用該方法得到的隨機信號y（n）的功率譜是直接由傅立葉變換得到的。傅立葉級數是對周期信號求解頻域特性，傅立葉變換則是對非周期信號求解其頻域信息。一個周期信號的傅立葉級數的實質是：把所要研究時域的周期波形分解成許多不同頻率的正弦波的疊加和。傅立葉變換可以看作是時間函數在頻率域上的表示。由傅立葉變換給出的頻率域包含的信息和原函數時間域內包含的完全相同，不同的僅是信息的表示形式。由于計算機的離散性，對一個時間連續信號進行分析要在遵守抽樣定理的前提下，進行抽樣。同樣，對一個時域信號進行分析時，也要在頻域呈離散性，離散傅立葉變換應運而生。綜上所述：周期圖法是把隨機序列y（n）的N個觀測數據視為一個能量有限的序列，直接計算y（n）的離散傅立葉變換得Y（k），然后再取其幅值的平方，并除以N，作為序列y（n）真實功率譜的估計。以下是利用直接法進行功率譜估計的MATLAB代碼：clear;

fs=1000;%采樣頻率

n=0:1/Fs:1;%產生含有噪聲的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

window=boxcar(length(xn));%矩形窗

nfft=1024;

[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法

plot(f,10*log10(Pxx));下面利用直接法對其信號進行分析，圖3.1為原始仿真信號，圖3.2為直接法求得的功率譜圖：圖3.1原始仿真信號圖3.2利用直接法求得的功率譜圖這種將周期圖作為功率譜估計的一個主要缺點就是頻率分辨率低，這是由于周期圖在計算中把觀察到的有限長的N個數據以外的數據認為是零。該方法的優點是計算效率高。在實際的快速傅立葉變換算法的使用中還需要注意數據截斷帶來的譜泄漏問題以及頻率離散帶來的柵欄效應。（2）間接法間接法先由序列x（n）估計出自相關函數R（n），然后對R（n）進行傅立葉變換，便得到x（n）的功率譜估計。以下是利用間接法進行功率譜的估計的MATLAB代碼：

clear;

fs=1000;%采樣頻率

n=0:1/Fs:1;%產生含有噪聲的序列

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;

cxn=xcorr(xn,'unbiased');%計算序列的自相關函數

CXk=fft(cxn,nfft);

Pxx=abs(CXk);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));plot(k,plot_Pxx);下面利用間接法對其信號進行分析，圖3.3為原始仿真信號，圖3.4為間接法求得的功率譜圖：圖3.3原始仿真信號圖3.4利用間接法求得的功率譜圖（3）改進的直接法對于直接法的功率譜估計，當數據長度N太大時，譜曲線起伏加劇，若N太小，譜的分辨率又不好，因此需要改進。而改進的直接法主要有Bartlett算法和Welch算法[13]。

1）Bartlett平均周期圖法Bartlett平均周期圖法的基本思想很簡單，即將具有N個觀測點的可用樣本分成L=N/M個子樣本，每個子樣本有M個觀測點，然后在每個w值上對所有子樣本的周期圖進行平均，一次來減小周期圖中較大的波動[7]。它的數學描述如下，設：,t=1,…,M j=1,…,L（式3.1）表示第j個子樣本的觀測值，并設 （式3.2）表示相應的周期圖。則Bartlett方法為 （式3.3）Bartlett方法是對長度為M的分段數據進行計算，其分辨率大約為1/M。可見，與原始的周期圖方法相比，Bartlett方法的分辨率下降了L倍。可以證明Bartlett方法的方差減小了L倍。選擇M（或L）時可以根據經驗對分辨率和方差進行折中。以下是利用Bartlett進行功率譜的估計的MATLAB代碼：Fs=1000;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));

nfft=1024;

window=boxcar(length(n));%矩形窗

noverlap=0;%數據無重疊

p=0.9;%置信概率

[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);

pause;

figure(2)

plot(k,[plot_Pxxplot_Pxx-plot_Pxxcplot_Pxx+plot_Pxxc]);下面利用Bartlett對其信號進行分析，圖3.5為原始仿真信號，圖3.6為Bartlett求得的功率譜圖：圖3.5原始仿真信號圖3.6利用Bartlett法求得的一信號的功率譜圖基于上述討論有下面的結論：Bartlett與基本周期圖相比，它的分辨率降低了，方差也下降了。2）Welch平滑平均周期圖法目前在工程實際中，經典譜估計獲得廣泛應用的是1967年由Welch提出的修正周期圖法[14]。該方法綜合了Bartlett改進周期圖法的優點，是通過先對數據分段加窗處理，然后再求平均的方法進行的。Welch法對Bartlett法進行了兩方面的修正，一是選擇適當的窗函數w（n），并在周期圖計算前直接加進去。二是在分段時，可使各段之間有重疊,這樣會使方差減小[15]。為了用數學形式描述welch方法[8]，可以先假設：,t=1,...,Mj=1,...,S（式3.4）式（3.4）表示第j個數據段，其中的是第j個觀測序列的起始點。如果K=M，則序列不重疊（但是相鄰連接的），可以采用Bartlett方法對樣本分段（即產生S=L=N/M個數據子樣本）。但是Welch方法中建議K值選為M/2，此時，有個數據段（在連續子段之間有50%的重疊）。以下是利用Welch法進行功率譜的估計的MATLAB代碼：Fs=1000;n=0:1/Fs:1;yn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));nfft=1024;window=boxcar(100);%矩形窗window1=hamming(100);%海明窗window2=blackman(100);%blackman窗noverlap=20;%數據無重疊range='half';%頻率間隔為[0Fs/2],計算一半的頻率[Pxx,f]=pwelch(yn,window,noverlap,nfft,Fs,range);[Pxx1,f]=pwelch(yn,window1,noverlap,nfft,Fs,range);[Pxx2,f]=pwelch(yn,window2,noverlap,nfft,Fs,range);plot_Pxx=10*log10(Pxx);plot_Pxx1=10*log10(Pxx1);plot_Pxx2=10*log10(Pxx2);subplot(311);plot(f,plot_Pxx);title('矩形窗');subplot(312);plot(f,plot_Pxx1);title('海明窗');subplot(313);plot(f,plot_Pxx2);title('blackman窗');下面利用Welch對其信號進行分析，圖3.7為原始仿真信號，圖3.8為Welch得的功率譜圖：圖3.7原始仿真信號圖3.8利用Welch法不同加窗函數的頻譜分析結果通過仿真結果，很容易解釋Welch方法對Bartlett方法改進的原因，因為希望降低被估計的功率譜的方差，所以在數據分段時允許數據之間有重疊，這樣，就增加了（式3.4）式中被平均的周期圖數。同時，因為希望能夠更好地控制被估計的功率譜的偏差/分辨率特性，所以在計算周期圖時引入了時間窗。另外，時間窗對每一個子樣本序列末端數據的加權較小，所以即使相鄰子樣本序列之間有重疊，彼此之間的相關性也較小。這種“去相關”的主要作用是通過式（式3.4）的平均運算完成的。在經典譜估計中，無論是周期圖法還是其改進的方法,都存在著頻率分辨率低、方差性能不好的問題，原因是進行頻譜估計時需要對數據加窗截斷，用有限個數據或其自相關函數來估計無限個數據的功率譜，這其實是假定了窗以外的數據或自相關函數全為零，這種假定是不符合實際的，正是由于這些不符合實際的假設造成了經典譜估計分辨率較差[15]。另外，經典譜估計的功率譜定義中既無求均值運算又無求極限運算，因而使得譜估計的方差性能較差，當數據很短時，這個問題更為突出，所以大家又熱衷于模型參數估計法的研究。3.2.2現代功率譜估計分析方法用模型法進行譜估計的方法稱為現代功率譜估計方法，它主要是針對經典譜估計的分辨率低和方差性能不好的問題而提出的，它不是直接地進行功率譜的計算，而是假設隨機信號服從一個模型，通過模型參數得到信號的功率譜。現代功率譜估計的內容極其豐富,涉及的學科和領域也相當廣泛,大致可分為參數模型估計和非參數模型估計兩大類，前者有AR模型、MA模型、ARMA模型、PRONY指數模型等；后者有最小方差方法、多分量的MUSIC方法等。本小節針對AR模型參數的幾種典型求解算法原理進行簡單介紹[17]。AR模型的物理意義是認為序列y(n)是白噪聲作用于線性時不變系統時的系統響應。假設一個隨機過程可以由AR(p)模型刻畫[18]，即 （式3.5）式中：u(n)、x(n)為實平穩的隨機信號；u(n)為白噪聲；R為方差。經過推導，可得：=+（式3.6）式中：為的自相關函數。上式是AR模型的正則方程Yule-Walke。解Yule-Walke方程是一個復雜的矩陣求逆數學問題，當N大時運算量會很大。因此就有了各種簡便的算法，如Levinson-Durbin遞推算法、協方差算法、Burg遞推算法、自相關算法等[19]。（1）AR模型譜估計的協方差算法協方差即數據矩陣的協方差方法。觀察數據窗[0，N-1]，為了不對觀察數據窗之外的數據任何假設，將這組數據分成初始數據并記為，其它數據記為。計算協方差算法可以分為三步：第一，確定模型階p,由觀察數據共軛轉置構成數據矩陣；第二，對進行特征分解，得到W個不為零的特征值和部分特征矩陣；第三，將估計參數帶入功率譜密度公式(式3.7)。 （式3.7）以下是利用改進協方差算法進行功率譜的估計的MATLAB代碼：fs=1000;%取樣頻率n=0:1/fs:1;xn=4.1*cos(2*pi*300*n)+3.8*cos(2*pi*305*n)+0.5*randn(size(n));nfft=1024;%FFT的長度subplot(411);plot(n,xn);title('仿真信號x(n)')[px3,f3]=pmcov(xn,30,nfft,fs);subplot(412);plot(f3,px3);title('改進協方差算法功率譜估計')下面利用改進協方差算法對其信號進行分析，圖3.9為原始仿真信號和改進協方差算法求得的功率譜圖：圖3.9原始仿真信號和改進協方差算法功率譜圖（2）AR模型譜估計的自相關算法與協方差算法不同，由觀測數據對信號在盡可能多的時刻進行預測，產生盡可能多的預測誤差，是自相關算法[20]。計算協方差算法：第一，利用觀測數據，得到預測的最大時間范圍是[0，N+P-1]；第二，用觀測數據以外的數據，令，則（式3.8）計算該范圍的預測誤差和；第三，由得到自相關矩陣，即Yule-walker方程，利用levinson遞歸求解Yule-walker方程得到AR模型的參數等效于前向預測器的系數。以下是利用自相關法進行功率譜的估計的MATLAB代碼：fs=1000;%取樣頻率n=0:1/fs:1;xn=4.1*cos(2*pi*300*n)+3.8*cos(2*pi*305*n)+0.5*randn(size(n));nfft=2048;%FFT的長度subplot(411);plot(n,xn);title('仿真信號x(n)')[px1,f1]=pyulear(xn,30,nfft,fs);subplot(413);plot(f1,px1);title('自相關法功率譜估計')下面利用自相關算法對其信號進行分析，圖3.10為原始仿真信號和自相關算法求得的功率譜圖：圖3.10原始仿真信號和自相關算法功率譜圖（3）AR模型譜估計的Burg算法在現代功率譜估計的歷史中，Burg提出的最大熵方法產生了很大的影響，刺激了現代譜估計方法的發展。Burg算法進行功率譜估計時是令前后向預測誤差功率之和最小,即對,前后都不加窗,使用Levinson-Durbin遞推可快速的求解AR系數。Burg算法是建立在數據基礎之上的，避免了先計算自相關函數，從而提高計算速度，計算不太復雜，且分辨率優于自相關法，是較為通用的方法。以下是利用Burg法進行功率譜的估計的MATLAB代碼：fs=1000;%取樣頻率n=0:1/fs:1;xn=4.1*cos(2*pi*300*n)+3.8*cos(2*pi*305*n)+0.5*randn(size(n));nfft=2048;%FFT的長度subplot(411);plot(n,xn);title('仿真信號x(n)')[px2,f2]=pburg(xn,30,nfft,fs);subplot(414);plot(f2,px2);title('Burg算法功率譜估計')下面利用Burg算法對其信號進行分析，圖3.11為原始仿真信號和Burg算法求得的功率譜圖：圖3.11原始信號和Burg算法功率譜圖仿真結果直觀地說明了自相關算法、Burg算法和改進協方差算法各自的優缺點。對相隔僅有5HZ的2個混合信號進行譜估計時，自相關法不容易看出其頻率成份，而Burg算法和改進協方差算法提高了參數估計的精度和頻率分辨率。參數模型譜估計方法是現代譜估計的重要內容，AR模型譜估計隱含著數據和自相關函數的外推，其長度可能超過給定的長度，分辨率不受信源信號長度的限制，這是經典譜估計無法做到的。在AR模型自相關算法中由于對ef（n）前后都加窗，使得自相關法的分辨率降低，數據越短分辨率越差。3.3用戶界面的設計圖形用戶界面是指由窗口、菜單、對話框等各種元素組成的用戶界面。在這種用戶界面中，用戶的操作既形象生動，又方便靈活，所以當今絕大部分開發環境與應用程序都采用圖形用戶界面，許多流行的開發工具也都可以進行圖形用戶界面的設計。MATLAB作為功能強大的科學計算軟件，同樣也提供了一套可視化的創建圖形窗口的工具，本設計利用MATLAB所提供的圖像用戶界面（GUI），實現一個可視的面向對象的操作界面[21]。3.31窗口界面的實現過程該系統界面的設計主要是利用MATLAB所提供的GUI（GraphUserInterface）向導設計控件而完成的，該向導可以實現多種控件的設計，給用戶提供了一種友好的交互方式，同時也給操作帶來很多方便。圖形用戶界面GUI是包含圖形對象（如圖形窗口、菜單、控件、文本）的用戶界面，用戶以某種方式選擇或者激活這些對象會發生變化或引起動作。（1）啟動GUI的方法啟動圖形用戶界面的方法有很多種，例如可以利用工具欄上的命令按鈕，也可以采用菜單和命令的方式。在本人設計的案例中采用的是命令方式：在命令窗口直接輸入guide命令，彈出的窗口如圖3.12所示：圖3.12為啟動GUI界面彈出的窗口在彈出的guidequickstart窗口中選擇createnewgui選項卡里面地BlankGUI選項，這樣就可以進入圖形用戶界面來設計我們的系統靜態界面，如下圖3.13所示：圖3.13為系統靜態界面上圖就是Guide提供的圖形界面設計工具集，在此界面下就可以利用控件組件、文本菜單、排列工具等對系統的界面進行設計。當靜態界面設計完成以后，對該界面進行保存，此時Guidie將自動生成兩個發布文件，分布是.fig文件和.m文件。Fig文件：該文件包括圖形窗口及其所有后裔的完全描述，包括所有對象的屬性值。Fig文件是一個二進制文件，調用hgsave命令或界面設計編輯器的file菜單save選項保存圖像窗口時將產生該文件。Fig文件最有用的地方之一就是對象句柄的保存和引用，可以使用open、openfig和hgload命令打開一個后綴為.fig的文件。M文件：該文件包括GUI設計、控制函數以及為子函數的用戶控件回調函數，主要用于控制GUI展開時的各種特征。這個M文件可以分為GUI初始化和回調函數兩個部分，用戶控件的回調函數根據用戶與GUI的具體交互方式分別調用。3.32用戶界面本設計中用戶界面，大體分為信號輸入模塊，原始信號顯示模塊，頻譜分析模塊，經典譜分析模塊和現代譜分析模塊四部分。經典譜分析模塊包含直接法、間接法、Welch法和Bartlett法四個小模塊。現代譜分析模塊則包括Burg法和自相關法兩個小模塊。根據不同的信號來源，界面主要采用兩種不同輸入方式。一種信號來源是由串口實時采集到的信號，把該信號作為一種輸入方式。然后點擊打開串口，串口開始讀取數據同時把數據保存在D盤下的一個文件中。接著關閉串口，點擊打開按鈕，即可讀回由串口實時采集到的聲信號數據。打開文件的同時首先在坐標軸中會顯示原始信號波形。其次點擊頻譜分析按鈕，坐標軸會顯示信號的幅度譜。最后就是功率譜分析部分，用到了八種方法，分別是經典功率譜估計中的直接法、間接法、welch法、barlett法，以及現代功率譜估計中AR模型法的Burg算法和自相關算法。圖3.14為輸入信號為串口接收的數據的界面。另一種信號來源是打開mat文件中已存的兩個不同頻率疊加的聲音信號，把該信號作為一種輸入方式。然后點打開按鈕，在打開文件名的同時首先會在坐標軸中顯示原始信號波形。其他的操作如上第一種方式。圖3.15為輸入信號為mat文件中已存的信號。圖3.14輸入信號為串口接收的數據界面圖3.15輸入信號為mat文件已存的數據通過實驗仿真可以直觀地看出以下特性：（1）經典功率譜估計中的周期圖法所得到的結果特點是離散性大，曲線粗糙，方差較大，但是分辨率較高。（2）Bartlett平均周期圖法和Welch平滑平均周期圖法的收斂性較好，曲線平滑，估計的結果方差較小，但是功率譜主瓣較寬，分辨率低。這是由于對隨機序列的分段處理引起了長度有限所帶來的Gibbs現象而造成的。（3）與Bartlett法相比，Welch法的估計曲線比較粗糙，但是分辨率較好，原因是Welch法中對數據進行截斷時加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相對于矩形窗，窗的主瓣包含更多的能量，因而使功率譜的主瓣較窄，分辨率較高。經典功率譜估計的分辨率反比于有效信號的長度，但現代譜估計的分辨率可以不受此限制。這是因為對于給定的N點有限長序列x（n），雖然其估計出的自相關函數也是有限長的，但是現代譜估計的一些隱含著數據和自相關函數的外推，使其可能的長度超過給定的長度，不象經典譜估計那樣受窗函數的影響。因而現代譜的分辨率比較高，而且現代譜線要平滑得多，從實驗仿真圖可以清楚看出[22]。4總結本設計首先通過傳感器采集目標聲信號，聲信號進入DSP芯片TMS320F2812內部的AD轉換器實現模數轉換，之后數據通過串口實時上傳至計算機運用Matlab軟件進行頻譜和功率譜分析。主要是對于信號的頻率分析算法進行研究，所謂的頻率分析即是在頻域范圍內研究信號的特征，通常用求信號的功率譜來研究信號的頻率特性。首先是介紹經典功率譜估計和現代功率譜估計幾種常見的頻譜分析算法的原理，并在Matlab軟件上利用用戶界面對采集的原始聲信號數據進行頻譜分析，則直觀的對比出每種算法各自的優缺點。而有待改進的地方是在MATLAB環境下如何采用更多的算法進行頻譜分析，如何編寫更美觀實用的界面，也可以進一步研究如何提高頻譜分析的分辨率，降低估計方差等。附錄A基于MATLAB用戶界面的譜分析的程序functionvarargout=untitledzuizhong(varargin)gui_Singleton=1;gui_State=struct('gui_Name',mfilename,...'gui_Singleton',gui_Singleton,...'gui_OpeningFcn',@bishe_OpeningFcn,...'gui_OutputFcn',@bishe_OutputFcn,...'gui_LayoutFcn',[],...'gui_Callback',[]);ifnargin&&ischar(varargin{1})gui_State.gui_Callback=str2func(varargin{1});endifnargout[varargout{1:nargout}]=gui_mainfcn(gui_State,varargin{:});elsegui_mainfcn(gui_State,varargin{:});endfunctionbishe_OpeningFcn(hObject,eventdata,handles,varargin)handles.output=hObject;guidata(hObject,handles);functionvarargout=untitledzuizhong_OutputFcn(hObject,eventdata,handles)varargout{1}=handles.output;functiontxtfile_Callback(hObject,eventdata,handles)%打開txt文件，單擊單選按鈕，相應參數設置部分可用if~get(handles.txtfile,'value')set(handles.filename_01,'enable','off');set(handles.open_txt,'enable','off');elseset(handles.filename_01,'enable','on');set(handles.open_txt,'enable','on');set(handles.wavefile,'value',0);set(handles.filename_02,'enable','off');set(handles.open_wave,'enable','off');endfunctionwavefile_Callback(hObject,eventdata,handles)%打開wav文件，單擊單選按鈕，相應參數設置部分可用if~get(handles.wavefile,'value')set(handles.filename_02,'enable','off');set(handles.open_wave,'enable','off');elseset(handles.filename_02,'enable','on');set(handles.open_wave,'enable','on');set(handles.txtfile,'value',0);set(handles.filename_01,'enable','off');set(handles.open_txt,'enable','off');endfunctionopen_wave_Callback(hObject,eventdata,handles)[filename,filepath]=uigetfile('.wav','Openwavfile');[y,fs,nbits]=wavread([filepath,filename]);set(handles.filename_02,'string',filename);handles.y=y;handles.fs=fs;plot(handles.axes1,handles.y);guidata(hObject,handles);functionfilename_02_Callback(hObject,eventdata,handles)functionfilename_02_CreateFcn(hObject,eventdata,handles)ifispc&&isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunctionpinpu_Callback(hObject,eventdata,handles)%進行FFT變換并做頻譜圖Y=fft(handles.y);%進行fft變換mag=abs(Y);%求幅值f=(0:length(Y)-1)'*handles.fs/length(Y);%進行對應的頻率轉換plot(handles.xianshi2,f,fuzhi);%做幅值譜圖functionzhijie_Callback(hObject,eventdata,handles)functionradiobutton3_Callback(hObject,eventdata,handles)functionradiobutton4_Callback(hObject,eventdata,handles)functionpushbutton4_Callback(hObject,eventdata,handles)Fs=22050;window=boxcar(length(handles.y));%矩形窗nfft=1024;[Pxx,f]=periodogram(handles.y,window,nfft,Fs);%直接法plot(handles.xianshi4,f,10*log10(Pxx));functionjianjie_Callback(hObject,eventdata,handles)Fs=1000;nfft=1024;cxn=xcorr(handles.y,'unbiased');CXK=fft(cxn,nfft);Pxx=abs(CXK);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot(handles.xianshi10,k,plot_Pxx);functionwelch_Callback(hObject,eventdata,handles)Fs=22050;nfft=1024;window=boxcar(100);%矩形窗noverlap=20;%數據無重疊range='half';%頻率間隔為[0Fs/2]，只計算一半的頻率[Pxx,f]=pwelch(handles.y,window,noverlap,nfft,Fs,range);plot_Pxx=10*log10(Pxx);plot(handles.xianshi5,f,plot_Pxx);functionbartlett_Callback(hObject,eventdata,handles)Fs=22050;nfft=1024;window=boxcar(100);%矩形窗noverlap=0;%數據無重疊p=0.9;%置信概率[Pxx,Pxxc]=psd(handles.y,nfft,Fs,window,noverlap,p);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot(handles.xinashi7,k,plot_Pxx);functionAR_Callback(hObject,eventdata,handles)fs=1000;nfft=2048;[px1,f1]=pburg(handles.y,50,nfft,fs);plot(handles.xianshi8,f1,px1);functionzixinagguan_Callback(hObject,eventdata,handles)fs=1000;nfft=2048;[px2,f2]=pyulear(handles.y,30,nfft,fs);plot(handles.xianshi11,f2,px2);functionopen_txt_Callback(hObject,eventdata,handles)filepath=get(handles.filename_01,'string')fid=fopen(filepath,'rt');[a,count]=fscanf(fid,'%X',inf);fori=1:count/2-2b(i)=a(2*i-1)*256+a(2*i);endhandles.y=b;handles.fs=2000;plot(handles.axes1,b);guidata(hObject,handles);functionfilename_01_Callback(hObject,eventdata,handles)functionfilename_01_CreateFcn(hObject,eventdata,handles)ifispc&&isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunctionxianshi_Callback(hObject,eventdata,handles)functionxianshi_CreateFcn(hObject,eventdata,handles)ifispc&&isequal(get(hObject,'BackgroundColor'),get(0,'defaultUicontrolBackgroundColor'))set(hObject,'BackgroundColor','white');endfunctionpushbutton10_Callback(hObject,eventdata,handles)functionstart_serial_Callback(hObject,eventdata,handles)globalscom;globalfid;ifget(hObject,'value')set(handles.start_serial,'string','關閉串口');scom=serial('com1');set(s