Fast Hadamard Transform

Fortran Code

C
	SUBROUTINE HADAM(A,M)
C
C	ORDERED HADAMARD TRANSFORM FOR N=2**M, M INTEGER
C
C	AUTHOR: B. H. SUITS, MICHIGAN TECHNOLOGICAL UNIVERSITY, 1988.
C
C	INPUTS: A = ARRAY TO TRANSFORM (REAL)
C		   (REPLACED BY TRANSFORM)
C		M = LOG BASE TWO OF NUMBER OF DATA POINTS TO TRANSFORM
C			(INTEGER, M > -1)	
C	OUTPUT: A= TRANSFORMED DATA
C
C	SEE ALSO: M. KUNT, IEEE TRANSACTIONS ON COMPUTERS, VOL 24,
C		PAGE 1120 (1975).
C
C	GENERAL REFS: W. K. PRATT, DIGITAL IMAGE PROCESSING, 
C		(WILEY AND SONS, N.Y. 1978), PAGE 250 FF.
C	     W. K. PRATT, J. KANE, H. C. ANDREWS, PROC. IEEE, VOL 57,
C		NO. 1, PAGE 58 (1969).
C	REF. FOR BIT INVERSE PORTION: COOLEY, LEWIS, AND WELCH, IEEE
C		TRANSACTIONS E-12, NO. 1 (1965) AS QUOTED IN NUMERICAL
C		METHODS, BY DAHLQUIST AND BJORCK, PRENTICE HALL (1974),
C		PAGE 416.
C
C	COMMENT: TRANSFORM IS NOT NORMALIZED IN THE SENSE DESCRIBED IN
C		PRATT'S BOOK OR KUNT'S ROUTINE.  TO NORMALIZE, NEED
C		TO MULTIPLY EACH ELEMENT BY 1/SQRT(N).
C
C	FOR INTEGER, DOUBLE PRECISION, OR OTHER DATA TYPES, CHANGE
C	THE TYPE DECLARATION IN THE FOLLOWING LINE.
C
	REAL T,A(*)
C
	N=2**M
C
C	REORDER INTO GRAY SCALE ORDER
C
	K1=1
	K2=2
	DO 10 L=1,M-1
	  K0=K1
	  K1=K2
	  K2=2*K2
	  DO 10 J=K1,N,K2
	    DO 10 I=1,K0
	      IP=I+J
	      IPP=IP+K0
	      T=A(IP)
	      A(IP)=A(IPP)
	      A(IPP)=T
   10	CONTINUE
C
C	NOW DO BIT INVERSE SWAP AS IN FOURIER TRANSFORM
C
	NV2=N/2
	NM1=N-1
	J=1
	DO 17 I=1,NM1
	IF(I.LT.J) THEN
		T=A(J)
		A(J)=A(I)
		A(I)=T
	ENDIF
	K=NV2
   16	IF(K.GE.J) GOTO 17
	J=J-K
	K=K/2
	GOTO 16
   17	J=J+K
C
C	DO THE UNORDERED TRANSFORM ON THE REORDERED DATA TO GET
C	AN ORDERED TRANSFORM.
C
	L0=1
	DO 20 L=0,M-1
	  L1=L0
	  L0=2*L0
	  DO 20 J=1,L1
	    DO 20 I=J,N,L0
	      IP=I+L1
	      T=A(IP)
	      A(IP)=A(I)-T
   20	      A(I)=A(I)+T
C
	RETURN
	END