C----------------------------------------------------------------------------- C SUBROUTINE FHTFOR TAKES INPUT F() AND RETURNS THE DHT IN THE SAME F() C LENGTH OF F() IS 2**P C----------------------------------------------------------------------------- C THIS IS THE CORRECTED VERSION OF FHTFOR.FOR C RN BRACEWELL, THE HARTLEY TRANSFORM, Oxford Univ. Press (1986) C pp. 140-144 C----------------------------------------------------------------------------- C P >= 2 C----------------------------------------------------------------------------- SUBROUTINE FHTFR(P,F) DIMENSION F(0:255) DIMENSION S9(256), T9(256) DIMENSION V9(0:10) INTEGER P, D9, E9, Q9, SX9, U9, A9(0:64), M9(0:10), C9(0:10) INTEGER POLD DATA POLD /-99/ PI = 3.1415926536 IF (P.EQ.1) THEN J = F(0) + F(1) F(1) = F(0) - F(1) F(0) = J GOTO 9636 ! RETURN ENDIF IF(P .EQ. POLD) GO TO 9400 N9 = 2**(P-2) N = 4*N9 C9(5) = N-1 C9(6) = P-1 C GET POWERS OF 2 I = 1 M9(0) = 1 M9(1) = 2 9204 M9(I+1) = M9(I) + M9(I) I = I + 1 IF(I.LE.P) GO TO 9204 C TRIGONOMETRIC FUNCTION TABLES IF(N.EQ.4) GO TO 9400 IF(N.EQ.8) THEN S9(2) = 1. S9(1) = SIN(PI/4.) GO TO 9330 ENDIF C GET SINES S9(N9) = 1. C COARSE SEED TABLE FOR SINES DO I = 1,3 S9(I*N9/4) = SIN(I*PI/8.) ENDDO C INITIAL HALF SECANT H9 = 1./2/COS(PI/16) C FILL SINE TABLE C9(4) = P-4 DO I = 1,P-4 C9(4) = C9(4)-1 V9(0) = 0 DO J = M9(C9(4)),N9-M9(C9(4)),M9(C9(4)+1) V9(1) = J + M9(C9(4)) S9(J) = H9*(S9(V9(1))+V9(0)) V9(0) = S9(V9(1)) ENDDO C HALF SECANT RECURSION H9 = 1/SQRT(2+1/H9) ENDDO C GET TANGENTS 9330 CONTINUE C9(0) = N9 - 1 DO I = 1,N9-1 T9(I) = (1-S9(C9(0)))/S9(I) C9(0) = C9(0) - 1 ENDDO T9(N9) = 1 9400 CONTINUE C FAST PERMUTE C FOR P = 2,3 PERMUTE DIRECTLY IF(P.EQ.2) THEN V9(9) = F(1) F(1) = F(2) F(2) = V9(9) GO TO 9500 ENDIF IF(P.EQ.3) THEN V9(9) = F(1) F(1) = F(4) F(4) = V9(9) V9(9) = F(3) F(3) = F(6) F(6) = V9(9) GO TO 9500 ENDIF IF(P.EQ.POLD) GOTO 9420 C FOR P=4, 5, 6 (Q9=2,3), SKIP STRUCTURE TABLE Q9 = P/2 C9(2) = M9(Q9) Q9 = Q9 + MOD(P,2) IF(Q9.EQ.2) THEN A9(0) = 0 A9(1) = 2 A9(2) = 1 A9(3) = 3 GO TO 9420 ENDIF IF(Q9.EQ.3) THEN A9(0) = 0 A9(1) = 4 A9(2) = 2 A9(3) = 6 A9(4) = 1 A9(5) = 5 A9(6) = 3 A9(7) = 7 GO TO 9420 ENDIF C SET UP STRUCTURE TABLE A9(0) = 0 A9(1) = 1 DO I = 2,Q9 DO J = 0,M9(I-1)-1 A9(J) = A9(J) + A9(J) A9(J+M9(I-1))=A9(J)+1 ENDDO ENDDO 9420 CONTINUE C PERMUTE DO I = 1,C9(2)-1 V9(4) = C9(2)*A9(I) V9(5) = I V9(6) = V9(4) V9(7) = F(V9(5)) F(V9(5)) = F(V9(6)) F(V9(6)) = V9(7) DO J =1,A9(I)-1 V9(5) = V9(5) + C9(2) V9(6) = V9(4) + A9(J) V9(7) = F(V9(5)) F(V9(5)) = F(V9(6)) F(V9(6)) = V9(7) ENDDO ENDDO C STAGES 1 & 2 C GET TWO-ELEMENT DHTs 9500 CONTINUE DO I = 0,N-2,2 T = F(I) - F(I+1) F(I) = F(I) + F(I+1) F(I+1) = T ENDDO C GET FOUR-ELEMENT DHTs DO I = 0,N-4,4 T = F(I) - F(I+2) F(I) = F(I) + F(I+2) F(I+2) = T T = F(I+1)-F(I+3) F(I+1) = F(I+1)+F(I+3) F(I+3) = T ENDDO IF(P.EQ.2) GOTO 9636 ! RETURN C STAGES 3, 4, ... U9 = C9(6) SX9 = 4 DO L9 = 2,C9(6) V9(2) = SX9 + SX9 U9 = U9-1 V9(3) = M9(U9-1) DO Q9 = 0,C9(5),V9(2) I = Q9 D9 = I+SX9 T = F(I) - F(D9) F(I) = F(I) + F(D9) F(D9) = T K9 = D9-1 DO J=V9(3),N9,V9(3) I = I + 1 D9 = I + SX9 E9 = K9 + SX9 T = F(D9) + F(E9)*T9(J) X9 = F(E9) - T*S9(J) Y9 = X9*T9(J) + T T = F(I) + Y9 F(D9) = F(I) - Y9 F(E9) = F(K9) + X9 F(K9) = F(K9) - X9 F(I) = T K9 = K9 -1 ENDDO ENDDO SX9= V9(2) ENDDO C SAVE P SO TABLE GENERATION CAN BE SKIPPED ON NEXT CALL 9636 CONTINUE POLD = P RETURN END