; This batch file creates a plot of the impulse response and the 
; frequency response of a notch filter, as described in the section
; on Infinite Impulse Response filters in Chapter 13, "Signal Processing",
; of _Using IDL_.

@sigprc13     ; Load the coefficients for the notch filter.

na = N_ELEMENTS(A)-1 ; degree of denominator polynomial
nb = N_ELEMENTS(B)-1 ; degree of numerator polynomial
N = 1024L

; set the input u to an impulse

U = FLTARR(N)
U(0) = FLOAT(N)

Y = FLTARR(N)
Y(0) = B(2)*U(0) / A(na)

FOR K = 1, N-1 DO $ ; compute filtered signal recursively
   Y(K) = ( TOTAL( B(nb-K>0:nb  )*U(K-nb>0:K  ) ) $
          - TOTAL( A(na-K>0:na-1)*Y(K-na>0:K-1) ) ) / A(na)

V = FFT(Y) ; compute spectrum v

F = FINDGEN(N/2+1) / (N*delt) ; f = [0.0, 1.0/(N*delt), ... , 1.0/(2.0*delt)]
mag = ABS(V(0:N/2)); magnitude of first half of v
phi = ATAN(V(0:N/2)) ; phase of first half of v

; log plots of magnitude in dB and phase in degrees

!P.MULTI = [0, 1, 2, 0, 0] ; set up for two plots in window

PLOT, F, 20*ALOG10(mag), YTITLE='Magnitude in dB', $
    XTITLE='Frequency in cycles / second', $
    /XLOG, XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1, $
    TITLE='Frequency Response Function of b(z)/a(z)'

PLOT, F, phi/!DTOR, YTITLE='Phase in degrees', $
    YRANGE=[-180,180], YSTYLE=1, YTICKS=4, YMINOR=3, $
    XTITLE='Frequency in cycles / second', /XLOG, $
    XRANGE=[1.0,1.0/(2.0*delt)], XSTYLE=1

!P.MULTI = 0