PRINT		OFF
			INCLUDE		'Traps.a'
			PRINT		ON
			
			TITLE	'FHTutils.a'
			STRING	ASIS					; to make DC.B work w/ normal 'strings'
			MC68881							; Needed for routines using FPU
			
; Included in this file are several routines needed in the computation of 1 and 2
; dimensional FFT's, such as bit-reversal and power spectrum calculation.  They
; may be freely duplicated and used (at your own risk).
;
; Questions about the implementation of these routines can be directed to me at:
; 
; Until 4/15/90:	Arlo Reeves
;					Thayer School of Engineering
;					Dartmouth College
;					Hanover, NH   03755
;					(603) 643 9076
;					BITNET: arlo@mac.dartmouth.edu
;
; Permanently:		Box 345
;					Mendocino, CA   95460
;					(707) 937 5686
; 

BitRevEZ	PROC	EXPORT

; function BitRevEZ(x: ptr; 
;				    length: integer):ErrCode;
;
; function BitRevEasy bitreverses the integer array of 'length' elements pointed
; to by x.  It does a getresource call to get a bitreversal lookup table of type
; 'BREV' of the appropriate size.  If BitRevEZ finds the resource, ErrCode is returned
; as 0, otherwise, it is returned as 1 (Nil resource handle).
; NOTE: Length must be a power of 2.

ErrCode		EQU		14						; ErrCode result A6 offset
x			EQU		10						;
length		EQU		8						;
NilHErr		EQU		1						;

UsedRegs	REG		A2/A4/D2-D4				;

			LINK	A6, #0					; No Locals
			MOVEM.L	UsedRegs, -(SP)			; Save registers
			MOVE.W	#0, ErrCode(A6)			; Default ErrCode is 0
			MOVE.L	x(A6), A1				; A1: x
			MOVE.W	length(A6), D0			; D0: length
			MOVE.W	D0, D1					; Duplicate D0
			ADD.W	#124, D1				; Add 124 to D1 to get Rsrc ID
			
			SUBQ.L	#4,SP					; Reserve space for Rsrc handle
			MOVE.L	#'BREV', -(SP)			; Pass RSRC Type and			
			MOVE.W	D1, -(SP)				; RSRC ID
			_GetResource					; Get the resource
			MOVE.L	(SP)+, D0				; Put Resource Handle in DO to set CCR
			BNE.S	continue				; just in case handle is nil
			MOVE.W	NilHErr, ErrCode(A6)	; set ErrCode = 1 => NilRHandle
			BRA.S	exit					; and exit
			
continue	MOVE.L	D0, A2					; use the handle
			MOVEA.L	A2, A4					; Copy the handle into A4
			MOVE.L	(A4), A2				; Dereference the Rsrc Handle in A2
			MOVE.W	(A2)+, D1				; Move numSwaps into D1
			SUBQ.W	#1, D1					; numSwaps -1 needed for DBRA
loop		MOVE.L	(A2)+, D2				; Displacement of first swap add
			MOVE.W	D2,D0					; move loword disp into D0
			SWAP	D2						; hiword disp in D2
			MOVE.W	0(A1,D0.W), D3			; 1st integer in D3
			MOVE.W	0(A1,D2.W), D4			; 2nd integer in D4
			MOVE.W	D3, 0(A1,D2.W)			; Store back in reversed order
			MOVE.W	D4, 0(A1,D0.W)			; "
			DBRA	D1, loop				; until numSwaps (left) = 0
			
			MOVE.L	A4, -(SP)				; push resource handle
			_ReleaseResource				; release Resource
			
exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDQ.L	#6, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'BITREVEZ'				; MacsBug Name
			
			ENDPROC							; BitRevEZ


BitRevQ		PROC	EXPORT

; procedure BitRevQ(x: ptr; 
;				    length: integer;
;				    BitRevTbl: ptr);
;
; function BitRevQuick bitreverses the integer array of 'length' elements pointed
; to by x.  It relies on the caller to provide a pointer to the appropriate
; resource of type 'BREV'.  If this pointer is nil, BitRevQ returns w/o having
; done anything to the input data.  By eliminating the getresource call, BitRevQ
; is quicker than BitRevEZ.
; NOTE: Length must be a power of 2.

data		EQU		14						; A6 offsets
length		EQU		12						;
BitRevTbl	EQU		8						;

UsedRegs	REG		A2/A4/D2-D4				;

			LINK	A6, #0					; No Locals
			MOVEM.L	UsedRegs, -(SP)			; Save registers
			MOVE.L	BitRevTbl(A6), D0		; to set condition codes
			BEQ.S	exit					; if nil then quit
			MOVEA.L	D0, A2					; A2: BitRevTbl
			MOVE.L	data(A6), A1			; A1: data
			MOVE.W	length(A6), D0			; D0: length
						
continue	MOVE.W	(A2)+, D1				; Move numSwaps into D1
			SUBQ.W	#1, D1					; numSwaps - 1 needed for DBRA
loop		MOVE.L	(A2)+, D2				; Displacement of first swap add
			MOVE.W	D2,D0					; move loword disp into D0
			SWAP	D2						; hiword disp in D2
			MOVE.W	0(A1,D0.W), D3			; 1st integer in D3
			MOVE.W	0(A1,D2.W), D4			; 2nd integer in D4
			MOVE.W	D3, 0(A1,D2.W)			; Store back in reversed order
			MOVE.W	D4, 0(A1,D0.W)			; "
			DBRA	D1, loop				; until numSwaps (left) = 0
			
exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#10, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'BITREVQ '				; MacsBug Name
			
			ENDPROC							; BitRevQ
			

BitRevRows	PROC	EXPORT

; procedure BitRevRows(iMatPtr: ptr;
;					   rowWords: integer;
;					   BitRevTbl: ptr);
;
; BitRevRows is similar to BitRevQ, but it BitReverses EACH row of the 
; rowWords x rowWords matrix pointed to by iMatPtr.  Like BitRevQ, it also relies on the
; calling routine to supply a pointer to the lookup table (resource BREV),
; instead of loading it itself (as in BitRev).
; A NIL BitRevTbl will cause this routine to return w/o having done anything.
			
iMatPtr		EQU		14						; A6 Offsets
rowWords	EQU		12						;
BitRevTbl	EQU		8						;

UsedRegs	REG		A2/D2-D6				;
			
			LINK	A6, #0					; No Locals
			MOVEM.L	UsedRegs, -(SP)			; Save registers
			MOVE.L	BitRevTbl(A6), D0		;
			BEQ.S	exit					; if BitRevTbl is NIL then exit
			MOVEA.L	D0, A1					; A1: BitRevTbl
			MOVEA.L	iMatPtr(A6), A0			; A0: iMatPtr
			MOVE.W	rowWords(A6), D0		; D0: rowWords
			
			MOVE.W	(A1)+, D1				; D1: numSwaps
			MOVEA.L	A1, A2					; A2: swap start
			SUBQ.W	#1, D1					; D1: numSwaps cntr.
			MOVE.W	D1, D2					; D2: numSwaps cntr store		
			MOVE.W	D0, D3					;
			EXT.L	D3						;
			ADD.L	D3, D3					; D3: rowOffset
			SUBQ.W	#1, D0					; D0: rowWords cntr
			
rowLoop		MOVE.L	(A1)+, D4				; D4: swap indices
			MOVE.W	(A0, D4.W), D5			;
			SWAP	D4						;
			MOVE.W	(A0, D4.W), D6			;
			MOVE.W	D5, (A0, D4.W)			;
			SWAP	D4						;
			MOVE.W	D6, (A0, D4.W)			;
			DBRA.W	D1, rowLoop				; for each Swap
			
			MOVE.W	D2, D1					; restore numSwaps cntr
			ADDA.L	D3, A0					; increment addr for next row
			MOVEA.L	A2, A1					; restore Tbl Addr
			DBRA.W	D0, rowLoop				; for each row

exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#10, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'BTRVROWS'				; MacsBug Name

			ENDPROC							; BitRevRows
			

ShiftRowsR	PROC	EXPORT

; procedure ShiftRowsR(iMatPtr: ptr;
;					   rowWords: integer;
;					   maxShift: integer;
;					   RowShfArr: ptr);
;
; Each row of the rowWords x rowWords integer matrix pointed to by iMatPtr is right
; shifted by a number of bits specified by the difference maxShift - RowShfArr[r]
; where RowShfArr (= array[0..rowWords -1] of integer) contains the number of
; bits each row has shifted (during transformation).
; During the right shift operation ShiftRowsR rounds the numbers instead of
; throwing out the bits shifted out.  It does this in a very simple way by
; checking the carry bit (the value of the last bit shifted out).  If it is
; set, then the number is rounded up; if it is not, the number is left as is.
; Because each word must be manipulated in this routine, there is no fast loop
; for the common 1 bit shift case.
; maxShift is the maximum number of bits by which any row was shifted during the
; transform (in FHT).  RowShfArr contains the number of bits by which each row
; was shifted.  To make all rows shifted by maxShift, each row is shifted here
; by maxShift - RowShfArr[rowIndex] bits.

iMatPtr		EQU		16						; Ptr to data
rowWords	EQU		14						; length of row
maxShift	EQU		12						; maximum of array
RowShfArr	EQU		8						; Array w/number of bits to shift right 

UsedRegs	REG		D2-D6/A2

start		LINK	A6, #0					; No Locals
			MOVEM.L	UsedRegs, -(SP)			; Save registers
			MOVE.L	iMatPtr(A6), A0			; A0: iMatPtr
			MOVEA.L	A0, A2					; A2: iMatPtr (row modified) 
			MOVE.L	RowShfArr(A6), A1		; A1: RowShfArr
			MOVE.W	rowWords(A6), D0		; D0: rowWords
			
			MOVE.W	D0, D4					;
			EXT.L	D4						;
			ADD.L	D4, D4					; D4: rowOffset
			SUBQ.W	#1, D0					; D0: rowWords cntr store
			MOVE.W	D0, D1					; D1: row cntr
			MOVE.W	D1, D2					; D2: col cntr
			MOVE.W	maxShift(A6), D3		; D3: maxShift
			
rowLoop		MOVEA.L	A2, A0					; A0: points to current row
			MOVE.W	(A1)+, D5				; D5: numBits already shifted
			SUB.W	D3, D5					;
			NEG.W	D5						; D5: maxShift - BitShfArr[r]
			BLE.S	ShfDone					; if numBits <= 0 then No Shift necc.
			
GenShift	MOVE.W	(A0), D6				; D4: data
			ASR.W	D5, D6					; shift it
			BCC.S	NoRound					;
			ADDQ.W	#1, D6					;
NoRound		MOVE.W	D6, (A0)+				; put it back
			DBRA.W	D2, GenShift			; for each col

