基于MATLAB的信號(hào)分析

1、摘 要本文首先介紹了三種典型數(shù)字信號(hào)，對(duì)離散信號(hào)的均值、方差、相關(guān)和高斯隨機(jī)信號(hào)的統(tǒng)計(jì)特性用MATLAB仿真和分析，用MATLAB實(shí)現(xiàn)離散信號(hào)的加減運(yùn)算。其次編程實(shí)現(xiàn)了三種典型離散信號(hào)的離散傅里葉變換，顯示時(shí)域信號(hào)和頻譜圖形（幅度值和相位譜），最后用經(jīng)典功率譜估計(jì)中的周期圖估計(jì)法、Bartlett譜估計(jì)法及Welch譜估計(jì)法，對(duì)正弦序列加高斯隨機(jī)序列進(jìn)行功率譜估計(jì)，并且用時(shí)域提取法進(jìn)行提取。關(guān)鍵詞：離散信號(hào) 功率譜估計(jì) 提取目 錄摘 要1前 言1一 基本原理1二 系統(tǒng)分析3三 詳細(xì)設(shè)計(jì)4總 結(jié)16附 錄17參考文獻(xiàn)27致 謝280前 言 數(shù)字信號(hào)處理無論在理論上還是技術(shù)上都有了突破性的發(fā)展。它

2、在工業(yè)中的成功是由于開發(fā)和利用了廉價(jià)的硬件和軟件，現(xiàn)在各個(gè)領(lǐng)域的新工藝和新應(yīng)用中現(xiàn)在都想利用DSP算法，這就導(dǎo)致了對(duì)具有DSP基礎(chǔ)的電氣工程師的巨大需求，同時(shí)也對(duì)DSP的電子仿真提出了更高的要求。 隨著MATLAB的出現(xiàn)和不斷完善，尤其是MATLAB的信號(hào)分析工具箱的提出，越來越多的電氣工程師門已經(jīng)意識(shí)到用MATLAB來解決電信工程的實(shí)際問題是一種省事又省力的選擇，由于MATLAB在信號(hào)處理與分析等方面有很大的優(yōu)點(diǎn)，在過去的10年中，MATLAB已經(jīng)成為DSP應(yīng)用中分析和設(shè)計(jì)的主要仿真工具。MATLAB是一個(gè)高性能的科學(xué)計(jì)算軟件，廣泛應(yīng)用于數(shù)學(xué)計(jì)算、算法開發(fā)、數(shù)學(xué)建模、系統(tǒng)仿真、數(shù)據(jù)分析處理及

3、可視化、科學(xué)和工程測繪、應(yīng)用系統(tǒng)開發(fā)等。當(dāng)前它的使用范圍涵蓋了工業(yè)、電子、醫(yī)學(xué)、醫(yī)療、建筑等領(lǐng)域。MATLAB用于數(shù)字信號(hào)處理與分析仿真實(shí)驗(yàn)平臺(tái)的設(shè)計(jì)已經(jīng)十分普遍。一 基本原理1.1離散時(shí)間信號(hào)（離散序列的時(shí)域分析）信號(hào)可以表示為一個(gè)或多個(gè)自變量的函數(shù)。在信號(hào)的分類中，自變量連續(xù)變化的信號(hào)稱為連續(xù)時(shí)間（模擬）信號(hào)，自變量離散變化的信號(hào)稱為離散時(shí)間信號(hào)。一個(gè)連續(xù)（模擬）信號(hào)通常用表示其中變量t可代表任何物理量，一般都假定它代表時(shí)間，以秒計(jì)。一個(gè)離散信號(hào)通常用表示，其中變量n是整數(shù)值并在時(shí)間上代表一些離散時(shí)刻；因此，它也稱作離散時(shí)間信號(hào)。離散時(shí)間信號(hào)通常稱為(離散)序列。在通用的計(jì)算機(jī)或?qū)Ｓ玫臄?shù)字

