SUBROUTINE FOURIE(X,N) C C PURPOSE--THIS SUBROUTINE PERFORMS A FOURIER ANALYSIS C OF THE DATA IN THE INPUT VECTOR X. C THE ANALYSIS CONSISTS OF THE FOLLOWING-- C 1) COMPUTING (AND PRINTING) C (FOR EACH OF THE HARMONIC FREQUENCIES C 1/N, 2/N, 3/N, ..., 1/2) C THE CORRESPONDING FOURIER COEFICIENTS, C THE AMPLITUDE, THE PHASE, C THE CONTRIBUTION TO THE TOTAL VARIANCE, C AND THE RELATIVE CONTRIBUTION TO THE TOTAL C VARIANCE. C 2) PLOTTING OUT A FOURIER LINE SPECTRUM = C THE PERIODOGRAM = THE PLOT OF RELATIVE C CONTRIBUTION TO TOTAL VARIANCE C (AT EACH FOURIER FREQUENCY) VERSUS C THE FOURIER FREQUENCY. C C IN ORDER THAT THE RESULTS OF THE FOURIER ANALYSIS 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 C THE HORIZONTAL AXIS OF THE SPECTRA PRODUCED C BY THIS SUBROUTINE IS FREQUENCY. C THIS FREQUENCY IS MEASURED IN UNITS OF C CYCLES PER 'DATA POINT' OR, MORE PRECISELY, IN C CYCLES PER UNIT TIME WHERE C 'UNIT TIME' IS DEFINED AS THE C ELAPSED TIME BETWEEN ADJACENT OBSERVATIONS. C THE RANGE OF THE FREQUENCY AXIS IS 0.0 TO 0.5. 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 OUTPUT--2 TO 10 PAGES (DEPENDING ON C THE INPUT SAMPLE SIZE) OF C AUTOMATIC PRINTOUT-- C 1) A LISTING OF THE AMPLITUDE, C PHASE, CONTRIBUTION TO THE C TOTAL VARIANCE, AND RELATIVE C CONTRIBUTION TO THE TOTAL C VARIANCE FOR EACH OF THE C FOURIER FREQUENCIES C (1/N, 2/N, 3/N, ..., 1/2). C THIS LISTING MAY TAKE AS LITTLE AS 1 C PAGE OR AS MANY AS N/100 PAGES C (THE EXACT NUMBER DEPENDING ON C THE INPUT SAMPLE SIZE N). C THIS LISTING IS TERMINATED C AFTER AT MOST 8 COMPUTER PAGES. C IF MORE PAGES ARE DESIRED, C CHANGE THE VALUE OF THE C VARIABLE MAXPAG C WITHIN THIS SUBROUTINE C FROM 8 TO WHATEVER DESIRED. C 2) A PLOT OF THE RELATIVE C CONTRIBUTION TO THE C TOTAL VARIANCE VERSUS FREQUENCY. C PRINTING--YES. C RESTRICTIONS--THE MAXIMUM ALLOWABLE VALUE OF N C FOR THIS SUBROUTINE IS 15000. C --THE SAMPLE SIZE N MUST BE GREATER C THAN OR EQUAL TO 3. C OTHER DATAPAC SUBROUTINES NEEDED--PLOTSP AND CHSPPF. C FORTRAN LIBRARY SUBROUTINES NEEDED--SQRT, SIN, COS, ATAN. C MODE OF INTERNAL OPERATIONS--SINGLE PRECISION. C LANGUAGE--ANSI FORTRAN. C COMMENT--FOURIER ANALYSIS DIFFERS FROM SPECTRAL ANALYSIS C (AS, FOR EXAMPLE, PRODUCED BY THE DATAPAC C TIMESE SUBROUTINE) IN THAT A C FOURIER ANALYSIS DOES NO SMOOTHING ON C THE SPECTRAL ESTIMATES WHEREAS A SPECRRAL C ANALYSIS DOES SMOOTH THE SPECTRAL ESTIMATES. C THE NET RESULT IS THAT THE SPECTRAL C ESTIMATES OBTAINED FROM A FOURIER C ANALYSIS ARE ALMOST ALWAYS MORE C VARIABLE THAN THOSE OBTAINED IN A C SPECTRAL ANALYSIS. C THE PRACTICAL CONCLUSION IS THAT C WHEN THE DATA ANALYST HAS A CHOICE C OF WHETHER TO PERFORM A FOURIER C ANALYSIS OR A SPECTRAL ANALYSIS, C THE SPECTRAL ANALYSIS SHOULD C ALMOST ALWAYS BE PREFERRED. C --THE MAXIMUM NUMBER OF FOURIER FREQUENCIES C FOR WHICH THE FOURIER COEFFICIENTS IS C COMPUTED (AND LISTED) IS N/2 WHERE N IS C THE SAMPLE SIZE (LENGTH OF THE C DATA RECORD IN THE VECTOR X). C THIS RULE IS OVERRIDDEN C (FOR LISTING PURPOSES ONLY) C IN LARGE DATA SETS AND IS REPLACED C BY THE RULE THAT THE MAXIMUM C NUMBER OF LAGS LISTED = 800 C (WHICH CORRESPONDS TO AN C 8-PAGE LISTING OF FOURIER COEFFICIENTS. C IF MORE PAGES ARE DESIRED, C CHANGE THE VALUE OF THE C VARIABLE MAXPAG C WITHIN THIS SUBROUTINE C FROM 8 TO WHATEVER DESIRED. C --IF THE INPUT OBSERVATIONS IN X ARE CONSIDERED C TO HAVE BEEN COLLECTED 1 SECOND APART IN TIME, C THEN THE FREQUENCY AXIS OF THE RESULTING C SPECTRA WOULD BE IN UNITS OF HERTZ C (= CYCLES PER SECOND). C --THE FREQUENCY OF 0.0 CORRESPONDS TO A CYCLE C IN THE DATA OF INFINITE (= 1/(0.0)) C LENGTH OR PERIOD. C THE FREQUENCY OF 0.5 CORRESPONDS TO A CYCLE C IN THE DATA OF LENGTH = 1/(0.5) = 2 DATA POINTS. C --ANY EQUI-SPACED FOURIER ANALYSIS IS C INTRINSICALLY LIMITED TO DETECTING FREQUENCIES C NO LARGER THAN 0.5 CYCLES PER DATA POINT; C THIS CORRESPONDS TO THE FACT THAT THE C SMALLEST DETECTABLE CYCLE IN THE DATA C IS 2 DATA POINTS PER CYCLE. C REFERENCES--JENKINS AND WATTS, ESPECIALLY PAGE 290. 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 CHARACTER*4 ALPERC DIMENSION X(1) DIMENSION A(7500),B(7500) COMMON /BLOCK2/ WS(15000) EQUIVALENCE (A(1),WS(1)),(B(1),WS(7501)) DATA PI/3.14159265358979/ DATA ALPERC/'%'/ C IPR=6 ILOWER=3 IUPPER=15000 MAXPAG=8 C C CHECK THE INPUT ARGUMENTS FOR ERRORS C IF(N.LT.ILOWER.OR.N.GT.IUPPER)GOTO50 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 90 CONTINUE 9 FORMAT(1H ,109H***** NON-FATAL DIAGNOSTIC--THE FIRST INPUT ARGUME 1NT (A VECTOR) TO THE FOURIE SUBROUTINE HAS ALL ELEMENTS = ,E15.8,6 1H *****) 17 FORMAT(1H , 96H***** FATAL ERROR--THE SECOND INPUT ARGUMENT TO THE 1 FOURIE SUBROUTINE IS OUTSIDE THE ALLOWABLE (,I6,1H,,I6,16H) INTER 1VAL *****) 47 FORMAT(1H , 35H***** THE VALUE OF THE ARGUMENT IS ,I8 ,6H *****) C C-----START POINT----------------------------------------------------- C AN=N C C DETERMINE IF N IS ODD OR EVEN C IEVODD=N-2*(N/2) DEL=(AN+1.0)/2.0 IF(IEVODD.EQ.0)DEL=(AN+2.0)/2.0 C C COMPUTE THE SAMPLE MEAN C SUM=0.0 DO100I=1,N SUM=SUM+X(I) 100 CONTINUE XBAR=SUM/AN C C COMPUTE THE BIASED SAMPLE VARIANCE C SUM=0.0 DO200I=1,N SUM=SUM+(X(I)-XBAR)**2 200 CONTINUE VBIAS=SUM/AN C C COMPUTE THE FOURIER COSINE AND SINE COEFFICIENTS--THEY ARE PLACED C IN VECTORS A AND B, RESPECTIVELY. C NHALF=N/2 DO400I=1,NHALF AI=I SUMA=0.0 SUMB=0.0 DO500J=1,N T=J SUMA=SUMA+X(J)*COS(2.0*PI*(AI/AN)*(T-DEL)) SUMB=SUMB+X(J)*SIN(2.0*PI*(AI/AN)*(T-DEL)) 500 CONTINUE A(I)=SUMA/AN B(I)=SUMB/AN 400 CONTINUE C C WRITE OUT THE SAMPLE SIZE, THE SAMPLE MEAN, C AND THE (BIASED) SAMPLE VARIANCE. C WRITE(IPR,998) WRITE(IPR,801) WRITE(IPR,999) WRITE(IPR,999) WRITE(IPR,805)N WRITE(IPR,810)XBAR WRITE(IPR,815)VBIAS WRITE(IPR,999) C C COMPUTE THE HARMONIC CONTRIBUTION C AT EACH OF THE FOURIER FREQUENCIES. C THE FUNDAMENTAL FOURIER FREQUENCY C IS 1/N CYCLES PER DATA POINT C (WHERE N = THE INPUT SAMPLE SIZE). C THE OTHER FOURIER FREQUENCIES C ARE MULTIPLES OR HARMONICS C (2/N, 3/N, 4/N, ...1/2) OF THE FUNDAMENTAL. C COMPUTE AMPLITUDES, PHASES, AND C CONTRIBUTIONS TO THE VARIANCE AT EACH C OF THE FOURIER FREQUENCIES. C COMPUTE THE PERCENTAGE CONTRIBUTION C TO THE TOTAL VARIANCE AT EACH C OF THE FOURIER FREQUENCIES. C NOTE--TO SAVE STORAGE, ALSO COPY C THE PERCENTAGE CONTRIBUTIONS TO THE VARIANCE) C (WHICH WILL LATER BE PLOTTED OUT LIKE A SPECTRUM) C INTO THE VECTOR A; THIS WILL DESTROY C THE PREVIOUS CONTENTS OF THE VECTOR A. C WRITE OUT ALL OF THE ABOVE. C NNPAGE=50 I=0 DO600IPAGE=1,MAXPAG WRITE(IPR,998) WRITE(IPR,820) WRITE(IPR,821) WRITE(IPR,822) DO700J=1,NNPAGE I=I+1 AI=I FFREQ =AI/AN PERIOD=1.0/FFREQ ANGRAD =(AI/AN)*2.0*PI ANGDEG =(AI/AN)*360.0 AMP =SQRT(A(I)*A(I)+B(I)*B(I)) PHASE1 =ATAN(-B(I)/A(I)) PHASE2 =PHASE1 *360.0/(2.0*PI) CONMSQ =2.0*AMP *AMP IF(I.EQ.NHALF.AND.IEVODD.EQ.0)CONMSQ=CONMSQ/2.0 PERCON =100.0*(CONMSQ /VBIAS) WRITE(IPR,825)I,FFREQ,PERIOD,A(I),B(I),AMP,PHASE1,PHASE2, 1CONMSQ,PERCON,ALPERC A(I)=PERCON IF(I.GE.NHALF)GOTO750 ISKIP=I-10*(I/10) IF(ISKIP.EQ.0)WRITE(IPR,999) 700 CONTINUE 600 CONTINUE 750 CONTINUE C C PLOT OUT THE PERCENTAGE CONTRIBUTIONS C TO THE TOTAL VARIANCE AT C EACH OF THE FOURIER FREQUENCIES C (1/N, 2/N, 3/N, ..., 1/2). C THIS WILL CORRESPOND TO A SPECTRAL C PLOT IN SPECTRAL ANALYSIS. C CALL PLOTSP(A,NHALF,0) WRITE(IPR,855) C 801 FORMAT(1H ,44X,16HFOURIER ANALYSIS) 805 FORMAT(1H ,40X,41HTHE SAMPLE SIZE N = ,I8) 810 FORMAT(1H ,40X,41HTHE SAMPLE MEAN = ,F20.8) 815 FORMAT(1H ,40X,41HTHE SAMPLE VARIANCE (WITH DIVISOR N-1) = ,F20.8) 820 FORMAT(1H ,40H I FOURIER PERIOD FOURIER , 1 30H FOURIER AMPLITUDE , 1 46H PHASE PHASE VARIANCE , 1 10H RELATIVE) 821 FORMAT(1H ,44H FREQUENCY COEFFICIENT , 1 11HCOEFFICIENT , 1 59H RADIANS DEGREES COMPONENT, 1 12H VARIANCE) 822 FORMAT(1H ,43H (CYCLES/POINT) A(I) , 1 14H B(I) , 1 50H , 1 22H COMPONENT (%)) 825 FORMAT(1H ,I6,2X,F8.6,1X,F8.2,6(1X,E14.7),2X,F6.2,A1) 855 FORMAT(1H ,40X,56HPERIODOGRAM = FOURIER LINE SPECTRUM OF THE ORIGI 1NAL DATA) 998 FORMAT(1H1) 999 FORMAT(1H ) C RETURN END