現(xiàn)代信號(hào)處理經(jīng)典的功率譜估計(jì)

..《現(xiàn)代信號(hào)處理》__李建強(qiáng)__201512172087專業(yè)：電子科學(xué)與技術(shù)作業(yè)內(nèi)容：在MATLAB平臺(tái)上對(duì)一個(gè)特定的平穩(wěn)隨機(jī)信號(hào)進(jìn)行經(jīng)典功率譜估計(jì)和現(xiàn)代功率譜估計(jì)的比較前言功率譜估計(jì)是信息學(xué)科中的研究熱點(diǎn),在過去的30多年里取得了飛速的發(fā)展。在許多工程應(yīng)用中,它能給出被分析對(duì)象的能量隨頻率的分布情況。平滑周期圖是一種計(jì)算簡(jiǎn)單的經(jīng)典方法,它的主要特點(diǎn)是與任何模型參數(shù)無(wú)關(guān),但估計(jì)出來的功率譜很難與信號(hào)的真是功率譜相匹配。與周期圖方法不同,現(xiàn)代譜估計(jì)主要是針對(duì)經(jīng)典譜估計(jì)<周期圖和自相關(guān)法>的分辨率低和方差性能不好的問題而提出的。其使用參數(shù)化的模型,能夠給出比周期圖方法高得多的頻率分辨率。其內(nèi)容極其豐富,涉及的學(xué)科和領(lǐng)域也相當(dāng)廣泛,按是否有參數(shù)大致可分為參數(shù)模型估計(jì)和非參數(shù)模型估計(jì),前者有AR模型、MA模型、ARMA模型、PRONY指數(shù)模型等;后者有最小方差方法、多分量的MUSIC方法等?？傮w概述本次實(shí)驗(yàn)分別使用經(jīng)典的功率譜估計(jì)〔如周期圖法與AR模型法對(duì)某一特定的平穩(wěn)隨機(jī)信號(hào)進(jìn)行其功率譜估計(jì),由圖像得到信號(hào)的頻率。利用MATLAB平臺(tái),直觀形象地觀察并比較二者估計(jì)效果的區(qū)別,以便于加深對(duì)功率譜估計(jì)的理解和掌握。具體的實(shí)現(xiàn)步驟經(jīng)典法功率譜估計(jì)周期圖法又稱直接法,它是從隨機(jī)信號(hào)x<n>中截取N長(zhǎng)的一段,把它視為能量有限的真實(shí)功率譜的估計(jì)的一個(gè)抽樣。實(shí)現(xiàn)步驟〔1、模擬系統(tǒng)輸出參數(shù)x<n>=A*sin<2πf1*n>+B*sin<2πf2*n>,包括序列長(zhǎng)度N〔128或512或1024,加性高斯白噪聲〔AGWN功率一定,設(shè)置A,B,f1,f2,n的值?！?、應(yīng)用周期圖法〔不加窗對(duì)信號(hào)的功率譜密度進(jìn)行估計(jì),使用直接法在MATLAB平臺(tái)上進(jìn)行編程實(shí)現(xiàn)?！?、輸出相應(yīng)波形圖,進(jìn)行觀察,記錄。1.2MATLAB源代碼實(shí)現(xiàn)clearall;%清除工作空間所有之前的變量closeall; %關(guān)閉之前的所有的figureclc;%清除命令行之前所有的文字n=1:1:128;%設(shè)定采樣點(diǎn)n=1-128f1=0.2;%設(shè)定f1頻率的值0.2f2=0.213;%設(shè)定f2頻率的值0.213A=1;%取定第一個(gè)正弦函數(shù)的振幅B=1;%取定第一個(gè)正弦函數(shù)的振幅a=0; %設(shè)定相位為0x1=A*sin<2*pi*f1*n+a>+B*sin<2*pi*f2*n+a>; %定義x1函數(shù),不添加高斯白噪聲x2=awgn<x1,3>;%在x1基礎(chǔ)上添加加性高斯白噪聲,信噪比為3,定義x2函數(shù)temp=0;%定義臨時(shí)值,并規(guī)定初始值為0temp=fft<x2,128>;%對(duì)x2做快速傅里葉變換pw1=abs<temp>.^2/128;%對(duì)temp做經(jīng)典功率估計(jì)k=0:length<temp>-1;w=2*pi*k/128;figure<1>;%輸出x1函數(shù)圖像plot<w/pi/2,pw1>%輸出功率譜函數(shù)pw1圖像xlabel<'信號(hào)頻率/Hz'>;ylabel<'PSD/傅立葉功率譜估計(jì)'>;title<'正弦信號(hào)x<n>添加高斯白噪聲后的,周期圖法功率頻譜分析'>;grid;%-------------------------------------------------------------------------pw2=temp.*conj<temp>/128;%對(duì)temp做向量的共軛乘積k=0:length<temp>-1;w=2*pi*k/128;figure<2>;plot<w/pi/2,pw2>;%輸出功率譜函數(shù)pw2圖像xlabel<'信號(hào)頻率/Hz'>;ylabel<'PSD/傅立葉功率譜估計(jì)'>;title<'正弦信號(hào)x〔n自相關(guān)法功率譜估計(jì)'>;grid;1.3matlab仿真圖形〔1、用直接法,功率譜圖像,采樣點(diǎn)N=128。 〔2用直接法,功率譜圖像,采樣點(diǎn)N=512。1.4、經(jīng)典功率譜估計(jì)分析當(dāng)采樣的點(diǎn)數(shù)為N=128時(shí),此時(shí)采樣的得到的圖像分辨力很低,并且分辨率也比較低,這就導(dǎo)致了功率譜圖像只能看到一個(gè)峰值點(diǎn)。采樣點(diǎn)數(shù)為N=512時(shí),此時(shí),分辨力和分辨率比較高,可以清楚的區(qū)分到兩個(gè)峰值點(diǎn)的橫坐標(biāo),此時(shí)的橫坐標(biāo)就是信號(hào)的頻率。但是這是以犧牲效率為代價(jià)的,采樣的點(diǎn)數(shù)越多,所花的時(shí)間越長(zhǎng),這在實(shí)際的工程中是不切合實(shí)際的,因此,在我們估計(jì)隨機(jī)信號(hào)的頻率的時(shí)候,要合理的采取樣本點(diǎn)數(shù),盡可能的采取多的樣點(diǎn),來接近真實(shí)的信號(hào)頻率,也要考慮實(shí)際的效率問題。2、AR模型一般最小二乘法譜分析方法要求ARMA模型的階數(shù)和參數(shù)以及噪聲的方差已知.然而這類要求在實(shí)際中是不可能提供的,即除了一組樣本值x<1>,x<2>,…,x<T>以供利用<有時(shí)會(huì)有一定的先驗(yàn)知識(shí)>外,再?zèng)]有其它可用的數(shù)據(jù).因此必須估計(jì)有關(guān)的階數(shù)和參數(shù),以便獲得譜密度的估計(jì)。2.1實(shí)現(xiàn)步驟〔1、模擬系統(tǒng)輸出參數(shù)y=A*sin<2πf1*n>+B*sin<2πf2*n>,包括序列長(zhǎng)度N,加性高斯白噪聲〔AGWN,設(shè)置A,B,f1,f2,n的值?！?、應(yīng)用AR模型一般最小二乘法對(duì)信號(hào)進(jìn)行功率譜估計(jì),編寫程序。取定|B<z>|=1,構(gòu)造AR模型,然后不斷變換p的值,觀察不同p值下功率譜密度波形的分辨率高低?！?、輸出相應(yīng)波形圖,進(jìn)行觀察,記錄。2.2源代碼%AR模型的一般最小二乘估計(jì)%-----------------------------------------clearall;%清除workspace之前的變量closeall;%關(guān)閉之前的圖像clc;%清除命令行之前的文字n=[1:128];%取定采樣點(diǎn)n=1至128f1=0.2;%取定f1頻率的值f2=0.213;%取定f2頻率的值<根據(jù)f1與f2之差=2*pi/n=0.0491A=sqrt<20>;%取定第一個(gè)正弦函數(shù)的振幅B=sqrt<2>;%取定第一個(gè)正弦函數(shù)的振幅x=A*sin<2*pi*f1*n>+B*sin<2*pi*f2*n>;%定義x函數(shù)noise=0+1*randn<1,length<n>>;%添加均值為0、方差為1的高斯白噪聲xn=x+noise;%在x1基礎(chǔ)上添加加性高斯白噪聲,定義xn函數(shù)m=xcorr<xn>;%m為xn的自相關(guān)函數(shù)〔序列>%-----------------------------------------p=100;%取定R的階數(shù),更改p的值,觀察相對(duì)應(yīng)的譜估計(jì)q=125;%此處一定要滿足q>=pfori=1:pforj=1:pR<i,j>=m<q+i+j-1-p>;%構(gòu)造一個(gè)p*p階的自相關(guān)矩陣〔Hankel矩陣%<課本P883.4.33aendendRlegnth=size<R>%輸出驗(yàn)證R矩陣的行列數(shù)的值fori=1:p%i=1~pr<i>=m<q+i>;%定義一個(gè)1*p的向量,對(duì)應(yīng)課本P883.4.22endr=-r';%對(duì)應(yīng)課本P88Ra=-ra=<inv<R'*R>*R'>*r;%用LS方法求解aa1=fliplr<a>%對(duì)應(yīng)課本P88b,將a進(jìn)行元素對(duì)調(diào),使a1=[ap,...,a1]'figure<1>;freqz<1,a1,128,1>;title<'AR模型的一般最小二乘估計(jì)'>;legend<strcat<'AR階數(shù)=',int2str<p>>>;gridon;2.3matlab仿真圖形〔1、當(dāng)p=4時(shí),信號(hào)的功率譜密度波形：〔2、當(dāng)p=64時(shí),信號(hào)的功率譜密度波形：〔3、當(dāng)p=100時(shí),信號(hào)的功率譜密度波形：3、AR模型的總體最小二乘法一般最小二乘法會(huì)帶來兩個(gè)問題：其一,必須重新列出方程組,使它只包含p個(gè)未知數(shù)；其二,求解Ax=b的最小二乘方法只認(rèn)為b含有誤差,但實(shí)際上系數(shù)矩陣A也含有誤差。因此,引入總體最小二乘法〔SVD-TLS可以比一般最小二乘法更合理地同時(shí)考慮A和b的誤差或擾動(dòng)。3.1實(shí)現(xiàn)步驟〔1、模擬系統(tǒng)輸出參數(shù)y=Asin<2πf1*n>+Bsin<2πf2*n>,包括序列長(zhǎng)度N,加性高斯白噪聲〔AGWN,設(shè)置A,B,f1,f2,n的值。〔2、應(yīng)用AR模型的總體最小二乘法對(duì)信號(hào)進(jìn)行功率譜估計(jì),按照課本P96求總體最小二乘解的算法步驟,編寫程序。取定|B<z>|=1,構(gòu)造AR模型,然后不斷變換p的值,觀察不同p值下功率譜密度波形的分辨率高低?！?、輸出相應(yīng)波形圖,進(jìn)行觀察,記錄?！?、算法步驟〔SVD-TLS算法步驟1：計(jì)算增廣矩陣B的SVD,并存儲(chǔ)奇異值和矩陣V；步驟2：確定增廣矩陣B的有效秩p；步驟3：利用式〔3.4.56和式〔3.4.53計(jì)算矩陣S〔p；步驟4：求S〔p的逆矩陣S-<p>,并由式〔3.4.59計(jì)算未知參數(shù)的總體最小二乘估計(jì)。3.2源代碼%實(shí)現(xiàn)AR模型的總體最小二乘估計(jì)〔SVD-TLS算法%-----------------------------------------clearall;%清除workspace之前的變量closeall;%關(guān)閉之前的圖像clc;%清除命令行之前的文字n=[1:128];%取定采樣點(diǎn)n=1至128f1=0.2;%取定f1頻率的值f2=0.213;%取定f2頻率的值A(chǔ)=sqrt<20>;%取定第一個(gè)正弦函數(shù)的振幅B=sqrt<2>;%取定第二個(gè)正弦函數(shù)的振幅x=A*sin<2*pi*f1*n>+B*sin<2*pi*f2*n>;%定義x函數(shù)noise=0+1*randn<1,length<n>>;%添加均值為0、方差為1的高斯白噪聲xn=x+noise;%在x1基礎(chǔ)上添加加性高斯白噪聲,定義xn函數(shù)m=xcorr<xn>;%m為xn的自相關(guān)函數(shù)〔序列>%-----------------------------------------p=100;%取定R的階數(shù),更改p=4,64,100,的值,觀察%相對(duì)應(yīng)的譜估計(jì)q=125;%此處一定要滿足q>=pfori=1:pforj=1:pR<i,j>=m<q+i+j-1-p>;%構(gòu)造一個(gè)pxp階的自相關(guān)矩陣〔Hankel矩陣%<課本P883.4.33aendendRlegnth=size<R>%輸出驗(yàn)證R矩陣的行列數(shù)的值fori=1:p%i=1~pr<i>=m<q+i>;%定義一個(gè)1*p的向量,對(duì)應(yīng)課本P883.4.22endB=[-r',R];%對(duì)應(yīng)P94b中的B[U,K,V]=svd<B>;%由P96算法步驟1求得增廣矩陣B的%SVD,并存儲(chǔ)奇異值和矩陣VP=rank<B>;%由P96算法步驟2求得增廣矩陣B的有%效秩,定義為PS=zeros<P+1>;%構(gòu)造一個(gè)<p+1>*<p+1>維的矩陣S,對(duì)應(yīng)課本P95forj=1:pfori=1:p+1-Pdjj=K<j,j>*K<j,j>;%對(duì)應(yīng)課本P96,構(gòu)造djj,并求其平方vij=V<i:i+p,j>;%對(duì)應(yīng)課本P96和課本P953.4.53,構(gòu)造vijS=S+djj*vij*vij';%對(duì)應(yīng)課本P96,計(jì)算矩陣S的二重級(jí)數(shù)求和endendSni=inv<S>;%對(duì)應(yīng)課本P96算法步驟4,求S逆矩陣a=zeros<1,P>;%對(duì)應(yīng)課本P88b,構(gòu)造a矩陣fori=1:Pa<1,i>=Sni<i+1,1>/Sni<1,1>;%對(duì)應(yīng)課本P96,求出矩陣a=[a1,...,ap]'enda1=fliplr<a>%對(duì)應(yīng)課本P88b,將a進(jìn)行元素對(duì)調(diào),使a1=[ap,...,a1]'figure<1>;freqz<1,a1,128,1>;%求出信號(hào)的幅頻響應(yīng)