功率譜密度估計(jì)方法的MATLAB實(shí)現(xiàn)(共12頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上功率譜密度估計(jì)方法的MATLAB實(shí)現(xiàn)在應(yīng)用數(shù)學(xué)和物理學(xué)中，譜密度、功率譜密度和能量譜密度是一個(gè)用于信號(hào)的通用概念，它表示每赫茲的功率、每赫茲的能量這樣的物理量綱。在物理學(xué)中，信號(hào)通常是波的形式，例如電磁波、隨機(jī)振動(dòng)或者聲波。當(dāng)波的頻譜密度乘以一個(gè)適當(dāng)?shù)南禂?shù)后將得到每單位頻率波攜帶的功率，這被稱為信號(hào)的功率譜密度（power spectral density, PSD）或者譜功率分布（spectral power distribution, SPD）。功率譜密度的單位通常用每赫茲的瓦特?cái)?shù)（W/Hz）表示，或者使用波長(zhǎng)而不是頻率，即每納米的瓦特?cái)?shù)（W/nm）來(lái)表示。信號(hào)的

2、功率譜密度當(dāng)且僅當(dāng)信號(hào)是廣義的平穩(wěn)過(guò)程的時(shí)候才存在。如果信號(hào)不是平穩(wěn)過(guò)程，那么自相關(guān)函數(shù)一定是兩個(gè)變量的函數(shù)，這樣就不存在功率譜密度，但是可以使用類(lèi)似的技術(shù)估計(jì)時(shí)變譜密度。信號(hào)功率譜的概念和應(yīng)用是電子工程的基礎(chǔ)，尤其是在電子通信系統(tǒng)中，例如無(wú)線電和微波通信、雷達(dá)以及相關(guān)系統(tǒng)。因此學(xué)習(xí)如何進(jìn)行功率譜密度估計(jì)十分重要，借助于Matlab工具可以實(shí)現(xiàn)各種譜估計(jì)方法的模擬仿真并輸出結(jié)果。下面對(duì)周期圖法、修正周期圖法、最大熵法、Levinson遞推法和Burg法的功率譜密度估計(jì)方法進(jìn)行程序設(shè)計(jì)及仿真并給出仿真結(jié)果。以下程序運(yùn)行平臺(tái)：Matlab R2015a（8.5.0.）一、 周期圖法譜估計(jì)程序1、

