rls_AR_pred.m
上传用户:kendun0711
上传日期:2007-06-03
资源大小:32k
文件大小:2k
源码类别:

技术管理

开发平台:

Matlab

  1. function [Wo, xp, alpha, e] = rls_AR_pred(Xi, Y, verbose, lambda, delta)
  2. % function [Wo, xp, alpha, e] = rls_AR_pred(Xi, Y, verbose, lambda, delta)
  3. %
  4. % rls_AR_pred.m - use basic RLS algorithm to predict real-valued AR process
  5. % written for MATLAB 4.0
  6. %
  7. % Reference: Haykin, _Adaptive Filter Theory_, 2nd (corr.) ed., 1991
  8. %
  9. % Note that we use the algorithm in Table 13.2, i.e.,
  10. % we do not exploit the Hermitian property of P(n), to
  11. % minimize the possibility of numerical instability.
  12. %
  13. %
  14. % Input parameters:
  15. %       Xi     : matrix of training/test points - each row is
  16. %                considered a datum
  17. %       y      : vector of corresponding desired outputs for
  18. %                predictor
  19. %       verbose: set to 1 for interactive processing
  20. %       lambda : the initial value of the forgetting factor
  21. %       delta  : the initial value for the diagonal P matrix
  22. %
  23. % Output parameters:
  24. %       Wo     : column-wise matrix of weights at each iteration
  25. %       xp     : row vector of predicted outputs
  26. %       alpha  : row vector of a priori prediction errors
  27. %       e      : row vector of a posteriori prediction errors Y - xp
  28. Nout = size(Xi, 2);
  29. % length of maximum number of timesteps that can be predicted
  30. N    = size(Xi, 1);
  31. % order of predictor
  32. d = size(Xi, 2);
  33. % initialize weight matrix and associated parameters for RLS predictor
  34. W  = zeros(Nout, 1);
  35. Wo = []; 
  36. P  = eye(d) / delta;
  37. for n = 1:N,
  38.    % save weights
  39. Wo = [Wo W];
  40. % adapt weight matrix
  41. u = Xi(n, :)';
  42. p = u' * P;
  43. kappa = lambda + p * u;
  44. k = P * u / kappa;
  45. alpha(n) = Y(n) - W' * u;
  46. W  = W + k * alpha(n);
  47. Pp = k * p;
  48. P  = (P - Pp) / lambda;
  49. % predict next sample and compute error
  50. xp(n) = W' * u;
  51. e(n)  = Y(n) - xp(n);
  52. if (verbose ~= 0)
  53.   disp(['time step ', int2str(n), ': mag. pred. err. = ', num2str(abs(e(n)))]);
  54. end;
  55. end % for n