4、信號(hào)處理器中，只能用有限位二進(jìn)制數(shù)來表示離散時(shí)間信號(hào)的幅度，即離散時(shí)間信號(hào)的幅度只能取有限個(gè)離散值。這類在幅度上只能取有限個(gè)離散值的離散時(shí)間信號(hào)稱為數(shù)字信號(hào)。在分析離散時(shí)間信號(hào)時(shí)，常將復(fù)雜的離散時(shí)間信號(hào)表示為基本離散序列的線性組合，這些基本離散時(shí)間信號(hào)也是實(shí)際應(yīng)用中常用的序列。因此基本離散序列在離散時(shí)間信號(hào)分析中起著非常重要的作用。1.2傅里葉變換（頻域分析）離散傅里葉變換（Discrete Fourier Transform，縮寫為DFT），是傅里葉變換在時(shí)域和頻域上都呈離散的形式，將信號(hào)的時(shí)域采樣變換為其DTFT的頻域采樣。在形式上，變換兩端（時(shí)域和頻域上）的序列是有限長的，而實(shí)際上這兩組

5、序列都應(yīng)當(dāng)被認(rèn)為是離散周期信號(hào)的主值序列。即使對(duì)有限長的離散信號(hào)作DFT，也應(yīng)當(dāng)將其看作其周期延拓的變換。在實(shí)際應(yīng)用中通常采用快速傅里葉變換計(jì)算DFT。在各種信號(hào)序列中，有限長序列占重要地位。對(duì)有限長序列可以利用離散傅立葉變換(DFT)進(jìn)行分析。DFT不但可以很好的反映序列的頻譜特性，而且易于用快速算法(FFT)在計(jì)算機(jī)上進(jìn)行分析。有限長序列的DFT是其z變換在單位圓上的等距離采樣，或者說是序列傅立葉的等距離采樣，因此可以用于序列的譜分析。1.3數(shù)字濾波器及濾波濾波器的作用是利用離散時(shí)間系統(tǒng)的特性對(duì)輸入信號(hào)波形(或頻譜)進(jìn)行加工處理，或者說利用數(shù)字方法按預(yù)定的要求對(duì)信號(hào)進(jìn)行變換。數(shù)字濾波器,

6、是數(shù)字信號(hào)處理中及其重要的一部分。隨著信息時(shí)代和數(shù)字技術(shù)的發(fā)展，受到人們?cè)絹碓蕉嗟闹匾?。?shù)字濾波器可以通過數(shù)值運(yùn)算實(shí)現(xiàn)濾波，所以數(shù)字濾波器處理精度高、穩(wěn)定、體積小、重量輕、靈活不存在阻抗匹配問題，可以實(shí)現(xiàn)模擬濾波器無法實(shí)現(xiàn)的特殊功能。數(shù)字濾波器種類很多，根據(jù)其實(shí)現(xiàn)的網(wǎng)絡(luò)結(jié)構(gòu)或者其沖激響應(yīng)函數(shù)的時(shí)域特性，可分為兩種，即有限沖激響應(yīng)( FIR，F(xiàn)inite Impulse Response)濾波器和無限沖激響應(yīng)( IIR，Infinite Impulse Response)濾波器。28二 系統(tǒng)分析2.1 系統(tǒng)設(shè)計(jì)流程：（1） 明確典型數(shù)字信號(hào)在MATLAB中的產(chǎn)生方法，產(chǎn)生正弦信號(hào)、周期信號(hào)和隨機(jī)

