C PROGRAM WINDOW.FOR C Modified in 1994 by S.A. Tretter to add bandpass Hilbert transform C Modified 10/4/95 by S.A. Tretter to correctly handle the case of C less than 3 taps C Modified 4/23/98 by S.A. Tretter to allow full-band Hilbert C transforms, i.e., FL = 0 and FH = fs/2. C Modified 11/12/99 by S.A. Tretter to include bandpass differentiator C----------------------------------------------------------------------- C MAIN PROGRAM: WINDOW DESIGN OF LINEAR PHASE, LOWPASS, HIGHPASS C BANDPASS, AND BANDSTOP FIR DIGITAL FILTERS C AUTHOR: LAWRENCE R. RABINER AND CAROL A. MCGONEGAL C BELL LABORATORIES, MURRAY HILL, NEW JERSEY, 07974 C MODIFIED JAN. 1978 BY DOUG PAUL, MIT LINCOLN LABORATORIES C TO INCLUDE SUBROUTINES FOR OBTAINING FILTER BAND EDGES AND RIPPLES C C INPUT: NF IS THE FILTER LENGTH IN SAMPLES C 3 <= NF <= 1024 C C ITYPE IS THE WINDOW TYPE C ITYPE = 1 RECTANGULAR WINDOW C ITYPE = 2 TRIANGULAR WINDOW C ITYPE = 3 HAMMING WINDOW C ITYPE = 4 GENERALIZED HAMMING WINDOW C ITYPE = 5 HANNING WINDOW C ITYPE = 6 KAISER (I0-SINH) WINDOW C ITYPE = 7 CHEBYSHEV WINDOW C C JTYPE IS THE FILTER TYPE C JTYPE = 1 LOWPASS FILTER C JTYPE = 2 HIGHPASS FILTER C JTYPE = 3 BANDPASS FILTER C JTYPE = 4 BANDSTOP FILTER C JTYPE = 5 BANDPASS HILBERT TRANSFORM C H(w) = - j sign w for FL < |f/fs| < FH C 0 elsewhere C h(n) = 2 ((sin n pi B)/ (n pi)) sin 2 pi n Fo C where B = (FH_actual - FL_actual)/fs, C Fo = (FH_actual + FL_actual)/(2 fs) C JTYPE = 6 BANDPASS DIFFERENTIATOR C C H(w) = jw for FL < |w/ws| < FH C 0 elsewhere C C h(n) = 2(fs/n)[ FH cos(2 pi n FH) - FL cos(2 pi n FL)] C -2(fs/n**2 pi) sin(pi n B) cos(2 pi n Fo) C C FC IS THE NORMALIZED CUTOFF FREQUENCY C 0 <= FC <= 0.5 C FL AND FH ARE THE NORMALIZED FILTER CUTOFF FREQUENCIES C 0 <= FL <= FH <= 0.5 C IWP OPTIONALLY PRINTS OUT THE WINDOW VALUES C IWP = 0 DO NOT PRINT C IWP = 1 PRINT C IMD REQUESTS ADDITIONAL RUNS C IMD = 1 NEW RUN C IMD = 0 TERMINATES PROGRAM C----------------------------------------------------------------------- C DIMENSION W(512), G(512),H(1025) INTEGER OTCD1, OTCD2 CHARACTER BELL, COEFFILE*72,KPLT COMPLEX RESP C BELL = CHAR(7) PI = 4.0*ATAN(1.0) TWOPI = 2.0*PI C C DEFINE I/O DEVICE CODES C INPUT: INPUT TO THIS PROGRAM IS USER-INTERACTIVE C THAT IS - A QUESTION IS WRITTEN ON THE USER C TERMINAL (OTCD1) AND THE USER TYPES IN THE ANSWER. C C INCOD = 5 OTCD1 = 8 OTCD2 = 7 C C INPUT THE FILTER LENGTH(NF), WINDOW TYPE(ITYPE) AND FILTER TYPE(JTYPE) C 5540 WRITE(*,'(A\)') ' ENTER NAME OF LISTING FILE: ' READ(*,'(A)') COEFFILE OPEN(8,FILE = COEFFILE) 5550 WRITE(*,'(A,\)') ' ENTER FILENAME FOR COEFFICIENTS: ' READ(*,'(A)') COEFFILE OPEN(7,FILE=COEFFILE) WRITE(*,'(A\)') ' ENTER SAMPLING FREQUENCY IN HZ: ' READ(*,*) FSAMP WRITE(8,*) ' SAMPLING FREQUENCY = ',FSAMP WRITE(*,9998) 9998 FORMAT( 1' WINDOW TYPES'/ 2' 1 RECTANGULAR WINDOW'/ 3' 2 TRIANGULAR WINDOW'/ 4' 3 HAMMING WINDOW 0.54+0.46 cos(theta)'/ 5' 4 GENERALIZED HAMMING WINDOW alpha + (1-alpha *) cos(theta)'/ 6' 5 HANNING WINDOW 0.5+0.5 cos(theta)'/ 7' 6 KAISER (I0-SINH) WINDOW'/ 8' 7 CHEBYSHEV WINDOW'/ 9' '/ A' FILTER TYPES '/ B' 1 LOWPASS FILTER '/ C' 2 HIGHPASS FILTER '/ D' 3 BANDPASS FILTER '/ E' 4 BANDSTOP FILTER '/ F' 5 BANDPASS HILBERT TRANSFORM '/ G' 6 BANDPASS DIFFERENTIATOR '/ H' ' ) 10 WRITE (*,9999) 9999 FORMAT (' ENTER FILTER LENGTH, WINDOW TYPE, FILTER TYPE: '\) READ (INCOD,*) NF, ITYPE, JTYPE IF (NF.GT.1024) THEN WRITE(*,*) ' ' WRITE (*,9997) NF 9997 FORMAT(4H NF=,I4,' IS TOO LARGE. MUST BE NO GREATER THAN 1024.') WRITE(*,*) BELL GO TO 10 ENDIF IF (NF.LT.3) THEN WRITE(*,*) ' ' WRITE(*,6667) NF 6667 FORMAT(' NF=',I4,' IS TOO SMALL. MUST BE NO LESS THAN 3.') WRITE(*,*) BELL GOTO 10 ENDIF C C N IS HALF THE LENGTH OF THE SYMMETRIC FILTER C N = (NF+1)/2 IF (JTYPE.NE.1 .AND. JTYPE.NE.2) GO TO 50 C C FOR THE IDEAL LOWPASS OR HIGHPASS DESIGN - INPUT FC C 40 WRITE (*,9996) 9996 FORMAT(' ENTER IDEAL CUTOFF FREQUENCY IN HZ: '\) READ (INCOD,*) FCACT FC = FCACT/FSAMP IF (FC.GT.0.0 .AND. FC.LT.0.5) GO TO 60 WRITE(*,*) BELL WRITE (*,9995) FCACT 9995 FORMAT(' FC=',F14.7,'GREATER THAN FS/2, REENTER DATA.') GO TO 40 C C FOR BANDPASS, BANDSTOP, HILBERT, OR DIFFERENTIATOR - INPUT FL AND FH C 50 WRITE (*,9994) 9994 FORMAT (' SPECIFY LOWER, UPPER CUTOFF IN HZ: '\) READ (INCOD,*) FLACT, FHACT FL = FLACT/FSAMP FH = FHACT/FSAMP IF (FL.GE.0.0 .AND. FL.LT.0.5 .AND. FH.GT.0.0 .AND. FH.LE.0.5 * .AND. FH.GT.FL) GO TO 62 WRITE(*,*) BELL IF (FL.LT.0. .OR. FL.GT.0.5) WRITE (*,8885) FLACT 8885 FORMAT(' FL=',F14.7,' < 0 OR > Fs/2, REENTER DATA') IF (FH.LT.0. .OR. FH.GT.0.5) WRITE (*,8886) FHACT 8886 FORMAT(' FH=',F14.7,' < 0 OR > Fs/2, REENTER DATA') IF (FH.LT.FL) WRITE (*,9992) FHACT, FLACT 9992 FORMAT (4H FH=, F14.7, 20H IS SMALLER THAN FL=, F14.7, 8H REENTER, * 5H DATA) GO TO 50 C INPUT SLOPE FOR DIFFERENTIATOR 62 IF(JTYPE.EQ.6) THEN WRITE(*,*) ' ' WRITE(*,*) 'The amplitude response for a differentiator will' WRITE(*,*) 'approximate A(f) = SLOPE*f .' WRITE(*,1000) 1000 FORMAT(' ENTER SLOPE: '\) READ(INCOD,*) SLOPE ENDIF 60 IF (ITYPE.NE.7) GO TO 70 C C INPUT FOR CHEBYSHEV WINDOW--2 OF THE 3 PARAMETERS NF, DPLOG, AND DF C MUST BE SPECIFIED, WHERE DPLOG IS THE DESIRED FILTER RIPPLE(DB SCALE), C DF IS THE TRANSITION WIDTH (NORMALIZED) OF THE FILTER, C AND NF IS THE FILTER LENGTH. THE UNSPECIFIED PARAMETER C IS READ IN WITH THE ZERO VALUE. C WRITE (*,9991) 9991 FORMAT('0SPECIFY CHEBYSHEV SIDELOBE PEAK RIPPLE ATTENUATION RELATI *VE TO MAIN LOBE PEAK'/' IN +DB (DPLOG), AND TRANSITION WIDTH IN HZ * (DF). ONE OF THE THREE PARAMETERS'/' (NF, DPLOG, DF) SHOULD BE 0: * ',\) READ (INCOD,*) DPLOG, DF DF = DF/FSAMP DP = 10.0**(-DPLOG/20.0) CALL CHEBC(NF, DP, DF, N, X0, XN) C C IEO IS AN EVEN, ODD INDICATOR, IEO = 0 FOR EVEN, IEO = 1 FOR ODD C 70 IEO = MOD(NF,2) IF (IEO.EQ.1 .OR. JTYPE.EQ.1 .OR. JTYPE.EQ.3) GO TO 80 WRITE (*,9990) 9990 FORMAT(' NF MUST BE ODD FOR HP, BS, HILBERT, or DIFFERENTIATOR.'/ * ' NF IS BEING INCREASED BY 1') NF = NF + 1 N = (1+NF)/2 IEO = 1 80 CONTINUE C C COMPUTE IDEAL (UNWINDOWED) IMPULSE RESPONSE FOR FILTER C C1 = FC IF (JTYPE.EQ.3 .OR. JTYPE.EQ.4 .OR. JTYPE.EQ.5) C1 = FH - FL IF (IEO.EQ.1) G(1) = 2.*C1 IF (JTYPE.EQ.5 .OR. JTYPE.EQ.6) G(1) = 0. I1 = IEO + 1 DO 90 I=I1,N XN = I - 1 IF (IEO.EQ.0) XN = XN + 0.5 C = PI*XN C3 = C*C1 IF (JTYPE.EQ.1 .OR. JTYPE.EQ.2) C3 = 2.*C3 G(I) = SIN(C3)/C IF (JTYPE.EQ.3 .OR. JTYPE.EQ.4) G(I) = G(I)*2.*COS(C*(FL+FH)) IF (JTYPE.EQ.5) G(I) = G(I)*2.*SIN(C*(FL+FH)) F0 = (FH + FL)/2. B = FH - FL IF(JTYPE.EQ.6) THEN G(I) = 2.*(FSAMP/XN)*(FH*COS(2*PI*XN*FH) - * FL*COS(2.*PI*XN*FL)) - 2.*(FSAMP/(PI*XN**2))*SIN(PI*XN*B)* * COS(2.*PI*XN*F0) G(I) = G(I)*SLOPE*0.5/PI ENDIF 90 CONTINUE C C COMPUTE A RECTANGULAR WINDOW C IF (ITYPE.EQ.1) WRITE (OTCD1,9989) NF 9989 FORMAT (' RECTANGULAR WINDOW -- TAPS = ', I4) DO 100 I=1,N W(I) = 1. 100 CONTINUE C C DISPATCH ON WINDOW TYPE C GO TO (200, 110, 120, 140, 150, 160, 170), ITYPE C C TRIANGULAR WINDOW C 110 CALL TRIANG(W, N, IEO) WRITE (OTCD1,9988) NF 9988 FORMAT (' TRIANGULAR WINDOW -- TAPS = ', I4) GO TO 180 C C HAMMING WINDOW C 120 ALPHA = 0.54 WRITE (OTCD1,9987) NF 9987 FORMAT (' HAMMING WINDOW -- TAPS = ', I4) 130 BETA = 1. - ALPHA CALL HAMMIN(NF, W, N, IEO, ALPHA, BETA) GO TO 180 C C GENERALIZED HAMMING WINDOW C FORM OF WINDOW IS W(M)=ALPHA+BETA*COS((TWOPI*M)/(NF-1)) C BETA IS AUTOMATICALLY SET TO 1.-ALPHA C READ IN ALPHA C 140 WRITE (*,'(A\)')' SPECIFY ALPHA FOR GENERALIZED HAMMING WINDOW: ' READ (*,*) ALPHA WRITE (OTCD1,9984) NF,ALPHA 9984 FORMAT (' GENERALIZED HAMMING WINDOW -- TAPS = ', I4/' alpha = ' * ,F10.4) GO TO 130 C C HANNING WINDOW C 150 ALPHA = 0.5 WRITE (OTCD1,9983) NF 9983 FORMAT (' HANNING WINDOW -- TAPS = ', I4) C C INCREASE NF BY 2 AND N BY 1 FOR HANNING WINDOW SO ZERO C ENDPOINTS ARE NOT PART OF WINDOW C NF = NF + 2 N = N + 1 GO TO 130 C C KAISER (I0-SINH) WINDOW C NEED TO SPECIFY PARAMETER ATT=STOPBAND ATTENUATION IN DB C 160 WRITE (*,'(A\)') ' SPECIFY OUT-OF-BAND ATTENUATION IN +DB FOR KAIS *ER WINDOW: ' READ (*,*) ATT IF (ATT.GT.50.) BETA = 0.1102*(ATT-8.7) IF (ATT.GE.20.96 .AND. ATT.LE.50.) BETA = 0.58417*(ATT-20.96)** * 0.4 + 0.07886*(ATT-20.96) IF (ATT.LT.20.96) BETA = 0. CALL KAISER(NF, W, N, IEO, BETA) WRITE (OTCD1,9981) NF 9981 FORMAT (' KAISER WINDOW -- TAPS = ', I4) WRITE (OTCD1,9980) ATT, BETA 9980 FORMAT (6H ATT=, F14.7, 7H BETA=, F14.7) GO TO 180 C C CHEBYSHEV WINDOW C 170 CALL CHEBY(NF, W, N, IEO, DP, X0, XN) WRITE (OTCD1,9979) NF 9979 FORMAT (' CHEBYSHEV WINDOW -- TAPS = ', I4) WRITE (OTCD1,9978) -20.*LOG10(DP), DF*FSAMP 9978 FORMAT (4H DP=, F14.7, 5H DF=, F14.7) C C WINDOW IDEAL FILTER RESPONSE C CHANGE BACK NF AND N FOR HANNING WINDOW C 180 IF (ITYPE.EQ.5) NF = NF - 2 IF (ITYPE.EQ.5) N = N - 1 DO 190 I=1,N G(I) = G(I)*W(I) 190 CONTINUE C C PRINT OUT RESULTS C 200 CONTINUE C WRITE (OTCD1,9977) C9977 FORMAT (36H PRINT OUT WINDOW VALUES(1=YES,0=NO)) C READ (INCOD,9976) IWP C9976 FORMAT (I1) IWP = 0 IF (IWP.EQ.0) GO TO 220 WRITE (OTCD1,9975) 9975 FORMAT (14H WINDOW VALUES) DO 210 I=1,N J = N + 1 - I K = NF + 1 - I WRITE (OTCD1,9974) I, W(J), K 9974 FORMAT (10X, 3H W(, I3, 2H)=, E15.8, 4H =W(, I4, 1H)) 210 CONTINUE 220 IF (JTYPE.EQ.1) WRITE (OTCD1,9973) 9973 FORMAT (26H **LOWPASS FILTER DESIGN**) IF (JTYPE.EQ.2) WRITE (OTCD1,9972) 9972 FORMAT (27H **HIGHPASS FILTER DESIGN**) IF (JTYPE.EQ.3) WRITE (OTCD1,9971) 9971 FORMAT (27H **BANDPASS FILTER DESIGN**) IF (JTYPE.EQ.4) WRITE (OTCD1,9970) 9970 FORMAT (27H **BANDSTOP FILTER DESIGN**) IF (JTYPE.EQ.5) WRITE (OTCD1,8887) 8887 FORMAT(' **BANDPASS HILBERT DESIGN**') IF (JTYPE.EQ.6) WRITE (OTCD1,8888) 8888 FORMAT(' **BANDPASS DIFFERENTIATOR**') IF(JTYPE.EQ.1) WRITE(OTCD1,*)' IDEAL LOWPASS CUTOFF = ',FCACT IF(JTYPE.EQ.2) WRITE(OTCD1,*)' IDEAL HIGHPASS CUTOFF = ',FCACT IF(JTYPE.EQ.3 .OR. JTYPE.EQ.4 .OR. JTYPE.EQ.5 .OR. JTYPE.EQ.6) * WRITE (OTCD1,*) ' IDEAL CUTOFF FREQUENCIES = ',FLACT, FHACT IF (JTYPE.EQ.1 .OR. JTYPE.EQ.3 .OR. JTYPE.EQ.5 .OR. JTYPE.EQ.6) * GO TO 240 DO 230 I=2,N G(I) = -G(I) 230 CONTINUE G(1) = 1.0 - G(1) C C WRITE OUT IMPULSE RESPONSE C 240 DO 250 I=1,N J = N + 1 - I K = NF + 1 - I H(I) = G(J) IF(JTYPE.EQ.5 .OR. JTYPE.EQ.6) H(I) = -H(I) H(K) = G(J) IF(JTYPE.NE.5.AND. JTYPE.NE.6) WRITE(OTCD1,9966)I-1, G(J), K-1 9966 FORMAT (10X, 3H H(, I3, 2H)=, E15.8, 4H =H(, I4, 1H)) IF(JTYPE.EQ.5 .OR. JTYPE.EQ.6) WRITE (OTCD1,8866) * I-1,G(J),K-1 8866 FORMAT(10X,4H -H(,I3, 2H)=, E15.8, 4H =H(, I4, 1H)) 250 CONTINUE WRITE(7,*) NF DO 255 I=1,NF 255 WRITE(7,*) H(I) C WRITE(*,*) BELL WRITE(*,'(A\)') ' CREATE (FREQ,RESPONSE) FILE (Y OR N)? ' READ(*,'(A)') KPLT IF(KPLT.NE.'Y'.AND.KPLT.NE.'y') STOP 5530 WRITE(*,'(A,\)') ' ENTER FILENAME: ' READ(*,'(A)') COEFFILE OPEN(9,FILE=COEFFILE) WRITE(*,'(A\)') ' LINEAR (L) OR DB (D) SCALE ?: ' READ(*,'(A)') KPLT DELF = 0.5/(5.*NF) DO 260 FNORM = 0, 0.5,DELF RESP = 0. DO 261 K = 0,NF-1 ANG = K*TWOPI*FNORM 261 RESP = RESP + H(K+1)*CMPLX(COS(ANG),-SIN(ANG)) IF(KPLT.EQ.'L' .OR. KPLT.EQ.'l') THEN AMP = CABS(RESP) ELSEIF(KPLT.EQ.'d' .OR. KPLT.EQ.'D') THEN AMP = CABS(RESP) AMP = MAX(1.E-20,AMP) AMP = 20.*ALOG10(AMP) ENDIF WRITE(9,*) FNORM*FSAMP,AMP 260 CONTINUE C=============================================================== C CALL FLCHAR(NF,ITYPE,JTYPE,FC,FL,FH,N,IEO,G,OTCD1,FSAMP) C WRITE (OTCD1,9965) C9965 FORMAT (1H /1H /1H /1H ) C WRITE (OTCD1,9964) C9964 FORMAT (1H1) C WRITE (OTCD1,9963) C9963 FORMAT (26H MORE DATA(1=YES,0=NO)(I1)) C READ (INCOD,9962) IMD C9962 FORMAT (I1) C IF (IMD.EQ.1) GO TO 10 END C C----------------------------------------------------------------------- C SUBROUTINE: TRIANG C TRIANGULAR WINDOW C----------------------------------------------------------------------- C SUBROUTINE TRIANG( W, N, IEO) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW COEFFICIENTS FOR HALF THE WINDOW C N = HALF WINDOW LENGTH=(NF+1)/2 C IEO = EVEN - ODD INDICATION--IEO=0 FOR NF EVEN C DIMENSION W(1) FN = N DO 10 I=1,N XI = I - 1 IF (IEO.EQ.0) XI = XI + 0.5 W(I) = 1. - XI/FN 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: HAMMIN C GENERALIZED HAMMING WINDOW ROUTINE C WINDOW IS W(N) = ALPHA + BETA * COS( TWOPI*(N-1) / (NF-1) ) C----------------------------------------------------------------------- C SUBROUTINE HAMMIN(NF, W, N, IEO, ALPHA, BETA) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = HALF LENGTH OF FILTER=(NF+1)/2 C IEO = EVEN ODD INDICATOR--IEO=0 IF NF EVEN C ALPHA = CONSTANT OF WINDOW C BETA = CONSTANT OF WINDOW--GENERALLY BETA=1-ALPHA C DIMENSION W(1) PI2 = 8.0*ATAN(1.0) FN = NF - 1 DO 10 I=1,N FI = I - 1 IF (IEO.EQ.0) FI = FI + 0.5 W(I) = ALPHA + BETA*COS((PI2*FI)/FN) 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: KAISER C KAISER WINDOW C----------------------------------------------------------------------- C SUBROUTINE KAISER(NF, W, N, IEO, BETA) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = FILTER HALF LENGTH=(NF+1)/2 C IEO = EVEN ODD INDICATOR--IEO=0 IF NF EVEN C BETA = PARAMETER OF KAISER WINDOW C DIMENSION W(1) REAL INO BES = INO(BETA) XIND = FLOAT(NF-1)*FLOAT(NF-1) DO 10 I=1,N XI = I - 1 IF (IEO.EQ.0) XI = XI + 0.5 XI = 4.*XI*XI W(I) = INO(BETA*SQRT(1.-XI/XIND)) W(I) = W(I)/BES 10 CONTINUE RETURN END C C----------------------------------------------------------------------- C FUNCTION: INO C BESSEL FUNCTION FOR KAISER WINDOW C----------------------------------------------------------------------- C REAL FUNCTION INO(X) Y = X/2. T = 1.E-08 E = 1. DE = 1. DO 10 I=1,25 XI = I DE = DE*Y/XI SDE = DE*DE E = E + SDE IF (E*T-SDE) 10, 10, 20 10 CONTINUE 20 INO = E RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: CHEBC C SUBROUTINE TO GENERATE CHEBYSHEV WINDOW PARAMETERS WHEN C ONE OF THE THREE PARAMETERS NF,DP AND DF IS UNSPECIFIED C----------------------------------------------------------------------- C SUBROUTINE CHEBC(NF, DP, DF, N, X0, XN) C C NF = FILTER LENGTH (IN SAMPLES) C DP = FILTER RIPPLE (ABSOLUTE SCALE) C DF = NORMALIZED TRANSITION WIDTH OF FILTER C N = (NF+1)/2 = FILTER HALF LENGTH C X0 = (3-C0)/(1+C0) WITH C0=COS(PI*DF) = CHEBYSHEV WINDOW CONSTANT C XN = NF-1 C PI = 4.*ATAN(1.0) IF (NF.NE.0) GO TO 10 C C DP,DF SPECIFIED, DETERMINE NF C C1 = COSHIN((1.+DP)/DP) C0 = COS(PI*DF) X = 1. + C1/COSHIN(1./C0) C C INCREMENT BY 1 TO GIVE NF WHICH MEETS OR EXCEEDS SPECS ON DP AND DF C NF = X + 1.0 N = (NF+1)/2 XN = NF - 1 GO TO 30 10 IF (DF.NE.0.0) GO TO 20 C C NF,DP SPECIFIED, DETERMINE DF C XN = NF - 1 C1 = COSHIN((1.+DP)/DP) C2 = COSH(C1/XN) DF = ARCCOS(1./C2)/PI GO TO 30 C C NF,DF SPECIFIED, DETERMINE DP C 20 XN = NF - 1 C0 = COS(PI*DF) C1 = XN*COSHIN(1./C0) DP = 1./(COSH(C1)-1.) 30 X0 = (3.-COS(2.*PI*DF))/(1.+COS(2.*PI*DF)) RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: CHEBY C DOLPH CHEBYSHEV WINDOW DESIGN C----------------------------------------------------------------------- C SUBROUTINE CHEBY(NF, W, N, IEO, DP,X0, XN) C C NF = FILTER LENGTH IN SAMPLES C W = WINDOW ARRAY OF SIZE N C N = HALF LENGTH OF FILTER = (NF+1)/2 C IEO = EVEN-ODD INDICATOR--IEO=0 FOR NF EVEN C DP = WINDOW RIPPLE ON AN ABSOLUTE SCALE C DF = NORMALIZED TRANSITION WIDTH OF WINDOW C X0 = WINDOW PARAMETER RELATED TO TRANSITION WIDTH C XN = NF-1 C DIMENSION W(1) DIMENSION PR(1024), PI(1024) PIE = 4.*ATAN(1.0) XN = NF - 1 FNF = NF ALPHA = (X0+1.)/2. BETA = (X0-1.)/2. TWOPI = 2.*PIE C2 = XN/2. DO 40 I=1,NF XI = I - 1 F = XI/FNF X = ALPHA*COS(TWOPI*F) + BETA IF (ABS(X)-1.) 10, 10, 20 10 P = DP*COS(C2*ARCCOS(X)) GO TO 30 20 P = DP*COSH(C2*COSHIN(X)) 30 PI(I) = 0. PR(I) = P C C FOR EVEN LENGTH FILTERS USE A ONE-HALF SAMPLE DELAY C ALSO THE FREQUENCY RESPONSE IS ANTISYMMETRIC IN FREQUENCY C IF (IEO.EQ.1) GO TO 40 PR(I) = P*COS(PIE*F) PI(I) = -P*SIN(PIE*F) IF (I.GT.(NF/2+1)) PR(I) = -PR(I) IF (I.GT.(NF/2+1)) PI(I) = -PI(I) 40 CONTINUE C C USE DFT TO GIVE WINDOW C TWN = TWOPI/FNF DO 60 I=1,N XI = I - 1 SUM = 0. DO 50 J=1,NF XJ = J - 1 SUM = SUM + PR(J)*COS(TWN*XJ*XI) + PI(J)*SIN(TWN*XJ*XI) 50 CONTINUE W(I) = SUM 60 CONTINUE C1 = W(1) DO 70 I=1,N W(I) = W(I)/C1 70 CONTINUE RETURN END C C----------------------------------------------------------------------- C FUNCTION: COSHIN C FUNCTION FOR HYPERBOLIC INVERSE COSINE OF X C----------------------------------------------------------------------- C REAL FUNCTION COSHIN(X) COSHIN = ALOG(X+SQRT(X*X-1.)) RETURN END C C----------------------------------------------------------------------- C FUNCTION: ARCCOS C FUNCTION FOR INVERSE COSINE OF X C----------------------------------------------------------------------- C FUNCTION ARCCOS(X) IF (X) 30, 20, 10 10 A = SQRT(1.-X*X)/X ARCCOS = ATAN(A) RETURN 20 ARCCOS = 2.*ATAN(1.0) RETURN 30 A = SQRT(1.-X*X)/X ARCCOS = ATAN(A) + 4.*ATAN(1.0) RETURN END C C----------------------------------------------------------------------- C FUNCTION: COSH C FUNCTION FOR HYPERBOLIC COSINE OF X C----------------------------------------------------------------------- C REAL FUNCTION COSH(X) COSH = (EXP(X)+EXP(-X))/2. RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: FLCHAR C SUBROUTINE TO DETERMINE FILTER CHARACTERISTICS C----------------------------------------------------------------------- C SUBROUTINE FLCHAR(NF,ITYPE,JTYPE,FC,FL,FH,N,IEO,G,OTCD1,FSAMP) C C NF = FILTER LENGTH IN SAMPLES C ITYPE = WINDOW TYPE C JTYPE = FILTER TYPE C FC = IDEAL CUTOFF OF LP OR HP FILTER C FL = LOWER CUTOFF OF BP OR BS FILTER C FH = UPPER CUTOFF OF BP OR BS FILTER C N = FILTER HALF LENGTH = (NF+1) / 2 C IEO = EVEN ODD INDICATOR C G = FILTER ARRAY OF SIZE N C OTCD1 = OUTPUT CODE FOR LINE PRINTER USED IN WRITE STATEMENTS C DIMENSION G(1) DIMENSION RESP(2048) INTEGER OTCD1 C C NOT FOR FOR TRIANGULAR WINDOW C IF (ITYPE.EQ.2) RETURN C C DFT TO GET FREQ RESP C PI = 4.*ATAN(1.0) C C UP TO 4096 PT DFT C NR = 8*NF IF (NR.GT.2048) NR = 2048 XNR = NR TWN = PI/XNR SUMI = -G(1)/2. IF (IEO.EQ.0) SUMI = 0. DO 20 I=1,NR XI = I - 1 TWNI = TWN*XI SUM = SUMI DO 10 J=1,N XJ = J - 1 IF (IEO.EQ.0) XJ = XJ + .5 SUM = SUM + G(J)*COS(XJ*TWNI) 10 CONTINUE RESP(I) = 2.*SUM 20 CONTINUE C C DISPATCH ON FILTER TYPE C GO TO (30, 40, 50, 60), JTYPE C C LOWPASS C 30 CALL RIPPLE(NR, 1., 0., FC, RESP, F1, F2, DB) WRITE (OTCD1,9999) F2*FSAMP, DB 9999 FORMAT (17H PASSBAND CUTOFF , F20.4, 9H RIPPLE , F8.3, 3H DB) CALL RIPPLE(NR, 0., FC, .5, RESP, F1, F2, DB) WRITE (OTCD1,9998) F1*FSAMP, DB 9998 FORMAT (17H STOPBAND CUTOFF , F20.4, 9H RIPPLE , F8.3, 3H DB) RETURN C C HIGHPASS C 40 CALL RIPPLE(NR, 0., 0., FC, RESP, F1, F2, DB) WRITE (OTCD1,9998) F2*FSAMP, DB CALL RIPPLE(NR, 1., FC, .5, RESP, F1, F2, DB) WRITE (OTCD1,9999) F1*FSAMP, DB RETURN C C BANDPASS C 50 CALL RIPPLE(NR, 0., 0., FL, RESP, F1, F2, DB) WRITE (OTCD1,9998) F2*FSAMP, DB CALL RIPPLE(NR, 1., FL, FH, RESP, F1, F2, DB) WRITE (OTCD1,9997) F1*FSAMP, F2*FSAMP, DB 9997 FORMAT (18H PASSBAND CUTOFFS , F6.4, 2X, F6.4, 8H RIPPLE, F9.3, * 3H DB) CALL RIPPLE(NR, 0., FH, .5, RESP, F1, F2, DB) WRITE (OTCD1,9998) F1*FSAMP, DB RETURN C C STOPBAND C 60 CALL RIPPLE(NR, 1., 0., FL, RESP, F1, F2, DB) WRITE (OTCD1,9999) F2*FSAMP, DB CALL RIPPLE(NR, 0., FL, FH, RESP, F1, F2, DB) WRITE (OTCD1,9996) F1*FSAMP, F2*FSAMP, DB 9996 FORMAT (18H STOPBAND CUTOFFS , F6.4, 2X, F6.4, 8H RIPPLE, F9.3, * 3H DB) CALL RIPPLE(NR, 1., FH, .5, RESP, F1, F2, DB) WRITE (OTCD1,9999) F1*FSAMP, DB RETURN END C C----------------------------------------------------------------------- C SUBROUTINE: RIPPLE C FINDS LARGEST RIPPLE IN BAND AND LOCATES BAND EDGES BASED ON THE C POINT WHERE THE TRANSITION REGION CROSSES THE MEASURED RIPPLE BOUND C----------------------------------------------------------------------- C SUBROUTINE RIPPLE(NR, RIDEAL, FLOW, FHI, RESP, F1, F2, DB) C C NR = SIZE OF RESP C RIDEAL = IDEAL FREQUENCY RESPONSE C FLOW = LOW EDGE OF IDEAL BAND C FHI = HIGH EDGE OF IDEAL BAND C RESP = FREQUENCY RESPONSE OF SIZE NR C F1 = COMPUTED LOWER BAND EDGE C F2 = COMPUTED UPPER BAND EDGE C DB = DEVIATION FROM IDEAL RESPONSE IN DB C DIMENSION RESP(1) XNR = NR C C BAND LIMITS C IFLOW = 2.*XNR*FLOW + 1.5 IFHI = 2.*XNR*FHI + 1.5 IF (IFLOW.EQ.0) IFLOW = 1 IF (IFHI.GE.NR) IFHI = NR - 1 C C FIND MAX AND MIN PEAKS IN BAND C RMIN = RIDEAL RMAX = RIDEAL DO 20 I=IFLOW,IFHI IF(I.EQ.1) THEN IF(RESP(I).LE.RMAX.OR.RESP(I).LT.RESP(I+1)) GOTO 10 ELSE IF (RESP(I).LE.RMAX .OR. RESP(I).LT.RESP(I-1) .OR. * RESP(I).LT.RESP(I+1)) GO TO 10 RMAX = RESP(I) ENDIF 10 IF(I.EQ.1) THEN IF(RESP(I).GE.RMIN.OR.RESP(I).GT.RESP(I+1))GOTO 20 ELSE IF (RESP(I).GE.RMIN .OR. RESP(I).GT.RESP(I-1) .OR. * RESP(I).GT.RESP(I+1)) GO TO 20 RMIN = RESP(I) ENDIF 20 CONTINUE C C PEAK DEVIATION FROM IDEAL C RIPL = AMAX1(RMAX-RIDEAL,RIDEAL-RMIN) C C SEARCH FOR LOWER BAND EDGE C F1 = FLOW IF (FLOW.EQ.0.0) GO TO 50 DO 30 I=IFLOW,IFHI IF (ABS(RESP(I)-RIDEAL).LE.RIPL) GO TO 40 30 CONTINUE 40 XI = I - 1 C C LINEAR INTERPOLATION OF BAND EDGE FREQUENCY TO IMPROVE ACCURACY C X1 = .5*XI/XNR X0 = .5*(XI-1.)/XNR Y1 = ABS(RESP(I)-RIDEAL) Y0 = ABS(RESP(I-1)-RIDEAL) IF(Y1.NE.Y0) F1 = (X1-X0)/(Y1-Y0)*(RIPL-Y0) + X0 C C SEARCH FOR UPPER BAND EDGE C 50 F2 = FHI IF (FHI.EQ.0.5) GO TO 80 DO 60 I=IFLOW,IFHI J = IFHI + IFLOW - I IF (ABS(RESP(J)-RIDEAL).LE.RIPL) GO TO 70 60 CONTINUE 70 XI = J - 1 C C LINEAR INTERPOLATION OF BAND EDGE FREQUENCY TO IMPROVE ACCURACY C X1 = .5*XI/XNR X0 = .5*(XI+1.)/XNR Y1 = ABS(RESP(J)-RIDEAL) Y0 = ABS(RESP(J+1)-RIDEAL) IF(Y1.NE.Y0) F2 = (X1-X0)/(Y1-Y0)*(RIPL-Y0) + X0 C C DEVIATION FROM IDEAL IN DB C 80 DB = 20.*ALOG10(RIPL+RIDEAL) RETURN END