; A MICROSOFT-C PROCEDURE         R.Tervo    March 20, 1993
; Copyright (c) 1985,1997 SCCON Scientific Research Inc
; For information, contact  < tervo@sccon.ca >
;*************************************************************
; adapted from....
; A TURBO PASCAL PROCEDURE        R.Tervo  December 1, 1985
;*************************************************************
; 02 FEB 98 (rt) changed label to CONSTx (assembler errors)
;*************************************************************
; CALLING CONVENTION:
;
; fft( vector, N, direction )
; struct complex *vector;
; int N, direction
; *** USES 8087 ***
;-------------------------------------------------------------
; DESCRIPTION:
;
; Calculate Fourier transform of input vector xi,xr
; N is log(base2) of the number of points in the vector
; IF direction = 0 THEN do forward FFT, ELSE do reverse FFT
;-------------------------------------------------------------
; ASSUMES:  (Standard math.h definition of type complex)
;
; struct complex { double x;
;                  double y;
;                };
;-------------------------------------------------------------
 COMPLEX STRUC   ; define corresponding assembler structure
  R DQ ?         ; 8087 longreals are 8 bytes each (QuadWord)
  I DQ ?
 COMPLEX ENDS
;-------------------------------------------------------------
	TITLE   fft.asm
	NAME    fft

_TEXT   SEGMENT  WORD PUBLIC 'CODE'
_TEXT   ENDS

_DATA   SEGMENT  WORD PUBLIC 'DATA'
_DATA   ENDS

CONSTx  SEGMENT  WORD PUBLIC 'CONST'
CONSTx  ENDS

_BSS    SEGMENT  WORD PUBLIC 'BSS'
_BSS    ENDS

DGROUP  GROUP  CONSTx, _BSS, _DATA
	ASSUME  CS: _TEXT, DS: DGROUP, SS: DGROUP

; EXTRN   __chkstk : NEAR

_TEXT   SEGMENT
	ASSUME  CS: _TEXT

	PUBLIC _fft5
_fft5   PROC NEAR

	PUSH BP
	MOV  BP,SP
	SUB  SP,50
	PUSH DS
	PUSH SI
	CALL begin
	POP  SI
	POP  DS
	MOV  SP,BP
	POP  BP
	RET
;--------------------------------------------------------------
; The C-language convention PUSHes variables from right-to-left
; FFT3.ASM does not need a segment register for the data buffer...
;
xone   EQU  word ptr [BP+ 4] ; base address of input array
N      EQU  word ptr [BP+ 6] ; log_2 the number of values
dir    EQU  word ptr [BP+ 8] ; forward or reverse transform
mag    EQU  word ptr [BP+10] ; # of points to do magnitude with
y0     EQU  word ptr [BP+12] ; centreline of output plot
ymax   EQU  word ptr [BP+14] ; range of output plot
buf    EQU  word ptr [BP+16] ; base address of output plotting array

; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Xj1    EQU  word ptr [BP-4]  ; offset pointer to X[j1]
Xj2    EQU  word ptr [BP-6]  ; offset pointer to X[j2]

status EQU  word ptr [BP-8]  ; used for trig calculation
sign   EQU  word ptr [BP-10] ; keep sign of input angle
four   EQU  word ptr [BP-12] ; (constant = 4)

NP     EQU  word ptr [BP-14] ; NP = Number of Points = 2 to the N
L      EQU  word ptr [BP-16]
MEX2   EQU  word ptr [BP-18]
NMEX2  EQU  word ptr [BP-20]
HALF   EQU  word ptr [BP-22]
MP1    EQU  word ptr [BP-24]
j0     EQU  word ptr [BP-26]
j      EQU  word ptr [BP-28]
CH2PI  EQU qword ptr [BP-36] ; (constant = 2 * pi / NP) (-ve for Forward FFT)

rvalue EQU qword ptr [DI].R  ; implied DS:[DI]
ivalue EQU qword ptr [DI].I  ; implied DS:[DI+8]

rvalu1 EQU qword ptr [SI].R  ; implied DS:[SI]
ivalu1 EQU qword ptr [SI].I  ; implied DS:[SI+8]

rvalu2 EQU qword ptr [DI].R  ; implied DS:[DI]
ivalu2 EQU qword ptr [DI].I  ; implied DS:[DI+8]

magx   EQU qword ptr [DI].R  ;
bufy   EQU  word ptr [SI]    ; pointer to plotting buffer
;---------------------- main program ----------------------------

begin:  MOV   CX,N
	MOV   AX,1
	SHL   AX,CL
	MOV   NP,AX      ; NP := 2^N   (note: NP = FNP)
	
	FINIT            ;--0------1------
	FILD  NP         ; FNP
	FLDPI            ;  pi     FNP
	FADD  ST,ST(0)   ; 2pi     FNP
	FDIVP ST(1),ST   ; 2pi/FNP
	CMP   dir,0
	JNZ   noneg
	FCHS             ; -ve for forward FFT
 noneg: FSTP  ch2pi      ; store result as constant
	MOV   four,4     ; store another constant
;;      MOV   DS,xseg    ; DS points to input array

