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.
| CODE |
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 |