ShfDone		ADDA.L	D4, A2					;
			MOVE.W	D0, D2					; restore col cntr
			DBRA.W	D1, rowLoop				; for each row

exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#12, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'SFTROWSR'				; MacsBug Name

			ENDPROC							; ShiftRowsR
			
MeanZero	PROC	EXPORT

; procedure MeanZero(srcPtr: ptr;
;					 bufSize: longint);
;
; Procedure MeanZero traverses the bufSize integer array pointed to by srcPtr
; and calculates the average value of the matrix elements.  It the traverses the
; array again, subtracting the average from each element.  This removes the DC
; component of the 2D transform of this matrix.
;
srcPtr		EQU		12						; A6 offsets
bufSize		EQU		8						;

UsedRegs	REG		D2-D6

			LINK	A6, #0					; no locals
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	srcPtr(A6), A0			; A0: srcPtr
			MOVE.L	bufSize(A6), D0			;
			MOVE.L	D0, D1					; D1: counter
			MOVEQ.L	#0, D2					; D2: lo longword of accumulator
			MOVEQ.L	#0, D3					; D3: hi longword of accumulator
			MOVEQ.L	#0, D4					; D4: zero needed for addx
			MOVEQ.L	#0, D5					; D5: Hi word must be zero

sumLoop		MOVE.W	(A0)+, D5				;
			ADD.L	D5, D2					;
			ADDX.L	D4, D3					; if carry bit, carry on
			SUBQ.L	#1, D1					;
			BNE.S	sumLoop					;
			
; have accumulated 64 bit sum in D3:D2; now simulate divide by subtracting DO repeatedly
			
			MOVEQ.L	#0, D1					; D1: zero needed for subx
divLoop		SUB.L	D0, D2					;
			SUBX.L	D1, D3					;
			BLT.S	divDone					;
			ADDQ.L	#1, D4					; D4: quotient
			BRA.S	divLoop					;
			
divDone		MOVE.L	srcPtr(A6), A0			;
subtLoop	SUB.W	D4, (A0)+				;
			SUBQ.L	#1, D0					;
			BNE.S	subtLoop				;
			
exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDQ.L	#8, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'MEANZERO'				; MacsBug Name

			ENDPROC							; MeanZero
			
			
FHTtoFFT1D	PROC	EXPORT

; procedure HtoF1D(H: iArrPtr;
;				   F: icArrPtr;
;				   numPts: integer);
;
; procedure FHTtoFFT1D takes a numPts long Hartley transform pointed to by x and
; converts it into a complex Fourier transform pointed to F.  F is of type cArrPtr:
;  icPoint = record
;    re, im: integer;
;   end;
;  icArr = array[0..maxStoreLess1] of icPoint;
;  icArrPtr = ^icArr;
; Both pointers H & F must be valid (not nil) and numPts must be a power of 2.
; The memory pointed to by F needn't be initialized.

H			EQU		14						; A6 offsets
F			EQU		10						;
numPts		EQU		8						;

UsedRegs	REG		D2-D6

			LINK	A6, #0					; no locals
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.W	numPts(A6), D1			; D1: numPts
			MOVEA.L	H(A6), A0				; A0: H
			MOVEA.L	F(A6), A1				; A1: F
			
; Do zero'th element of array first - special case

			MOVE.W	(A0), D0				; D0: ??? | H(0)
			SWAP	D0						;
			CLR.W	D0						; D0: H(0) | 0 = F(0)
			MOVE.L	D0, (A1)				; Store result

; Do N div 2 + 1'th element next - also a special case

			MOVE.W	0(A0, D1.W), D0			; N div 2 + 1 th element - middle of spectrum
			SWAP	D0						;
			CLR.W	D0						; D0: H(N div 2 + 1) | 0 = F(N div 2 + 1)
			MOVE.W	D1, D2					;
			ASL.W	#1, D2					;
			MOVE.L	D0, 0(A1, D2.W)			; Store result
			
; Now do rest of the points in a loop

			MOVE.W	D1, D0					; D0: N
			ASR.W	#1, D1					; 
			SUB.W	#2, D1					; D1: counter
			MOVE.W	#2, D2					; D2: f	    (integer index)
			MOVE.W	D0, D3					;
			ASL.W	#1, D3					;
			SUB.W	#2, D3					; D3: N - f (integer index)

loop		MOVE.W	0(A0, D2.W), D4			; D4: H(f)
			MOVE.W	0(A0, D3.W), D5			; D5: H(N - f)
			MOVE.W	D5, D6					; D6: H(N - f)
			ADD.W	D4, D5					; D5:   H(f) + H(N - f)  = re
			SUB.W	D4, D6					; D6: -(H(f) - H(N - f)) = im
			ASR.W	#1, D5					;
			ASR.W	#1, D6					;
			SWAP	D5						;
			MOVE.W	D6, D5					; D5: re | im
			
			MOVE.W	D2, D6					;
			ASL.W	#1, D6					;
			MOVE.L	D5, 0(A1, D6.W)			; store F(f)
			MOVE.W	D3, D6					;
			ASL.W	#1, D6					;
			NEG.W	D5						; CHS of im part
			MOVE.L	D5, 0(A1, D6.W)			; store F(N - f)
			
			ADDQ.W	#2, D2					;
			SUBQ.W	#2, D3					;			
			
			DBRA	D1, loop				;

exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#10, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'HtoF1D  '				; MacsBug Name

			ENDPROC							; HtoF1D
			
			
FFTtoFHT1D	PROC	EXPORT

; procedure FtoH1D(F: icArrPtr;
;				   H:	iArrPtr;
;				   numPts: integer);
;
; procedure FtoH1D converts an integer Fourier transform of length numPts into 
; a Hartley transform.  icP points to the complex integer array holding the Fourier
; transform and ixP points to the integer array Hartley transform result.
; For a definition of icArrPtr type, see HtoF1D.
; Again, no checking is done for nil pointers and numPts must be a power of 2.

F			EQU		14						; A6 offsets
H			EQU		10						;
numPts		EQU		8						;

UsedRegs	REG		D2-D3

			LINK	A6, #0					; no locals
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.W	numPts(A6), D1			; D1: numPts
			MOVEA.L	H(A6), A0				; A0: H
			MOVEA.L	F(A6), A1				; A1: F
			
			SUBQ.W	#1, D1					; D1: counter
loop		MOVE.L	(A1)+, D2				; D2: re | im
			MOVE.L	D2, D3					;
			SWAP	D3						; D3: im | re
			SUB.W	D2, D3					; D3: im | re - im 
			MOVE.W	D3, (A0)+				;
			DBRA	D1, loop				;

exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#10, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'FtoH1D  '				; MacsBug Name

			ENDPROC							; FtoH1D
			
psFHT1D		PROC	EXPORT

; procedure psFHT1D(h: iArrPtr;
;					ps: iArrPtr;
;					numPts: integer;
;					scaleType: integer);
;
; procedure psFHT1D computes the power spectrum from the numPts long Hartley transform
; pointed to by h and puts the result in the integer array ps.  The multiplication
; inherent in the ps calculations gives 32 bit results from 16 bit words, so these values
; are mapped into a their 16 bit destinations.  The mapping chosen is based on the
; scaleType argument as follows:
; scaleType			scaling implemented:
;		<=0				linear
;		1-9				nth root		(1st root = linear)
;		>=10			log
; In order to first determine the scale factor and offset of the linear mapping, 
; h is first scanned to determine the range of output by actually calculating 
; the power spectrum.  During the second scan, the computation again takes 
; place and the scaled data is stored in the dest.  This redundancy in computation
; could only be eliminated by an increase in storage requirements.
;
; WARNING: the use of this routine requires the existance of an FPU!

h			EQU		16						; A6 offsets
ps			EQU		12						;
numPts		EQU		10						;
scaleType	EQU		8						;

UsedRegs	REG		D2-D7

			LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			FMOVE.X	FP2, -(SP)				; save FP2
			MOVE.W	numPts(A6), D0			; D0: numPts
			MOVEA.L	h(A6), A0				; A0: h
			MOVEA.L	ps(A6), A1				; A1: ps
			
; first find the maximum and minimum of the range
			
			MOVE.L	#$7FFFFFFF, D6			; D6: min
			MOVE.L	#$80000000, D7			; D7: max
			
; first special case: H(0)
			MOVE.W	(A0), D1				; D1: H(0)
			MULS.W	D1, D1					;
			ADD.L	D1, D1					;
			CMP.L	D6, D1					;
			BGE.S	notMin1					;
			MOVE.L	D1, D6					;
notMin1		CMP.L	D7, D1					;
			BLE.S	notMax1					;
			MOVE.L	D1, D7					;

; second special case: H(N div 2 + 1)
notMax1		MOVE.W	0(A0, D0.W), D1			; D1: H(N Div 2 + 1)
			MULS.W	D1, D1					;
			ADD.L	D1, D1					;
			CMP.L	D6, D1					;
			BGE.S	notMin2					;
			MOVE.L	D1, D6					;
notMin2		CMP.L	D7, D1					;
			BLE.S	notMax2					;
			MOVE.L	D1, D7					;

; now for the rest of the spectrum
notMax2		MOVE.W	D0, D1					; D0: N
			ASR.W	#1, D1					; 
			SUB.W	#2, D1					; D1: counter
			MOVE.W	#2, D2					; D2: f	    (integer index)
			MOVE.W	D0, D3					;
			ASL.W	#1, D3					;
			SUB.W	#2, D3					; D3: N - f (integer index)
			
loop1		MOVE.W	0(A0, D2.W), D4			;
			MOVE.W	0(A0, D3.W), D5			;
			MULS.W	D4, D4					;
			MULS.W	D5, D5
			ADD.W	D5, D4					;
			CMP.L	D6, D4					;
			BGE.S	notMin3					;
			MOVE.L	D4, D6					;
notMin3		CMP.L	D7, D4					;
			BLE.S	notMax3					;
			MOVE.L	D4, D7					;

notMax3		ADDQ.W	#2, D2					;
			SUBQ.W	#2, D3					;
			DBRA	D1, loop1				;

; Now the values of max and min are recorded in D7 and D6 resp.
;
; Find which scaling to use : 	D6: 	min
;								FP1: 	mulFact
;								FP2:	1/n (for nth root calculation)
; for any x the mapping is: scaleType(x - xmin) * mulFact
; where scaleType produces the scaling desired and mulFact maps that range
; onto 0..maxint.

			SUB.L	D6, D7					; D7: max - min :: range
			BNE.S	findScale				;
			MOVEQ.L	#1, D7					; don't want a range of 0!
		
