+PATCH,MYWV0LIB. +KEEP,WVLTDATA. C************************************************************************ C Main common block for storing WAVELET transformed data * C * C now in YPATCHY @Erik Kats. Feb 1998 * C************************************************************************ C Here we define the pointer assignment scheme for a WAVELET * C (HAAR) transformation: * C Original signal: N samples = 2**J numbered from 1..N or 0..(N-1)* C Resolution scales: J numbered from 1..J * C Scale Signal (Approximation) Details (Haar coef) * C 1 N=2**J 2**(J-1) * C 2 N/2=2**(J-1) 2**(J-2) * C ... * C j N/[2**(j-1)]=2**(J-j+1) 2**(J-j) * C ... * C J-1 N/[2**(J-2)]=2**2 2**1=2 * C J N/[2**(J-1)]=2 2**0=1 * C 1 (DC) * C Total number of resulting coefficients after Haar transform: * C 1(DC)+1+2+2**2+...+2**(J-2)+2**(J-1)=2**J * C************************************************************************ C Notice that the number below should (suggested to) be equal to NSECS of C the MYSPELIB patch. DOUBLE PRECISION PLSTSQR ! the sum of squares for RAW waveform INTEGER NSMPL ! that's the maximum number of samples PARAMETER (NSMPL=262144)! i.e., 262 microsecs... ASSUMING TICK=1ns C or for 5ns sampling 1.31 ms C Notice: all samples, start from 1 (NOT zero) REAL PLSTARRAY(NSMPL) ! this is the RAW waveform REAL WVLTARRAY(NSMPL) ! this is the WAVELET transformed waveform C INTEGER PLSTDIM ! total number of samples C we are not ready to process PLST and WVLT arrays of different dimension INTEGER WVLTDIM ! total number of samples ==> WVLT coef's REAL S ! this is the square root of 2... INTEGER LEVEL ! total number of resolution scales REAL PLSTSUM ! integral for RAW waveform COMMON/WVLTDATA/PLSTSQR,PLSTARRAY,WVLTARRAY,WVLTDIM,S, + LEVEL,PLSTSUM +KEEP,WVLTFEAT. C************************************************************************ C Main common block for storing features of the WAVELET * C transformed data * C PRETRG1 is a number that reflects the alignment of the * C wavelet maxima over all the resolution scales * C PRETRG3 is the scale where the maximum of the wavelet maxima * C occurs * C GMAX is the maximum of ALL wavelet coefficient among all * C scales * C DENER is the difference in energy content between High and * C Low scales * C LOCGMAX is the location of the location of the GMAX * C INFTRG is an array with FOUR variables for each resolution * C scale * C * C The above are the "traditional" wavelet features from the * C original Haar code. * C * C New features added: * C RS_DENER is the number of scales used for constructing DENER * C (this is actually used for DENER) * C WVLTENER is the sum of squares of wavelet coefs for a give * C scale; (1,I) is the raw sqr sum, (2,I) is "normalized" * C WVLTSQR the sum of squares for WAVELET waveform * C WVLTSUM wave-DC * C * C New features added: * C DECO_CPU is the cpu in 10's of miliseconds used for WAVELET * C decomposition * C RECO_CPU is the cpu in 10's of miliseconds used for WAVELET * C recomposition * C MSE is mean square error between the original WAVEFORM and * C the one reconstructed with the inverse WAVELET transform * C * C @Erik Kats. 1994 * C now in YPATCHY @Erik Kats. Feb 1998 * C * C We now introduced the "second maxima" processing: * C PRETRG12 is a number that reflects the alignment of the * C SECOND wavelet maxima over all the resolution scales * C PRETRG32 is the scale where the maximum of the SECOND wavelet * C maxima occurs * C GMAX2 is the maximum of ALL SECOND MAX wavelet coefficients * C among all scales * C DENER2 is the difference in energy content between High and * C Low scales constructed by the SECOND WAVELET max * C LOCGMAX2 is the location of the location of the GMAX2 * C * C @Erik Kats. Mar 1998 * C************************************************************************ INTEGER MAXRES ! maximum number of resolution scales PARAMETER (MAXRES=100) DOUBLE PRECISION WVLTENER(2,MAXRES) ! energy for each scale DOUBLE PRECISION WVLTSQR,MSE ! the sum of squares... INTEGER PRETRG1,PRETRG3,LOCGMAX REAL GMAX,DENER REAL INFTRG(4,MAXRES) ! 4 features for each resolution scale INTEGER RS_DENER REAL WVLTSUM ! wave-DC INTEGER*4 DECO_CPU,RECO_CPU ! CPU usage INTEGER PRETRG12,PRETRG32,LOCGMAX2 REAL GMAX2,DENER2 COMMON/WVLTFEAT/WVLTENER,WVLTSQR,MSE, + PRETRG1,PRETRG3,LOCGMAX,GMAX,DENER,INFTRG, + RS_DENER,WVLTSUM, + DECO_CPU,RECO_CPU, + PRETRG12,PRETRG32,LOCGMAX2,GMAX2,DENER2 C************************************************************************ C SUBROUTINE WVLTDATA_INIT * C SUBROUTINE WVLTDATA_DUMP(ILUN) * C SUBROUTINE WVLTFEAT_INIT * C SUBROUTINE WVLTFEAT_DUMP(ILUN) * C SUBROUTINE PREP_WVLTD * C SUBROUTINE EVLTCOMP(WVLTARRAY,WVLTDIM,WVPED,WVSAMPLES,PEDINF) * C SUBROUTINE WVLTCOMP(WVLTARRAY,WVLTDIM,WVTHRES,WVPED,WVREJ) * C * C now in YPATCHY @Erik Kats. Feb 1998 * C************************************************************************ +DECK,WVLTDATA_INIT. SUBROUTINE WVLTDATA_INIT C************************************************************************ C It initializes common block WVLTDATA * C * C @Erik Kats. Feb 1998 * C************************************************************************ IMPLICIT NONE +SEQ,WVLTDATA. +SEQ,GLOBALDT. INTEGER I C INTEGER LEVEL0 C 5 WRITE(*,*) 'Enter number of samples for WAVELET array (power of 2!):' C READ(*,*) WVLTDIM C C IF (WVLTDIM.GT.NSMPL) THEN C WRITE(*,7) NSMPL,INT(LOG10(FLOAT(NSMPL))/LOG10(2.0)+0.5) C GOTO 5 C END IF C C 7 FORMAT(1X,'Maximum Haar array allocation:',I7, C + ' i.e., over',I3,' resolution scales.') C LEVEL=INT(LOG10(FLOAT(WVLTDIM))/LOG10(2.0)+0.5) C LEVEL0=INT(LOG10(FLOAT(WVLTDIM))/LOG10(2.0)) C check if input HAAR dimension if a power of two... C IF (WVLTDIM.NE.2**LEVEL) THEN C WRITE(*,10) WVLTDIM C WRITE(*,20) WVLTDIM,2**LEVEL0,2**(LEVEL0+1) C GOTO 5 C END IF C C 10 FORMAT(1X,'The number you entered (',I10,') is NOT', C + ' a power of 2.') C 20 FORMAT(1x,'The closest to (',I10,') powers of 2 are:',2I11) C we grab the WAVELET dimension from GLOBAL common block... WVLTDIM=GWVLTDIM LEVEL=INT(LOG10(FLOAT(WVLTDIM))/LOG10(2.0)+0.5) S=0.707106781186547524400844362104849039 C S=SQRT(2.0)/2.0 DO I=1,WVLTDIM ! loop over all USED samples PLSTARRAY(I)=0. WVLTARRAY(I)=0. END DO ! close loop over all USED samples PLSTSQR=0.0 ! RAW waveform sum of squares PLSTSUM=0.0 ! RAW waveform integral RETURN END +DECK,WVLTDATA_DUMP. SUBROUTINE WVLTDATA_DUMP(ILUN) C************************************************************************ C It dumps to logical unit common block WVLTDATA * C * C @Erik Kats. Feb 1998 * C************************************************************************ IMPLICIT NONE +SEQ,WVLTDATA. INTEGER ILUN C INTEGER I,TMPLEVEL INTEGER I,J IF (ILUN.LE.0) RETURN ! return if no valid LUN WRITE(ILUN,10) WVLTDIM,LEVEL 10 FORMAT(1X,'Total processed samples:',I7,' over ',I3, + ' resolution scales') WRITE(ILUN,20) PLSTSQR,PLSTSUM 20 FORMAT(1X,'Sum of squares/Integral:',2F15.2) WRITE(ILUN,*) ' Sample/ RAW-WFD/', + ' HAAR-COEFF/ ThisLevel/' C DO I=1,WVLTDIM ! loop over processed samples C we loop over the scales rather the buckets... I=0 J=1 WRITE(ILUN,*) J,PLSTARRAY(J),WVLTARRAY(J),I DO I=1,LEVEL ! loop over scales... C TMPLEVEL=INT(LOG10(FLOAT(I))/LOG10(2.0)) C TMPLEVEL=INT(LOG10(FLOAT(I))/LOG10(2.0)+0.5) DO J=1+2**(I-1),2**I ! loop over samples of each scale WRITE(ILUN,*) J,PLSTARRAY(J),WVLTARRAY(J),I END DO ! close loop over samples of each scale END DO ! close loop over processed samples RETURN END +DECK,WVLTFEAT_INIT. SUBROUTINE WVLTFEAT_INIT C************************************************************************ C It INITIALIZES the content of sequence WVLTFEAT * C * C @Erik Kats. Feb 1998 * C************************************************************************ IMPLICIT NONE +SEQ,WVLTDATA. +SEQ,WVLTFEAT. INTEGER I PRETRG1=0 ! initialize output variables PRETRG3=0 ! initialize output variables GMAX =0.0 ! initialize output variables DENER =0.0 ! initialize output variables LOCGMAX=0 ! initialize output variables DO I=1,LEVEL ! here we initialize the INFTRG array INFTRG(1,I)=0. ! wavelet max INFTRG(2,I)=0. ! wavelet max position INFTRG(3,I)=0. ! wavelet max2 INFTRG(4,I)=0. ! wavelet max2 position END DO ! close initialize INFTRG array C Here we found out how many resolution scales to use for constructing C the energy contrast. Out of the LEVEL resolution scales, the last one C j=LEVEL is just a single coef and is ignored. That leaves us with LEVEL-1 C coefficients. RS_DENER=INT(((LEVEL-1)/2)) DO I=1,LEVEL ! here we initialize the WVLTENER array WVLTENER(1,I)=0. WVLTENER(2,I)=0. END DO ! close init WVLTENER array WVLTSQR=0. WVLTSUM=0. DECO_CPU=0 RECO_CPU=0 MSE=0. C features of the "second MAX" PRETRG12=0 ! initialize output variables PRETRG32=0 ! initialize output variables GMAX2 =0.0 ! initialize output variables DENER2 =0.0 ! initialize output variables LOCGMAX2=0 ! initialize output variables RETURN END +DECK,WVLTFEAT_DUMP. SUBROUTINE WVLTFEAT_DUMP(ILUN) C************************************************************************ C It dumps to logical unit ILUN the content of sequence WVLTFEAT * C It actually adds some info from WVLTDATA (PLSTSQRT,PLSTSUM...) * C * C @Erik Kats. Feb 1998 * C************************************************************************ IMPLICIT NONE +SEQ,WVLTDATA. +SEQ,WVLTFEAT. INTEGER ILUN INTEGER I IF (ILUN.LE.0) RETURN ! return if no valid LUN WRITE(ILUN,10) WVLTDIM,LEVEL 10 FORMAT(1X,'Total processed samples:',I7,' over ',I3, + ' resolution scales') WRITE(ILUN,89) DO I=1,LEVEL WRITE(ILUN,90) I,INFTRG(1,I),INT(INFTRG(2,I)), + INFTRG(3,I),INT(INFTRG(4,I)), + WVLTENER(1,I),WVLTENER(2,I) END DO WRITE(ILUN,91) PRETRG1,PRETRG1 WRITE(ILUN,92) PRETRG3,GMAX,LOCGMAX WRITE(ILUN,93) DENER,RS_DENER,LEVEL-RS_DENER,LEVEL-1 C 89 FORMAT(1X,'Level First MAX Location Second MAX Location') C Level First MAX/Location Second MAX/Location Total Energy Tot NormEner C 1 18.41 1234567 18.41 1234567x123456789012x123456789012 89 FORMAT(1X,'Level First MAX/Location Second MAX/' + 'Location Total Energy Tot NormEner') C 90 FORMAT(1X,I7,1X,F9.2,1X,I8,1X,F10.2,1X,I8) 90 FORMAT(1X,I5,F13.4,I7,F13.4,I7,F15.2,F15.2) 91 FORMAT(1X,' PRETRG1=',I6,' octal:',O) 92 FORMAT(1X,' PRETRG3=',I4,' GMAX=',F13.4, + ' LOCGMAX=',I7) 93 FORMAT(1X,' DENER=',F15.2,' scales 1 to',I3,' and',I3,' to',I3) WRITE(ILUN,191) PRETRG12,PRETRG12 WRITE(ILUN,192) PRETRG32,GMAX2,LOCGMAX2 WRITE(ILUN,193) DENER2,RS_DENER,LEVEL-RS_DENER,LEVEL-1 191 FORMAT(1X,' PRETRG12=',I6,' octal:',O) 192 FORMAT(1X,' PRETRG32=',I4,' GMAX2=',F13.4, + ' LOCGMAX2=',I7) 193 FORMAT(1X,' DENER2=',F15.2,' scales 1 to',I3,' and',I3,' to',I3) WRITE(ILUN,20) WVLTSQR,PLSTSQR 20 FORMAT(1X,'Sum of squares: WAVELET/SOURCE:',2F15.2) WRITE(ILUN,30) WVLTSUM,PLSTSUM 30 FORMAT(1X,'WAVELET-DC/SOURCE-INTEGRAL: ',2F15.2) WRITE(ILUN,40) FLOAT(DECO_CPU)/100.,FLOAT(RECO_CPU)/100. 40 FORMAT(1X,'Wavelet Deco/Reco CPU (in secs):',2F10.2) WRITE(ILUN,*) 'With a mean square error: ',MSE RETURN END +DECK,PREP_WVLTD. SUBROUTINE PREP_WVLTD C************************************************************************ C Prepares the array to feed into the wavelet algorithm. * C It practically prepares the array PLSTARRAY with the data * C before wavelet-transform it. * C * C @Erik Kats. Feb 1998 * C************************************************************************ IMPLICIT NONE +SEQ,PLSTRAIN. +SEQ,PLSTCON1. +SEQ,WVLTDATA. INTEGER I INTEGER SBIN ! start of digitized buckets in HAAR IF (FLOAT(WVLTDIM)*TICK*FLOAT(NBIN).LT.DURATION) THEN WRITE(*,*) 'Warning: pulse width is greater than WV window!', + ' Pulse will be trancated!' SBIN=0 DO I=1,WVLTDIM PLSTARRAY(I)=DIGIPULSE(I) WVLTARRAY(I)=PLSTARRAY(I) END DO RETURN ENDIF SBIN=INT((FLOAT(WVLTDIM)*TICK*FLOAT(NBIN)-DURATION) + /(2.0*TICK*FLOAT(NBIN))) DO I=1,SBIN ! loop over pedestal presamples PLSTARRAY(I)=0. END DO ! close loop over pedestal presamples DO I=SBIN+1,SBIN+DLASTSMPL ! loop over digitized samples PLSTARRAY(I)=DIGIPULSE(I-SBIN) END DO ! close over digitized samples DO I=SBIN+DLASTSMPL+1,WVLTDIM ! loop beyond end of digitized pulse PLSTARRAY(I)=0. END DO ! close loop to the Haar limit C Since all WAVELET transforms return the WAVELET vector on the same C array, we INITIALIZE the WAVELET array with the RAW array... C we copy the RAW data to the WAVELET vector... DO I=1,WVLTDIM WVLTARRAY(I)=PLSTARRAY(I) END DO RETURN END C+QUIT. +DECK,WVLTCOMP. SUBROUTINE WVLTCOMP(WVLTARRAY,WVLTDIM,WVTHRES,WVPED,WVREJ) C************************************************************************ C It accepts a wavelet array WVLTARRAY of dimension WVLTDIM * C and sets to WVPED all the WAVELET coefficients with ABSOLUTE * C values LESS than WVTHRES. * C If by any chance a NEGATIVE WVTHRES is passed, its ABS value * C is used. * C WVREJ is the number of WV coefficients set to WVPED * C * C @Erik Kats. Mar 1998 * C************************************************************************ IMPLICIT NONE INTEGER NSMPL ! that's the maximum number of samples PARAMETER (NSMPL=262144)! i.e., 262 microsecs... ASSUMING TICK=1ns C or for 5ns sampling 1.31 ms C Notice: all samples, start from 1 (NOT zero) REAL WVLTARRAY(NSMPL) ! this is the WAVELET transformed waveform INTEGER WVLTDIM ! total number of samples ==> WVLT coef's REAL WVTHRES ! WAVELET threshold REAL WVPED ! WAVELET pedestal INTEGER WVREJ ! number of PEDESTAL points C local auxiliary variable INTEGER I C initialize returned variable WVREJ=0 ! i.e., no sample compressed DO I=1,WVLTDIM ! loop over the WAVELET dimension IF (ABS(WVLTARRAY(I)).LT.ABS(WVTHRES)) THEN WVLTARRAY(I)=WVPED ! set to PEDESTAL rejected samples WVREJ=WVREJ+1 ! update counter for rejected points END IF END DO ! close loop over the WAVELET dimension RETURN END +DECK,EVLTCOMP. SUBROUTINE EVLTCOMP(WVLTARRAY,WVLTDIM,WVPED,WVSAMPLES,PEDINF) C************************************************************************ C It accepts a wavelet array WVLTARRAY of dimension WVLTDIM * C and calculates the number of words needed in order to store it. * C We assume that when we find WVPED in the WVLTARRAY we store the * C start and the number of WVPED points following. * C This routine effectively evaluates compression factor. * C PEDINF is the number of words we assume we need to save for * C every pedestal region, e.g. START FLAG + PED LENGTH + PED VALUE * C * C @Erik Kats. Mar 1998 * C************************************************************************ IMPLICIT NONE INTEGER NSMPL ! that's the maximum number of samples PARAMETER (NSMPL=262144)! i.e., 262 microsecs... ASSUMING TICK=1ns C or for 5ns sampling 1.31 ms C Notice: all samples, start from 1 (NOT zero) REAL WVLTARRAY(NSMPL) ! this is the WAVELET transformed waveform INTEGER WVLTDIM ! total number of samples ==> WVLT coef's REAL WVPED ! WAVELET pedestal INTEGER WVSAMPLES ! number of SAMPLES needed to record frame INTEGER PEDINF ! number of WORDS needed for PED suppression C local auxiliary variable INTEGER I,XPED C initialize returned variable WVSAMPLES=0 ! i.e., no sample compressed C initialize PED region control variable XPED=0 ! i.e., not in PED region DO I=1,WVLTDIM ! loop over the WAVELET dimension IF (ABS(WVLTARRAY(I)).EQ.ABS(WVPED)) THEN C We have previously set WVLTARRAY(I)=WVPED for those points that C were less than a given threshold, so they should satisfy this condition. C Be cautious, although, of floating point arithmetic. IF (XPED.EQ.0) THEN ! first time we get in this PED region WVSAMPLES=WVSAMPLES+PEDINF ! increase sample counter... XPED=1 ! i.e., we're in PED END IF ! close check on PED region ELSE ! this is a NON-ped point... XPED=0 ! i.e., we're out of PED WVSAMPLES=WVSAMPLES+1 ! increase sample counter... END IF END DO ! close loop over the WAVELET dimension RETURN END +QUIT.