CTITLES2CSRJ -- DEVELOP BANDPASS FILTER IN TIME DOMAIN 00000100 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA AUTHOR HANS HOOGSTRAAT 00000200 CA DESIGNER HANS HOOGSTRAAT 00000300 CA LANGUAGE FORTRAN 00000400 CA SYSTEM IBM AND CRAY 00000400 CA WRITTEN 00000500 C REVISED 05-26-77 BY R. MCMILLAN. CHANGE NAME TO S2CSRJ FROM 00000600 C S1CSRJ AND ADD IPR TO PARAMETER LIST. 00000700 C REVISED 10-17-77 MENDEKE , NYMAN MADE APPROXIMATION OF 00000900 C XN1, XN2 FOR D1, D2 VERY SMALL. 00001000 C REVISED 07-29-83 PARKER, NYMAN INCREASE APPROX. OF XN1, XN2. 00001100 C REVISED 08-01-83 PARKER, NYMAN FORCE NPTS TO BE ODDSET 00001300 C D1,D2 = MINVAL IF THEY ARE 0. 00001400 C REVISED 12-10-84 KNIGHT; FOR THE CRAY. 00001410 C REVISED 10-02-85 RDK; BROUGHT CRAY/IBM VERSIONS INTO ACCORD. 00001410 C REVISED 11-13-89 RDK; REMOVE EXTERNAL FOR S1ATP. 00001410 C 00001500 C 00001600 CA 00001700 CA 00001800 CA CALL S2CSRJ (FILT,NPTS,SAMPR,LOWCUT,LOWPAS,HIGPAS,HIGCUT,IPR) 00001900 CA OUTPUT FILT = RETURNED BANDPASS FILTER COEFFIEIENTS R4 00002000 CA IN/OUT NPTS = FILTER LENGTH IN SAMPLES (MAY BE I4 00002100 CA REDUCED BY THIS ROUTINE. SHOULD BE ODD. 00002200 CA IF EVEN, THIS ROUTINE WILL MAKE IT ODD.) 00002300 CA INPUT SAMPR = FILTER SAMPLE RATE IN MS I4 00002400 CA INPUT LOWCUT = LOW CUT FILTER SETTING IN HZ I4 00002500 CA INPUT LOWPAS = LOW PASS FILTER SETTING IN HZ I4 00002600 CA INPUT HIGPAS = HIGH PASS FILTER SETTING IN HZ I4 00002700 CA INPUT HIGCUT = HIGH CUT FILTER SETTING IN HZ I4 00002800 CA INPUT IPR = PRINTER UNIT NUMBER I4 00002900 CA 00003000 CA 00003100 CA S2CSRJ DEVELOPS A COSINE REJECTION BANDPASS FILTER. 00003200 C 00003300 C ====================================================== 00003400 C * PROGRAM DEVELOPED BY 00003500 C * H * HOOGSTRAAT PROGRAMMING SERVICES LTD. 00003600 C * P S * BOX 20, SITE 7, SS #1 00003700 C * * * * * * * CALGARY, ALTA. CANADA PH 288-8088 00003800 C ====================================================== 00003900 C 00004000 SUBROUTINE S2CSRJ (FILT, NPTS, SAMPR, LOWCUT, LOWPAS, HIGPAS, 00004100 * HIGCUT, IPR) 00004200 C 00004300 IMPLICIT INTEGER (A-Z) 00004400 C EXTERNAL S1ATP 00004500 C 00004600 C=================================================================== 00004700 C 00004800 C REAL ARRAYS IN PARAMETER LIST 00004900 C 00005000 REAL FILT (1) 00005100 C 00005200 C REAL VARIABLES -- LOCAL 00005300 C 00005400 REAL AL 00005500 DOUBLE PRECISION D1 00005600 DOUBLE PRECISION D2 00005700 DOUBLE PRECISION MINVAL 00005800 REAL HC 00005900 DOUBLE PRECISION HHC 00006000 REAL HP 00006100 REAL LC 00006200 DOUBLE PRECISION LLC 00006300 REAL LP 00006400 DOUBLE PRECISION PI 00006500 REAL POW 00006600 REAL SCF 00006610 DOUBLE PRECISION SS 00006700 DOUBLE PRECISION SSS 00006800 DOUBLE PRECISION T 00006900 DOUBLE PRECISION XK 00007000 DOUBLE PRECISION XN1 00007100 DOUBLE PRECISION XN2 00007200 DOUBLE PRECISION Y 00007300 DOUBLE PRECISION YSQ 00007400 C 00007411 C INITIALIZATION 00007412 DATA MINVAL /1.D-4/ 00007420 DATA PI /6.283185307179586/ 00007430 C 00007500 IF (1.EQ.2) CALL S1ATP C T = SAMPR / 1000. 00007600 M = NPTS / 2 00007700 L = M * 2 + 1 00007800 C 00007900 LC = LOWCUT 00008000 LP = LOWPAS 00008100 HP = HIGPAS 00008200 HC = HIGCUT 00008300 C 00008400 HHC = HC - HP 00008500 HHC = HHC * HHC 00008600 LLC = LP - LC 00008700 LLC = LLC * LLC 00008800 C 00008900 C SET MIDPOINT TO THE VALUE 1./NPTS 00009000 C 00009100 FILT(M+1) = 1./NPTS 00009200 SSS = 1./((HP+HC-LC-LP)*T*NPTS) 00009300 C 00009400 HP = HP * PI 00009500 HC = HC * PI 00009600 LP = LP * PI 00009700 LC = LC * PI 00009800 C 00009900 DO 20 00010000 * K = 1, M 00010100 XK = K - M - 1 00010200 C 00010300 Y = XK * T 00010400 YSQ = Y * Y 00010500 D1 = 1. - 4. * YSQ * HHC 00010600 D2 = 1. - 4. * YSQ * LLC 00010700 C 00010800 XN1 = DSIN ( Y * HP ) + DSIN ( Y * HC ) 00010900 XN2 = DSIN ( Y * LC ) + DSIN ( Y * LP ) 00011000 C 00011100 CC FOR CASES IN WHICH D1 IS SMALL: (DOUG NYMAN OCT 17, 1977) 00011200 CC 00011300 CC D1 = 1 - (2*Y (HC-HP))**2 = @ << 1 00011400 CC THIS IMPLIES, 00011500 CC HC = HP - (1-@)/(2*Y) WHICH IS APPROX.= HP - 1 / (2*Y) FOR(Y<0) 00011600 CC SO, 00011700 CC XN1 = SIN(2*PI*Y*HP) + SIN(2*PI*Y*HC) 00011800 CC = SIN(2*PI*Y*HP) + SIN(2*PI*Y*(HP-1/(2*Y)+@/(4*Y)) 00011900 CC WHICH IS APPROX. 00012000 CC = SIN(2*PI*Y*HP) + SIN(2*PI*Y*HP - PI + PI*@/2) 00012100 CC = SIN(2*PI*Y*HP) - SIN(2*PI*Y*HP + PI*@/2) 00012200 CC = SIN(2*PI*Y*HP) - SIN(S*PI*Y*HP) * COS(PI*@) 00012300 CC -COS(2*PI*Y*HP) * SIN(PI*@/2) 00012400 CC WHICH IS APPROX. 00012500 CC = - PI*@/2 * COS(2*PI*Y*HP) 00012600 CC OR, 00012700 CC XN1 IS APPROX. 00012800 CC = - 2*PI*D1/4 * COS(2*PI*Y*HP) 00012900 CC 00013000 CC AND SIMILARLY, 00013100 CC XN2 IS APPROX. 00013200 CC = - 2*PI*D2/4 * COS(2*PI*Y*LC) 00013300 CC 00013400 CC 00013500 CC 00013600 IF (DABS(D1) .GT. MINVAL) GO TO 10 00013700 IF (D1 .EQ. 0.0) D1 = MINVAL 00013800 XN1 = - 0.25 * PI * D1 * DCOS(Y*HP) 00013900 C 00014000 10 IF (DABS(D2) .GT. MINVAL) GO TO 15 00014100 IF (D2 .EQ. 0.0) D2 = MINVAL 00014200 XN2 = - 0.25 * PI * D2 * DCOS(Y*LC) 00014300 C 00014400 C 00014500 15 SS = SSS / (XK * PI) 00014600 FILT(K) = (XN1/D1-XN2/D2) * SS 00014700 FILT(L+1-K) = FILT(K) 00014800 C 00014900 20 CONTINUE 00015000 C 00015100 C NORMALIZE FILTER TO A POWER UNITY. 00015200 C 00015300 CCCC CALL ARSMFS (FILT,L,POW) 00015301 C 00015302 POW = 0.0 00015303 C 00015304 DO 25 I = 1, L 00015310 POW = POW + FILT(I)*FILT(I) 00015320 25 CONTINUE 00015330 C 00015340 AL = L 00015500 POW = POW * SQRT(AL) 00015600 CCCC CALL ARMPFC (FILT,FILT,SQRT(1./POW),L) 00015700 C 00015710 SCF = SQRT(1./POW) 00015711 C 00015712 DO 30 I = 1, L 00015720 FILT(I) = FILT(I)*SCF 00015730 30 CONTINUE 00015740 C 00015800 C 00015900 NPTS = L 00016000 RETURN 00016100 C 00016200 C 00016300 9000 FORMAT (/' *** S2CSRJ BANDPASS DESIGN. TOO HIGH REJECT/LENGTH', 00016400 *' RATIO -- ADJUSTMENTS MADE') 00016500 C 00016600 END 00016700