findScale	MOVE.W	scaleType(A6), D4		;
			BLE.S	linMap1					; scaleType is linear
			CMP.W	#10, D4					;
			BGE.S	logMap1					; scaleType is log

rootMap1	FMOVE.L	#1, FP2					; scaleType is nth root
			FDIV.W	D4, FP2					; FP1: 1/n
			FLOGN.L	D7, FP0					; FP0: ln(range)
			FMUL.X	FP2, FP0				; FP0: 1/n*ln(range)
			FETOX.X	FP0, FP0				; FP0: (range)^(1/n)
			MOVE.W	#2, D5					; D5 : scaleType flag
			BRA.S	GoOn1					;
			
linMap1		FMOVE.L	D7, FP0					;
			MOVE.W	#1, D5					; D5: scaleType flag
			BRA.S	GoOn1					;

logMap1		SUBQ.L	#1, D6					; min := min - 1 ;; to prevent log(0)
			ADDQ.L	#1, D7					; increases range by 1
			FLOG2.L	D7, FP0					; FP0: range !!! 6888x
			MOVE.W	#3, D5					; D5: scaleType flag
			
; Calculate the mulFact and put it in FP1

GoOn1		FMOVE.L	#$7FFF, FP1				;
			FDIV.X	FP0, FP1				; FP1: mulFact
			MOVE.W	D5, D7					; D7: scaleType flag
;
; D7 contains a scaleType flag: 1 -> linear; 2 -> nth root; 3 -> log
;
; now calculate actual ps values
; special case 1: H(0)

			MOVE.W	(A0), D2				;
			MULS.W	D2, D2					;
			ADD.L	D2, D2					;
			SUB.L	D6, D2					;
			MOVE.W	D7, D5					;
			SUBQ.W	#1, D5					;
			BEQ.S	linMap2					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap2				;
logMap2		FLOG2.L	D2, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn2					;
			
linMap2		FMOVE.L	D2, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn2					;
			
rootMap2	FLOGN.L	D2, FP0					;
			FMUL.X	FP2, FP0				; FP0: 1/n * ln(x)
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn2		FMOVE.W	FP0, D2					;
			MOVE.W	D2, (A1)				;
			
; special case 2: H(N div 2 + 1)

			MOVE.W	0(A0, D0.W), D2			;
			MULS.W	D2, D2					;
			ADD.L	D2, D2					;
			SUB.L	D6, D2					;
			MOVE.W	D7, D5					;
			SUBQ.W	#1, D5					;
			BEQ.S	linMap3					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap3				;
logMap3		FLOG2.L	D2, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn3					;
			
linMap3		FMOVE.L	D2, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn3					;
			
rootMap3	FLOGN.L	D2, FP0					;
			FMUL.X	FP2, FP0				; FP0: 1/n * ln(x)
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn3		FMOVE.W	FP0, D2					;
			MOVE.W	D2, 0(A1, D0.W)			;
			
; now do it for the rest
			MOVE.W	D0, D1					; D0: N
			ASR.W	#1, D1					; 
			SUB.W	#2, D1					; D1: counter
			MOVE.W	#2, D2					; D2: f	    (integer index)
			MOVE.W	D0, D3					;
			ASL.W	#1, D3					;
			SUB.W	#2, D3					; D3: N - f (integer index)
			
loop2		MOVE.W	0(A0, D2.W), D4			;
			MOVE.W	0(A0, D3.W), D5			;
			MULS.W	D4, D4					;
			MULS.W	D5, D5					;
			ADD.L	D4, D5					;
			SUB.L	D6, D5					;
			MOVE.W	D7, D4					;
			SUBQ.W	#1, D4					;
			BEQ.S	linMap4					;
			SUBQ.W	#1, D4					;
			BEQ.S	rootMap4				;
logMap4		FLOG2.L	D5, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn4					;
			
linMap4		FMOVE.L	D5, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn4					;
			
rootMap4	FLOGN.L	D5, FP0					;
			FMUL.X	FP2, FP0				; FP0: 1/n * ln(x)
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn4		FMOVE.W	FP0, D5					;
			MOVE.W	D5, 0(A1, D2.W)			;
			MOVE.W	D5, 0(A1, D3.W)			;
			
			ADDQ.W	#2, D2					;
			SUBQ.W	#2, D3					;
			DBRA	D1, loop2				;
			
			FMOVE.X	(SP)+, FP2				; restore FP2
exit		MOVEM.L	(SP)+,UsedRegs			; restore reg's
			UNLK	A6						; restore old A6
			MOVEA.L	(SP)+, A0				; pop return address
			ADDA.W	#12, SP					; discard parameters
			JMP	(A0)						; return to caller
			
			DC.B	'psFHT1D '				; MacsBug Name

			ENDPROC							; psFHT1D
			
			
ByteAvgRect	PROC	EXPORT

