*****************************************************
* fft.c11 - fast fourier transform for the MC68HC11
*        written by:
*
*        Ron Williams
*        Department of Chemistry
*        Ohio University
*        Athens, OH 45701
*
* This is a modification of the 6800 FFT presented by:
*      Richard Lord
*      Byte Magazine, pp. 108-119
*      February 1979
*
* My version is written in ROMable code for the HC11. 
* It uses a sine look-up table for speed and can only 
* transform 256 8-bit data points.  The program 
* assumes that the address of the real data is pushed 
* on the stack prior to the call and that a 256 byte 
* imaginary buffer is at data+256 therefore you must 
* declare a 512 byte data array in the calling routine 
* and load the lower 256 bytes with data.  The FFT 
* will zero out the imaginary portion.  Also note that 
* the FFT uses memory in the stack RAM for its dynamic 
* variables and the FFT returns a value on the stack 
* which contains the number of times the data was 
* divided by 2 during transform. 
* 
* As mentioned in Lord's article, "power" spectra can 
* be computed using the sum of absolute values routine 
* included at the end of the FFT. Simply change the 
* beq done at the end of the FFT to beq smsq. 
* 
* note - this copy has been modified to use $DD00 for 
* data because this is easiest with BUFFALO
* I have timed this transform on some test data.  The 
* results are an impressive 350 milliseconds per 
* transform including the "power" spectra computation. 
* 
* Please note that the origin of this code may require 
* adjustment for your specific memory map. 
* 
* Please let me know of any bugs you find.

RETURN  EQU     0
REAL    EQU     2
CELNM   EQU     4
CELCT   EQU     5
PAIRNM  EQU     6
CELDIS  EQU     7
DELTA   EQU     8
SCLFCT  EQU     9
COSA    EQU     $0A
SINA    EQU     $0B
SINPT   EQU     $0C
REAL1   EQU     $0E
REAL2   EQU     $10
TREAL   EQU     $12
TIMAG   EQU     $13
TMP     EQU     $14
TMP2    EQU     $15
DATA    EQU     $DD00
 
        ORG     $C000

        TSX             ; top of stack for frame pointer
        XGDX            ;     to be placed in x
        SUBD    #$18    ; subtract offset to make room
        XGDX            ; x now has frame pointer
        PULY            ; get return address
        STY     RETURN,X ; save it
*        PULY            ; get data address
*        STY     REAL,X  ; save it
        LDY    #DATA
        STY    REAL,X
        CLR     SCLFCT,X ; zero scale factor
        INY             ; inc y for imag data
        CLRB
ZERO    CLR     $FF,Y   ; note special place of imag
        INY             ;  256 above data
        DECB
        BNE     ZERO
*
* MUST DO BIT SORTING BEFORE TRANSFORMING
*
        LDAB    #$FE    ; setup start for bit reversal
REVBIT  LDAA    #08     ; get # of bits to reverse
        PSHB            ; save address offset
REV1    RORB            ; rotate b right - bit to carry
        ROL     TMP,X   ; rotate left - carry bit in
        DECA            ; decrement counter
        BNE     REV1    ; go back if not done
        PULB            ; get unshifted address
        PSHB            ; save copy
        CMPB    TMP,X   ; check to see if already done
        BCS     NOSWAP  ; if so don't swap bytes
SWAP    LDY     REAL,X  ; get data address
        ABY             ; add to base address
        LDAA    0,Y     ; get value
        PSHY            ; store away
        LDY     REAL,X  ; get base again
        LDAB    TMP,X   ; get shifted address
        ABY             ; add to base
        LDAB    0,Y     ; get second member
        STAA    0,Y     ; put away first member
        PULY            ; get first address
        STAB    0,Y     ; put second member in first slot
NOSWAP  PULB            ; get current address back
        DECB            ; decrement it
        BNE     REVBIT  ; do next if not done
*
* SPECIAL CASE OF FIRST PASS OF FFT
*
        JSR     SCALE
        LDY     REAL,X  ; set up data pointer
        LDAA    #128    ; get number of cells
        STAA    TMP,X   ; store in temp
FPSS    LDAA    0,Y     ; get rm
        LDAB    1,Y     ; get rn
        PSHA            ; make copy
        ABA             ; rm'=rm+rn
        STAA    0,Y     ; save back in data array
        PULA            ; get rm again
        SBA             ; rn'=rm-rn
        STAA    1,Y     ; put away
        INY             ; point to next pair
        INY
        DEC     TMP,X   ; decrement # cells
        BNE     FPSS    ; go back if not done
*
* NOW THE FFT PROPER FOR PASSES 2 THRU N
*
FOUR    LDAA    #64     ; # of cells is now 64
        STAA    CELNM,X ; store
        STAA    DELTA,X ; so is delta
        LDAA    #02     ; number of pairs is 2
        STAA    PAIRNM,X
        STAA    CELDIS,X ; so is distance between