;------------------ FOR MP1:= 1 to N --------------------------------

	MOV   MP1,1      ; MP1 counts this loop

MP1top: MOV   CX,MP1
	DEC   CX         ; CX = (MP1-1)
	MOV   AX,1
	SHL   AX,CL
	MOV   MEX2,AX    ; MEX2 := 2 to the (MP1-1)
	MOV   AX,N
	SUB   AX,CX
	MOV   CX,AX
	MOV   AX,1       ; NMEX2 := NP  div MEX2
	SHL   AX,CL      ;        = 2^N div 2^(MP-1)
	MOV   NMEX2,AX   ;        = 2^( N - (MP-1) )
	SHR   AX,1
	MOV   HALF,AX    ; HALF := NMEX2 div 2
	MOV   L ,0
	MOV   j0,0

;------------- FOR  j00:= 1 to MEX2 ------------------------------

	MOV   DX,MEX2    ; DX counts this loop  (j00 is a dummy)

j00top: MOV   AX,xone   ; get pointer to X[array]
	MOV   Xj1,AX
	MOV   Xj2,AX

	MOV   AX,j0      ; get J0
	MOV   CL,4
	SHL   AX,CL      ; j1 := J0  (multiplied by 16)
	ADD   Xj1,AX
	MOV   BX,half    ; get half
	SHL   BX,CL
	ADD   AX,BX      ; j2 := j1 + half  (multiplied by 16)
	ADD   Xj2,AX

;---------------- calculate sinV , cosV -------------------------
;
QT100  EQU 0100h  ;C0       ; ** STATUS WORD BIT MASKS **
QT010  EQU 4000h  ;C3       ; C0 C3 C1 = HalfQuadrant of theta
QT001  EQU 0200h  ;C1       ;          = QT in 4,2,1 order...
QT011  EQU 4200h            ; ...as set by FPREM instruction
QT110  EQU 4100h            ;
C0     EQU 0100h            ; sign bit from FTST instruction

	FINIT               ;--0-------1-------2-------3-
	FILD   four         ;  4
	FLDPI               ;  pi      4
	FDIVP  ST(1),ST     ; pi/4
	FILD   L            ;  L      pi/4
	FLD    ch2pi        ; ch2pi    L      pi/4
	FMULP  ST(1),ST     ; theta   pi/4
	FTST                ;
	FSTSW  sign         ; save original sign of theta for later...
	FABS                ; ...and work only on positive angles
	FPREM               ; theta  pi/4  (scaling: 0 <= theta <= pi/4 )
	FSTSW  status       ; theta  pi/4
	TEST   status,QT001 ; IF odd HalfQuadrant
	JZ     ok           ; THEN theta :=  pi/4 - theta;
	FSUBR  ST,ST(1)     ;
   ok:  FPTAN               ;  X       Y
	FLD    ST           ;  X       X       Y
	FMUL   ST,ST(0)     ; X*X      X       Y
	FLD    ST(2)        ;  Y      X*X      X       Y
	FMUL   ST,ST(0)     ; Y*Y     X*X      X       Y
	FADDP  ST(1),ST     ; Y*Y+X*X  X       Y
	FSQRT               ; hyp      X       Y
	FDIVR  ST(2),ST     ; hyp      X      sin
	FDIVRP ST(1),ST     ; cos     sin    (pi/4)
	FFREE  ST(2)        ; NOTE: these are still 'raw' values
	MOV    AX,status
	TEST   AX,QT011
	JZ     noswap
	XOR    AX,QT011
	TEST   AX,QT011     ; true sin and cos are swapped
	JZ     noswap       ; when in HalfQuadrants 1,2,5,6
	FXCH   ST(1)        ; cos     sin
noswap: MOV    AX,status
	TEST   AX,QT110
	JZ     cosign
	XOR    AX,QT110     ; COSIGN(x) FIXUP
	TEST   AX,QT110     ; cosign is -ve
	JZ     cosign       ; when in HalfQuadrants 2,3,4,5
	FCHS                ; cos     sin
cosign: FXCH   ST(1)        ; sin     cos
	TEST   status,QT100
	JZ     cont         ; SIGN(x) FIXUP
	FCHS                ; sine is -ve when in 4,5,6,7
 cont:  TEST   sign,C0
	JZ     suite        ; check original sign of theta
	FCHS                ; sin(-x) := -sin(x)

;---------------- FOR i:= 1 to HALF ------------------------------
; this is a critical loop, done N * 2^(N-1) times = 2304 for (N=9)
;
suite:  MOV   SI,Xj1        ; CX counts this loop
	MOV   DI,Xj2        ;--0------1------2------3------4------5---
	MOV   CX,half       ; sinv   cosv