; function ByteAvgRect(baseAddr: ptr;
;					   pixelsPerLine: integer;
;					   roiWidth: integer;
;					   roiHeight: integer): integer;
;
; function ByteAvgRect returns the average byte of a rectangular part of a PixMap.
; baseAddr is the address of the first byte in the rectangle, pixelsPerLine is the
; width of the PixMap (right out of the Image's PicInfo record), roiWidth is the
; width of the rectangle (osRoiRect's width) and roiHeight is its height (osRoiRect's
; height).
; NOTE: lots of odd addresses generated here.

result		EQU		18						; A6 Offsets
baseAddr	EQU		14						;
pixelsPL	EQU		12						;
roiWidth	EQU		10						;
roiHeight	EQU		8						;

UsedRegs	REG		D2-D5/A2

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVEA.L	baseAddr(A6), A0		; A0: baseAddr 
			MOVEA.L	A0, A1					; A1: lineAddr stored
			MOVE.W	roiHeight(A6), D0		; D0: numLines (row)
			MOVE.W	roiWidth(A6), D1		; D0: rowBytes (col)
			MOVE.W	D1, D5					;
			MULS.W	D0, D5					; D5: numBytes
			SUBQ.W	#1, D0					; D0: row counter
			SUBQ.W	#1, D1					; D1: col counter
			MOVE.W	D1, D2					; D2: column count (for restore)
			MOVE.W	pixelsPL(A6), A2		; A2: pixelsPerLine (bytes offset)
			MOVEQ.L	#0, D3					; temp byte store
			MOVEQ.L	#0, D4					; sum of bytes
			
loop		MOVE.B	(A0)+, D3				;
			ADD.L	D3, D4					;
			DBRA.W	D1, loop				;
			
			MOVE.W	D2, D1					; restore column counter
			ADDA.L	A2, A1					;
			MOVEA.L	A1, A0					; update address
			DBRA.W	D0, loop				;
			
			CLR.W	D3						; D3: quotient
divLoop		SUB.L	D5, D4					;
			BLT.S	done					;
			ADDQ.W	#1, D3					;
			BRA.S	divLoop					;
			
done		MOVE.W	D3, result(A6)			;

exit		MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVEA.L	(SP)+, A0				;
			ADDA.L	#10, SP					;
			JMP		(A0)					;
			
			DC.B	'BTAVGRCT'				; Macsbug name
			
			ENDPROC
			
			
ByteRange	PROC	EXPORT

; function ByteRange(baseAddr: ptr;
;					 bufSize: longint):longint;
;
; function ByteRange scans the array of bytes starting at baseAddr that is bufSize
; bytes in length, recording max and min.  The max is returned in the high word, 
; the min in the low word of the result.
;
;
result		EQU		16						;
baseAddr	EQU		12						;
bufSize		EQU		8						;

UsedRegs	REG		D2-D3/D6-D7				;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	baseAddr(A6), A0		;
			MOVE.L	bufSize(A6), D0			;
			ASR.L	#1, D0					; fetching in bytes
			MOVE.W	#$00FF, D3				; D3: lo byte mask
			MOVE.W	D3, D6					; D6: min
			MOVE.W	#0, D7					; D7: max			
			
loop		MOVE.W	(A0)+, D1				;
			MOVE.W	D1, D2					;
			LSR.W	#8, D1					; D1: hibyte
			AND.W	D3, D2					; D2: lobyte
			
			CMP.W	D6, D2					; check hibyte
			BGE.S	notMin1					;
			MOVE.W	D2, D6					;
notMin1		CMP.W	D7, D2					;
			BLE.S	notMax1					;
			MOVE.W	D2, D7					;

notMax1		CMP.W	D6, D1					; check lobyte
			BGE.S	notMin2					;
			MOVE.W	D1, D6					;
notMin2		CMP.W	D7, D1					;
			BLE.S	notMax2					;
			MOVE.W	D1, D7					;
			
notMax2		SUBQ.L	#1, D0					;
			BGT.S	loop					;
			
			SWAP	D7						;
			MOVE.W	D6, D7					;
			MOVE.L	D7, result(A6)			;

exit		MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#8, SP					;
			JMP		(A0)					;
			
			DC.B	'BYTERNGE'				; MacsBug Name

			ENDPROC				
			

IntRange	PROC	EXPORT

; function IntRange(baseAddr: ptr;
;					bufSize: longint): result;
;
; function IntRange scans an array of integers bufSize long starting at baseAddr
; and returns the max in the high word, the min in the low word of the result
;
result		EQU		16						;
baseAddr	EQU		12						;
bufSize		EQU		8						;

UsedRegs	REG		D2/D6-D7				;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	baseAddr(A6), A0		;
			MOVE.L	bufSize(A6), D0			;
			ASR.L	#1, D0					;
			MOVE.W	#$7FFF, D6				; D6: min
			MOVE.W	#$8000, D7				; D7: max
loop		MOVE.L	(A0)+, D1				;
			MOVE.L	D1, D2					;
			SWAP	D2						;
			
			CMP.W	D6, D2					; check 1st #
			BGE.S	notMin1					;
			MOVE.W	D2, D6					;
notMin1		CMP.W	D7, D2					;
			BLE.S	notMax1					;
			MOVE.W	D2, D7					;

notMax1		CMP.W	D6, D1					; check 2nd #
			BGE.S	notMin2					;
			MOVE.W	D1, D6					;
notMin2		CMP.W	D7, D1					;
			BLE.S	notMax2					;
			MOVE.W	D1, D7					;
			
notMax2		SUBQ.L	#1, D0					;
			BGT.S	loop					;
			
			SWAP	D7						; D7: max | ?
			MOVE.W	D6, D7					; D7: max | min
			MOVE.L	D7, result(A6)			; 

exit		MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#8, SP					;
			JMP		(A0)					;
			
			DC.B	'INTRANGE'				; MacsBug Name

			ENDPROC						


AddSub2BufsF	PROC	EXPORT

;	procedure AddSub2BufsF(srcPtr1: 		ptr;
;						   srcPtr2: 		ptr;
;						   destPtr: 		ptr;
;						   srcScale1: 		integer;
;						   srcScale2: 		integer;
;					var	   destScale: 		integer;  (ptr)
;						   srcSize: 		longint;
;						   operation:		integer);
;
; procedure AddSub2BufsF adds/subtracts the contents of the FHTBuf pointed
; to by srcPtr2 to/from the FHTBuf pointed to by srcPtr1 and puts the result
; in the FHTBuf pointed to by destPtr (all Bufs are srcSize words in length).  
; If operation is >= 0 then the two arrays are added;
; If operation is <  0 then the two arrays are subtracted.
; srcScale is the power of 2 by which FHTBuf is multiplied to give a 
; 'correctly' scaled FHTBuf.  Since the scalings for the two buffers may 
; differ, one must be scaled according to the difference of srcScale1
; and srcScale2 before adding/subtracting.
; NOTE: FPU used here.

srcPtr1		EQU		30						; A6 Offsets
srcPtr2		EQU		26						;
destPtr		EQU		22						;
srcScale1	EQU		20						;
srcScale2	EQU		18						;
destScale	EQU		14						;
srcSize		EQU		10						;
operation	EQU		8						;

UsedRegs	REG		A2-A3					;
UsedFRegs	FREG	FP2-FP7					;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			FMOVEM.X	UsedFRegs, -(SP)	;
			
			FMOVE.W	srcScale1(A6), FP6		;
			FTWOTOX.X	FP6					; FP6: scale1
			FMOVE.W	srcScale2(A6), FP7		;
			FTWOTOX.X	FP7					; FP7: scale2
			
			MOVEA.L	destScale(A6), A0		; @destScale
			MOVEA.L	srcPtr1(A6), A1			;
			MOVEA.L	srcPtr2(A6), A2			;
			MOVE.L	srcSize(A6), D0			;
			FMOVE.L #$7FFFFFFF, FP4			; FP4: min
			FMOVE.L	#$80000000, FP5			; FP5: max
		
			TST.W	operation(A6)			;
			BLT.S	SubLoop1				;
			
AddLoop1	FMOVE.W	(A1)+, FP1				;
			FMOVE.W	(A2)+, FP2				;
			FMUL.X	FP6, FP1				;
			FMUL.X	FP7, FP2				;
			FADD.X	FP2, FP1				;
			FCMP.X	FP4, FP1				;
			FBGE.W	notMin1					;
			FMOVE.X	FP1, FP4				;
notMin1		FCMP.X	FP5, FP1				;
			FBLE.W	notMax1					;
			FMOVE.X	FP1, FP5				;
notMax1		SUBQ.L	#1, D0					;
			BGT.S	AddLoop1				;
			BRA.S	GoOn1					;
			
SubLoop1	FMOVE.W	(A1)+, FP1				;
			FMOVE.W	(A2)+, FP2				;
			FMUL.X	FP6, FP1				;
			FMUL.X	FP7, FP2				;
			FSUB.X	FP2, FP1				;
			FCMP.X	FP4, FP1				;
			FBGE.W	notMin2					;
			FMOVE.X	FP1, FP4				;
notMin2		FCMP.X	FP5, FP1				;
			FBLE.W	notMax2					;
			FMOVE.X	FP1, FP5				;
notMax2		SUBQ.L	#1, D0					;
			BGT.B	SubLoop1				;

GoOn1		FABS.X	FP4						; FP4: abs(min)
			FABS.X	FP5						; FP5: abs(max)
			FCMP.X	FP4, FP5				;
			FBGT.W	Bigger1					;
			FMOVE.X	FP4, FP5				;
Bigger1		FMOVE.W	#0, FP4					; FTST wouldn't work in MPW 2.02's asm!
			FCMP.X	FP4, FP5				;
			FBNE.W	notZero1				;
			FMOVE.W	#1, FP5					;
notZero1	CLR.W	D0						; D0: destScale
			FLOG2.X	FP5						;
			FMOVE.W	#13, FP4				;
			FCMP.X	FP4, FP5				;
			FBLT.W	GoOn2					;
			FSUB.X	FP4, FP5				;
			FADD.W	#1, FP5					;
			FINTRZ.X	FP5					; Round Up to next integer
			FMOVE.W	FP5, D0					; D0: destScale
GoOn2		MOVE.W	D0, (A0)				; destScale stored
			FMOVE.W	D0, FP4					; FP4: destScale
			FTWOTOX.X	FP4					; FP4: 2^destScale (divFactor)
			FMOVE.W	#1, FP5					;
			FDIV.X	FP4, FP5				; FP5: 1/(2^destScale (mulFactor)
			
			MOVEA.L	srcPtr1(A6), A1			;
			MOVEA.L	srcPtr2(A6), A2			;
			MOVEA.L	destPtr(a6), A3			;
			MOVE.L	srcSize(A6), D0			;
			TST.W	operation(A6)			;
			BLT.S	SubLoop2				;
			
AddLoop2	FMOVE.W	(A1)+, FP1				;
			FMOVE.W	(A2)+, FP2				;
			FMUL.X	FP6, FP1				;
			FMUL.X	FP7, FP2				;
			FADD.X	FP2, FP1				;
			FMUL.X	FP5, FP1				;
			FMOVE.W	FP1, (A3)+				;
			SUBQ.L	#1, D0					;
			BGT.B	AddLoop2				;
			BRA.S	exit					;
			
SubLoop2	FMOVE.W	(A1)+, FP1				;
			FMOVE.W	(A2)+, FP2				;
			FMUL.X	FP6, FP1				;
			FMUL.X	FP7, FP2				;
			FSUB.X	FP2, FP1				;
			FMUL.X	FP5, FP1				;
			FMOVE.W	FP1, (A3)+				;
			SUBQ.L	#1, D0					;
			BGT.B	SubLoop2				;
			
exit		FMOVEM.X	(SP)+, UsedFRegs	;
			MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDA.W	#26, SP					;
			JMP		(A0)					;
			
			DC.B	'ADSB2BFF'				; MacsBug Name

			ENDPROC


Hcdc2BufsF	PROC	EXPORT					

;	procedure Hcdc2BufsF(srcPtr1: 	ptr;
;						 srcPtr2: 	ptr;
;						 destPtr: 	ptr;
;						 srcScale1:	integer;
;						 srcScale2:	integer;
;					var	 destScale:	integer;	(ptr)
;						 rowWords: 	integer;
;						 operation:	integer);
;
; procedure Hcdc2BufsF (for Hartley Convolve or Deconvolve or Correlate two Buffers w/ FPU)
; is a feature packed routine that performs the Hartley analog of complex multiplication
; or division or conjugate multiplication on 2 rowWords x rowWords integer Hartley
; transform buffers pointed to by srcPtr1 and srcPtr2, placing the result in
; the Hartley transform buffer pointed to by destPtr.  All three of these computations
; are so similar (and the addressing so complex) that they have been lumped togehter
; in one routine.
; The operations on H1 and H2 (src1 and src2) can be expressed using Bracewell's notation
; (p 43, 'The Hartley Transform') as
; (Convolution):  	 H1(+f)*H2e + H1(-f)*H2o
; (Deconvolution):	(H1(+f)*H2e - H1(-f)*H2o)/(H2(+f)^2 + H2(-f)^2)
; (Correlation):	 H1(+f)*H2e - H1(-f)*H2o
; The operation desired is specified by the operation paramter as follows:
;				operation			function performed
;					< 0					Hartley Multiplication			 (for Convolution)
;					= 0					Hartley Division       			 (for Deconvolution)
;					> 0					Hartley Conjugate Multiplication (for Correlation)
; As in PSFHT2D and other routines, this calculation must be performed twice.  The first
; time through, only the maximum and minimum values of the output are recorded.  From the
; extrema, a linear mapping from the 32 bit output into the 16 bit output buffer is 
; devised.  On the second time through, this mapping is used to store the computed result
; into the output buffer.
; The order of arguments is not important for multiplication (convolution) as it commutes,
; however, for division and conjugate multiplication the operations are
; src1/src2 and src1(src2)* respectivley, where * indicates conjugation.
; NOTE: this accomplishes with one loop what PSFHT2D does with two loops and two special
; 		cases; the addressing calculation involved here does not take long in comparison
;		with PSFHT2D, so that routine should be modified to use this more compact format.
; NOTE: rowWords MUST be a power of 2; this routine is designed to work exclusively on
;	    square, two-dimensional integer matrices.
; NOTE: Autocorrelation needs only half the operations that correlation needs; this is
; 		not taken advantage of here; THIS ROUTINE IS NOT OPTIMIZED AT ALL.
; NOTE: FPU used here.

srcPtr1		EQU		28						; A6 Offsets
srcPtr2		EQU		24						;
destPtr		EQU		20						;
srcScale1	EQU		18						;
srcScale2	EQU		16						;
destScale	EQU		12						;
rowWords	EQU		10						;
operation	EQU		8						;

UsedRegs	REG		D2-D7/A2-A4
UsedFRegs	FREG	FP2-FP7

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			FMOVEM.X	UsedFRegs, -(SP)	;
			
			FMOVE.W	srcScale1(A6), FP6		;
			FTWOTOX.X	FP6					; FP6: scale1
			FMOVE.W	srcScale2(A6), FP7		;
			FTWOTOX.X	FP7					; FP7: scale2
						
			MOVE.W	operation(A6), D0		; D0: * | operation
			SWAP	D0						;
			MOVE.W	rowWords(A6), D0		; D0: operation | rowWords
			MOVEA.L	srcPtr1(A6), A1			; A1: srcPtr1
			MOVEA.L	srcPtr2(A6), A2			; A2: srcPtr2
			
			MOVE.W	#15, D2					;
logLoop		BTST.L	D2, D0					;
			DBNE.W	D2, logLoop				; D2: log2(rowWords)
			MOVE.W	D2, D1					;
			SWAP	D1						;
			MOVE.W	D0, D1					;
			SUBQ.W	#1, D1					; D1: log2(rowWords) | mod mask & column count
			MOVE.W	D1, D2					; D2: column counter
			MOVE.W	D1, D3					; D3: row counter		
			FMOVE.L	#$7FFFFFFF, FP4			; FP4: min
			FMOVE.L	#$80000000, FP5			; FP5: max
			
RCLoop1		MOVE.W	D3, D4					;
			EXT.L	D4						;
			SWAP	D1						;
			ASL.L	D1, D4					; D4: N * row
			SWAP	D1						;
			MOVE.W	D2, D5					;
			EXT.L	D5						;
			ADD.L	D5, D4					; D4: N * row + col
			ASL.L	#1, D4					; D4: even (word) address of H(r, c)
			MOVEA.L	D4, A3					; A3: address of H(row, col)
			
			MOVE.L	D3, D4					;
			MOVE.W	D0, D5					; D5: N
			SUB.W	D4, D5					; D5: N - row
			AND.W	D1, D5					; D4: (N - row) mod N = modRow
			EXT.L	D5						;
			SWAP	D1						;
			ASL.L	D1, D5					; D4: N * modRow
			SWAP	D1						;
			MOVE.W	D0, D6					; D6: N
			SUB.W	D2, D6					; D6: N - col
			AND.W	D1, D6					; D6: (N - col) mod N = modCol
			EXT.L	D6						;
			ADD.L	D6, D5					; D4: N * modRow + modCol
			ASL.L	#1, D5					; D4: even (word) address of H(modRow, modCol)
			MOVEA.L	D5, A4					; A4: address of H(modRow, modCol)
			
; address calculation done

			FMOVE.W	0(A2, A3.L), FP0		; FP0: H2(r,c)
			FMUL.X	FP7, FP0				;
			FMOVE.X	FP0, FP1				;
			FMOVE.W	0(A2, A4.L), FP2		; FP2: H2(modRow, modCol)
			FMUL.X	FP7, FP2				;
			FADD.X	FP2, FP0				;
			FSUB.X	FP2, FP1				;
			FDIV.W	#2, FP0					; FP0: H2even
			FDIV.W	#2, FP1					; FP1: H2odd
			
			MOVE.L	D0, D7					; because a swap will change CCR's Z & N bits
			SWAP	D7						;
			TST.W	D7						;
			BGT.W	Corr1					;
			BLT.B	Conv1					;
			
Deconv1		FMOVE.X	FP0, FP2				;
			FMUL.X	FP2, FP2				; FP2: sqr(H2e)
			FMOVE.X	FP1, FP3				;
			FMUL.X	FP3, FP3				; FP3: sqr(H2o)
			FADD.X	FP3, FP2				; FP2: sqr(H2e) + sqr(H2o) [= (H2^2 + H2(-)^2)/2]
			FMUL.W	#2, FP2					; FP2: = sqr(H2) + sqr(H2(-)) = denominator
			FBNE.W	notZero1				;
			FMOVE.W	#1, FP2					; to prevent divide by zero
notZero1	FMOVE.W	0(A1, A3.L), FP3		;
			FMUL.X	FP6, FP3				; FP3: H1(r,c)
			FMUL.X	FP3, FP0				; FP0: H1(r,c)* H2even
			FMOVE.W	0(A1, A4.L), FP3		;
			FMUL.X	FP6, FP3				; FP3: H1(modRow, modCol)
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)* H2odd
			FSUB.X	FP1, FP0				; FP0: deconv numerator
			FDIV.X	FP2, FP0				; FP0: ratio
			FCMP.X	FP4, FP0				;
			FBGE.W	notMin1					;
			FMOVE.X	FP0, FP4				;
notMin1		FCMP.X	FP5, FP0				;
			FBLE.W	notMax1					;
			FMOVE.X	FP0, FP5				;
notMax1		BRA.S	GoOn1					;

Conv1		FMOVE.W	0(A1, A3.L), FP2		; FP2: H1(r,c)
			FMOVE.W	0(A1, A4.L), FP3		; FP3: H1(modRow, modCol)
			FMUL.X	FP6, FP2				;
			FMUL.X	FP6, FP3				;
			FMUL.X	FP2, FP0				; FP0: H1(r,c)*H2even
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)*H2odd
			FADD.X	FP1, FP0				;
			FCMP.X	FP4, FP0				;
			FBGE.W	notMin2					;
			FMOVE.X	FP0, FP4				;
