功率譜密度估計方法的MATLAB實現(xiàn)(共13頁)

1、 15/15功率譜密度估計方法(fngf)的MATLAB實現(xiàn)在應(yīng)用數(shù)學(xué)和物理學(xué)中，譜密度、功率譜密度和能量譜密度是一個用于信號的通用概念，它表示每赫茲的功率、每赫茲的能量這樣的物理量綱。在物理學(xué)中，信號通常是波的形式，例如電磁波、隨機振動或者聲波。當(dāng)波的頻譜密度乘以一個適當(dāng)?shù)南禂?shù)后將得到每單位頻率波攜帶的功率，這被稱為信號的功率譜密度（power spectral density, PSD）或者譜功率分布（spectral power distribution, SPD）。功率譜密度的單位通常用每赫茲的瓦特數(shù)（W/Hz）表示，或者使用波長而不是頻率，即每納米的瓦特數(shù)（W/nm）來表示。信號的功

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

3、3）周期圖法譜估計程序源程序Fs=100000; %采樣頻率100kHzN=1024; %數(shù)據(jù)長度N=1024n=0:N-1;t=n/Fs;xn=sin(2000*2*pi*t); %正弦波,f=2000HzY=awgn(xn,10); %加入(jir)信噪比為10db的高斯白噪聲subplot(2,1,1);plot(n,Y) title(信號(xnho)xlabel(時間(shjin);ylabel(幅度);grid on;window=boxcar(length(xn); %矩形窗nfft=N/4; %采樣點數(shù)Pxx f=periodogram(Y,window,nfft,Fs); %直

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

5、);plot(n,Y) title(信號)xlabel(時間);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(修正周期圖法譜估計 N=,int2str(N), M=,int2str(M);xlabel(頻率（Hz）);ylabel(功率譜密度);仿真結(jié)果最大熵法譜估計程序(chngx)1、源程序fs=1; %設(shè)采樣(ci yn)頻率N=128; %數(shù)據(jù)長度(chngd) 改變數(shù)據(jù)長度會導(dǎo)致分

6、辨率的變化；f1=0.2*fs; %第一個sin信號的頻率,f1/fs=0.2f2=0.3*fs; %第二個sin信號的頻率，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為原始信號x=awgn(s,10); %x為觀測信號，即對原始信號加入白噪聲，信噪比10dBfigure(1); %畫出原始信號和觀測信號subplot(2,1,1);plot(s,b),xlabel(時間),ylabel(幅度),title(原始信號s);grid;subplot(2,1,2);plot(x,r),xla

7、bel(時間),ylabel(幅度),title(觀測信號x);Pxx1,f=pmem(x,P,N,fs); %最大熵譜估計figure(2);plot(f,10*log10(Pxx1);xlabel(頻率(pnl)(Hz) );ylabel(功率譜(dB) );title(最大熵法譜估計 模型(mxng)階數(shù)P=,int2str(P), 數(shù)據(jù)長度N=,int2str(N);仿真(fn zhn)結(jié)果Levinson遞推法譜估計程序(chngx)源程序fs=1; %設(shè)采樣(ci yn)頻率為1N=1000; %數(shù)據(jù)長度 改變數(shù)據(jù)長度會導(dǎo)致(dozh)分辨率的變化；f1=0.2*fs; %第一個s

8、in信號的頻率,f1/fs=0.2f2=0.3*fs; %第二個sin信號的頻率，f1/fs=0.2或者0.3M=16; %濾波器階數(shù)的最大取值,超過則認(rèn)為代價太大而放棄L=2*N; %有限長序列進行離散傅里葉變換前，序列補零的長度n=1:N; s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s為原始信號x=awgn(s,10);%x為觀測信號，即對原始信號加入白噪聲，信噪比10dBfigure(1); %畫出原始信號和觀測信號subplot(2,1,1);plot(s,b),axis(0 100 -3 3),xlabel(時間),ylabel(幅度),title

9、(原始信號s);grid;subplot(2,1,2);plot(x,r),axis(0 100 -3 3),xlabel(時間(shjin),ylabel(幅度),title(觀測信號x);grid;%計算自相關(guān)(xinggun)函數(shù)rxx = xcorr(x,x,M,biased);%計算有偏估計自相關(guān)函數(shù)(hnsh)，長度為-M到M，%共2M+1r0 = rxx(M+1); %r0為零點上的自相關(guān)函數(shù)，相對于-M，第M+1個點為零點R = rxx(M+2:2*M+1);% R為從1到第M個點的自相關(guān)函數(shù)矩陣%確定矩陣大小a = zeros(M,M);FPE = zeros(1,M);%F

10、PE：最終預(yù)測誤差，用來估計模型的階次var = zeros(1,M);%求初值a(1,1) = -R(1)/r0;%一階模型參數(shù)var(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

11、)*var(p-1); %求方差 FPE(p)=var(p)*(M+1+p)/(M+1-p);%求最終預(yù)測誤差end %確定(qudng)AR模型的最佳階數(shù)min=FPE(1); %求出FPE最小時(xiosh)對應(yīng)的階數(shù)p = 1;for k=2:M if FPE(k)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);% 求解(qi ji)前向預(yù)測誤差 for n=p+1:N ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1); end%求解后向預(yù)測誤差 for

12、n=p:N-1 eb(p,n)=eb(p-1,n-1)+k(p-1)*ef(p-1,n); endend % 計算功率譜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)*sin(2*pi*i*j/N); end pxx=sqrt(sum3*sum3+sum4*sum4); pxx=q(M)/pxx; pxx=10*log10(pxx); pp(j)=pxx;end%畫出功率(gngl)譜ff=1:N;ff=ff/N;figure;plot(ff,pp),axis(0 0.5 -20 10),xlabel(頻率),ylabel(幅度(fd)（dB）)