CTITLECEPNEW -- 3838 ARRAY PROCESSOR CEPSTRUM OPERATOR ROUTINE 00000010 C*********************************************************************** C COPYRIGHT ATLANTIC RICHFIELD COMPANY 1991 * C*********************************************************************** CA 00000020 CA AUTHOR R. C. DECKER 00000030 CA DESIGNER R. C. DECKER 00000040 CA LANGUAGE S/370 FORTRAN H EXTENDED 00000050 CA DATE 03/15/82 00000060 CA MODIFIED 10/1/82 ADDED PREDICTION DISTANCE 00000070 C 07/01/86 JMP. CHANGE ISTATE TO 0. CA CA REVISED FOR VS FORTRAN TRA 10-03-8400004400 CA CA RELEASED FOR PRODUCTION SPARC TRA 01-21-8500004400 CA CA 00000080 CA CALLING SEQUENCE: 00000090 CA 00000100 CA CALL CEPNEW(CEP, K, TRACE, N, P, LOP, AUTO, OPR, WAV, JPD) 00000110 CA 00000120 CA INPUT PARAMETERS: 00000130 CA 00000140 CA CEP - INPUT CEPSTRUM R4 00000150 CA K - FFT LENGTH (POWER OF 2--MIN=128 MAX=8192) I4 00000160 CA TRACE - INPUT TRACE / OUTPUT DECONVOLVED TRACE R4 00000170 CA N - TRACE LENGTH (MAX=8192) I4 00000180 CA P - WHITE NOISE R4 00000190 CA LOP - LENGTH OF THE OPERATOR I4 00000200 CA AUTO - OUTPUT AUTOCORRELATION R4 00000210 CA OPR - OUTPUT OPERATOR R4 00000220 CA WAV - OUTPUT WAVELET R4 00000230 CA JPD - PREDICTION DISTANCE IN SAMPLES I4 00000240 CA 00000250 CA 00000260 CA 00000270 CA PURPOSE: THIS ROUTINE DEVELOPS A MIN-PHASE OPERATOR FROM 00000280 CA THE CEPSTRUM AND APPLIES THE OPERATOR TO THE TRACE. 00000290 CA 00000300 CA 00000310 CA 00000320 CA NOTE: THE TRACE LENGTH PLUS THE OPERATOR LENGTH CANNOT 00000330 CA BE GREATER THAN 8192. 00000340 CA 00000341 CA AUTOCORRELATION LENGTH < OR = K/2 00000342 CA OPERATOR LENGTH < OR = K/2 00000343 CA WAVELET LENGTH < OR = K/2 00000344 CA 00000345 CA 00000350 CA 00000360 C PROCEDURE: 00000370 C 00000380 C 00000390 C 1. FFT CEPSTRUM C (LENGTH K) 00000400 C 00000410 C 2. DOUBLE THE REAL(C) 00000420 C 00000430 C 3. SET IMAG(C) TO 0 00000440 C 00000450 C 4. C = EXP(C) 00000460 C 00000470 C 5. INVERSE FFT(C) 00000480 C 00000490 C 6. REPLACE C(1) = C(1)*(1+P) 00000500 C 00000510 C 7. DEVELOP THE OPERATOR WITH AUTOCORRELATION IN C 00000520 C 00000530 C 8. MODIFY THE OPERATOR FOR PREDICTION DISTANCE 00000540 C 00000550 C 9. SCALE THE OPERATOR 00000560 C 00000570 C 10. CONVOLVE THE TRACE WITH THE OPERATOR (USE FFT) 00000580 C 00000590 C 00000600 C BULK STORAGE: 00000610 C ______________________________ 00000620 C |-----0|____________________________|WHITE NOISE (1+P) 00000630 C | 1| | 00000640 C | | | 00000650 C | | TRACE SEGMENT | 00000660 C | | (N) |R5=LOP 00000670 C | | |R6=K 00000680 C | | |R7=K/2+1 00000690 C | | |R8=1.0/FLOAT(K) 00000700 C | | |R9=1.0 00000710 C | | | 00000720 C | |____________________________|___ 00000730 C | |R1=N+1 | | 00000740 C | | | | 00000750 C | | AUTOCORRELATION | | 00000760 C OUTPUT| | (K/2) | | 00000770 C AREA| |____________________________| CEPSTRUM INPUT 00000780 C | |R3=R1+K/2 |(K) 00000790 C | | OPERATOR | | 00000800 C | | (K/2) | | 00000810 C | | | | 00000820 C | |____________________________|___ 00000830 C | |R13=R1+K | 00000840 C | | WAVELET | 00000850 C | | (K/2) | 00000860 C | | | 00000870 C |______|____________________________| 00000880 C |R2=R1+K*2-K/2 | 00000890 C | | 00000900 C | CEPSTRUM & TRACE | 00000910 C | FFT AREA AREA | 00000920 C | (K+2) (NN) | 00000930 C | | 00000940 C | | 00000950 C | | 00000960 C | | 00000970 C |____________________________| 00000980 C |R4=R2+NN | 00000990 C | |R11=NN 00001000 C | |R12=NN/2+1 00001010 C | |R15=1.0/FLOAT(NN) 00001020 C | | 00001030 C | OPERATOR AREA | 00001040 C | (NN) | 00001050 C | | 00001060 C | | 00001070 C |____________________________| 00001080 C |R13=R4+NN | 00001090 C | | 00001100 C | | 00001110 C | | 00001120 C | FFT AREA FOR TRACE | 00001130 C | (NN+2) | 00001140 C | | 00001150 C | | 00001160 C | | 00001170 C |____________________________| 00001180 C |R14=R13+NN+2 | 00001190 C | | 00001200 C | | 00001210 C | | 00001220 C | FFT AREA FOR OPR | 00001230 C | (NN+2) | 00001240 C | | 00001250 C | | 00001260 C | | 00001270 C |____________________________| 00001280 C 00001290 C 00001300 C 00001310 C 00001320 C 00001330 C 00001340 SUBROUTINE CEPNEW(CEP,K,TRACE,N,P,LOP,AUTO,OPR,WAV,JPD) 00001350 C 00001360 DIMENSION APBUF(16385), AUTO(1), OPR(1), CEP(1), WAV(1) 00001370 DIMENSION TRACE(1), CIT(400), IREG(15), REGS(15) 00001380 C 00001390 REAL*8 CCW(200) 00001400 C 00001410 INTEGER R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14,R15 00001420 DATA R1/1/,R2/2/,R3/3/,R4/4/,R5/5/,R6/6/,R7/7/,R8/8/,R9/9/,R10/10/00001430 DATA R11/11/,R12/12/,R13/13/,R14/14/,R15/15/ 00001430 DATA ISTAT/3/, ISTATE/0/, L150/150/ C 00001460 C 00001470 C 00001480 EQUIVALENCE(IREG(1),REGS(1)) 00001490 EQUIVALENCE (IREG(1), NR1 ) 00001500 EQUIVALENCE (IREG(2), NR2 ) 00001510 EQUIVALENCE (IREG(3), NR3 ) 00001520 EQUIVALENCE (IREG(4), NR4 ) 00001530 EQUIVALENCE (IREG(5), NR5 ) 00001540 EQUIVALENCE (IREG(6), NR6 ) 00001550 EQUIVALENCE (IREG(7), NR7 ) 00001560 EQUIVALENCE (IREG(8), XR8 ) 00001570 EQUIVALENCE (IREG(9), XR9 ) 00001580 EQUIVALENCE (IREG(10), XR10 ) 00001590 EQUIVALENCE (IREG(11), NR11 ) 00001600 EQUIVALENCE (IREG(12), NR12 ) 00001610 EQUIVALENCE (IREG(13), NR13 ) 00001620 EQUIVALENCE (IREG(14), NR14 ) 00001630 EQUIVALENCE (IREG(15), XR15 ) 00001640 C 00001650 DATA IFIRST/0/ 00001660 C 00001670 IC = ALOG(FLOAT(N+LOP+JPD)) / ALOG(2.) 00001680 NN = 2**(IC+1) 00001690 C 00001700 IF(IFIRST.NE.0) GO TO 100 00001710 C 00001720 C GET A UNIT NUMBER FOR THE ARRAY PROCESSOR 00001730 C 00001740 C 00001750 KP1 = 2*8192 + 1 00001760 NWDS = 5*NN + 2*K + 5 00001770 IF(NWDS.LT.16385) NWDS = 16385 00001780 CALL CSAPUN(NWDS,IU) 00001790 C 00001800 C 00001810 C START CONSTRUCTION OF THE 3838 PROGRAM 00001820 C 00001830 CALL VPSS(IU,'BLD ',ISTAT,CCW,200,CIT,400) 00001840 CALL VPSS(IU,'XWR ',REGS,15,1) 00001850 CALL VPSS(IU,'VPUT',APBUF,KP1,0,0) 00001860 C 00001870 C FFT CEPSTRUM 00001880 C 00001890 CALL VPSS(IU,'FFTR',ISTATE, 00001900 * 96, 0, 0, R2, R7, 00001910 * 96, 0, 0, R1, R6) 00001920 C 00001930 C FFT OUTPUT IN ALREADY DOUBLED SO 00001940 C IT IS NOT NECESSARY TO DOUBLE IT. 00001950 C 00001960 C EXPONENTIATE THE REAL COMPONENTS 00001970 C 00001980 CALL VPSS(IU,'EXP ',ISTATE, 00001990 * 96, 0, 0, 2, R2, R7, 00002000 * 32, 0, 2, R2) 00002010 C 00002020 C ZERO OUT THE IMAG COMPONENT 00002030 C 00002040 CALL VPSS(IU,'ZMV ', ISTATE, 00002050 * 96, 1, 0, 2, R2, R7) 00002060 C 00002070 C SCALE AND CONJUGATE 00002080 C 00002090 CALL VPSS(IU,'XMVX', 0, 1, R1, R8) 00002100 CALL VPSS(IU,'ZMV ', ISTATE, 00002110 * 32, 1, 1, 1, R1) 00002120 C 00002130 C NORMALIZE 00002140 C 00002150 CALL VPSS(IU, 'SCMC', ISTATE, 00002160 * 96, 0, 0, R2, R7, 00002170 * 32, 0, R2, 00002180 * 32, 0, R1) 00002190 C 00002200 C INVERSE FFT 00002210 C 00002220 CALL VPSS(IU, 'IFTR', ISTATE, 00002230 * 96, 0, 0, R1, R6, 00002240 * 96, 0, 0, R2, R7) 00002250 C 00002260 C SCALE ZERO LAG VALUE BY 1.0 + P 00002270 C 00002280 CALL VPSS(IU, 'SMY ', ISTATE, 00002290 * 32, 0, 1, 1, R1, 00002300 * 32, 0, 1, R1, 00002310 * 0, 0) 00002320 C 00002322 C 00002324 C 00002330 C CALCULATE INDEX FOR PREDICTION DISTANCEINTO AC 00002340 C 00002360 CALL VPSS(IU, 'XMV ', R8, R13) 00002370 CALL VPSS(IU, 'XAD ', R13, R1) 00002380 C 00002400 C MOVE ZEROS INTO AREA FOLLOWING THE AC SO 00002410 C WHEN THE PREDICTION DISTANCE ADJUSTMENT IS 00002420 C MADE THE END OF THE AC WILL BE ZEROS. 00002430 C 00002440 CALL VPSS(IU, 'ZMV ', ISTATE, 00002450 * 96, 0, -1, 1, R3, R7) 00002460 C 00002470 C CALL THE WEINER-LEVISON TO DEVELOP THE OPERATOR 00002480 C 00002500 00002510 C DEVELOP THE OPERATOR 00002520 C 00002530 CALL VPSS(IU, 'WLEV', ISTATE, 00002540 * 96, 0, 0, R2, R5, 00002550 * 32, 0, 1, R1, 00002560 * 32, 0, 1, R13) 00002570 C 00002580 C CONSTRUCT THE NEW OPERATOR AND 00002590 C MOVE THE CONSTANT 1.0 TO 0 OFF R3 00002600 C 00002610 CALL VPSS(IU, 'XMVX', 0, 1, R3, R9) 00002620 C 00002630 C 00002690 C CALCULATE ADDRESS IN OPERATOR TO MOVE WLEV 00002700 C OPERATOR ADJUSTED FOR THE PREDICTION DISTANCE 00002710 C 00002720 CALL VPSS(IU, 'XMV ',R13, R8) 00002730 CALL VPSS(IU, 'XMV ', R6, R5) 00002740 CALL VPSS(IU, 'XAD ', R6, R13) 00002750 CALL VPSS(IU, 'XAD ', R13, R3) 00002760 CALL VPSS(IU, 'VMV ', ISTATE, 00002770 * 96, 0, 0, 1, R13, R5, 00002780 * 40, 0, 1, R2) 00002790 C 00002800 C 00002802 C 00002804 C 00002808 C THE OPERATOR IS NOW MODIFIED AND STORED AT R3 00002810 C 00002811 C 00002813 C 00002820 C 00002830 C==================================================================== 00002840 C 00002850 C 00002860 C 00002880 C 00002890 C==================================================================== 00002900 C 00002910 C 00002920 C SET UP REGISTERS FOR CONVOLUTION 00002930 C 00002931 C MOVE CONSTANT 1.0 TO ADDRESS 0 00002932 C 00002933 CALL VPSS(IU, 'XMVX', 0, 1, 0, R9) 00002936 C 00002940 C 00002950 CALL VPSS(IU, 'XMV ', R9, R7) 00002960 CALL VPSS(IU, 'XSB ', R9, R6) 00002970 CALL VPSS(IU, 'XMV ', R13, R3) 00002980 CALL VPSS(IU, 'XAD ', R13, R7) 00002990 CALL VPSS(IU, 'XSBI', R13, 1) 00003000 C 00003001 C ZERO OUT THE WAVELET ARRAY SO IF CVM CALCULATES 00003002 C FEWER THAN (LOP+JPD) SAMPLES THE REMAINING 00003003 C SAMPLES WILL BE ZERO. 00003004 C 00003005 CALL VPSS(IU, 'ZMV ', ISTATE, 00003006 * 96, 0, -1, 1, R13, R7) 00003007 C 00003010 C CONVOLVE THE TIME REVERSED OPERATOR WITH THE AC 00003020 C TO PRODUCE THE WAVELET 00003030 C 00003040 C 00003048 CALL VPSS(IU, 'CVM ', ISTATE, 00003050 * 96, 0, 0, 1, R13, R9 , 00003060 * 96, 0, -1, 1, R1, R7 , 00003070 * 96, 0, 0, 1, R3, R6 ) 00003080 C 00003081 C SKIP SCALING OF OPERATOR IF PREDICTION 00003082 C 00003083 CALL VPSS(IU, 'XCI ', R8, 1, L150, 'GT ') 00003085 C 00003090 C 00003100 C DOT PRODUCT KW 00003110 C 00003120 CALL VPSS(IU, 'SSQ ', ISTATE, 00003130 * 32, 0, R2, 00003140 * 96, 0, 0, 1, R13, R6) 00003150 C 00003160 C DIVIDE BY THE ZERO LAG OF AC 00003170 C 00003180 CALL VPSS(IU, 'SDIV', ISTATE, 00003190 * 32, 0, 1, 1, R2, 00003200 * 32, 0, 1, R1, 00003210 * 32, 0, R2) 00003220 C 00003230 C TAKE SQRT TO GET K 00003240 C 00003250 CALL VPSS(IU, 'SQRT', ISTATE, 00003260 * 32, 1, 1, 1, R2, 00003270 * 32, 0, 1, R2) 00003280 C 00003290 C 00003300 C 00003301 CALL VPSS(IU, 'SDIV', ISTATE, 00003302 * 32, 1, 1, 1, R2, 00003303 * 32, 1, 1, R2, 00003304 * 0, 0) 00003305 C 00003306 C 00003307 C 00003310 C SCALE THE OPERATOR BY K 00003320 C 00003330 CALL VPSS(IU, 'SMY ', ISTATE, 00003340 * 96, 0, 0, 1, R3, R6, 00003350 * 32, 0, 1, R3, 00003360 * 32, 1, R2) 00003370 C 00003372 C 00003374 C 00003380 C 00003390 CALL VPSS(IU, 'XID ', L150) 00003400 C 00003410 C 00003420 C RESET R13 FOR CONVOLVING TRACE AND OPERATOR 00003430 C 00003440 CALL VPSS(IU, 'XMV ', R13, R4) 00003450 CALL VPSS(IU, 'XAD ', R13, R11) 00003460 C 00003470 C 00003480 C 00003490 C CONVOLVE THE OPERATOR WITH THE TRACE 00003500 C USE FFTS TO DO THE CONVOLUTION 00003510 C 00003520 C 00003530 C ZERO OUT THE FFT AREA FOR THE TRACE 00003540 C 00003550 CALL VPSS(IU, 'ZMV ', ISTATE, 00003560 * 96, 0, 0, 1, R2, R11) 00003570 C 00003580 C MOVE THE TRACE TO THE FFT AREA 00003590 C 00003600 CALL VPSS(IU, 'VMV ', ISTATE, 00003610 * 96, 0, -1, 1, R2, R1, 00003620 * 0, 1) 00003630 C 00003640 C ZERO OUT THE FFT AREA FOR THE OPERATOR 00003650 C 00003660 CALL VPSS(IU, 'ZMV ', ISTATE, 00003670 * 96, 0, 0, 1, R4, R11) 00003680 C 00003690 C MOVE THE OPERATOR TO THE FFT AREA 00003700 C 00003710 CALL VPSS(IU, 'VMV ', ISTATE, 00003720 * 96, 0, 0, 1, R4, R6, 00003730 * 32, 0, 1, R3) 00003740 C 00003750 C FFT THE TRACE 00003760 C 00003770 CALL VPSS(IU, 'FFTR', ISTATE, 00003780 * 96, 0, 0, R13, R12, 00003790 * 96, 0, 0, R2, R11) 00003800 C 00003810 C FFT THE OPERATOR 00003820 C 00003830 CALL VPSS(IU, 'FFTR', ISTATE, 00003840 * 96, 0, 0, R14, R12, 00003850 * 96, 0, 0, R4, R11) 00003860 C 00003870 CALL VPSS(IU, 'XMVX', 0, 1, R2, R10) 00003880 C 00003881 C SCALE THE FFT BY 0.5 00003882 C 00003885 C 00003890 CALL VPSS(IU, 'SMY ', ISTATE, 00003900 * 96, 0, 2, 1, R13, R11, 00003910 * 32, 0, 1, R13, 00003920 * 32, 0, R2) 00003930 C 00003940 CALL VPSS(IU, 'SMY ', ISTATE, 00003950 * 96, 0, 2, 1, R14, R11, 00003960 * 32, 0, 1, R14, 00003970 * 32, 0, R2) 00003980 C 00003990 C 00004000 C 00004020 C MULTIPLY THE SPECTRUM 00004030 C 00004040 CALL VPSS(IU, 'CMCO', ISTATE, 00004050 * 96, 0, 0, R13, R12, 00004060 * 32, 0, R13, 00004070 * 32, 0, R14) 00004080 C 00004090 C 00004100 C 00004110 C INVERSE FFT 00004120 C 00004130 CALL VPSS(IU, 'IFTR', ISTATE, 00004140 * 96, 0, 0, R2, R11, 00004150 * 96, 0, 0, R13, R12) 00004160 C 00004170 C 00004180 C 00004190 C 00004200 C NORMALIZE THE FFT - MOVE SCALE FACTOR TO MEMORY 00004210 C 00004220 CALL VPSS(IU, 'XMVX', 0, 1, R13, R15) 00004230 C 00004240 C NORMALIZE THE FFT 00004250 C 00004260 CALL VPSS(IU, 'SMY ', ISTATE, 00004270 * 96, 0, 0, 1, R2, R11, 00004280 * 32, 0, 1, R2, 00004290 * 32, 0, R13) 00004300 C 00004310 C MOVE OUTPUT TRACE INTO OUTPUT AREA 00004320 C 00004330 CALL VPSS(IU, 'VMV ', ISTATE, 00004340 * 64, 1, -1, 1, R1, 00004350 * 32, 0, 1, R2) 00004360 C 00004370 C 00004380 C 00004390 C 00004400 C 00004410 50 CONTINUE 00004420 C 00004430 CALL VPSS(IU, 'XRD ', REGS, 15, 1) 00004440 CALL VPSS(IU, 'VGET', APBUF, KP1, 1, 0) 00004450 CALL VPSS(IU, 'XLTE', CEP2) 00004460 C 00004470 C 00004480 C 00004490 C 00004500 C 00004510 IFIRST = 1 00004520 C 00004530 C SET UP REGISTERS FOR 3838 00004540 C 00004550 100 NR1 = N + 1 00004560 NR2 = NR1 + K*2 00004570 NR3 = NR1 + K/2 00004580 NR4 = NR2 + NN 00004590 NR5 = LOP 00004600 NR6 = K 00004610 NR7 = K/2 + 1 00004620 XR8 = 1.0/FLOAT(K) 00004630 XR9 = 1.0 00004640 XR10 = 0.5 00004650 NR11 = NN 00004660 NR12 = NN/2 + 1 00004670 NR13 = JPD 00004680 NR14 = NR4 + 2*NN + 2 00004690 XR15 = 1.0/FLOAT(NN) 00004700 C 00004710 C 00004720 C 00004730 C 00004740 C 00004750 C SET UP THE APBUF 00004760 C 00004770 APBUF(1) = 1. + P 00004780 CALL ARMVE(TRACE,APBUF(2),N) 00004790 CALL ARMVE(CEP,APBUF(N+2),K) 00004800 C 00004810 C NPK = N + K 00004820 C WRITE(6,8010) 00004830 C8010 FORMAT(//,' APBUF BEFORE CEPNEW ') 00004840 C WRITE(6,8020) (APBUF(I),I=1,NPK) 00004850 C 00004860 C 00004870 C 00004880 C EXECUTE THE 3838 PROGRAM 00004890 C 00004900 CALL VPSS(IU,'EXCW',CEP2) 00004910 C 00004920 C 00004930 C 00004940 C WRITE(6,9006) (IREG(I),I=1,15) 00004941 C9006 FORMAT(5X, 15I8) 00004942 C 00004943 C 00004944 C 00004945 C 00004946 C NREND = NR14 + NN + 2 00004950 C WRITE(6,8030) 00004960 C8030 FORMAT(' APBUF AFTER CEPNEW ') 00004970 C WRITE(6,8020) (APBUF(I),I=1,NREND) 00004980 C8020 FORMAT(5X,10E12.3) 00004990 C 00005000 C 00005010 C 00005020 C 00005030 C MOVE THE OUTPUT INTO ARRAYS 00005040 C 00005050 C 00005060 NP = LOP+JPD 00005070 CALL ARMVE(APBUF,TRACE,N) 00005080 CALL ARMVE(APBUF(N+1),AUTO,NP) 00005090 CALL ARMVE(APBUF(N+K/2+1),OPR,NP) 00005100 CALL ARMVE(APBUF(N+K+1),WAV,NP) 00005110 C 00005130 C 00005140 C 00005150 C 00005160 RETURN 00005170 END 00005180