NPASS   JSR     SCALE   ; check for over-range
        LDAA    CELNM,X ; get current cell #
        STAA    CELCT,X ; store at cell counter
        LDY     REAL,X
        STY     REAL1,X ; get copy of data
NCELL   LDY     #SINTAB ; get address of sines
        STY     SINPT,X ; save copy
        LDAA    PAIRNM,X ; get current pairnm
NP1     PSHA            ; save pair counter
        LDAA    0,Y     ; get cosine
        LDAB    64,Y    ; get sine
        STAA    COSA,X  ; save copy
        STAB    SINA,X  ; ditto
        LDY     REAL1,X ; point to top of data
        LDAB    CELDIS,X ; get current offset
        ABY             ; add to y for current 
        STY     REAL2,X ; copy it
        LDAA    0,Y     ; get data point rn
        PSHA            ; copy it
        LDAB    COSA,X  ; get cosine
        JSR     SMUL    ; rn*cos(a)
        STAA    TREAL,X
        PULA            ; get copy of rn
        LDAB    SINA,X  ; get sin(a)
        JSR     SMUL    ; rn*sin(a)
        STAA    TIMAG,X ; store imaginary tmp
        INY
        LDAA    $FF,Y   ; get imaginary data
        PSHA            ; save it
        LDAB    SINA,X  ; get sin(a)
        JSR     SMUL    ; in*sin(a)
        ADDA    TREAL,X
        STAA    TREAL,X  ; tr=rn*cos + in*sin
        PULA            ; get data back
        LDAB    COSA,X  ; get cosine
        JSR     SMUL    ; in*cos(a)
        SUBA    TIMAG,X  ; ti=in*cos-rn*sin
        STAA    TIMAG,X
        LDY     REAL1,X
        LDAA    00,Y    ; get rm 
        TAB             ; save a copy
        ADDA    TREAL,X ; rm'=rm+tr
        STAA    00,Y    ; store new rm
        SUBB    TREAL,X ; rn'=rm-tr
        LDY     REAL2,X
        STAB    00,Y    ; store new rn
        LDY     REAL1,X
        INY
        STY     REAL1,X ; save real1 for nxt
        LDAA    $FF,Y   ; get im
        TAB             ; save copy
        ADDA    TIMAG,X ; im'=im+ti
        STAA    $FF,Y   ; put back in array
        LDY     REAL2,X
        INY
        SUBB    TIMAG,X ; in'=im-ti
        STAB    $FF,Y   ; put back in array
        LDY     SINPT,X
        LDAB    DELTA,X ; increment sine pntr
        ABY
        STY     SINPT,X ; save away
        PULA
        DECA            ; dec pair counter
        BNE     NP1
AR1     LDY     REAL1,X
        LDAB    CELDIS,X
        ABY
        STY     REAL1,X
        DEC     CELCT,X
        BEQ     AR3
        JMP     NCELL
AR3     LSR     CELNM,X  ; half cells
        BEQ     SMSQ    ; done when all cells
        ASL     PAIRNM,X ; double pairs
        ASL     CELDIS,X ; twice as far apart
        LSR     DELTA,X  ; delta is half
        JMP     NPASS    ; one more time!
DONE    LDAA    SCLFCT,X ; get scale factor
        PSHA             ; save on stack
        LDY     RETURN,X
        PSHY
        RTS
*
* SUM OF ABSOLUTE VALUES INSTEAD OF SUM OF SQUARES
*
SMSQ    LDY     REAL,X   ; compute sum of "sqrs"
        CLRA             ; clear byte counter
SUM     PSHA             ; save on stack
        LDAA    0,Y      ; get real data point
        BPL     SM1      ; force positive
        NEGA
        BVC     SM1      ; watch for $80
        CLRA             ;   which is really 0
SM1     INY              ; get imaginary data
        LDAB    $FF,Y
        BPL     SM2      ; force positive
        NEGB
        BVC     SM2      ; watch for $80 again
        CLRB
SM2     DEY              ; correct data pointer
        ABA              ; compute sum
        STAA    0,Y      ; save back in real
        INY              ; inc y for next round
        PULA             ; get byte counter
        DECA             ; done when zero
        BNE     SUM
        BRA     DONE     ; let's get out of here
*
* SUBROUTINE FOR CATCHING OVERSCALED DATA
*
SCALE   LDY     REAL,X   ; start at top of data
        LDAB    #$FF
        ABY              ; top of data
        ABY              ; top or imag
        INY              ; need two more
        INY
        LDAA    #$C0     ; -64
        LDAB    #$40     ; +64
TOP     CMPA    0,Y      ; check for minimum
        BLO     NXT      ; if more negative fix
        CMPB    0,Y      ; check for too big
        BCS     SCL      ; go fix it
NXT     DEY              ; bump pointer
        CPY  REAL,X      ; done when both
        BNE     TOP      ; imag and data done
        RTS
