C XAPIIR -- SUBROUTINE: IIR FILTER DESIGN AND IMPLEMENTATION C C ARGUMENTS: C ---------- C C DATA REAL ARRAY CONTAINING SEQUENCE TO BE FILTERED C ORIGINAL DATA DESTROYED, REPLACED BY FILTERED DATA C C NSAMPS NUMBER OF SAMPLES IN DATA C C C APROTO CHARACTER*8 VARIABLE, CONTAINS TYPE OF ANALOG C PROTOTYPE FILTER C '(BU)TTER ' -- BUTTERWORTH FILTER C '(BE)SSEL ' -- BESSEL FILTER C 'C1 ' -- CHEBYSHEV TYPE I C 'C2 ' -- CHEBYSHEV TYPE II C C TRBNDW TRANSITION BANDWIDTH AS FRACTION OF LOWPASS C PROTOTYPE FILTER CUTOFF FREQUENCY. USED C ONLY BY CHEBYSHEV FILTERS. C C A ATTENUATION FACTOR. EQUALS AMPLITUDE C REACHED AT STOPBAND EDGE. USED ONLY BY C CHEBYSHEV FILTERS. C C IORD ORDER (#POLES) OF ANALOG PROTOTYPE C NOT TO EXCEED 10 IN THIS CONFIGURATION. 4 - 5 C SHOULD BE AMPLE. C C TYPE CHARACTER*8 VARIABLE CONTAINING FILTER TYPE C 'LP' -- LOW PASS C 'HP' -- HIGH PASS C 'BP' -- BAND PASS C 'BR' -- BAND REJECT C C FLO LOW FREQUENCY CUTOFF OF FILTER (HERTZ) C IGNORED IF TYPE = 'LP' C C FHI HIGH FREQUENCY CUTOFF OF FILTER (HERTZ) C IGNORED IF TYPE = 'HP' C C TS SAMPLING INTERVAL (SECONDS) C C PASSES INTEGER VARIABLE CONTAINING THE NUMBER OF PASSES C 1 -- FORWARD FILTERING ONLY C 2 -- FORWARD AND REVERSE (I.E. ZERO PHASE) FILTERING C C C SUBPROGRAMS REFERENCED: BILIN2, BUROOTS, WARP, CUTOFFS, LPTHP, LPTBP C LP, LPTBR, BEROOTS, C1ROOTS, C2ROOTS, CHEBPARM, DESIGN, APPLY C SUBROUTINE XAPIIR( DATA, NSAMPS, APROTO, TRBNDW, A, IORD, + TYPE, FLO, FHI, TS, PASSES ) include 'mach' C DIMENSION DATA(1) CHARACTER*8 TYPE, APROTO INTEGER NSAMPS, PASSES, IORD c REAL*4 TRBNDW, A, FLO, FHI, TS, SN(30), SD(30) REAL TRBNDW, A, FLO, FHI, TS, SN(30), SD(30) LOGICAL ZP C C Filter designed C CALL DESIGN( IORD, TYPE(1:2), APROTO(1:2), A, TRBNDW, & FLO, FHI, TS, SN, SD, NSECTS ) C C Filter data C IF ( PASSES .EQ. 1 ) THEN ZP = .FALSE. ELSE ZP = .TRUE. END IF CALL APPLY( DATA, NSAMPS, ZP, SN, SD, NSECTS ) C RETURN END C C Subroutine to design IIR digital filters from analog prototypes. C C Input Arguments: C ---------------- C C IORD Filter order (10 MAXIMUM) C C TYPE Character*2 variable containing filter type C LOWPASS (LP) C HIGHPASS (HP) C BANDPASS (BP) C BANDREJECT (BR) C C APROTO Character*2 variable designating analog prototype C Butterworth (BU) C Bessel (BE) C Chebyshev Type I (C1) C Chebyshev Type II (C2) C C A Chebyshev stopband attenuation factor C C TRBNDW Chebyshev transition bandwidth (fraction of C lowpass prototype passband width) C C FL Low-frequency cutoff C C FH High-frequency cutoff C C TS Sampling interval (in seconds) C C Output Arguments: C ----------------- C C SN Array containing numerator coefficients of C second-order sections packed head-to-tail. C C SD Array containing denominator coefficients C of second-order sections packed head-to-tail. C C NSECTS Number of second-order sections. C C SUBROUTINE DESIGN( IORD, TYPE, APROTO, A, TRBNDW, & FL, FH, TS, SN, SD, NSECTS ) C COMPLEX P(10), Z(10) CHARACTER*2 TYPE, APROTO CHARACTER*3 STYPE(10) c REAL*4 SN(1), SD(1) REAL SN(1), SD(1),FHW,FH,TS,DCVALUE C C Analog prototype selection C IF ( APROTO .EQ. 'BU' ) THEN C CALL BUROOTS( P, STYPE, DCVALUE, NSECTS, IORD ) C END IF C C Analog mapping selection C C IF (TYPE.EQ.'LP') THEN C FHW = WARP(FH*TS/2.,2.) CALL LP( P, Z, STYPE, DCVALUE, NSECTS, SN, SD ) CALL CUTOFFS( SN, SD, NSECTS, FHW ) C END IF C C Bilinear analog to digital transformation C CALL BILIN2( SN, SD, NSECTS ) C RETURN END C C WARP -- FUNCTION, APPLIES TANGENT FREQUENCY WARPING TO COMPENSATE C FOR BILINEAR ANALOG -> DIGITAL TRANSFORMATION C C ARGUMENTS: C ---------- C C F ORIGINAL DESIGN FREQUENCY SPECIFICATION (HERTZ) C TS SAMPLING INTERVAL (SECONDS) C C LAST MODIFIED: SEPTEMBER 20, 1990 C REAL FUNCTION WARP( F , TS ) C REAL F,TS,TWOPI,ANGLE TWOPI = 6.2831853 ANGLE = TWOPI*F*TS/2. WARP = 2.*TAN(ANGLE)/TS WARP = WARP/TWOPI C RETURN END C C Subroutine to generate second order section parameterization C from an pole-zero description for lowpass filters. C C Input Arguments: C ---------------- C C P Array containing poles C C Z Array containing zeros C C RTYPE Character array containing root type information C (SP) Single real pole or C (CP) Complex conjugate pole pair C (CPZ) Complex conjugate pole and zero pairs C C DCVALUE Zero-frequency value of prototype filter C C NSECTS Number of second-order sections C C Output Arguments: C ----------------- C C SN Numerator polynomials for second order C sections. C C SD Denominator polynomials for second order C sections. C C SUBROUTINE LP( P, Z, RTYPE, DCVALUE, NSECTS, SN, SD ) C COMPLEX P(1), Z(1) CHARACTER*3 RTYPE(1) c REAL*4 SN(1), SD(1), DCVALUE REAL SN(1), SD(1), DCVALUE C IPTR = 1 DO 1 I = 1, NSECTS C IF ( RTYPE(I) .EQ. 'CPZ' ) THEN C SCALE = REAL( P(I) * CONJG( P(I) ) ) & / REAL( Z(I) * CONJG( Z(I) ) ) SN( IPTR ) = REAL( Z(I) * CONJG( Z(I) ) ) * SCALE SN( IPTR + 1 ) = -2. * REAL( Z(I) ) * SCALE SN( IPTR + 2 ) = 1. * SCALE SD( IPTR ) = REAL( P(I) * CONJG( P(I) ) ) SD( IPTR + 1 ) = -2. * REAL( P(I) ) SD( IPTR + 2 ) = 1. IPTR = IPTR + 3 C ELSE IF ( RTYPE(I) .EQ. 'CP' ) THEN C SCALE = REAL( P(I) * CONJG( P(I) ) ) SN( IPTR ) = SCALE SN( IPTR + 1 ) = 0. SN( IPTR + 2 ) = 0. SD( IPTR ) = REAL( P(I) * CONJG( P(I) ) ) SD( IPTR + 1 ) = -2. * REAL( P(I) ) SD( IPTR + 2 ) = 1. IPTR = IPTR + 3 C ELSE IF ( RTYPE(I) .EQ. 'SP' ) THEN C SCALE = -REAL( P(I) ) SN( IPTR ) = SCALE SN( IPTR + 1 ) = 0. SN( IPTR + 2 ) = 0. SD( IPTR ) = -REAL( P(I) ) SD( IPTR + 1 ) = 1. SD( IPTR + 2 ) = 0. IPTR = IPTR + 3 C END IF C 1 CONTINUE C SN(1) = DCVALUE * SN(1) SN(2) = DCVALUE * SN(2) SN(3) = DCVALUE * SN(3) C RETURN END C CUTOFFS C C Subroutine to alter the cutoff of a filter. Assumes that the C filter is structured as second order sections. Changes C the cutoffs of normalized lowpass and highpass filters through C a simple polynomial transformation. C C C Input Arguments: C ---------------- C C F New cutoff frequency C C Input/Output Arguments: C ----------------------- C C SN Numerator polynomials for second order C sections. C C SD Denominator polynomials for second order C sections. C C NSECTS Number of second order sectionsects C C SUBROUTINE CUTOFFS( SN, SD, NSECTS, F ) C c REAL*4 SN(1), SD(1) REAL SN(1), SD(1),SCALE,F C SCALE = 2.*3.14159265*F C IPTR = 1 DO 1 I = 1, NSECTS C SN( IPTR + 1 ) = SN( IPTR + 1 ) / SCALE SN( IPTR + 2 ) = SN( IPTR + 2 ) / (SCALE*SCALE) SD( IPTR + 1 ) = SD( IPTR + 1 ) / SCALE SD( IPTR + 2 ) = SD( IPTR + 2 ) / (SCALE*SCALE) IPTR = IPTR + 3 C 1 CONTINUE C RETURN END C C C Transforms an analog filter to a digital filter via the bilinear transformati C Assumes both are stored as second order sections. The transformation is C done in-place. C C Input Arguments: C ---------------- C C SN Array containing numerator polynomial coefficients for C second order sections. Packed head-to-tail. C C SD Array containing denominator polynomial coefficients f C second order sections. Packed head-to-tail. C C NSECTS Number of second order sections. C C SUBROUTINE BILIN2( SN, SD, NSECTS ) C c REAL*4 SN(1), SD(1) REAL SN(1), SD(1),A0,A1,A2,SCALE C IPTR = 1 DO 1 I = 1, NSECTS C A0 = SD(IPTR) A1 = SD(IPTR+1) A2 = SD(IPTR+2) C SCALE = A2 + A1 + A0 SD(IPTR) = 1. SD(IPTR+1) = (2.*(A0 - A2)) / SCALE SD(IPTR+2) = (A2 - A1 + A0) / SCALE C A0 = SN(IPTR) A1 = SN(IPTR+1) A2 = SN(IPTR+2) C SN(IPTR) = (A2 + A1 + A0) / SCALE SN(IPTR+1) = (2.*(A0 - A2)) / SCALE SN(IPTR+2) = (A2 - A1 + A0) / SCALE C IPTR = IPTR + 3 C 1 CONTINUE C RETURN END C C BUROOTS -- SUBROUTINE TO COMPUTE BUTTERWORTH POLES FOR C NORMALIZED LOWPASS FILTER C C LAST MODIFIED: SEPTEMBER 7, 1990 C C OUTPUT ARGUMENTS: C ----------------- C P COMPLEX ARRAY CONTAINING POLES C CONTAINS ONLY ONE FROM EACH C COMPLEX CONJUGATE PAIR, AND C ALL REAL POLES C C RTYPE CHARACTER ARRAY INDICATING 2ND ORDER SECTION C TYPE: C (SP) SINGLE REAL POLE C (CP) COMPLEX CONJUGATE POLE PAIR C (CPZ) COMPLEX CONJUGATE POLE-ZERO PAIRS C C DCVALUE MAGNITUDE OF FILTER AT ZERO FREQUENCY C C NSECTS NUMBER OF SECOND ORDER SECTIONS C C INPUT ARGUMENTS: C ---------------- C C IORD DESIRED FILTER ORDER C C SUBROUTINE BUROOTS( P, RTYPE, DCVALUE, NSECTS, IORD ) C COMPLEX P(1) INTEGER HALF CHARACTER*3 RTYPE(1) C PI=3.14159265 C HALF = IORD/2 C C TEST FOR ODD ORDER, AND ADD POLE AT -1 C NSECTS = 0 IF ( 2*HALF .LT. IORD ) THEN P(1) = CMPLX( -1., 0. ) RTYPE(1) = 'SP' NSECTS = 1 END IF C DO 1 K = 1, HALF ANGLE = PI * ( .5 + FLOAT(2*K-1) / FLOAT(2*IORD) ) NSECTS = NSECTS + 1 P(NSECTS) = CMPLX( COS(ANGLE), SIN(ANGLE) ) RTYPE(NSECTS) = 'CP' 1 CONTINUE C DCVALUE = 1.0 C RETURN END C APPLY C Subroutine to apply an iir filter to a data sequence. C The filter is assumed to be stored as second order sections. C Filtering is in-place. C Zero-phase (forward and reverse) is an option. C C Input Arguments: C ---------------- C C DATA Array containing data C C NSAMPS Number of data samples C C ZP Logical variable, true for C zero phase filtering, false C for single pass filtering C C SN Numerator polynomials for second C order sections. C C SD Denominator polynomials for second C order sections. C C NSECTS Number of second-order sections C C Output Arguments: C ----------------- C C DATA Data array (same as input) C C SUBROUTINE APPLY( DATA, NSAMPS, ZP, SN, SD, NSECTS ) C C REAL*4 SN(1), SD(1), DATA(1) REAL SN(1), SD(1), DATA(1) C REAL*4 OUTPUT REAL OUTPUT,X1,X2,Y1,Y2,B0,B1,B2,A1,A2 LOGICAL ZP C JPTR = 1 DO 1 J = 1, NSECTS C X1 = 0.0 X2 = 0.0 Y1 = 0.0 Y2 = 0.0 B0 = SN(JPTR) B1 = SN(JPTR+1) B2 = SN(JPTR+2) A1 = SD(JPTR+1) A2 = SD(JPTR+2) C DO 2 I = 1, NSAMPS C OUTPUT = B0*DATA(I) + B1*X1 + B2*X2 OUTPUT = OUTPUT - ( A1*Y1 + A2*Y2 ) Y2 = Y1 Y1 = OUTPUT X2 = X1 X1 = DATA(I) DATA(I) = OUTPUT C 2 CONTINUE C JPTR = JPTR + 3 C 1 CONTINUE C IF ( ZP ) THEN C JPTR = 1 DO 3 J = 1, NSECTS C X1 = 0.0 X2 = 0.0 Y1 = 0.0 Y2 = 0.0 B0 = SN(JPTR) B1 = SN(JPTR+1) B2 = SN(JPTR+2) A1 = SD(JPTR+1) A2 = SD(JPTR+2) C DO 4 I = NSAMPS, 1, -1 C OUTPUT = B0*DATA(I) + B1*X1 + B2*X2 OUTPUT = OUTPUT - ( A1*Y1 + A2*Y2 ) Y2 = Y1 Y1 = OUTPUT X2 = X1 X1 = DATA(I) DATA(I) = OUTPUT C 4 CONTINUE C JPTR = JPTR + 3 C 3 CONTINUE C END IF C RETURN END c subroutine zero(x,lx) c dimension x(1) do 12 i=1,lx 12 x(i)=0.0 return end c subroutine decimate(n,ndec,datain,dt,samdec,nd) dimension datain(1) c c D e c i m a t i o n c samdec=dt*float(ndec) nd=(n-1)/ndec+1 temp=0.0 do 4200 i=1,nd temp=datain(ndec*(i-1)+1) datain(i)=temp 4200 continue return end c subroutine rmavg(rdata,npts) c dimension rdata(1) c c ... average removal c sum = 0.0 c do 22 j=1,npts 22 sum = sum + rdata(j) c avg = sum / float(npts) do 23 k=1,npts 23 rdata(k) = rdata(k) - avg 99 return end