信號(hào)分析與處理實(shí)驗(yàn)：ch1_6 實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析

1、實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 一、實(shí)驗(yàn)?zāi)康囊?、?shí)驗(yàn)?zāi)康?深刻理解隨機(jī)信號(hào)的特性，掌握隨機(jī)信號(hào)功率譜估計(jì)的基本原理，靈活運(yùn)用各種隨機(jī)信號(hào)功率譜估計(jì)的基本方法。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 功率譜估計(jì)是隨機(jī)信號(hào)處理中的一個(gè)重要的研究和應(yīng)用領(lǐng)域。功率譜估計(jì)基本上可以分為非參數(shù)估計(jì)的經(jīng)典方法和參數(shù)估計(jì)的近代方法。 經(jīng)典功率譜估計(jì)是基于FFT算法的非參數(shù)估計(jì)，對(duì)足夠長的記錄數(shù)據(jù)效果較好。近代功率譜估計(jì)是建立在隨機(jī)信號(hào)參數(shù)模型的基礎(chǔ)上，通過信號(hào)參數(shù)模型或預(yù)測(cè)誤差濾波器參數(shù)的估計(jì)，實(shí)現(xiàn)功率譜估計(jì)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功

2、率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 經(jīng)典功率譜估計(jì)法獲得廣泛應(yīng)用的是修正周期圖法。該方法將數(shù)據(jù)分段和加窗處理后，先求各段數(shù)據(jù)的功率譜，然后再對(duì)各段功率譜進(jìn)行平均。修正周期圖法有效地減少了方差和偏差，提高了估計(jì)質(zhì)量，得到普遍應(yīng)用。經(jīng)典功率譜估計(jì)法存在以下主要缺陷：1. 利用有限的N個(gè)觀察數(shù)據(jù)進(jìn)行自相關(guān)估計(jì)，隱含著在已知N個(gè)數(shù)據(jù)之外的全部數(shù)據(jù)均為零的假設(shè)；2. 假定數(shù)據(jù)是由N個(gè)觀察數(shù)據(jù)以N為周期的周期性延拓。由于在計(jì)算過程中采用了加窗處理，因而功率泄漏較大，頻率分辨率較低，不適用于短序列的譜分析和微弱信號(hào)的檢測(cè)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 近代功率

3、譜估計(jì)既不需要加窗，又不需要對(duì)相關(guān)函數(shù)的估計(jì)進(jìn)行如經(jīng)典法那樣的假設(shè)，從而減少了功率泄漏，提高了頻譜分辨率。常用的參數(shù)模型有自回歸(AR)模型、滑動(dòng)平均(MA)模型、自回歸滑動(dòng)平均(ARMA)模型。其中AR模型是基本模型，求解AR模型的參數(shù)主要有L-D算法和Burg算法。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理已知某隨機(jī)信號(hào)由兩個(gè)余弦信號(hào)和噪聲構(gòu)成：s(t)是均值為0且方差為1的白噪聲。試用周期圖法和修正周期圖法（Welch法）估計(jì)信號(hào)的功率譜。% 產(chǎn)生并顯示隨機(jī)信號(hào)x(t) Fs = 1000; % sampling frequency t = 0:1/Fs

