RLS、QR-RLS算法的譜估計(jì)

1、課程（論文）題目：仿真報(bào)告3內(nèi)容：1 算法原理估計(jì)誤差定義：可取濾波器的實(shí)際輸入d*(i)作為期望響應(yīng)d(i)。將誤差代入代價(jià)函數(shù)得到加權(quán)誤差平方和的完整表達(dá)式抽頭權(quán)向量取的是n時(shí)刻的w(n)而不是i時(shí)刻的w(i)。為了使代價(jià)函數(shù)取得最小值,可通過(guò)對(duì)權(quán)向量求導(dǎo)解得：其中：由此可見(jiàn)指數(shù)加權(quán)最小二乘法的解轉(zhuǎn)化為Wiener濾波器的形式:由令：由矩陣求逆引理得令其中k(n)為增益向量又由式中：指數(shù)加權(quán)RLS算法的步驟。1、初始化：w(0)=0，R(0)=I,2、更新：對(duì)于n=1、2計(jì)算：2 RLS算法譜估計(jì)仿真設(shè)計(jì)一個(gè)觀測(cè)信號(hào)（至少包含三個(gè)主頻分量），分別用RLS和QR-RLS完成信號(hào)的頻譜估計(jì)，并

2、比較兩種算法。%產(chǎn)生零均值、方差為1的復(fù)高斯白噪聲序列v(n)N=1000;noise=0.005*(randn(1,N)+j*randn(1,N)/sqrt(2);%產(chǎn)生帶噪聲的信號(hào)樣本u(n)sig1=exp(j*0.31*2*pi*(0:N-1)+j*2*pi*rand) ;%產(chǎn)生第一個(gè)信號(hào)sig2=exp(-j*0.2*2*pi*(0:N-1)+j*2*pi*rand) ;%產(chǎn)生第二個(gè)信號(hào)sig3=exp(j*0.5*2*pi*(0:N-1)+j*2*pi*rand);%產(chǎn)生第三個(gè)信號(hào)Un=sig1+sig2+sig3+noise;%產(chǎn)生帶噪聲的信號(hào)unun=zeros(1,M-1),

3、Un.;%A=zeros(M,N);M=4;%濾波器抽頭數(shù)N=1000;%樣本數(shù)f=0.1 0.25 0.27;%歸一化頻率SNR=30 30 27;%信噪比sigma=1;Am=sqrt(sigma*10.(SNR/10);%信號(hào)幅度t=linspace(0,1,N);phi=2*pi*rand(size(f);%隨機(jī)相位vn=sqrt(sigma/2)*randn(size(t)+j*sqrt(sigma/2)*randn(size(t);Un=vn;%加高斯白噪聲for k=1:length(f) s=Am(k)*exp(j*2*pi*N*f(k).*t+j*phi(k); Un=Un+

4、s;endUn=Un.;%構(gòu)建矩陣A=zeros(M,N-M+1);%構(gòu)建觀測(cè)矩陣for n=1:N-M+1 A(:,n)=Un(M+n-1:-1:n);endU,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);for k=1:Pfor m=1:M a(m,k)=exp(-j*omega(k)*(m-1);endendun=zeros(1,M-1)Un;%擴(kuò)展數(shù)據(jù)A=zeros(M,N);%構(gòu)建樣本矩陣for n=1:NA(:,n)=u

5、n(M+n-1 -1.n);enddelta=0.004;%調(diào)整參數(shù)lambda=0.98;%遺忘因子dn=Un(2:end);%一步預(yù)測(cè)期望信號(hào)w=zeros(M,N);epsilon=zeros(N-1,1);%先驗(yàn)估計(jì)誤差P1=eye(M)/delta;for k=1:N-1 %RLS算法迭代過(guò)程 PIn=P1*A(:,k); deno=lambda+A(:,k)*A(:,k); w(:,k+1)=w(:,k)+kn*conj(epsilon(k); P1=P1/lambda-kn*U(:,k)*P1/lambda;endMSE=abs(epsilon).2;%均方誤差w=zeros(1

6、,M);for k=1:M %取后500個(gè)點(diǎn)的平均值 w(k)=sum(wopt(k,501:end)/500;enda=-conj(w);%AR模型的參數(shù)向量sigma=sum(MSE(501:end)/500;%AR模型輸出白噪聲方差%構(gòu)建頻率矩陣P=1024;%將【0 2*pi】采樣1024點(diǎn)f=linspace(-0.5,0.5,P);%歸一化頻率omega=2*pi*f;%相對(duì)角頻率aw=zeros(M,P);for k=1:M aw(m,k)=exp(-j*omega(k)*(m);end%計(jì)算功率譜Sx=zeros(size(f);for m=1:P %計(jì)算功率譜過(guò)程 deno=

7、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)過(guò)上述的步驟，即可估計(jì)出信號(hào)的頻率和功率譜密度，如下圖所示：考慮一個(gè)一階AR模型：式中v(n)是方差為0.997的白噪聲信號(hào)，使用兩抽頭FIR濾波器，分別采用RLS和QR-RLS完成信號(hào)u(n)的線性預(yù)測(cè)，并比較結(jié)果。a1=0.99;%AR模型系數(shù)sigma=0.997;%白噪聲方差N=1000;

8、%數(shù)據(jù)個(gè)數(shù)vn=sqrt(sigma)*randn(N,1);%產(chǎn)生白噪聲樣本nume=1;%分子系數(shù)deno=1,a1;%分母系數(shù)u0=zeros(length(deno)-1,1);%初始數(shù)據(jù) xic=filtic(nume,deno,u0);%初始條件un=filter(nume,deno,vn,xic);%產(chǎn)生數(shù)據(jù)%產(chǎn)生期望響應(yīng)和觀測(cè)數(shù)據(jù)矩陣n0=1;%需要實(shí)現(xiàn)n0步線性預(yù)測(cè)M=2;%濾波器階數(shù)b=un(n0+1:N);%預(yù)測(cè)期望響應(yīng)L=length(b);un1=zeros(M-1,1).,un.;%擴(kuò)展數(shù)據(jù)A=zeros(M,L);for k=1:L %構(gòu)建觀測(cè)數(shù)據(jù)矩陣 A(:,k

9、)=un1(M-1+k:-1:k);end%應(yīng)用RLS算法進(jìn)行迭代求最優(yōu)權(quán)向量delta=0.004;%調(diào)整參數(shù)lambda=0.98;%遺忘因子w=zeros(M,L+1);%存取權(quán)向量epsilon=zeros(L,1);P1=eye(M)/delta;for k=1:L %RLS算法迭代過(guò)程 PIn=P1*A(:,k); denok=lambda+A(:,k)*A(:,k); kn=PIn/denok; epsilon(k)=b(k)-w(:,k)*A(:,k); w(:,k+1)=w(:,k)+kn*conj(epsilon(k); PI=PI/lambda-kn*A(:,k)*PI/lambda;endMSE=abs(epsilon).2;w_mean=w_mean/trials;%500次獨(dú)立實(shí)驗(yàn)權(quán)向量的均值?MSE_mean=MSE_mean/trials;%500次獨(dú)立實(shí)驗(yàn)的MSE?t=1:1000;figure(1)plot(t,w(1,:,