%最小均方差算法
% RLS 算法 
randn('seed', 0) ;
rand('seed', 0) ;
NoOfData = 8000 ;  % 設置序列的數據點數
Order = 32 ;    % 設置自適應濾波器的階數
Lambda = 0.98 ;    % 設定因數
Delta = 0.001 ;    % 已初始化的系數
x = randn(NoOfData, 1) ;% 假定輸入
h = rand(Order, 1) ;  % 系統任意選擇的序列
d = filter(h, 1, x) ;  % 產生輸出
% RLS初始化
P = Delta * eye ( Order, Order ) ;
w = zeros ( Order, 1 ) ;
% RLS 自適應
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))) ;
title('Learning Curve') ;%學習曲線
xlabel('Iteration Number') ;
ylabel('Output Estimation Error in dB') ;
pause;
figure ;
semilogy(w_err) ;
title('Weight Estimation Error') ;%誤差函數
xlabel('Iteration Number') ;
ylabel('Weight Error in dB') ;
pause;

clear all
close all
hold off
%系統頻道數
sysorder = 2 ; snr=30 ;
% 系統點數
N=2000;
inp = randn(N,1);
n = randn(N,1)/10.^(snr/10);
b=[1,0,0]; a=[1,-1.558,0.81];
Gz = tf(b,a,1);
h= [1.558;-0.81];
y = lsim(Gz, inp);% inp通過gz的輸出響應
%加噪聲
d = y + n;
totallength=size(d,1);
%算法開始
M=totallength;
%初始化各系數
f0=zeros(1,M);f1=zeros(1,M);f2=zeros(1,M);
k1=zeros(1,M);k2=zeros(1,M);
b0=zeros(1,M);b1=zeros(1,M);b2=zeros(1,M);
a1=zeros(1,M);a2=zeros(1,M);
f0=inp;b0=inp;
f1(1)=inp(1);
k1(1)=0;k2(1)=0;b1(1)=0;
k1(2)=0;f2(1)=inp(1);
b2(1)=0;k2(2)=0;
u=0.05;
for n = sysorder : M-1
 i=n+2-sysorder;
 f1(i)=f0(i)+k1(i)*b0(i-1);
 b1(i)=b0(i-1)+k1(i)*f0(i);
 k1(i+1)=k1(i)-2*u*(f1(i)*b0(i-1)+b1(i)*f0(i));
 f2(i)=f1(i)+k2(i)*b1(i-1);
 b2(i)=b1(i-1)+k2(i)*f1(i);
 k2(i+1)=k2(i)-2*u*(f2(i)*b1(i-1)+b2(i)*f1(i));
end 
for n = sysorder:M
i=n+1-sysorder;
a1(i)=(-k1(i)*(1+k2(i)));a2(i)=-k2(i);
end
hold on
figure(1); subplot(3,1,1)
plot(d)
plot(y, 'r');
title('System output') ;%系統輸出
xlabel('Samples')
ylabel('True and estimated output')
subplot(3,1,3)
plot(h, 'k+')
hold on
plot(a,'r*')
legend('Actual weights',' Estimated weights')
title('Comparison of the actual weights and the estimated weights') ;%真實加權值與估計加權值的比較
axis([0 3 -2.35 2.35])