SCL     INC     SCLFCT,X ; keep track of scale
        LDY     REAL,X   ; set up pointer
        LDAB    #$FF
        ABY
        ABY
        INY
        INY
SCL1    LDAA    0,Y      ; get data
        ADDA    #$80     ; make positive
        LSRA             ; divide by two
        SUBA    #$40     ; put back
        STAA    0,Y      ; store away
        DEY              ; bump pointer
        CPY     REAL,X   ; done when both
        BNE     SCL1     ; imag and data done
        RTS
*
* THE HC11 MULTIPLY MUST BE MODIFIED TO HANDLE
* NEGATIVE DATA
*
SMUL    STAA    TMP,X   ; copy multiplier
        STAB    TMP2,X  ; ditto multiplicand
        TSTA            ; check sign of multiplier
        BPL     SK1     ; skip negation
        NEGA
        BVS     SKO     ; check for $80
        BEQ     SKO     ; check for zero
SK1     TSTB            ; check multiplier sign
        BPL     SK2
        NEGB
        BVS     SKO     ; check for $80
        BEQ     SKO
SK2     MUL             ; do multiplication
        ADCA    #0      ; 8 bit conversion
        ASLA            ; and correct for sine
        LDAB    TMP2,X  ; get original multiplicand
        EORB    TMP,X   ; check for result
        BPL     OUT
        NEGA            ; result is negative
OUT     RTS
SKO     CLRA            ; return zero to main
        RTS
*
* NOW FOR THE SINE LOOK UP TABLE
*
SINTAB            
 FCB  127, 127, 127, 127, 126, 126, 126, 125, 125, 124
 FCB  123, 122, 122, 121, 120, 118, 117, 116, 115, 113
 FCB  112, 111, 109, 107, 106, 104, 102, 100,  98,  96
 FCB   94,  92,  90,  88,  85,  83,  81,  78,  76,  73
 FCB   71,  68,  65,  63,  60,  57,  54,  51,  49,  46
 FCB   43,  40,  37,  34,  31,  28,  25,  22,  19,  16
 FCB   12,   9,   6,   3,   0,  -3,  -6,  -9, -12, -16
 FCB  -19, -22, -25, -28, -31, -34, -37, -40, -43, -46
 FCB  -49, -51, -54, -57, -60, -63, -65, -68, -71, -73
 FCB  -76, -78, -81, -83, -85, -88, -90, -92, -94, -96
 FCB  -98,-100,-102,-104,-106,-107,-109,-111,-112,-113
 FCB -115,-116,-117,-118,-120,-121,-122,-122,-123,-124
 FCB -125,-125,-126,-126,-126,-127,-127,-127,-127,-127
 FCB -127,-127,-126,-126,-126,-125,-125,-124,-123,-122
 FCB -122,-121,-120,-118,-117,-116,-115,-113,-112,-111
 FCB -109,-107,-106,-104,-102,-100, -98, -96, -94, -92
 FCB  -90, -88, -85, -83, -81, -78, -76, -73, -71, -68
 FCB  -65, -63, -60, -57, -54, -51, -49, -46, -43, -40
 FCB  -37, -34, -31, -28, -25, -22, -19, -16, -12,  -9
 FCB   -6,  -3,   0,   3,   6,   9,  12,  16,  19,  22
 FCB   25,  28,  31,  34,  37,  40,  43,  46,  49,  51
 FCB   54,  57,  60,  63,  65,  68,  71,  73,  76,  78
 FCB   81,  83,  85,  88,  90,  92,  94,  96,  98, 100
 FCB  102, 104, 106, 107, 109, 111, 112, 113, 115, 116
 FCB  117, 118, 120, 121, 122, 122, 123, 124, 125, 125
 FCB  126, 126, 126, 127, 127, 127
*
* SOME ROUTINES FOR FAST A2D AND FOR SENDING DATA TO 
* HOST
*   THESE HAVE NOT BEEN DEBUGGED
*
CHAN4   EQU     $04
ATDCTR  EQU     $1030
PORTC   EQU     $1003
CCONTR  EQU     $1002
LPORTC  EQU     $1005
STRT    EQU     $02
        LDY     #DATA
        LDAA    #CHAN4
        LDAB    #STRT
WAIT:   BITB    PORTC
        BEQ     WAIT
CKLP:   LDAB    CCONTR
        BPL     CKLP
        STAA    ATDCTR
        LDAB    LPORTC
        INY
        CPY     #DATA+256
        BEQ     ADONE
        NOP
        NOP
        NOP
        NOP
        LDAB    ATDCTR+1
        STAB    $00,Y
        BRA     CKLP
ADONE:  RTS

OUTPUT  EQU    $E3B3
SENDAT  LDX    #DATA
        CLRB
S1      JSR    OUTPUT
        DECB
        BNE    S1
        RTS
