RLSQR-RLS算法的谱估计..docx

课程（论文）题目：仿真报告3内容：1 算法原理估计误差定义：可取滤波器的实际输入d*(i)作为期望响应d(i)。将误差代入代价函数得到加权误差平方和的完整表达式抽头权向量取的是n时刻的w(n)而不是i时刻的w(i)。为了使代价函数取得最小值,可通过对权向量求导解得：其中：由此可见指数加权最小二乘法的解转化为Wiener滤波器的形式:由令：由矩阵求逆引理得令其中k(n)为增益向量又由式中：指数加权RLS算法的步骤。1、初始化：w(0)=0，R(0)=σI,2、更新：对于n=1、2···计算：2 RLS算法谱估计仿真设计一个观测信号（至少包含三个主频分量），分别用RLS和QR-RLS完成信号的频谱估计，并比较两种算法。%产生零均值、方差为1的复高斯白噪声序列v(n)N=1000;noise=0.005*(randn(1,N)+j*randn(1,N))/sqrt(2);%产生带噪声的信号样本u(n)sig1=exp(j*0.31*2*pi*(0:N-1)+j*2*pi*rand) ;%产生第一个信号sig2=exp(-j*0.2*2*pi*(0:N-1)+j*2*pi*rand) ;%产生第二个信号sig3=exp(j*0.5*2*pi*(0:N-1)+j*2*pi*rand);%产生第三个信号Un=sig1+sig2+sig3+noise;%产生带噪声的信号unun=[zeros(1,M-1),Un].;%A=zeros(M,N);M=4;%滤波器抽头数N=1000;%样本数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+s;endUn=Un.;%构建矩阵A=zeros(M,N-M+1);%构建观测矩阵for n=1:N-M+1 A(:,n)=Un(M+n-1:-1:n);end[U,S,V]=svd(A);invphi=V*inv(S*S)*V;%构建矩阵phi%构建向量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];%扩展数据A=zeros(M,N);%构建样本矩阵for n=1:NA(:,n)=un(M+n-1 -1.n);enddelta=0.004;%调整参数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,M);for k=1:M %取后500个点的平均值 w(k)=sum(wopt(k,501:end))/500;enda=-conj(w);%AR模型的参数向量sigma=sum(MSE(501:end))/500;%AR模型输出白噪声方差%构建频率矩阵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=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);经过上述的步骤，即可估计出信号的