递归最小二乘法RLS功率谱估计MATLAB完整程序分享

递归最小二乘法RLS简介：

        RLS算法的主要问题之一是每次迭代中的计算量与阶数M的平方成正比。虽然比之最小二乘法（M的三次方成正比)好，但比LMS算法(M成正比)要差。

       RLS算法的优点：RLS算法与LMS算法的计算量实际上是一样的，因为RLS算法一次迭代M*M次，总共需要迭代2M次；RLS算法通过自相关矩阵的迭代计算，克服了非平稳随机信号对自相关矩阵估计误差的影响，因为RLS算法来一个样本递推一次自相关矩阵比你用若干个样本一次估计自相关矩阵的精确度要高。一个广泛的共识是RLS 算法的收敛速度和跟踪性能都优于 LMS 算法，所付出的代价是需要更复杂的计算 。

功率谱简介：

        随机信号是时域无限信号，不具备可积分条件，因此不能直接进行傅氏变换。一般用具有统计特性的功率谱来作为谱分析的依据。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲。所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。

        功率谱密度，从名字分解来看就是说，观察对象是功率，观察域是谱域，通常指频域，密度，就是指观察对象在观察域上的分布情况。一般我们讲的功率谱密度都是针对平稳随机过程的，由于平稳随机过程的样本函数一般不是绝对可积的，因此不能直接对它进行傅立叶分析。

1、产生数据

N = 1001;                           %样本点数
snr=[20 25 30];                     %信噪比
g=100;                              %蒙特卡诺仿真次数
n=0:N-1;                            %数据轴
MSE=zeros(N-1,3);                   %存放误差
M=10;                               %滤波器抽头数
signal1 = exp(1i*0.15*2*pi*n+1i*2*pi*rand);
signal2 = exp(1i*0.25*2*pi*n+1i*2*pi*rand);
signal3 = exp(1i*0.30*2*pi*n+1i*2*pi*rand);
UN = signal1 + signal2 + signal3;   %输入信号产生
w_RLS=zeros(M,N,3);                 %存放权系数

2、构造数据矩阵(扩展的数据矩阵)及嵌入RLS

for i=1:3
    for kk=1:g
        un=awgn(UN,snr(i),'measured');
       %% 构造数据矩阵(扩展的数据矩阵)
        Un=[zeros(1,M-1) un];       %数据扩展
        A=zeros(M,N-M+1);
        for nn=1:N
            A(:,nn)=Un(M+nn-1:-1:nn);
        end
       %% RLS算法
        delta=0.004;                %调整因子
        lamda=0.99;                 %遗忘因子
        dn=un(2:end);               %期望信号
        w=zeros(M,N);               %权向量初始化
        epsilon=zeros(N-1,1);       %误差初始化
        P1=eye(M)/delta;            %相关矩阵逆阵初始化
        for k=1:N-1
            pin=P1*A(:,k);                      %增益向量分子
            deno=lamda+A(:,k)'*pin;             %增益向量分母
            kn=pin/deno;                        %增益向量计算
            epsilon(k)=dn(k)-w(:,k)'*A(:,k);    %估计误差
            w(:,k+1)=w(:,k)+kn*conj(epsilon(k));%权向量更新
            P1=P1/lamda-kn*A(:,k)'*P1/lamda;    %相关矩阵逆阵更新
        end
        w_RLS(:,:,i)=w_RLS(:,:,i)+w;
        MSE(:,i)=MSE(:,i)+abs(epsilon).^2;
    end
end
w_RLS=w_RLS/g;
MSE=MSE/g;

3、绘图

figure;
hold on
plot(abs(w_RLS(1,:,1)));
plot(abs(w_RLS(2,:,1)));
hold off
figure
hold on
plot(MSE(1:200,1));
plot(MSE(1:200,2));
plot(MSE(1:200,3));
hold off
title('遗忘因子0.99下的学习曲线');
legend('SNR=20','SNR=25','SNR=30');
xlabel('迭代次数n');ylabel('均方误差MSE');
grid on

4、计算功率谱

sigma=sum(MSE(700:800,1)/100);          %白噪声方差
for k=1:M
    w_G(k)=sum(w_RLS(k,700:800,1))/100; %权值取平均
end
a=-conj(w_G);               %权值转AR参数
P=1024;                     %搜索点数
f=linspace(-0.5,0.5,P);     %频率轴
omega=2*pi*f;               %相对角频率
aw=zeros(M,P);              %构建频率矩阵
for k=1:P
    for m=1:M
        aw(m,k)=exp(-1j*omega(k)*m);
    end
end
Sx=zeros(size(f));          %功率谱
for m=1:P
    deno=abs(1+a*aw(:,m))^2;
    Sx(m)=sigma/deno;
end
Sx=10*log10(abs(Sx/max(abs(Sx))));
figure;
plot(f,Sx);