notMin2		FCMP.X	FP5, FP0				;
			FBLE.W	notMax2					;
			FMOVE.X	FP0, FP5				;
notMax2		BRA.S	GoOn1					;

Corr1		FMOVE.W	0(A1, A3.L), FP2		; FP2: H1(r,c)
			FMOVE.W	0(A1, A4.L), FP3		; FP3: H1(modRow, modCol)
			FMUL.X	FP6, FP2				;
			FMUL.X	FP6, FP3				;
			FMUL.X	FP2, FP0				; FP0: H1(r,c)*H2even
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)*H2odd
			FSUB.X	FP1, FP0				; 
			FCMP.X	FP4, FP0				;
			FBGE.W	notMin3					;
			FMOVE.X	FP0, FP4				;
notMin3		FCMP.X	FP5, FP0				;
			FBLE.W	GoOn1					;
			FMOVE.X	FP0, FP5				;
			
; calculation complete

GoOn1		DBRA.W	D2, RCLoop1				;

			MOVE.W	D1, D2					; Restore column counter
			DBRA.W	D3, RCLoop1				;
			
; First traversal done - max and min found
			
			MOVEA.L	destScale(A6), A0		; A0: @destScale
			FABS.X	FP4						; FP4: abs(min)
			FABS.X	FP5						; FP5: abs(max)
			FCMP.X	FP4, FP5				;
			FBGT.W	Bigger2					;
			FMOVE.X	FP4, FP5				;
Bigger2		FMOVE.W	#0, FP4					; FTST wouldn't work in MPW 2.02's asm!
			FCMP.X	FP4, FP5				;
			FBNE.W	notZero2				;
			FMOVE.W	#1, FP5					;
notZero2	CLR.W	D5						; D5: destScale
			FLOG2.X	FP5						;
			FMOVE.W	#13, FP4				;
			FCMP.X	FP4, FP5				;
			FBLT.W	GoOn2					;
			FSUB.X	FP4, FP5				;
			FADD.W	#1, FP5					;
			FINTRZ.X	FP5					; Round Up to next integer
			FMOVE.W	FP5, D5					; D5: destScale
GoOn2		MOVE.W	D5, (A0)				; destScale stored
			FMOVE.W	D5, FP4					; FP4: destScale
			FTWOTOX.X	FP4					; FP4: 2^destScale (divFactor)
			FMOVE.W	#1, FP5					;
			FDIV.X	FP4, FP5				; FP5: 1/(2^destScale (mulFactor)
			
; mapping set up - restore row and column counters, etc. for next traversal

			MOVEA.L	srcPtr1(A6), A1			;
			MOVEA.L	srcPtr2(A6), A2			;
			MOVEA.L	destPtr(A6), A0			;
			MOVE.W	D1, D2					; D2: col counter
			MOVE.W	D1, D3					; D3: row counter
			
RCLoop2		MOVE.W	D3, D4					;
			EXT.L	D4						;
			SWAP	D1						;
			ASL.L	D1, D4					; D4: N * row
			SWAP	D1						;
			MOVE.W	D2, D5					;
			EXT.L	D5						;
			ADD.L	D5, D4					; D4: N * row + col
			ASL.L	#1, D4					; D4: even (word) address of H(r, c)
			MOVEA.L	D4, A3					; A3: address of H(row, col)
			
			MOVE.L	D3, D4					;
			MOVE.W	D0, D5					; D5: N
			SUB.W	D4, D5					; D5: N - row
			AND.W	D1, D5					; D4: (N - row) mod N = modRow
			EXT.L	D5						;
			SWAP	D1						;
			ASL.L	D1, D5					; D4: N * modRow
			SWAP	D1						;
			MOVE.W	D0, D6					; D6: N
			SUB.W	D2, D6					; D6: N - col
			AND.W	D1, D6					; D6: (N - col) mod N = modCol
			EXT.L	D6						;
			ADD.L	D6, D5					; D4: N * modRow + modCol
			ASL.L	#1, D5					; D4: even (word) address of H(modRow, modCol)
			MOVEA.L	D5, A4					; A4: address of H(modRow, modCol)
			
; address calculation done

			FMOVE.W	0(A2, A3.L), FP0		; FP0: H2(r,c)
			FMUL.X	FP7, FP0				;
			FMOVE.X	FP0, FP1				;
			FMOVE.W	0(A2, A4.L), FP2		; FP2: H2(modRow, modCol)
			FMUL.X	FP7, FP2				;
			FADD.X	FP2, FP0				;
			FSUB.X	FP2, FP1				;
			FDIV.W	#2, FP0					; FP0: H2even
			FDIV.W	#2, FP1					; FP1: H2odd
			
			MOVE.L	D0, D7					; because a swap will change CCR's Z & N bits
			SWAP	D7						;
			TST.W	D7						;
			BGT.S	Corr2					;
			BLT.S	Conv2					;
			
Deconv2		FMOVE.X	FP0, FP2				;
			FMUL.X	FP2, FP2				; FP2: sqr(H2e)
			FMOVE.X	FP1, FP3				;
			FMUL.X	FP3, FP3				; FP3: sqr(H2o)
			FADD.X	FP3, FP2				; FP2: sqr(H2e) + sqr(H2o) [= (H2^2 + H2(-)^2)/2]
			FMUL.W	#2, FP2					; FP2: = sqr(H2) + sqr(H2(-)) = denominator
			FBNE.W	notZero3				;
			FMOVE.W	#1, FP2					; to prevent divide by zero
notZero3	FMOVE.W	0(A1, A3.L), FP3		;
			FMUL.X	FP6, FP3				; FP3: H1(r,c)
			FMUL.X	FP3, FP0				; FP0: H1(r,c)* H2even
			FMOVE.W	0(A1, A4.L), FP3		;
			FMUL.X	FP6, FP3				; FP3: H1(modRow, modCol)
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)* H2odd
			FSUB.X	FP1, FP0				; FP0: deconv numerator
			FDIV.X	FP2, FP0				; FP0: ratio
			BRA.S	GoOn3					;

Conv2		FMOVE.W	0(A1, A3.L), FP2		; FP2: H1(r,c)
			FMOVE.W	0(A1, A4.L), FP3		; FP3: H1(modRow, modCol)
			FMUL.X	FP6, FP2				;
			FMUL.X	FP6, FP3				;
			FMUL.X	FP2, FP0				; FP0: H1(r,c)*H2even
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)*H2odd
			FADD.X	FP1, FP0				;
			BRA.S	GoOn3					;

Corr2		FMOVE.W	0(A1, A3.L), FP2		; FP2: H1(r,c)
			FMOVE.W	0(A1, A4.L), FP3		; FP3: H1(modRow, modCol)
			FMUL.X	FP6, FP2				;
			FMUL.X	FP6, FP3				;
			FMUL.X	FP2, FP0				; FP0: H1(r,c)*H2even
			FMUL.X	FP3, FP1				; FP1: H1(modRow, modCol)*H2odd
			FSUB.X	FP1, FP0				; 
			
; calculation complete

GoOn3		FMUL.X	FP5, FP0				; FP0: result scaled down
			FMOVE.W	FP0, 0(A0, A3.L)		; result stored
			DBRA.W	D2, RCLoop2				;

			MOVE.W	D1, D2					; Restore column counter
			DBRA.W	D3, RCLoop2				;

