FUNCTION filter_ss,b,a,u,STATES=xhist,X0=x0
;+
;y=filter_ss(b,a,u) uses a state-space model to compute the
;output sequence of a filter that is described by the
;rational function Y(z)=(B(z)/A(z))U(z).
;
;See documentation to filter.
;
;KEYWORDS
;y=filter_ss(b,a,u,STATES=x) returns an array in which column
;n is the state at step n. The first column of x is the
;initial state. y(n)=H##x(n+1) since the state is updated
;prior to computing the output.
;
;X0=x0 is the initial state, which must be a column vector
;of length N, where N is the order of polynomial a.
;
;This implementation assumes N=n_elements(a) GE
;M=n_elements(b).
;
;H. Rhody
;February 4, 2005
;April 24, 2005 Revised to properly handle b longer than a
;-

;Construct the transition matrix F such that x(n+1)=F##x(n)+u(n)
;where x(n) is the state at step n.

P=n_elements(a)
Q=n_elements(b)
N=P>(Q+1)-1 ;Order of the filter
F=fltarr(N,N)
IF P GE 2 THEN F[0:P-2,0]=-a[1:*]
IF N GE 2 THEN F[0:N-2,1:N-1]=Identity(N-1)

;Construct the output matrix H such that y=H##x
H=(Q GE N) ? b : [b,fltarr(N-Q)]

;Initialize the state vector.
x=Keyword_Set(X0) ? x0 : fltarr(1,N)
xhist=x

;Do the iteration over the input sequence.
for k=0L,n_elements(u)-1 do begin
	x=F##x                         ;Step the state vector
	x[0]=x[0]+u[k]                 ;Add the input value.
	y=(k EQ 0) ? [H##x] : [y,H##x] ;Compute output from current state
	xhist=[xhist,x]
end

return,y
end