C....*...1.........2.........3.........4.........5.........6.........7.*.......8 C FASTF 1/28/71 C C PURPOSE C COMPUTES ALL FOURIER COEFFICIENTS AND PERIODOGRAM ORDINATES FOR A C REAL TIME SERIES OF LENGTH N. ON RETURN THE ARRANGEMENT IS: C C=A(0),A(1),B(1),...,A((N-2)/2),B((N-2)/2),A(N/2) N COEFFICIENTS C P=P(W0),P(W1),P(W2),...,P(WN/2) NP ORDINATES C C DPER is probably a better choice. C C USAGE C CALL FASTF(X,N,C,P,NP) C C SUBROUTINES CALLED C FFT, REALTR C C ARGUMENTS C X - INPUT REAL TIME SERIES. C N - LENGTH OF X. SEE RESTRICTIONS ON N BELOW. C C - VECTOR OF COMPUTED FOURIER COEFFICIENTS, DIMENSION C .GT. N+2 C P - VECTOR OF COMPUTED PERIODOGRAM ORDINATES, DIMENSION P .GT. N+2 C NP- COMPUTED NUMBER OF PERIODOGRAM ORDINATES C C RESTRICTIONS ON N C FACTOR N INTO C N = N(1)*N(2)*N(3)*...*N(K) C WHERE EACH N(I) IS A PRIME NUMBER. C 1. N MUST BE EVEN AND LARGER THAN 2. C 2. THE LARGEST N(I) MUST BE LESS THAN OR EQUAL TO 23. C 3. K MUST BE LESS THAN OR EQUAL TO 11. C 4. THE PRODUCT OF THOSE PRIME NUMBERS WHICH OCCUR AN ODD NUMBER C OF TIMES IN THE FACTORIZATION OF N MUST BE LESS THAN OR EQUAL C TO 210. C C EXAMPLE C 10200 = 2*2*2*3*5*5*17 C 1. 10200 IS EVEN AND LARGER THAN 2. C 2. 17 IS LESS THAN 23. C 3. K=7 IS LESS THAN 11. C 4. 2*3*17=102 IS LESS THAN 210. C C REMARK C IF THE FOURIER COEFFICIENTS ARE NOT NEEDED THE USAGE C CALL FASTF(X,N,P,P,NP) C IS PERMISSIBLE. IF X CAN BE DESTROYED THE USAGE C CALL FASTF(X,N,X,X,NP) C IS PERMISSIBLE. X WILL CONTAIN THE PERIODOGRAM ORDINATES ON C RETURN. DIMENSION X .GT. N+2 WITH THIS USAGE. C C REFERENCE C R. C. SINGLETON. AN ALGORITHM FOR COMPUTING THE MIXED RADIX FAST C FOURIER TRANSFORM. IEEE TRANSACTIONS ON AUDIO AND ELECTROACOUSTICS C AU-17. JUNE 1969. C SUBROUTINE FASTF(X,N,C,P,NP) implicit real*8 (a-h,o-z) save REAL*8 X(1),C(2),P(1) DO 10 I=1,N 10 C(I)=X(I) M=N/2 NP=M+1 ISN=2 CALL FFT(C,C(2),M,M,M,ISN) CALL REALTR(C,C(2),M,ISN) A=1.d0/N B=N*.5d0 C(1)=A*C(1)*.5d0 P(1)=B*C(1)*C(1)*4.d0 DO 20 I=2,M K=2*I-2 C(K)=A*C(K+1) C(K+1)=A*C(K+2) P(I)=B*C(K)**2+B*C(K+1)**2 20 CONTINUE C(N)=A*C(N+1)*.5d0 P(NP)=B*C(N)*C(N)*4.d0 RETURN END