注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Mr.Right

不顾一切的去想,于是我们有了梦想。脚踏实地的去做,于是梦想成了现实。

 
 
 

日志

 
 
关于我

人生一年又一年,只要每年都有所积累,有所成长,都有那么一次自己认为满意的花开时刻就好。即使一时不顺,也要敞开胸怀。生命的荣枯并不是简单的重复,一时的得失不是成败的尺度。花开不是荣耀,而是一个美丽的结束,花谢也不是耻辱,而是一个低调的开始。

网易考拉推荐

RLS 自适应算法 -- IIR结构  

2012-03-23 12:58:40|  分类: 学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rlsalgo : IIR RLS algorithm demo
% Reference : S. Haykin, Adaptive Filter Theory. 3rd edition, Upper Saddle River, NJ: Prentice-Hall, 1996. 
% Note : The author doesn't take any responsibility for any harm caused by the use of this file
clear all
close all
hold off
% Number of system points
N=2000;
inp = randn(N,1);
n = randn(N,1);
[b,a] = butter(2,0.25);
Gz = tf(b,a,-1);

%If you don't have access to Control toolbox use the sample data
%load IIRsampledata;
%inp= IIRsampledata(1:2000);
%d= IIRsampledata(1:2000);

% use ldiv to get the approximate IIR weights of the filter ( a function only in the input)
% y=h*u
%ldiv is a function submitted to get inverse Z-transform (Matlab central file exchange)
%The first sysorder weight value
%use h=ldiv(b,a,sysorder)'; ==> here we use sysorder == 10

%channel system order (you can change the sysorder value and you don't need to change anything in the algorithm )
sysorder = 10 ;
%h= [0.0976   ; 0.2873  ;  0.3360   ; 0.2210   ; 0.0964   ; 0.0172 ;  -0.0159 ;  -0.0207  ; -0.0142  ; -0.0065 ; -0.0014 ;  0.0009 ;   0.0013   ; 0.0009   ; 0.0004  ;  0.0001  ; -0.0000  ; -0.0001  ; -0.0001  ; -0.0000];
h=[0.097631   0.287310   0.335965   0.220981 0.096354 0.017183  -0.015917 -0.020735  -0.014243  -0.006517 -0.001396   0.000856   0.001272  0.000914 0.000438 0.000108 -0.000044  -0.00008  -0.000058 -0.000029];
h=h(1:sysorder);
y = lsim(Gz,inp);
%add some noise
n = n * std(y)/(10*std(n));
d = y + n;
totallength=size(d,1);
%Take only 70 points for training ( N - systorder 70 = 80 - 10 )
N=80 ;
%begin of the algorithm
%forgetting factor
lamda = 0.9995 ;
%initial P matrix
delta = 1e10 ;  
P = delta * eye (sysorder ) ;
w = zeros ( sysorder  , 1 ) ;
for n = sysorder : N 
u = inp(n:-1:n-sysorder+1) ;
    phi = u' * P ;
k = phi'/(lamda + phi * u );
    y(n)=w' * u;
    e(n) = d(n) - y(n) ;
w = w + k * e(n) ;
P = ( P - k * phi ) / lamda ;
    % Just for plotting
    Recordedw(1:sysorder,n)=w;
end 
%check of results
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('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
figure
semilogy((abs(e))) ;
title('Error curve') ;
xlabel('Samples');
ylabel('Error value');
figure
plot(h, 'r+')
hold on
plot(w, '.')
legend('filter weights','Estimated filter weights');
title('Comparison of the filter weights and estimated weights') ;
figure
plot(Recordedw(1:sysorder,sysorder:N)');
title('Estimated weights convergence') ;
xlabel('Samples');
ylabel('Weights value');
axis([1 N-sysorder min(min(Recordedw(1:sysorder,sysorder:N)')) max(max(Recordedw(1:sysorder,sysorder:N)')) ]);
hold off
  评论这张
 
阅读(516)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016