C C This source-code file contains (in alphabetical order by name) C the source code of the following 22 subroutines: C C burgdp.cl C bwlow.cl C bwring.cl C bwtdf.cl C dftftr.cl C errfor.cl C extend.cl C factr.cl C fejer.cl C hamm3.cl C hann3.cl C kais3.cl C lsqply.cl C rclow.cl C riesz.cl C trendr.cl C tukey.cl C unwin0.cl C varfqd.cl C window.cl C * winnow.cl C zmbes5.cl C C * - The entry subroutine. All others are called from or C through this subroutine. C C All subroutines were compiled and tested under: C C UNIX Fortran 77, HP Fortran/9000, HP-UX release 11.0 C C Some of them were compiled with the following flags. It is not C known whether these flags are required for successful C compilation or execution. The +E7 flags probably are NOT C required for any of the subroutines. The -K flags probably ARE C required. All 22 subroutines may be compiled with both the +E7 C and -K flags without any know problems: C C burgdp.f +E7 C bwlow.f +E7 -K C bwtdf.f +E7 -K C dftftr.f -K C rclow.f +E7 -K C unwin0.f -K C varfqd.f +E7 -K C winnow.f -K C C C -K Automatically SAVE all local variables in all C subprograms. This option forces static storage C for these variables in order to provide a C convenient path for importing FORTRAN 66 and C FORTRAN 77 programs that were written to depend on C static allocation of memory (that is, variables C retaining their values between invocations of C their containing program units). This option has C two side-effects: C C + All non-initialized variables are initialized C to zero. C + The DATA statement can appear among executable C statements. C C C +E7 Use static storage for arguments passed to C subroutines and functions. C C C C C C C________________________________________________________________ C C SUBROUTINE B U R G D P C________________________________________________________________ C C SUBROUTINE BURGDP FINDS THE DOUBLE-PRECISION COEFFICIENTS OF A C BURG ERROR-PREDICTION FILTER (UNIT SPAN). THE TIME SERIES MUST C BE OF DOUBLE-PRECISION DATA TYPE. IF PROCESSING COMPLEX DATA, C USE SUBROUTINE BURGDC. C C THE BURG ERROR-PREDICTION FILTER IS A REVERSABLE MINIMUM-PHASE C TIME-DOMAIN FILTER WITH A ZERO OUTPUT. IT'S PREDICTION C CHARACTERISTIC RESULTS FROM HOLDING THE FIRST COEFFICIENT TO C ONE DURING DERIVATION. ALL OF THE OTHER COEFFICIENTS OPERATE C ON THE DATA TO CANCEL THE FIRST POINT. C C THERE ARE THREE PRIMARY USES OF THE BURG FILTER. C C A FIRST USE IS TO EXTEND A SHORT DATA PROFILE WITH ARTIFICIAL C DATA OF A SIMILAR FREQUENCY SPECTRUM. TO DO THIS, EXTEND ONE C POINT AT A TIME ACCORDING TO: X(N+1) = -(X(N)*A(2)+X(N-1)*A(3)+ C X(N-2)*A(4)+ ... ). THEN SLIDE THE FILTER COEFFICIENTS ONE C POINT TO X(N+1) AND ITERATE TO THE DESIRED LENGTH. THIS C OPERATION MAY BE DONE IN BOTH DIRECTIONS. C C NOTE THAT THE BURG ALGORITHM DISCOVERS THE TRUE FREQUENCY C SPECTRUM SO ACCURATELY THAT IT CAN BE USED FOR EXTENDING C PROFILES TO MANY TIMES THEIR ORIGINAL LENGTH BEFORE FOURIER C TRANSFORMING. IN FACT, IF USED IN CONJUNCTION WITH WINDOWING, C EXTREMELY CLOSE SPECTRAL PEAKS CAN BE SEPARATED THAT OTHERWISE C WOULD HAVE BEEN BLURRED BY THE WINDOW OR IT'S SIDE LOBES. C C A SECOND USE IS TO FIND THE SPECTRAL CONTENT OF A DATA PROFILE C WITHOUT THE CLUTTERING END EFFECTS. THE FILTER IS DERIVED SUCH C THAT IT ASSUMES THE FUNCTION TO CONTINUE "IN A REASONABLE WAY" C INSTEAD OF TRUNCATING WITH ZEROS OR REPEATING THE INTERVAL. C THE SPECTRUM OF THE FUNCTION IS THE INVERSE OF THE FILTER C SPECTRUM. TO FIND THE FUNCTION SPECTRUM, EXTEND THE FILTER C COEFFICIENTS WITH ZEROS TO A REASONABLE LENGTH; THEN FOURIER C TRANSFORM THE EXTENDED FILTER COEFFICIENTS INTO THE FREQUENCY C DOMAIN; FINALLY, INVERT EACH FREQUENCY. TO PERFORM THE C FREQUENCY INVERSION AND PRODUCE A NATURAL-LOG POWER SPECTRUM IN C ONE EASY STEP, SIMPLY TAKE THE COMPLEX LOG AND MULTIPLY BY C NEGATIVE 2. (NOTE: THIS USE IS NOT EASILY INTERPRETED BECAUSE C THE SPECTRUM DOES NOT EVEN RESEMBLE A SPECTRUM FROM A FOURIER C TRANSFORM.) C C A THIRD USE IS TO PREDICT FUTURE VALUES IN A DATA STREAM. IF C THE DATA STREAM IS STATIONARY, ANOMALOUS VALUES WILL NOT BE C PREDICTED CORRECTLY AND CAN THEREFORE BE ISOLATED. IT MAY BE C REASONABLE TO AUGMENT THIS TECHNIQUE BY THE USE OF AN ADAPTIVE C FILTER SUCH AS THE WIDROW ALGORITHM (SEE CLAERBOUT CHAPTER 7-3, C PAGES 136-139). C C A COMPLETE DESCRIPTION OF THE BURG FILTER (INCLUDING A DETAILED C DERIVATION) MAY BE FOUND IN NOTES ENTITLED BURG SPECTRAL C ESTIMATION, 6NOV96. THE FOUNDATIONAL MATERIAL WAS OBTAINED C FROM "FUNDAMENTALS OF GEOPHYSICAL DATA PROCESSING" BY JON F. C CLAERBOUT (ISBN 0-86542-305-9), CHAPTER 7-2, PAGES 133-137. C C NOTE: TESTING HAS SHOWN THAT THE OPTIMUM FILTER LENGTH (NF) IS C ABOUT .707*ND. A SHORTER LENGTH CANNOT SUFFICIENTLY DETERMINE C ALL OF THE FREQUENCIES AND THEIR AMPLITUDES (ALTHOUGH FAIRLY C GOOD FREQUENCY DETERMINATION OCCURS DOWN TO ABOUT .250*ND; IF C SHORTER THAN THAT, IT RAPIDLY DETERIORATES). A GREATER LENGTH C BEGINS INTRODUCING MINOR NOISE. C C UNDER CERTAIN CIRCUMSTANCES THE BASIC BURG PREDICTION FILTER C CAN BECOME UNSTABLE. THIS OCCURS WHEN A BAND OF FREQUENCIES C (PARTICULARLY THE HIGH FREQUENCIES) IS MISSING FROM THE C SPECTRUM OF THE DATA TO BE FILTERED. TO AVERT THE POSSIBILITY C OF THIS HAPPENING, THIS SUBROUTINE ARTIFICIALLY INJECTS A C RANDOM (ALMOST WHITE) NOISE SIGNAL INTO THE COEFFICIENT C GENERATING PROCESS WITH AN AMPLITUDE OF -235 DB BELOW THE C VARIANCE OF THE PROFILE. THIS SIGNAL LEVEL IS SO SMALL THAT IT C WOULD PROBABLY NEVER BE SEEN IN SUBSEQUENT SIGNAL PROCESSING C (EVEN IF ONE IS LOOKING SPECIFICALLY FOR IT!). C C THE PROCESS THAT GENERATES THE RANDOM NOISE CAN BE MADE TO C GENERATE IDENTICAL NOISE PROFILES IN SUBSEQUENT SUBROUTINE C CALLS OR TO GENERATE DIFFERING NOISE PROFILES. TO CHANGE, THE C TWO COMMENTED STATEMENTS IN THE CODE MUST BE SWAPPED WITH THEIR C UNCOMMENTED COUNTERPARTS (IQX, ISEED AND RANDM4, RAN). THE C IDENTICAL NOISE PROFILES WILL RESULT IN IDENTICAL BURG C EXTENSIONS AND COEFFICIENTS FROM CALL TO CALL, BUT IF STACKING, C THEY COULD ADD CONSTRUCTIVELY. THE DIFFERING NOISE WILL CAUSE C SLIGHT VARIATIONS IN BURG EXTENSIONS IF THE ORIGINAL DATA C PROFILE HAD AN INCOMPLETE SPECTRUM (WAS UNSTABLE). ON BALANCE, C THE INJECTED NOISE IS SO SMALL THAT IT PROBABLY WOULD NEVER C STACK TO SIGNIFICANT LEVELS; THEREFORE, THE IDENTICAL NOISE C PROFILES ARE LIKELY THE BEST SELECTION. C C C INPUT ARGUMENTS: C ND - INTEGER*4. THE NUMBER OF POINTS IN XDAT. ND MUST C NOT EXCEED MXARR GIVEN AS A PARAMETER IN THE C SUBROUTINE. THE CALCULATION TIME WILL BE C PROPORTIONAL TO ND*NF. C XDAT - REAL*8. ARRAY OF DIMENSION ND CONTAINING THE DATA C FOR FILTER DERIVATION. IDEALLY, XDAT SHOULD BE SOME C KIND OF ERGODIC STATIONARY DATA IF EXTENSION IS C CONTEMPLATED FOR PREDICTIVE PURPOSES. C NF - INTEGER*4. THE DESIRED NUMBER OF FILTER C COEFFICIENTS. NF MUST NOT EXCEED ND. FOR A C MEANINGFUL FILTER, NF SHOULD NOT BE SMALLER THAN 2. C C OUTPUT ARGUMENT: C APF - REAL*8. ARRAY OF DIMENSION NF CONTAINING THE C TIME-DOMAIN COEFFICIENTS OF THE BURG ERROR- C PREDICTION FILTER DERIVED FROM THE DATA IN XDAT. C THE FIRST COEFFICIENT, APF(1), ALWAYS HAS A VALUE OF C ONE. C C SUBROUTINE BURGC GIVEN BY JON CLAERBOUT, STANFORD UNIVERSITY. C SUBROUTINE BURGDP ADAPTED BY ROB BRACKEN, USGS, 11NOV96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 11NOV96. C C subroutine burgdp(nd,xdat,nf,apf) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 nd real*8 xdat(0:*) integer*4 nf C C OUTPUT ARGUMENT real*8 apf(0:*) C C NUMBERS OF ELEMENTS integer*4 nx,mf C C INTERNAL ARRAYS integer*4 mxarr parameter(mxarr=200000) real*8 ep(0:mxarr),em(0:mxarr),bpf(0:mxarr) common /work1/ ep,em,bpf C C MISC VARIABLES real*8 ck,top,bot,epj C C NOISE INJECTION real*8 acoef(2),varf,varu,dbnoi,ampnoi parameter(dbnoi=-235.d0) integer*4 iseed c integer*4 iqx c real*4 randm4 C C CHECK ARRAY SIZES C if(nd.gt.mxarr+1) call errfor(mxarr, & '(burgdp) Array xdat too large') if(nf.gt.nd) call errfor(nf, & '(burgdp) Too many filter elements') C C INITIALIZE ARRAYS AND VARIABLES C C ADJUST UPPER LIMITS TO A ZERO TO N-1 SYSTEM (FROM ONE TO N) nx=nd-1 mf=nf-1 C C LOCATION OF CURRENT FILTER COEFS (ISW=0 APF, ISW=1 BPF) isw=0 C C ZEROTH FILTER COEFFICIENT apf(0)=1.d0 C C FIND AMPLITUDE OF STABILIZING RANDOM NOISE SIGNAL TO INJECT call trendr(1,2,acoef, 0,-nd,xdat(0)) call varfqd(0.d0,0,nd,xdat(0),acoef(2),acoef(1), varf,varu) ampnoi=dsqrt(varu*10.d0**(dbnoi/10.d0)) if(ampnoi.eq.0.d0) then C C PROF IS ALL 0; FLTR WILL BE UNSTBLE; MAKE ARTIFICIAL FLTR do i=1,mf apf(i)=-1.d0/dflotj(mf) enddo goto 990 endif C C PUT DATA VALUES INTO THE EP AND EM ARRAYS AND INJECT NOISE C (NOTE: USE OF FUNCTION RANDM4 AUTOMATICALLY SELECTS A C DIFFERENT SEED VALUE FOR EACH CALL TO ISEED(). THEREFORE, C IF SEVERAL PROFILES ARE EXTENDED USING THIS SUBROUTINE FOR C THE BURG FILTER COEFS AND THEN STACKED, THE NOISE WILL TEND C TO BE SUMMED OUT.) c iqx=iseed(0) iseed=15377317 do i=0,nx c ep(i)=xdat(i) c ep(i)=xdat(i)+ampnoi*dble(randm4()-0.5) ep(i)=xdat(i)+ampnoi*dble(ran(iseed)-0.5) em(i)=ep(i) enddo C C ITERATE FILTER COEFS OVER LENGTHS K=1,MF C do 701 k=1,mf C C FIND THE REFLECTION COEFFICIENT, CK top=0.d0 bot=0.d0 do j=k,nx jm=j-k top=top+em(jm)*ep(j) bot=bot+em(jm)*em(jm)+ep(j)*ep(j) enddo ck=2.d0*top/bot C C UPDATE EP AND EM do j=k,nx jm=j-k epj=ep(j) ep(j)=epj-ck*em(jm) em(jm)=em(jm)-ck*epj enddo C C FIND THE NEXT FILTER COEFS USING CK & LEVINSON RECURSION if(isw.eq.0) then C C CURRENT FILTER IN APF; PUT NEW IN BPF isw=1 apf(k)=0.d0 do i=0,k bpf(i)=apf(i)-ck*apf(k-i) enddo else C C CURRENT FILTER IN BPF; PUT NEW IN APF isw=0 bpf(k)=0.d0 do i=0,k apf(i)=bpf(i)-ck*bpf(k-i) enddo endif 701 continue C C MOVE NEW FILTER BACK INTO APF if(isw.eq.0) goto 990 do i=0,mf apf(i)=bpf(i) enddo C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE B W L O W C________________________________________________________________ C C SUBROUTINE BWLOW FILTERS A SIGNAL ARRAY USING A LOW-PASS C BUTTERWORTH TIME-DOMAIN FILTER WITH ANY NUMBER OF POLES. C C THE TRANSFER FUNCTION FOR EACH 2-POLE SECTION IS: C C H(F) = ( (ZZ)C0 + (Z)C1 + C2 ) / ( (ZZ)B0 + (Z)B1 + B2 ) C WHERE: C Z=EXP(2PIF), P=PI, I=SQRT(-1), F=FREQ IN RAD/SEC C C0=WW, C1=2WW, C2=WW C B0=WW+1-2WCOS(X), B1=2WW-2, B2=WW+1+2WCOS(X) C W=TAN(PFCDT), FC=CUTOFF FQ, DT=SAMPLE RATE C X=ANGLE OF THE UPPER POLE IN THE CONJUGATE POLE-PAIR C C THE RECURSIVE EQUATION FOR EACH 2-POLE SECTION IS: C C Y(N)B0 + Y(N-1)B1 + Y(N-2)B2 = C X(N)C0 + X(N-1)C1 + X(N-2)C2 C C WHICH BECOMES C C Y(N) = C ( X(N)C0 + X(N-1)C1 + X(N-2)C2 - Y(N-1)B1 - Y(N-2)B2 ) / B0 C C C INPUT ARGUMENTS: C FCUT - REAL*8. THE CUT-OFF FREQUENCY IN NORMALIZED UNITS C OF 1/SAMPLES. FCUT = FREQ (HZ) * SAMPLE-RATE (SEC). C NPOLE - INTEGER*4. THE NUMBER OF FILTER POLES IN THE LEFT C HALF OF THE COMPLEX S-PLANE. IF NPOLE IS NEGATIVE, C ALL POLES WILL BE MADE R-C WITH THE GIVEN FCUT. C (I.E. IT WILL BECOME AN R-C FILTER RATHER THAN A C BUTTERWORTH). C NREV - INTEGER*4. CAUSES SELECTED 2-POLE SECTIONS TO BE C FILTERED IN REVERSE FOR THE PURPOSE OF PHASE C COMPENSATION. EACH 2-POLE SECTION CONSISTS OF TWO C BUTTERWORTH POLES WHICH ARE COMPLEX CONJUGATES; ONE C ABOVE THE REAL AXIS AND THE OTHER REFLECTED BELOW. C THE SECTIONS ARE NUMBERED 1 THROUGH NPOLE/2 WITH THE C FIRST SECTION BEING FURTHEST FROM THE REAL AXIS. C THE POSSIBLE VALUES OF NREV ARE AS FOLLOWS: C 0 = NO SECTIONS FILTERED IN REVERSE. C 1 = ODD NUMBERED SECTIONS FILTERED IN REVERSE. C 2 = EVEN NUMBERED SECTIONS FILTERED IN REVERSE. C 3 = ALL SECTIONS FILTERED IN REVERSE. C WHEN NREV IS 1 OR 2, FILTERS WITH MULTIPLES OF C 4-POLES (NPOLE= 4, 8, 12, ...) WILL EXHIBIT NO PHASE C SHIFT. C NSAMP - INTEGER*4. THE NUMBER OF SAMPLES IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENTS: C YIO - REAL*8. AN ARRAY CONTAINING THE UNFILTERED INPUT C SIGNAL. UPON RETURN FROM THE SUBROUTINE, IT WILL C CONTAIN THE FILTERED OUTPUT SIGNAL. C C BUTTERWORTH EQUATIONS DERIVED FROM: C C HAMMING, R.W., 1977, DIGITAL FILTERS: PRENTICE-HALL, INC., C ENGLEWOOD CLIFFS, NEW JERSEY, PAGES 180-195. C C DAVID E. RYERSON, 1985, THE DESIGN AND USE OF DIGITAL C FILTERS WITH PERSONAL COMPUTER SPREADSHEET PROGRAMS, SANDIA C NATIONAL LABORATORIES, ALBUQUERQUE, NM, NTIS REPORT C DE85013805, SAND-85-1073, 55P. C C C SUBROUTINE BWLOW WRITTEN BY ROB BRACKEN, USGS, 23MAY91. C VAX/VMS VERSION 2.0 29MAY91. C HP-9000 SERIES 700/800 UNIX VERSION 2.0, 2JUN95. C subroutine bwlow(fcut,npole,nrev,nsamp,yio) C C DECLARATIONS C C PASSED ARGUMENTS real*8 fcut integer*4 npole,nrev,nsamp real*8 yio(*) C C OTHER VARIABLES c real*8 ax,ax0,ax1,ax2, ay1,ay2, b0,b0inv,b1,b2, c0,c1,c2 real*8 ax, ay1,ay2, b0,b0inv,b1,b2, c0,c1,c2 real*8 offset,polarc,polx, wc,wcsq,wcsq1,wctwo real*8 xin0,xin1,xin2 C C CONSTANTS real*8 pie parameter( pie = 3.1415 92653 58979 32384 62643 ) C C CHECK FOR NUMBERS OF POLES LESS THAN TWO C if(npole.eq.0) then goto 999 else if(npole.le.1) then mpole=jiabs(npole) krev=nrev goto 103 endif C C WARP THE CUTOFF FREQUENCY C wc=dtan(pie*fcut) wctwo=wc*2.d0 wcsq=wc*wc wcsq1=wcsq+1.d0 C C FIND COEFFICIENTS IN NUMERATOR OF TRANSFER FUNCTION C c0=wcsq c1=wcsq*2.d0 c2=wcsq C C FIND THE MIDDLE COEFFICIENT IN DENOMINATOR OF TRANSFER FUNCTION C b1=2.d0*(wcsq-1.d0) C C SET UP CONSTANTS FOR COMPUTING POLE PAIR LOCATIONS C polarc=pie/npole offset=(pie-polarc)/2 C C FILTER EACH 2-POLE SECTION C do 701 k=1,npole/2 C C COMPUTE THE OTHER COEFFICIENTS IN DENOMINATOR OF XFER FCN C polx=wctwo*dcos(polarc*k+offset) b0=(wcsq1-polx) b2=wcsq1+polx C C TIME-DOMAIN FILTER COEFFICIENTS FOR EACH BUTTERWORTH POLE-PAIR C b0inv=1.d0/b0 C AX0=C0*B0INV C AX1=C1*B0INV C AX2=C2*B0INV ax=c0*b0inv ay1=b1*b0inv ay2=b2*b0inv C C DETERMINE THE DIRECTION OF FILTERING C if(nrev.eq.0) then krev=0 else if(nrev.eq.3) then krev=1 else krev=k-(k/2)*2 if(nrev.eq.2) krev=1-krev endif C if(krev.eq.0) then C C SET INDICIES TO FILTER FORWARD FROM BEGINNING OF ARRAY i2=1 i1=2 i0=3 ie=nsamp istep=1 istep2=2 else C C SET INDICIES TO FILTER IN REVERSE FROM END OF ARRAY i2=nsamp i1=nsamp-1 i0=nsamp-2 ie=1 istep=-1 istep2=-2 endif C C FILTER THE STARTING SAMPLE POINT C if(nsamp.ge.1) then C KEEP UNFILTERED VALUE OF STARTING SAMPLE POINT xin2=yio(i2) C CALCULATE FILTERED VALUE OF STARTING SAMPLE POINT yio(i2)=ax*xin2 endif C C FILTER THE NEXT SAMPLE POINT C if(nsamp.ge.2) then C KEEP UNFILTERED VALUE OF NEXT SAMPLE POINT xin1=yio(i1) C CALCULATE FILTERED VALUE OF NEXT SAMPLE POINT yio(i1)= ax*(xin1 + 2.D0*xin2) 1 - ay1*yio(i2) endif C C FILTER THE REST OF THE SAMPLE POINTS C do 702 m=i0,ie,istep C KEEP UNFILTERED VALUE OF CURRENT SAMPLE POINT xin0=yio(m) C CALCULATE FILTERED VALUE OF CURRENT SAMPLE POINT yio(m)= ax*(xin0 + 2.d0*xin1 + xin2) 1 - ay1*yio(m-istep) - ay2*yio(m-istep2) C ROTATE ADDRESSES OF PREVIOUS UNFILTERED VALUES xin2=xin1 xin1=xin0 702 continue C 701 continue C C WHEN EVEN NUMBER OF POLES, SKIP THE R-C SECTION C if((npole/2)*2.eq.npole) goto 999 C C DETERMINE THE NUMBER OF POLES IN THE R-C SECTION C mpole=1 C C DETERMINE THE DIRECTION TO FILTER THE R-C SECTION C k=npole/2+1 if(nrev.eq.0) then krev=0 else if(nrev.eq.3) then krev=1 else krev=k-(k/2)*2 if(nrev.eq.2) krev=1-krev endif C C FILTER THE R-C SECTION C 103 call rclow(fcut,mpole,krev,nsamp,yio) C C EXIT PROCEDURE C 999 return end C C________________________________________________________________ C C SUBROUTINE B W R I N G C________________________________________________________________ C C SUBROUTINE BWRING CALCULATES THE NUMBER OF SAMPLES IN A TIME- C DOMAIN PROFILE NEEDED TO DAMP RINGING FROM A BUTTERWORTH C LOW-PASS FILTER TO A VALUE BELOW A GIVEN THRESHOLD. NUMBERS C MAY BE CALCULATED FOR EITHER UNSQUARED OR SQUARED TRANSFER C FUNCTIONS. (UNSQUARED ARE FOR A SINGLE PASS OF THE FILTER; C SQUARED ARE FOR A FORWARD AND REVERSE PASS WHICH ELIMINATES C PHASE SHIFTS). C C THE EQUATIONS USED FOR THIS CALCULATION HAVE BEEN DERIVED C EMPIRICALLY FROM SAMPLE DATA USING A TIME-DOMAIN SQUARED C BUTTERWORTH FILTER. CALCULATIONS FOR A STANDARD BUTTERWORTH C HAVE BEEN ESTIMATED FROM THE SQUARED DATA. C C INPUT ARGUMENTS: C THRESH - REAL*8. THE THRESHOLD BELOW WHICH RINGING MAY BE C TOLERATED. A TYPICAL VALUE OF THRESH MIGHT BE C 1/4096 WHICH WOULD BE AN ABSOLUTE POWER OF -72 DB. C THRESH MUST NOT BE ZERO. C STEP - REAL*8. THE SIZE OF THE STEP CAUSING THE RINGING. C THIS VALUE IS THE EXPECTED PEAK VALUE OF THE WHOLE C PROFILE CALCULATED AS THE SQUARE ROOT OF TWO TIMES C THE PROFILE VARIANCE -> SQRT(2*VAR) = SQRT(2)*SIGMA. C FQC - REAL*8. THE CUTOFF FREQUENCY (-3 DB FOR UNSQUARED C AND -6 DB FOR SQUARED) WHERE 0.5 IS NYQUIST. FQC C MUST NOT BE ZERO. C NPOLE - INTEGER*4. THE NUMBER OF POLES IN THE UNSQUARED C FILTER. FOR AN UNSQUARED BUTTERWORTH, THE ROLL-OFF C IS ABOUT -6 DB PER OCTIVE PER POLE; FOR A SQUARED C BUTTERWORTH, ABOUT -12 DB. C KSQ - INTEGER*4. IF KSQ IS 0, THE NUMBERS FOR THE C UNSQUARED FILTER WILL BE CALCULATED. IF KSQ IS 1, C THE SQUARED FILTER NUMBERS WILL BE FOUND. C C OUTPUT ARGUMENTS: C JLEFT - INTEGER*4. THE NUMBER OF DATA SAMPLES FROM THE LEFT C TO WHERE THE RINGING DAMPS BELOW THE SIZE OF THE C THRESHOLD. C KRIGHT - INTEGER*4. THE NUMBER OF DATA SAMPLES FROM THE C RIGHT TO WHERE THE RINGING DAMPS BELOW THE SIZE OF C THE THRESHOLD. C C INPUT/OUTPUT ARGUMENTS: C C SUBROUTINE BWRING WRITTEN BY ROB BRACKEN, USGS. C HP-9000 SERIES 700/800 UNIX FORTRAN 77. C SUBROUTINE BWRING VERSION 1.0, 19970204. C C subroutine bwring(thresh,step, fqc,npole,ksq, jleft,kright) C C DECLARATIONS C C INPUT ARGUMENTS real*8 thresh,step,fqc integer*4 npole,ksq C C OUTPUT ARGUMENTS integer*4 jleft,kright C C INTERNAL VARIABLES real*8 tos,fqc2,pole,ipl,ipr C C CHECK THRESHOLD AND STEP C if(thresh.eq.0.d0) & call errfor(0,'(bwring) thresh = 0') tos=1.d0 if(step.ne.0.d0) then tos=dabs(thresh/step) if(tos.gt.1.d0) tos=1.d0 endif C C CHECK LIMITS OF FREQUENCY C if(fqc.eq.0.d0) & call errfor(0,'(bwring) fqc = 0') fqc2=dabs(fqc) if(fqc2.gt.0.5d0) fqc2=0.5d0 C C FIND WHETHER TO SQUARE THE XFER FUNCTION C pole=dflotj(jiabs(npole)) if(ksq.eq.0) pole=pole/2.d0 C C FIND DAMPING VALUES C jl=jidnnt((pole**1.4d0)/(10.58d0*fqc2)) kright=jidnnt((-pole*dlog10(tos))/(4.23d0*fqc2)) if(kright.lt.jl) kright=jl jleft=kright-jl C if(tos.lt.1.d0) then ipl=jidnnt(1.d0/(2.d0*fqc2)) ipr=jidnnt(pole/(9.52d0*fqc2))+ipl if(jleft.lt.ipl) jleft=ipl if(kright.lt.ipr) kright=ipr endif C C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE B W T D F C________________________________________________________________ C C SUBROUTINE BWTDF LOWPASS FILTERS A PROFILE USING A TIME-DOMAIN C BUTTERWORTH FILTER. OPTIONALLY, THE FILTER MAY BE SET TO MAKE C A FORWARD PASS (UNSQUARED) OR A FORWARD AND REVERSE PASS C (SQUARED). IF THE FILTER IS SQUARED, THE TRANSFER FUNCTION IS C SQUARED AS WELL. BWTDF PERFORMS AUTOMATICALLY A NUMBER OF DATA C PREPARATIONS WHICH REDUCE END EFFECTS (RINGING). THESE STEPS C INCLUDE REMOVAL OF A LEAST SQUARES LINE, EXTENSION WITH BURG C PREDICTION, AND APPLICATION OF A SPLIT RIESZ WINDOW TO TAPER C THE EXTENSIONS TO ZERO. UPON RETURN, THE EXTENSIONS ARE C REMOVED AND THE LSQ LINE IS ADDED BACK. C C FOLLOWING IS A BRIEF DESCRIPTION OF THE BUTTERWORTH FILTER C CALLED BY THIS SUBROUTINE: C C THE TRANSFER FUNCTION OF THE BUTTERWORTH FILTER HAS POLE PAIRS C ON THE LEFT COMPLEX PLANE STARTING NEAR THE POSITIVE AND C NEGATIVE IMAGINARY AXIS AND PROGRESSING TOWARD THE LEFT REAL C AXIS. IF THE NUMBER OF POLES IS ODD, THE FINAL (ODD) POLE IS C ON THE LEFT REAL AXIS AND BECOMES A SIMPLE RC POLE. THE TIME C DOMAIN OPERATION INVOLVES RUNNING THE ENTIRE LENGTH OF THE C PROFILE WITH EACH 2-POLE SECTION (POLE PAIR) STARTING NEAR THE C IMAGINARY AXIS AND PROGRESSING TO THE LEFT REAL AXIS JUST AS IN C THE TRANSFER FUNCTION. C C THE TIME-DOMAIN BUTTERWORTH FILTER HAS A SIGNAL DELAY RESULTING C FROM VARIOUS PHASE SHIFTS AT VARIOUS FREQUENCIES. THESE DELAYS C SEEM TO SETTLE OUT TO A CONSTANT TIME DELAY BELOW CUTOFF C FREQUENCY WHICH MAY BE DESCRIBED (EMPIRICALLY) AS 37.071 C DEGREES PER POLE AT CUTOFF. NOTE THAT FOR PURPOSES OF C CALCULATING DELAY, THE CUTOFF FREQUENCY IS WARPED INTO THE C INTERVAL (ZERO,INFINITY) BY TAN(PI*FCUT)/PI. THE DELAY MAY BE C ACCOUNTED FOR BY APPLYING AN APPROPRIATE CORRECTION BASED ON C THIS INFORMATION. ALTERNATIVELY, THE PHASE SHIFTS (AND HENCE C THE DELAY) MAY BE ELIMINATED ENTIRELY BY RUNNING THE IDENTICAL C FILTER BOTH DIRECTIONS, THUS SQUARING THE TRANSFER FUNCTION. C THE PRIMARY EFFECTS OF THIS OPERATION ARE: REMOVE PHASE C SHIFTS, DOUBLE THE EFFECTIVE NUMBER OF POLES (DOUBLE THE SLOPE C IN THE STOP BAND), AND DECREASE THE POWER AT CUT-OFF TO -6 DB. C C THE TRANSFER FUNCTION FOR EACH 2-POLE SECTION IS: C C H(F) = ( (ZZ)C0 + (Z)C1 + C2 ) / ( (ZZ)B0 + (Z)B1 + B2 ) C WHERE: C Z=EXP(2PIF), P=PI, I=SQRT(-1), F=FREQ IN RAD/SEC C C0=WW, C1=2WW, C2=WW C B0=WW+1-2WCOS(X), B1=2WW-2, B2=WW+1+2WCOS(X) C W=TAN(PFCDT), FC=CUTOFF FQ, DT=SAMPLE RATE C X=ANGLE OF THE UPPER POLE IN THE CONJUGATE POLE-PAIR C C THE RECURSIVE EQUATION FOR EACH 2-POLE SECTION IS: C C Y(N)B0 + Y(N-1)B1 + Y(N-2)B2 = C X(N)C0 + X(N-1)C1 + X(N-2)C2 C C WHICH BECOMES C C Y(N) = C ( X(N)C0 + X(N-1)C1 + X(N-2)C2 - Y(N-1)B1 - Y(N-2)B2 ) / B0 C C BUTTERWORTH EQUATIONS DERIVED FROM: C C HAMMING, R.W., 1977, DIGITAL FILTERS: PRENTICE-HALL, INC., C ENGLEWOOD CLIFFS, NEW JERSEY, PAGES 180-195. C C DAVID E. RYERSON, 1985, THE DESIGN AND USE OF DIGITAL C FILTERS WITH PERSONAL COMPUTER SPREADSHEET PROGRAMS, SANDIA C NATIONAL LABORATORIES, ALBUQUERQUE, NM, NTIS REPORT C DE85013805, SAND-85-1073, 55P. C C C INPUT ARGUMENTS: C FCUT - REAL*8. THE CUTOFF FREQUENCY IN UNITS OF 1/SAMPLES C AND THE PASSBAND FREQUENCY. THE MAXIMUM FREQUENCY C IS NYQUIST (FCUT=0.5D0); THE MINIMUM FREQUENCY MAY C APPROACH D.C. (FCUT=0.D0) IN THE LIMIT BUT SHOULD C NEVER BE D.C. (AS THE FREQUENCY DECREASES, THE C PROCESSING TIME INCREASES DUE TO AN INCREASED C EXTENSION LENGTH WHICH IS CALCULATED BY AN N*N C ALGORITHM). FCUT IS ALSO THE PASSBAND FREQUENCY C USED FOR AN INTERNAL CALCULATION OF THE NUMBER OF C FILTER POLES. THE POWER OF FREQUENCIES AT FCUT WILL C BE REDUCED BY -3 DB (AMP=1/SQRT(2)) FOR THE C UNSQUARED FILTER AND -6 DB (AMP=1/2) IF SQUARED. AT C FREQUENCIES BELOW FCUT, THE POWER REDUCTION C APPROACHES 0 DB IN A MAXIMALLY FLAT CURVE AS IS C CHARACTERISTIC OF A BUTTERWORTH FILTER. AT C FREQUENCIES ABOVE FCUT, THE POWER REDUCES AT A RATE C OF ABOUT 6 DB PER OCTIVE PER POLE FOR THE UNSQUARED C FILTER AND 12 DB IF SQUARED. THE ACTUAL RATE OF C REDUCTION CAN BE SET USING FSTOP. C FSTOP - REAL*8. THE STOPBAND FREQUENCY IN THE UNITS OF C 1/SAMPLES (SAME AS FCUT). FSTOP IS THE FREQUENCY AT C WHICH THE POWER HAS BEEN REDUCED BY -36 DB C (AMP=1/SQRT(4097)) FOR THE UNSQUARED FILTER AND -72 C DB (AMP=1/4097) IF SQUARED. BECAUSE FSTOP IS USED C ONLY FOR AN INTERNAL CALCULATION OF THE NUMBER OF C POLES, NYQUIST FREQUENCY DOES NOT POSE A LIMITATION. C THEREFORE FSTOP MAY BE ANY VALUE ABOVE FCUT WHICH C PRODUCES THE DESIRED ROLLOFF. HOWEVER, IT SHOULD BE C STRESSED THAT A VERY STEEP ROLLOFF MAY CAUSE RINGING C WITHIN THE PROFILE FROM NON-END-EFFECT STEPS AND MAY C INCREASE COMPUTATION TIME. SOME RINGING MAY BE C NOTICED ABOVE 8 POLES AND IT MAY BECOME SERIOUS C ABOVE 16 POLES. THE NUMBER OF POLES MAY BE FOUND C BY: C C NPOLE=LOG(64)/LOG(FSTOP/FCUT) C C NOTE THAT THIS EQUATION HOLDS REGARDLESS OF WHETHER C THE FILTER IS SQUARED BECAUSE THE PASSBAND AND C STOPBAND POWERS HAVE BEEN REDEFINED FOR THE SQUARED C FILTER. THE UNADAPTED EQUATION FOR THE NUMBER OF C POLES IS: C C N = LOG( SQRT(1/PP-1)/SQRT(1/PS-1) )/LOG( WP/WS ) C C WHERE: C N = NUMBER OF POLES C PP = POWER AT PASSBAND FREQUENCY, WP C PS = POWER AT STOPBAND FREQUENCY, WS C WP = PASSBAND FREQUENCY FREQUENCY C WS = STOPBAND FREQUENCY C C KSQ - INTEGER*4. DETERMINES WHETHER THE FILTER SHOULD BE C SQUARED BY A SECOND REVERSE PASS. IF KSQ=0, NO C SQUARING WILL OCCUR. IF KSQ=1, SQUARING WILL BE C PERFORMED. IF SQUARING IS SELECTED, THE POWER WILL C BE -6 DB AT FCUT AND -72 DB AT FSTOP; THERE WILL BE C NO PHASE-SHIFTS NOR SIGNAL DELAYS. THIS IS THE C RECOMMENDED OPTION. HOWEVER, IF SQUARING IS NOT C SELECTED, THE POWER WILL BE -3 DB AT FCUT AND -36 DB C AT FSTOP; THERE WILL BE PHASE SHIFTS AND SIGNAL C DELAYS WITH SOME OF THE PROFILE BEING WASHED OFF THE C RIGHT END. AT FREQUENCIES BELOW FCUT A DELAY WILL C OCCUR WHICH CAN BE DESCRIBED AS 37.071 DEGREES PER C POLE AT A CUTOFF FREQUENCY WARPED INTO THE INTERVAL C (ZERO,INFINITY) BY TAN(PI*FCUT)/PI. C RINGDB - REAL*8. THE AMOUNT OF RINGING FROM END EFFECTS THAT C CAN BE TOLERATED IN THE PROFILE. RINGDB IS A C NEGATIVE POWER EXPRESSED IN DECIBELS RELATIVE TO C TWICE THE VARIANCE OF THE PROFILE (TWICE THE C VARIANCE IS CONSIDERED 0 DB). A REASONABLE VALUE C FOR RINGDB WOULD BE -72 DB. THIS MEANS THAT AFTER C FILTERING, RINGING AT EITHER END OF THE PROFILE WILL C HAVE AMPLITUDE ON THE ORDER OF SQRT(2)*STDDEV/4096. C IF DESIRED, BY MAKING RINGDB A POSITIVE VALUE, THE C RINGING CAN BE SPECIFIED ABSOLUTELY WHERE 0 DB MEANS C A NORMALIZED VARIANCE OF 1/2 (RATHER THE RELATIVE C VALUE AS DESCRIBED ABOVE WHERE 0 DB IS TWICE THE C PROFILE VARIANCE). IN THIS CASE, RINGDB=+72 WOULD C MEAN RINGING WITH AMPLITUDE ON THE ORDER OF 1/4096. C THE RINGING WILL BE AS SPECIFIED BY RINGDB ONLY AT C THE ENDS OF THE PROFILE; IT WILL DECREASE TOWARD THE C CENTER. C NPTS - INTEGER*4. THE NUMBER OF DATA POINTS IN YFNC. NPTS C MUST NOT EXCEED NDIM UNLESS NDIM = 0. ALSO, IF NPTS C IS 2 OR SMALLER, BWTDF WILL NOT DO ANYTHING; C FILTERING IN THAT CONTEXT IS NOT MEANINGFUL. C C INPUT/OUTPUT ARGUMENTS: C YFNC - REAL*8. ARRAY OF DIMENSION NDIM CONTAINING THE C PROFILE OF NPTS DATA POINTS TO BE FILTERED. THE C EXTENSION BEYOND NPTS TO NDIM WILL BE USED ONLY AS A C WORK SPACE AND UPON RETURN WILL CONTAIN NO USEFUL C INFORMATION. C NDIM - INTEGER*4. THE DIMENSION OF YFNC. YFNC MUST BE C LARGER THAN THE NUMBER OF DATA POINTS TO ACCOMODATE C PROFILE EXTENSIONS. IF, DURING COMPUTATIONS, THE C PROFILE MUST BE EXTENDED FARTHER THAN NDIM, AN ERROR C WILL BE GENERATED. HOWEVER, IF NDIM IS GIVEN A C VALUE OF ZERO DURING THE CALL, THE MAXIMUM NEEDED C DIMENSION WILL BE CALCULATED AND RETURNED IN NDIM; C THE PROFILE WILL NOT BE FILTERED. C C SUBROUTINE BWTDF WRITTEN BY ROB BRACKEN, USGS, 19970604. C HP9000 SERIES 700/800 VERSION 1.0, 19970604. C C subroutine bwtdf(fcut,fstop,ksq,ringdb,npts, yfnc,ndim) C C DECLARATIONS C C INPUT ARGUMENTS real*8 fcut,fstop,ringdb integer*4 ksq,npts C C INPUT/OUTPUT ARGUMENT real*8 yfnc(*) C C PASS AND STOP-BAND POWERS real*8 ppdb,psdb parameter(ppdb=-3.0103d0,psdb=-36.125d0) real*8 pp,ps,ppx,psx integer*4 npole C C LEAST SQUARES LINE REMOVAL integer*4 ktype,ncoef,nverse real*8 acoef(2) C C FREQUENCY DEPENDENT VARIANCE integer*4 npolev real*8 slope,rcept real*8 varf,varu C C EXTENSION-AMOUNT CALCULATION real*8 rdb,profdb,thresh,step integer*4 jleft,kright C C EXTENSION PROCEDURES real*8 tparms(3) integer*4 just,nout1,nout2 C C WINDOWING FUNCTION real*8 alpha,wfrac integer*4 kdft C C BUTTERWORTH FILTER integer*4 nrev C C CHECK FREQUENCIES C if(fcut.le.0.d0) call errfor(0, & '(bwtdf) fcut is d.c. or negative') if(fcut.gt.0.5d0) call errfor(0, & '(bwtdf) fcut is above nyquist (0.5)') C C FIND NUMBER OF POLES C if(fstop.gt.fcut) then pp=10.d0**(ppdb/10.d0) ppx=dsqrt(1.d0/pp-1.d0) ps=10.d0**(psdb/10.d0) psx=dsqrt(1.d0/ps-1.d0) npole=jidnnt(dlog10(ppx/psx)/dlog10(fcut/fstop)) if(npole.lt.1) npole=1 if(npole.gt.64) npole=64 else npole=8 endif C C CHECK THE NUMBER OF POINTS C if(npts.le.2) then if(ndim.le.0) ndim=npts goto 990 endif C C FIND COEFS OF A LEAST SQUARES LINE LATER TO BE REMOVED C if(ndim.ge.1.or.ringdb.ge.0.d0) then ktype=1 ncoef=2 nverse=0 call trendr(ktype,ncoef,acoef,nverse,-npts,yfnc) endif C C FIND THE FREQ DEPENDENT VARIANCE (PWR IN PROF BELOW CUTOFF FQ) C if(ringdb.lt.0.d0) then rdb=ringdb C C RINGDB IS A RELATIVE THRESHOLD; EXTENSION LENGTHS WILL BE C CALCULATED BY BWRING SUCH THAT POWER IS RELATIVE TO THE C PROFILE VARIANCE. ADJUSTMENT OF THE REQUESTED RINGING C POWER IS UNNECESSARY. profdb=0.d0 else rdb=-ringdb C C RINGDB IS AN ABSOLUTE THRESHOLD; BECAUSE EXTENSION C LENGTHS ARE CALCULATED BY BWRING RELATIVE TO PROFILE C VARIANCE, THE EFFECT OF VARIANCE MUST BE REMOVED FROM THE C REQUESTED RINGING POWER IN ORDER TO TIE IT TO THE C ABSOLUTE STANDARD (VARIANCE = 1/2). slope=acoef(2) rcept=acoef(1) npolev=4 call varfqd(fcut,npolev,npts,yfnc,slope,rcept,varf,varu) profdb=10.d0*dlog10(2.d0*varf) endif C C DETERMINE AN AMOUNT TO EXTEND THE PROFILE C C FIND A THE RINGING THRESHOLD thresh=rdb-profdb C C CONVERT THRESH TO AN AMPLITUDE FOR SUBR BWRING thresh=10.d0**(thresh/20.d0) step=1.d0 call bwring(thresh,step, fcut,npole,ksq, jleft,kright) C C FIND THE LENGTHS AFTER THE TWO NEEDED EXTENSIONS nout1=jleft+npts+jleft nout2=jleft+npts+kright C C CHECK NDIM if(ndim.le.0) then C C CALCULATE THE VALUE OF NDIM AND RETURN TO CALLING ROUTINE ndim=nout2 goto 990 endif C C REMOVE A LSQ LINE USING THE PREVIOUSLY CALCULATED COEFS C ktype=1 ncoef=2 nverse=0 call trendr(ktype,-ncoef,acoef,nverse,npts,yfnc) C C EXTEND THE PROFILE (WITH A BURG PREDICTION FILTER) C C NOTE: THE MOST ACCURATE EXTENSION CAN BE OBTAINED BY USING THE C LENGTH-OF-PROFILE DEPTH OF COEFS RATHER THAN THE C DEPTH-OF-EXTENSION. HOWEVER, THE TRADE-OFF IS THAT THE GREATER C DEPTH REQUIRES MORE COMPUTATION TIME. C C BE SURE THE EXTENSIONS WILL FIT IN THE ALOTTED DIMENSION if(ndim.lt.nout2) & call errfor(nout2, & '(bwtdf) dimension ndim of array yfnc is too small') C C SELECT BURG PREDICTION FILTER FOR EXTENSION ktype=7 C C DEPTH OF BURG FILTER COEFS = DEPTH OF EXTENSION tparms(1)=-1.d0 C DEPTH OF BURG FILTER COEFS = LENGTH OF PROFILE c tparms(1)=0.d0 C DEPTH OF BURG FILTER COEFS = VALUE OF TPARMS(2) AND (3) c tparms(1)=-3.d0 c tparms(2)=2048.d0 c tparms(3)=2048.d0 just=0 call extend(ktype,tparms,just,npts,nout1,yfnc) C C WINDOW THE PROFILE WITH A SPLIT RIESZ (ALPHA=2.5) WINDOW FCN C C NOTE: DO NOT ATTEMPT WINDOWING WITH A FULL LENGTH WINDOW C BECAUSE THE NECESSARY SUBSEQUENT UNWINDOWING WILL RE-AMPLIFY C THE RINGING (AND MAY ADD PHASE SHIFTS) MAKING IT WORSE THAN IF C NO WINDOW WAS USED! C C NOTE2: THE SPLIT RIESZ WITH ALPHA=2.5 AND WFRAC SUCH THAT THE C WINDOWING TAPER IS COMPLETED AT HALF EACH EXTENSION IS FOUND TO C BE THE BEST; BUT IT IS INCONSISTANT. SOMETIMES IT DOES NOT C IMPROVE THE RINGING AT ALL OVER AN UNWINDOWED PROFILE. C ktype=4 alpha=2.5d0 wfrac=dflotj(jleft)/dflotj(nout1) kdft=0 nverse=0 call window(ktype,alpha,wfrac,kdft,nverse,nout1,yfnc) C C EXTND RIGHT END FURTHER TO MAKE ROOM FOR WAVES BEING PUSHED OFF C ktype=0 tparms=0.d0 just=1 call extend(ktype,tparms,just,nout1,nout2,yfnc) C C PERFORM BUTTERWORTH LOW-PASS TIME-DOMAIN FILTERING C nrev=0 call bwlow(fcut,npole,nrev,nout2,yfnc) if(ksq.ge.1) then nrev=3 call bwlow(fcut,npole,nrev,nout2,yfnc) endif C C C REMOVE EXTENSIONS C ktype=0 tparms=0.d0 just=1 call extend(ktype,tparms,just,nout2,nout1,yfnc) just=0 call extend(ktype,tparms,just,nout1,npts,yfnc) C C REPLACE LSQ LINE C ktype=1 ncoef=2 nverse=1 call trendr(ktype,ncoef,acoef,nverse,npts,yfnc) C C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE D F T F T R C________________________________________________________________ C C SUBROUTINE DFTFTR PERFORMS A DISCRETE FOURIER TRANSFORM ON ONE- C DIMENSIONAL DATA. THE TIME DOMAIN IS DOUBLE-PRECISION REAL AND C THE FREQUENCY DOMAIN IS DOUBLE-COMPLEX. THE TIME AND FREQUENCY C DOMAINS ARE ALLOWED TO DIFFER IN LENGTH. DFTFTR IS SIMILAR TO C DFTFDR BUT HAS AN ADDITIONAL FEATURE THAT ALLOWS THE TIME C DOMAIN TO BE TRUNCATED IMPLICITLY AT THE ENDS THUS GIVING A C FASTER TRANSFORM WHEN TRUNCATION IS CONTEMPLATED. C C DIFFERING LENGTH DOMAINS ARE INTERPRETED BY DFTFTR AS REQUIRING C A HIGH FREQUENCY TRUNCATION OR EXTENSION WITH ZEROS (SEE C INTERPRETATION 2, BELOW). UPON FORWARD AND INVERSE C TRANSFORMING HOWEVER, RESTORATION OF THE ORIGINAL TIME-DOMAIN C PROFILE CAN ONLY BE ASSURED IF THE FREQUENCY DOMAIN HAS THE C SAME OR MORE POINTS THAN THE TIME DOMAIN. C C IF THE NUMBERS OF POINTS IN THE TIME AND FREQUENCIES DOMAINS C ARE EQUAL, THIS SUBROUTINE PRODUCES A TRUE DFT (SIGNI=-1 FOR C FORWARD AND SIGNI=+1 FOR INVERSE). THE SCALE FACTORS ARE 1/NRX C FOR FORWARD TRANSFORM AND 1 FOR INVERSE TRANSFORM. C C THE FOLLOWING TABLE SHOWS TIMES REQUIRED FOR VARIOUS ROUTINES C TO PERFORM A FORWARD AND INVERSE TRANSFORM WHEN THE NUMBERS OF C TIME AND FREQUENCY DOMAIN POINTS ARE EQUAL. (TIMES TAKEN ON C MUSETTE, HP-9000 SERIES 700/800 COMPUTER, ON 14MAY97. THIS C COMPUTER RUNS ROUGHLY 30 TIMES FASTER THAN A VAX 11/750.) C C ROUTINE TYPE RULE LENGTH TIME (SEC) C C FORKDC FFT N*LOGN 4096 0.80 C DFTTDR DFT N*N 4096 31.69 C DFTFDR DFT N*N 4096 31.69 C DFTTDC DFT N*N 4096 66.98 C C (THE DFT IS FASTER THAN THE FFT FOR N < 64) C (THE BREAK-EVEN POINT IS ACTUALLY N = 48 BUT THE FFT C CANNOT ACCEPT A 48 POINT PROFILE) C C C C C A DFT CAN INTERPRET DIFFERING LENGTH DOMAINS IN TWO WAYS: 1) C THE BEGINNING AND END OF THE TIME-DOMAIN ARE TO BE TRUNCATED OR C EXTENDED WITH ZEROS, 2) THE HIGH FREQUENCIES IN THE FREQUENCY C DOMAIN ARE TO BE TRUNCATED OR EXTENDED WITH ZEROS C C 1) TIME-DOMAIN TRUNCATION/EXTENSION IS THE SAME AS CHANGING THE C DENSITY OF FREQUENCIES WITHIN THE FREQUENCY DOMAIN. THE C PROCEDURE IS DONE BY SUBROUTINE DFTTDR. BY EXTENDING THE TIME C DOMAIN WITH ZEROS, THE NUMBER OF POINTS IN THE FREQUENCY DOMAIN C INCREASES BUT THE OVERALL SHAPE OF THE FREQUENCY-DOMAIN PROFILE C REMAINS THE SAME. SIMILARLY WITH TRUNCATION. (NOTE THAT C TRUNCATION OF THE TIME DOMAIN IS SO SIMPLE THAT IT SHOULD BE C DONE IN THE SUBROUTINE CALL. IF IN THE FORWARD CALL THE C TIME-DOMAIN HAS MORE POINTS THAN THE FREQUENCY DOMAIN, DFTTDR C WILL PRODUCE A DEGENERATE FREQUENCY DOMAIN RATHER THAN THE C FREQUENCY DOMAIN OF A TRUNCATED TIME DOMAIN PROFILE.) C C 2) FREQUENCY-DOMAIN TRUNCATION/EXTENSION IS THE SAME AS C CHANGING THE DENSITY OF SAMPLES WITHIN THE TIME DOMAIN. THIS C PROCEDURE IS DONE BY THIS SUBROUTINE, DFTFTR. BY EXTENDING THE C HIGH FREQUENCIES IN THE FREQUENCY DOMAIN WITH ZEROS, THE NUMBER C OF POINTS IN THE TIME DOMAIN (AFTER INVERSE TRANSFORM) C INCREASES BUT THE OVERALL SHAPE OF THE TIME-DOMAIN PROFILE C REMAINS THE SAME. SIMILARLY WITH TRUNCATION. (TRUNCATION OF C THE FREQUENCY DOMAIN IS NOT AS SIMPLE AS IN THE TIME DOMAIN; C THEREFORE DFTFTR DOES NOT ALLOW THE POSSIBILITY OF A DEGENERATE C TIME DOMAIN. RATHER, IF ON INVERSE TRANSFORM THE NUMBER OF C TIME-DOMAIN POINTS IS LESS THAN THE NUMBER OF FREQUENCY-DOMAIN C POINTS, DFTFTR PRODUCES THE FREQUENCY TRUNCATION INTERNALLY.) C C C C C DFTFTR PROVIDES THE OPTION OF PRODUCING A NEW TIME DOMAIN WHICH C HAS BEEN ACCORDIONED-IN OR -OUT TO HAVE FEWER OR MORE POINTS C THAN THE ORIGINAL TIME DOMAIN. THIS IS DONE IMPLICITLY BY C TRUNCATING OR EXTENDING THE HIGH FREQUENCIES IN THE FREQUENCY C DOMAIN. THE IMPLICIT TRUNCATION/EXTENSION IS FASTER THAN USING C A STANDARD DFT TO TRANSFORM INTO THE FREQUENCY DOMAIN AND THEN C MANUALLY TRUNCATE OR EXTEND THE FREQUENCIES. TO ACCORDION A C TIME DOMAIN, THE RECOMMENDED PROCEDURE IS TO FORWARD TRANSFORM C THE ORIGINAL TIME DOMAIN INTO A FREQUENCY DOMAIN WITH EQUAL OR C FEWER POINTS; THEN INVERSE TRANSFORM INTO THE NEW TIME DOMAIN C WITH EQUAL OR MORE POINTS THAN THE FREQUENCY DOMAIN. THE C FREQUENCY DOMAIN PRODUCED BY DFTFTR IS NOT PARTICULARLY C INTERESTING BECAUSE IT SIMPLY HAS THE HIGH FREQUENCIES C TRUNCATED OR EXTENDED WITH ZEROS. C C THE FOLLOWING TABLE GIVES THE EFFECTS DFTFTR HAS ON VARIOUS C EXAMPLE-LENGTH TIME AND FREQUENCY DOMAIN PROFILES: C C TD1 -1 FD +1 TD2 DESCRIPTION C NRX NCY NRX2 C C 100 100 PT ORIGINAL TIME DOMAIN (TD1) C C 100 100 100 PT FREQ DOMAIN (FD) NORMAL DFT-STYLE C 100 100 100 100 PT NEW TIME DOM (TD2) IDENTICAL TO TD1 C 100 100 200 TD2 = TD1 ACCORDIONED TO 200 PTS (FASTER) C C 100 200 100 PT FD WITH 100 ADDITIONAL ZERO HF PTS C 100 200 100 TD2 = TD1 C 100 200 200 TD2 = TD1 ACCORDIONED TO 200 PTS (SLOWER) C C 200 200 PT TD1 C C 200 100 200 PT FD WITH 100 HF POINTS TRUNCATED C 200 100 100 TD2 = TD1 ACCORDIONED TO 100 PTS (FASTER) C 200 100 200 TD2 = TD1 WITH HIGH FREQUENCIES MISSING C C 200 200 200 PT FD NORMAL DFT-STYLE C 200 200 100 TD2 = TD1 ACCORDIONED TO 100 PTS (SLOWER) C 200 200 200 TD2 = TD1 C C C AN ADDITIONAL FUNCTION OF DFTFTR IS THAT IT ALLOWS THE TIME C DOMAIN TO BE TRUNCATED IMPLICITLY USING THE TWO ARGUMENTS JB C AND JE. THIS TRUNCATION IS INTENDED PURELY AS A TIME SAVINGS C WHICH IS USEFUL ONLY IF THE TIME DOMAIN IS TO BE TRUNCATED C IMMEDIATELY AFTER INVERSE TRANSFORM OR EXTENDED WITH ZEROS C IMMEDIATELY BEFORE FORWARD TRANSFORM. IF THIS FEATURE WAS NOT C AVAILABLE THE DFT WOULD HAVE TO SPEND TIME FINDING OR USING C VALUES THAT DO NOT CONTRIBUTE TO THE RESULT. C C A SPECIAL FEATURE OF DFTFTR IS THAT THE FREQUENCY DOMAIN DOES C NOT INCLUDE THE NEGATIVE FREQUENCIES AND ZERO & NYQUIST C FREQUIENCIES ARE PASSED THROUGH DEDICATED ARGUMENTS. BECAUSE C THE TIME DOMAIN IS REAL, THE NEGATIVE FREQUENCIES DO NOT ADD C ANY NEW INFORMATION. THEREFORE, DROPPING THE NEGATIVE C FREQUENCIES DOES NOT ALTER THE FREQUENCY DOMAIN WHILE ALLOWING C IT TO TAKE THE SAME ARRAY SPACE AS THE TIME DOMAIN (RATHER THAN C TWICE). C C C INPUT ARGUMENTS: C KSCALE - INTEGER*4. THE TYPE OF SCALE FACTOR TO BE USED: C KSCALE DESCRIPTION C 0 FFT SCALE FACTOR C 1 DFT SCALE FACTOR C THE FFT SCALE FACTOR TRIES TO MAKE THE TIME AND C FREQUENCY DOMAINS MATHEMATICALLY EQUIVALENT BY C SCALING THE SQUARE ROOT OF THE NUMBER OF POINTS ONCE C ON FORWARD TRANSFORM AND AGAIN ON INVERSE TRANSFORM. C THE DFT SCALE FACTOR PRODUCES A FREQUENCY DOMAIN C WITH TRUE AMPLITUDES ON FORWARD TRANSFORM AND C THEREFORE DOES NO SCALING ON INVERSE TRANSFORM. THE C DFT SCALE FACTOR IS GENERALLY MORE PHYSICALLY C MEANINGFUL AND IS MORE LIKELY TO REMAIN MEANINGFUL C WHEN TRUNCATIONS AND EXTENSIONS ARE BEING USED. C JB - INTEGER*4. THE BEGINNING LOCATION IN THE TIME C DOMAIN OF THE INTERVAL (JB,JE) FROM WHICH C CALCULATIONS WILL BE MADE. ANY VALUES OUTSIDE OF C THIS INTERVAL WILL HAVE AN EFFECTIVE VALUE OF ZERO. C THE VALUE OF JB MUST BE WITHIN THE RANGE 1 TO NRX C WHERE NRX IS THE NUMBER OF POINTS IN THE TIME C DOMAIN. UPON INVERSE TRANSFORM ALL POINTS BEFORE JB C WILL BE IMPLICITLY TRUNCATED MAKING THE TRUNCATED C ARRAY LOCATION, RX(1) HAVE THE VALUE OF THE C UNTRUNCATION ARRAY LOCATION, RX(JB). C JE - INTEGER*4. THE ENDING LOCATION IN THE TIME DOMAIN C OF THE INTERVAL (JB,JE) FROM WHICH CALCULATIONS WILL C BE MADE. ANY VALUES OUTSIDE OF THIS INTERVAL WILL C HAVE AN EFFECTIVE VALUE OF ZERO. THE VALUE OF JE C MUST BE GREATER THAN OR EQUAL TO JB AND WITHIN THE C RANGE 1 TO NRX WHERE NRX IS THE NUMBER OF POINTS IN C THE TIME DOMAIN. UPON INVERSE TRANSFORM ALL POINTS C AFTER JE WILL BE IMPLICITLY TRUNCATED. THE C TRUNCATED ARRAY LOCATION, RX(JE-JB+1) HAVE THE VALUE C OF THE UNTRUNCATION ARRAY LOCATION, RX(JE). C NRX - INTEGER*4. THE NUMBER OF TIME-DOMAIN PTS. NRX IS C ALWAYS USED AS THE EFFECTIVE LENGTH OF THE C TIME-DOMAIN REGARDLESS OF THE THE VALUES IN JB AND C JE. BUT, UPON FORWARD TRANSFORM, NRX IS THE ACTUAL C NUMBER OF POINTS IN AND THE DIMENSION OF ARRAY RX; C UPON INVERSE TRANSFORM THE DIMENSION IS REDUCED TO C JE-JB+1 BECAUSE OF IMPLICIT TRUNCATION. C NCY - INTEGER*4. THE NUMBER OF FREQUENCY-DOMAIN PTS. NCY C IS ALWAYS USED AS THE EFFECTIVE LENGTH OF THE C FREQUENCY-DOMAIN. BUT, BECAUSE NEGATIVE FREQUENCIES C ARE NOT USED, NCY IS ROUGHLY TWICE THE ACTUAL NUMBER C OF POINTS IN ARRAY CYH. THE ACTUAL DIMENSION OF CYH C IS (NCY-1)/2. C SIGNI - REAL*8. DIRECTION OF TRANSFORM: C -1.D0 = FORWARD TRANSFORM (TIME TO FQ) (E^-IW). C +1.D0 = INVERSE TRANSFORM (FQ TO TIME) (E^+IW). C C INPUT/OUTPUT ARGUMENT: C RX - REAL*8. ARRAY CONTAINING REAL TIME-DOMAIN DATA. IF C SIGNI IS ZERO OR NEG (FORWARD TRANSFORM), RX WILL BE C AN INPUT ARGUMENT. IF SIGNI IS GREATER THAN ZERO C (INVERSE TRANSFORM), RX WILL BE AN OUTPUT ARGUMENT. C THE DIMENSION OF RX IS NRX DURING FORWARD TRANSFORM; C IT IS REDUCED TO JE-JB+1 DURING INVERSE. C CYH - COMPLEX*16. ARRAY OF DIMENSION (NCY-1)/2 CONTAINING C A HALF SPECTRUM OF THE COMPLEX FREQUENCY-DOMAIN C DATA. IF SIGNI IS ZERO OR NEG (FORWARD TRANSFORM), C CYH WILL BE AN OUTPUT ARGUMENT. IF SIGNI IS GREATER C THAN ZERO (INVERSE TRANSFORM), CYH WILL BE AN INPUT C ARGUMENT. THE HALF SPECTRUM IS STRUCTURED AS C FOLLOWS: CYH(1)=WAVE#1, CYH(2)=WAVE#2, ... C CYH((NCY-1)/2)=HIGHEST POSITIVE FREQUENCY THAT IS C NOT NYQUIST. ZERO FREQUENCY (WAVE#0) AND NYQUIST C (WAVE#(NCY/2)) ARE BOTH REAL AND SEPARATED TO THE C ARGUMENTS RYA (ZERO) AND RYN (NYQUIST). NOTE THAT C IF NCY IS ODD, NYQUIST DOES NOT EXIST. (A FULL C SPECTRUM WOULD CONTAIN NCY EVENLY-SPACED FREQUENCIES C IN AN ASCENDING SEQUENCE FROM ZERO TO NYQUIST AND ON C UP TO 1 POINT BEFORE TWICE NYQUIST.) C RYA - REAL*8. THE ZERO FREQUENCY (WAVE#0 OR THE AVERAGE C VALUE). THE VALUE IN RYA NORMALLY APPEARS AS A C COMPLEX NUMBER IN ARRAY LOCATION 1 OF A FULL C SPECTRUM. HOWEVER, HERE IT IS A SEPARATE ARGUMENT C TO INSURE THAT THE HALF SPECTRUM HAVE A DIMENSION NO C LARGER THAN HALF OF THE FULL SPECTRUM. C RYN - REAL*8. THE NYQUIST FREQUENCY (WAVE#(NCY/2) OR THE C HIGHEST POSSIBLE FREQUENCY). THE VALUE IN RYN C NORMALLY APPEARS AS A COMPLEX NUMBER IN ARRAY C LOCATION NCY/2+1 OF A FULL SPECTRUM. HOWEVER, HERE C IT IS A SEPARATE ARGUMENT TO INSURE THAT THE HALF C SPECTRUM HAVE A DIMENSION NO LARGER THAN HALF OF THE C FULL SPECTRUM. NOTE THAT NYQUIST ONLY EXISTS IF NCY C IS EVEN. IF NCY IS ODD OR NCY IS GREATER THAN NRX, C THE VALUE IN RYN WILL BE ZERO. C C SUBROUTINE DFTFTR WRITTEN BY ROB BRACKEN, USGS, 19970613. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 19970613. C C subroutine dftftr(kscale, & jb,je, nrx,ncy,signi,rx,cyh,rya,ryn) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 kscale integer*4 jb,je,nrx,ncy real*8 signi C C INPUT/OUTPUT ARGUMENTS real*8 rx(*) complex*16 cyh(*) real*8 rya,ryn C C INTERNAL VARIABLES AND CONSTANTS complex*16 ci,cw0,cwk,cwj,csum real*8 rsum real*8 pi,dscale real*8 powmin parameter(pi = 3.1415 92653 58979 32384 62643) integer*4 ksc2 C C C SET UP INTERNAL SCALE TYPE C C FFT TYPE (KSC2=0) OR DFT TYPE (KSC2=1) ksc2=0 if(kscale.ge.1) ksc2=1 C C CALCULATE DISCRETE FOURIER TRANSFORM WITH REAL TIME DOMAIN, C COMPLEX HALF SPECTRUM, HIGH FREQUENCY TRUNCATION/EXTENSION, AND C IMPLICIT SPECIFIABLE TIME-DOMAIN TRUNCATION. C C NOTE: IF NCY IS ODD, NYQUIST DOES NOT EXIST; IF NCY IS EVEN, C NYQUIST EXISTS, IS REAL, AND IS DOUBLE AMPLITUDE (AS IS THE C ZERO FREQUENCY) C C CALCULATE THE DISCRETE KERNEL OF INCREMENTAL NORMALIZED C ANGULAR FREQ, CW0 FOR FREQUENCY DOMAIN TRUNCATION/EXTENSION ci=dcmplx(0.d0,1.d0) cw0=cdexp(-ci*2.d0*pi/dflotj(nrx)) C C ADJ BEG & END TIM-DOM LOCS (JB, JE) INTO INTERVAL (1,NRX) jb2=jb je2=je if(jb2.lt.1) jb2=1 if(je2.lt.1) je2=1 if(jb2.gt.nrx) jb2=nrx if(je2.gt.nrx) je2=nrx if(je2.lt.jb2) then jhold=je2 je2=jb2 jb2=jhold endif C C DETERMINE WHETHER TO PERFORM FORWARD OR INVERSE TRANSFORM if(signi.gt.0.d0) goto 101 C C C FORWARD TRANSFORM (TIME DOMAIN TO FREQUENCY DOMAIN) C C INITIALIZE MINIMUM POWER CHECKER c powmin=1.1d+38 C C SET UP INTRINSIC DFT TYPE OF SCALE FACTOR c dscale=1.d0/dflotj(je2-jb2+1) dscale=1.d0/dflotj(nrx) C C CHANGE SCALE FACTOR IF FFT TYPE IS DESIRED if(ksc2.eq.0) dscale=dsqrt(dscale) C C FIND THE AVERAGE VALUE (ZERO FREQUENCY OR WAVENUMBER ZERO) rsum=0.d0 do j=jb2,je2 rsum=rsum+rx(j) enddo rya=rsum*dscale C C INITIALIZE NYQUIST VARIABLE ryn=0.d0 C C FIND THE NUMBER OF FREQS TO CALCULATE (INTERVAL (ZERO,NYQ)) nfq=jmin0(nrx,ncy)/2+1 C C FIND FREQUENCIES: WAVENUMBER 1 UP THROUGH NYQUIST cwk=cw0 do k=2,nfq csum=dcmplx(0.d0,0.d0) cwj=cwk**dflotj(jb2-1) do j=jb2,je2 csum=csum+rx(j)*cwj cwj=cwj*cwk enddo cwk=cwk*cw0 c if(cdabs(csum).lt.powmin) powmin=cdabs(csum) C C SCALE AND PUT IN LEFT HALF (AVG IS NOT IN CYH(1)) if(k.lt.nfq) cyh(k-1)=csum*dscale enddo k=k-1 C C MAKE ADJUSTMENTS AT AND ABOVE THE NYQUIST FREQUENCY nyqck=(nfq-1)*2 if(ncy.gt.nrx) then C C IF CSUM IS THE OLD NYQUIST, SPLIT IT if(nyqck.eq.nrx) csum=csum/2.d0 cyh(k-1)=csum*dscale C C FILL EXTRA FREQUENCIES WITH ZEROS powmin=0.d0 C (OTHER OPTIONS ARE PICK A POWER LEVEL BELOW POWMIN C OR CREATE WHITE NOISE AT SOME SPECIFIED POWER) c powmin=powmin*dscale/4096.d0 c powmin=powmin*dscale/2.d0 do i=k,(ncy-1)/2 cyh(i)=dcmplx(powmin,0.d0) enddo else C C IF CSUM IS NEW NYQUIST, COMBINE WITH ITS REFLECTION if(nyqck.eq.ncy) then if(ncy.lt.nrx) csum=csum+dconjg(csum) ryn=dreal(csum)*dscale else cyh(k-1)=csum*dscale endif endif C goto 990 C C C INVERSE TRANSFORM (FREQUENCY DOMAIN TO TIME DOMAIN) C C INVERT THE KERNEL 101 cw0=dconjg(cw0) C C SET UP INTRINSIC DFT TYPE OF INVERSE SCALE FACTOR ( = 1 ) dscale=1.d0 C C CHANGE INVERSE SCALE FACTOR IF FFT TYPE IS DESIRED c if(ksc2.eq.0) dscale=1.d0/dsqrt(dflotj(je2-jb2+1)) if(ksc2.eq.0) dscale=1.d0/dsqrt(dflotj(nrx)) C C FIND LOC OF FREQ BEFORE NYQ & FIND WHETHER NYQUIST EXISTS kbn=(jmin0(nrx,ncy)-1)/2+1 nyq=0 if(kbn*2.eq.jmin0(nrx,ncy)) nyq=1 C C FIND TIME-DOM POINTS; UPPER HALF OF FREQS CAN BE IGNORED cwj=cw0**dflotj(jb2-1) do j=jb2,je2 rsum=0.d0 C C ADD IN FREQS WITH WAVENUMBERS 1 THROUGH ALMOST NYQUIST cwk=cwj do k=1,kbn-1 rsum=rsum+dreal(cyh(k)*cwk) cwk=cwk*cwj enddo rsum=2.d0*rsum cwj=cwj*cw0 C C ADD IN NYQUIST IF IT EXISTS (MODIFY IT IF NECESSARY) if(nyq.eq.1) then if(ncy.gt.nrx) then rsum=rsum+dreal((cyh(k)+dconjg(cyh(k)))*cwk) else rsum=rsum+dreal(ryn*cwk) endif endif C C ADD IN AVERAGE rsum=rsum+rya C C SCALE THE SUM rx(j-jb2+1)=rsum*dscale enddo C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE E R R F O R C________________________________________________________________ C C SUBROUTINE ERRFOR CAUSES A FORTRAN PROGRAM TO STOP WITH ERROR C STATUS AND SENDS AN ERROR NUMBER AND MESSAGE TO THE TERMINAL. C C IERR - INTEGER*4. THE ERROR NUMBER TO BE DISPLAYED BEFORE C THE PROGRAM EXECUTION IS TERMINATED. C MSG - CHARACTER*(*). THE MESSAGE TEXT TO BE DISPLAYED. C C IN CASE YOU ARE WONDERING, THIS SUBROUTINE IS USEFUL FOR C DELIBERATE PROGRAM TERMINATION WHICH PIN POINTS A PROGRAMMER C DETERMINED ERROR (RATHER THAN A SYSTEM DETERMINED ERROR). C IT CAUSES A FATAL SYSTEM ERROR WHICH STOPS BATCH JOBS SO THAT C THEY CAN'T GO ON AND MUNCH OR DESTROY GOOD DATA. C C NOTE: IN UNIX SHELL (SH), A SCRIPT WILL CONTINUE TO EXECUTE C AFTER THE PROGRAM ERRFOR GENERATES A FATAL ERROR. TO CAUSE C THE SCRIPT TO TERMINATE, INSERT THE COMMAND set -e IN THE C THE SCRIPT FILE. C C SUBROUTINE ERRFOR WRITTEN BY ROB BRACKEN, USGS, 20SEP87. C VAX 11/780 VERSION 1.0, 20SEP87. C HP9000 SERIES 300/400 VERSION 1.1, 30MAR92. C HP9000 SERIES 300/400 VERSION 1.2, 19970324. C C subroutine errfor(ierr,msg) C C DECLARATIONS C C ARGUMENTS integer*4 ierr character*(*) msg C C OTHER VARIABLES integer*4 ttyo parameter(ttyo=6) character*80 blank c equivalence(iarmed,blast) data blank /' '/ C C PRODUCE A MESSAGE AT THE TERMINAL C write(ttyo,801) ierr,msg//blank 801 format(/x,'TO TERMINATE SCRIPT, USE set -e', & /x,'%PGM-F-ERRORFOR, ERROR #',i12,/2x,a79) C C CLOSE UNIT TTYO TO FLUSH BUFFER C (NECESSARY ON UNIX MACHINES WITH HUGE OUTPUT BUFFERS) C close(unit=ttyo) C C PRODUCE A FATAL ERROR C c iarmed=-1 c iarmed=blast c pause 'do NOT restart program' call exit(2) C C EXIT PROCEDURE (NOT USED) C 990 return end C C________________________________________________________________ C C SUBROUTINE E X T E N D C________________________________________________________________ C C SUBROUTINE EXTEND EXTENDS A FUNCTION PROFILE IN ONE OR BOTH C DIRECTIONS USING VARIOUS OPTIONS FOR TYPE, LENGTH, AND C DIRECTION OF EXTENSION. IT WILL ALSO PERFORM THE INVERSE BY C TRUNCATION. C C INPUT ARGUMENTS: C KTYPE - INTEGER*4. THE TYPE OF EXTENSION TO PRODUCE. C FOLLOWING ARE POSSIBLE VALUES. A STAR (*) INDICATES C THAT THE OPTION HAS BEEN INSTALLED: C * 0 = ZEROS C * 1 = AVERAGE VALUE C * 2 = FIRST/LAST VALUE C * 3 = REPEAT C * 4 = REFLECT C * 5 = REFLECT AND INVERT (TPARMS NEEDED) C TPARMS(1) IS THE VALUE ABOUT WHICH TO C INVERT THE LEFT EXTENSION. TPARMS(2) IS C THE VALUE ABOUT WHICH TO INVERT THE RIGHT C EXTENSION. C 6 = SPLINE TO ZERO (TPARMS NEEDED) C * 7 = BURG PREDICTION WITH LINEAR FUNCTION C REMOVED BEFORE PREDICTION AND ADDED BACK C AFTER PREDICTION. (TPARMS NEEDED) C TPARMS(1) IS THE "ROOT" LENGTH OR THE C METHOD TO USE FOR OBTAINING THE LENGTH OF C THE PROFILE FROM WHICH FILTER COEFS ARE C DERIVED: C -3) INT(TPARMS(2))=LENGTH OF LEFT ROOT C INT(TPARMS(3))=LENGTH OF RIGHT ROOT C -2) TPARMS(2)*(L AND R EXTENSION AMOUNTS) C -1) LEFT AND RIGHT EXTENSION AMOUNTS C 0) INPUT ARRAY LENGTH (NIN) C 1) HALF THE TOT EXTND AMT (NOUT-NIN)/2 C >1) INT(TPARMS(1)) C * 8 = BURG PREDICTION AS IN KTYPE=7 ONLY WITHOUT C THE LINEAR FUNCTION. (TPARMS SAME AS 7). C C TPARMS - REAL*8. ARRAY CONTAINING ANY PARAMETERS ASSOCIATED C WITH ONE OF THE KTYPE NUMBERS. THIS ARGUMENT WILL C NOT BE USED UNLESS DESCRIBED UNDER KTYPE. C JUST - INTEGER*4. JUSTIFICATION OF THE ORIGINAL PROFILE C WITHIN THE EXTENDED PROFILE. FOLLOWING ARE POSSIBLE C VALUES: C <-1 = EXACT LOCATION OF LAST POINT FROM RIGHT C -1 = RIGHT JUSTIFY (LAST POINT AT END), C 0 = CENTER JUSTIFY, C 1 = LEFT JUSTIFY (FIRST POINT AT LOCATION 1) C >1 = EXACT LOCATION OF FIRST POINT FROM LEFT C NIN - INTEGER*4. THE NUMBER OF POINTS INPUT. IF NIN IS C LESS THAN NOUT, PROFILE WILL BE EXTENDED. C NOUT - INTEGER*4. THE NUMBER OF POINTS TO OUTPUT. IF NOUT C IS LESS THAN NIN, PROFILE WILL BE TRUNCATED. C C INPUT/OUTPUT ARGUMENT: C YFNC - REAL*8. ARRAY OF DIMENSION NIN OR NOUT (WHICH EVER C IS GREATER) CONTAINING THE DATA TO BE EXTENDED (OR C TRUNCATED). C C SUBROUTINE EXTEND WRITTEN BY ROB BRACKEN, USGS, 12NOV96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 12NOV96. C VERSION 1.1, 20010608. (FIXED BUG IN KTYPE=4) C C subroutine extend(ktype,tparms,just,nin,nout,yfnc) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 ktype real*8 tparms(*) integer*4 just,nin,nout C C INPUT/OUTPUT ARGUMENT real*8 yfnc(*) C C INTERNAL VARIABLES integer*4 mxburg parameter(mxburg=200000) real*8 avg,apf(mxburg),acoef(2) C C DO NOTHING IF EQUAL NUMBERS OF INPUT AND OUTPUT POINTS C if(nout.eq.nin) goto 990 C C DETERMINE NUMBER OF EXTENDED ELEMENTS ON THE LEFT SIDE C C FIND NPTS AND NPTSX npts=jmin0(nin,nout) nptsx=jmax0(nin,nout) C C MAKE DETERMINATION if(just.ge.1) left=just-1 if(just.eq.0) left=(nptsx-npts)/2 if(just.le.-1) left=nptsx-npts+just+1 if(left.lt.0.or.left.gt.nptsx-npts) call errfor(left, & '(extend) npts and nptsx inconsistant w/r just') C C PERFORM TRUNCATION IF FEWER OUTPUT THAN INPUT POINTS C if(nout.lt.nin) then C C PERFORM TRUNCATION if(left.gt.0) then do i=1,npts yfnc(i)=yfnc(i+left) enddo endif goto 990 endif C C PERFORM EXTENSION C C MOVE PROFILE TO RIGHT TO MAKE ROOM FOR THE EXTENSION C if(left.gt.0) then do i=npts,1,-1 yfnc(left+i)=yfnc(i) enddo endif C C FIND LIMITS OF EXTENSION SECTIONS C C LEFT SIDE il=1 jl=left C C RIGHT SIDE ir=left+npts+1 jr=nptsx C C DETERMINE THE TYPE OF EXTENSION C goto (100,101,102,103,104,105,106,107,108),ktype+1 call errfor(ktype, & '(extend) Non-existant extension-type number') C C ZERO EXTENSION C 100 do i=il,jl yfnc(i)=0.d0 enddo do i=ir,jr yfnc(i)=0.d0 enddo goto 990 C C AVERAGE-VALUE EXTENSION C 101 avg=0.d0 do i=jl+1,ir-1 avg=avg+yfnc(i) enddo avg=avg/dflotj(npts) do i=il,jl yfnc(i)=avg enddo do i=ir,jr yfnc(i)=avg enddo goto 990 C C FIRST AND LAST VALUE EXTENSION C 102 do i=il,jl yfnc(i)=yfnc(jl+1) enddo do i=ir,jr yfnc(i)=yfnc(ir-1) enddo goto 990 C C REPEAT EXTENSION C 103 do i=jl,il,-1 yfnc(i)=yfnc(i+npts) enddo do i=ir,jr yfnc(i)=yfnc(i-npts) enddo goto 990 C C REFLECT EXTENSION C 104 j=jl+1 kdir=-1 do i=jl,il,-1 if(j.ge.ir-1.or.j.le.jl+1) kdir=-kdir j=j+kdir yfnc(i)=yfnc(j) enddo C j=ir-1 kdir=1 do i=ir,jr if(j.ge.ir-1.or.j.le.jl+1) kdir=-kdir j=j+kdir if(j.le.0) then yfnc(i)=yfnc(1) else yfnc(i)=yfnc(j) endif enddo goto 990 C C REFLECT-AND-INVERT EXTENSION C 105 j=jl+1 kdir=-1 do i=jl,il,-1 if(j.ge.ir-1.or.j.le.jl+1) kdir=-kdir j=j+kdir yfnc(i)=tparms(1)+(tparms(1)-yfnc(j)) enddo C j=ir-1 kdir=1 do i=ir,jr if(j.ge.ir-1.or.j.le.jl+1) kdir=-kdir j=j+kdir yfnc(i)=tparms(2)+(tparms(2)-yfnc(j)) enddo goto 990 C C SPLINE-TO-ZERO EXTENSION C 106 call errfor(ktype, & '(extend) Spline-to-zero extension not installed') C C BURG-PREDICTION EXTENSION C NOTE: FILTER LENGTH IS NOW .707*(FILTER "ROOT" LENGTH) C C FIND LEFT AND RIGHT FILTER "ROOT" LENGTHS 108 continue 107 jp1=jidnnt(tparms(1)) if(jp1.lt.-3) jp1=0 C if(jp1.lt.0) then lenl=jl-il+1 lenr=jr-ir+1 if(jp1.eq.-2) then lenl=jidnnt(dflotj(lenl)*tparms(2)) lenr=jidnnt(dflotj(lenr)*tparms(2)) else if(jp1.eq.-3) then lenl=jidnnt(tparms(2)) lenr=jidnnt(tparms(3)) endif else lenl=npts if(jp1.eq.1) lenl=(nptsx-npts)/2 if(jp1.gt.1) lenl=jp1 lenr=lenl endif C if(lenl.lt.3) lenl=3 if(lenr.lt.3) lenr=3 if(lenl.gt.npts) lenl=npts if(lenr.gt.npts) lenr=npts c lenlf=lenl c lenrf=lenr lenlf=jnint(floatj(lenl-2)/sqrt(2.))+2 lenrf=jnint(floatj(lenr-2)/sqrt(2.))+2 if(npts.le.1) goto 102 C C REMOVE LSQ LINE if(ktype.eq.7) call trendr(1,2,acoef,0,npts,yfnc(jl+1)) C C EXTEND TO THE LEFT if(lenlf.gt.mxburg) call errfor(lenlf, & '(extend) lenlf exceeded upper dimension of array apf()') call burgdp(lenl,yfnc(jl+1),lenlf,apf) do j=jl,il,-1 c yfnc(j)=0.d0 jm1=j-1 avg=0.d0 do i=2,lenlf c yfnc(j)=yfnc(j)-yfnc(j+i-1)*apf(i) avg=avg-yfnc(jm1+i)*apf(i) enddo yfnc(j)=avg enddo C C EXTEND TO THE RIGHT if(lenrf.gt.mxburg) call errfor(lenrf, & '(extend) lenrf exceeded upper dimension of array apf()') if(lenr.ne.lenl.or.lenr.ne.npts) & call burgdp(lenr,yfnc(ir-lenr),lenrf,apf) do j=ir,jr c yfnc(j)=0.d0 jp1=j+1 avg=0.d0 do i=2,lenrf c yfnc(j)=yfnc(j)-yfnc(j-i+1)*apf(i) avg=avg-yfnc(jp1-i)*apf(i) enddo yfnc(j)=avg enddo C C ADD-BACK LSQ LINE if(ktype.eq.7) then acoef(1)=acoef(1)-acoef(2)*dflotj(left) call trendr(1,2,acoef,1,nptsx,yfnc) endif goto 990 C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE F A C T R C________________________________________________________________ C C SUBROUTINE FACTR FINDS ALL OF THE PRIME FACTORS OF A GIVEN C INTEGER VALUE. C C INPUT ARGUMENTS: C FVAL - REAL*8. THE VALUE TO FACTOR. IF FVAL IS NOT AN C INTEGER, THE NEAREST INTEGER WILL BE ASSUMED FOR C FACTORING. FVAL MAY BE POSITIVE OR NEGATIVE. IF C FVAL IS ZERO, THE FACTOR ARRAY WILL CONTAIN 0 1 0. C THE MAXIMUM FACTORABLE VALUE IS APPROXIMATELY C 7.21E+16 (EXACTLY 2**56-1) ON THE VAX 11/780 AND C 9.01E+15 (EXACTLY 2**53-1) ON THE IEEE MACHINES. C ANY VALUE LARGER THAN THESE WILL GENERATE AN ERROR. C NOTE THAT THESE VALUES ARE GREATER THAN THE MAXIMUM C 32 BIT INTEGER (I*4) AND CONTAIN MORE PRECISION THAN C THE SINGLE PRECISION VARIABLES (R*4). NOTE ALSO C THAT FACTORING VERY LARGE PRIME NUMBERS MAY BE TIME C CONSUMING. C C OUTPUT ARGUMENTS: C FACTRS - REAL*8. ARRAY CONTAINING THE INTEGRAL FACTORS. C THE FIRST ARRAY ELEMENT ALWAYS CONTAINS THE PRODUCT C OF THE FACTORS. THE SECOND ARRAY ELEMENT ALWAYS C CONTAINS THE SIGN FACTOR, 1 OR -1. SUBSEQUENT ARRAY C ELEMENTS CONTAIN THE LIST OF FACTORS IN ORDER OF C INCREASING VALUE. FACTORS OCCURING MORE THAN ONCE C WILL BE REPEATED. IN ORDER TO ACCOMODATE ALL C POSSIBLE FACTORING COMBINATIONS, ARRAY FACTRS SHOULD C BE DIMENSIONED 57 OR MORE ELEMENTS. C NFACT - INTEGER*4. THE NUMBER OF ARRAY ELEMENTS USED. C C SUBROUTINE FACTR WRITTEN BY ROB BRACKEN, USGS, 25FEB92. C VAX/VMS 11/780 VERSION 2.0, 14MAR92. C HP9000 SERIES 300/400, VERSION 2.0, 14MAR92. C C subroutine factr(fval,factrs,nfact) C C DECLARATIONS C C ARGUMENTS real*8 fval,factrs(*) C C FACTORING VARIABLES real*8 rval,sval,tval,uval real*8 sfact integer*4 mxchk C C MAXIMUM FACTORABLE NUMBER = 2**56-1 ON THE VAX 11/780 C MAXIMUM FACTORABLE NUMBER = 2**53-1 ON THE HP9000 real*8 xval c parameter(xval='ffffffffffff5c7f'x) parameter(xval='433fffffffffffff'X) C C MAXIMUM INTEGER VALUE parameter(idval='7fffffff'X) C C GET VALUE TO FACTOR C rval=fval C C COMPARE WITH MAXIMUM VALUE C if(dabs(rval).gt.xval) & call errfor(0,'Exceeded maximum factorable value') C C FIND THE NEAREST INTEGER VALUE C sval=rval if(dabs(sval).le.(xval-1.d0)/2.d0) sval=dnint(sval) if(dabs(sval-rval).gt.0.5d0) & call errfor(0,'Round off error') rval=sval C C CHECK FOR THE ZERO CASE C if(sval.eq.0.d0) then factrs(1)=0.d0 factrs(2)=1.d0 factrs(3)=0.d0 nfact=3 goto 990 endif C C FIND THE FIRST FACTOR (POSITVE OR NEGATIVE ONE) C nfact=2 if(sval.ge.0.d0) then factrs(nfact)=1.d0 else factrs(nfact)=-1.d0 sval=-sval endif C C FIND THE SECOND FACTOR (TWO) C 202 tval=sval/2.d0 if(tval.eq.dint(tval)) then nfact=nfact+1 factrs(nfact)=2.d0 sval=tval goto 202 else if(tval.lt.2.d0) then if(sval.gt.1.d0) goto 103 goto 104 endif C C FIND THE MAXIMUM VALUE TO CHECK C mxchk=jidint(dsqrt(sval)) C C FIND THE REST OF THE FACTORS C do 705 ifact=3,mxchk,2 C C ESTABLISH THE FACTOR AND THE NEW VALUE sfact=dflotj(ifact) 206 tval=sval/sfact C C TRUNCATE THE NEW VALUE TO AN INTEGER uval=tval if(uval.gt.idval) then jval=uval*1.d-8 uval=uval-dflotj(jval)*1.d+8 endif jval=uval C C CHECK DIVISIBILITY if(uval.eq.jval) then nfact=nfact+1 factrs(nfact)=dflotj(ifact) sval=tval goto 206 else if(tval.lt.sfact) then if(sval.gt.1.d0) goto 103 goto 104 endif 705 continue C C INCLUDE THE FINAL FACTOR C 103 nfact=nfact+1 factrs(nfact)=sval C C CHECK THE FACTORS AGAINST THE ORIGINAL NUMBER C 104 tval=1.d0 do 707 i=2,nfact tval=tval*dnint(factrs(i)) 707 continue C C INCLUDE THE ORIGINAL VALUE C factrs(1)=tval C C IF THE CALCULATED VALUE DOES NOT MATCH THE ORIGINAL VALUE, C GENERATE AN ERROR C if(tval.ne.rval) call errfor(0,'Failed internal check') C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE F E J E R C________________________________________________________________ C C SUBROUTINE FEJER REMOVES OR APPLIES A TRIANGLE (FEJER OR C BARTLETT) WINDOW IN AN EVENLY SPACED ARRAY OF REAL DATA. C C THE FEJER WINDOW IS MERELY A TRIANGULAR SHAPE WITH A STRAIGHT C LINE FROM 0 VALUE AT THE BEGINNING TO 1 IN THE MIDDLE AND THE C SAME REFLECTED ACROSS THE MIDDLE INTO THE SECOND HALF. C CHANGING THE INTENSITY VARIABLE, ALPHA HAS NO EFFECT IN THIS C SUBROUTINE. FOR VARIATIONS ON THE FEJER, SEE THE GENERALIZED C RIESZ WINDOW. C C THE BASIC EQUATION IS C C W(J)=1-(2J/M) C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE INTENSITY OF THE WINDOW. C ALPHA HAS NO EFFECT IN THIS SUBROUTINE. IF USING C SUBROUTINE RIESZ, ALPHA=1 PRODUCES A FEJER WINDOW C BUT NOT AS QUICKLY AS THIS SUBROUTINE. C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE FEJER WRITTEN BY ROB BRACKEN, USGS, 19970609. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 19970609. C C subroutine fejer(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp real*8 winval,xr,xctr,scale,xj C C INITIALIZATIONS C C FIND POINT WITH HIGHEST INDEX ON INTERVAL 1 TO N C DEPENDING ON WHETHER SYMMETRIC OR DFT EVEN nwp=nio if(kdft.ge.1) nwp=nio+1 C C SLIDE INTERVAL LEFT 1 POINT TO BE 0 TO N-1 nwp=nwp-1 xr=dflotj(nwp) C C FIND VALUE OF X CORRESPONDING TO THE CENTER OF THE PROFILE xctr=0.5d0*xr C C FIND SCALE FACTOR WHICH EQUALS 1 WHEN MULTIPLIED BY XCTR. scale=1.d0/xctr C C FIND AND APPLY (OR REMOVE) THE WINDOW VALUES C C APPLY (OR REMOVE) ALL EXCEPT BEG PT (& END PT IF SYMMETRIC) do j=1,jidint(xctr) xj=dflotj(j) C C FEJER WINDOW winval=scale*xj C C DETERMINE WHETHER TO APPLY OR REMOVE THE WINDOW if(nrem.ge.1) winval=1.d0/winval C C APPLY (OR REMOVE) THE LEFT HALF OF THE WINDOW yio(j+1)=yio(j+1)*winval C C APPLY (OR REMOVE) THE RIGHT HALF OF THE WINDOW yio(nwp-j+1)=yio(nwp-j+1)*winval enddo C C APPLY (OR REMOVE) FIRST POINT (AND LAST POINT IF SYMMETRIC) if(nrem.le.0) then yio(1)=0.d0 if(kdft.le.0) yio(nio)=0.d0 else yio(1)=yio(2) if(kdft.le.0) yio(nio)=yio(nio-1) endif C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE H A M M 3 C________________________________________________________________ C C SUBROUTINE HAMM3 REMOVES OR APPIES A HAMMING WINDOW IN AN C EVENLY SPACED ARRAY OF REAL DATA. C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE INTENSITY OF THE WINDOW. C IN THE CASE OF THE HAMMING WINDOW, ALPHA HAS NO USE; C IT MAY BE SET TO ANY VALUE. C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE HAMM3 WRITTEN BY ROB BRACKEN, USGS, 22AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 3.0, 22AUG96. C C subroutine hamm3(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp,jdmax real*8 half,pi,xd,cosp,winval parameter(pi = 3.1415 92653 58979 32384 62643) C C INITIALIZATIONS C C EFFECTIVE NUMBER OF WINDOW POINTS nwp=nio if(kdft.ge.1) nwp=nwp+1 C THE MAXIMUM DISTANCE FROM THE LEFT MOST PROFILE POINT jdmax=nwp/2-1 C THE EXACT LENGTH OF HALF THE PROFILE half=dflotj(nwp-1)/2.d0 C C WINDOW ONE POINT AND IT'S REFLECTION ON EACH LOOP C do 701 jd=0,jdmax C C THE DISTANCE TRANSLATED TO XD=(n/(N/2)) xd=1.d0-dflotj(jd)/half C C THE WINDOW FUNCTION cosp=0.54d0+0.46d0*dcos(xd*pi) C if(nrem.le.0) then C C THE WINDOW VALUE winval=cosp else C C THE WINDOW INVERSE VALUE winval=1.d0/cosp endif C C WINDOW POINT NUMBERS JD+1 AND NWP-JD yio(jd+1)=yio(jd+1)*winval if(jd.gt.0.or.kdft.le.0) yio(nwp-jd)=yio(nwp-jd)*winval 701 continue C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE H A N N 3 C________________________________________________________________ C C SUBROUTINE HANN3 REMOVES OR APPIES A HANNING WINDOW IN AN C EVENLY SPACED ARRAY OF REAL DATA. C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE INTENSITY OF THE WINDOW. C TYPICAL VALUES ARE 1.0 TO 4.0. THE HIGHER NUMBERS C DECREASE THE NEAREST SIDE-LOBE BUT INCREASE THE BAND C WIDTH (INCREASE LOSS). IF ALPHA IS 1, THE HANNING C BECOMES A SIMPLE COSINE TAPER EQUIVIALENT TO A TUKEY C WITH ALPHA=1. IF ALPHA IS LESS THAN 1, A C DISCONTINUITY IN THE FIRST DERIVATIVE DEVELOPS AND C INCREASES UNTIL AS ALPHA APPROACHES ZERO, IT BECOMES C A RECTANGULAR WINDOW. NOTE THAT THE HANNING DOES C NOT MIMIC THE TUKEY FOR ALPHA LESS THAN 1. C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE HANN3 WRITTEN BY ROB BRACKEN, USGS, 22AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 3.0, 22AUG96. C C subroutine hann3(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp,jdmax real*8 half,pi,xd,cosp,winval parameter(pi = 3.1415 92653 58979 32384 62643) C C INITIALIZATIONS C C EFFECTIVE NUMBER OF WINDOW POINTS nwp=nio if(kdft.ge.1) nwp=nwp+1 C THE MAXIMUM DISTANCE FROM THE LEFT MOST PROFILE POINT jdmax=nwp/2-1 C THE EXACT LENGTH OF HALF THE PROFILE half=dflotj(nwp-1)/2.d0 C C WINDOW ONE POINT AND IT'S REFLECTION ON EACH LOOP C do 701 jd=0,jdmax C C THE DISTANCE TRANSLATED TO XD=(n/(N/2)) xd=1.d0-dflotj(jd)/half C C THE WINDOW FUNCTION cosp=(dcos(xd*pi/2.d0))**alpha C if(nrem.le.0) then C C THE HANNING WINDOW VALUE winval=cosp else C C THE HANNING WINDOW INVERSE VALUE if(jd.ne.0) winval=1.d0/cosp endif C C WINDOW POINT NUMBERS JD+1 AND NWP-JD if(jd.gt.0.or.nrem.le.0) then yio(jd+1)=yio(jd+1)*winval if(jd.gt.0.or.kdft.le.0) yio(nwp-jd)=yio(nwp-jd)*winval endif 701 continue C C EXTRAPOLATION TO BEG AND END POINTS IF INVERSE WINDOW C if(nrem.ge.1) then yio(1)=yio(2) if(kdft.le.0) yio(nwp)=yio(nwp-1) endif C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE K A I S 3 C________________________________________________________________ C C SUBROUTINE KAIS3 REMOVES OR APPIES A KAISER-BESSEL WINDOW IN AN C EVENLY SPACED ARRAY OF REAL DATA. C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE INTENSITY OF THE WINDOW. C ALPHA = 0.0 DOES NO WINDOWING; THE BEST VALUE IS C ALPHA = 3.0; VALUES OF ALPHA > 3.5 ARE NOT C RECOMMENDED. C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE KAISER WRITTEN BY MIKE WEBRING, USGS, DATE? C MODIFIED TO KAIS2 BY ROB BRACKEN, USGS, 11OCT91. C MODIFIED TO KAIS3 BY ROB BRACKEN, USGS, 20AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 2.1, 9AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 3.0,20AUG96. C C subroutine kais3(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp,jdmax real*8 half,pi,pial,pialx,xd,winval,zmbes5,zmpial parameter(pi = 3.1415 92653 58979 32384 62643) C C INITIALIZATIONS C C EFFECTIVE NUMBER OF WINDOW POINTS nwp=nio if(kdft.ge.1) nwp=nwp+1 C THE MAXIMUM DISTANCE FROM THE LEFT MOST PROFILE POINT jdmax=nwp/2-1 C THE EXACT LENGTH OF HALF THE PROFILE half=dflotj(nwp-1)/2.d0 C PI * ALPHA pial=pi*alpha zmpial=zmbes5(pial) C C WINDOW ONE POINT AND IT'S REFLECTION ON EACH LOOP C do 701 jd=0,jdmax C C THE DISTANCE TRANSLATED TO XD=(n/(N/2)) xd=1.d0-dflotj(jd)/half C C THE ARGUMENT OF THE BESSEL FUNCTION IN THE NUMERATOR pialx=pial*dsqrt(1.d0-xd*xd) C if(nrem.le.0) then C C THE KAISER-BESSEL WINDOW VALUE winval=zmbes5(pialx)/zmpial else C C THE KAISER-BESSEL WINDOW INVERSE VALUE winval=zmpial/zmbes5(pialx) endif C C WINDOW POINT NUMBERS JD+1 AND NWP-JD yio(jd+1)=yio(jd+1)*winval if(jd.gt.0.or.kdft.le.0) yio(nwp-jd)=yio(nwp-jd)*winval 701 continue C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE L S Q P L Y C________________________________________________________________ C C SUBROUTINE LSQPLY REMOVES OR ADDS BACK A POLYNOMIAL OF ORDER 0, C 1, OR 2 FROM THE INPUT FUNCTION DATA ARRAY. THE METHOD IS C GIVEN IN "RANDOM DATA ANALYSIS AND MEASUREMENT PROCEDURES", C SECOND EDITION, BY BENDAT AND PIERSOL, WILEY & SONS, PAGE 363. C C NOTE: A MORE GENERAL (HIGHER ORDER CAPABILITIES) VERSION OF C LSQPLY IS MAY BE FOUND IN DIRECTORY: C /DATA/RBRACKEN/ C TMGS.DIR/ C HI.DIR/ C PROC.DIR/ C FFT_PGMS_TABLED_G0914.DIR/ C SUMNXK_PRCSUM_BITSTUFF.DIR C LSQPLY0* C C C C INPUT ARGUMENTS: C NCOEF - INTEGER*4. NUMBER OF COEFS IN ACOEF. THE ORDER OF C THE POLYNOMIAL FUNCTION WILL BE NCOEF-1. MAXIMUM C VALUE IS 2. IF NCOEF IS MADE NEGATIVE AND NVERSE <= C 0, INSTEAD OF FINDING NEW COEFS, THE COEFS GIVEN IN C ACOEF WILL BE USED FOR REMOVING THE TREND. IF NCOEF C IS NEGATIVE AND NVERSE => 1, THERE WILL BE NO EFFECT C DIFFERENT FROM THAT OF POSITIVE NCOEF (I.E. THE C COEFS IN ACOEF WILL BE USED FOR ADDING BACK THE C TREND). C ACOEF - REAL*8. THE COEFS OF THE POLYNOMIAL IN SEQUENCE C ZEROTH, FIRST ORDER, SECOND ORDER, ETC. IF THE C TREND IS BEING REMOVED, ACOEF IS AN OUTPUT ARGUMENT. C IF THE TREND IS BEING ADDED BACK, ACOEF IS AN INPUT C ARGUMENT. C NVERSE - INTEGER*4. DETERMINES WHETHER TO REMOVE OR ADD BACK C THE TREND: C NVERSE <= 0 MEANS TO REMOVE THE TREND. C NVERSE => 1 MEANS TO ADD BACK THE TREND. C NPTS - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YFNC. C IF NPTS IS MADE NEGATIVE AND NVERSE <= 0, THE TREND C WILL NOT BE REMOVED (THE COEFS CAN THEN BE FOUND C WITHOUT CHANGING THE DATA ARRAY). IF NPTS IS C NEGATIVE AND NVERSE => 1, THE TREND WILL NOT BE C ADDED BACK (THIS HAS NO PARTICULAR USE BUT IS C PROVIDED FOR FUNCTIONAL BALANCE). C C INPUT/OUTPUT ARGUMENT: C YFNC - REAL*8. ARRAY OF DIMENSION NPTS CONTAINING THE C FUNCTION. C C SUBROUTINE LSQPLY WRITTEN BY ROB BRACKEN, USGS, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 2.0, 08SEP96. C HP-9000 SERIES 700/800 UNIX VERSION 2.1, 19970611. C ADDITION OF OPTION TO FIND COEFS BUT NOT REMOVE THE TREND C C subroutine lsqply(ncoef,acoef, nverse,npts,yfnc) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 ncoef real*8 acoef(*) integer*4 nverse,npts C C INPUT/OUTPUT ARGUMENT real*8 yfnc(*) C C INTERNAL VARIABLES real*8 tverse real*8 den,yn(3),pwr real*8 x0,y0,x1,x2,x3,x4,y1,y2 C C C CHECK NUMBER OF COEFS C if(jiabs(ncoef).gt.3) & call errfor(jiabs(ncoef),'number of coefs out of range') if(ncoef.eq.0) goto 990 C C FIND DIRECTION OF CALCULATION C if(nverse.ge.1) then C C ADD BACK THE TREND GIVEN BY ACOEF tverse=1.d0 goto 101 else C C FIND ACOEF AND REMOVE THE TREND tverse=-1.d0 if(ncoef.lt.0) goto 101 endif C C C FIND COEFFICIENTS C C SET UP SUMMATIONS do 704 j=1,ncoef yn(j)=0.d0 pwr=dflotj(j-1) do 705 i=1,jiabs(npts) yn(j)=yn(j)+yfnc(i)*dflotj(i)**pwr 705 continue 704 continue C C SOLVE SYSTEM OF EQUATIONS FOR THE COEFS if(ncoef.eq.1) then C C ONE COEF (ZEROTH ORDER = AVERAGE) x0=dflotj(jiabs(npts)) acoef(1)=yn(1)/x0 else if(ncoef.eq.2) then C C TWO COEFS (FIRST ORDER = LSQ LINE) x0=dflotj(jiabs(npts)) acoef(1)=(2.d0*(2.d0*x0+1.d0)*yn(1)-6.d0*yn(2))/ & (x0*(x0-1.d0)) acoef(2)=(12.d0*yn(2)-6.d0*(x0+1.d0)*yn(1))/ & (x0*(x0-1.d0)*(x0+1.d0)) else if(ncoef.eq.3) then C C THREE COEFS (SECOND ORDER = LSQ PARABOLA) x0=dflotj(jiabs(npts)) x1=x0*(x0+1.d0)/2.d0 x2=x0*(x0+1.d0)*(2.d0*x0+1.d0)/6.d0 x3=x0*(x0+1.d0)*x0*(x0+1.d0)/4.d0 x4=x0*(x0+1.d0)*(2.d0*x0+1.d0)*(3.d0*x0*x0+3.d0*x0-1.d0)/ & 30.d0 y0=yn(1) y1=yn(2) y2=yn(3) C den= & x0*(x2*x4-x3*x3)-x1*(x1*x4-x3*x2)+x2*(x1*x3-x2*x2) acoef(1)= & (y0*(x2*x4-x3*x3)-x1*(y1*x4-x3*y2)+x2*(y1*x3-x2*y2))/den acoef(2)= & (x0*(y1*x4-x3*y2)-y0*(x1*x4-x3*x2)+x2*(x1*y2-y1*x2))/den acoef(3)= & (x0*(x2*y2-y1*x3)-x1*(x1*y2-y1*x2)+y0*(x1*x3-x2*x2))/den endif C C C REMOVE OR ADD THE POLYNOMIAL TREND C 101 if(npts.le.0) goto 990 do 702 i=1,npts x0=dflotj(i) y0=0.d0 do 703 j=1,jiabs(ncoef) pwr=dflotj(j-1) y0=y0+(acoef(j)*(x0**pwr)) 703 continue yfnc(i)=yfnc(i)+tverse*y0 702 continue C C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE R C L O W C________________________________________________________________ C C SUBROUTINE RCLOW FILTERS A SIGNAL ARRAY USING A LOW-PASS C R-C TIME-DOMAIN FILTER WITH ANY NUMBER OF POLES. C C THE TRANSFER FUNCTION FOR EACH R-C SECTION IS: C C H(F) = ( (Z)C0 + C1 ) / ( (Z)B0 + B1 ) C WHERE: C Z=EXP(2PIF), P=PI, I=SQRT(-1), F=FREQ IN RAD/SEC C C0=W, C1=W, C B0=W+1, B1=W-1, C W=TAN(PFCDT), FC=CUTOFF FQ, DT=SAMPLE RATE C C THE RECURSIVE EQUATION FOR EACH 2-POLE SECTION IS: C C Y(N)B0 + Y(N-1)B1 = X(N)C0 + X(N-1)C1 C C WHICH BECOMES C C Y(N) = ( X(N)C0 + X(N-1)C1 - Y(N-1)B1 ) / B0 C C C INPUT ARGUMENTS: C FCUT - REAL*8. THE CUT-OFF FREQUENCY IN NORMALIZED UNITS C OF 1/SAMPLES. FCUT = FREQ (HZ) * SAMPLE-RATE (SEC). C NPOLE - INTEGER*4. THE NUMBER OF FILTER POLES IN THE LEFT C HALF OF THE COMPLEX S-PLANE. ALL OF THE R-C POLES C ARE ON THE NEGATIVE REAL AXIS. C NREV - INTEGER*4. CAUSES SELECTED POLES TO BE FILTERED IN C REVERSE FOR THE PURPOSE OF PHASE COMPENSATION. THE C POSSIBLE VALUES OF NREV ARE AS FOLLOWS: C 0 = NO POLES FILTERED IN REVERSE. C 1 = ODD NUMBERED POLES FILTERED IN REVERSE. C 2 = EVEN NUMBERED POLES FILTERED IN REVERSE. C 3 = ALL POLES FILTERED IN REVERSE. C NSAMP - INTEGER*4. THE NUMBER OF SAMPLES IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENTS: C YIO - REAL*8. AN ARRAY CONTAINING THE UNFILTERED INPUT C SIGNAL. UPON RETURN FROM THE SUBROUTINE, IT WILL C CONTAIN THE FILTERED OUTPUT SIGNAL. C C R-C EQUATIONS DERIVED FROM: C HAMMING, R.W., 1977, DIGITAL FILTERS: PRENTICE-HALL, INC., C ENGLEWOOD CLIFFS, NEW JERSEY, PAGES 180-195. C C SUBROUTINE RCLOW WRITTEN BY ROB BRACKEN, USGS, 23MAY91. C VAX/VMS VERSION 2.0 29MAY91. C HP-9000 SERIES 700/800 UNIX VERSION 2.0, 2JUN95. C subroutine rclow(fcut,npole,nrev,nsamp,yio) C C DECLARATIONS C C PASSED ARGUMENTS real*8 fcut integer*4 npole,nrev,nsamp real*8 yio(*) C C OTHER VARIABLES c real*8 ax,ax0,ax1,ay1, b0,b1, c0,c1 real*8 ax,ay1, b0,b1, c0,c1 real*8 wc, xin0,xin1 C C CONSTANTS real*8 pie parameter( pie = 3.1415 92653 58979 32384 62643 ) C C CHECK FOR NUMBERS OF POLES LESS THAN ONE C if(npole.le.0) goto 999 C C COMPUTE FILTER COEFFICIENTS C C WARP THE CUTOFF FREQUENCY wc=dtan(pie*fcut) C C FIND COEFFICIENTS IN NUMERATOR OF TRANSFER FUNCTION c0=wc c1=wc C C FIND COEFFICIENTS IN DENOMINATOR OF TRANSFER FUNCTION b0=wc+1.d0 b1=wc-1.d0 C C FIND TIME-DOMAIN FILTER COEFFICIENTS FOR EACH R-C POLE C AX0=C0/B0 C AX1=C1/B0 ax=c0/b0 ay1=b1/b0 C C FILTER EACH POLE C do 701 k=1,npole C C DETERMINE THE DIRECTION OF FILTERING C if(nrev.eq.0) then krev=0 else if(nrev.eq.3) then krev=1 else krev=k-(k/2)*2 if(nrev.eq.2) krev=1-krev endif C if(krev.eq.0) then C C SET INDICIES TO FILTER FORWARD FROM BEGINNING OF ARRAY i1=1 i0=2 ie=nsamp istep=1 else C C SET INDICIES TO FILTER IN REVERSE FROM END OF ARRAY i1=nsamp i0=nsamp-1 ie=1 istep=-1 endif C C FILTER THE STARTING SAMPLE POINT C if(nsamp.ge.1) then C KEEP UNFILTERED VALUE OF STARTING SAMPLE POINT xin1=yio(i1) C CALCULATE FILTERED VALUE OF STARTING SAMPLE POINT yio(i1)=ax*xin1 endif C C FILTER THE REST OF THE SAMPLE POINTS C do 702 m=i0,ie,istep C KEEP UNFILTERED VALUE OF CURRENT SAMPLE POINT xin0=yio(m) C CALCULATE FILTERED VALUE OF CURRENT SAMPLE POINT yio(m)= ax*(xin0 + xin1) - ay1*yio(m-istep) C ROTATE ADDRESSES OF PREVIOUS UNFILTERED VALUES xin1=xin0 702 continue C 701 continue C C EXIT PROCEDURE C 999 return end C C________________________________________________________________ C C SUBROUTINE R I E S Z C________________________________________________________________ C C SUBROUTINE RIESZ REMOVES OR APPLIES A GENERALIZED RIESZ WINDOW C IN AN EVENLY SPACED ARRAY OF REAL DATA. C C THE RIESZ WINDOW IS THE SIMPLEST CONTINUOUS POLYNOMIAL WINDOW. C IT EXHIBITS A DISCONTINUOUS FIRST DERIVATIVE AT THE BOUNDARIES C AND APPROACHES THE CENTER WITH A CONVEX BULGE AS A NORMAL C WINDOW SHOULD. (THE GENERALIZED RIESZ WINDOW IS CALCULATED C WITH THE SAME EQUATION AS A RIESZ BUT THE INTENSITY, ALPHA IS C ALLOWED TO BE A VALUE OTHER THAN THAT REQUIRED FOR A STRICT C RIESZ.) AS LONG AS ALPHA IS GREATER THAN 1, THE BULGE REMAINS C CONVEX BECOMING MORE PRONOUNCED AS ALPHA INCREASES. IF ALPHA C IS LESS THAN 1, THE BULGE BECOMES CONCAVE AND THE FIRST C DERIVATIVE BECOMES CONTINOUS AT THE BOUNDARIES AND DISCONTINOUS C IN THE CENTER, A FAIRLY USELESS CONDITION. TECHNICALLY, IT C IS NOT A RIESZ WINDOW UNLESS ALPHA=2. FOR EXAMPLE, IF ALPHA=1, C THE RIESZ BECOMES A TRIANGLE (FEJER OR BARTLET) WINDOW. C C THE BASIC EQUATION IS C C W(J)=1-(2J/M)**ALPHA C C C IF ALPHA=0, THE ENTIRE WINDOW BECOMES ZERO. BUT IF THE WINDOW C IS BEING SPLIT (SEE WFRAC IN SUBROUTINE WINDOW) WITH ALPHA=0, A C RECTANGLE (BOX CAR) WITH ZERO EXTENSIONS IS PRODUCED. REMOVING C THE RIESZ WINDOW WITH ALPHA=0 IS NOT MEANINGFUL AND MAY RESULT C IN ERRORS. C C IF ALPHA IS MADE NEGATIVE, THIS SUBROUTINE PRODUCES A REFLECTED C RIESZ WINDOW. A REFLECTED RIESZ IS A SIMPLE RIESZ WINDOW C MODIFIED AS FOLLOWS: 1) THE LEFT HALF IS REFLECTED THROUGH THE C ORIGIN, 2) THE RIGHT HALF IS REFLECTED THROUGH THE LAST POINT, C 3) THE NEW WINDOW IS MOVED RIGHT HALF OF THE OLD LENGTH TO PUT C THE NEW BEGINNING POINT AT ZERO, 4) THE NEW LENGTH IS SCALED BY C HALF TO PUT THE NEW ENDING POINT AT ITS OLD LOCATION, 5) A C VALUE OF ONE (1) IS ADDED TO THE NEW WINDOW TO MAKE THE C BEGINNING AND ENDING VALUES EQUAL TO ZERO, AND 6) THE NEW C WINDOW IS SCALED BY HALF TO MAKE THE CENTER VALUE EQUAL TO ONE C (1). THE REFLECTION PROCESS RESULTS IN A CONTINUOUS FIRST C DERIVATIVE AT THE BOUNDARIES. C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE INTENSITY OF THE WINDOW. C IF ALPHA IS 2, A RIESZ WINDOW IS PRODUCED; ALPHA C GREATER THAN 1 PRODUCES A SIMILAR WINDOW; ALPHA C EQUAL TO 1 PRODUCES A TRIANGLE (FEJER OR BARTLETT) C WINDOW; ALPHA LESS THAN 1 PRODUCES A USELESS WINDOW; C ALPHA=0 IS A ZERO WINDOW; -ALPHA PRODUCES A C REFLECTED RIESZ WINDOW THE INTENSITY OF WHICH IS C CONTROLLED BY ABS(-ALPHA). C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE RIESZ WRITTEN BY ROB BRACKEN, USGS, 19970602. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 19970602. C C subroutine riesz(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp real*8 winval,xr,xctr,scale,xj C C INITIALIZATIONS C C FIND POINT WITH HIGHEST INDEX ON INTERVAL 1 TO N C DEPENDING ON WHETHER SYMMETRIC OR DFT EVEN nwp=nio if(kdft.ge.1) nwp=nio+1 C C SLIDE INTERVAL LEFT 1 POINT TO BE 0 TO N-1 nwp=nwp-1 xr=dflotj(nwp) C C FIND VALUE OF X CORRESPONDING TO THE CENTER OF THE PROFILE xctr=0.5d0*xr C C FIND SCALE FACTOR WHICH EQUALS 1 WHEN MULTIPLIED BY XCTR. scale=1.d0/xctr C C FIND AND APPLY (OR REMOVE) THE WINDOW VALUES C C APPLY (OR REMOVE) ALL EXCEPT BEG PT (& END PT IF SYMMETRIC) do j=1,jidint(xctr) xj=dflotj(j) C if(alpha.ge.0.d0) then C C GENERALIZED RIESZ WINDOW winval=1.d0-(scale*(xctr-xj))**alpha else C C REFLECTED GENERALIZED RIESZ WINDOW if(xj.gt.xctr/2.d0) then winval=1.d0-0.5d0*(2.d0*scale*(xctr-xj))**-alpha else winval=0.5d0*(2.d0*scale*xj)**-alpha endif endif C C DETERMINE WHETHER TO APPLY OR REMOVE THE WINDOW if(nrem.ge.1) winval=1.d0/winval C C APPLY (OR REMOVE) THE LEFT HALF OF THE WINDOW yio(j+1)=yio(j+1)*winval C C APPLY (OR REMOVE) THE RIGHT HALF OF THE WINDOW yio(nwp-j+1)=yio(nwp-j+1)*winval enddo C C APPLY (OR REMOVE) FIRST POINT (AND LAST POINT IF SYMMETRIC) if(nrem.le.0) then yio(1)=0.d0 if(kdft.le.0) yio(nio)=0.d0 else yio(1)=yio(2) if(kdft.le.0) yio(nio)=yio(nio-1) endif C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE T R E N D R C________________________________________________________________ C C SUBROUTINE TRENDR REMOVES OR ADDS BACK A TREND FROM AN C EVENLY-SPACED DOUBLE-PRECISION FUNCTION ARRAY. THE TYPE OF C TREND MAY BE SELECTED USING KTYPE. C C INPUT ARGUMENTS: C KTYPE - INTEGER*4. THE TYPE OF TREND. OPTIONS ARE: C * 0 = NO TREND C * 1 = POLYNOMIAL OF ORDER NCOEF-1 C 2 = COS FNC WITH AMPS, FQS, AND PHASES IN ACOEF C NCOEF - INTEGER*4. NUMBER OF COEFS IN ACOEF. IF NCOEF IS C MADE NEGATIVE AND NVERSE <= 0, INSTEAD OF FINDING C NEW COEFS, THE COEFS GIVEN IN ACOEF WILL BE USED FOR C REMOVING THE TREND. IF NCOEF IS NEGATIVE AND NVERSE C => 1, THERE WILL BE NO EFFECT DIFFERENT FROM THAT OF C POSITIVE NCOEF (I.E. THE COEFS IN ACOEF WILL BE C USED FOR ADDING BACK THE TREND). C C NOTE: CURRENTLY, NCOEF MAY NOT EXCEED 3 WHEN KTYPE C IS 1. C C ACOEF - REAL*8. THE COEFS DESCRIBING THE TREND FUNCTION. C IF THE TREND IS BEING REMOVED, ACOEF IS AN C INPUT/OUTPUT ARGUMENT. IF THE TREND IS BEING ADDED C BACK, ACOEF IS AN INPUT ARGUMENT. C C NVERSE - INTEGER*4. DETERMINES WHETHER TO REMOVE OR ADD BACK C THE TREND: C NVERSE <= 0 MEANS TO REMOVE THE TREND. C NVERSE => 1 MEANS TO ADD BACK THE TREND. C NPTS - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YFNC. C IF NPTS IS MADE NEGATIVE AND NVERSE <= 0, THE TREND C WILL NOT BE REMOVED (THE COEFS CAN THEN BE FOUND C WITHOUT CHANGING THE DATA ARRAY). IF NPTS IS C NEGATIVE AND NVERSE => 1, THE TREND WILL NOT BE C ADDED BACK (THIS HAS NO PARTICULAR USE BUT IS C PROVIDED FOR FUNCTIONAL BALANCE). C INPUT/OUTPUT ARGUMENT: C YFNC - REAL*8. ARRAY OF DIMENSION NPTS CONTAINING THE C FUNCTION TO HAVE THE TREND REMOVED OR ADDED BACK. C C C SUBROUTINE TRENDR WRITTEN BY ROB BRACKEN, USGS, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.1, 19970611. C ADDITION OF OPTION TO FIND COEFS BUT NOT REMOVE THE TREND C C subroutine trendr(ktype,ncoef,acoef, & nverse,npts,yfnc) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 ktype,ncoef,nverse,npts C C INPUT/OUTPUT ARGUMENTS real*8 acoef(*),yfnc(*) C C SELECT THE TYPE OF TREND C goto (101,102),ktype goto 990 C C TRENDS C C POLYNOMIAL 101 continue call lsqply(ncoef,acoef, nverse,npts,yfnc) goto 990 C C COSINE 102 continue call errfor(ktype,'cosine trend not installed') goto 990 C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE T U K E Y C________________________________________________________________ C C SUBROUTINE TUKEY REMOVES OR APPIES A TUKEY WINDOW IN AN EVENLY C SPACED ARRAY OF REAL DATA. THE TUKEY WINDOW IS THE SAME AS A C RECTANGLE THAT HAS BOTH ENDS TAPERED WITH A COSINE. IT IS C SOMETIMES CALLED A COSINE TAPER. C C THE TRANSFORM OF THE TUKEY WINDOW HAS A "CONFUSING ARRAY OF C LARGE SIDE LOBES" AND HENCE, IT IS A VERY POOR WINDOWING C FUNCTION; ITS USE IS NOT RECOMMENDED FOR DISCERNING PEAKS IN C THE FREQUENCY DOMAIN. C C HOWEVER, THE TUKEY CAN BE HELPFUL FOR KNOCKING DOWN STEPS AT C THE ENDS OF FUNCTIONS WHICH ARE GOING TO HAVE A TIME-DOMAIN C FILTER APPLIED. TESTS HAVE SHOWN THAT IF THE PROFILE IS FIRST C EXTENDED (WITH A BURG PREDICTION) AND THEN A TUKEY WINDOW IS C APPLIED IN WHICH THE TAPER IS THE SAME LENGTH AS THE EXTENSION, C RINGING CAN BE REDUCED BY UP TO 20 DB. (NOTE THAT WITH THIS C WINDOW, RINGING MAY BE PARTIALLY RE-ENERGIZED NEAR THE C BEGINNING AND END OF EACH TAPER). IF THE TAPER IS HALF THE C LENGTH OF THE EXTENSION, SOMETIMES A GREATER REDUCTION IN C RINGING CAN BE ACHIEVED. C C KEEP IN MIND THAT A STEP AT EITHER END OF A PROFILE WILL RING C IN PROPORTION TO ITS SIZE; UPON UNWINDOWING, THE RINGING FROM A C STEP WILL BE AMPLIFIED BACK TO THE SIZE (OR NEARLY SO) IT WOULD C HAVE BEEN IF NO WINDOWING HAD BEEN PERFORMED. THEREFORE, IF A C WINDOW IS USED TO REDUCE RINGING, THE WINDOWING SHOULD BE C RESTRICTED TO EXTENDED PORTIONS WHICH ARE TO BE REMOVED BEFORE C FURTHER PROCESSING. C C C INPUT ARGUMENTS: C ALPHA - REAL*8. DETERMINES THE PERCENTAGE OF THE RECTANGLE C OVERPRINTED WITH COSINE TAPERS. THE RANGE OF VALUES C IS 0.0 TO 1.0. FOR EXAMPLE, IF ALPHA=0.4, 6/10 OF C THE WINDOW WIDTH CENTERED WILL HAVE A VALUE OF 1; C THE LEFT 2/10 WILL TAPER TO ZERO AND THE RIGHT 2/10 C WILL TAPER TO ZERO. IF ALPHA IS GREATER THAN 1.0, C SUBROUTINE HANN3 WILL BE CALLED AUTOMATICALLY WITH C THE SAME ARGUMENTS AS TUKEY. THIS IS DONE BECAUSE, C A HANNING WINDOW IS THE LOGICAL ANALOGUE TO A TUKEY C WINDOW WITH TAPERS MEETING IN THE MIDDLE. C NREM - INTEGER*4. DETERMINES WHETHER TO REMOVE (NREM = 1) C OR APPLY (NREM = 0) THE WINDOW. C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT = 1) OR SYMMETRIC (KDFT = 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NIO - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YIO. C C INPUT/OUTPUT ARGUMENT: C YIO - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE TUKEY WRITTEN BY ROB BRACKEN, USGS, 30MAY97. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 30MAY97. C C subroutine tukey(alpha,nrem,kdft,nio,yio) C C DECLARATIONS C C INPUT ARGUMENTS real*8, alpha integer*4 nrem,kdft,nio C C INPUT/OUTPUT ARGUMENT real*8 yio(*) C C INTERNAL VARIABLES integer*4 nwp real*8 pi,winval,xr,xpi,scale parameter(pi = 3.1415 92653 58979 32384 62643) C C CHECK VALUE OF ALPHA C if(alpha.gt.1.d0) then C C THE HANNING WINDOW IS THE NATURAL EXTENSION OF A TUKEY call hann3(alpha,nrem,kdft,nio,yio) goto 990 endif C C INITIALIZATIONS C C FIND POINT WITH HIGHEST INDEX ON INTERVAL 1 TO N C DEPENDING ON WHETHER SYMMETRIC OR DFT EVEN nwp=nio if(kdft.ge.1) nwp=nio+1 C C SLIDE INTERVAL LEFT 1 POINT TO BE 0 TO N-1 nwp=nwp-1 xr=dflotj(nwp) C C FIND THE VALUE OF X CORRESPONDING TO COS(PI) IN LEFT TAPER xpi=(0.5d0*alpha)*xr C C FIND SCALE FACTOR WHICH EQUALS PI WHEN MULTIPLIED BY XPI. scale=pi/xpi C C FIND AND APPLY THE WINDOW VALUES C C APPLY ALL EXCEPT FIRST POINT (AND LAST POINT IF SYMMETRIC) do j=1,jidint(xpi) winval=0.5d0*(1.d0-dcos(scale*dflotj(j))) if(nrem.ge.1) winval=1.d0/winval yio(j+1)=yio(j+1)*winval yio(nwp-j+1)=yio(nwp-j+1)*winval enddo C C APPLY FIRST POINT (AND LAST POINT IF SYMMETRIC) if(nrem.le.0) then yio(1)=0.d0 if(kdft.le.0) yio(nio)=0.d0 else yio(1)=yio(2) if(kdft.le.0) yio(nio)=yio(nio-1) endif C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE U N W I N 0 C________________________________________________________________ C C SUBROUTINE UNWIN0 "UNWINNOWS" A PROFILE OF DATA. C C WINNOWING IS SIMPLY THE PROCESS OF REMOVING INTERVENING DATA C POINTS AND COMPACTING THE DATA PROFILE. SO TECHNICALLY, C "UNWINNOWING" SHOULD PERFORM THE EXACT INVERSE OF WINNOWING, C THAT IS, RESTORE THE ORIGINAL PROFILE. HOWEVER, THE WINNOWING C PROCESS ESSENTIALLY ERASES THE HIGHER FREQUENCIES MAKING THEM C UNAVAILABLE FOR RESTORING. THEREFORE, THE "UNWINNOWING" AS C PERFORMED BY THIS SUBROUTINE RESTORES THE ORIGINAL SAMPLING C INTERVAL AND NUMBER OF POINTS BUT WITH THE SPECTRUM OF ONLY THE C WINNOWED PROFILE. C C IN THE UNWINNOWING PROCESS, ALL OF THE WINNOWED POINTS ARE C RESTORED TO THE ORIGINAL PROFILE WITH THEIR EXACT WINNOWED C VALUES. THE FOURIER TRANSFORM METHOD (WHICH IS EMPLOYED BY C THIS SUBROUTINE) INSURES THROUGH ORTHOGONALITY THAT THESE C POINTS WILL IN FACT BE RESTORED. HOWEVER, THE INTERVENING C POINTS MUST BE INTERPOLATED AND THEIR VALUES DEPEND UPON C SELECTION OF EXTENSION PROCESSES AND WINDOWING METHODS. THESE C HAVE BEEN OPTIMIZED FOR THIS SUBROUTINE BUT VERY SMALL ERRORS C REMAIN. TYPICALLY, MID-PROFILE ERRORS SHOULD BE WELL BELOW -72 C DB REFERENCED TO THE PROFILE VARIANCE. END POINT ERRORS C (WITHIN THE FIRST AND LAST 2*KWIN POINTS) HAVE BEEN FOUND AS C LARGE AS -32 DB. C C THE FREQUENCY SPECTRUM OF THE UNWINNOWED PROFILE MATCHES VERY C WELL TO THE ORIGINAL PROFILE IN THE LOWER HALF OF THE WINNOWED C NYQUIST INTERVAL. OF COURSE, THE UPPER HALF OF THE WINNOWED C NYQUIST INTERVAL MAY DEPART SIGNIFICANTLY FROM THE ORIGINAL C BECAUSE OF THE LOW-PASS FILTERING WHICH WOULD HAVE OCCURED C BEFORE WINNOWING. C C C INPUT ARGUMENTS: C KWIN - INTEGER*4. THE WINNOW FACTOR, AN INTEGER FACTOR BY C WHICH THE ORIGINAL SAMPLING INTERVAL WAS MULTIPLIED C VIA WINNOWING. THE FIRST POINT OF THE ORIGINAL C PROFILE WAS KEPT, AS WELL AS EVERY KWINTH POINT C AFTERWARDS. FOR EXAMPLE, IF KWIN=3, POINTS 1,2,3, C ... N OF THE WINNOWED PROFILE WILL HAVE COME FROM C POINTS 1,4,7, ... 3N+1 OF THE ORIGINAL PROFILE. THE C TOTAL NUMBER OF POINTS IN THE WINNOWED PROFILE WILL C EQUAL 1+(NPTS-1)/KWIN (INTEGER DIVISION ASSUMED). C IF UNWINNOWING IS TO BE PERFORMED, ABS(KWIN) MUST BE C GREATER THAN 1. C NPTS - INTEGER*4. THE ORIGINAL, UNWINNOWED NUMBER OF C POINTS IN ARRAY YFNC. C C NOTE: THE NUMBER OF WINNOWED POINTS IS ALWAYS C 1+(NPTS-1)/KWIN. THIS VALUE IS NOT GIVEN AS AN C ARGUMENT BECAUSE ITS VALUE CAN ALWAYS BE CALCULATED C EXACTLY FROM KWIN AND NPTS; BUT, THE VALUE OF NPTS C CANNOT BE CALCULATED EXACTLY FROM KWIN AND THE C WINNOWED NUMBER OF POINTS. C C INPUT/OUTPUT ARGUMENTS: C YFNC - REAL*8. A SINGLE DIMENSIONED ARRAY CONTAINING C 1+(NPTS-1)/KWIN DATA POINTS TO BE UNWINNOWED. UPON C OUTPUT, YFNC WILL CONTAIN NPTS OF UNWINNOWED DATA. C YFNC DOUBLES AS A WORKSPACE DURING THE UNWINNOWING C PROCESS AND THE REQUIRED WORKSPACE (THE DIMENSION OF C YFNC) MAY EXCEED NPTS. TO CALCULATE THE NECESSARY C DIMENSION, SEE NDIM. C NDIM - INTEGER*4. THE DIMENSION OF YFNC. NDIM CAN BE C GIVEN EITHER OF TWO VALUES, ZERO (0) OR A POSITIVE C VALUE THAT IS EQUAL TO OR GREATER THAN THE LARGEST C DIMENSION REQUIRED FOR YFNC. IF THE POSITIVE VALUE C IS GIVEN, PROCESSING WILL OCCUR. BUT IF THE VALUE C IS LESS THAN THE REQUIRE SPACE, AN ERROR WILL BE C GENERATED. IF THE ZERO VALUE IS GIVEN, NO C PROCESSING WILL OCCUR DURING THAT CALL, BUT RATHER C THE MINIMUM REQUIRED DIMENSION OF YFNC WILL BE C CALCULATED AND RETURNED IN NDIM. C C OUTPUT ARGUMENTS: C FDH - COMPLEX*16. WORK ARRAY OF OF DIMENSION C (NDIM/KWIN-1)/2 WHEN NDIM HAS BEEN GIVEN A VALUE OF C 0 AND UNWIN0 ALLOWED TO CALCULATE A LENGTH. FDH IS C A HALF SPECTRUM (SEE SUBROUTINES DFTFTR.F AND C CNVSPC.F) OF THE WINNOWED, EXTENDED PROFILE. THE C COMPLEX*16 HALF SPECTRUM IS DESIGNED TO REQUIRE NO C MORE REAL*8 ELEMENTS THAN ITS REAL*8 TIME-DOMAIN C COUNTERPART. IF NECESSARY, FDH MAY BE COMBINED WITH C FDAVG AND FDNYQ TO RECONSTRUCT A FULL SPECTRUM OF C THE WINNOWED EXTENDED PROFILE. HOWEVER, FOR THE C C PURPOSES OF UNWINNOWING, FDH IS SIMPLY A WORKSPACE. C FDAVG - REAL*8. WORK VARIABLE WHICH IS THE ZERO FREQUENCY C PART OF THE HALF SPECTRUM AS DESCRIBED ABOVE AND IN C SUBROUTINE DFTFTR.F. C FDNYQ - REAL*8. WORK VARIABLE WHICH IS THE NYQUIST C FREQUENCY PART OF THE HALF SPECTRUM AS DESCRIBED C ABOVE AND IN SUBROUTINE DFTFTR.F. C C C SUBROUTINE UNWIN0 WRITTEN BY ROB BRACKEN, USGS, 19970506. C HP-9000 SERIES 700/800 VERSION 1.0, 19970506 C C subroutine unwin0(kwin,npts,yfnc,ndim,fdh,fdavg,fdnyq) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 kwin,npts C C INPUT/OUTPUT ARGUMENTS real*8 yfnc(*) integer*4 ndim C C OUTPUT ARGUMENTS complex*16 fdh(*) real*8 fdavg,fdnyq C C UNWINNOWING AND DFT VARIABLES integer*4 nptw,kfac real*8 signi C C LEAST SQUARES LINE REMOVAL integer*4 ktype,ncoef,nverse real*8 slope,rcept real*8 acoef(2) C C FREQUENCY DEPENDENT VARIANCE real*8 varf C C EXTENSION-AMOUNT CALCULATION real*8 rdb,profdb,thresh,step real*8 tenpct integer*4 jleft,kright,jmn,jlefw C C EXTENSION PROCEDURES real*8 tparms(3) integer*4 just,npxt,npxw C C WINDOWING FUNCTION real*8 alpha,wfrac integer*4 kdft C C BUTTERWORTH-FILTER PARAMETERS C C KSQ=SQUARED FILTER (TO MAKE ZERO PHASE), 0=OFF, 1=ON C RINGDB=FILTER END-EFFECT RINGING IN DECIBELS C WCRNYQ=CUTOFF FREQUENCY RELATIVE TO NYQUIST C WSRWC=STOP BAND (-36 DB UNSQRD, -72 DB SQRD) FQ REL TO FCUT C (STOP BAND FREQ AT 1.29*FCUT PRODUCES A 16-POLE FILTER) integer*4 ksq real*8 ringdb,wcrnyq,wsrwc parameter(ringdb=-72.d0,ksq=1) parameter(wcrnyq=0.8053d0,wsrwc=1.2968d0) C real*8 fcut integer*4 npole C C SELECT SPLIT RIESZ WINDOW (KRIESZ=1) OR KAISER (KRIESZ=0) integer*4 kriesz parameter(kriesz=0) C C DETERMINE WHETHER TO UNWINNOW C if(jiabs(kwin).le.1) goto 990 C C C UNWINNOW C C CONVERT THE WINNOW FACTOR kfac=jiabs(kwin) C C FIND THE NUMBER OF POINTS IN THE WINNOWED PROFILE nptw=1+(npts-1)/kfac C C FIND WINNOWED CUTOFF FREQ REFERENCED TO UNWINNOWED NYQUIST fcut=0.5d0*wcrnyq/dble(kfac) npole=16 C C FIND COEFS OF A LEAST SQUARES LINE LATER TO BE REMOVED C if(ndim.ge.1.or.ringdb.ge.0.d0) then ktype=1 ncoef=2 nverse=0 call trendr(ktype,ncoef,acoef,nverse,-nptw,yfnc) endif C C FIND THE VARIANCE (POWER) OF THE WINNOWED PROFILE C (SAME AS THE FREQ-DEPENDENT VARIANCE OF THE UNWINNOWED PROFILE) C if(ringdb.lt.0.d0) then rdb=ringdb C C RINGDB IS A RELATIVE THRESHOLD; EXTENSION LENGTHS WILL BE C CALCULATED BY BWRING SUCH THAT POWER IS RELATIVE TO THE C PROFILE VARIANCE. ADJUSTMENT OF THE REQUESTED RINGING C POWER IS UNNECESSARY. profdb=0.d0 else rdb=-ringdb C C RINGDB IS AN ABSOLUTE THRESHOLD; BECAUSE EXTENSION C LENGTHS ARE CALCULATED BY BWRING RELATIVE TO PROFILE C VARIANCE, THE EFFECT OF VARIANCE MUST BE REMOVED FROM THE C REQUESTED RINGING POWER IN ORDER TO TIE IT TO THE C ABSOLUTE STANDARD (VARIANCE = 1/2). C C (NOTE: INTERCEPT OCCURS AT I=0 NOT I=1) slope=acoef(2) rcept=acoef(1) varf=0.d0 do i=1,nptw varf=varf+( yfnc(i)-(slope*dflotj(i)+rcept) )**2.d0 enddo varf=varf/dflotj(nptw) profdb=10.d0*dlog10(2.d0*varf) endif C C FIND AMOUNT THE UNWINNOWED PROFILE SHOULD HAVE BEEN EXTENDED C C FIND A THE RINGING THRESHOLD thresh=rdb-profdb C C CONVERT THRESH TO AN AMPLITUDE FOR SUBR BWRING thresh=10.d0**(thresh/20.d0) step=1.d0 call bwring(thresh,step, fcut,npole,ksq, jleft,kright) C C ADJUST JLEFT TO EXTEND FAR ENOUGH THAT THE ENDS OF THE C PROFILE ARE HIGHER THAN 10% ON A KAISER WINDOW WITH ALPHA=3 tenpct=0.1632d0 jmn=dflotj(3*kfac+npts)*tenpct/(1.d0-2.d0*tenpct)+1.d0 if(jleft.lt.jmn) jleft=jmn C C ADJUST JLEFT TO THE WINNOWED & UNWINNOWED EXTENSION AMOUNTS jlefw=1+(jleft-1)/kfac jleft=jlefw*kfac C C FIND THE WINNOWED AND UNWINNOWED NUMBER OF POINTS EXTENDED npxw=1+(jleft+npts+jleft-1)/kfac npxt=npxw*kfac C C CHECK NDIM if(ndim.le.0) then C C CALCULATE THE VALUE OF NDIM AND RETURN TO CALLING ROUTINE if(kriesz.eq.1) then C C SPLIT RIESZ WINDOW INTERNALLY SELECTED ndim=jmax0(npts,npxw) else C C KAISER WINDOW INTERNALLY SELECTED ndim=npxt endif goto 990 endif C C REMOVE AN LSQ LINE FROM THE WINNOWED PROFILE USING COEFS C CALCULATED EARLIER IN THIS SUBR C ktype=1 ncoef=2 nverse=0 call trendr(ktype,-ncoef,acoef,nverse,nptw,yfnc) C C EXTEND THE PROFILE WITH A BURG PREDICTION FILTER C C NOTE: THE EXTENSION USED HAS EVERYTHING TO DO WITH HOW WELL C THE RESTORED WINNOWED POINTS NEAR THE PROFILE ENDS WILL MATCH C THE ORIGINAL. HOWEVER, MATCHING THE ORIGINAL IS A BIT C SUBJECTIVE BECAUSE FREQUENCIES ARE LOST GOING THROUGH THE C WINNOWING STAGE. GENERALLY, THE LONGER ROOT IN THE BURG C EXTENSION IS BETTER. BUT, THERE IS A TRADE-OFF WITH TIME ... C AND THE BURG EXTENSION IS VERY SLOW. BECAUSE OF THIS, THE C COMPROMISE HAS BEEN TAKEN TO MAKE THE ROOT LENGTH EQUAL THE C EXTENSION LENGTH. C C BE SURE THE EXTENSIONS WILL FIT IN THE ALOTTED DIMENSION if(ndim.lt.npxw) & call errfor(npxw, & '(unwin0) dimension ndim of array yfnc is too small') C C SELECT BURG PREDICTION FILTER FOR EXTENSION ktype=7 C C DEPTH OF BURG FILTER COEFS = DEPTH OF EXTENSION tparms(1)=-1.d0 C DEPTH OF BURG FILTER COEFS = LENGTH OF PROFILE c tparms(1)=0.d0 C DEPTH OF BURG FILTER COEFS = VALUE OF TPARMS(2) AND (3) c tparms(1)=-3.d0 c tparms(2)=2048.d0 c tparms(3)=2048.d0 just=jlefw+1 call extend(ktype,tparms,just,nptw,npxw,yfnc) C C WINDOW THE PROFILE WITH EITHER A SPLIT RIESZ OR KAISER WINDOW C C NOTE: WINDOWING WITH A SPLIT RIESZ PRODUCES ABOUT 10 DB MORE C NOISE THAN THE KAISER THROUGHOUT THE PROFILE. HOWEVER, THE C SPLIT RIESZ WINDOW REMOVES THE DANGER OF DIFFERING WEIGHTS IN C VARIOUS PARTS OF THE PROFILE. SELECTION OF THE WINDOWING C FUNCTION DOES NOT EFFECT THE TIME-DOMAIN ENDS. C if(kriesz.eq.1) then C C WINDOWING WITH A SPLIT RIESZ WINDOW FUNCTION (ALPHA=3.0) ktype=4 alpha=2.5d0 wfrac=dflotj(jleft)/dflotj(npxt) else C C WINDOWING WITH A KAISER WINDOW FUNCTION (ALPHA=3.0) ktype=14 alpha=3.0d0 wfrac=1.d0 endif kdft=1 nverse=0 call window(ktype,alpha,wfrac,kdft,nverse,npxw,yfnc) C C EXPAND THE WINNOWED EXTENDED PROFILE BACK TO THE UNWINNOWED C EXTENDED PROFILE WHILE IMPLICITLY REMOVING EXTENSIONS C C NOTE: DFTFDR IS THE SAME AS DFTFTR EXCEPT THAT IT CALCULATES C THE WHOLE TIME-DOMAIN PROFILE INCLUDING THE PARTS THAT ARE TO C BE REMOVED. WHEREAS DFTFTR SAVES TIME, IMPLICITLY REMOVING C THE EXTENSIONS BY NOT CALCULATING THEM. C C BE SURE THE UNWINNOWED PROFILE WILL FIT IN THE ALOTTED DIM if(ndim.lt.npts) & call errfor(npts, & '(unwin0) dimension ndim of array yfnc is too small') C kscale=1 signi=-1.d0 call dftftr(kscale, & 1,npxw,npxw,npxw,signi,yfnc,fdh,fdavg,fdnyq) c call dftfdr(npxw,npxw,signi,yfnc,fdh) jb=jleft+1 je=jleft+npts signi=+1.d0 call dftftr(kscale, & jb,je,npxt,npxw,signi,yfnc,fdh,fdavg,fdnyq) c call dftfdr(npxt,npxw,signi,yfnc,fdh) C C REMOVE THE FULL-LENGTH KAISER WINDOW C if(kriesz.ne.1) then if(ndim.lt.npxt) call errfor(ndim, & '(unwin0) ndim is less than npxt') ktype=0 just=jleft+1 call extend(ktype,tparms,just,npts,npxt,yfnc) C ktype=14 alpha=3.0d0 wfrac=1.d0 kdft=1 nverse=1 call window(ktype,alpha,wfrac,kdft,nverse,npxt,yfnc) ktype=0 just=jleft+1 call extend(ktype,tparms,just,npxt,npts,yfnc) endif C C UNWINNOW THE LSQ LINE AND REPLACE INTO UNWINNOWED PROFILE C C IMPLICITLY UNWINNOW THE LSQ LINE BY CHANGING THE COEFS slope=acoef(2) rcept=acoef(1) acoef(2)=slope/dble(kfac) acoef(1)=rcept+slope-acoef(2) C C REPLACE THE LSQ LINE ktype=1 ncoef=2 nverse=1 call trendr(ktype,ncoef,acoef,nverse,npts,yfnc) C C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE V A R F Q D C________________________________________________________________ C C SUBROUTINE VARFQD FINDS THE VARIANCE OF A PROFILE OF DATA ABOUT C A LINE (SUPPLIED BY THE CALLING PROGRAM). THE VARIANCE IS C CALCULATED FOR BOTH FILTERED (VARF) AND UNFILTERED (VARU) C VERSIONS OF THE DATA. FILTERING IS DONE IMPLICITLY (DOES NOT C AFFECT THE DATA PROFILE) BY A LOW-PASS BUTTERWORTH TIME-DOMAIN C FILTER WITH SELECTABLE CUTOFF FREQUENCY AND NUMBER OF POLES. C C THE VARIANCE OF THE FILTERED DATA IS THE TOTAL POWER CONTAINED C IN ALL FREQUENCIES BELOW CUTOFF; THE VARIANCE OF THE UNFILTERED C DATA IS THE TOTAL POWER IN THE PROFILE. THE POWER ABOVE C CUTOFF MAY BE FOUND BY SUBTRACTING THE POWER BELOW CUTOFF FROM C THE TOTAL POWER. C C IF A DATA PROFILE IS LATER TO BE EXTENDED AND THEN FILTERED, C THE FILTERED VARIANCE FROM THIS SUBROUTINE GIVES A STATISTICAL C MEASURE OF THE EXPECTED PEAK AMPLITUDE WHICH MAY OCCUR AT THE C ENDS OF THE PROFILE. THE EXPECTED PEAK AMPLITUDE IS NECESSARY C TO DETERMINE THE LENGTH THAT THE PROFILE MUST BE EXTENDED C BEFORE SUBSEQUENT FILTERING IN ORDER TO REDUCE RINGING IN THE C UNEXTENDED PART. THIS LENGTH CAN BE FOUND BY SUBROUTINE C BWRING. IF USING BWRING, THE STEP ARGUMENT CAN BE GIVEN THE C VALUE SQRT(2*VARF) WHICH IS THE EXPECTED PEAK AMPLITUDE C ESTIMATED FROM THE FILTERED VARIANCE. THE CUTOFF FREQUENCIES C SHOULD MATCH IN VARFQD, BWRING, AND THE FILTER SUBROUTINE. C (THE NUMBERS OF POLES SHOULD MATCH IN BWRING AND THE FILTER C SUBROUTINE.) C C C FOLLOWING IS A BRIEF DESCRIPTION OF THE BUTTERWORTH FILTER USED C BY THIS SUBROUTINE: C C THE TRANSFER FUNCTION OF THE BUTTERWORTH FILTER HAS POLE PAIRS C ON THE LEFT COMPLEX PLANE STARTING NEAR THE POSITIVE AND C NEGATIVE IMAGINARY AXIS AND PROGRESSING TOWARD THE LEFT REAL C AXIS. IF THE NUMBER OF POLES IS ODD, THE FINAL (ODD) POLE IS C ON THE LEFT REAL AXIS AND BECOMES A SIMPLE RC POLE. THE TIME C DOMAIN OPERATION INVOLVES RUNNING THE ENTIRE LENGTH OF THE C PROFILE WITH EACH 2-POLE SECTION (POLE PAIR) STARTING NEAR THE C IMAGINARY AXIS AND PROGRESSING TO THE LEFT REAL AXIS JUST AS IN C THE TRANSFER FUNCTION. C C THE TIME-DOMAIN BUTTERWORTH FILTER HAS A SIGNAL DELAY RESULTING C FROM VARIOUS PHASE SHIFTS AT VARIOUS FREQUENCIES. THESE DELAYS C SEEM TO SETTLE OUT TO A CONSTANT TIME DELAY BELOW CUTOFF C FREQUENCY WHICH MAY BE DESCRIBED (EMPIRICALLY) AS 37.071 C DEGREES PER POLE AT CUTOFF. NOTE THAT FOR PURPOSES OF C CALCULATING DELAY, THE CUTOFF FREQUENCY IS WARPED INTO THE C INTERVAL (ZERO,INFINITY) BY TAN(PI*FCUT)/PI. THE DELAY MAY BE C ACCOUNTED FOR BY APPLYING AN APPROPRIATE CORRECTION BASED ON C THIS INFORMATION. ALTERNATIVELY, THE PHASE SHIFTS (AND HENCE C THE DELAY) MAY BE ELIMINATED ENTIRELY BY RUNNING THE IDENTICAL C FILTER BOTH DIRECTIONS, THUS SQUARING THE TRANSFER FUNCTION. C THE PRIMARY EFFECTS OF THIS OPERATION ARE: REMOVE PHASE C SHIFTS, DOUBLE THE EFFECTIVE NUMBER OF POLES (DOUBLE THE SLOPE C IN THE STOP BAND), AND DECREASE THE POWER AT CUT-OFF TO -6 DB. C C THE TRANSFER FUNCTION FOR EACH 2-POLE SECTION IS: C C H(F) = ( (ZZ)C0 + (Z)C1 + C2 ) / ( (ZZ)B0 + (Z)B1 + B2 ) C WHERE: C Z=EXP(2PIF), P=PI, I=SQRT(-1), F=FREQ IN RAD/SEC C C0=WW, C1=2WW, C2=WW C B0=WW+1-2WCOS(X), B1=2WW-2, B2=WW+1+2WCOS(X) C W=TAN(PFCDT), FC=CUTOFF FQ, DT=SAMPLE RATE C X=ANGLE OF THE UPPER POLE IN THE CONJUGATE POLE-PAIR C C THE RECURSIVE EQUATION FOR EACH 2-POLE SECTION IS: C C Y(N)B0 + Y(N-1)B1 + Y(N-2)B2 = C X(N)C0 + X(N-1)C1 + X(N-2)C2 C C WHICH BECOMES C C Y(N) = C ( X(N)C0 + X(N-1)C1 + X(N-2)C2 - Y(N-1)B1 - Y(N-2)B2 ) / B0 C C BUTTERWORTH EQUATIONS DERIVED FROM: C C HAMMING, R.W., 1977, DIGITAL FILTERS: PRENTICE-HALL, INC., C ENGLEWOOD CLIFFS, NEW JERSEY, PAGES 180-195. C C DAVID E. RYERSON, 1985, THE DESIGN AND USE OF DIGITAL C FILTERS WITH PERSONAL COMPUTER SPREADSHEET PROGRAMS, SANDIA C NATIONAL LABORATORIES, ALBUQUERQUE, NM, NTIS REPORT C DE85013805, SAND-85-1073, 55P. C C C INPUT ARGUMENTS: C FCUT - REAL*8. THE CUT-OFF FREQUENCY IN NORMALIZED UNITS C OF 1/SAMPLES WHERE NYQUIST IS 0.5. FCUT = FREQ (HZ) C * SAMPLE-RATE (SEC). FCUT MAY RANGE FROM 0. TO 0.5 C (1/SAMPLES). C NPOLE - INTEGER*4. THE NUMBER OF BUTTERWORTH FILTER POLES C TO APPLY. IF NPOLE IS 0 OR NEG, NO FILTERING WILL C BE DONE. IF THE VARIANCE OF A SQUARED BUTTERWORTH C TRANSFER FUNCTION IS DESIRED, SIMPLY DOUBLE THE C NUMBER OF POLES. HOWEVER, FOR THIS CALCULATION AN C ABRUPT CUTOFF IS NOT REALLY NECESSARY AS IT HAS C LITTLE EFFECT ON THE VARIANCE. THEREFORE, IT IS C BETTER TO MAXIMIZE SPEED BY USING FEWER POLES. C ALTHOUGH THIS ROUTINE IS SET TO HANDLE UP TO 64, C EIGHT POLES SHOULD BE SUFFICIENT FOR ANY C APPLICATION. C NSAMP - INTEGER*4. THE NUMBER OF SAMPLES IN ARRAY YIO. C YIO - REAL*8. AN ARRAY CONTAINING THE UNFILTERED INPUT C SIGNAL. C SLOPE - REAL*8. THE SLOPE OF THE LINE TO REMOVE IMPLICITLY C FROM YIO BEFORE FILTERING. IF THE VARIANCE ABOUT C ZERO OR A CONSTANT IS TO BE CALCULATED, SLOPE SHOULD C BE SET TO ZERO. C RCEPT - REAL*8. THE Y-INTERCEPT OF THE LINE TO REMOVE C IMPLICITLY FROM YIO BEFORE FILTERING. IF THE C VARIANCE ABOUT ZERO IS TO BE CALCULATED, RCEPT C SHOULD BE SET TO ZERO. C C OUTPUT ARGUMENTS: C VARF - REAL*8. THE VARIANCE OF THE FILTERED DATA ABOUT THE C LINE DEFINED BY SLOPE AND RCEPT. THE RMS AMPLITUDE C OF THE FILTERED DATA IS SQRT(VARF); THE EXPECTED C PEAK AMPLITUDE IS SQRT(2*VARF); AND THE EXPECTED C PEAK-TO-PEAK AMPLITUDE IS 2*SQRT(2*VARF). VARF IS C THE POWER CONTAINED IN ALL FREQUENCIES BELOW FCUT. C VARU - REAL*8. THE VARIANCE OF THE UNFILTERED DATA ABOUT C THE LINE DEFINED BY SLOPE AND RCEPT. VARU IS THE C TOTAL POWER CONTAINED IN THE PROFILE, YIO REFERENCED C TO THE LINE DETERMINED BY SLOPE AND RECEPT. C C C SUBROUTINE VARFQD WRITTEN BY ROB BRACKEN, USGS, 19970529. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 19970529. C C subroutine varfqd(fcut,npole,nsamp,yio,slope,rcept, & varf,varu) C C DECLARATIONS C C INPUT ARGUMENTS real*8 fcut integer*4 npole,nsamp real*8 yio(*) real*8 slope,rcept C C OUTPUT ARGUMENTS real*8 varf,varu C C POLE-COEF CONSTANTS real*8 b0,b0inv,b1,b2, c0,c1,c2 real*8 offset,polarc,polx, wc,wcsq,wcsq1,wctwo C C POLE-COEF ARRAYS integer*4 mxp1 parameter(mxp1=65) real*8 ax(mxp1),ay1(mxp1),ay2(mxp1),yx(mxp1,2) C C PI CONSTANT real*8 pi parameter( pi = 3.1415 92653 58979 32384 62643 ) C C C CHECK FOR MAXIMUM ALLOWED NUMBER OF POLES C if(npole.gt.mxp1-1) call errfor(mxp1-1, & '(varfqd) Max number of poles exceeded') C C INITIALIZE VARIANCE SUMS C varf=0.d0 varu=0.d0 C C ASSUME NO FILTERING IF ZERO OR NEGATIVE NUMBER OF POLES C if(npole.le.0) then do i=1,nsamp varu=varu+( yio(i)-(slope*dflotj(i)+rcept) )**2.d0 enddo varf=varu goto 102 endif C C C SET UP BUTTERWORTH FILTER COEFS, AX(), AY1(), AND AY2() C C WARP THE CUTOFF FREQUENCY wc=dtan(pi*fcut) wctwo=wc*2.d0 wcsq=wc*wc wcsq1=wcsq+1.d0 C C FIND COEFFICIENTS IN NUMERATOR OF TRANSFER FUNCTION c0=wcsq c1=wcsq*2.d0 c2=wcsq C C FIND THE MIDDLE COEFFICIENT IN DENOMINATOR OF TRANSFER FUNC b1=2.d0*(wcsq-1.d0) C C SET UP CONSTANTS FOR COMPUTING POLE PAIR LOCATIONS polarc=pi/npole offset=(pi-polarc)/2 C C FIND NUMBER OF BUTTERWORTH POLE PAIRS nppair=npole/2 c C FIND NUMBER OF RC POLES (1 OR 0) nrc=npole-nppair*2 C C FIND BUTTERWORTH POLE PAIR COEFS do k=1,nppair C C COMPUTE THE OTHER COEFFICIENTS IN DENOMINATOR OF XFER FCN polx=wctwo*dcos(polarc*k+offset) b0=(wcsq1-polx) b2=wcsq1+polx C C TIME-DOMAIN FILTER COEFS FOR EACH BUTTERWORTH POLE-PAIR b0inv=1.d0/b0 ax(k)=c0*b0inv ay1(k)=b1*b0inv ay2(k)=b2*b0inv enddo C C FIND RC-POLE COEFS, IF ANY if(nrc.eq.1) then ax(k)=wc/(wc+1.d0) ay1(k)=(wc-1.d0)/(wc+1.d0) endif C C C FILTER THE SAMPLE POINTS AND SYMULTANEOUSLY FIND THE VARIANCE C C DETERMINE THE NUMBER OF YX-ARRAY LOCATIONS TO BE USED nyx=npole+1 C C LOCATION IN YX ARRAY OF VALUE FILTERED BY ALL POLES jy=3-nrc C C NUMBER OF POINTS IN YX ARRAY TO THE RIGHT OF FILTERED VALUE ky=nyx-jy C C INIT YX-ARRAY (WHICH KEEPS TRACK OF POLE INPUT AND OUTPUT) do i=1,nyx yx(i,1)=0.d0 yx(i,2)=0.d0 enddo C C FIND THE AMOUNT OF SHIFT CAUSED BY THE FILTER C (9.711=360/37.071 WAS DERIVED EMPIRICALLY) msh=jidnnt(pi*dflotj(npole)/(wc*9.711d0))+ky C C FILTER AND SUM THE VARIANCES do 701 m=1,nsamp+msh C C SLIDE VALUES IN THE YX-ARRAY ONE POINT TO THE LEFT do i=2,nyx yx(i-1,1)=yx(i,1) yx(i-1,2)=yx(i,2) enddo C C COPY A NEW VALUE INTO RIGHT END OF YX ARRAY if(m.le.nsamp) then yx(nyx,1)=yio(m)-(slope*dflotj(m)+rcept) C C ADD VALUE TO UNFILTERED VARIANCE SUM varu=varu+(yx(nyx,1))**2.d0 else yx(nyx,1)=0.d0 endif C C APPLY BUTTERWORTH POLE PAIRS j=nyx+2 do k=1,nppair j=j-2 yx(j,2)=yx(j,1) yx(j,1)= & ax(k)*( yx(j,2) +2.d0*yx(j-1,2) +yx(j-2,2) ) & -ay1(k)*yx(j-1,1) -ay2(k)*yx(j-2,1) enddo C C APPLY RC POLE, IF ANY if(nrc.eq.1) then j=j-2 yx(j,2)=yx(j,1) yx(j,1)=ax(k)*(yx(j,2)+yx(j-1,2))-ay1(k)*yx(j-1,1) endif C C ADD VALUE TO FILTERED VARIANCE SUM if(m.gt.msh) varf=varf+(yx(jy,1))**2.d0 C C PUT FILTERED ARRAY BACK INTO INPUT ARRAY C (THIS LINE IS FOR TESTING PURPOSES ONLY; THE SUBROUTINE C SHOULD NOT CHANGE ANY VALUES IN ARRAY YIO) c if(m.gt.msh) yio(m-msh)=(yx(jy,1)) c & +(slope*dflotj(m-msh)+rcept) 701 continue C C DIVIDE OUT VARIANCE SUMS C 102 varf=varf/dflotj(nsamp) varu=varu/dflotj(nsamp) C C EXIT PROCEDURE C 999 return end C C________________________________________________________________ C C SUBROUTINE W I N D O W C________________________________________________________________ C C SUBROUTINE WINDOW APPLIES OR REMOVES A SELECTED WINDOW FROM AN C EVENLY-SPACED DOUBLE-PRECISION FUNCTION ARRAY. WINDOWS MAY BE C SPECIFIED AS DFT-EVEN OR SYMMETRIC. ALSO, A NUMBER OF C WINDOWING FUNCTIONS ARE AVAILABLE AND MAY BE SELECTED WITH C KTYPE. C C KTYPE - INTEGER*4. NUMBER INDICATING THE TYPE OF WINDOW. C FOLLOWING IS A LIST OF WINDOWS. A STAR (*) C INDICATES THAT THE OPTION HAS BEEN INSTALLED: C C KTYPE NAME (OTHER NAMES OR DESC) ALPHA C * 0 = RECTANGLE (BOXCAR, NO WINDOWING) - C * 1 = TRIANGLE (FEJER OR BARTLET) - C * 2 = COS**ALPHA (HANNING) 1.0 - 4.0 C * 3 = HAMMING (MODIFIED HANNING) - C * 4 = RIESZ (SIMPLEST CONT POLYNOM) 2.0 C 5 = RIEMANN - C 6 = DE LA VALLE-POUSSIN (CUBIC) - C * 7 = TUKEY (COSINE TAPERED) 0.0 - 1.0 C 8 = BOHMAN - C 9 = POISSON (2-SIDED EXPONENTIAL) 2.0 - 4.0 C 10 = HANNING-POISSON 0.5 - 2.0 C 11 = CAUCHY 3.0 - 5.0 C 12 = GAUSSIAN (MIN TIME/BANDWID PROD) 2.5 - 3.5 C 13 = DOLPH-TCHEBYSHEV 2.5 - 4.0 C * 14 = KAISER-BESSEL (OVERALL BEST) 2.0 - 3.5 C 15 = BARCILON-TEMES 3.0 - 4.0 C C ALPHA - REAL*8. THE DEGREE OF WINDOWING FOR WINDOWS WHICH C ALLOW A SELECTION. IF ALPHA=0.D0, THE BEST VALUE C FOR THAT WINDOW WILL BE SELECTED AUTOMATICALLY C (ALPHA WILL NOT BE CHANGED HOWEVER). C WFRAC - REAL*8. WINDOW FRACTION FOR SPLITTING WINDOWS. IF C WFRAC IS GREATER THAN ZERO AND LESS THAN ONE, A C SPLIT WINDOW WILL BE INDICATED. A SPLIT WINDOW IS C SIMPLY THE SPECIFIED WINDOW, KTYPE, APPLIED TO ONLY C THE ENDS OF THE PROFILE SUCH THAT THE PROFILE TAPERS C TO ZERO AT EACH END AND HAS ITS VALUES UNCHANGED IN C THE MIDDLE. WFRAC SPECIFIES THE FRACTION OF THE C ORIGINAL PROFILE LENGTH TAKEN BY WINDOW TAPERING. C NOTE THAT A BOX CAR WINDOW (KTYPE=0) CANNOT BE C SHORTENED USING WFRAC; THAT PROCEDURE MUST BE DONE C WITH ZERO EXTENSION (SEE SUBROUTINE EXTEND OR USE C KTYPE=4, RIESZ WINDOW AND ALPHA NEARLY 0). C KDFT - INTEGER*4. DETERMINES WHETHER THE WINDOW IS C DFT-EVEN (KDFT => 1) OR SYMMETRIC (KDFT <= 0). C DFT-EVEN WINDOWING SHOULD BE USED WITH DESCRETE C FOURIER TRANSFORMS. HOWEVER, IF THE DATA ARE TO BE C PADDED WITH ZEROS AFTER WINDOWING, SYMMETRIC C WINDOWING IS MORE APPROPRIATE. C NVERSE - INTEGER*4. DETERMINES WHETHER TO APPLY OR REMOVE C THE WINDOW: C NVERSE <= 0 MEANS TO APPLY THE WINDOW. C NVERSE => 1 MEANS TO REMOVE THE WINDOW. C NOTE: SOME WINDOWS GO TO ZERO AT THE END POINTS. IN C THESE CASES, WHEN THE WINDOW IS REMOVED, EACH END C POINT OF THE FUNCTION WILL BE SET EQUAL TO THE C ADJACENT FUNCTION VALUE. C NPTS - INTEGER*4. THE NUMBER OF DATA POINTS IN ARRAY YFNC. C C INPUT/OUTPUT ARGUMENT: C YFNC - REAL*8. ARRAY CONTAINING THE DATA TO BE WINDOWED. C C C SUBROUTINE WINDOW WRITTEN BY ROB BRACKEN, USGS, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 26AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.1, 19970325. C HP-9000 SERIES 700/800 UNIX VERSION 1.2, 19970602. C HP-9000 SERIES 700/800 UNIX VERSION 2.0, 19970602. C HP-9000 SERIES 700/800 UNIX VERSION 2.1, 19970609. C C subroutine window(ktype,alpha,wfrac,kdft,nverse,npts,yfnc) SAVE C C DECLARATIONS C C INPUT ARGUMENTS integer*4 ktype real*8 alpha,wfrac integer*4 kdft,nverse,npts C C INPUT/OUTPUT ARGUMENT real*8 yfnc(*) C C INTERNAL VARIABLES real*8 alpha2,alpha3,wfrac2 integer*4 ktype2,kdft2,npts2 data ktype2,alpha2,wfrac2,kdft2,npts2 & / -999, 1.7d+37, 1.7d+37, -999, -999 / C C INTERNAL WINDOW integer*4 mxw parameter(mxw= 800 000 ) real*8 dwindo(mxw) C C C CHECK WINDOW TYPE C if(ktype.ge.1.and.ktype.le.15) goto 254 if(ktype.eq.0) goto 990 call errfor(ktype,'(WINDOW) Non-existant window type') C C C USE THE PREVIOUSLY CALCULATED WINDOW C 254 if(ktype.eq.ktype2.and.alpha.eq.alpha2.and.wfrac.eq.wfrac2 & .and.kdft.eq.kdft2.and.npts.eq.npts2) then C C DWINDO CONTAINS THE DESIRED WINDOW if(nverse.le.0) then C C APPLY THE PREVIOUSLY CALCULATED WINDOW do 751 i=1,npts2 yfnc(i)=yfnc(i)*dwindo(i) 751 continue else C C REMOVE THE PREVIOUSLY CALCULATED WINDOW do 752 i=2,npts2-1 yfnc(i)=yfnc(i)/dwindo(i) 752 continue C C FIRST WINDOW POINT MAY BE ZERO if(dwindo(1).eq.0.d0) then yfnc(1)=yfnc(2) else yfnc(1)=yfnc(1)/dwindo(1) endif C C LAST WINDOW POINT MAY BE ZERO if(dwindo(npts2).eq.0.d0) then yfnc(npts2)=yfnc(npts2-1) else yfnc(npts2)=yfnc(npts2)/dwindo(npts2) endif endif goto 990 endif C C C CALCULATE A NEW WINDOW C C CHECK SIZE OF WINDOW AGAINST INTERNAL ARRAY SIZE if(npts.gt.mxw) call errfor(mxw, & '(WINDOW) Internal array dwindo too small for data') C C UPDATE INTERNAL WINDOW CHARACTERISTICS ktype2=ktype alpha2=alpha wfrac2=wfrac kdft2=kdft npts2=npts C C CHECK VALUE OF WFRAC TO SET UP FOR POSSIBLE SPLIT WINDOW if(wfrac2.gt.0.d0.and.wfrac2.lt.1.d0) then C C WINDOW TO BE SPLIT WITH ONES FILLING THE MIDDLE ksym=0 if(kdft2.le.0) ksym=1 npts3=jidnnt(wfrac2*dflotj(npts2-ksym))+ksym ncopy=(npts3-ksym)/2+ksym nfill=npts2-npts3 else C C NORMAL WINDOW npts3=npts2 ncopy=0 nfill=0 endif C C INITIALIZE INTERNAL ARRAY, DWINDO do 753 i=1,npts3 dwindo(i)=1.d0 753 continue C C SET INTERNAL INVERSE TO APPLY THE WINDOW nversx=0 C C SELECT THE DESIRED WINDOW goto (101,102,103,104,105,106,107,108,109, & 110,111,112,113,114,115),ktype2 goto 990 C C C ================= BEGIN WINDOWS ================= C C TRIANGLE (FEJER OR BARTLET) ALPHA = - 101 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=1.d0 call fejer(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C COS**ALPHA (HANNING) ALPHA = 1.0 - 4.0 102 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=2.d0 call hann3(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C HAMMING (MODIFIED HANNING) ALPHA = - 103 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=0.d0 call hamm3(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C RIESZ (SIMPLEST CONT POLYNOM) ALPHA = 2.0 104 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=2.d0 call riesz(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C RIEMANN ALPHA = - 105 continue call errfor(ktype2, & '(WINDOW) riemann window not installed') goto 150 C C DE LA VALLE-POUSSIN (CUBIC) ALPHA = - 106 continue call errfor(ktype2, & '(WINDOW) de la valle-poussin window not installed') goto 150 C C TUKEY (COSINE TAPERED) ALPHA = 0.0 - 1.0 107 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=0.5d0 call tukey(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C BOHMAN ALPHA = - 108 continue call errfor(ktype2, & '(WINDOW) bohman window not installed') goto 150 C C POISSON (2-SIDED EXPONENTIAL) ALPHA = 2.0 - 4.0 109 continue call errfor(ktype2, & '(WINDOW) poisson window not installed') goto 150 C C HANNING-POISSON ALPHA = 0.5 - 2.0 110 continue call errfor(ktype2, & '(WINDOW) hanning-poisson window not installed') goto 150 C C CAUCHY ALPHA = 3.0 - 5.0 111 continue call errfor(ktype2, & '(WINDOW) cauchy window not installed') goto 150 C C GAUSSIAN (MIN TIME/BANDWID PROD) ALPHA = 2.5 - 3.5 112 continue call errfor(ktype2, & '(WINDOW) gaussian window not installed') goto 150 C C DOLPH-TCHEBYSHEV ALPHA = 2.5 - 4.0 113 continue call errfor(ktype2, & '(WINDOW) dolph-tchebyshev window not installed') goto 150 C C KAISER-BESSEL (OVERALL BEST) ALPHA = 2.0 - 3.5 114 continue alpha3=alpha2 if(alpha3.eq.0.d0) alpha3=3.d0 call kais3(alpha3,nversx,kdft2,npts3,dwindo) goto 150 C C BARCILON-TEMES ALPHA = 3.0 - 4.0 115 continue call errfor(ktype2, & '(WINDOW) barcilon-temes window not installed') goto 150 C 150 continue C C ================= END WINDOWS ================= C C C IF SPLIT WINDOW, SLIDE RIGHT HALF RIGHT & FILL SPLIT W/ONES if(nfill.gt.0) then do i=npts2,npts2-ncopy+1,-1 dwindo(i)=dwindo(i-nfill) enddo do j=i,i-nfill+1,-1 dwindo(j)=1.d0 enddo endif C C RETURN TO APPLY OR REMOVE THE WINDOW FROM THE DATA goto 254 C C EXIT PROCEDURE C 990 continue end C C________________________________________________________________ C C SUBROUTINE W I N N O W C________________________________________________________________ C C SUBROUTINE WINNOW WINNOWS (DECIMATES) AND "UNWINNOWS" A PROFILE C OF DATA. C C WINNOWING IS SIMPLY THE PROCESS OF REMOVING INTERVENING DATA C POINTS AND COMPACTING THE DATA PROFILE. AFTER WINNOWING, THE C NEW PROFILE WILL CONTAIN ONLY EVERY NTH DATA POINT OF THE C ORIGINAL. THE EFFECT OF WINNOWING IS TO LOWER THE NYQUIST C FREQUENCY OF THE DATA AND THEREFORE IMPOSE ALIASING ON THE NEW C NYQUIST INTERVAL. CONSEQUENTLY, BEFORE WINNOWING, A LOW-PASS C FILTER IS APPLIED TO THE ORIGINAL DATA PROFILE TO REMOVE OR C MINIMIZE ANY FREQUENCIES THAT ARE ABOVE THE NEW NYQUIST. C C TECHNICALLY, "UNWINNOWING" SHOULD PERFORM THE EXACT INVERSE OF C WINNOWING, THAT IS, RESTORE THE ORIGINAL PROFILE. HOWEVER, THE C WINNOWING PROCESS ESSENTIALLY ERASES THE HIGHER FREQUENCIES C MAKING THEM UNAVAILABLE FOR RESTORING. THEREFORE, THE C "UNWINNOWING" AS PERFORMED BY THIS SUBROUTINE RESTORES THE C ORIGINAL SAMPLING INTERVAL AND NUMBER OF POINTS BUT WITH THE C SPECTRUM OF ONLY THE WINNOWED PROFILE. C C C THIS SUBROUTINE USES A LOW-PASS TIME-DOMAIN FILTER WITH A C SQUARED BUTTERWORTH CHARACTERISTIC AND CUTOFF FREQUENCY 0.8053 C OF THE WINNOWED NYQUIST FREQUENCY (WC=0.8053*PI/KWIN). THE C CHARACTERISTICS OF THIS FILTER INCLUDE ZERO PHASE SHIFT, AND A C TRANSIENT BAND WHOSE ROLLOFF CAN BE SPECIFIED WITH ARGUMENT C NPOLE. IF NPOLE=16, THE FOLLOWING ROLLOFF WILL OCCUR: C C FREQUENCY POWER(DB) POWER(RATIO) COMMENTS C 0.4686*PI/KWIN 0 DB = 1-(1/2**24) REAL*4 PREC LIMIT C 0.5000*PI/KWIN -0 DB = 1-(1/2**21) HALF NYQUIST C 0.5110*PI/KWIN -0 DB = 1-(1/2**20) C 0.7834*PI/KWIN -3 DB = 1/2 SQUARED HALF POWER C 0.8053*PI/KWIN -6 DB = 1/2**2 UNSQUARED CUTOFF FQ C 1.0000*PI/KWIN -60 DB = 1/2**20 NYQUIST C 1.0444*PI/KWIN -72 DB - 1/2**24 REAL*4 PREC LIMIT C C THE SQUARED BUTTERWORTH WITH NPOLE=16 (EFFECTIVE NUMBER OF C POLES IS 32) HAS A ROLLOFF OF 12 DB PER POLE PER OCTIVE. C HOWEVER, WITH THIS LARGE NUMBER OF POLES, RINGING CAN BECOME A C SIGNIFICANT ISSUE. SEE DISCUSSION UNDER ARGUMENT NPOLE. C C BEFORE FILTERING, THE ORIGINAL DATA PROFILE WILL BE EXTENDED C WITH A BURG PREDICTION ALGORITHM. THIS PROCESS DIMINISHES THE C END EFFECTS (RINGING) LIKELY TO OCCUR WITH THE TIME-DOMAIN C FILTER. HOWEVER, IT SHOULD BE NOTED THAT IF STEPS OCCUR WITHIN C THE PROFILE, RESULTANT RINGING MAY BE VISIBLE (AS RINGING) C WITHIN THE WINNOWED PROFILE. C C C INPUT ARGUMENTS: C KWIN - INTEGER*4. THE WINNOW FACTOR, AN INTEGER FACTOR BY C WHICH THE ORIGINAL SAMPLING INTERVAL IS TO BE C MULTIPLIED VIA REMOVAL OF INTERVENING DATA POINTS. C THE FIRST POINT OF THE ORIGINAL PROFILE WILL ALWAYS C BE KEPT, AS WELL AS EVERY KWINTH POINT AFTERWARDS. C FOR EXAMPLE, IF KWIN=3, POINTS 1,2,3, ... N OF THE C WINNOWED PROFILE WILL HAVE COME FROM POINTS 1,4,7, C ... 3N+1 OF THE ORIGINAL PROFILE. THE TOTAL NUMBER C OF POINTS IN THE WINNOWED PROFILE WILL EQUAL C 1+(NPTS-1)/KWIN (INTEGER DIVISION ASSUMED). C C IF KWIN IS GREATER THAN 1, WINNOWING WILL BE DONE. C IF KWIN IS LESS THAN -1, UNWINNOWING WILL BE DONE. C IF KWIN IS -1, 0, OR 1, NOTHING WILL BE DONE. C (IF YOU SIMPLY NEED A PROFILE LOW-PASS FILTERED, SEE C SUBROUTINE BWTDF) C C NOTE: TO MINIMIZE WINNOWING COMPUTATION TIME FOR C LARGE PROFILES, KWIN SHOULD BE FACTORABLE WITH C SEVERAL SMALL FACTORS SUCH AS 2, 3, AND 5 AND LSTAGE C SHOULD BE SET TO 1. CONVERSELY, IF KWIN IS ALLOWED C TO BE A LARGE PRIME NUMBER, THE TIME REQUIRED FOR C THE EXTENSION PART OF THE FILTERING PROCESS MAY C BECOME UNREASONABLY LARGE AS WILL THE NECESSARY C EXTENSION LENGTH. C C LSTAGE - INTEGER*4. DISABLES OR ENABLES STAGING DURING C WINNOWING. IF LSTAGE IS 0, STAGING WILL NOT BE C USED; IF LSTAGE IS 1, STAGING WILL BE USED. LSTAGE C IS NOT USED DURING UNWINNOWING. WINNOWING IN STAGES C INVOLVES COMPOUNDED WINNOWINGS BY EACH OF THE C FACTORS OF KWIN. BY WINNOWING IN STAGES, TIME IS C SAVED BUT ACCURACY IS COMPROMISED BECAUSE THE C REPEATED FILTERING CAUSED BY STAGING OVERPRINTS THE C UPPER HALF OF THE NYQUIST INTERVAL WITH VARIOUS C PARTS OF THE FILTER CURVE. THE GREATEST ACCURACY IS C OBTAINED BY CHOOSING NO STAGING OR BY CHOOSING C WINNOW FACTORS THAT ARE PRIME NUMBERS. C C NPOLE - INTEGER*4. THE NUMBER OF POLES DESIRED FOR THE C BUTTERWORTH FILTER DURING WINNOWING. IF NPOLE=0, C THE FILTERING STEP WILL BE SKIPPED. NPOLE IS NOT C USED DURING UNWINNOWING. THE FILTER HAS BEEN C PRE-DEFINED AS A SQUARED TIME-DOMAIN BUTTERWORTH (NO C PHASE SHIFT) WITH A CUT-OFF FREQUENCY OF 0.8053 * C (THE NEW NYQUIST FREQ). BECAUSE THE FILTER IS C SQUARED, THE EFFECTIVE NUMBER OF POLES IS TWICE C NPOLE. THE FOLLOWING TABLE GIVES A -72 dB STOP C FREQUENCY IN OCTAVES RELATIVE TO THE CUTOFF C FREQ, DERIVED FROM NPOLE=LOG(64)/LOG(FSTOP/FCUT): C C NPOLE STOP FQ SUPPRESSION AT NYQUIST C (OCTAVES) (dB) (AMPLITUDE) C 0 FILTERING IS SKIPPED FOR NPOLE=0 C 1 6 -3.25 2/3 C 2 3 -7.5 1/2 C 4 1.5 -15 1/5 C 8 0.75 -30 1/31 C 16 0.375 -60 1/1000 C 32 0.1875 -120 1/1000000 C C THE PRE-DEFINED ASPECTS OF THE FILTER WERE C ORIGINALLY DESIGNED TO PRODUCE A VERY SHARP CUTOFF C USING NPOLE=16. HOWEVER, THIS PRODUCES SIGNIFICANT C RINGING IF THE FUNCTION CONTAINS ANY "SHARP C CORNERS". BY REDUCING TO NPOLE=4 RINGING DECREASES C SUBSTANTIALLY, BUT THE SUPPRESSION OF FREQUENCIES C ABOVE NYQUIST IS COMPROMISED. NEVERTHELESS, IN C THE NORMALLY-PINK POTENTIAL-FIELD FUNCTIONS, NPOLE=4 C IS USUALLY SUFFICIENT. C C NPTS - INTEGER*4. THE ORIGINAL, UNWINNOWED NUMBER OF C POINTS IN ARRAY YFNC. NPTS SHOULD ALWAYS BE THE C UNWINNOWED NUMBER OF POINTS REGARDLESS OF WHETHER C THE PROFILE IS BEING WINNOWED (POSITIVE KWIN) OR C UNWINNOWED (NEGATIVE KWIN). IF NPTS IS NEGATIVE C DURING WINNOWING, ONLY THE LOW-PASS FILTERING C STEP(S) OF THE WINNOWING PROCESS WILL BE PERFORMED; C THE ACTUAL WINNOWING WILL BE SKIPPED (IF UNDER THESE C CIRCUMSTANCES NDIM=0, THE NEW DIMENSION RETURNED C WILL BE GREATER THAN ITS VALUE HAD NPTS BEEN C POSITIVE). SIMILARLY, IF NPTS IS NEGATIVE DURING C UNWINNOWING, THE UNWINNOWING WILL BE SKIPPED. C (HOWEVER, BECAUSE UNFILTERING IS NOT PERFORMED C DURING UNWINNOWING, NEGATIVE NPTS DURING UNWINNOWING C ACCOMPLISHES NOTHING; THE OPTION IS INCLUDED MERELY C FOR SYMMETRY). C C INPUT/OUTPUT ARGUMENTS: C YFNC - REAL*8. ARRAY OF DIMENSION NDIM CONTAINING NPTS OF C LEFT-JUSTIFIED DATA TO BE WINNOWED OR C 1+(NPTS-1)/ABS(KWIN) POINTS TO BE UNWINNOWED. C NDIM MUST BE LARGER THAN NPTS TO ACCOMOMDATE C SUBROUTINE COMPUTATIONS. THE MINIMUM VALUE OF NDIM C MAY BE FOUND AS DESCRIBED IN NDIM. C C NOTE: THE NUMBER OF WINNOWED POINTS IS NEVER TAKEN C AS AN INPUT ARGUMENT. SEE DESCRIPTION OF NPTW. C C NDIM - INTEGER*4. THE DIMENSION OF YFNC. DURING WINNOWING C AND UNWINNOWING NDIM MUST BE GREATER THAN NPTS TO C ACCOMMODATE PROFILE EXTENSIONS NEEDED FOR FILTERING C AND TRANSFORMATIONS. IF THE SPACE REQUIREMENTS IN C YFNC ARE GREATER THAN NDIM, AN ERROR WILL BE C GENERATED. HOWEVER, IF NDIM IS GIVEN A VALUE OF C ZERO DURING THE CALL, THE MINIMUM NEEDED DIMENSION C OF YFNC WILL BE CALCULATED AND RETURNED IN NDIM C WITHOUT ALTERING YFNC. THIS WILL BE DONE REGARDLESS C OF THE SIGN OF NPTS. FOLLOWING IS AN EXAMPLE OF THE C USAGE OF NDIM: C C integer*4 MXYFNC C parameter(MXYFNC = 200 000 ) C real*8 yfnc(MXYFNC) C integer*4 ndim C . . . C . . . C . . . C C FIND A MINIMUM VALUE FOR NDIM C ndim=0 C call winnow(kwin,lstage,npole,npts, yfnc,ndim, nptw) C if(ndim.gt.MXYFNC) then C write(6,888) ndim C 888 format('Size requirements of array yfnc are ', C & ,'greater than its allocation.',/,'Please modify' C & ,' parameter MXYFNC to be greater than ndim=',i7) C stop C endif C C C C THE VALUE OF NDIM IS NOW ASSURED TO BE LARGE ENOUGH C call winnow(kwin,lstage,npole,npts, yfnc,ndim, nptw) C C OUTPUT ARGUMENT: C NPTW - INTEGER*4. THE NUMBER OF POINTS IN THE WINNOWED C PROFILE. THE NUMBER OF WINNOWED POINTS IS C 1+(ABS(NPTS)-1)/ABS(KWIN). THIS VALUE CAN ALWAYS BE C CALCULATED EXACTLY FROM KWIN AND NPTS; BUT, THE C VALUE OF NPTS CANNOT BE CALCULATED EXACTLY FROM KWIN C AND THE WINNOWED NUMBER OF POINTS. THEREFORE, THE C NUMBER OF WINNOWED POINTS IS NEVER TAKEN AS AN INPUT C ARGUMENT. RATHER, IT IS PROVIDED HERE AS AN OUTPUT C ARGUMENT STRICTLY FOR INFORMATIONAL PURPOSES. A C VALUE FOR NPTW IS PROVIDED WHENEVER A CALL IS MADE C TO THIS SUBROUTINE REGARDLESS OF THE VALUE OF NDIM C OR KWIN (ALTHOUGH THE VALUE OF KWIN INFLUENCES C NPTW). C C C SUBROUTINE WINNOW WRITTEN BY ROB BRACKEN, USGS. C FORTRAN 77, HP FORTRAN/9000, HP-UX RELEASE 11.0 C VERSION 1.0, 19970506, (INITIAL ROUTINE). C VERSION 2.0, 20000201, (ADDED ARGUMENTS, LSTAGE AND NPOLE). C C subroutine winnow(kwin,lstage,npole,npts, yfnc,ndim, nptw) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 kwin,lstage,npole,npts C C INPUT/OUTPUT ARGUMENTS real*8 yfnc(*) integer*4 ndim C C OUTPUT ARGUMENT integer*4 nptw C C WINNOWING STAGES integer*4 mxfac parameter(mxfac=60) real*8 winfac(mxfac) integer*4 nstage C C BUTTERWORTH-FILTER PARAMETERS C C KSQ=SQUARED FILTER (TO MAKE ZERO PHASE), 0=OFF, 1=ON C RINGDB=FILTER END-EFFECT RINGING IN DECIBELS C WCRNYQ=CUTOFF FREQUENCY RELATIVE TO NYQUIST C WSRWC=STOP BAND (-36 DB UNSQRD, -72 DB SQRD) FQ REL TO FCUT C (STOP BAND FREQ AT 1.29*FCUT PRODUCES A 16-POLE FILTER) integer*4 ksq real*8 ringdb,wcrnyq,wsrwc parameter(ringdb=-72.d0,ksq=1) c parameter(wcrnyq=0.8053d0,wsrwc=1.2968d0) parameter(wcrnyq=0.8053d0) real*8 fcut,fstop C C UNWINNOWING real*8 fdavg,fdnyq C C C FIND THE WINNOWED NUMBER OF POINTS, NPTW C if(kwin.ne.0) then nptw=1+(jiabs(npts)-1)/jiabs(kwin) else nptw=jiabs(npts) endif C C DETERMINE WHETHER TO WINNOW OR UNWINNOW C if(kwin.ge.2) goto 101 if(kwin.le.-2) goto 102 if(ndim.eq.0) ndim=jiabs(npts) goto 990 C C WINNOW C C DETERMINE THE FACTORS OF WINNOWING STAGES 101 if(lstage.ge.1) then C C STAGING IS DESIRED call factr(dflotj(kwin),winfac,nstage) if(nstage.gt.mxfac) call errfor(nstage, & '(winnow) nstage exceeds mxfac') else C C STAGING IS NOT DESIRED nstage=3 winfac(3)=dflotj(kwin) endif C C DO WINNOWING IN STAGES FROM SMALLEST FACTORS TO LARGEST ndim0=ndim npts2=jiabs(npts) do j=3,nstage C C GET THE NEXT LARGEST FACTOR IN KWIN kfac=jidnnt(winfac(j)) C C CHECK WHETHER TO FILTER if(npole.ge.1) then C C FILTER THIS STAGE if(npts.gt.0.or.j.eq.3) then fcut=0.5d0*wcrnyq/dflotj(kfac) else fcut=fcut/dflotj(kfac) endif wsrwc=dexp(dlog(64.d0)/dble(npole)) fstop=wsrwc*fcut ndim2=ndim0 call bwtdf(fcut,fstop,ksq,ringdb,npts2, yfnc,ndim2) else C C DO NOT FILTER THIS STAGE ndim2=ndim0 if(ndim2.le.0) then ndim2=npts2 endif C if(ndim0.le.0) then C C ADJUST THE DIMENSION if(ndim2.gt.ndim) ndim=ndim2 if(npts.gt.0) npts2=1+(npts2-1)/kfac else if(npts.gt.0) then C C WINNOW THE PROFILE C C IF ONLY 2 POINTS LEFT IN PROFILE, TAKE THE AVERAGE if(npts2.eq.2.and.kfac.gt.1) & yfnc(1)=(yfnc(1)+yfnc(2))/2.d0 C i=1 do k=1+kfac,npts2,kfac i=i+1 yfnc(i)=yfnc(k) enddo npts2=i endif enddo C C CHECK THE NUMBER OF POINTS IN THE WINNOWED PROFILE if(ndim0.gt.0.and.npts.gt.0 & .and.npts2.ne.1+(jiabs(npts)-1)/kwin) & call errfor(0,'(winnow) Algorithm error') goto 990 C C UNWINNOW C C CHECK THE VALUE OF NDIM 102 ndim2=0 call unwin0(jiabs(kwin),jiabs(npts),yfnc(1),ndim2, & yfnc(2),fdavg,fdnyq) ndim0=jmax0(jiabs(npts),ndim2) C IF KRIESZ=1 IN UNWIN0, THIS STATEMENT MUST BE MODIFIED ndim2=((ndim2/jiabs(kwin)-1)/2)*2 if(ndim.le.0) then ndim=ndim0+ndim2 goto 990 else if(ndim.lt.ndim0+ndim2) call errfor(ndim0+ndim2, & '(winnow) ndim is too small') endif C UNWINNOW THE PROFILE if(npts.gt.0) & call unwin0(jiabs(kwin),jiabs(npts),yfnc(1),ndim, & yfnc(ndim0+1),fdavg,fdnyq) C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C FUNCTION Z M B E S 5 C________________________________________________________________ C C FUNCTION ZMBES5 IS A ZERO-ORDER MODIFIED BESSEL FUNCTION OF THE C FIRST KIND. ZMBES5 HAS BEEN WRITTEN SO AS TO MAINTAIN MAXIMUM C POSSIBLE DOUBLE PRECISION ACCURACY WHILE EXTENDING THE DOMAIN. C THE FUNCTION IS AS FOLLOWS: C C inf C Io[Z] = sum[ ( ((Z/2)**k)/k! )**2 ] C k=0 C C FUNCTION: C ZMBES5 - REAL*8. THE VALUE OF THE MODIFIED BESSEL FUNCTION. C C INPUT ARGUMENT: C Z - REAL*8. THE INDEPENDANT VARIABLE. Z MUST BE C A SMALL POSITIVE NUMBER WHERE 0 <= Z <= 713.8 C ACCURACY SHOULD BE MAINTAINED TO WITHIN 1.D-16 FOR C Z < 130. ACCURACY BEGINS DECREASING ABOVE Z = 130. C THE LARGEST VALUE OF Z FOR WHICH THE BESSEL FUNCTION C CAN BE CONVERTED TO A REAL*4 VALUE IS 94.2 C C FUNCTION ZMBES5 WRITTEN BY ROB BRACKEN, USGS, 19AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 19AUG96. C C function zmbes5(z) C C DECLARATIONS C C FUNCTION AND ARGUMENT real*8 zmbes5,z C C OTHER VARIABLES real*8 x2,dlogx2,fact integer*4 mxterm parameter(mxterm=1000) real*8 term0,term2,term(0:mxterm) real*8 peak,dk real*8 zmbesl,zmbesr C C CHECK FOR Z VALUE WITHIN BOUNDS C if(z.lt.0.d0) then call errfor(0,'bessel fcn z is negative') else if(z.eq.0.d0) then zmbes5=1.d0 goto 990 else if(z.gt.713.8d0) then call errfor(0,'bessel fcn z is larger than 713.8') endif C C INITIALIZE VALUES C C NUMERATOR x2=z/2.d0 dlogx2=dlog10(x2) C INITIAL DENOMINATOR VALUE (FACTORIAL VALUE) fact=1.d0 C ASSIGN VALUE OF THE 0 TERM term0=1.d0 term2=term0*term0 term(0)=term2 C INITIALIZE PEAK TERM VALUE AND LOCATION peak=1.d0 kpeak=0 C C CALCULATE THE FUNCTION TERMS (1-NTERM) AND FIND PEAK TERM VALUE C C INDEXING k=0 201 k=k+1 if(k.gt.mxterm) call errfor(k, & '(zmbes5) Index exceeded upper boundary of array term()') dk=dflotj(k) C if(k.le.170.and.dlogx2*dk.lt.308.d0) then C C IF FACTORIAL AND NUMERATOR**K OKAY, THEN RECALC NEW TERM fact=fact*dk term0=((x2**dk)/fact) else C C FACTORIAL OR NUMERATOR**K TOO LARGE; ADJUST FM PREV TERM term0=term0*(x2/dk) endif C C SQUARE TERM AND STORE IN ARRAY FOR FUTURE REFERENCE term2=term0*term0 term(k)=term2 C C FIND LOCATION OF PEAK TERM if(term0.gt.peak) then peak=term0 kpeak=k endif C C IF CURRENT TERM CAN AFFECT PEAK TERM, CALCULATE NEXT TERM if(peak-term0.ne.peak) goto 201 nterm=k C C SUM THE FUNCTION TERMS FROM THE LEFT AND RIGHT C C ADD THE TERMS FROM THE LEFT TO THE PEAK zmbesl=0.d0 do 702 k=0,kpeak zmbesl=zmbesl+term(k) 702 continue C C ADD THE TERMS FROM THE RIGHT TO THE PEAK zmbesr=0.d0 do 703 k=nterm,kpeak+1,-1 zmbesr=zmbesr+term(k) 703 continue C C ADD THE LEFT AND RIGHT SUMS zmbes5=zmbesl+zmbesr C C EXIT PROCEDURE C 990 return end