RLS、QR-RLS算法的谱估计

1、课程（论文）题目：仿真报告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完成信号的频谱估计，并

2、比较两种算法。%产生零均值、方差为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),

3、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+

4、s;endUn=Un.;%构建矩阵A=zeros(M,N-M+1);%构建观测矩阵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;%构建矩阵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)=u

5、n(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

6、,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=

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);经过上述的步骤，即可估计出信号的频率和功率谱密度，如下图所示：考虑一个一阶AR模型：式中v(n)是方差为0.997的白噪声信号，使用两抽头FIR滤波器，分别采用RLS和QR-RLS完成信号u(n)的线性预测，并比较结果。a1=0.99;%AR模型系数sigma=0.997;%白噪声方差N=1000;

8、%数据个数vn=sqrt(sigma)*randn(N,1);%产生白噪声样本nume=1;%分子系数deno=1,a1;%分母系数u0=zeros(length(deno)-1,1);%初始数据 xic=filtic(nume,deno,u0);%初始条件un=filter(nume,deno,vn,xic);%产生数据%产生期望响应和观测数据矩阵n0=1;%需要实现n0步线性预测M=2;%滤波器阶数b=un(n0+1:N);%预测期望响应L=length(b);un1=zeros(M-1,1).,un.;%扩展数据A=zeros(M,L);for k=1:L %构建观测数据矩阵 A(:,k

9、)=un1(M-1+k:-1:k);end%应用RLS算法进行迭代求最优权向量delta=0.004;%调整参数lambda=0.98;%遗忘因子w=zeros(M,L+1);%存取权向量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次独立实验权向量的均值?MSE_mean=MSE_mean/trials;%500次独立实验的MSE?t=1:1000;figure(1)plot(t,w(1,:,