資源描述:
《RLS、QR-RLS算法的譜估計.docx》由會員上傳分享,免費在線閱讀,更多相關(guān)內(nèi)容在教育資源-天天文庫。
1、課程(論文)題目:仿真報告3內(nèi)容:1算法原理估計誤差定義:可取濾波器的實際輸入d*(i)作為期望響應(yīng)d(i)。將誤差代入代價函數(shù)得到加權(quán)誤差平方和的完整表達式抽頭權(quán)向量取的是n時刻的w(n)而不是i時刻的w(i)。為了使代價函數(shù)取得最小值,可通過對權(quán)向量求導(dǎo)解得:其中:由此可見指數(shù)加權(quán)最小二乘法的解轉(zhuǎn)化為Wiener濾波器的形式:由令:由矩陣求逆引理得令其中k(n)為增益向量又由式中:指數(shù)加權(quán)RLS算法的步驟。1、初始化:w(0)=0,R(0)=σI,2、更新:對于n=1、2···計算:2RLS算法譜估計仿真設(shè)計一個觀測信號(至少包含三個主頻分量),分別用R
2、LS和QR-RLS完成信號的頻譜估計,并比較兩種算法。%產(chǎn)生零均值、方差為1的復(fù)高斯白噪聲序列v(n)N=1000;noise=0.005*(randn(1,N)+j*randn(1,N))/sqrt(2);%產(chǎn)生帶噪聲的信號樣本u(n)sig1=exp(j*0.31*2*pi*(0:N-1)+j*2*pi*rand);%產(chǎn)生第一個信號sig2=exp(-j*0.2*2*pi*(0:N-1)+j*2*pi*rand);%產(chǎn)生第二個信號sig3=exp(j*0.5*2*pi*(0:N-1)+j*2*pi*rand);%產(chǎn)生第三個信號Un=sig1+sig2+s
3、ig3+noise;%產(chǎn)生帶噪聲的信號unun=[zeros(1,M-1),Un].';%A=zeros(M,N);M=4;%濾波器抽頭數(shù)N=1000;%樣本數(shù)f=[0.10.250.27];%歸一化頻率SNR=[303027];%信噪比sigma=1;Am=sqrt(sigma*10.^(SNR/10));%信號幅度t=linspace(0,1,N);phi=2*pi*rand(size(f));%隨機相位vn=sqrt(sigma/2)*randn(size(t))+j*sqrt(sigma/2)*randn(size(t));Un=vn;%加高斯白噪聲
4、fork=1:length(f)s=Am(k)*exp(j*2*pi*N*f(k).*t+j*phi(k));Un=Un+s;endUn=Un.';%構(gòu)建矩陣A=zeros(M,N-M+1);%構(gòu)建觀測矩陣forn=1:N-M+1A(:,n)=Un(M+n-1:-1:n);end[U,S,V]=svd(A');invphi=V*inv(S'*S)*V';%構(gòu)建矩陣phi%構(gòu)建向量a(w)P=1024;f=linspace(-0.5,0.5,P);omega=2*pi*f;a=zeros(M,P);fork=1:Pform=1:Ma(m,k)=exp(-j*o
5、mega(k)*(m-1));endendun=[zeros(1,M-1)Un']';%擴展數(shù)據(jù)A=zeros(M,N);%構(gòu)建樣本矩陣forn=1:NA(:,n)=un(M+n-1-1.n);enddelta=0.004;%調(diào)整參數(shù)lambda=0.98;%遺忘因子dn=Un(2:end);%一步預(yù)測期望信號w=zeros(M,N);epsilon=zeros(N-1,1);%先驗估計誤差P1=eye(M)/delta;fork=1:N-1%RLS算法迭代過程PIn=P1*A(:,k);deno=lambda+A(:,k)'*A(:,k);w(:,k+1)
6、=w(:,k)+kn*conj(epsilon(k));P1=P1/lambda-kn*U(:,k)'*P1/lambda;endMSE=abs(epsilon).^2;%均方誤差w=zeros(1,M);fork=1:M%取后500個點的平均值w(k)=sum(wopt(k,501:end))/500;enda=-conj(w);%AR模型的參數(shù)向量sigma=sum(MSE(501:end))/500;%AR模型輸出白噪聲方差%構(gòu)建頻率矩陣P=1024;%將【02*pi】采樣1024點f=linspace(-0.5,0.5,P);%歸一化頻率omega=
7、2*pi*f;%相對角頻率aw=zeros(M,P);fork=1:Maw(m,k)=exp(-j*omega(k)*(m));end%計算功率譜Sx=zeros(size(f));form=1:P%計算功率譜過程deno=abs(1+a*aw(:,m))^2;Sx(m)=sigma/deno;endSx=abs(Sx/max(abs(Sx)));%功率譜歸一化Sx=10*log10(Sx);Sx=abs(Sx/max(abs(Sx)));Sx=10*log10(Sx);kk=-511:512;plot(kk/1024,Sx);經(jīng)過上述的步驟,即可估計出信號
8、的頻率和功率譜密度,如下圖所示:考慮一個一階AR模型