SUBROUTINE ENVEL (DATA,NP) 

 C THIS ROUTINE RECEIVES A TIME SERIES OF LENGTH NP 



C IN ARRAY DATA, AND RETURNS ENVELOPE ESTIMATES BASED 



C ON A HILBERT TRANSFORM METHOD AS DESCRIBED IN KANASEWICH 



C (1981), P. 362-368. 



C C.G. FOX, ADVANCED TECHNOLOGY STAFF, NAVOCEANO-5/12/83 



DIMENSION OATA(l) 

 COMPLEX XDATA(2048) 

 NSTART=1 

 NPSEG=NP 

 10 IF(NPSEG.LE.O) RETURN 

 C PERFORM ENVELOPE CALCULATIONS IN SEGMENTS OF 2048 POINTS 



IF(NPSEG.GT.2048) THEN 

 NPSEG=NPSEG-2048 

 LENGTH =2048 

 GO TO 100 

 ELSE 



LENGTH=NPSEG 

 NPSEG=0 

 END IF 

 C TRANSFER INPUT DATA TO REAL PORTION OF COMPLEX ARRAY XDATA 



100 DO 110 I=1,LENGTH 



110 XDATA(I)=(1.0,0.0)*DATA(NSTART+(I-1)) 

 C ZERO FILL ARRAY IF LESS THAN 2048 POINTS 



IF ( LENGTH. EQ. 2048) GO TO 200 

 DO 120 I=LENGTH+1,2048 

 120 XDATA(I)=(0.0,0.0) 

 C FOURIER TRANSFORM TO FREQUENCY DOMAIN USING FFT 



200 CALL NLOGNdl, XDATA, -1.0) 

 C COMPUTE HILBERT TRANSFORM IN THE FREQUENCY DOMAIN. THE 



C TRANSFORM IS COMPUTED AS I*SGN (OMEGA )*F (OMEGA), WHERE 



C OMEGA=FREQUENCY F(OMEGA)=FOURIER TRANSFORM 



C SGN(X)=1,(X>0),=0,(X=0),=-1(X<0) 



C THIS OPERATION INDUCES A NINETY DEGREE PHASE SHIFT ON 



C COMPLEX PLANE. THE CALCULATION IS DONE BY ZEROING X(l), 



C (OMEGA=0), TRANSFERRING REAL TO IMAGINARY, MINUS IMAG TO REAL 



C FOR X(2)-X(N/2),(OMEGA>0), AND TRANSFERRING MINUS REAL TO 



C IMAGINARY, AND IMAGINARY TO REAL FOR X (N/2+1) TO X(N), 



C (OMEGA<0) 



XDATA(1)=(0.0,0.0) 

 DO 250 1=2,1024 

 TEMPR=REAL(XDATA(I)) 

 TEMPI=AIMAG(XDATA(I)) 

 ■250 XDATA(I)=CMPLX(-TEMPI,TEMPR) 

 DO 260 1=1025,2048 

 TEMPR=REAL(XDATA(I)) 

 TEMPI=AIMAG(XDATA(I)) 

 260 XDATA(I)=CMPLX(TEMPI,-TEMPR) 

 C INVERSE TRANSFORM SHIFTED DATA 



CALL NLOGNdl, XDATA, +1.0) 

 C CALCULATE ENVELOPE AFTER EQUATION 21.3-1 OF KANASEWICH 



C DUE TO SYMMETRIES IN THE TRANSFORM, IMAGINARY PORTION 



C OF XDATA IS NEAR ZERO AND CAN BE IGNORED 



172 



