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.

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