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;%系統輸出
r(n)=u'*u;%自相關矩陣
e(n) = d(n) - y(n) ;%誤差
fai=.0001;%修正參數，為防止u'*u過小導致步長值太大而設置的
if n < 20
mu=0.32;
else
mu=0.15;
end
w = w + mu * u * e(n)/(r(n)+fai) ;%迭代方程
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])