3、源程序Fs=; %采樣頻率100kHzN=1024; %數(shù)據(jù)長(zhǎng)度N=1024n=0:N-1;t=n/Fs;xn=sin(2000*2*pi*t); %正弦波,f=2000HzY=awgn(xn,10); %加入信噪比為10db的高斯白噪聲subplot(2,1,1);plot(n,Y) title('信號(hào)')xlabel('時(shí)間');ylabel('幅度');grid on;window=boxcar(length(xn); %矩形窗nfft=N/4; %采樣點(diǎn)數(shù)Pxx f=periodogram(Y,window,nfft,Fs); %直接法s

4、ubplot(2,1,2);plot(f,10*log10(Pxx);grid on;title('周期圖法譜估計(jì),',int2str(N),'點(diǎn)');xlabel('頻率（Hz）');ylabel('功率譜密度');2、 仿真結(jié)果二、 修正周期圖法（加窗）譜估計(jì)程序1、源程序Fs=; %采樣頻率100kHzN=512; %數(shù)據(jù)長(zhǎng)度M=32; %漢明窗寬度n=0:N-1;t=n/Fs;xn=sin(2000*2*pi*t); %正弦波,f=2000HzY=awgn(xn,10); %加入信噪比為10db的高斯白噪聲subplot(

5、2,1,1);subplot(2,1,1);plot(n,Y) title('信號(hào)')xlabel('時(shí)間');ylabel('幅度');grid on;window=hamming(M); %漢明窗Pxx f=pwelch(Y,window,10,256,Fs); subplot(2,1,2);plot(f,10*log10(Pxx);grid on;title('修正周期圖法譜估計(jì) N=',int2str(N),' M=',int2str(M);xlabel('頻率（Hz）');ylabel(&

6、#39;功率譜密度');2、 仿真結(jié)果三、 最大熵法譜估計(jì)程序1、源程序fs=1; %設(shè)采樣頻率N=128; %數(shù)據(jù)長(zhǎng)度 改變數(shù)據(jù)長(zhǎng)度會(huì)導(dǎo)致分辨率的變化；f1=0.2*fs; %第一個(gè)sin信號(hào)的頻率,f1/fs=0.2f2=0.3*fs; %第二個(gè)sin信號(hào)的頻率，f2/fs=0.2或者0.3P=10; %濾波器階數(shù) n=1:N; s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs); %s為原始信號(hào)x=awgn(s,10); %x為觀測(cè)信號(hào)，即對(duì)原始信號(hào)加入白噪聲，信噪比10dBfigure(1); %畫(huà)出原始信號(hào)和觀測(cè)信號(hào)subplot(2,1,1);plo

7、t(s,'b'),xlabel('時(shí)間'),ylabel('幅度'),title('原始信號(hào)s');grid;subplot(2,1,2);plot(x,'r'),xlabel('時(shí)間'),ylabel('幅度'),title('觀測(cè)信號(hào)x');Pxx1,f=pmem(x,P,N,fs); %最大熵譜估計(jì)figure(2);plot(f,10*log10(Pxx1);xlabel('頻率(Hz) ');ylabel('功率譜(dB) '

8、);title('最大熵法譜估計(jì) 模型階數(shù)P=',int2str(P),' 數(shù)據(jù)長(zhǎng)度N=',int2str(N);2、 仿真結(jié)果四、 Levinson遞推法譜估計(jì)程序1、 源程序fs=1; %設(shè)采樣頻率為1N=1000; %數(shù)據(jù)長(zhǎng)度 改變數(shù)據(jù)長(zhǎng)度會(huì)導(dǎo)致分辨率的變化；f1=0.2*fs; %第一個(gè)sin信號(hào)的頻率,f1/fs=0.2f2=0.3*fs; %第二個(gè)sin信號(hào)的頻率，f1/fs=0.2或者0.3M=16; %濾波器階數(shù)的最大取值,超過(guò)則認(rèn)為代價(jià)太大而放棄L=2*N; %有限長(zhǎng)序列進(jìn)行離散傅里葉變換前，序列補(bǔ)零的長(zhǎng)度n=1:N; s=sin(2*pi*f

9、1*n/fs)+sin(2*pi*f2*n/fs);%s為原始信號(hào)x=awgn(s,10);%x為觀測(cè)信號(hào)，即對(duì)原始信號(hào)加入白噪聲，信噪比10dBfigure(1); %畫(huà)出原始信號(hào)和觀測(cè)信號(hào)subplot(2,1,1);plot(s,'b'),axis(0 100 -3 3),xlabel('時(shí)間'),ylabel('幅度'),title('原始信號(hào)s');grid;subplot(2,1,2);plot(x,'r'),axis(0 100 -3 3),xlabel('時(shí)間'),ylabel(&#

10、39;幅度'),title('觀測(cè)信號(hào)x');grid;%計(jì)算自相關(guān)函數(shù)rxx = xcorr(x,x,M,'biased');%計(jì)算有偏估計(jì)自相關(guān)函數(shù)，長(zhǎng)度為-M到M，%共2M+1r0 = rxx(M+1); %r0為零點(diǎn)上的自相關(guān)函數(shù)，相對(duì)于-M，第M+1個(gè)點(diǎn)為零點(diǎn)R = rxx(M+2:2*M+1);% R為從1到第M個(gè)點(diǎn)的自相關(guān)函數(shù)矩陣%確定矩陣大小a = zeros(M,M);FPE = zeros(1,M);%FPE：最終預(yù)測(cè)誤差，用來(lái)估計(jì)模型的階次var = zeros(1,M);%求初值a(1,1) = -R(1)/r0;%一階模型參數(shù)v

11、ar(1) = (1-(abs(a(1,1)2)*r0;%一階方差FPE(1) = var(1)*(M+2)/(M);%遞推for p=2:M sum=0; for k=1:p-1%求a(p,p) sum=sum+a(p-1,k)*R(p-k); end a(p,p)=-(R(p)+sum)/var(p-1); for k=1:p-1 %求a(p,k) a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); end var(p)=(1-a(p,p)2)*var(p-1); %求方差 FPE(p)=var(p)*(M+1+p)/(M+1-p);%求最終預(yù)測(cè)誤差end %確定AR模型

12、的最佳階數(shù)min=FPE(1); %求出FPE最小時(shí)對(duì)應(yīng)的階數(shù)p = 1;for k=2:M if FPE(k)<min min=FPE(k); p=k; endend %功率譜估計(jì)W=0.01:0.01:pi; %功率譜以2*pi為周期，又信號(hào)為實(shí)信號(hào)，只需輸出0到PI即可；he=ones(1,length(W); %length()求向量的長(zhǎng)度f(wàn)or k=1:p he=he+(a(p,k).*exp(-j*k*W); endPxx=var(p)./(abs(he).2); %功率譜函數(shù)；F=W*fs/(pi*2); %將角頻率坐標(biāo)換算成HZ坐標(biāo)，便于觀察；重要！figure;plot

13、(F,abs(Pxx) xlabel('頻率/Hz'),ylabel('功率譜P'),title(' AR模型的最佳階數(shù)p=' int2str(p) ); grid;2、 仿真結(jié)果五、 Burg法譜估計(jì)程序1、 源程序fs=1;%設(shè)采樣頻率為1N=900;%數(shù)據(jù)長(zhǎng)度 改變數(shù)據(jù)長(zhǎng)度會(huì)導(dǎo)致分辨率的變化；f1=0.2*fs;%第一個(gè)sin信號(hào)的頻率,f1/fs=0.2f2=0.3*fs;%第二個(gè)sin信號(hào)的頻率，f1/fs=0.2或者0.3M=512;%濾波器階數(shù)的最大取值,超過(guò)則認(rèn)為代價(jià)太大而放棄n=1:N; s = sin(2*pi*f1*n/fs

14、)+sin(2*pi*f2*n/fs);%s為原始信號(hào)x = awgn(s,10);%x為觀測(cè)信號(hào)，即對(duì)原始信號(hào)加入白噪聲，信噪比10dBfor i=1:N ef(1,i)=x(i); eb(1,i)=x(i);endsum=0;for i=1:N sum=sum+x(i)*x(i);endr(1)=sum/N;% Burg遞推for p=2:M% 求解第p個(gè)反射系數(shù) sum1=0; for n=p:N sum1=sum1+ef(p-1,n)*eb(p-1,n-1); end sum1=-2*sum1; sum2=0; for n=p:N sum2=sum2+ef(p-1,n)*ef(p-1,

15、n)+eb(p-1,n-1)*eb(p-1,n-1); end k(p-1)=sum1/sum2;% 求解預(yù)測(cè)誤差平均功率 r(p)=(1-k(p-1)*k(p-1)*r(p-1);% 求解p階白噪聲方差 q(p)=r(p);% 系數(shù)a if p>2 for i=1:p-2 a(p-1,i)=a(p-2,i)+k(p-1)*a(p-2,p-1-i); end end a(p-1,p-1)=k(p-1);% 求解前向預(yù)測(cè)誤差 for n=p+1:N ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1); end%求解后向預(yù)測(cè)誤差 for n=p:N-1 eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n); endend % 計(jì)算功率譜for j=1:N sum3=0; sum4=0; for i=1:p-1 sum3=sum3+a(p-1,i)*cos(2*pi*i*j/N); end sum3=1+sum3; for i=1:p-1 sum4=sum4+a(p-1,i)*