4、:1; % one second of samples x = sin(2*pi*50*t) + 2*sin(2*pi*120*t) + randn(size(t); plot(t,x);)()220cos()210cos()(tstttx00.20.40.60.81-6-4-20246Time (s)x(t)waveform of signal實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理% 直接利用周期圖法估計(jì)信號(hào)的功率譜。% 選取DFT的點(diǎn)數(shù)N=1024 Pxx = abs(fft(x,1024).2/1001; subplot(2,1,1); plot(0

5、:1023*Fs/1024,10*log10(Pxx); axis(0,1023,-30,30);01002003004005006007008009001000-30-20-100102030power spectrum (dB) periodgram01002003004005006007008009001000-30-20-100102030Frequency (Hz)power spectrum (dB)256 sample section, nonovelapping實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理直接周期圖法的不足是方差大且方差不隨采樣點(diǎn)

6、數(shù)增加而減小。減小方差的方法是把數(shù)據(jù)x分段，每段做周期圖，然后相加取平均值。分段時(shí)不重疊的平均周期圖法稱為Bartlett法。 % Bartlett法估計(jì)隨機(jī)信號(hào)的功率譜 Pxx1 = (abs(fft(x (1:256).2+abs(fft(x(257:512).2+ abs(fft(x(513:768).2)/(256*3); subplot(2,1,2); plot(0:255*Fs/256,10*log10(Pxx1); axis(0,1023,-30,30);01002003004005006007008009001000-30-20-100102030power spectrum

7、(dB) periodgram01002003004005006007008009001000-30-20-100102030Frequency (Hz)power spectrum (dB)256 sample section, nonovelapping實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理Bartlett法雖然可以通過分段而減少方差，但其分段數(shù)不能太多，因?yàn)榉侄螖?shù)增加會(huì)使得每段數(shù)據(jù)點(diǎn)數(shù)減少，造成頻率分辨率下降。為了減小偏差，提高頻率分辨率，Welch對(duì)Bartlett法進(jìn)行了改進(jìn)，提出了重疊平均周期圖法(Welch法)。% Welch 法估計(jì)隨機(jī)信號(hào)

8、的功率譜，選取128點(diǎn)重疊分段 Pxx2 = (abs(fft(x( 1:256).2 + abs(fft(x(129:384).2 + abs(fft(x(257:512).2 + abs(fft(x(385:640).2 + abs(fft(x(513:768).2 +abs(fft(x(641:896).2 ) / (256*6); subplot(2,1,1); plot(0:255*Fs/256,10*log10(Pxx2); axis(0,1023,-30,30);01002003004005006007008009001000-30-20-100102030power spect

9、rum (dB)256 sample section,128 sample overlap01002003004005006007008009001000-30-20-100102030Frequency (Hz)power spectrum (dB)Averaged modified periodogram,256 sample section,128sample overlap,hanning window實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理為了進(jìn)一步提高Welch法估計(jì)質(zhì)量，可對(duì)每段數(shù)據(jù)加非矩形窗，然后再求各段的功率譜，即“修正周期圖法”。非矩形窗

10、由于通過降低旁瓣的幅度而減小旁瓣造成的“譜泄露”。%修正周期圖法估計(jì)功率譜，選取256點(diǎn)哈寧窗重新計(jì)算 w = hanning(256); Pxx3 = (abs(fft(w.*x( 1:256).2 + abs(fft(w.*x(129:384).2 + abs(fft(w.*x(257:512).2 + abs(fft(w.*x(385:640).2 + abs(fft(w.*x(513:768).2 + abs(fft(w.*x(641:896).2 ) ； Pxx3=Pxx3/ (norm(w)2*6); subplot(2,1,2); plot(0:255*Fs/256,10*log

11、10(Pxx3); axis(0,1023,-30,30);01002003004005006007008009001000-30-20-100102030power spectrum (dB)256 sample section,128 sample overlap01002003004005006007008009001000-30-20-100102030Frequency (Hz)power spectrum (dB)Averaged modified periodogram,256 sample section,128sample overlap,hanning window實(shí)驗(yàn)六實(shí)

12、驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 MATLAB提供了經(jīng)典功率譜估計(jì)的函數(shù)。函數(shù)periodogram可以實(shí)現(xiàn)周期圖法的功率譜估計(jì)，其調(diào)用格式為 Pxx,F = PERIODOGRAM(x, WINDOW,NFFT,Fs)其中： x為進(jìn)行功率譜估計(jì)的輸入有限長序列； WINDOW用于指定采用的窗函數(shù)，默認(rèn)值為矩形窗； NFFT為DFT的點(diǎn)數(shù)，一般取大于輸入序列x的長度，默認(rèn)值為256； Fs 是繪制功率譜曲線的抽樣頻率，默認(rèn)值為1； Pxx為功率譜估計(jì)值； F為Pxx值所對(duì)應(yīng)的頻率點(diǎn)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理

13、 MATLAB提供了經(jīng)典功率譜估計(jì)的函數(shù)。Welch-Bartlett平均周期圖法可以利用函數(shù)psd實(shí)現(xiàn)，其調(diào)用格式為 Pxx,F = PSD(x, NFFT, Fs, WINDOW, NOVERLAP)其中： 參數(shù)x,NFFT,Fs用法同periodogram 函數(shù)； WINDOW用于指定采用的窗函數(shù)，默認(rèn)值為hanning窗； NOVERLAP指定分段重疊的樣本數(shù)。 如果使用boxcar窗，且NOVERLAP=0，則為Bartlett法。 如果NOVERLAP=length(x)/2，則為重疊50%的Welch法。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理

14、 MATLAB提供了經(jīng)典功率譜估計(jì)的函數(shù)。Welch法還可以利用函數(shù)pwelch實(shí)現(xiàn)，其調(diào)用格式為Pxx,F = PWELCH(x, WINDOW, NOVERLAP, NFFT, Fs)函數(shù)pwelch和函數(shù)psd進(jìn)行功率譜估計(jì)的方法相同，只是參數(shù)設(shè)置略有不同。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 MATLAB也提供了近代功率譜估計(jì)的函數(shù)。函數(shù)levinson可以實(shí)現(xiàn)L-D遞推算法，其求解AR模型參數(shù)及白噪聲序列的方差，其調(diào)用格式為A = levinson (R, ORDER)其中：ORDER為AR模型的階數(shù)； R為觀測(cè)序列的自相關(guān)函數(shù)； 返回值A(chǔ)即為

15、白噪聲序列的方差和AR模型參數(shù)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 MATLAB也提供了近代功率譜估計(jì)的函數(shù)。函數(shù)aryule也可以實(shí)現(xiàn)L-D遞推算法，它與函數(shù)levinson的區(qū)別僅在于輸入?yún)?shù)不同，其調(diào)用格式為 A,E,K= aryule (x, ORDER)其中： ORDER為AR模型的階數(shù)； x為觀測(cè)序列； E為預(yù)測(cè)誤差；K為反射系數(shù)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 二、實(shí)驗(yàn)原理二、實(shí)驗(yàn)原理 MATLAB也提供了近代功率譜估計(jì)的函數(shù)。Burg算法可以利用函數(shù)arburg實(shí)現(xiàn)，其調(diào)用格式為 A,E,K= ARBURG(x, O

16、RDER)其中： ORDER為AR模型的階數(shù)； x為觀測(cè)序列； E為預(yù)測(cè)誤差；K為反射系數(shù)。實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 三、實(shí)驗(yàn)內(nèi)容三、實(shí)驗(yàn)內(nèi)容1. 某隨機(jī)信號(hào)由兩余弦信號(hào)及噪聲構(gòu)成：其中，s(t)是均值為0且方差為1的高斯白噪聲。 (1) 生成信號(hào)x(t) ，繪出其時(shí)域波形； (2) 試分別用周期圖法、Bartlett法和Welch法估計(jì)該序列的功率譜， 并對(duì)上述各種估計(jì)的效果進(jìn)行比較； (3) 選用不同窗函數(shù)，利用修正周期圖法估計(jì)其功率譜，有何結(jié)論？2. 使用MATLAB函數(shù)load mtlb 調(diào)用MATLAB中的一段語音信號(hào)，利用AR模型估計(jì)其功率譜，并討論不同階

17、數(shù)的AR模型對(duì)隨機(jī)信號(hào)功率譜估計(jì)的影響。)()40cos()20cos()(tstttx實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 三、實(shí)驗(yàn)內(nèi)容三、實(shí)驗(yàn)內(nèi)容3.在某個(gè)音樂廳中，原始音頻信號(hào)由于墻壁和天花板的反射而產(chǎn)生回音。觀眾感受到的音頻信號(hào)是xk與其回音的合成，即若 ，a=0.1, n=50， 產(chǎn)生200個(gè)yk的樣本值，求出其自相關(guān)函數(shù)Ryn， 根據(jù)Ryn能否得到a和n？4. 若對(duì)某地區(qū)每隔6小時(shí)采集其溫度達(dá)一年，得到序列xk，假設(shè)該序列為平穩(wěn)序列，試通過該序列預(yù)測(cè)此地區(qū)下一年的溫度。nkaxkxkykkkx6 . 0cos2 . 0cos實(shí)驗(yàn)六實(shí)驗(yàn)六 隨機(jī)信號(hào)功率譜分析隨機(jī)信號(hào)功率譜分析 三、實(shí)驗(yàn)內(nèi)容三、實(shí)驗(yàn)內(nèi)容5. 已知某含有噪聲的信號(hào)其中h (t)是均值為0且方差為1的高斯白噪聲。 (1) 生成信號(hào) x(t)，繪出其時(shí)域