Here's a PIC routine that will take the sine of a 14-bit frequency variable. It was originally written for the C64, but that's easy enough to change. Also, there is some code illustrating it's use.

        list    p=16C84,t=ON,c=132,n=80,st=off
        radix   dec


        include "p16c84.inc"

  cblock  0x20

        f_hi    ;Low byte of frequency variable
        f_lo
        x1
        x2
        interp
        index
        product

        fstep_lo
        fstep_hi
  endc

ODD_QUADRANT    equ     4       ;In f_hi
NEG_QUADRANT    equ     5       ;In f_hi
INDEX_MASK      equ     0x0f    ;In f_hi
INTERP_MASK     equ     0xf0    ;In f_lo


        ORG     0               ;Reset Vector

        GOTO    Main

        ORG     4               ;Interrupt Vector


Main

        BCF     STATUS,RP0      ;Point to BANK 0

    ;Initialize the step in frequency
        MOVLW   0x10
        MOVWF   fstep_lo        ;One unit of fstep corresponds to 360/16384 degrees
        CLRF    fstep_hi        ;i.e. 0.02197 degrees. Set fstep = 455 (1c7) to step
                                ;10 degrees for example

        CLRF    f_hi            ;Start off at 0 degrees
        CLRF    f_lo

xxx
        CALL    sine
        NOP

    ;Now that the sine has been calculated, we could do something with it other than
    ;throwing it away.

    ;Make a frequency step. Note, the frequency step doesn't necessarily have to be
    ;constant. E.g. it could be changed as part of a PLL algorithm.

        MOVF    fstep_lo,W
        ADDWF   f_lo,F
        SKPNC
          INCF  f_hi,F
        MOVF    fstep_hi,W
        ADDWF   f_hi,F

        GOTO    xxx


;--------------------------------------------------------
;sine  -- update 02JAN00
;
;  The purpose of this routine is to take the sine of the
;16-bit frequency variable f_hi:f_lo to produce a signed
;8-bit result that is within +/- 1 count of the true sine
;value.
;  Only the lower 14 bits of the frequency variable are
;actually used. The frequency variable maps into degrees
;by the following transfer function:
;  degrees = (f & 0x3fff) * 360 / 0x4000
;  The sine output is approximately
;sine = int(127*sin( (f & 0x3fff) * 360 / 0x4000) )
; where
; sin() is the true sine function
; int() is the nearest integer function
;
;The technique used to obtain the sine value is a combination
;of table look-up and first order linear interpolation. Sixteen
;equidistant frequency samples of the first quadrant of sin(x) 
;are stored in sine_table.
;
; The frequency variable is broken down as follows:
; xxQQTTTT IIIIPPPP
; where 
;  xx - don't care
;  QQ - Quadrant: 00 = quadrant 1, 01 = quadrant 2, etc.
;  TTTT - Index into the sine_table.
;  IIII - Interpolation between successive entries in the table
;  PPPP - Phase accumulation (not needed in this function, it's
;         only used to increase the dynamic range in frequency
;         steps).
;Once the frequency has been broken down in to these parts, the
;sine function for the first quadrant can be calculated as follows:
;  x1 = sine_table[index]
;  x2 = sine_table[index+1]
;  sine = x1 + ((x2-x1) * interp) / 16
;The first term, x1, is seen to be a first order approximation to
;sine function. The second term improves on that approximation by
;interpolating between x1 and x2. The interpolation variable interp
;is 0 <= interp <= 15 and it is divided by 16. Consequently, the
;interpolation factor ranges between 0 and 15/16.
;
;The sine function in the other three quadrants can be obtained
;from calculations based on the first quadrant by using the following
;trig identities:
; first, let 0 <= f <= 90, i.e. f covers the first quadrant.
;  quadrant 2:  u = 90 + f,   90 < u < 180
;     sin(u-90) = sin(f)
;     x1 = sine_table(16-index), x2 = sine_table(15-index)
;  quadrant 3:  u = 180 + f,  180 < u < 270
;     sin(u) = sin(f+180) = -sin(f)
;     x1 = -sine_table(index), x2 = -sine_table(index+1)
;  quadrant 4:  u = 270 + f,  270 < u < 360
;     sin(u-90) = sin(f+180) = -sin(f)
;     x1 = -sine_table(16-index), x2 = -sine_table(15-index)
;
;Thus, for quadrants 2 and 4, the sine table is indexed in reverse
;order and for quadrants 3 and 4 the values from the sine table 
;are negated. A slight change is made on this indexing and negation
;scheme so that the operation (x2-x1) * interp / 16 only deals with
;positive numbers. This significantly simplifies the multiplication.
;The modification changes the formula for each quadrant as follows:
; quadrant 1:   (no change)
;   x1 = sine_table[index],  x2 = sine_table[index+1]
;   sine = x1 + ((x2-x1) * interp) / 16
; quadrant 2:
;   x1 = sine_table[15-index],  x2 = sine_table[16-index]
;   sine = x2 - ((x2-x1) * interp) / 16
; quadrant 3:
;   x1 = sine_table[index],  x2 = sine_table[index+1]
;   sine = -(x1 + ((x2-x1) * interp) / 16)
; quadrant 4:
;   x1 = sine_table[15-index],  x2 = sine_table[16-index]
;   sine = -(x2 - ((x2-x1) * interp) / 16)
;
;Input
; f_hi:f_lo  -  16-bit frequency variable
;Output
; W = int(127*sin( (f & 0x3fff) * 360 / 0x4000) )
;
;Execution time: 48 Cycles (for all cases)
;

