求基于RLS算法和LMS的自适应均衡系统的MATLAB程序``

求基于RLS算法和LMS的自适应均衡系统的MATLAB程序``

里面有些代码有问题，可以参考，代码还是自己写：！
%基于RLS算法的自适应线性预测
clc;
clear all;
N=300;
M=100;%计算的次数
w1=zeros(N,M);w2=zeros(N,M);I=eye(2);e1=zeros(N,M);

for k=1:M
%产生白噪声
Pv=0.008;%定义白噪声方差
a1=-0.195;a2=0.95;o=0.02;r=0.95;
m=5000;%产生5000个随机数
v=randn(1,m);
v=v*sqrt(Pv);%产生均值为0，方差为Pv的白噪声
%m=1:N;
v=v(1:N);%取出前1000个
%plot(m,v);title('均值为0，方差为0.0965的白噪声');ylabel('v(n)');xlabel('n');
v=v';

%向量初使化
x=zeros(1,N);
x(1)=v(1);%x(0)=v(0)
x(2)=v(2)-a1*v(1);%x(1)=v(1)-a1*v(0)
w=zeros(2,N);
w(:,1)=[0 0]';%w(0)=[0 0]';
X=zeros(2,N);
X(:,2)=[v(1) 0]';%X(0)=[0 0]';X(1)=[v(0) 0]'
C=zeros(2,2*N);
C(:,1:2)=1/o.*I;%C(0)=1/o*I
e=zeros(1,N)';%定义误差向量
u=zeros(1,N);
g=zeros(2,N);

%根据RLS算法进行递推

for n=1:N-2
x(n+2)=v(n+2)-a1*x(n+1)-a2*x(n);
X(:,n+2)=[x(n+1) x(n)]';
u(n)=X(:,n+1)'*C(:,2*n-1:2*n)*X(:,n+1);
g(:,n)=(C(:,2*n-1:2*n)*X(:,n+1))./(r+u(n));
w(:,n+1)=w(:,n)+g(:,n)*(x(n+1)-X(:,n+1)'*w(:,n));
C(:,2*n+1:2*(n+1))=1/r.*(C(:,2*n-1:2*n)-g(:,n)*X(:,n+1)'*C(:,2*n-1:2*n));
e(n)=x(n+1)-X(:,n+1)'*w(:,n);

w1(:,k)=w(1,:)'; w2(:,k)=w(2,:)';%将每次计算得到的权矢量值储存
e1(:,k)=e(:,1);%将每次计算得到的误差储存
end

end
%求权矢量和误差的M次的平均值
wa1=zeros(N,1);wa2=zeros(N,1);en=zeros(N,1);
for k=1:M
wa1(:,1)=wa1(:,1)+w1(:,k);
wa2(:,1)=wa2(:,1)+w2(:,k);
en(:,1)=en(:,1)+e1(:,k);
end
n=1:N;
subplot(221)
plot(n,w(1,n),n,w(2,n));%作出单次计算权矢量的变化曲线
xlabel('n');ylabel('w(n)');title('w1(n)和w2(n)的单次变化曲线(线性预测，RLS)')
subplot(222)
plot(n,wa1(n,1)./M,n,wa2(n,1)./M);%作出100次计算权矢量的平均变化曲线
xlabel('n');ylabel('w(n)');title('w1(n)和w2(n)的100次平均变化曲线')
subplot(223)
plot(n,e(n,1).^2);%作出单次计算e^2的变化曲线
xlabel('n');ylabel('e^2');title('单次计算e^2的变化曲线');
subplot(224)
plot(n,(en(n,1)/M).^2);%作出M次计算e^2的平均变化曲线
xlabel('n');ylabel('e^2');title('100次计算e^2的平均变化曲线');

% rls算法 randn('seed', 0) ; rand('seed', 0) ; noofdata = 8000 ; % set no of data points used for training order = 32 ; % 自适应滤波权数 lambda = 0.98 ; % 遗忘因子 delta = 0.001 ; % 相关矩阵r的初始化 x = randn(noofdata, 1) ;%高斯随机系列 h = rand(order, 1) ; % 系统随机抽样 d = filter(h, 1, x) ; % 期望输出 % rls算法的初始化 p = delta * eye ( order, order ) ;%相关矩阵 w = zeros ( order, 1 ) ;%滤波系数矢量的初始化 % rls adaptation for n = order : noofdata ; u = x(n:-1:n-order+1) ;%延时函数 pi_ = u' * p ;%互相关函数 k = lambda + pi_ * u ; k = pi_'/k;%增益矢量 e(n) = d(n) - w' * u ;%误差函数 w = w + k * e(n) ;%递归公式 pprime = k * pi_ ; p = ( p - pprime ) / lambda ;%误差相关矩阵 w_err(n) = norm(h - w) ;%真实估计误差 end ; % 作图表示结果 figure ; plot(20*log10(abs(e))) ;%| e |的误差曲线 title('学习曲线') ; xlabel('迭代次数') ; ylabel('输出误差估计') ; figure ; semilogy(w_err) ;%作实际估计误差图 title('矢量估计误差') ; xlabel('迭代次数') ; ylabel('误差权矢量') ; %lms 算法 clear all close all hold off%系统信道权数 sysorder = 5 ;%抽头数 n=1000;%总采样次数 inp = randn(n,1);%产生高斯随机系列 n = randn(n,1); [b,a] = butter(2,0.25); gz = tf(b,a,-1);%逆变换函数 h= [0.0976;0.2873;0.3360;0.2210;0.0964;];%信道特性向量 y = lsim(gz,inp);%加入噪声 n = n * std(y)/(10*std(n));%噪声信号 d = y + n;%期望输出信号 totallength=size(d,1);%步长 n=60 ; %60节点作为训练序列 %算法的开始 w = zeros ( sysorder , 1 ) ;%初始化 for n = sysorder : n u = inp(n:-1:n-sysorder+1) ;% u的矩阵 y(n)= w' * u;%系统输出 e(n) = d(n) - y(n) ;%误差 if n < 20 mu=0.32; else mu=0.15; end w = w + mu * u * e(n) ;%迭代方程 end %检验结果 for n = n+1 : totallength u = inp(n:-1:n-sysorder+1) ; y(n) = w' * u ; e(n) = d(n) - y(n) ;%误差 end hold on plot(d) plot(y,'r'); title('系统输出') ; xlabel('样本') ylabel('实际输出') figure semilogy((abs(e))) ;% e的绝对值坐标 title('误差曲线') ; xlabel('样本') ylabel('误差矢量') figure%作图 plot(h, 'k+') hold on plot(w, 'r*') legend('实际权矢量','估计权矢量') title('比较实际和估计权矢量') ; axis([0 6 0.05 0.35])