7、信號(hào)。對(duì)確知信號(hào)（正弦、周期）進(jìn)行頻域分析，對(duì)隨機(jī)信號(hào)進(jìn)行統(tǒng)計(jì)分析（均值、方差、相關(guān)）和頻域分析。（2） 分別進(jìn)行正弦信號(hào)與隨機(jī)信號(hào)、周期信號(hào)與隨機(jī)信號(hào)的加減運(yùn)算，并對(duì)正弦信號(hào)與隨機(jī)信號(hào)組成的混合信號(hào)進(jìn)行功率譜分析。系統(tǒng)初始化（3）設(shè)計(jì)數(shù)字低通濾波器對(duì)上述產(chǎn)生的混合信號(hào)進(jìn)行濾噪處理，恢復(fù)原始正弦信號(hào)產(chǎn)生數(shù)字信號(hào)是否為隨機(jī)信號(hào)？統(tǒng)計(jì)分析、頻域分析頻域分析數(shù)字信號(hào)的加減運(yùn)算對(duì)正弦余隨機(jī)組成的混合信號(hào)進(jìn)行功率譜分析設(shè)計(jì)濾波器對(duì)混合信號(hào)進(jìn)行濾波恢復(fù)原始正弦信號(hào)三 詳細(xì)設(shè)計(jì)3.1典型數(shù)字信號(hào)的產(chǎn)生與分析3.1.1正弦序列 （3-1） 其中是幅度,而是相位（弧度rad）?？捎肕ATLAB函數(shù)cos（或si

8、n）產(chǎn)生正弦序列。 下圖是正弦信號(hào)的正弦序列、正弦序列時(shí)域信號(hào)及正弦序列幅頻、相頻特性圖，運(yùn)行MATLAB得到其正弦信號(hào)的均值為0.0000 ,方差為2.0328圖3-1正弦序列 正弦序列在點(diǎn)DFT，正弦序列的幅頻、相頻特性圖圖3-2正弦序列幅頻、相頻特性3.1.2 周期序列 如果有序列則序列是周期的。滿足上面關(guān)系的最小整數(shù)成為基波周期?，F(xiàn)用表示一個(gè)周期序列。為了從一個(gè)周期的產(chǎn)生的個(gè)周期，可將重復(fù)次。 下圖是周期序列及幅頻、相頻特性圖，其均值為1.5000 方差為1.2500圖3-3 周期序列 圖3-4周期序列的幅頻、相頻特性3.1.3 高斯隨機(jī)序列 利用MATLAB很容易產(chǎn)生兩類隨機(jī)信號(hào)：

9、Rand(1,N)在0.1區(qū)間上產(chǎn)生N點(diǎn)均勻分布的隨機(jī)序列； Randn(1,N)產(chǎn)生均值為0，方差為1的高斯隨機(jī)序列，也就是白噪聲序列。 圖3-5正弦信號(hào)加高斯隨機(jī)信號(hào)自相關(guān)3.2 離散信號(hào)進(jìn)行加減運(yùn)算 若信號(hào)： （3-2）值得注意的是，當(dāng)序列和的長度不等或位置不對(duì)應(yīng)時(shí)，首先應(yīng)使兩者的位置對(duì)齊，然后通過zeros函數(shù)左右補(bǔ)零使其長度相等后再相加。下圖是正弦序列加減隨機(jī)序列圖：圖3-6離散信號(hào)加減運(yùn)算3.3 經(jīng)典法功率譜估計(jì) （3-3）經(jīng)典功率譜估計(jì)主要包括Tukey提出的相關(guān)函數(shù)估計(jì)，Schuster提出的周期圖估計(jì)，間接法（BT法）、Bartlett法譜估計(jì)、Welch法譜估計(jì)、Butta

10、ll法譜估計(jì)，以及Brurg提出的最大熵譜估計(jì)。本節(jié)主要運(yùn)用基于直接法的功率譜估計(jì)、基于間接法的功率譜估計(jì)、Welch法功率譜估計(jì)對(duì)正弦信號(hào)加高斯噪聲信號(hào)的功率譜估計(jì)。3.3.1基于直接法的功率譜估計(jì) 直接法是Schuster于1899年提出的，又稱周期圖法。它將平穩(wěn)信號(hào)序列的N點(diǎn)觀察數(shù)據(jù)視為能量有限的信號(hào)，直接取的傅里葉變換，然后再取其幅值的平方，并除以數(shù)據(jù)長度N，作為序列功率譜的估計(jì)。 假設(shè)直接估計(jì)出的功率譜為，則： 直接法之所以得到廣泛使用，是由于它于序列的頻譜有對(duì)應(yīng)關(guān)系，可以采用FFT算法來快速計(jì)算。但是在直接法功率譜估計(jì)中，對(duì)無限長的平穩(wěn)信號(hào)序列進(jìn)行截?cái)?，這等于對(duì)無限長的序列加以矩形

