function filter,d,c,u,XSTATE=xstate,YSTATE=ystate,FREQ_RESP=freq_resp
;+
; y=filter(b,a,x) where a and b are the ARMA filter functions
; and x is a number sequence. The filter output is generated by
;
; y(n)=b(0)x(n)+b(1)x(n-1)+...+b(Q)x(n-Q)-a(1)y(n-1)-a(2)y(n-2)-...
;										  -a(P)y(n-P)
; This construction assumes that a(0)=1.
;
; INPUTS
; a=[1,a[1],...,a[P]] "autoregressive" coefficients
; b=[b[0]...b[Q]], "moving average" coefficients
; x=[x[0]...x[N-1]] input sequence of length N.
;
; The element a[0] is required but ignored. The value is set
; to a[0]=1 by the program. But don't omit it because then
; the program would, in effect, shorten the a vector.
;
; Pure Moving Average Filter: a=1 and b=[b[0]...b[Q]]
; Pure Autoregressive Filter: b=b0 and a=[1,a[1],...,a[P]]
;
; OUTPUTS
; The sequence y=[y[0]...y[N-1]].
;
; KEYWORD PARAMETERS
; XSTATE - set this parameter to a set of values that will
;          serve as the past input values. This is a vector
;          of length Q or less that represents the "past" values
;          of the input.
;
; YSTATE - A vector of up to P past output values used to
;          initialize the system.
;
; HISTORY
; Harvey E. Rhody  January, 1997, Initial version of program.
;                  July, 2000, Added XSTATE and YSTATE keywords.
;
;-

; Change the names of the parameters to protect their values.
b=d[*] & a=c[*] & x=u[*]
sb=size(b) & IF SB[0] EQ 0 THEN b=[b]
sa=size(a) & IF SA[0] eq 0 THEN a=[a]
P=N_ELEMENTS(a)-1
n=N_ELEMENTS(x)
Q=N_ELEMENTS(b)

IF NOT KEYWORD_SET(freq_resp) THEN BEGIN
;Generate the output sequence for the given input sequence.
IF (P GT 0) THEN BEGIN
	z=FLTARR(P) ; History of past outputs array.
	IF KEYWORD_SET(ystate) THEN BEGIN
		IF N_ELEMENTS(ystate) GT P THEN ystate=ystate[0:P-1]
		z[0:N_ELEMENTS(ystate)-1]=ystate
	ENDIF
	a=a[1:P]    ; Drop the a[0] element, assumed to be 1.
ENDIF

IF Q GT 1 THEN  BEGIN
	x=[FLTARR(Q-1),x] 	; Augment the input array with zeros
						; so the full output can be generated.
    IF KEYWORD_SET(xstate) THEN BEGIN ;Set MA Initial Conditions
       IF N_ELEMENTS(xstate) GE Q THEN xstate=xstate[0:Q-2]
       x[0:N_ELEMENTS(xstate)-1]=xstate
    ENDIF
ENDIF

y=fltarr(n) 	; Create an array for the output values.
;b=ROTATE(b,2)	; Flip b end for end
b=reverse(b)
b=TRANSPOSE(b)	; Make b a column vector
IF (P GT 1) THEN BEGIN
a=TRANSPOSE(a)  ; Make a into a column vector.
for k=0L,n-1 do begin
	y[k]=x[k:k+Q-1]##b-z##a
	z=[y[k],z[0:P-2]]
endfor
ENDIF ELSE $
IF P EQ 1 THEN BEGIN
	FOR k=0L,n-1 DO BEGIN
		y[k]=x[k:k+Q-1]##b-z*a
		z=y(k)
	ENDFOR
ENDIF ELSE BEGIN
FOR k=0L,n-1 DO BEGIN
	y[k]=x[k:k+Q-1]##b
ENDFOR
ENDELSE
return,y
ENDIF ELSE BEGIN
;Generate the frequency response. Now u contains the frequencies
;of interest for the plot.
nu=n_elements(u)
i=complex(0,1)
h=complexarr(nu)
for k=0L,nu-1 do begin
	z=exp(-i*u[k]*2*!pi)
	bsum=complex(0,0)
	asum=complex(0,0)
	for m=Q-1,0,-1 do bsum=bsum*z+b[m]
	for m=P,0,-1 do asum=asum*z+a[m]
	h[k]=bsum/asum
end
return,h
endelse
end