; This batch file creates a plot a bandstop filter which suppresses
; frequencies between 7 cycles per second and 15 cycles per second for
; data sampled every 0.02 seconds, using the Hanning window. It is used
; in Chapter 13, "Signal Processing", of _Using IDL_.

delt = 0.02 ; sampling period in seconds

f_low = 15. ; frequencies above f_low will be passed

f_high = 7. ; frequencies below f_high will be passed

nfilt = 81 ; the length of the filter

f_filt = FINDGEN(nfilt/2+1) / (nfilt*delt)

ideal_fr = (f_filt GT f_low) $ ; pass frequencies greater than f_low
        OR (f_filt LT F_high)  ; pass frequencies less than f_high

ideal_fr = FLOAT(ideal_fr) ; convert from byte to floating point

; replicate to obtain values for negative frequencies:

ideal_fr = [ideal_fr, REVERSE(ideal_fr(1:*))]

; now use an inverse FFT to get the impulse response of the ideal filter

ideal_ir = FLOAT(FFT(ideal_fr, /INVERSE)) ; ideal_fr is an even function,
                                          ; so the result is real
ideal_ir = ideal_ir / nfilt ; scale by the # of points

ideal_ir = SHIFT(ideal_ir, nfilt/2) ; shift it before applying the window

; apply a Hanning window to the shifted ideal impulse response

bs_ir_n = ideal_ir*HANNING(nfilt) ; these are the coefficients of the filter

; The frequency response of the filter is the FFT of its impulse response: 

bs_fr_n = FFT(bs_ir_n) * nfilt ; scale by the number of points

; log plot of magnitude in dB

mag = ABS(bs_fr_n(0:nfilt/2)) ; mag of Hanning bandstop filter x'fer f'n

PLOT, f_filt, 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 for Bandstop!CFIR Filter (Hanning)'