11、窗，使之變成有限長的數(shù)據(jù)。這也意味著對(duì)自相關(guān)函數(shù)的加窗，使得功率譜與窗函數(shù)的卷積。這種頻域卷積會(huì)產(chǎn)生頻譜泄露，容易使弱信號(hào)的主瓣被強(qiáng)信號(hào)的旁瓣淹沒，造成頻譜的模糊和失真，使得周期圖功率譜的分辨率較低。 在MATLAB中可以利用periodogram實(shí)現(xiàn)直接法的功率譜估計(jì)，其調(diào)用格式為： Pxx,w=periodogram(x) Pxx,w=periodogram(x,window) Pxx,w=periodogram(x,window,nfft,fs)Pxx,=periodogram(x,range)Periodogram()Periodogram函數(shù)參數(shù)的說明如下：x為進(jìn)行功率譜估計(jì)的輸入有

12、限長序列。Window用于指定采用的窗函數(shù)，缺省時(shí)采用矩形窗，窗長等于輸入序列x的長度。對(duì)于周期圖法，算法未指定窗函數(shù)，但有限長序列相當(dāng)于采用矩形窗。nfft設(shè)定FFT算法的長度，一般取2的整次冪，缺省值為256。fs為采樣頻率，缺省值為1.range用于指定頻率的取值范圍。若range=twosided,它是復(fù)序列的默認(rèn)選項(xiàng)，其頻率取值范圍為0,fs,但是如果fs為空向量時(shí)，則頻率取值范圍為0,1,如果沒有指定fs,則頻率取值范圍為0,2;若range=onesided，它是實(shí)序列的默認(rèn)選項(xiàng)，其頻率范圍為0,fs/2。Pxx為輸出的功率譜估計(jì)值。f為頻率向量。w為歸一化的頻率向量。 下圖用公

13、式法和periodogram實(shí)現(xiàn)直接法的功率譜估計(jì)曲線，由圖可以看出利用公式直接計(jì)算的序列的功率譜與利用函數(shù)periodogram計(jì)算得到的功率譜完全一致。 圖3-7正弦序列加高斯隨機(jī)序列圖3-8公式直接計(jì)算的功率譜曲線圖3-9 periodogram函數(shù)計(jì)算的功率譜估計(jì)曲線3.4 基于改進(jìn)直接法的功率譜估計(jì) 對(duì)于功率譜直接法估計(jì)，它是通過有限長樣本序列的FFT變換來計(jì)算平穩(wěn)信號(hào)序列功率譜的，當(dāng)數(shù)據(jù)長度N太大時(shí)，譜曲線起伏加劇，當(dāng)數(shù)據(jù)N長度太小，譜的分辨率又不好。因此必須對(duì)周期圖法的功率譜估計(jì)方法進(jìn)行改進(jìn)，這里的改進(jìn)指的是數(shù)據(jù)方差的改進(jìn)。當(dāng)前主要有兩中改進(jìn)方法：一種是周期圖的平滑，即采用間接法

