C C This source-code file contains (in alphabetical order by name) C the source code of the following 18 subroutines: C C burgdp.f C dfttdc.f C dftfdr.f C errfor.f C extend.f C fejer.f C * fftdc2.f C forkdc.f C hamm3.f C hann3.f C kais3.f C lsqply.f C riesz.f C trendr.f C tukey.f C varfqd.f C window.f C zmbes5.f 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 Two of them were compiled with the following flags. It is not C known whether these flags are actually required for successful C compilation or execution. The +E7 flags probably are NOT C required for either of the subroutines. The -K flag probably C IS required. All 22 subroutines may be compiled with both the C +E7 and -K flags without any known problems: C C burgdp.f +E7 C varfqd.f +E7 -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 D F T T D C C________________________________________________________________ C C SUBROUTINE DFTTDC PERFORMS A DISCRETE FOURIER TRANSFORM ON ONE- C DIMENSIONAL DATA. THE TIME DOMAIN AND FREQUENCY DOMAIN ARE C BOTH PASSED IN ONE DOUBLE-COMPLEX ARRAY. HOWEVER, THE TIME AND C FREQUENCY DOMAINS ARE ALLOWED TO DIFFER IN LENGTH. THE MAXIMUM C PROFILE LENGTH IS DETERMINED INTERNALLY BY THE PARAMETER C MXWORK. IF THE ARGUMENT NIN EXCEEDS MXWORK, AN ERROR WILL BE C GENERATED. C C DIFFERING LENGTH DOMAINS ARE INTERPRETED BY DFTTDC AS REQUIRING C A TIME-DOMAIN TRUNCATION OR EXTENSION WITH ZEROS (SEE C INTERPRETATION 1, 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 FREQUENCY DOMAINS C ARE EQUAL, THIS SUBROUTINE PRODUCES A TRUE DFT (SIGNI=-1 FOR C FORWARD AND SIGNI=+1 FOR INVERSE). THE SCALE FACTORS ARE C 1/NOUT 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. THIS C PROCEDURE IS DONE BY THIS SUBROUTINE, DFTTDC. BY EXTENDING THE C TIME DOMAIN WITH ZEROS, THE NUMBER OF POINTS IN THE FREQUENCY C DOMAIN INCREASES BUT THE OVERALL SHAPE OF THE FREQUENCY-DOMAIN C PROFILE REMAINS THE SAME. SIMILARLY WITH TRUNCATION. (NOTE C THAT TRUNCATION OF THE TIME DOMAIN IS SO SIMPLE THAT IT SHOULD C BE DONE IN THE SUBROUTINE CALL. IF, IN THE FORWARD CALL, THE C TIME-DOMAIN HAS MORE POINTS THAN THE FREQUENCY DOMAIN, DFTTDC 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. THE C PROCEDURE IS DONE BY SUBROUTINE DFTFDC. BY EXTENDING THE HIGH C FREQUENCIES IN THE FREQUENCY DOMAIN WITH ZEROS, THE NUMBER OF C POINTS IN THE TIME DOMAIN (AFTER INVERSE TRANSFORM) INCREASES C BUT THE OVERALL SHAPE OF THE TIME-DOMAIN PROFILE REMAINS THE C SAME. SIMILARLY WITH TRUNCATION. (TRUNCATION OF THE FREQUENCY C DOMAIN IS NOT AS SIMPLE AS IN THE TIME DOMAIN; THEREFORE DFTFDC C DOES NOT ALLOW THE POSSIBILITY OF A DEGENERATE TIME DOMAIN. C RATHER, IF ON INVERSE TRANSFORM THE NUMBER OF TIME-DOMAIN C POINTS IS LESS THAN THE NUMBER OF FREQUENCY-DOMAIN POINTS, C DFTFDC PRODUCES THE FREQUENCY TRUNCATION INTERNALLY.) C C C C C DFTTDC PROVIDES THE OPTION OF PRODUCING A FREQUENCY DOMAIN WITH C MORE DETAIL IN THE FREQUENCIES BY IMPLICITLY EXTENDING THE TIME C DOMAIN TO THE RIGHT WITH ZEROS. THE IMPLICIT EXTENSION IS C FASTER THAN FIRST EXTENDING THE TIME DOMAIN AND THEN C TRANSFORMING A LARGER NUMBER OF POINTS. IF A FREQUENCY DOMAIN C WITH LESS DETAIL IS DESIRED, THE TIME-DOMAIN SHOULD BE C TRUNCATED SIMPLY BY REDUCING THE NUMBER OF TIME-DOMAIN POINTS C (NIN) IN THE SUBROUTINE CALL TO EQUAL A SMALLER NUMBER OF C FREQUENCY DOMAIN POINTS (NOUT). IF, IN THE FORWARD TRANSFORM, C NIN IS GREATER THAN NOUT, A DEGENERATE FREQUENCY DOMAIN WILL BE C PRODUCED. IT WILL LOOK SIMILAR TO A PROPER REDUCED-LENGTH C FREQUENCY DOMAIN BUT WILL HAVE INCORRECT PHASES. FOR BALANCE C THE INVERSE TRANSFORM PROVIDES AN INVERSE EQUIVALENT TO C TIME-DOMAIN TRUNCATION/EXTENSION; BUT IT ACTUALLY HAS NO C PARTICULAR UTILITY OTHER THAN MAKING AN EXPENSIVE ZERO C EXTENSION OR REPEATING THE TIME-DOMAIN INTERVAL. C C THE FOLLOWING TABLE GIVES THE EFFECTS DFTTDC HAS ON VARIOUS C EXAMPLE-LENGTH TIME AND FREQUENCY DOMAIN PROFILES: C C TD1 -1 FD +1 TD2 DESCRIPTION C NIN NOUT NOUT 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 REPEATED C C 100 200 100 PT FD ACCORDIONED-OUT TO 200 PTS C 100 200 100 TD2 = TD1 C 100 200 200 TD2 = TD1 + 100 PT ZERO EXTENSION C C 200 200 PT TD1 C C 100 100 200 PT FD ACCORDIONED-IN TO 100 PTS C 100 100 100 TD2 = FIRST HALF OF 200 PT TD1 C 100 100 200 TD2 = FIRST HALF OF 200 PT TD1 REPEATED C C (200) (100) (200 PT FD ACCORDIONED-IN TO DEGEN 100 PTS) C (200) (100) (100) (TD2 = DEGENERATE HALF-LENGTH TD1) C (200) (100) (200) (TD2 = DEGENERATE HALF-LENGTH TD1 REPEATED) C C 200 200 200 PT FD NORMAL DFT-STYLE C 200 200 100 TD2 = FIRST HALF OF TD1 C 200 200 200 TD2 = TD1 C C C C C INPUT ARGUMENTS: C NIN - INTEGER*4. NUMBER OF PTS IN CX BEFORE TRANSFORM. C NOUT - INTEGER*4. NUMBER OF PTS IN CX AFTER TRANSFORM. 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 CX - COMPLEX*16. ARRAY CONTAINING THE FUNCTION TO BE C TRANSFORMED. CX MUST HAVE DIMENSION THE LARGER OF C NIN OR NOUT. C C SUBROUTINE DFTTDC WRITTEN BY ROB BRACKEN, USGS, 13NOV96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 13NOV96. C C subroutine dfttdc(nin,nout,signi,cx) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 nin,nout real*8 signi C C INPUT ARGUMENTS complex*16 cx(*) C C INTERNAL WORK ARRAY integer*4 mxwork parameter(mxwork=200000) complex*16 cy(mxwork) common /work1/ cy C C INTERNAL VARIABLES AND CONSTANTS complex*16 ci,cw0,cwk,cwj,cxk real*8 pi,scale parameter(pi = 3.1415 92653 58979 32384 62643) C C DETERMINE WHETHER NIN IS WITHIN LIMITS C if(nin.gt.mxwork) & call errfor(mxwork,'(dfttdc) nin exceeds mxwork') C C COPY INPUT ARRAY INTO WORK ARRAY C do i=1,nin cy(i)=cx(i) enddo C C CALCULATE TRANSFORM C C INITIALIZE VARIABLES cc scale=1.d0/dflotj(nout) cc if(signi.gt.0.d0) scale=dflotj(nin)*scale cc ci=dcmplx(0.d0,1.d0) cc cw0=cdexp(ci*signi*2.d0*pi/dflotj(nout)) cc cwk=dcmplx(1.d0,0.d0) C ci=dcmplx(0.d0,1.d0) if(signi.le.0.d0) then C C FORWARD TRANSFORM c scale=1.d0/dflotj(nin) scale=1.d0/dflotj(nout) cw0=cdexp(ci*-2.d0*pi/dflotj(nout)) else C C INVERSE TRANSFORM c scale=dflotj(nout)/dflotj(nin) scale=1.d0 cw0=cdexp(ci*2.d0*pi/dflotj(nin)) endif cwk=dcmplx(1.d0,0.d0) C C DISCRETE FOURIER TRANSFORM do k=1,nout c cx(k)=dcmplx(0.d0,0.d0) cxk=dcmplx(0.d0,0.d0) cwj=dcmplx(1.d0,0.d0) do j=1,nin c cx(k)=cx(k)+cy(j)*cwj cxk=cxk+cy(j)*cwj cwj=cwj*cwk enddo c cx(k)=cx(k)*scale cx(k)=cxk*scale cwk=cwk*cw0 enddo C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE D F T F D R C________________________________________________________________ C C SUBROUTINE DFTFDR 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. C C DIFFERING LENGTH DOMAINS ARE INTERPRETED BY DFTFDR 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, DFTFDR. 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 DFTFDR 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, DFTFDR PRODUCES THE FREQUENCY TRUNCATION INTERNALLY.) C C C C C DFTFDR 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 DFTFDR 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 DFTFDR 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 C C INPUT ARGUMENTS: C NRX - INTEGER*4. NUMBER OF TIME-DOMAIN PTS IN ARRAY RX. C NCY - INTEGER*4. NUMBER OF FREQUENCY-DOMAIN PTS IN CY. 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 CY - COMPLEX*16. ARRAY CONTAINING COMPLEX FREQUENCY- C DOMAIN DATA. IF SIGNI IS ZERO OR NEG (FORWARD C TRANSFORM), CY WILL BE AN OUTPUT ARGUMENT. IF SIGNI C IS GREATER THAN ZERO (INVERSE TRANSFORM), CY WILL BE C AN INPUT ARGUMENT. C C SUBROUTINE DFTFDR WRITTEN BY ROB BRACKEN, USGS, 13MAY97. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 13MAY97. C C subroutine dftfdr(nrx,ncy,signi,rx,cy) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 nrx,ncy real*8 signi C C INPUT/OUTPUT ARGUMENTS real*8 rx(*) complex*16 cy(*) C C INTERNAL VARIABLES AND CONSTANTS complex*16 ci,cw0,cwk,cwj,csum real*8 rsum real*8 pi,scale real*8 powmin parameter(pi = 3.1415 92653 58979 32384 62643) C C C CALCULATE DISCRETE FOURIER TRANSFORM 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 DISCRETE KERNEL OF INCREMENTAL NORMALIZED ANGULAR FREQ, CW0 C FOR FREQUENCY DOMAIN TRUNCATION/EXTENSION ci=dcmplx(0.d0,1.d0) cw0=cdexp(-ci*2.d0*pi/dflotj(nrx)) C C DETERMINE WHETHER TO PERFORM FORWARD OR INVERSE TRANSFORM if(signi.le.0.d0) then C C FORWARD TRANSFORM (TIME DOMAIN TO FREQUENCY DOMAIN) C C USE INTRINSIC DFT SCALING FACTOR scale=1.d0/dflotj(nrx) C C FIND THE NUMBER OF FREQUENCIES TO CALCULATE nfq=jmin0(nrx,ncy)/2+1 C C FIND AVERAGE VALUE (ZERO FREQUENCY, WAVENUMBER = 0) rsum=0.d0 do j=1,nrx rsum=rsum+rx(j) enddo cy(1)=dcmplx(rsum*scale,0.d0) C C INITIALIZE MINIMUM POWER CHECKER c powmin=1.1d+38 C C FIND NON-ZERO FREQUENCIES INCLUDING NYQUIST cwk=cw0 do k=2,nfq csum=dcmplx(0.d0,0.d0) cwj=dcmplx(1.d0,0.d0) do j=1,nrx csum=csum+rx(j)*cwj cwj=cwj*cwk enddo cwk=cwk*cw0 C C MAKE ADJUSTMENTS AT AND ABOVE THE NYQUIST FREQUENCY c if(cdabs(csum).lt.powmin) powmin=cdabs(csum) if(k.eq.nfq) then if(ncy.gt.nrx) then C C IF CSUM IS NYQUIST, SPLIT IT if((k-1)*2.eq.nrx) csum=csum/2.d0 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*scale/4096.d0 c powmin=powmin*scale/2.d0 do i=k+1,ncy-k+1 cy(i)=dcmplx(powmin,0.d0) enddo else if(ncy.lt.nrx) then C C IF CSUM IS NEW NYQUIST, COMBINE WITH ITS REFLECTION if((k-1)*2.eq.ncy) csum=csum+dconjg(csum) endif endif C C SCALE & REFLECT THE CONJG OF THE NEWLY CALCULATED FREQ cy(k)=csum*scale cy(ncy-k+2)=dconjg(cy(k)) enddo else C C INVERSE TRANSFORM (FREQUENCY DOMAIN TO TIME DOMAIN) C C C INVERT THE KERNEL cw0=dconjg(cw0) C C INVERSE SCALING FACTOR IS 1 scale=1.d0 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=dcmplx(1.d0,0.d0) do j=1,nrx rsum=0.d0 C C ADD IN FREQS WITH WAVENUMBERS 1 THROUGH ALMOST NYQUIST cwk=cwj do k=2,kbn rsum=rsum+dreal(cy(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((cy(k)+dconjg(cy(k)))*cwk) else rsum=rsum+dreal(cy(k)*cwk) endif endif C C ADD IN AVERAGE rsum=rsum+dreal(cy(1)) C C SCALE THE SUM rx(j)=rsum*scale enddo endif 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 A LINEAR FUNCTION C REMOVED BEFORE PREDICTION AND ADDED BACK C AFTER PREDICTION. ONE OR MORE VALUES MUST C BE SPECIFIED IN THE ARGUMENT, TPARMS(). C C FOR KTYPE=7, TPARMS(1) SPECIFIES THE ROOT C LENGTH OF THE BURG FILTER. THE ROOT LENGTH C IS THE NUMBER OF DATA POINTS COUNTED FROM C THE FIRST PREDICTED POINT BACKWARDS INTO C THE EXISTING PROFILE. THESE POINTS ARE C USED TO PRODUCE THE PREDICTION-FILTER C COEFFICIENTS. FOR EXAMPLE, IF TPARMS(1)=3, C AND THE PROFILE TO BE EXTENDED HAS NIN=20 C POINTS, THEN POINTS 18, 19, AND 20 WILL BE C USED TO PREDICT THE NEXT POINT, 21. THEN, C POINT 22 WILL BE PREDICTED FROM POINTS 19, C 20, AND 21; AND SO FORTH. THIS PROCESS IS C REFLECTED AT THE LEFT END OF THE PROFILE, C IF EXTENSION IS OCCURRING THERE. C C BECAUSE THE ROOT LENGTH CANNOT BE LESS THAN C 2, VALUES OF TPARMS(1) THAT ARE LESS THAN 2 C ARE USED TO INDICATE SPECIFIC METHODS OF C DERIVING THE ROOT LENGTH. THESE METHODS C ARE GIVEN IN THE FOLLOWING TABLE, AND SOME C REQUIRE ADDITIONAL VALUES TO BE SPECIFIED C IN TPARMS(2) AND TPARMS(3): C C TPARMS(1) ROOT-LENGTH CALCULATION METHOD C -3 LEFT ROOT = TPARMS(2) C RIGHT ROOT = TPARMS(3) C -2 LEFT ROOT = C TPARMS(2)*(LEFT EXTENSION) C RIGHT ROOT = C TPARMS(2)*(RIGHT EXTENSION) C -1 LEFT ROOT = LEFT EXTENSION C RIGHT ROOT = RIGHT EXTENSION C 0 EACH ROOT = INPUT ARRAY LEN (NIN) C 1 EACH ROOT = C HALF THE TOTAL EXTENDED LENGTH C WHICH IS (NOUT-NIN)/2 C >1 EACH ROOT = TPARMS(1) C 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 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 DISCRETE 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 F F T D C 2 C________________________________________________________________ C C SUBROUTINE FFTDC2 PERFORMS FORWARD AND INVERSE FOURIER C TRANSFORMS (FT), WINDOWS, EXTENDS, AND REMOVES TRENDS IN A C DOUBLE PRECISION, 1-DIMENSIONAL ARRAY. C C RECOMMENDED PROCEDURE C C THE BEST SEPARATION OF SPECTRAL PEAKS, AND DISCERNMENT OF C SPECTRAL CONTENT CAN BE OBTAINED BY THE FOLLOWING PROCEDURE: C C 1) REMOVE THE AVERAGE OR LEAST SQUARES LINE C KTREND=1, NLSQ=1 OR 2 C 2) EXTEND THE FUNCTION IN BOTH DIRECTIONS 2 TO 4 TIMES IT'S C ORIGINAL LENGTH USING BURG PREDICTION C KEXTEN=7, XPARMS=0, JUSTX=0 C NFQD = 2*NTMD TO 4*NTMD (NFQD = 2 ** INTEGER) C 3) WINDOW THE FUNCTION WITH THE KAISER-BESSEL WINDOW C KWINDO=14, ALPHA=3 OR 3.5 C 4) PERFORM AN FFT C KDXFRM=0 C C C DESCRIPTION C C THE PURPOSE OF THIS SUBROUTINE IS TO PERFORM A FOURIER C TRANSFORM AND PROVIDE ALL REASONABLE DATA-PREPARATION STEPS IN C A SINGLE PACKAGE. ALTHOUGH THE FOURIER TRANSFORM IS A GENERAL C ABSTRACT MATHEMATICAL OPERATION, THIS SUBROUTINE PLACES ON IT A C FEW PRACTICAL CONSTRAINTS: 1) -IW IS DEFINED AS FORWARD C TRANSFORM FROM THE TIME-DOMAIN TO THE FREQUENCY-DOMAIN (+IW IS C INVERSE); 2) THE TIME DOMAIN IS CONSIDERED TO BE REAL AND THE C FREQUENCY DOMAIN COMPLEX; 3) ALL DATA PREPARATION OCCURS IN THE C TIME-DOMAIN BEFORE FORWARD TRANSFORM OR AFTER INVERSE C TRANSFORM. C C IF FORWARD TRANSFORM IS SELECTED: 1) A LEAST SQUARES TREND C WILL BE REMOVED; 2) THE PROFILE WILL THEN BE EXTENDED OR C TRUNCATED; 3) IT WILL BE WINDOWED; AND FINALLY 4) IT WILL BE C TRANSFORMED FROM THE TIME DOMAIN TO THE FREQUENCY DOMAIN VIA AN C FFT OR DFT. IF INVERSE IS SELECTED, THE INVERSE OF THE ABOVE C STEPS WILL BE TAKEN IN THE REVERSE SEQUENCE. THE ORDER OF C ARGUMENTS ROUGHLY REFLECTS THE ORDER OF STEPS TAKEN UPON C FORWARD TRANSFORM. C C IF THE FAST FOURIER TRANSFORM (FFT: TIME REQUIRED IS C PROPORTIONAL TO N*LOG(N)/LOG(2)) IS DESIRED, THE C FREQUENCY-DOMAIN ARRAY IS RESTRICTED IN LENGTH TO POWERS OF C TWO. OTHERWISE, THE DISCRETE FOURIER TRANSFORM (DFT: TIME C REQUIRED IS PROPORTIONAL TO N*N) WILL BE USED. SPIKES AND C DVALS MUST HAVE BEEN REMOVED BEFORE CALLING FFTDC2. C C THE FOURIER TRANSFORM SUBROUTINES (FORKDC AND DFTTDC) CALLED BY C THIS ROUTINE BOTH INCLUDE A SCALING FACTOR EQUAL TO 1/SQRT(N). C (WHERE N IS THE TOTAL NUMBER OF POINTS IN THE FUNCTION). THE C FOURIER SUMS IN EITHER DIRECTION (-IW OR +IW) ARE MULTIPLIED BY C THIS SCALING FACTOR. THE PURPOSE IS TO MAKE THE FORWARD AND C INVERSE TRANSFORMS MATHEMATICALLY INDISTIQUISHABLE (EXCEPT FOR C PHASE) AS DESCRIBED IN CLAERBOUT. THE FORWARD AND INVERSE C FUNCTIONS ARE ILLUSTRATED AS FOLLOWS: C C F.D. = (1/SQRT(N))*FSUM(-IW) T.D. = (1/SQRT(N))*FSUM(+IW) C C IN CONTRAST, THE TRUE DISCRETE FOURIER TRANSFORM APPEARS AS: C C F.D. = (1/N)*FSUM(-IW) T.D. = 1*FSUM(+IW) C C THEREFORE, IF ACTUAL FREQUENCY DOMAIN AMPLITUDES ARE REQUIRED C AFTER TRANSFORMING WITH THIS ROUTINE, ALL VALUES IN FQD (REAL C AND IMAGINARY) MUST BE MULTIPLIED BY 1/SQRT(NFQD). IF HOWEVER C A FORWARD TRANSFORM IS TO BE FOLLOWED BY AN INVERSE TRANSFORM, C THE SCALING FACTORS WILL PROPERLY MULTIPLY TO 1/N AFTER BOTH C TRANSFORMS, LEAVING THE AMPLITUDES IN THE FINAL TIME DOMAIN C UNCHANGED. (NOTE THAT IF THE DFT IS USED WITH DISPARATE C NUMBERS OF INPUT AND OUTPUT POINTS, THE SCALING FACTORS WILL C NOT BEHAVE IN A MEANINGFUL WAY). C C C C THE TIME DOMAIN C C IN THE CONTEXT OF A DISCRETE FOURIER TRANSFORM, THE TIME DOMAIN C IS AN ARRAY OF EVENLY-SPACED REAL (IMPOSED BY THIS SUBROUTINE) C VALUES COMPRISING A 1-DIMENSIONAL FUNCTION OF TIME. IN C FORTRAN, ARRAYS BEGIN WITH THE FIRST ELEMENT NUMBERED 1 (NOT C 0). THEREFORE, THE NUMBER OF POINTS IN THE ARRAY EQUALS THE C DIMENSION OF THE ARRAY. HOWEVER, IN THIS DESCRIPTION, IT WILL C BE ASSUMED THAT ARRAYS OF N ELEMENTS BEGIN WITH ELEMENT NUMBER C 0 (NOT 1) AND EXTEND TO ELEMENT NUMBER M = N-1. THE TIMES OF C THE POINTS ARE IMPLIED TO EXTEND FROM T0 TO TM: C C ARRAY OF FUNCTION VALUES X(T) = X0, X1, X2, ..., XM C TIME T = T0, T1, T2, ..., TM C C BECAUSE THE FOURIER TRANSFORM IS BASED ON PERIODIC FUNCTIONS, C IT ASSUMES THE INTERVAL OF DATA [X0...XM] TO REPEAT: C C FUNC VALS ... XM, X0, X1, X2, ..., XM, X0, X1, X2, ... C TIME ..., T-1, T0, T1, T2, ..., TM, TM+1, TM+2, TM+3, ... C C THEREFORE, AS FAR AS THE FOURIER TRANSFORM IS CONCERNED, THE C VALUE OF THE POINT IMMEDIATELY AFTER THE LAST POINT IS KNOWN C IMPLICITLY TO EQUAL THE VALUE OF THE FIRST POINT IN THE TIME C SERIES. (THERE IS OFTEN A LARGE STEP BETWEEN THESE POINTS C WHICH WOULD NOT REALLY EXIST IN THE TRUE EXTENDED FUNCTION. C HENCE, VARIOUS DATA PREPARATIONS ARE NEEDED TO CALM DOWN THE C IMPLIED DISCONTINUITY). C C THE IMPLIED POINT MUST BE INCLUDED WHENEVER SYMMETRIES ARE C CONSIDERED. FOR EXAMPLE, THE CENTER TIME IS NOT TM/2 BUT C RATHER (TM+1)/2. THIS LEADS TO THE CONCEPT OF DFT-EVEN C SYMMETRY WHICH REQUIRES THE CENTER LOCATION IN A SERIES TO BE C 1/2 OF AN INTERVAL TO THE RIGHT OF THE ACTUAL CENTER; IN OTHER C WORDS, THE END-POINT OF THE SERIES IS NOT M BUT RATHER M+1. C (IF THE THE ACTUAL DIMENSION IS N, THE IMPLIED DFT-EVEN C DIMENSION IS N+1.) C C THE DFT-EVEN CONCEPT IS USED HERE FOR TIME-DOMAIN WINDOWING BUT C IS ALSO HELPFUL IN UNDERSTANDING THE SYMMETRIES APPEARING IN C THE FREQUENCY DOMAIN. BECAUSE THE FOURIER TRANSFORM DOES NOT C DISTINGUISH BETWEEN FREQUENCY AND TIME, THE FREQUENCY DOMAIN C ALSO CONTAINS DFT-EVEN SYMMETRY. C C C C FREQUENCY DOMAIN C C JUST AS IN THE TIME DOMAIN, THE FREQUENCY DOMAIN IS PERIODIC. C BECAUSE OF THE PERIODICITY, IT IS USEFUL TO THINK OF THE C FREQUENCIES IN THE FREQUENCY DOMAIN AS ANGLES OR LOCATIONS ON A C UNIT CIRCLE. THE FULL RANGE OF TRANSFORMED VALUES WOULD THEN C BE REPRESENTED AFTER ONE TURN AROUND THE CIRCLE. IF THE ANGLE C W REPRESENTS A GIVEN FREQUENCY, THEN W, W+2PI, W+4PI, ETC ALL C REPRESENT THE SAME FREQUENCY BECAUSE THEY ALL ARE LOCATED AT C THE SAME PLACE ON THE CIRCLE. C C TO SIMPLIFY MATTERS, THE DOMAIN CAN BE CONFINED TO ONE COMPLETE C CIRCLE OR 2PI INTERVAL. IT IS CONVENIENT TO PICK THE INTERVAL C -PI < W <= +PI. IN THIS INTERVAL THE ABSOLUTE VALUES OF C FREQUENCIES ARE CONFINED TO 0 <= |W| <= +PI; THE MAXIMUM C FREQUENCY IS +PI. +PI IS ACTUALLY THE NYQUIST FREQUENCY; ANY C FREQUENCY ABOVE NYQUIST FOLDS BACK TO BECOME A LOWER FREQUENCY. C C THE TRANSFORM OF A REAL (TIME-DOMAIN) FUNCTION WILL CONSIST OF C A REAL EVEN FUNCTION AND AN IMAGINARY ODD FUNCTION; THE REAL C EVEN FUNCTION IS MIRRORED ABOUT 0 AND THE IMAGINARY ODD C FUNCTION IS MIRRORED AND INVERTED. THIS IMPLIES THAT WHILE THE C NEGATIVE FREQUENCIES (-PI TO 0) ARE IMPORTANT IN THE TRANSFORM, C THEY DO NOT CONTAIN ANY INFORMATION THAT CANNOT BE OBTAINED C FROM THE POSITIVE FREQUENCIES (0 TO PI). C C BECAUSE IT IS EASIER TO HANDLE ALL-POSITIVE ARRAY INDICIES, THE C VALUES IN THE NEGATIVE FREQUENCIES CAN BE COPIED INTO THE C POSITIVE INTERVAL, +PI < W < +2PI. THIS WOULD BE THE SAME AS C ROTATING THE INTERVAL A HALF CIRCLE (PI RADIANS) RE-DEFINING IT C 0 <= W < 2PI. IN THE NEW INTERVAL, NYQUIST FREQUENCY STILL C OCCURS AT PI. BUT BECAUSE OF FOLDING, THE EFFECTIVE FREQUENCY C DECREASES TO ZERO IN BOTH DIRECTIONS FROM NYQUIST WITH ZERO C FREQUENCY EFFECTIVELY OCCURING AT BOTH ZERO AND 2PI. C C THE DISCRETE FOURIER TRANSFORM (DFT) DOES NOT PRODUCE A C CONTINUUM OVER THE DOMAIN; BUT RATHER, IT PRODUCES A SET OF C VALUES AT DISCRETE FREQUENCIES EVENLY SPACED AROUND THE CIRCLE. C FURTHERMORE, IF ALL OF THE INFORMATION IN THE TIME-DOMAIN IS TO C BE TRANSFORMED TO THE FREQUENCY DOMAIN, IT IS NECESSARY AND C SUFFICIENT THAT THE FREQUENCY DOMAIN HAVE EXACTLY THE SAME C NUMBER OF POINTS AS THE TIME DOMAIN. IN OTHER WORDS, IF THE C TIME DOMAIN FUNCTION HAD N EVENLY-SPACED POINTS OVER THE TIME C INTERVAL 0 <= T < TN, THE FREQUENCY DOMAIN SHOULD HAVE N C POINTS EVENLY SPACED OVER THE FREQUENCY INTERVAL 0 <= W < 2PI. C C IT SHOULD BE NOTED THAT THE DFT-EVEN SYMMETRY EXISTS IN BOTH C DOMAINS. PARTICULARLY, THERE IS NO EXPLICIT POINT IN THE C FREQUENCY DOMAIN AT 2PI; THE VALUE AT W=2PI IS IMPLIED EQUAL TO C THE VALUE AT W=0. THE LEFT-MOST POINT IS AT ZERO FREQUENCY; C THE REST OF THE POINTS ARE EVENLY SPACED OVER THE INTERVAL 0 TO C 2PI WITH THE POINT AT 2PI MERELY IMPLIED. C C THE ANGULAR FREQUENCY (NOT ACCOUNTING FOR FOLDING) OF THE KTH C POINT IS: C C WK = 2PI*K/N K = 0, 1, 2, ..., N-1. C C THE TEMPORAL FREQUENCY (ACCOUNTING FOR FOLDING) IS: C C FK = WK /(2PI*DT) = K /(N*DT) K=0,1,2,...N/2 C FK = (2PI-WK)/(2PI*DT) = (N-K)/(N*DT) K=N/2+1,...N-1 C C WHERE DT = THE TEMPORAL SAMPLING INTERVAL C C IF THERE IS AN EVEN NUMBER OF POINTS IN THE FUNCTION, THE C NYQUIST FREQUENCY WILL HAVE A CORROSPONDING POINT AT K=N/2+1. C IF HOWEVER, THERE IS AN ODD NUMBER OF POINTS, THE NYQUIST C FREQUENCY WILL FALL HALF WAY BETWEEN K=(N+1)/2 AND K=(N+1)/2+1 C C C C REFERENCES C C A DISCUSSION OF THE FFT MAY BE FOUND IN "FUNDAMENTALS OF C GEOPHYSICAL DATA PROCESSING" BY JON F. CLAERBOUT (ISBN C 0-86542-305-9), CHAPTER 1, PAGES 6-12. (THE OPERATION OF THE C DFT MAY ALSO BE UNDERSTOOD FROM THIS DISCUSSION.) C C A DISCUSSION OF THE DFT MAY BE FOUND IN "FOURIER ANALYSIS OF C TIME SERIES: AN INTRODUCTION" BY PETER BLOOMFIELD, JOHN WILEY C & SONS, PAGES 46-49. C C AN EXCELLENT IN-DEPTH DISCUSSION OF THE FOURIER TRANSFORM MAY C BE FOUND IN "THE FOURIER TRANSFORM AND ITS APPLICATIONS", BY C RONALD N. BRACEWELL, MCGRAW-HILL. C C A RATHER INTENSE DISCUSSION OF THE BURG ERROR-PREDICTION FILTER C USED IN EXTENSION MAY BE FOUND IN CLAERBOUT CHAPTER 7-2, PAGES C 133-137. C C AN IN-DEPTH DISCUSSION OF WINDOWING FUNCTIONS MAY BE FOUND IN C "WINDOWS, HARMONIC ANALYSIS, AND THE DISCRETE FOURIER C TRANSFORM" BY FREDRIC J. HARRIS, UNDERSEA SURVEILLANCE C DEPARTMENT, SEPTEMBER 1976. C C C C DIRECTION OF TRANSFORM C C INPUT ARGUMENT: C C NVERSE - INTEGER*4. DIRECTION OF TRANSFORM. NVERSE <= 0 IS C FORWARD TRANSFORM (DEFINED HERE AS TIME-DOMAIN TO C FREQUENCY DOMAIN). NVERSE => 1 IS INVERSE TRANSFORM C (FREQUENCY-DOMAIN TO TIME-DOMAIN). C C MOST ARGUMENTS TO THIS SUBROUTINE ARE AT ALL TIMES C INPUT ONLY. ONE ARGUMENT (YWIND) IS AT ALL TIMES C OUTPUT ONLY. THREE OF THE ARGUMENTS (ALSQ, TMD, AND C FQD) CHANGE AMONG INPUT ONLY, OUTPUT ONLY, OR C INPUT/OUTPUT DEPENDING UPON THE VALUE OF NVERSE. C THIS IS ILLUSTRATED IN THE FOLLOWING TABLE: C C ARGUMENT FORWARD INVERSE C (NVERSE<=0) (NVERSE=>1) C C (MOST ARGS) INPUT ONLY INPUT ONLY C ALSQ OUTPUT ONLY INPUT ONLY C YWIND OUTPUT ONLY OUTPUT ONLY C TMD INPUT/OUTPUT OUTPUT ONLY C FQD OUTPUT ONLY INPUT/OUTPUT C C C C TREND REMOVAL C C INPUT ARGUMENTS: C C KTREND - INTEGER*4. NUMBER INDICATING THE TYPE OF TREND TO C REMOVE. A STAR (*) INDICATES THAT THE OPTION HAS C BEEN INSTALLED: C C KTREND TYPE OF TREND C * 0 = NO TREND REMOVAL C * 1 = POLYNOMIAL OF ORDER NLSQ-1 C 2 = COSINE CURVE C C NLSQ - INTEGER*4. NUMBER OF CONSTANTS ASSOCIATED WITH THE C DESIRED LEAST-SQUARES-TREND REMOVAL. POLYNOMIAL C FUNCTIONS TO ORDER 9 ARE INSTALLED; HOWEVER GREATER C THAN ORDER 3 IS NOT RECOMMENDED. TYPICAL VALUES ARE C AS FOLLOWS: C C NLSQ ORDER OF TREND C 0 = (NO TREND REMOVAL IF POLYNOMIAL) C 1 = ZEROTH ORDER (AVERAGE IF POLYNOMIAL) C 2 = FIRST ORDER OR (LINEAR IF POLYNOMIAL) C 3 = SECOND ORDER (PARABOLIC IF POLYNOMIAL) C 4 = THIRD ORDER C C OUTPUT OR INPUT ARGUMENT: C C ALSQ - REAL*8. ARRAY OF DIMENSION NLSQ CONTAINING THE C COEFFICIENTS OF THE TREND. IF FORWARD TRANSFORM, C ALSQ WILL BECOME AN OUTPUT ARGUMENT; IF INVERSE C TRANSFORM, ALSQ WILL BE AN INPUT ARGUMENT. C C C C EXTENSION AND TRUNCATION C C INPUT ARGUMENTS: C C KEXTEN - INTEGER*4. NUMBER INDICATING THE TYPE OF EXTENSION C THAT IS DESIRED. IF EXTENSION WITH ZEROS IS C SELECTED, WINDOWING (THE DATA PREPARATION STEP C FOLLOWING EXTENSION) WILL BE PERFORMED ONLY OVER THE C ORIGINAL FUNCTION. FOR ALL OTHER EXTENSION TYPES, C THE WINDOWING WILL BE PERFORMED OVER THE ENTIRE C EXTENDED LENGTH. A STAR (*) INDICATES THAT THE C OPTION HAS BEEN INSTALLED: C C KEXTEN TYPE OF EXTENSION C * -1 = EXTENSION/TRUNCATION DISALLOWED C * 0 = ZEROS C * 1 = AVERAGE VALUE C * 2 = FIRST/LAST VALUE C * 3 = REPEAT ORIGINAL FUNCTION C * 4 = REFLECT ORIGINAL FUNCTION C 5 = REFLECT AND INVERT (XPARMS NEEDED) C 6 = SPLINE TO ZERO (XPARMS NEEDED) C * 7 = BURG PREDICTION WITH A LINEAR FUNCTION C REMOVED BEFORE PREDICTION AND ADDED BACK C AFTER PREDICTION. ONE OR MORE VALUES MUST C BE SPECIFIED IN THE ARGUMENT, XPARMS(). C C FOR KEXTEN=7, XPARMS(1) SPECIFIES THE ROOT C LENGTH OF THE BURG FILTER. THE ROOT LENGTH C IS THE NUMBER OF DATA POINTS COUNTED FROM C THE FIRST PREDICTED POINT BACKWARDS INTO C THE EXISTING PROFILE. THESE POINTS ARE C USED TO PRODUCE THE PREDICTION-FILTER C COEFFICIENTS. FOR EXAMPLE, IF XPARMS(1)=3, C AND THE TIME-DOMAIN PROFILE HAS NTMD=20 C POINTS, THEN POINTS 18, 19, AND 20 WILL BE C USED TO PREDICT THE NEXT POINT, 21. THEN, C POINT 22 WILL BE PREDICTED FROM POINTS 19, C 20, AND 21; AND SO FORTH. THIS PROCESS IS C REFLECTED AT THE LEFT END OF THE PROFILE, C IF EXTENSION IS OCCURRING THERE. C C BECAUSE THE ROOT LENGTH CANNOT BE LESS THAN C 2, VALUES OF XPARMS(1) THAT ARE LESS THAN 2 C ARE USED TO INDICATE SPECIFIC METHODS OF C DERIVING THE ROOT LENGTH. THESE METHODS C ARE GIVEN IN THE FOLLOWING TABLE, AND SOME C REQUIRE ADDITIONAL VALUES TO BE SPECIFIED C IN XPARMS(2) AND XPARMS(3): C C XPARMS(1) ROOT-LENGTH CALCULATION METHOD C -3 LEFT ROOT = XPARMS(2) C RIGHT ROOT = XPARMS(3) C -2 LEFT ROOT = C XPARMS(2)*(LEFT EXTENSION) C RIGHT ROOT = C XPARMS(2)*(RIGHT EXTENSION) C -1 LEFT ROOT = LEFT EXTENSION C RIGHT ROOT = RIGHT EXTENSION C 0 EACH ROOT = INPUT ARRAY LEN (NTMD) C 1 EACH ROOT = C HALF THE TOTAL EXTENDED LENGTH C WHICH IS (NFQD-NTMD)/2 C >1 EACH ROOT = XPARMS(1) C C * 8 = BURG PREDICTION AS IN KEXTEN=7 ONLY WITHOUT C THE LINEAR FUNCTION. C C NOTE: WHETHER THE PROFILE IS EXTENDED OR TRUNCATED C DEPENDS ON WHETHER NTMD IS LESS THAN OR GREATER C THAN NFQD; THE AMOUNT OF EXTENSION IS EQUAL TO THE C DIFFERENCE BETWEEN NTMD AND NFQD; AND THE LOCATION C WITHIN THE EXTENDED PROFILE IS GOVERNED BY THE C JUSTIFICATION, JUSTX. C C NOTE: USE OF THE BURG PREDICTION IS AN EXCELLENT C WAY OF EXTENDING THE PROFILE TO ANY ARBITRARY LENGTH C (UP TO 4 TIMES THE ORIGINAL PROFILE) WHILE C PRESERVING EXACTLY THE FREQUENCY CONTENT. UNLIKE C THE FOURIER TRANSFORM, BURG PREDICTION IS IMMUNE TO C THE REPETITION DISCONTINUITY ASSOCIATED WITH THE C BEGINNING AND END OF THE PROFILE. THE RESULT IS TO C NARROW SPECTRAL PEAKS TO AN ARBITRARILY THIN LINE. C WHEN USED IN CONJUNCTION WITH THE KAISER-BESSEL C WINDOW, RESULTS CAN APPROACH AN IDEAL NARROWING OF C THE PEAK AND REMOVAL OF SIDE LOBES. C C NOTE: THE BURG PREDICTION IS AN N*N PROCEDURE AND C HENCE IS QUITE SLOW. TO REDUCE THE TIME REQUIRED, A C SUITABLY SHORT "ROOT" LENGTH SHOULD BE CHOSEN USING C XPARMS. KEEP IN MIND HOWEVER, THAT SHORTENING THE C "ROOT" WILL CAUSE POORER FREQUENCY REPRODUCTION IN C THE EXTENDED SECTION. A RULE OF THUMB IS THAT THE C "ROOT" SHOULD NOT BE LESS THAN HALF OF THE AMOUNT TO C BE EXTENDED. C C XPARMS - REAL*8. ARRAY CONTAINING ANY PARAMETERS ASSOCIATED C WITH ONE OF THE KEXTEN NUMBERS. THIS ARGUMENT WILL C NOT BE USED UNLESS DESCRIBED UNDER KEXTEN. C C JUSTX - INTEGER*4. JUSTIFICATION OF THE ORIGINAL PROFILE C WITHIN THE EXTENDED PROFILE. FOLLOWING ARE POSSIBLE C VALUES: C 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 C C C WINDOWING C C INPUT ARGUMENTS: C C KWINDO - INTEGER*4. NUMBER INDICATING THE TYPE OF WINDOW. C EVERY WINDOW WILL BE APPLIED AS DFT-EVEN. THAT IS, C THE END-POINT SYMMETRIC WITH THE BEGINNING POINT IS C ASSUMED TO EXIST IMPLICITLY AFTER THE LAST C TIME-DOMAIN ARRAY POINT. AN EXCEPTION OCCURS WHEN C THE TIME-DOMAIN PROFILE HAS BEEN EXTENDED WITH C ZEROS; IN THIS CASE, A SYMMETRIC WINDOW IS APPLIED C ONLY TO THE NON-EXTENDED PART OF THE DATA. A STAR C (*) INDICATES THAT THE OPTION HAS BEEN INSTALLED: C C KWINDO 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, KWINDO, 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 (KWINDO=0) CANNOT BE C SHORTENED USING WFRAC; THAT PROCEDURE MUST BE DONE C WITH ZERO EXTENSION (KTYPE=0 IMPLIES ZERO EXTENSION, C OR USE KTYPE=4, RIESZ WINDOW, AND ALPHA NEARLY 0). C C OUTPUT ARGUMENT: C C YWIND - REAL*8. ARRAY OF DIMENSION NTMD CONTAINING THE C WINDOWING VALUES NORMALIZED TO 1 AT THE MAXIMUM. C YWIND IS INTENDED TO SERVE AS A WEIGHTING FUNCTION C FOR RESTORING CONTIGUOUS DATA THAT WAS STACKED AND C FILTERED USING OVERLAPPING WINDOWING. (IF EXTENSION C IS BEING PERFORMED AND NFQD IS GREATER THAN NTMD AND C A WINDOW IS BEING APPLIED OR UNAPPLIED, THEN C WORKSPACE MUST BE SUPPLIED IN YWIND FROM LOCATION C NTMD+1 THROUGH NFQD. THEREFORE, ALTHOUGH THE USEFUL C INFORMATION IN YWIND IS ALWAYS IN ELEMENTS 1-NTMD, C YWIND SHOULD BE DIMENSIONED THE GREATER OF NTMD AND C NFQD.) C C C C TRANSFORM AND DATA C C INPUT ARGUMENTS: C C KDXFRM - INTEGER*4. THE TYPE OF TRANSFORM DESIRED. C C KDXFRM DESCRIPTION C -1 FOURIER TRANSFORM NOT PERFORMED C 0 FFT WITH SCALE FACTOR (IF POSSIBLE) C 1 DFT WITH SCALE FACTOR C 2 FFT WITHOUT SCALE FACTOR (IF POSSIBLE) C 3 DFT WITHOUT SCALE FACTOR C C IF KDXFRM IS NEGATIVE, THE TRANSFORM STEP WILL NOT C BE PERFORMED. INSTEAD, THE TIME-DOMAIN DATA (AFTER C TREND REMOVAL, EXTENSION, AND WINDOWING) WILL BE C COPIED INTO THE REAL PART OF FQD. THE IMAGINARY C PART WILL BE SET TO ZERO. THIS IS EXACTLY AS IT C WOULD HAVE APPEARED AT THE INPUT TO A FOURIER C TRANSFORM ROUTINE. IF NVERSE IS GREATER THAN ZERO, C THE REAL PART OF FQD WILL BE COPIED DIRECTLY INTO C TMD AND THE IMAGINARY PART DROPPED. (THEN C UNWINDOWING, TRUNCATION, AND TREND ADD-BACK WILL BE C PERFORMED.) C C IF KDXFRM IS EQUAL TO ZERO, THE FFT WILL BE USED C WHEN POSSIBLE (NFQD MUST BE A POWER OF 2); WHEN NOT C POSSIBLE, THE DFT WILL BE SUBSTITUTED AUTOMATICALLY. C C IF KDXFRM IS GREATER THAN ZERO, THE DFT WILL BE C USED. THIS DFT IS DESIGNED TO HAVE THE SAME SCALING C FACTOR AS FFT. THEREFORE WHEN EITHER FFT OR DFT CAN C BE USED THEIR RESULTS ARE INDISTINGUISHABLE. NOTE, C HOWEVER THAT THE DFT IS CONSIDERABLY SLOWER THAN THE C FFT TAKING SEVERAL SECONDS FOR PROFILES AROUND 1000 C POINTS AND INCREASING BY THE SQUARE OF N. C C IF THE DFT IS SELECTED, IT IS POSSIBLE TO HAVE C DIFFERENT LENGTH TIME AND FREQUENCY DOMAINS C (OCCURING WHEN EXTENSION/TRUNCATION HAS BEEN C DISALLOWED AND NTMD IS NOT EQUAL TO NFQD). THIS C OPTION MAY HOWEVER BE OF LIMITED UTILITY BECAUSE THE C TRANSFORM IS RENDERED INCOMPLETE BY THE LENGTH C DISPARITY. C C NTMD - INTEGER*4. THE NUMBER OF POINTS IN THE TIME-DOMAIN C ARRAY, TMD. NTMD MAY BE ANY REASONABLE NUMBER AND C IS NOT RESTRICTED TO POWERS OF 2. (IF THE FFT IS C DESIRED FOR EITHER FORWARD OR INVERSE, EXTENSION C MUST BE ENABLED AND NFQD MUST BE A POWER OF 2). C C NOTE: WHETHER THE PROFILE IS EXTENDED OR TRUNCATED C DEPENDS ON WHETHER NTMD IS LESS THAN OR GREATER C THAN NFQD; THE AMOUNT OF EXTENSION IS EQUAL TO THE C DIFFERENCE BETWEEN NTMD AND NFQD; AND THE LOCATION C WITHIN THE EXTENDED PROFILE IS GOVERNED BY THE C JUSTIFICATION, JUSTX. C C INPUT/OUTPUT OR OUTPUT ARGUMENT: C C TMD - REAL*8. SINGLE DIMENSION ARRAY CONTAINING THE C TIME-DOMAIN FUNCTION. TMD MUST BE LARGE ENOUGH TO C ACCOMODATE THE GREATER OF NTMD AND NFQD. C C IF FORWARD TRANSFORM, AT THE TIME OF CALLING TMD C MUST CONTAIN A TIME-DOMAIN FUNCTION OF LENGTH NTMD. C THE CONTENTS OF THE ARRAY WILL BE MODIFIED DURING C PROCESSING TO INCLUDE: TREND REMOVAL, EXTENSION (OR C TRUNCATION), AND WINDOWING. THEREFORE, TMD CAN BE C VIEWED AS AN INPUT/OUTPUT ARRAY. C C IF INVERSE TRANSFORM, TMD WILL BE AN OUTPUT-ONLY C ARRAY CONTAINING THE FINAL TIME-DOMAIN FUNCTION OF C LENGTH NTMD. HOWEVER, DURING PROCESSING IT WILL BE C USED AS A WORK ARRAY FOR UNWINDOWING, TRUNCATION (OR C EXTENSION), AND TREND ADD-BACK. C C C INPUT ARGUMENT: C C NFQD - INTEGER*4. THE NUMBER OF POINTS IN THE C FREQUENCY-DOMAIN ARRAY, FQD. NFQD MAY BE ANY C REASONABLE NUMBER AND IS NOT RESTRICTED TO POWERS OF C 2. HOWEVER, NFQD MUST BE A POWER OF 2 IF THE FFT IS C DESIRED FOR EITHER FORWARD OR INVERSE TRANSFORM. IF C EXTENSION IS DISABLED AND NFQD DOES NOT EQUAL NTMD, C THE DFT WILL BE USED REGARDLESS OF WHETHER NFQD IS A C POWER OF 2. C C NORMALLY, THE FREQUENCY-DOMAIN PROFILE WOULD BE C EQUAL OR LONGER THAN THE TIME-DOMAIN. HOWEVER, C CERTAIN OPERATIONS MAY REQUIRE IT TO BE SHORTER. C THIS CASE CAN BE ACCOMODATED SIMPLY BY MAKING NFQD C SMALLER THAN NTMD. IF THIS OCCURS, TRUNCATION WILL C BE PERFORMED IN THE TIME-DOMAIN DURING FORWARD C TRANSFORM, AND EXTENSION DURING INVERSE TRANSFORM. C C NOTE: WHETHER THE PROFILE IS EXTENDED OR TRUNCATED C DEPENDS ON WHETHER NTMD IS LESS THAN OR GREATER C THAN NFQD; THE AMOUNT OF EXTENSION IS EQUAL TO THE C DIFFERENCE BETWEEN NTMD AND NFQD; AND THE LOCATION C WITHIN THE EXTENDED PROFILE IS GOVERNED BY THE C JUSTIFICATION, JUSTX. C C OUTPUT OR INPUT/OUTPUT ARGUMENT: C C FQD - COMPLEX*16. SINGLE DIMENSION ARRAY CONTAINING THE C FREQUENCY-DOMAIN FUNCTION. FQD MUST BE LARGE ENOUGH C TO ACCOMODATE THE GREATER OF NTMD AND NFQD. C C IF FORWARD TRANSFORM, FQD WILL BE AN OUTPUT-ONLY C ARRAY CONTAINING THE FINAL FREQUENCY-DOMAIN FUNCTION C OF LENGTH NFQD. C C IF INVERSE TRANSFORM, AT THE TIME OF CALLING FQD C MUST CONTAIN A FREQUENCY-DOMAIN FUNCTION OF LENGTH C NFQD. THE CONTENTS OF THE ARRAY WILL BE TRANSFORMED C BACK TO THE TIME DOMAIN DURING PROCESSING. C THEREFORE, FQD CAN BE VIEWED AS AN INPUT/OUTPUT C ARRAY. C C C SUBROUTINE FFTDCX WRITTEN BY ROB BRACKEN, USGS, 23AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 23AUG96. C SUBROUTINE FFTDC2 WRITTEN BY ROB BRACKEN, USGS, 18OCT96. C HP-9000 SERIES 700/800 UNIX VERSION 2.0, 18NOV96. C VERSION 2.1, 20020824 ( ADDED DFTFDR CAPABILITY & ARG WFRAC ). C C subroutine fftdc2(nverse, & ktrend,nlsq,alsq, & kexten,xparms,justx, & kwindo,alpha,wfrac, ywind, & kdxfrm, ntmd,tmd, nfqd,fqd) C C DECLARATIONS C C ARGUMENTS C C FORWARD (-1)/INVERSE (+1) integer*4 nverse C C TREND REMOVAL/ADDITION integer*4 ktrend,nlsq real*8 alsq(*) C C FUNCTION EXTENSION/TRUNCATION integer*4 kexten real*8 xparms(*) integer*4 justx C C WINDOW APPLICATION/REMOVAL integer*4 kwindo real*8 alpha,wfrac,ywind(*) C C TRANSFORM-TYPE SELECTION (FFT OR DFT) integer*4 kdxfrm C C TIME-DOMAIN FUNCTION ARRAY integer*4 ntmd real*8 tmd(*) C C FREQUENCY-DOMAIN FUNCTION ARRAY integer*4 nfqd complex*16 fqd(*) C C INTERNAL VARIABLES C real*8 signi,scale C c wfrac=1.d0 C C DETERMINE WHETHER EXTENSION/TRUNCATION IS ALLOWED C C DEFAULT IS EXTENSION/TUNCATION NOT ALLOWED jxt=0 C ALLOWED IF EXTENSION IS NOT BARRED & NTMD DIFFERS FROM NFQD if(kexten.ge.0.and.ntmd.ne.nfqd) jxt=1 C C DETERMINE IF SYMMETRIC WINDOWING IN MIDDLE OF PROFILE IS NEEDED C C DEFAULT IS DFT-EVEN WINDOWING OVER ENTIRE LENGTH OF PROFILE jw0=0 C SYMMETRIC WINDOWING IF EXTENSION WITH ZEROS AND NTMD < NFQD if(jxt.eq.1.and.kexten.eq.0.and.ntmd.lt.nfqd) jw0=1 C C DETERMINE WHETHER TO USE FFT OR DFT C C FFT IS DEFAULT jdx=0 C DFT IS NEEDED IF EXT NOT ALLOWED & NTMD DIFFERS FROM NFQD if(jxt.eq.0.and.ntmd.ne.nfqd) jdx=2 C DFT IS NEEDED IF NFQD IS NOT A POWER OF 2 nchk=1 205 nchk=nchk*2 if(nchk.lt.nfqd) goto 205 if(nchk.gt.nfqd.and.jdx.eq.0) jdx=1 C DFT IF SPECIFICALLY REQUESTED if((kdxfrm.eq.1.or.kdxfrm.eq.3).and.jdx.eq.0) jdx=1 C NO TRANSFORM IF SPECIFICALLY REQUESTED if(kdxfrm.le.-1) jdx=-1 C C FIND BEG LOCATION OF INPUT PROFILE WITHIN EXTENDED PROFILE C inloc=justx if(justx.lt.0) inloc=nfqd-ntmd+justx+2 if(justx.eq.0) inloc=(nfqd-ntmd)/2+1 C C INITIALIZE VALUES IN NORMALIZED WINDOW FUNCTION C do 704 i=1,ntmd ywind(i)=1.d0 704 continue C C CHECK DIRECTION OF TRANSFORM C if(nverse.ge.1) goto 101 C C FORWARD TRANSFORM C C REMOVE TREND if(ktrend.ge.1) & call trendr(ktrend,nlsq,alsq, nverse,ntmd,tmd) C C EXTEND if(jxt.eq.0) then npts=ntmd else call extend(kexten,xparms,justx,ntmd,nfqd,tmd) npts=nfqd endif C C APPLY WINDOW if(kwindo.ge.1) then if(jw0.eq.0) then C C WINDOW ENTIRE PROFILE kdftev=1 call window(kwindo,alpha,wfrac,kdftev, nverse,npts,tmd) C C SET UP THE NORMALIZED WINDOW FUNCTION do 706 i=ntmd+1,npts ywind(i)=1.d0 706 continue call window(kwindo,alpha,wfrac,kdftev, & nverse,npts,ywind) ixwind=0 call extend(ixwind,xparms,justx,npts,ntmd,ywind) else C C WINDOW NON-ZERO PORTION OF THE EXTENDED PROFILE kdftev=0 call window(kwindo,alpha,wfrac,kdftev, & nverse,ntmd,tmd(inloc)) C C SET UP THE NORMALIZED WINDOW FUNCTION call window(kwindo,alpha,wfrac,kdftev, & nverse,ntmd,ywind) endif endif C C FORWARD TRANSFORM signi=-1.d0 do i=1,npts fqd(i)=dcmplx(tmd(i),0.d0) enddo if(jdx.eq.0) then C C USE FFT (SCALE FACTOR INCLUDED IN FORKDC) call forkdc(npts,fqd,signi) if(kdxfrm.eq.2) then C C REMOVE SCALE FACTOR scale=1.d0/dsqrt(dflotj(npts)) do i=1,npts fqd(i)=fqd(i)*scale enddo endif else if(jdx.eq.1) then C C USE DFT (SCALE FACTOR NOT INCLUDED DFTTDC) call dfttdc(npts,nfqd,signi,fqd) if(kdxfrm.ne.2.and.kdxfrm.ne.3) then C C ATTACH SCALE FACTOR scale=dsqrt(dflotj(nfqd)) do i=1,nfqd fqd(i)=fqd(i)*scale enddo endif else if(jdx.eq.2) then C C USE DFT WITH FREQUENCY DOMAIN ZERO-EXTENSION/TRUNCATION C (SCALE FACTOR NOT INCLUDED DFTFDR) call dftfdr(npts,nfqd,signi,tmd,fqd) if(kdxfrm.ne.2.and.kdxfrm.ne.3) then C C ATTACH SCALE FACTOR scale=dsqrt(dflotj(nfqd)) do i=1,nfqd fqd(i)=fqd(i)*scale enddo endif else C C USE NO TRANSFORM if(nfqd.gt.npts) then do i=npts+1,nfqd fqd(i)=dcmplx(0.d0,0.d0) enddo endif endif goto 990 C C INVERSE TRANSFORM C C INVERSE FFT 101 signi=1.d0 npts=nfqd if(jxt.eq.0) npts=ntmd C if(jdx.eq.0) then C C USE FFT (SCALE FACTOR INCLUDED IN FORKDC) call forkdc(npts,fqd,signi) if(kdxfrm.eq.2) then C C REMOVE SCALE FACTOR scale=dsqrt(dflotj(npts)) do i=1,npts fqd(i)=fqd(i)*scale enddo endif else if(jdx.eq.1) then C C USE DFT (SCALE FACTOR NOT INCLUDED IN DFTTDC) call dfttdc(nfqd,npts,signi,fqd) if(kdxfrm.ne.2.and.kdxfrm.ne.3) then C C ATTACH SCALE FACTOR scale=1.d0/dsqrt(dflotj(nfqd)) do i=1,npts fqd(i)=fqd(i)*scale enddo endif else if(jdx.eq.2) then C C USE DFT WITH FREQUENCY DOMAIN ZERO-EXTENSION/TRUNCATION C (SCALE FACTOR NOT INCLUDED DFTFDR) call dftfdr(nfqd,npts,signi,tmd,fqd) do 708 i=1,npts fqd(i)=dcmplx(tmd(i),0.d0) 708 continue C if(kdxfrm.ne.2.and.kdxfrm.ne.3) then C C ATTACH SCALE FACTOR scale=1.d0/dsqrt(dflotj(nfqd)) do i=1,npts fqd(i)=fqd(i)*scale enddo endif else C C USE NO TRANSFORM if(nfqd.lt.npts) then do i=nfqd+1,npts fqd(i)=dcmplx(0.d0,0.d0) enddo endif endif do 703 i=1,npts tmd(i)=dreal(fqd(i)) 703 continue C C REMOVE WINDOW if(kwindo.ge.1) then if(jw0.eq.0) then C C UNWINDOW ENTIRE PROFILE kdftev=1 call window(kwindo,alpha,wfrac,kdftev, nverse,npts,tmd) C C SET UP THE NORMALIZED WINDOW FUNCTION do 707 i=ntmd+1,npts ywind(i)=1.d0 707 continue call window(kwindo,alpha,wfrac,kdftev, & -nverse,npts,ywind) ixwind=0 call extend(ixwind,xparms,justx,npts,ntmd,ywind) else C C UNWINDOW NON-ZERO PORTION OF THE EXTENDED PROFILE kdftev=0 call window(kwindo,alpha,wfrac,kdftev, & nverse,ntmd,tmd(inloc)) C C SET UP THE NORMALIZED WINDOW FUNCTION call window(kwindo,alpha,wfrac,kdftev, & -nverse,ntmd,ywind) endif endif C C TRUNCATE if(jxt.ne.0) & call extend(kexten,xparms,justx,nfqd,ntmd,tmd) C C ADD BACK TREND if(ktrend.ge.1) & call trendr(ktrend,nlsq,alsq, nverse,ntmd,tmd) C C EXIT PROCEDURE C 990 return end C C________________________________________________________________ C C SUBROUTINE F O R K D C C________________________________________________________________ C C SUBROUTINE FORKDC PERFORMS A FAST FOURIER TRANSFORM (FFT) ON A C SAMPLE FUNCTION (SEE NOTE BELOW ON TERMINOLOGY) IN A C DOUBLE-PRECISION COMPLEX ARRAY. THE TRANSFORM MAY BE PERFORMED C EITHER FORWARD (TIME TO FREQUENCEY DOMAIN) OR INVERSE C (FREQUENCY TO TIME DOMAIN). THE INPUT/OUTPUT ARRAY LENGTH IS C RESTRICTED TO POWERS OF 2 (SEE NOTE BELOW ON POWERS OF 2). C C THE FOLLOWING TABLE SHOWS TIMES REQUIRED FOR VARIOUS ROUTINES C TO PERFORM A FORWARD AND INVERSE TRANSFORM. (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 C C THE FFT ALGORITHM WAS FIRST PUBLICIZED BY COOLEY AND TUKEY IN C 1965. THE FORTRAN CODING, SUBROUTINE FORK (2/15/69), WAS C MODIFIED FROM BRENNER BY JON F. CLAERBOUT, DEPARTMENT OF C GEOPHYSICS, STANFORD UNIVERSITY, AND PUBLISHED IN "FUNDAMENTALS C OF GEOPHYSICAL DATA PROCESSING: WITH APPLICATIONS TO PETROLEUM C PROSPECTING", 1976 MCGRAW-HILL, PAGE 12. IT WAS LATER ADAPTED C TO A USGS COMPUTER BY MIKE WEBRING, THEN CONVERTED TO DOUBLE C COMPLEX ARITHMETIC (FORKDC) AND NOTES SUPPLIED BY ROB BRACKEN C IN 1996. THE DFT PROCEDURE IS SUMMARIZED AS FOLLOWS (THE FFT C IS A STREAMLINED VERSION WHICH TAKES ADVANTAGE OF SYMMETRIES IN C A MATRIX VERSION OF THE FOLLOWING EQUATION): C C NX C CX(K)=SQRT(1/NX)* SUM (CX(J)*EXP(2*PI*SIGNI*I*(J-1)*(K-1)/NX)) C J=1 C WHERE K = 1,2, ... (NX=2**INTEGER) C SIGNI = -1 (FORWARD), C = 1 (INVERSE) C C NOTE ON TERMINOLOGY: A SAMPLE FUNCTION (IN THIS CONTEXT) IS C THE TIME-DOMAIN FUNCTION OF INTEREST WHICH IS BEING INPUT TO C THE FFT ALGORITHM. IT IS USUALLY A FINITE-LENGTH PIECE OF A C LONGER FUNCTION AND IS, HOPEFULLY, REPRESENTATIVE OF THAT C FUNCTION IN TERMS OF SPECTRAL CONTENT. THE FOURIER TRANSFORM C CONVERTS A TIME-DOMAIN FUNCTION INTO A FREQUENCY DOMAIN C SPECTRUM AND VICE-VERSA. THE ACTUAL FOURIER TRANSFORM HAS C INFINITE LIMITS OF INTEGRATION; THAT IS, IT REQUIRES A SAMPLE C FUNCTION OF INFINITE LENGTH. A VERSION OF THE FOURIER C TRANSFORM WHICH HAS FINITE LIMITS (AND THUS USES A PHYSICALLY C REALIZABLE SAMPLE FUNCTION) IS CALLED THE FINITE FOURIER C TRANSFORM. THE FINITE FOURIER TRANSFORM REQUIRES A C CONTINUOUSLY VARYING (OR ANALOG) SAMPLE FUNCTION WHICH CANNOT C BE SUPPLIED IN A DIGITAL COMPUTER. A DIGITAL VERSION IS THE C DISCRETE FOURIER TRANSFORM (DFT) WHICH USES AN EVENLY SAMPLED C (DISCRETE) SAMPLE FUNCTION AND PRODUCES SPECTRA WITH DISCRETE C FREQUENCIES. THE FAST FOURIER TRANSFORM IS A VERSION OF THE C DFT WHICH REDUCES COMPUTATION TIME. C C NOTE ON DATA PREPARATION: BEFORE TRANSFORMING A SAMPLE C FUNCTION, THE DATA SHOULD BE PREPARED SO AS TO MINIMIZE NOISE C AND FOCUS THE SPECTRUM APPROPRIATELY. THE PREPARATION SHOULD C INCLUDE THE FOLLOWING OPERATIONS: 1) REMOVE SPIKES, STEPS, AND C OTHER EXTRANEOUS NOISES WHICH ARE NOT OF INTEREST; 2) REMOVE A C TREND SUCH AS A LEAST-SQUARES-FIT LINE (KEEP THE COEFFICIENTS C IF THE LINE IS TO BE ADDED BACK LATER); AND 3) SHAPE THE DATA C WITH A WINDOWING FUNCTION. C C NOTE ON POWER OF 2 RESTRICTION: IF THE AVAILABLE FUNCTION IS C NOT AN INTEGER POWER OF 2 IN LENGTH AND IT CANNOT BE TRUNCATED, C IT MAY BE PADDED WITH ZEROS TO THE NEXT POWER OF 2. THIS C PROCEDURE IS GENERALLY ACCEPTED AS THE BEST THAT CAN BE DONE C WITH NON-POWER-OF-TWO DATA IN AN FFT. HOWEVER, BE WARNED THAT C IT MAY INCREASE BACKGROUND NOISE BY 10 dB OR MORE. THE NOISE C INCREASES IN PROPORTION TO THE PERCENTAGE OF PADDING. C C NOTE ON PERODIC EXTENSION: THE FINITE FOURIER TRANSFORM (AND C ALL SUBSETS) ASSUMES PERIODIC EXTENSION BEYOND IT'S LIMITS OF C INTEGRATION. IT EXTENDS THE SAMPLE RECORD IN BOTH DIRECTIONS C (POSITIVE AND NEGATIVE) BY REPLICATION. FOR EXAMPLE, AT THE C END OF THE ACTUAL SAMPLE FUNCTION, AN IDENTICAL COPY OF THE C ORIGINAL SAMPLE FUNCTION STARTS ALL OVER AGAIN. THIS C REPLICATION EXTENDS TO POSITIVE AND NEGATIVE INFINITY; AND C DISCONTINUITIES BETWEEN REPLICATIONS USUALLY REQUIRE WINDOWING C TO LESSEN THE EFFECTS. IN THE CASE OF A DFT, THE REPLICATION C BEGINS AT THE NEXT SAMPLE POINT AFTER THE END OF THE PREVIOUS C SAMPLE FUNCTION. THIS IS A SUBTLETY WHICH REQUIRES THE LAST C POINT OF AN OTHERWISE SYMMETRIC WINDOWING FUNCTION TO BE CUT C OFF. C C NOTE ON WINDOWING: ALL TIME-DOMAIN SAMPLE FUNCTIONS SHOULD BE C WINDOWED APPROPRIATELY BEFORE TRANSFORMING. BECAUSE THE C TRANSFORM ASSUMES PERIODIC EXTENSION OF THE FINITE SAMPLE C FUNCTION, DISCONTINUITIES MAY BE IMPLIED AT THE BEGINNING AND C END; THOSE DISCONTINUITIES ARE FOLDED BACK INTO THE FREQUENCY C DOMAIN AND PRODUCE SERIOUS SIDE LOBING ADJACENT TO EACH C DISCRETE FREQUENCY (BIN). A GOOD WINDOWING FUNCTION SLOWLY C DIMINISHES THE TIME-DOMAIN AMPLITUDE TOWARD THE END POINTS C (MINIMIZING REPLICATED DISCONTINUITIES); AND, IN SO DOING, IT C MINIMIZES SIDE LOBING AND FREQUENCY SMEARING WITHIN THE C FREQUENCY DOMAIN. C C NOTE ON THE KAISER-BESSEL WINDOW: ALTHOUGH THERE ARE MANY C WINDOWING FUNCTIONS AVAILABLE, DEMONSTRABLY, ONE OF THE BEST IS C THE KAISER-BESSEL WINDOW. IF THE RECOMMENDED ALPHA=3.0 IS C USED, IT HAS A MAXIMUM SIDE LOBE LEVEL OF -69 dB AND A 6-dB C BANDWIDTH OF 2.39 BINS. AN IN-DEPTH DISCUSSION OF THE C KAISER-BESSEL WINDOW (AND OTHER WINDOWS) MAY BE FOUND IN C "WINDOWS, HARMONIC ANALYSIS, AND THE DISCRETE FOURIER C TRANSFORM" BY FREDRIC J.HARRIS, UNDERSEA SURVEILLANCE C DEPARTMENT, SEPTEMBER 1976, NAVAL UNDERSEA CENTER, SAN DIEGO, C CALIFORNIA. THIS INFORMATION WAS SUPPLIED BY MIKE WEBRING, C USGS. THE KAISER-BESSEL-WINDOW EQUATION IS: C C w(n) = Io[pi*alpha*sqrt( 1-(n/(N/2))**2 )] / Io[pi*alpha] C C WHERE 0 <= |n| <= N/2 0 < alpha <= 3.5 C inf C Io[X] = sum[ ( ((X/2)**k)/k! )**2 ] C k=0 C C NOTE ON RESTORING WINDOWED DATA: A FULL LENGTH WINDOW (SUCH AS C THE KAISER-BESSEL) SIGNIFICANTLY REDUCES THE AMPLITUDE OF THE C ORIGINAL TIME DOMAIN SAMPLE FUNCTION ESPECIALLY NEAR THE END C POINTS. IF A SAMPLE FUNCTION HAS BEEN WINDOWED, FORWARD C TRANSFORMED, MODIFIED IN THE FREQUENCY DOMAIN, AND THEN C TRANSFORMED BACK TO THE TIME DOMAIN, OBVIOUSLY THE INVERSE C WINDOW MUST BE APPLIED TO RESTORE THE MODIFIED SAMPLE FUNCTION. C THIS PROCEDURE WORKS WELL EXCEPT FOR TWO MINOR PITFALLS: C 1) RESTORED FUNCTION VALUES NEAR THE ENDPOINTS MAY SUFFER C MACHINE PRECISION ERRORS; AND 2) CONTRIBUTIONS FROM WAVEFORMS C IN THE BEGINNING AND ENDING (31% FOR KAISER-BESSEL WITH C ALPHA=3.0) OF THE TIME DOMAIN FUNCTION WILL BE ATTENUATED BY C MORE THAN 6 dB. THIS MEANS THAT THE SPECTRUM WILL BE BIASED C TOWARD THE MIDDLE SECTION OF THE TIME DOMAIN FUNCTION AND THE C VERY LOW FREQUENCIES MAY BE COMPROMISED. HOWEVER, THESE C PITFALLS ARE MUCH LESS SEVERE THAN THOSE OF NOT USING A GOOD C WINDOW. C C NOTE ON OVERLAPPING WINDOWS: AS HAS BEEN EXPLAINED BY MIKE C WEBRING, THE CENTER BIAS PROBLEM MAY BE MINIMIZED BY WINDOWING C OVERLAPPING SECTIONS OF DATA AND THEN STACKING IN THE FREQUENCY C DOMAIN. OVERLAPPING AND STACKING OPTIMIZES THE EXPECTED (MEAN) C VALUE FOR THE SPECTRAL COMPONENTS. THE SEPARATION BETWEEN C CENTERS OF OVERLAPPING WINDOWS WILL VARY CONSIDERABLY DEPENDING C ON THE TYPE AND SPECIFICATIONS OF THE WINDOW USED. FOR C EXAMPLE, THE KAISER-BESSEL WINDOW WITH ALPHA=3.0 SHOULD HAVE C CENTER SEPARATIONS NO MORE THAN 38% OF EACH WINDOW LENGTH. C AMPLITUDES MAY INCREASE IN OVERLAPPING PARTS OF THE SAMPLE C FUNCTION. C C NOTE ON TESTS: MANY OF THE NOTES ABOVE WERE DERIVED FROM TESTS C DONE BY ROB BRACKEN, USGS. THE RESULTS OF THESE TESTS ARE IN C RING CORE NBF#6 9AUG96 - 16AUG96. C C C C INPUT ARGUMENTS: C NX - INTEGER*4. NUMBER OF POINTS IN ARRAY CX. NX MUST C BE AN INTEGER POWER OF 2. C SIGNI - REAL*8. DIRECTION OF TRANSFORM: C -1.D0 = FORWARD TRANSFORM (TIME TO FQ). C +1.D0 = INVERSE TRANSFORM (FQ TO TIME). C C INPUT/OUTPUT ARGUMENT: C CX - COMPLEX*16. ARRAY CONTAINING THE FUNCTION. CX MUST C HAVE DIMENSION NX (OR LARGER). ONLY ELEMENTS 1 C THROUGH NX WILL BE TRANSFORMED. C C SUBROUTINE FORK WRITTEN BY JON F. CLAERBOUT, STANFORD, 15FEB69. C ADAPTED TO USGS COMPUTER BY MIKE WEBRING, USGS, (DATE?) C CONVERSION TO DBLE CMPLX & NOTES BY ROB BRACKEN, USGS, 15AUG96. C HP-9000 SERIES 700/800 UNIX VERSION 1.0, 15AUG96. C C subroutine forkdc( nx, cx, signi ) C C DECLARATIONS C C INPUT ARGUMENTS integer*4 nx real*8 signi C C INPUT/OUTPUT ARGUMENT complex*16 cx(1) C C INTERNAL VARIABLES complex*16 carg,cw,ctemp real*8 spi,pi,sc,arg parameter(pi=3.1415 92653 58979 32384 62643) C C INITIALIZATIONS C spi=signi*pi sc=dsqrt( 1.d0/dflotj(nx) ) C C COMPUTATIONS C j=1 do 701 i=1,nx if(i.le.j) then ctemp=cx(j)*sc cx(j)=cx(i)*sc cx(i)=ctemp endif m=nx/2 202 if(j.gt.m) then j=j-m m=m/2 if(m.ge.1) goto 202 endif j=j+m 701 continue C ne=1 dowhile(ne.lt.nx) istep=2*ne do 703 m=1,ne arg=(spi*dflotj(m-1))/dflotj(ne) carg=dcmplx(0.d0,arg) cw=cdexp(carg) do 704 i=m,nx,istep k=i+ne ctemp=cw*cx(k) cx(k)=cx(i)-ctemp cx(i)=cx(i)+ctemp 704 continue 703 continue ne=istep enddo 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 DISCRETE 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 DISCRETE 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 DISCRETE 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 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 DISCRETE 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 DISCRETE 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 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 DISCRETE 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 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