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

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