14、估計(jì)功率譜；另一種是平均法，它將長度為N的數(shù)據(jù)分成L段，分別求出每段的功率譜，然后加以平均。下面介紹直接法的改進(jìn)方法Batlett法、Welch法。3.4.1 Batlett法Batlett平均周期圖法的方法是將N點(diǎn)的有限長序列，分成互不重疊的L段數(shù)據(jù)，每段有M個(gè)樣本，則有LM=N，即： （3-4）其中，為長度是M的矩形窗。然后求取每一段的M個(gè)樣本的功率譜估計(jì)，即： （3-5）最后求出所有M段數(shù)據(jù)功率譜的平均值，即： （3-6）并將上式作為整個(gè)序列的功率譜估計(jì)。隨著L的增大，平均周期圖的方差趨于零，因此它是功率譜的漸進(jìn)一致估計(jì)。雖然分段平均周期圖法功率譜估計(jì)可以減小估計(jì)誤差和波動(dòng)，但是由于這種

15、方法將長信號(hào)分段成短信號(hào)，從而使功率譜的分辨率下降。在MATLAB中可以利用函數(shù)psd來實(shí)現(xiàn)Barlett平均周期圖法的功率譜估計(jì)，其調(diào)用格式為：Pxx=psd(x,nfft,fs,window)Pxx,f=psd(x,nfft,fs,window,noverlap)Pxx,f,F=psd(x,nfft,fs,window,noverlap,p)Psd(x,.,dflag) 在使用psd函數(shù)計(jì)算Bartlett平均周期圖功率譜估計(jì)時(shí)，還應(yīng)注意以下幾點(diǎn)：(1) 由于Bartlett方法中沒有指定窗函數(shù)，但對(duì)于有限長序列來說，相當(dāng)于采用矩形窗，所以參數(shù)window應(yīng)設(shè)為boxcar。(2) Bar

16、tlett方法中沒有指定數(shù)據(jù)重疊，因而參數(shù)noverlap應(yīng)設(shè)為0；(3) 利用psd函數(shù)實(shí)現(xiàn)Bartlett方法，其分段數(shù)L默認(rèn)為： （3-7） 其中fix表示數(shù)值朝零方向取整； (4)參數(shù)p為置信概率； (5)參數(shù)Pxx為Pxx的p100%置信區(qū)間估計(jì)。圖3-10序列的Bartlett法功率譜估計(jì)曲線3.4.2 Welch法Welch法，又稱加權(quán)交疊平均法，它是對(duì)Bartlett法的改進(jìn)，主要體現(xiàn)在兩個(gè)方面： 在對(duì)序列分段時(shí)，允許每段數(shù)據(jù)有部分交疊； 每段數(shù)據(jù)可以選擇其他的窗函數(shù)，不一定要求是矩形窗。在每段數(shù)據(jù)的功率譜時(shí)，Welch法的算法與Bartlett法相同，即： （3-8）其中M為

17、每段數(shù)據(jù)樣本， 它保證由Welch方法得到的功率譜估計(jì)是無偏估計(jì)。由Welch法得到的平均功率譜為： （3-9）因?yàn)閃elch法允許各數(shù)據(jù)段交疊，因而數(shù)據(jù)段L增加，這樣方差得到更大的改善。但是數(shù)據(jù)的交疊又減小了每一段數(shù)據(jù)的不相關(guān)性，使方差的減小不會(huì)達(dá)到理論計(jì)算的程度。盡管如此，由Welch法估計(jì)出來的功率譜是漸進(jìn)無偏的。由于Welch法在估計(jì)序列的功率時(shí)，使用了窗函數(shù)，由窗函數(shù)基本知識(shí)可知，采用合適的非矩形窗可以減小信號(hào)“頻譜泄露”，同時(shí)也可以增加譜峰的寬度，從而提高頻譜的分辨率。但是在對(duì)窗函數(shù)進(jìn)行選擇時(shí)，有如下的要求： (1) 窗口寬度M遠(yuǎn)小于樣本序列長度N,以排除不可靠的自相關(guān)值； (2)

