function [err_norm, err_sq] = rls(input, output, N, hopt)
%	RLS adaptive filter for fig8_12.m
%	input	: input signal
%	output	: output signal
%	N	: no. of taps
%	hopt	: optimum set of weights
%	err_norm: normalised error of weights
%	err_sq	: squared estimation error
%
%	July 97

NN = max ( size (input));
h = zeros(size(hopt));
init = 100.0;
rinv = diag( ones(size(hopt)) * init );
for ii = 1:NN-N+1,
	x = input(ii:ii+N-1)';
	yhat = h*x;
	e = output(ii+N-1)-yhat;
	err_sq(ii+N-1) = e^2;
	c = rinv*x;
	rinv = rinv - c*c'/( 1+ c'*x);  
	h = h + (rinv*x*e)';
	herr = hopt-fliplr(h);
	err_norm(ii+N-1) = (herr*herr');
end