SUBROUTINE DEMOD(X,N,F) C C PURPOSE--THIS SUBROUTINE PERFORMS A COMPLEX DEMODULATION C ON THE DATA IN THE INPUT VECTOR X C AT THE INPUT DEMODULATION FREQUENCY = F. C THE COMPLEX DEMODULATION CONSISTS OF THE FOLLOWING-- C 1) AN AMPLITUDE VERSUS TIME PLOT; C 2) A PHASE VERSUS TIME PLOT; C 3) AN UPDATED DEMODULATION FREQUENCY ESTIMATE C TO ASSIST THE ANALYST IN DETERMINING A C MORE APPROPRIATE FREQUENCY AT WHICH C TO DEMODULATE IN CASE THE SPECIFIED C INPUT DEMODULATION FREQUENCY F C DOES NOT FLATTEN SUFFICIENTLY THE C PHASE PLOT. C C THE ALLOWABLE RANGE OF THE INPUT DEMODULATION C FREQUENCY F IS 0.0 TO 0.5 (EXCLUSIVELY). C THE INPUT DEMODULATION FREQUENCY F IS MEASURED OF C IN UNITS OF CYCLES PER 'DATA POINT' OR, C MORE PRECISELY, IN CYCLES PER UNIT TIME WHERE C 'UNIT TIME' IS DEFINED AS THE C ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS. C C INPUT ARGUMENTS--X = THE SINGLE PRECISION VECTOR OF C (UNSORTED) OBSERVATIONS. C N = THE INTEGER NUMBER OF OBSERVATIONS C IN THE VECTOR X. C F = THE SINGLE PRECISION C DEMODULATION FREQUENCY. C F IS IN UNITS OF CYCLES PER DATA POINT. C F IS BETWEEN 0.0 AND 0.5 (EXCLUSIVELY). C OUTPUT--2 PAGES OF AUTOMATIC PRINTOUT-- C 1) AN AMPLITUDE PLOT; C 2) A PHASE PLOT; AND C 3) AN UPDATED DEMODULATION FREQUENCY ESTIMATE. C PRINTING--YES. C RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N C FOR THIS SUBROUTINE IS 5000. C --THE SAMPLE SIZE N MUST BE GREATER C THAN OR EQUAL TO 3. C --THE INPUT FREQUENCY F MUST BE C GREATER THAN OR EQUAL TO 2/(N-2). C --THE INPUT FREQUENCY F MUST BE C SMALLER THAN 0.5. C OTHER DATAPAC SUBROUTINES NEEDED--PLOTX. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, SIN, COS, ATAN. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--IN ORDER THAT THE RESULTS OF THE COMPLEX DEMODULATION C BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA C IN X SHOULD BE EQUI-SPACED IN TIME C (OR WHATEVER VARIABLE CORRESPONDS TO TIME). C --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED C TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME, C THEN THE DEMODULATION FREQUENCY F C WOULD BE IN UNITS OF HERTZ C (= CYCLES PER SECOND). C --A FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE C IN THE DATA OF INFINITE (= 1/(0.0)) C LENGTH OR PERIOD. C A FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE C IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS. C --IN EXAMINING THE AMPLITUDE AND PHASE PLOTS, C ATTENTION SHOULD BE PAID NOT ONLY TO THE C STRUCTURE OF THE PHASE PLOT C (NEAR-ZERO SLOPE VERSUS NON-ZERO SLOPE) C BUT ALSO TO THE RANGE C OF VALUES ON THE VERTICAL AXIS. C A PLOT WITH MUCH STRUCTURE BUT C WITH A SMALL RANGE ON THE VERTICAL AXIS C IS USUALLY MORE INDICATIVE OF A C DEFINITE CYCLIC COMPONENT AT THE C SPECIFIED INPUT DEMODULATION FREQUENCY, C THAN IS A PLOT WITH LESS STRUCTURE BUT C A WIDER RANGE ON THE VERTICAL AXIS. C --INTERNAL TO THIS SUBROUTINE, 2 MOVING C AVERAGES ARE APPLIED, EACH OF LENGTH 1/F. C HENCE THE AMPLITUDE AND PHASE PLOTS C HAVE N - 2/F VALUES C (RATHER THAN N VALUES) ALONG THE C HORIZONTAL (TIME) AXIS. C IN ORDER THAT THE AMPLITUDE AND PHASE C PLOTS BE NON-EMPTY, AN INPUT C REQUIREMENT ON F FOR THIS SUBROUTINE C IS THAT THE SAMPLE SIZE N C AND THE DEMODULATION FREQUENCY F C MUST BE SUCH THAT C N - 2/F BE GREATER THAN ZERO. C FURTHER, SINCE A PLOT WITH BUT C 1 POINT IS MEANINGLESS C AND OUGHT ALSO BE EXCLUDED, C THE REQUIREMENT IS EXTENDED C SO THAT N - 2/F MUST BE GREATER THAN 1. C REFERENCES--GRANGER AND HATANAKA, PAGES 170 TO 189, C ESPECIALLY PAGES 173, 177, AND 182. C WRITTEN BY--JAMES J. FILLIBEN C STATISTICAL ENGINEERING LABORATORY (205.03) C NATIONAL BUREAU OF STANDARDS C WASHINGTON, D. C. 20234 C PHONE--301-921-2315 C ORIGINAL VERSION--NOVEMBER 1972. C UPDATED --NOVEMBER 1975. C UPDATED --FEBRUARY 1976. C C--------------------------------------------------------------------- C DIMENSION X(1) DIMENSION Y1(5000),Y2(5000),Z(5000) COMMON /BLOCK2/ WS(15000) EQUIVALENCE (Y1(1),WS(1)),(Y2(1),WS(5001)),(Z(1),WS(10001)) DATA PI/3.141592653/ C IPR=6 ILOWER=3 IUPPER=5000 AN=N FMIN=2.0/(AN-2.0) C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.ILOWER.OR.N.GT.IUPPER)GOTO50 IF(F.LE.FMIN.OR.F.GE.0.5)GOTO60 HOLD=X(1) DO65I=2,N IF(X(I).NE.HOLD)GOTO90 65 CONTINUE WRITE(IPR, 9)HOLD RETURN 50 WRITE(IPR,17)ILOWER,IUPPER WRITE(IPR,47)N RETURN 60 WRITE(IPR,27)FMIN WRITE(IPR,46)F WRITE(IPR,28)FMIN,N RETURN 90 CONTINUE 9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE DEMOD SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 DEMOD SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,1H,,I6,16H) INTER 1VAL *****) 27 FORMAT(1H , 96H***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE 1 DEMOD SUBROUTINE IS OUTSIDE THE ALLOWABLE (,F10.8,6H,0.5) , 1 14HINTERVAL *****) 28 FORMAT(1H ,42H THE ABOVE LOWER LIMIT (,F10.8, 1 46H) = 2/(N-2) WHERE N = THE INPUT SAMPLE SIZE = ,I8) 46 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,E15.8,6H *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C C FORM THE COSINE AND SINE SERIES C DO100I=1,N AI=I Y1(I)=X(I)*COS(6.2831853*F*AI) Y2(I)=X(I)*SIN(6.2831853*F*AI) 100 CONTINUE C C DEFINE THE LENGTH OF THE 2 MOVING AVERAGES C LENMA1=1.0/F LENMA2=1.0/F ALEN1=LENMA1 ALEN2=LENMA2 IMAX1=N-LENMA1 IMAX2=IMAX1-LENMA2 C C FORM THE FIRST MOVING AVERAGE FOR THE COSINE SERIES C DO200I=1,IMAX1 ISTART=I+1 IEND=I+LENMA1-1 IENDP1=I+LENMA1 SUM=0.0 DO300J=ISTART,IEND SUM=SUM+Y1(J) 300 CONTINUE SUM=SUM+Y1(I)/2.0+Y1(IENDP1)/2.0 Z(I)=SUM/ALEN1 200 CONTINUE C C FORM THE SECOND MOVING AVERAGE FOR THE COSINE SERIES C DO400I=1,IMAX2 ISTART=I+1 IEND=I+LENMA2-1 IENDP1=I+LENMA2 SUM=0.0 DO500J=ISTART,IEND SUM=SUM+Z(J) 500 CONTINUE SUM=SUM+Z(I)/2.0+Z(IENDP1)/2.0 Y1(I)=SUM/ALEN2 400 CONTINUE C C FORM THE FIRST MOVING AVERAGE FOR THE SINE SERIES C DO800I=1,IMAX1 ISTART=I+1 IEND=I+LENMA1-1 IENDP1=I+LENMA1 SUM=0.0 DO900J=ISTART,IEND SUM=SUM+Y2(J) 900 CONTINUE SUM=SUM+Y2(I)/2.0+Y2(IENDP1)/2.0 Z(I)=SUM/ALEN1 800 CONTINUE C C FORM THE SECOND MOVING AVERAGE FOR THE SINE SERIES C DO1000I=1,IMAX2 ISTART=I+1 IEND=I+LENMA1-1 IENDP1=I+LENMA1 SUM=0.0 DO1100J=ISTART,IEND SUM=SUM+Z(J) 1100 CONTINUE SUM=SUM+Z(I)/2.0+Z(IENDP1)/2.0 Y2(I)=SUM/ALEN2 1000 CONTINUE C C C FORM THE AMPLITUDES AND PLOT THEM C DO1500I=1,IMAX2 Z(I)=2.0*SQRT(Y1(I)*Y1(I)+Y2(I)*Y2(I)) 1500 CONTINUE CALL PLOTX(Z,IMAX2) WRITE(IPR,2005)F C C COMPUTE THE DIFFERENCE BETWEEN THE MAX AND MIN AMPLITUDES AND WRITE IT OUT C ZMIN=Z(1) ZMAX=Z(1) DO1600I=1,IMAX2 IF(Z(I).LT.ZMIN)ZMIN=Z(I) IF(Z(I).GT.ZMAX)ZMAX=Z(I) 1600 CONTINUE RANGE=ZMAX-ZMIN WRITE(IPR,2010)ZMIN,ZMAX,RANGE C C FORM THE PHASES AND PLOT THEM C DO1700I=1,IMAX2 Z(I)=ATAN(Y1(I)/Y2(I)) 1700 CONTINUE CALL PLOTX(Z,IMAX2) WRITE(IPR,2020)F C C COMPUTE A NEW ESTIMATE FOR THE DEMODULATION FREQUENCY AND WRITE IT OUT C AIMAX2=IMAX2 IMAX2M=IMAX2-1 IFLAG=0 ZMIN=Z(1) ZMAX=Z(1) DO1800I=1,IMAX2M IP1=I+1 DEL=Z(IP1)-Z(I) IF(DEL.GT.2.5)IFLAG=IFLAG-1 IF(DEL.LT.-2.5)IFLAG=IFLAG+1 AIFLAG=IFLAG ZNEW=Z(IP1)+AIFLAG*PI IF(ZNEW.LT.ZMIN)ZMIN=ZNEW IF(ZNEW.GT.ZMAX)ZMAX=ZNEW 1800 CONTINUE RANGE=ZMAX-ZMIN SLOPER=RANGE/AIMAX2 SLOPEH=SLOPER/(2.0*PI) FEST=F+SLOPEH WRITE(IPR,2025)ZMIN,ZMAX,RANGE WRITE(IPR,2030)SLOPER,SLOPEH,FEST C 2005 FORMAT(1H ,30X, 48HAMPLITUDE PLOT FOR THE DEMODULATION FREQUENCY = 1 ,F8.6,21H CYCLES PER UNIT TIME) 2010 FORMAT(1H ,9X,20HMINIMUM AMPLITUDE = ,E15.8,5X,20HMAXIMUM AMPLITUD 1E = ,E15.8,5X,22HRANGE OF AMPLITUDES = ,E15.8) 2020 FORMAT(1H ,32X, 44HPHASE PLOT FOR THE DEMODULATION FREQUENCY = ,F8 1.6,21H CYCLES PER UNIT TIME) 2025 FORMAT(1H ,3X,16HMINIMUM PHASE = ,E15.8,11H RADIANS ,16HMAXIMUM 1PHASE = ,E15.8,11H RADIANS ,18HRANGE OF PHASES = ,E15.8,8H RADIA 1NS) 2030 FORMAT(1H ,8HSLOPE = ,E14.8,11H RADIANS = ,E14.6,52H CYCLES PER UN 1IT TIME EST. OF NEW DEMOD. FREQ. = ,E15.8,15H CYC./UNIT TIME) C RETURN END