RLS、QR-RLS算法的譜估計

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

2、比較兩種算法。%產(chǎn)生零均值、方差為1的復高斯白噪聲序列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+sig3+noise;%產(chǎn)生帶噪聲的信號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);%信號幅度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;%加高斯白噪聲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)建觀測矩陣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;%擴展數(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);%一步預測期望信號w=zeros(M,N);epsilon=zeros(N-1,1);%先驗估計誤差P1=eye(M)/delta;for k=1:N-1 %RLS算法迭代過程 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個點的平均值 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點f=linspace(-0.5,0.5,P);%歸一化頻率omega=2*pi*f;%相對角頻率aw=zeros(M,P);for k=1:M aw(m,k)=exp(-j*omega(k)*(m);end%計算功率譜Sx=zeros(size(f);for m=1:P %計算功率譜過程 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)過上述的步驟，即可估計出信號的頻率和功率譜密度，如下圖所示：考慮一個一階AR模型：式中v(n)是方差為0.997的白噪聲信號，使用兩抽頭FIR濾波器，分別采用RLS和QR-RLS完成信號u(n)的線性預測，并比較結(jié)果。a1=0.99;%AR模型系數(shù)sigma=0.997;%白噪聲方差N=1000;

8、%數(shù)據(jù)個數(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)生期望響應和觀測數(shù)據(jù)矩陣n0=1;%需要實現(xiàn)n0步線性預測M=2;%濾波器階數(shù)b=un(n0+1:N);%預測期望響應L=length(b);un1=zeros(M-1,1).,un.;%擴展數(shù)據(jù)A=zeros(M,L);for k=1:L %構(gòu)建觀測數(shù)據(jù)矩陣 A(:,k

9、)=un1(M-1+k:-1:k);end%應用RLS算法進行迭代求最優(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算法迭代過程 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次獨立實驗權(quán)向量的均值?MSE_mean=MSE_mean/trials;%500次獨立實驗的MSE?t=1:1000;figure(1)plot(t,w(1,:,