exit		FMOVEM.X	(SP)+, UsedFRegs	;
			MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDA.L	#24, SP					;
			JMP		(A0)					;
			
			DC.B	'HCDC2BFF'				; MacsBug Name

			ENDPROC						


ClipMinMax	PROC	EXPORT

; procedure ClipMinMax(baseAddr: ptr	;
;					   bufSize: longint);
;
; procedure ClipMinMax is used on a PixMap to convert all of its 0 pixels to 1 and
; all of its 255 pixels to 254, thereby avoiding the sacred foreColor and bgColor
; extremeties.  baseAddr points to an array of bytes bufSize long.

baseAddr	EQU		12						;
bufSize		EQU		8						;

UsedRegs	REG		D2-D3					;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	baseAddr(A6), A0		; A0: baseAddr
			MOVE.L	bufSize(A6), D0			; D0: bufSize
			ASR.L	#1, D0					; D0: bufSize words
			MOVE.W	#$00FF, D1				; D1: low Byte mask
			
loop		MOVE.W	(A0), D2				;
			MOVE.W	D2, D3					;
			LSR.W	#8, D3					; D3: 1st byte
			AND.W	D1, D2					; D2: 2nd byte

			BNE.S	notZero1				; check 2nd byte
			MOVE.W	#1, D2					;
notZero1	CMPI.W	#254, D2				;
			BLE.S	notBig1					;			
			MOVE.W	#254, D2				;

notBig1		TST.W	D3						; check 1st byte
			BNE.S	notZero2				;
			MOVE.W	#1, D3					;
notZero2	CMPI.W	#254, D3				;
			BLE.S	notBig2					;			
			MOVE.W	#254, D3				;

notBig2		LSL.W	#8, D3					; store bytes back
			OR.W	D2, D3					;
			MOVE.W	D3, (A0)+				;
			SUBQ.L	#1, D0					;
			BGT.S	loop					;
			
exit		MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#8, SP					;
			JMP		(A0)					;
			
			DC.B	'CPMINMAX'				; MacsBug Name

			ENDPROC
			

DblMem		PROC	EXPORT

; procedure DblMem(srcPtr: 	ptr;
;				   destPtr:	ptr;
;				   srcSize:	longint);
;
; This procedure takes an input block of memory, a sequence of bytes (e.g. from
; an 8 bit deep pixel map), and copies it into a block of memory twice the size -
; a 16 bit deep pixel map. The byte with value 1 is mapped into 32, while the byte
; with value 255 is mapped into 8160; the values are simply left shifted by 5 bits.
; srcSize is measured in bytes.

srcPtr		EQU		16						; A6 offsets
destPtr		EQU		12						;
srcSize		EQU		8						;

UsedRegs	REG		D2-D4

			LINK	A6, #0					; no locals (disco surfer!)
			MOVEM.L	UsedRegs, -(SP)			;
			MOVEA.L	srcPtr(A6), A0			; A0: srcPtr
			MOVEA.L	destPtr(A6), A1			; A1: destPtr
			MOVE.L	srcSize(A6), D0			; D0: srcSize
			ASR.L	#1, D0					; D0: srcSize in words
			MOVE.W	#5, D1					; D1: num bits to shift
			MOVE.W	#$00FF, D2				; D2: low byte mask
			
loop		MOVE.W	(A0)+, D3				;
			MOVE.W	D3, D4					;
			LSR.W	#8, D3					; D3: 1st byte
			AND.W	D2, D4					; D4: 2nd byte
			ASL.W	D1, D3					;
			ASL.W	D1, D4					;
			SWAP	D3						;
			MOVE.W	D4, D3					;
			MOVE.L	D3, (A1)+				;
			SUBQ.L	#1, D0					;
			BGT.S	loop					;

exit		MOVEM.L	(SP)+, UsedRegs			; restore registers
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDA.W	#12, SP					; pop parameters
			JMP		(A0)					;
			
			DC.B	'DBLMEM  '				; MacsBug Name

			ENDPROC
			
			
IntToByteF	PROC	EXPORT

; procedure IntToByteF(srcPtr: ptr;
;				 	   destPtr: ptr;
;                      rowWords: integer);
;
; procedure IntToByte scans the rowWords^2 matrix of integers pointed 
; to by srcPtr recording the maximum and minimum values.
; These are then used to define a linear mapping onto the byte range 0-255.
; A second pass over the integer array is then made to map the integers
; into the rowWords^2 byte matrix pointed to by destPtr.
; IntToByteT differs from IntToByte in that it uses the FPU for the linear
; mapping and is therefore slower, but simpler to read, smaller and more accurate.

srcPtr		EQU		14						; A6 offsets
destPtr		EQU		10						;
rowWords	EQU		8						;

UsedRegs	REG		D2/D6-D7				;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	srcPtr(A6), A0			; A0: srcPtr
			MOVE.W	rowWords(A6), D0		; D0: rowWords
			
			MOVE.W	#$7FFF, D6				; D6: minVal
			MOVEQ.L	#0, D7					;
			BSET.L	#15, D7					; D7: maxVal

			MULS.W	D0, D0					; D0: srcSize
			MOVE.L	D0, D1					; D1: count
						
scanLoop	MOVE.W	(A0)+, D2				;
			CMP.W	D6, D2					;
			BGT.S	notMin					;
			MOVE.W	D2, D6					;
notMin		CMP.W	D7, D2					;
			BLE.S	notMax					;
			MOVE.W	D2, D7					;
notMax		SUBQ.L	#1, D1					;
			BGT.S	scanLoop				; max & min found
			SUB.W	D6, D7					; D7: range
			BNE.S	notZero					;
			MOVE.W	#1, D7					;
			
notZero		FMOVE.W	D7, FP0					; FPO: range
			FMOVE.W	#255, FP1				;
			FDIV.X	FP0, FP1				; FP1: mul Factor (scale)
			
			MOVE.L	srcPtr(A6), A0			; A0: srcPtr
			MOVE.L	destPtr(A6), A1			; A1: destPtr
			
			ASR.L	#1, D0					; D0: count
storeLoop	MOVE.L	(A0)+, D1				; D1: 2nd word (2 words at at time)
			MOVE.L	D1, D2					;
			SWAP	D2						; D2: 1st word
			SUB.W	D6, D1					;
			SUB.W	D6, D2					;
			FMOVE.W	D2, FP0					;
			FMUL.X	FP1, FP0				;
			FMOVE.W	FP0, D2					;
			
			FMOVE.W	D1, FP0					;
			FMUL.X	FP1, FP0				;
			FMOVE.W	FP0, D1					;
			
			LSL.W	#8, D2					;
			OR.W	D2, D1					;
			MOVE.W	D1, (A1)+				;
			
			SUBQ.L	#1, D0					;
			BGT.S	storeLoop				;
			
exit		MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDA.W	#10, SP					;
			JMP		(A0)					;
			
			DC.B	'INTOBYTF'				; MacsBug Name

			ENDPROC						
			
			
SwapBBlock	PROC	EXPORT

; procedure SwapBBlock(matPtr: ptr;
;                      rowBytes: integer);

; procedure SwapBBlock rearranges the square matrix of BYTES pointed to by matPtr.
; rowBytes must be a power of 2 >= 4. The moves are word sized.
; It swaps the matrix's 2nd and 4th quadrants and its 1st and 3rd quadrants, so
; that the transform is centered in the matrix.

matPtr		EQU		10						; A6 offsets
rowBytes	EQU		8						;

UsedRegs	REG		D2-D7/A2-A3				;

start		LINK	A6, #0					;
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	matPtr(A6), A0			;
			MOVE.W	rowBytes(A6), D0		;
			
			EXT.L	D0						; D0: rowBytes
			MOVE.L	D0, D1					;
			ASR.L	#1, D1					; D1: rowBytes div 2
			MOVE.W	#15, D2					;
logLoop		BTST.L	D2, D0					;
			DBNE.W	D2, logLoop				; D2: ?? | Log2(rowBytes)
			SUBQ.W	#1, D2					; 
			MOVE.L	D0, D3					;
			ASL.L	D2, D3					; D3: rowBytes^2 div 2
			
			MOVEA.L	A0, A1					; A0: Quad 2 ptr
			MOVEA.L	A0, A2					;
			ADDA.L	D1, A1					; A1: Quad 1 ptr
			ADDA.L	D3, A2					; A2: Quad 3 ptr
			MOVEA.L	A2, A3					;
			ADDA.L	D1, A3					; A3: Quad 4 ptr
			
			MOVE.W	D1, D3					; this must be pos, word sized.
			ASR.W	#1, D3					; D2: rowWords
			SUBQ.W	#1, D3					; D2: rowWords - 1
			MOVE.W	D3, D2					;
			SWAP	D2						;
			MOVE.W	D3, D2					; D2: col | col
			MOVE.W	D1, D3					; 
			SUBQ.W	#1, D3					; D3: row

loop		MOVE.W	(A0), D4				;
			MOVE.W	(A1), D5				;
			MOVE.W	(A2), D6				;
			MOVE.W	(A3), D7				;
			MOVE.W	D4, (A3)+				;
			MOVE.W	D5, (A2)+				;
			MOVE.W	D6, (A1)+				;
			MOVE.W	D7, (A0)+				;
			
			DBRA.W	D2, loop				;
			
			MOVE.L	D2, D4					;
			SWAP	D4						;
			MOVE.W	D4, D2					; restore col counter
			
			ADDA.L	D1, A0					;
			ADDA.L	D1, A1					;
			ADDA.L	D1, A2					;
			ADDA.L	D1, A3					; update the address ptrs

			DBRA.W	D3, loop				;
			
			MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#6, SP					;
			JMP		(A0)					;
			
			DC.B	'SWPBBLOK'				; MacsBug Name

			ENDPROC						
			
			
Transpose	PROC	EXPORT

; procedure transpose(matPtr: ptr; 
;					  rowWords: integer);

; procedure transpose transposes a rowWords x rowWords matrix of words.
; WARNING: the integer block pointed to by matPtr must have rowWords^2
; words allocated to it. (it must be a square matrix).  No checking for
; this is done.

matPtr		EQU		10						; A6 offsets
rowWords	EQU		8						;

UsedRegs	REG		D2-D7					;

			LINK	A6, #0					; no locals (my wave!)
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	matPtr(A6), A0			; Base Address
			MOVE.W	rowWords(A6), D0		; counter1
			EXT.L	D0						; make it long
			MOVE.L	D0, D7					;
			ASL.L	#1, D0					; integers = 2 bytes long -> all off's even
			MOVE.L	D0, D1					;

			SUBQ.L	#2, D1					; D1: rowWords - 1 == col
			MOVE.L	D1, D2					;
			SUBQ.L	#2, D2					; D2: rowWords - 2 == row
			MOVE.L	D1, D3					;
			MOVE.L	D2, D4					;
			MULU.W	D7, D3					; D3: col * N
			MULU.W	D7, D4					; D4: row * N