18、 當(dāng)平穩(wěn)信號(hào)為實(shí)過程時(shí)，為保證平滑周期圖和真實(shí)功率譜同樣也是偶函數(shù)，平滑窗函數(shù)必須是實(shí)偶對(duì)稱函數(shù)； (3) 平滑窗函數(shù)應(yīng)在m=0處有峰值，并隨m絕對(duì)值增加而單調(diào)下降，使可靠的自相關(guān)值有較大的權(quán)值； (4) 功率譜是頻率的非負(fù)函數(shù)，由于周期圖是非負(fù)的，因而要求窗函數(shù)的傅里葉變換是非負(fù)的。在MATLAB中利用函數(shù)psd與pwelch都可實(shí)現(xiàn)Welch法的功率譜估計(jì)，其方法是一樣的，只是部分參數(shù)設(shè)置有所不同，psd函數(shù)的調(diào)用格式與Bartlett法功率譜估計(jì)相似，這里僅說明函數(shù)pwelch的調(diào)用格式。Pwelch函數(shù)的調(diào)用格式如下：Pxx,w=pwelch(x)Pxx,w=pwelch(x,wind

19、ow)Pxx,w=pwelch(x,window,noverlap)Pxx,w=pwelch(x,window,noverlap,nfft)Pxx,f=pwelch(x,window,noverlap,nfft,fs).=pwelch(x,window,noverlap,.range) Pwelch(.) 在使用pwelch函數(shù)計(jì)算序列功率譜時(shí)，還應(yīng)注意以下兩點(diǎn)： (1) 序列的分段數(shù)L為： （3-10） (2) 參數(shù)range用來指定頻率的取值范圍。若range=half，頻率取值范圍為0,Fs/2;若range=whole，取值范圍為0,Fs; 在MATLAB中，還可以利用函數(shù)csd實(shí)現(xiàn)W

20、elch法的互功率譜估計(jì)，該函數(shù)的調(diào)用格式如下： Pxy=csd(x,y) Pxy=csd(x,y,nfft) Pxy,f=csd(x,y,nfft,fs) Pxy=csd(x,y,nfft,fs,window) Pxy=csd(x,y,nfft,fs,window,numoverlap) Pxy=csd(x,y,.dflag) Pxy,Pxyc,f=csd(x,y,nfft,fs,window,numoverlap,p) csd(x,y,.) 圖 3-11序列的Welch法功率譜估計(jì)曲線3.5 正弦信號(hào)的提取將隨機(jī)信號(hào)與正弦信號(hào)進(jìn)行相加，得到混合信號(hào)，正弦信號(hào)與隨機(jī)信號(hào)如圖所示：圖3-12正

21、弦信號(hào)與隨機(jī)信號(hào)總 結(jié)本次設(shè)計(jì)的核心內(nèi)容是利用MATLAB，實(shí)現(xiàn)三種典型的數(shù)字信號(hào)，并對(duì)實(shí)現(xiàn)它們的均值、方差、相關(guān)及高斯隨機(jī)序列高斯隨機(jī)信號(hào)的統(tǒng)計(jì)特性的編程，實(shí)現(xiàn)連續(xù)時(shí)間周期信號(hào)頻域分析的仿真波形，運(yùn)用經(jīng)典法功率譜估計(jì)對(duì)正弦信號(hào)加噪進(jìn)行頻譜估計(jì)，并用濾波器進(jìn)行信號(hào)的提取。整個(gè)設(shè)計(jì)過程中對(duì)信號(hào)與系統(tǒng)與數(shù)字信號(hào)處理有了更深的理解，比如傅立葉變換、信號(hào)頻譜分析等；其次，實(shí)現(xiàn)過程是通過MATLAB軟件設(shè)計(jì)完成的，MATLAB 的圖形功能強(qiáng)大，具有良好的人機(jī)界面，此次設(shè)計(jì)過程中熟練了MATLAB的編程，掌握了很多函數(shù)的作用及使用方法；最后，通過此次課程設(shè)計(jì)，我對(duì)設(shè)計(jì)所用到的軟件MATLAB有了更加深刻地