tight:  FLD   ST(1)         ; cosv   sinv   cosv
	FLD   ST(1)         ; sinv   cosv   sinv   cosv
	FLD   rvalu2        ; XRj2   sinv   cosv   sinv   cosv
	FMUL  ST(1),ST
	FMUL  ST,ST(2)      ; XRcos  XRsin  cosv   sinv   cosv
	FLD   ivalu2        ; XIj2   XRcos  XRsin  cosv   sinv   cosv
	FMUL  ST(3),ST
	FMUL  ST,ST(4)      ; XIsin  XRcos  XRsin  XIcos  sinv   cosv
	FSUBRP ST(1),ST     ; Q11    XRsin  XIcos  sinv   cosv
	FXCH  ST(2)         ; XIcos  XRsin  Q11    sinv   cosv
	FADDP ST(1),ST      ; Q12    Q11    sinv   cosv
			    ;--0------1------2------3------4------5---
	FLD   ST            ; Q12    Q12    Q11    sinv   cosv
	FLD   ivalu1        ; XIj1   Q12    Q12    Q11    sinv   cosv
	FSUB  ST(2),ST
	FADDP ST(1),ST
	FSTP  ivalu1        ;                          XIj1 := XIj1 + Q12
	FSTP  ivalu2        ; Q11                      XIj2 := XIj1 - Q12
			    ;--0------1------2------3------4------5---
	FLD   ST            ; Q11    Q11    sinv   cosv
	FLD   rvalu1        ; XRj1   Q11    Q11
	FSUB  ST(2),ST
	FADDP ST(1),ST
	FSTP  rvalu1        ;                          XRj1 := XRj1 + Q11
	FSTP  rvalu2        ; sinv   cosv              XRj2 := XRj1 - Q11

	ADD   SI,16         ; point to next value (16 bytes higher)...
	ADD   DI,16
	LOOP  tight
;----------------------- end {for i: 1 to half} ------------------------

done:   MOV   AX,NMEX2
	ADD   j0,AX           ; j0 := j0 + NMEX2

	MOV   AX,NP           ; UPDATE VALUE FOR L
	SHR   AX,1
	SHR   AX,1            ; AX := IOX[2]
	MOV   BX,L
	MOV   CX,N
	DEC   CX              ; FOR i:= 2 to N DO (CX is loop counter)
around: CMP   AX,BX
	JG    newel           ; AX := IOX[i]
	SUB   BX,AX           ; L  := L - IOX[i]
	SHR   AX,1            ; i  := i + 1
	LOOP around

newel:  ADD   BX,AX
	MOV   L,BX            ; L has been updated

	DEC   DX
	JZ    fin1
	JMP   j00top
;-------------------- end { FOR j00:= 1 to MEX2 } ------------------

fin1:   INC   MP1
	MOV   AX,MP1
	CMP   N,AX
	JS    fin2
	JMP   MP1top
;-------------------- end { FOR MP1:= 1 to N }----------------------

;---------------------- FOR J:=1 to NP -----------------------------
fin2:   MOV   L,0
	MOV   J,1             ; J is used directly as loop counter

jtop:   MOV   AX,J
	CMP   L,AX
	JS    VARYL           ; IF L < J THEN change L

	DEC   AX              ; (XR[1] is array element zero)
	MOV   CL,4            ; ELSE  swap    XR[J] <---> XR[L+1]
	SHL   AX,CL           ;               XI[J] <---> XI[L+1]
	ADD   AX,xone
	MOV   DI,AX           ; [DI] points to X[J]
	MOV   AX,L
	SHL   AX,CL
	ADD   AX,xone
	MOV   SI,AX           ; [SI] points to X[L+1]
	MOV   CX,8            ; sweep through real and imag parts
 swapr: MOV   AX,word ptr[SI]
	XCHG  AX,word ptr[DI]
	MOV   word ptr[SI],AX
	ADD   DI,2
	ADD   SI,2
	loop  swapr

 VARYL: MOV   AX,NP           ; UPDATE VALUE FOR L
	SHR   AX,1            ; AX := IOX[1]
	MOV   BX,L
	MOV   CX,N            ; FOR i:= 1 to N DO (CX is loop counter)
 about: CMP   AX,BX           ; IF ( IOX[i] > L ) THEN stop looping
	JG    newL            ; AX := IOX[i]
	SUB   BX,AX           ; L  := L - IOX[i]
	SHR   AX,1            ; i  := i + 1
	LOOP  about

 newl:  ADD   BX,AX
	MOV   L,BX            ; L has been updated

	INC   J
	MOV   AX,NP
	CMP   AX,J
	JS    final
	JMP   jtop

;------------ IF forward FFT THEN {scale output vector} --------------
final:  CMP   dir,0
	JNZ   normal

	MOV   DI,xone       ; X[i] := X[i]/FNP
	MOV   CX,NP         ; CX counts the loop

	FINIT               ;--0------1----
	FILD  NP            ; FNP
scalex: FLD   ivalue        ; XI      FNP
	FDIV  ST,ST(1)      ; XI/FNP  FNP
	FSTP  ivalue        ; FNP
	FLD   rvalue        ; XR      FNP
	FDIV  ST,ST(1)      ; XR/FNP  FNP
	FSTP  rvalue        ; FNP
	ADD   DI,16
	LOOP  scalex

;--------------------- end {IF forward } -------------------
normal: FINIT
	RET            ; that's all folks

_fft5   ENDP
_TEXT   ENDS
	END