sine

    ;Get the 4-bit index and add 1 to it.
        MOVF    f_hi,W
        ANDLW   INDEX_MASK
        ADDLW   1
        BTFSC   f_hi,ODD_QUADRANT
         SUBLW  17                      ;Odd quadrants, index = 16 - index
                                        ;Actually: (index + 1) = 17 - (index + 1)
                                        ;                      = 16 - index
        MOVWF   index
        CALL    sine_table              ;Get x2=sin(index+1)
        MOVWF   x2

        DECF    index,W
        CALL    sine_table              ;Get x1=sin(index)
        MOVWF   x1

        SUBWF   x2,W                    ;W=x2-x1, This is always positive. 

    ;Initialize the product of (x2-x1)*interp/16 to 1/2. Note 8/16 == 1/2
    ;(This rounds the product to the nearest integer.)
        CLRF    product
        BSF     product,3       ;(note, product and index could be aliased to
                                ; save one byte of ram).

    ;multiply interp and x2 - x1 and divide by 16. This is actually a 4 by 8
    ;bit multiplication. The division by 16 is implemented with a shift right
    ;one position for each of the four multiplication iterations.

	clrc
	btfsc	f_lo,4
         addwf  product,f       ;Then add (x2-x1) to the product
        rrf     product,f       ;Divide the product by two

	clrc
	btfsc	f_lo,5
         addwf  product,f
        rrf     product,f	;Divide the product by four

	clrc
	btfsc	f_lo,6
         addwf  product,f
        rrf     product,f	;Divide the product by eight

	clrc
	btfsc	f_lo,7
         addwf  product,f
        rrf     product,w	;Divide the product by sixteen


        BTFSS   f_hi,ODD_QUADRANT
          ADDWF   x1,W
        BTFSC   f_hi,ODD_QUADRANT
          SUBWF   x2,W
        BTFSC   f_hi,NEG_QUADRANT
          SUBLW 0
        RETURN

sine_table
        ADDWF   PCL,F
        RETLW   0       ;127*sin(0 * 90/16)
        RETLW   12      ;127*sin(1 * 90/16)
        RETLW   25      ;127*sin(2 * 90/16)
        RETLW   37      ;127*sin(3 * 90/16)
        RETLW   49      ;127*sin(4 * 90/16)
        RETLW   60      ;127*sin(5 * 90/16)
        RETLW   71      ;127*sin(6 * 90/16)
        RETLW   81      ;127*sin(7 * 90/16)
        RETLW   90      ;127*sin(8 * 90/16)
        RETLW   98      ;127*sin(9 * 90/16)
        RETLW   106     ;127*sin(10 * 90/16)
        RETLW   112     ;127*sin(11 * 90/16)
        RETLW   117     ;127*sin(12 * 90/16)
        RETLW   122     ;127*sin(13 * 90/16)
        RETLW   125     ;127*sin(14 * 90/16)
        RETLW   126     ;127*sin(15 * 90/16)
        RETLW   127     ;127*sin(16 * 90/16)

        END

Here's a similar routine for the 18Cxxx family (uses table reads and writes): sine18.asm

See Eric Smith's Implementing the Sine function on the PIC ) web page for a similar sine wave routine that does not use interpolation (which makes it really fast!).

Here's some more sine wave theory on this method and others.

And here's more software.

BACK HOME


This page is maintained by Scott Dattalo. You can reach me at: scott@dattalo.com
Last modified on 02JAN00.