22、了解，MATLAB不管在數(shù)值計(jì)算方面的功能很強(qiáng)大，而且其圖形仿真功能更能滿足各個(gè)領(lǐng)域的需要，因此我們以后更要經(jīng)常運(yùn)用MATLAB軟件，使其成為自己不可或缺的工具。 附 錄%正弦信號(hào)k1=-30;k2=30;k=k1:k2;w=pi/6;f=2*sin(k*0.5*w);figure(1);hold on;stem(k,f,'filled');title('正弦序列');xlabel('時(shí)間(k)');ylabel('幅值f(k)');fprintf('正弦信號(hào)的均值為%.4f 方差為%.4fn',mean(f),v

23、ar(f,1);N=64;k=0:N-1;w=pi/6;f=2*sin(k*0.5*w);f1=fft(f,N);figure(2);hold on;stem(k,f,'filled');xlabel('k');ylabel('幅值f(k)');grid on;title('時(shí)域信號(hào)');figure(3);hold on;subplot(2,1,1);stem(abs(f1),'filled');ylabel('幅值');title('幅頻特性');grid on;subplot(

24、2,1,2);stem(angle(f1),'filled');ylabel('相角');title('相頻特性');grid on;%周期信號(hào)統(tǒng)計(jì)分析k1=0;k2=3;k=k1:k2;Ts=1;f=k*Ts;xtilde=f'*ones(1,8);xtilde=xtilde(:);xtilde=xtilde'subplot(2,1,1);stem(k,f,'filled');title('一個(gè)周期');xlabel('時(shí)間(k)');ylabel('幅值f(k)'

25、);subplot(2,1,2);stem(xtilde,'filled');title('周期序列');xlabel('時(shí)間(k)');ylabel('幅值f(k)');fprintf('周期信號(hào)的均值為%.4f 方差為%.4fn',mean(xtilde),var(xtilde,1);k1=0;k2=3;k=k1:k2;Ts=1;f=k*Ts;xtilde=f'*ones(1,8);xtilde=xtilde(:);xtilde=xtilde'ff=fft(xtilde,32);figure(4

26、);hold on;stem(xtilde,'filled');title('周期序列');xlabel('時(shí)間(k)');ylabel('幅值f');grid on;figure(5);hold on;subplot(2,1,1);stem(abs(ff),'filled');ylabel('幅值');title('幅頻特性');grid on;subplot(2,1,2);stem(angle(ff),'filled');ylabel('相角');

27、title('相頻特性');grid on;.隨機(jī)信號(hào)及其分析.k=150;xn=randn(1,k);f_k=fft(xn,k);fprintf('隨機(jī)信號(hào)的均值為%.4f 方差為%.4fn',mean(xn),var(xn,1);figure(1);hold on;stem(xn,'filled');xlabel('時(shí)間（k）');ylabel('幅值xn(k)');title('高斯隨機(jī)信號(hào)');grid on;figure(2);hold on;subplot(2,1,1);stem(abs

28、(f_k),'filled');ylabel('幅值');title('幅頻特性');grid on;subplot(2,1,2);stem(angle(f_k),'filled');ylabel('相角');title('相頻特性');grid on;.離散信號(hào)的加減運(yùn)算.(1)正弦信號(hào)+/-隨機(jī)信號(hào)N=48;k=0:N-1;w=pi/10;f=sin(k*w);xn=randn(1,N);x=f+xn;y=f-xn;figure(1);hold on;subplot(2,2,1);stem(k,

29、f,'filled');xlabel('時(shí)間（k）');ylabel('幅值f(k)');title('f'); subplot(2,2,2);stem(k,xn,'filled');xlabel('時(shí)間（k）');ylabel('幅值xn(k)');title('xn');subplot(2,2,3);stem(k,x,'filled');xlabel('時(shí)間（k）');ylabel('幅值x(k)');title(&

30、#39;f+xn');subplot(2,2,4);stem(k,y,'filled');xlabel('時(shí)間（k）');ylabel('幅值y(k)');title('f-xn');(2)周期信號(hào)+/-隨機(jī)信號(hào)k=30;k1=1:4;Ts=1;f=k1*Ts;xtilde=f'*ones(1,6);xtilde=xtilde(:);xtilde=xtilde'xn=randn(1,k);x=f+xn;y=f-xn;figure(1);hold on;subplot(2,2,1);stem(xtilde,&

31、#39;filled');xlabel('時(shí)間（k）');ylabel('幅值f(k)');title('f'); subplot(2,2,2);stem(xn,'filled');xlabel('時(shí)間（k）');ylabel('幅值xn(k)');title('xn');subplot(2,2,3);stem(k,x,'filled');xlabel('時(shí)間（k）');ylabel('幅值x(k)');title('f

32、+xn');subplot(2,2,4);stem(k,y,'filled');xlabel('時(shí)間（k）');ylabel('幅值y(k)');title('f-xn');Fs=1000;%產(chǎn)生有噪信號(hào)n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*w1*n)+e;%繪制信號(hào)波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title(

33、9;x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%計(jì)算序列的DFTnfft=1024;Xk=fft(xn,nfft);%計(jì)算序列的PSDPxx1=abs(Xk).2/N;%繪制功率譜圖形index=0:round(nfft/2-1);k=index*N/nfft;plot_Pxx1=10*log10(Pxx1(index+1);figure(2)plot(k,plot_Pxx1)xlabel('頻率(Hz)')ylabel('功率（dB)')title('公式直接計(jì)算的功率譜')%periodogram函數(shù)計(jì)算的

34、功率譜window=boxcar(length(xn);Pxx2,f=periodogram(xn,window,nfft,Fs);plot_Pxx2=10*log10(Pxx2(index+1);figure(3)plot(k,plot_Pxx2)xlabel('頻率(Hz)')ylabel('功率(dB)')title('periodogram函數(shù)計(jì)算的功率譜')Fs=1000;%產(chǎn)生有噪信號(hào)n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*

35、w1*n)+e;%繪制信號(hào)波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title('x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%參數(shù)設(shè)計(jì)nfft=1024;window=boxcar(1001);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);pl

36、ot_Pxxc=10*log10(Pxxc(index+1);%功率譜對(duì)數(shù)曲線figure(2)plot(k,plot_Pxx)xlabel('頻率(Hz)')ylabel('功率（dB)')title('功率譜對(duì)數(shù)曲線')%置信區(qū)間中的功率譜曲線figure(3)plot(k,plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc)xlabel('頻率(Hz)')ylabel('功率(dB)')title('置信區(qū)間中的功率譜曲線')Fs=1000;%產(chǎn)生有噪

37、信號(hào)n=0:1/Fs:1;N=length(n);e=randn(1,N);w0=100*pi;w1=50*pi;xn=sin(2*w0*n)+sin(2*w1*n)+e;%繪制信號(hào)波形figure(1)plot(n,abs(xn)xlabel('n')ylabel('xn')title('x(n)=sin(2*wo*n)+sin(2*w1*n)+e(n)')%用Welch平均估計(jì)序列的功率譜nfft=1024;window=hamming(100);noverlap=20;range='half'Pxx,f=pwelch(xn,

38、window,noverlap,nfft,Fs,range);plot_Pxx=10*log10(Pxx);figure(2)plot(f,plot_Pxx)xlabel('頻率(Hz)')ylabel('功率（dB)')title('Welch法平均功率譜估計(jì)')fp=0.05;fs=3.14;rp=0.3;rs=64;wp=2*pi*fp;ws=2*pi*fs; %通帶、阻帶截止數(shù)字頻率n,wn=buttord(wp,ws,rp,rs,'s'); %確定巴特沃茲模擬濾波器的階次z,p,k=buttap(n); %設(shè)計(jì)歸一化巴特沃茲模擬低通濾波器bp,ap