rowLoop		MOVE.L	D1, D5					;
			ADD.L	D4, D5					; D5: Add1 = row*N + col
			MOVE.L	D2, D6					;
			ADD.L	D3, D6					; D6: Add2 = col*N + row
			
colLoop		MOVE.W	0(A0,D6.L), D7			;
			MOVE.W	0(A0,D5.L), 0(A0,D6.L)	;
			MOVE.W	D7, 0(A0,D5.L)			; swap Add1 & Add2
			
			SUB.L	D0, D6					;
			SUBQ.L	#2, D5					;			
			CMP.L	D5, D6					;
			BNE.S	colLoop					;
			
			SUBQ.L	#2, D2					;
			SUB.L	D0, D4					; 
			BGE.S	rowLoop					;
			
exit		MOVEM.L	(SP)+, UsedRegs			; restore registers
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#6, SP					; pop parameters
			JMP		(A0)					;
			
			DC.B	'TRANSPOS'				; MacsBug Name

			ENDPROC
			
			
ToRCFHT		PROC		EXPORT

; procedure toRCFHT(matPtr: Ptr; 
;				    rowWords: integer);
;
; procedure ToRCFHT takes a matrix of words which has been 1D FHT'd row
; by row and col by col and produces a real 2D FHT from the result
; From the pascal unit FFTnFHT.p, the algorithm is:

;  for row := 0 to maxN div 2 do 			{ Now calculate actual Hartley transform }
;   for col := 0 to maxN div 2 do begin
;     mRow := (maxN - row) mod maxN;
;     mCol := (maxN - col) mod maxN;
;     A := x[row, col];						{ see Bracewell, 'Fast 2D Hartley Transf.' IEEE Procs. 9/86 }
;     B := x[mRow, col];
;     C := x[row, mCol];
;     D := x[mRow, mCol];
;     E := ((A + D) - (B + C)) / 2;
;     x[row, col] := A - E;
;     x[mRow, col] := B + E;
;     x[row, mCol] := C + E;
;     x[mRow, mCol] := D - E;
;    end;

; WARNING: the integer block pointed to by matPtr must have rowWords^2
; words allocated to it. (it must be a square matrix).  No checking for
; this is done.  Furthermore, rowWords must be a power of 2; no checking is
; done for this either.

matPtr		EQU		10						; A6 Offsets
rowWords	EQU		8						;

UsedRegs	REG		D2-D7/A2-A4

			LINK	A6, #0					; no locals (shoulder-hopper!)
			MOVEM.L	UsedRegs, -(SP)			;
			MOVE.L	matPtr(A6), A0			; Base Address
			MOVE.W	rowWords(A6), D0		; counter
			
			MOVE.W	D0, D1					;
			SUBQ.W	#1, D0					;
			SWAP	D0						;
			MOVE.W	D1, D0					; D0: mask | maxN
			MOVE.W	#15, D2					;
logLoop		BTST.L	D2, D1					;
			DBNE.W	D2, logLoop				;
			MOVE.W	D2, D1					; 
			SWAP	D1						; 
			MOVE.W	D0, D1					;
			LSR.W	#1, D1					; D1: log2(maxN) | row 
			MOVE.W	D1, D2					; 
			SWAP	D2						;
			MOVE.W	D1, D2					; D2: col | col
			
rowLoop		MOVEQ.L	#0, D3					;
			MOVEQ.L	#0, D4					;
			MOVEQ.L	#0, D5					;
			MOVEQ.L	#0, D6					; clear the registersÉ
			MOVEQ.L	#0, D7					; for mixed w/l operations
			MOVE.W	D0, D3					;
			MOVE.W	D0, D4					;
			SUB.W	D1, D3					; D3: MaxN - row
			SUB.W	D2, D4					; D4: MaxN - col
			SWAP	D0						;
			AND.W	D0, D3					; D3: (MaxN - row) mod maxN = mRow
			AND.W	D0, D4					; D4: (MaxN - col) mod maxN = mCol
			SWAP	D0						;
			MOVE.W	D1, D5
			SWAP	D1						;
			ASL.L	D1, D3					; D3: maxN * mRow
			ASL.L	D1, D5					; D5: maxN * row
			SWAP	D1						;
			
			MOVE.W	D2, D6					;
			ADD.L	D5, D6					; maxN * row + col
			ASL.L	#1, D6					; integers -> double addr
			MOVEA.L	D6, A1					; A1: (A)
			
			MOVE.W	D2, D7					; 
			ADD.L	D3, D7					; maxN * mRow + col
			ASL.L	#1, D7					; integers take 2 bytes
			MOVEA.L	D7, A2					; A2: (B) 
			
			MOVEQ.L	#0, D6					; make space
			MOVE.W	D4, D6					;
			ADD.L	D5, D6					; maxN * row + mCol
			ASL.L	#1, D6					; integers take 2 bytes
			MOVEA.L	D6, A3					; A3: (C)
			
			MOVEQ.L	#0, D7					; make space
			MOVE.W	D4, D7					;
			ADD.L	D3, D7					; maxN * mRow + mCol
			ASL.L	#1, D7					; even addresses for int's
			MOVEA.L	D7, A4					; A4: (D)
			
			MOVE.W	0(A0, A1.L), D3			; D3: A
			MOVE.W	0(A0, A2.L), D4			; D4: B
			MOVE.W	0(A0, A3.L), D5			; D5: C
			MOVE.W	0(A0, A4.L), D6			; D6: D
			
			MOVE.W	D6, D7					;
			ADD.W	D3, D7					; D7: A + D
			SUB.W	D4, D7					; D7: A+D-B
			SUB.W	D5, D7					; D7: (A+D)-(B+C)
			ASR.W	#1, D7					; D7: E (no Rounding!!)
			BEQ.S	continue				; E=0 -> no changes
			
			SUB.W	D7, D3					; D3: A - E
			ADD.W	D7, D4					; D4: B + E
			ADD.W	D7, D5					; D5: C + E
			SUB.W	D7, D6					; D6: D - E
			
			MOVE.W	D3, 0(A0, A1.L)			;
			MOVE.W	D4, 0(A0, A2.L)			;
			MOVE.W	D5, 0(A0, A3.L)			;
			MOVE.W	D6, 0(A0, A4.L)			; all values stored
			
continue	DBRA	D2, rowLoop				;
			
			MOVE.L	D2, D3					;
			SWAP	D3						;
			MOVE.W	D3, D2					; restore col counter
			
			DBRA	D1, rowLoop				;

exit		MOVEM.L	(SP)+, UsedRegs			; restore registers
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADDQ.L	#6, SP					; pop parameters
			JMP		(A0)					;
			
			DC.B	'TORCFHT '				; MacsBug Name

			ENDPROC
			
			
psFHT2D		PROC		EXPORT

; procedure psFHT2D(srcPtr: ptr;
;					destPtr: ptr;
;					rowWords: integer;
;					scaleType: integer);
;
; function psFHT2D computes the power spectrum of a 2D FHT.
; srcPtr points to a matrix of integers containing the transform.
; destPtr points to a matrix of bytes in which the scaled power spectrum of
; the transform is written. the matrix pointed to by srcPtr is 
; rowWords x rowWords in size.  rowWords must be a power of 2.
;
; psFHT2D first scans the src matrix to determine the largest and
; smallest elements of the power spectrum.  A scaling function
; based on these extremes is developed according to 'scaleType'
; as follows:
; scaleType			scaling implemented:
;		<=0				linear
;		1-9				nth root		(1st root = linear)
;		>=10			log
; Next, the src matrix is scanned again, this time the power 
; spectrum at each point is again computed and scaled to a byte size
; before being stored in the dest matrix.  
; The necessity of scanning the matrix twice was originally avoided by 
; approximating the power spectrum element (a^2 + b^2)/2 by (|a| + |b|)/2, 
; but while this avoided many multiplications, it occasionally produced 
; spurious results.
; Here, some addressing calculations have been avoided by splitting up the
; matrix traversal into four separate sections: two loops and two elements.
; A more compact form of similar addressing may be found above in Hcdc2Bufs.

; register map
; D0		rowWords		A0		srcPtr
; D1		mask			A1		destPtr
; D2		loop counter	A2		longOffset
; D3		scratch			A3		shortOffset
; D4		scratch			
; D5		scratch			
;
; FP0		scratch			
; FP1		mulFactor
; FP2		1/n (for nth root mapping)

srcPtr		EQU		16						; A6 offsets
destPtr		EQU		12						;
rowWords	EQU		10						;
scaleType	EQU		8

UsedRegs	REG		D2-D7/A2-A3				;

start		LINK	A6, #0					; 
			MOVEM.L	UsedRegs, -(SP)			;
			FMOVE.X	FP2, -(SP)				;
			MOVE.L	srcPtr(A6), A0			; A0: srcPtr
			MOVE.W	rowWords(A6), D0		; D0: rowWords
			EXT.L	D0						;
			
			MOVE.W	D0, D1					;
			SUBQ.W	#1, D1					;
			EXT.L	D1						; D1: mod rowWords mask

			MOVE.W	D0, D2					;
			MULS.W	D2, D2					;
			MOVE.L	D2, -(SP)				; for later
			MOVEA.L	D2, A2					;
			ADDA.L	D2, A2					;
			SUBA.L	#2, A2					; A2: longOffset
			
			MOVE.W	D0, D3					;
			EXT.L	D3						;
			MOVEA.L	D3, A3					;
			ADDA.L	D3, A3					;
			ADDA.L	#2, A3					; A3: shortOffset
			
			MOVE	D0, D3					;
			EXT.L	D3						; 
			ASR.L	#1, D3					; D3: rowWords div 2
			ASR.L	#1, D2					; D2: rowWords^2 div 2
			SUB.L	D3, D2					; D2: count
			MOVE.L	#$7FFFFFFF, D6			; D6: minVal
			MOVEQ.L	#0, D7					;
			BSET.L	#31, D7					; D7: maxVal

loop1		MOVE.L	A3, D3					;
			MOVE.L	D3, D4					;
			ASR.L	#1, D4					;
			AND.L	D1, D4					; mod rowWords
			BNE.S	noMod1					;
			MOVE.W	D0, D4					;
			EXT.L	D4						;
			ADD.L	D4, D4					;
			SUB.L	D4, D3					;
