FUNCTION haar_matrix,N
;+
;Compute the Haar transform matrix H of size NxN
;using Gonzalez and Woods equation 7.1-26, page 361.
;
;N must be a power of 2.
;
;The Haar transform of an NxN array F is computed by
;T=H##F##H.
;
;-
m=Round(alog(N)/alog(2))
IF N NE 2^m THEN Message,'Argument must be a power of 2'

z=Findgen(N)/N
b=1.0 ;The value of 2^p in the loop below.
c=1.0 ;The value of 2^(p/2) in the loop below.
k=0; The row index in the loob below.
H=FltArr(N,N)
FOR p=0,m-1 DO BEGIN
	IF p EQ 0 THEN qlist=[0.,1.] ELSE qlist=Findgen(b)+1
	FOR j=0,N_Elements(qlist)-1 DO BEGIN
		IF p EQ 0 AND j EQ 0 THEN BEGIN
			z1=0.0 & z2=1.0 & z3=2.0
		ENDIF ELSE BEGIN
			z1=(qlist[j]-1)/b
			z2=(qlist[j]-0.5)/b
			z3=qlist[j]/b
		ENDELSE
		h1=c*(z1 LE z AND z LT z2)
		h2=-c*(z2 LE z AND z LT z3)
		H[*,k]=h1+h2
		k=k+1
	ENDFOR
	b=b*2
	c=sqrt(b)
ENDFOR
H=H/sqrt(N)
Return,H
END