demod(3f) - [M_datapac:STATISTICS] perform a complex demodulation
SUBROUTINE DEMOD(X,N,F)
demod(3f) performs a complex demodulation on the data in the input
vector x at the input demodulation frequency = f.
the complex demodulation consists of the following--
1. an amplitude versus time plot;
2. a phase versus time plot;
3. an updated demodulation frequency estimate
to assist the analyst in determining a
more appropriate frequency at which
to demodulate in case the specified
input demodulation frequency f
does not flatten sufficiently the
phase plot.
the allowable range of the input demodulation frequency f is 0.0 to
0.5 (exclusively).
the input demodulation frequency f is measured of in units of
cycles per 'data point' or, more precisely, in cycles per unit time
where 'unit time' is defined as the elapsed time between adjacent
observations.
X description of parameter
Y description of parameter
Sample program:
program demo_demod
use M_datapac, only : demod
implicit none
! call demod(x,y)
end program demo_demod
Results:
The original DATAPAC library was written by James Filliben of the Statistical
Engineering Division, National Institute of Standards and Technology.
John Urban, 2022.05.31
CC0-1.0
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=wp), | dimension(:) | :: | X | |||
| integer | :: | N | ||||
| real(kind=wp) | :: | F |
| Type | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|
| real | :: | WS(15000) |
SUBROUTINE DEMOD(X,N,F) REAL(kind=wp) :: ai, aiflag, aimax2, alen1, alen2, an, del, F, fest, fmin, hold, pi, range, slopeh, sloper, sum, WS, X, Y1, Y2 REAL(kind=wp) :: Z, zmax, zmin, znew INTEGER :: i, iend, iendp1, iflag, ilower, imax1, imax2, imax2m, ip1, istart, iupper, j, lenma1, lenma2, N ! ! INPUT ARGUMENTS--X = THE VECTOR OF ! (UNSORTED) OBSERVATIONS. ! N = THE INTEGER NUMBER OF OBSERVATIONS ! IN THE VECTOR X. ! F = THE ! DEMODULATION FREQUENCY. ! F IS IN UNITS OF CYCLES PER DATA POINT. ! F IS BETWEEN 0.0 AND 0.5 (EXCLUSIVELY). ! OUTPUT--2 pages OF AUTOMATIC PRINTOUT-- ! 1) AN AMPLITUDE PLOT; ! 2) A PHASE PLOT; AND ! 3) AN UPDATED DEMODULATION FREQUENCY ESTIMATE. ! PRINTING--YES. ! RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N ! FOR THIS SUBROUTINE IS 5000. ! --THE SAMPLE SIZE N MUST BE GREATER ! THAN OR EQUAL TO 3. ! --THE INPUT FREQUENCY F MUST BE ! GREATER THAN OR EQUAL TO 2/(N-2). ! --THE INPUT FREQUENCY F MUST BE ! SMALLER THAN 0.5. ! OTHER DATAPAC SUBROUTINES NEEDED--PLOTX. ! FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, SIN, COS, ATAN. ! MODE OF INTERNAL OPERATIONS--. ! COMMENT--IN ORDER THAT THE RESULTS OF THE COMPLEX DEMODULATION ! BE VALID AND PROPERLY INTERPRETED, THE INPUT DATA ! IN X SHOULD BE EQUI-SPACED IN TIME ! (OR WHATEVER VARIABLE CORRESPONDS TO TIME). ! --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED ! TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME, ! THEN THE DEMODULATION FREQUENCY F ! WOULD BE IN UNITS OF HERTZ ! (= CYCLES PER SECOND). ! --A FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE ! IN THE DATA OF INFINITE (= 1/(0.0)) ! LENGTH OR PERIOD. ! A FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE ! IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS. ! --IN EXAMINING THE AMPLITUDE AND PHASE PLOTS, ! ATTENTION SHOULD BE PAID NOT ONLY TO THE ! STRUCTURE OF THE PHASE PLOT ! (NEAR-ZERO SLOPE VERSUS NON-ZERO SLOPE) ! BUT ALSO TO THE RANGE ! OF VALUES ON THE VERTICAL AXIS. ! A PLOT WITH MUCH STRUCTURE BUT ! WITH A SMALL RANGE ON THE VERTICAL AXIS ! IS USUALLY MORE INDICATIVE OF A ! DEFINITE CYCLIC COMPONENT AT THE ! SPECIFIED INPUT DEMODULATION FREQUENCY, ! THAN IS A PLOT WITH LESS STRUCTURE BUT ! A WIDER RANGE ON THE VERTICAL AXIS. ! --INTERNAL TO THIS SUBROUTINE, 2 MOVING ! AVERAGES ARE APPLIED, EACH OF LENGTH 1/F. ! HENCE THE AMPLITUDE AND PHASE PLOTS ! HAVE N - 2/F VALUES ! (RATHER THAN N VALUES) ALONG THE ! HORIZONTAL (TIME) AXIS. ! IN ORDER THAT THE AMPLITUDE AND PHASE ! PLOTS BE NON-EMPTY, AN INPUT ! REQUIREMENT ON F FOR THIS SUBROUTINE ! IS THAT THE SAMPLE SIZE N ! AND THE DEMODULATION FREQUENCY F ! MUST BE SUCH THAT ! N - 2/F BE GREATER THAN ZERO. ! FURTHER, SINCE A PLOT WITH BUT ! 1 POINT IS MEANINGLESS ! AND OUGHT ALSO BE EXCLUDED, ! THE REQUIREMENT IS EXTENDED ! SO THAT N - 2/F MUST BE GREATER THAN 1. ! ORIGINAL VERSION--NOVEMBER 1972. ! UPDATED --NOVEMBER 1975. ! UPDATED --FEBRUARY 1976. ! !--------------------------------------------------------------------- ! DIMENSION X(:) DIMENSION Y1(5000) , Y2(5000) , Z(5000) COMMON /BLOCK2_real32/ WS(15000) EQUIVALENCE (Y1(1),WS(1)) EQUIVALENCE (Y2(1),WS(5001)) EQUIVALENCE (Z(1),WS(10001)) DATA pi/3.141592653_wp/ ! ilower = 3 iupper = 5000 an = N fmin = 2.0_wp/(an-2.0_wp) ! ! CHECK THE INPUT ARGUMENTS FOR ERRORS ! IF ( N<ilower .OR. N>iupper ) THEN WRITE (G_IO,99001) ilower , iupper 99001 FORMAT (' ', & &'***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE DEMOD SUBROU& &TINE IS OUTSIDE THE ALLOWABLE (',I0,',',I0,') INTERVAL *****') WRITE (G_IO,99002) N 99002 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',I0,' *****') RETURN ELSE IF ( F<=fmin .OR. F>=0.5_wp ) THEN WRITE (G_IO,99003) fmin 99003 FORMAT (' ', & &'***** FATAL ERROR--THE THIRD INPUT ARGUMENT TO THE DEMOD SUBROU& &TINE IS OUTSIDE THE ALLOWABLE (',F10.8,',0.5) ','INTERVAL *****') WRITE (G_IO,99004) F 99004 FORMAT (' ','***** THE VALUE OF THE ARGUMENT IS ',E15.8, & & ' *****') WRITE (G_IO,99005) fmin , N 99005 FORMAT (' ',' THE ABOVE LOWER LIMIT (', & & F10.8, & & ') = 2/(N-2) WHERE N = THE INPUT SAMPLE SIZE = ',I0) RETURN ELSE hold = X(1) DO i = 2 , N IF ( X(i)/=hold ) GOTO 50 ENDDO WRITE (G_IO,99006) hold 99006 FORMAT (' ', & &'***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUMENT (A VECTOR) & &TO THE DEMOD SUBROUTINE HAS ALL ELEMENTS = ',E15.8,' *****') RETURN ENDIF ! !-----START POINT----------------------------------------------------- ! ! FORM THE COSINE AND SINE SERIES ! 50 DO i = 1 , N ai = i Y1(i) = X(i)*COS(6.2831853_wp*F*ai) Y2(i) = X(i)*SIN(6.2831853_wp*F*ai) ENDDO ! ! DEFINE THE LENGTH OF THE 2 MOVING AVERAGES ! lenma1 = 1.0_wp/F lenma2 = 1.0_wp/F alen1 = lenma1 alen2 = lenma2 imax1 = N - lenma1 imax2 = imax1 - lenma2 ! ! FORM THE FIRST MOVING AVERAGE FOR THE COSINE SERIES ! DO i = 1 , imax1 istart = i + 1 iend = i + lenma1 - 1 iendp1 = i + lenma1 sum = 0.0_wp DO j = istart , iend sum = sum + Y1(j) ENDDO sum = sum + Y1(i)/2.0_wp + Y1(iendp1)/2.0_wp Z(i) = sum/alen1 ENDDO ! ! FORM THE SECOND MOVING AVERAGE FOR THE COSINE SERIES ! DO i = 1 , imax2 istart = i + 1 iend = i + lenma2 - 1 iendp1 = i + lenma2 sum = 0.0_wp DO j = istart , iend sum = sum + Z(j) ENDDO sum = sum + Z(i)/2.0_wp + Z(iendp1)/2.0_wp Y1(i) = sum/alen2 ENDDO ! ! FORM THE FIRST MOVING AVERAGE FOR THE SINE SERIES ! DO i = 1 , imax1 istart = i + 1 iend = i + lenma1 - 1 iendp1 = i + lenma1 sum = 0.0_wp DO j = istart , iend sum = sum + Y2(j) ENDDO sum = sum + Y2(i)/2.0_wp + Y2(iendp1)/2.0_wp Z(i) = sum/alen1 ENDDO ! ! FORM THE SECOND MOVING AVERAGE FOR THE SINE SERIES ! DO i = 1 , imax2 istart = i + 1 iend = i + lenma1 - 1 iendp1 = i + lenma1 sum = 0.0_wp DO j = istart , iend sum = sum + Z(j) ENDDO sum = sum + Z(i)/2.0_wp + Z(iendp1)/2.0_wp Y2(i) = sum/alen2 ENDDO ! ! ! FORM THE AMPLITUDES AND PLOT THEM ! DO i = 1 , imax2 Z(i) = 2.0_wp*SQRT(Y1(i)*Y1(i)+Y2(i)*Y2(i)) ENDDO CALL PLOTX(Z,imax2) WRITE (G_IO,99007) F ! 99007 FORMAT (' ',30X, & & 'AMPLITUDE PLOT FOR THE DEMODULATION FREQUENCY = ', & & F8.6,' CYCLES PER UNIT TIME') ! ! COMPUTE THE DIFFERENCE BETWEEN THE MAX AND MIN AMPLITUDES AND WRITE IT OUT ! zmin = Z(1) zmax = Z(1) DO i = 1 , imax2 IF ( Z(i)<zmin ) zmin = Z(i) IF ( Z(i)>zmax ) zmax = Z(i) ENDDO range = zmax - zmin WRITE (G_IO,99008) zmin , zmax , range 99008 FORMAT (' ',9X,'MINIMUM AMPLITUDE = ',E15.8,5X, & & 'MAXIMUM AMPLITUDE = ',E15.8,5X, & & 'RANGE OF AMPLITUDES = ',E15.8) ! ! FORM THE PHASES AND PLOT THEM ! DO i = 1 , imax2 Z(i) = ATAN(Y1(i)/Y2(i)) ENDDO CALL PLOTX(Z,imax2) WRITE (G_IO,99009) F 99009 FORMAT (' ',32X,'PHASE PLOT FOR THE DEMODULATION FREQUENCY = ',& & F8.6,' CYCLES PER UNIT TIME') ! ! COMPUTE A NEW ESTIMATE FOR THE DEMODULATION FREQUENCY AND WRITE IT OUT ! aimax2 = imax2 imax2m = imax2 - 1 iflag = 0 zmin = Z(1) zmax = Z(1) DO i = 1 , imax2m ip1 = i + 1 del = Z(ip1) - Z(i) IF ( del>2.5_wp ) iflag = iflag - 1 IF ( del<-2.5_wp ) iflag = iflag + 1 aiflag = iflag znew = Z(ip1) + aiflag*pi IF ( znew<zmin ) zmin = znew IF ( znew>zmax ) zmax = znew ENDDO range = zmax - zmin sloper = range/aimax2 slopeh = sloper/(2.0_wp*pi) fest = F + slopeh WRITE (G_IO,99010) zmin , zmax , range 99010 FORMAT (' ',3X,'MINIMUM PHASE = ',E15.8,' RADIANS ', & & 'MAXIMUM PHASE = ',E15.8,' RADIANS ', & & 'RANGE OF PHASES = ',E15.8,' RADIANS') WRITE (G_IO,99011) sloper , slopeh , fest 99011 FORMAT (' ','SLOPE = ',E14.8,' RADIANS = ',E14.6, & & ' CYCLES PER UNIT TIME EST. OF NEW DEMOD. FREQ. = ',& & E15.8,' CYC./UNIT TIME') ENDIF ! END SUBROUTINE DEMOD