noMod1		MOVE.W	0(A0, D3.L), D4			; 
			MOVE.W	0(A0, A2.L), D5			;
			MULS.W	D4, D4					; D4: a^2
			MULS.W	D5, D5					; D5: b^2
			ADD.L	D5, D4					; D4: a^2 + b^2
			CMP.L	D6, D4					;
			BGE.S	notMin1					;
			MOVE.L	D4, D6					;
notMin1		CMP.L	D7, D4					;
			BLE.S	notMax1					;
			MOVE.L	D4, D7					;

notMax1		SUBQ.L	#2, A2					;
			ADDQ.L	#2, A3					;
			SUBQ.L	#1, D2					;
			BGT.S	loop1					;
			
			MOVE.L	(SP), D3				; check middle element, 1st column
			MOVE.W	0(A0, D3.L), D4			;
			MULS.W	D4, D4					;
			ADD.L	D4, D4					;
			CMP.L	D6, D4					;
			BGE.S	notMin2					;
			MOVE.L	D4, D6					;
notMin2		CMP.L	D7, D4					;
			BLE.S	notMax2					;
			MOVE.L	D4, D7					;
			
notMax2		MOVE.W	(A0), D4				; check 0th matrix element
			MULS.W	D4, D4					;
			ADD.L	D4, D4					;
			CMP.L	D6, D4					;
			BGE.S	notMin3					;
			MOVE.L	D4, D6					;
notMin3		CMP.L	D7, D4					;
			BLE.S	notMax3					;
			MOVE.L	D4, D7					;
											; now check top Row
notMax3		MOVE.L	#2, A3					; A3: shortOffset
			MOVE.W	D0, D3					;
			EXT.L	D3						;
			ADD.L	D3, D3					;
			SUBQ.L	#2, D3					;
			MOVE.L	D3, A2					; A2: longOffset
			MOVE.L	D0, D2					;
			ASR.L	#1, D2					; D2: counter			
			
loop2		MOVE.W	0(A0, A2.L), D4			; 
			MOVE.W	0(A0, A3.L), D5			;
			MULS.W	D4, D4					;
			MULS.W	D5, D5					;
			ADD.L	D5, D4					; D4: a^2 + b^2
			CMP.L	D6, D4					;
			BGE.S	notMin4					;
			MOVE.L	D4, D6					;
notMin4		CMP.L	D7, D4					;
			BLE.S	notMax4					;
			MOVE.L	D4, D7					;
			
notMax4		SUBA.L	#2, A2					;
			ADDA.W	#2, A3					;
			SUBQ.L	#1, D2					;
			BGT.S	loop2					; done checking at last!
			
			SUB.L	D6, D7					; D7: range
			BNE.S	findScale				;
			MOVEQ.L	#1, D7					;

findScale	MOVE.W	scaleType(A6), D4		; 
			BLE.S	linMap1					; scaleType <= 0 -> linear Mapping
			CMP.W	#10, D4					;
			BGE.S	logMap1					; scaleType >= 10 -> log Mapping

rootMap1	FMOVE.L	#1, FP2					; 0 < scaleType < 10 -> nth root mapping
			FDIV.W	D4, FP2					; FP2: 1/n
			FLOGN.L	D7, FP0					; FP0: ln(range)
			FMUL.X	FP2, FP0				; FP0: 1/n * ln(range)
			FETOX.X	FP0, FP0				; FP0: e^(1/n * ln(range)) = range^(1/n)
			MOVE.W	#2, D5					; D5: ScaleType = 2 -> nth root
			BRA.S	GoOn1					;
			
linMap1		FMOVE.L	D7, FP0					;
			MOVE.W	#1, D5					; D5: ScaleType = 1 -> linear
			BRA.S	GoOn1					;
			
logMap1		SUBQ.L	#1, D6					; min := min - 1 ; to prevent log(0)
			ADDQ.L	#1, D7					; which increases range by 1
			FLOG2.L	D7, FP0					;
			MOVE.W	#3, D5					; D5: ScaleType = 3 -> log
			
GoOn1		FMOVE.L	#255, FP1				;
			FDIV.X	FP0, FP1				; FP1: mul factor
			MOVE.L	D5, D7					; D7: Copy of ScaleType
											; D6: minVal
			
doPS		MOVE.L	destPtr(A6), A1			; compute body of PS
			MOVE.L	(SP), D2				;
			MOVEA.L	D2, A2					;
			ADDA.L	D2, A2					;
			SUBA.L	#2, A2					; A2: longOffset
			
			MOVE.W	D0, D3					;
			EXT.L	D3						;
			MOVEA.L	D3, A3					;
			ADDA.L	D3, A3					;
			ADDA.W	#2, A3					; A3: shortOffset
			
			MOVE	D0, D3					;
			EXT.L	D3						; 
			ASR.L	#1, D3					; D3: rowWords div 2
			ASR.L	#1, D2					; D2: rowWords^2 div 2
			SUB.L	D3, D2					; D2: count

loop3		MOVE.L	A3, D3					;
			MOVE.L	D3, D4					;
			ASR.L	#1, D4					;
			AND.L	D1, D4					; mod rowWords
			BNE.S	noMod2					;
			MOVE.W	D0, D4					;
			EXT.L	D4						;
			ADD.L	D4, D4					;
			SUB.L	D4, D3					;
noMod2		MOVE.W	0(A0, D3.L), D4			; 
			MOVE.W	0(A0, A2.L), D5			;
			MULS.W	D4, D4					;
			MULS.W	D5, D5					;
			ADD.L	D5, D4					; D4: x^2 + y^2
			SUB.L	D6, D4					; D4: " - minVal
			MOVE.W	D7, D5					; D5: scaleType
			SUBQ.W	#1, D5					;
			BEQ.S	linMap2					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap2				;
			
logMap2		FLOG2.L	D4, FP0					; take log
			FMUL.X	FP1, FP0				; scale up to [0..255]
			BRA.S	GoOn2					;

linMap2		FMOVE.L	D4, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn2					;
			
rootMap2	FLOGN.L	D4, FP0					;
			FMUL.X	FP2, FP0				;
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn2		FMOVE.L	FP0, D4					; D4: scaled out: should be in [0,255]
			MOVE.L	D3, D5					;
			ASR.L	#1, D5					;
			MOVE.B	D4, 0(A1, D5.L)			; !!!doesn't work on 68000
			MOVE.L	A2, D5					; 
			ASR.L	#1, D5					;
			MOVE.B	D4, 0(A1, D5.L)			; stored result bytes
			
			SUBQ.L	#2, A2					;
			ADDQ.L	#2, A3					;
			SUBQ.L	#1, D2					;
			BGT.S	loop3					;
			
			MOVE.L	(SP)+, D3				; compute middle element, 1st column
			MOVE.W	0(A0, D3.L), D4			;
			MULS.W	D4, D4					; D4: (x^2 + y^2) div 2
			ADD.L	D4, D4					;
			SUB.L	D6, D4					; D4: " - minVal
			MOVE.W	D7, D5					; D6: scaleType
			SUBQ.W	#1, D5					;
			BEQ.S	linMap3					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap3				;
			
logMap3		FLOG2.L	D4, FP0					; take log
			FMUL.X	FP1, FP0				; scale up to [0..255]
			BRA.S	GoOn3					;

linMap3		FMOVE.L	D4, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn3					;
			
rootMap3	FLOGN.L	D4, FP0					;
			FMUL.X	FP2, FP0				;
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn3		FMOVE.L	FP0, D4					; D4: scaled out: should be in [0,255]
			MOVE.L	D3, D5					;
			ASR.L	#1, D5					;
			MOVE.B	D4, 0(A1, D5.L)			; !!!doesn't work on 68000
			
			MOVE.W	(A0), D4				; compute 0th matrix element
			MULS.W	D4, D4					; D4: (x^2 + y^2) div 2
			ADD.L	D4, D4
			SUB.L	D6, D4					; D4: " - minVal
			MOVE.W	D7, D5					; D6: scaleType
			SUBQ.W	#1, D5					;
			BEQ.S	linMap4					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap4				;
			
logMap4		FLOG2.L	D4, FP0					; take log
			FMUL.X	FP1, FP0				; scale up to [0..255]
			BRA.S	GoOn4					;

linMap4		FMOVE.L	D4, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn4					;
			
rootMap4	FLOGN.L	D4, FP0					;
			FMUL.X	FP2, FP0				;
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn4		FMOVE.L	FP0, D4					; D4: scaled out: should be in [0,255]
			MOVE.B	D4, (A1)				; !!!doesn't work on 68000
			
											; now compute top Row PS
			MOVE.L	#2, A3					; A3: shortOffset
			MOVE.W	D0, D3					;
			EXT.L	D3						;
			ADD.L	D3, D3					;
			SUBQ.L	#2, D3					;
			MOVE.L	D3, A2					; A2: longOffset
			MOVE.L	D0, D2					;
			SUBQ.L	#1, D2					; D2: counter
			
loop4		MOVE.W	0(A0, A2.L), D4			; 
			MOVE.W	0(A0, A3.L), D5			;
			MULS.W	D4, D4					;
			MULS.W	D5, D5					;
			ADD.L	D5, D4					; 
			SUB.L	D6, D4					; D4: " - minVal
			MOVE.W	D7, D5					; D6: scaleType
			SUBQ.W	#1, D5					;
			BEQ.S	linMap5					;
			SUBQ.W	#1, D5					;
			BEQ.S	rootMap5				;
			
logMap5		FLOG2.L	D4, FP0					; take log
			FMUL.X	FP1, FP0				; scale up to [0..255]
			BRA.S	GoOn5					;

linMap5		FMOVE.L	D4, FP0					;
			FMUL.X	FP1, FP0				;
			BRA.S	GoOn5					;
			
rootMap5	FLOGN.L	D4, FP0					;
			FMUL.X	FP2, FP0				;
			FETOX.X	FP0, FP0				;
			FMUL.X	FP1, FP0				;

GoOn5		FMOVE.L	FP0, D4					; D4: scaled out: should be in [0,255]
			MOVE.L	D3, D5					;
			ASR.L	#1, D5					;
			MOVE.B	D4, 0(A1, D5.L)			; !!!doesn't work on 68000
			MOVE.L	A2, D5					; 
			ASR.L	#1, D5					;
			MOVE.B	D4, 0(A1, D5.L)			; stored result bytes
			
			SUBA.L	#2, A2					;
			ADDA.W	#2, A3					;
			SUBQ.L	#1, D2					;
			BGT.S	loop4					; done at last!

exit		FMOVE.X	(SP)+, FP2				; restore FP2
			MOVEM.L	(SP)+, UsedRegs			;
			UNLK	A6						;
			MOVE.L	(SP)+, A0				;
			ADD.L	#12, SP					;
			JMP		(A0)					;
			
			DC.B	'PSFHT2D '				; MacsBug Name

			ENDPROC

			END