# CORDIC help!

#### throbscottle

##### Well-Known Member
Hello lovely people
I'm just getting my head around CORDIC multiplication and division. It sort of makes sense in a fuzzy kind of way.
What I actually want to do is 2 extra levels of complication:
1) 32 bit maths on an 8 bit pic, so I need 4 bytes for everything
2) the actual calculation is word = fout * 2^32 / clkin so I'm thinking is it ought to be possible to do a combined multiply and divide instead of 2 operations
3) I want to do this in assembler
(ok so I can't count!)

It's actually to generate the tuning word for an AD9850 I got ages ago pre-assembled into a module.

It doesn't have to be massively accurate, just better than the tempco of the cheap tcxo that gives clkin (and I've not idea what it's tempco actually is).

I think I can implement the algorithms individually, but don't understand enough to combine them. Any suggestions?
I also want to include a round up or down at the end.

#### jpanhalt

##### Well-Known Member
Have you considered what is sometimes called Kenyan division (a misnomer as the algorithm was not developed in Kenya and a similar algorithm also works for both multiplication)?

I compared integer division using that algorithm with a more classical algorithm, and it was considerably faster particularly if the divisor and dividend were not vastly different, e.g., 3 and a four-byte number. I was looking at duty cycles of 20% to 80%, so the added efficiency was appreciable. The analogous multiplication algorithm exists, but I have not tried it.

There are lots of examples on the Internet, including PICList. If you go that way and are using a PIC, I suggest using an enhanced mid-range device as the new math instructions will save many steps of checking STATUS and carry/borrow. In fact, on a PIC24F with multiple WREG's, it should really fly, but I have not made that jump yet.

John

#### throbscottle

##### Well-Known Member
Thanks John.
I hadn't heard of it, and until a few days ago hadn't heard of any of these methods at all. So it turns out the multiplication method I've found is already the one you suggest, I just didn't know it. It looks like it may be more adaptable to what I want to do if I can get the division to work.
Late now. Will have to try and play with it at the weekend

#### jpanhalt

##### Well-Known Member
Here is the crux of my other post:

For the standard method I chose Peter Hemsley's from Piclist (http://www.piclist.com/techref/microchip/math/div/32by16ph.htm). I have used it many times and simply adapted it to 24-bit X 16-bit for this experiment.

For the Kenyan method, I simply followed the algorithm developed by G. Reichborn-Kjennerud and presented here: http://www.piclist.com/techref/method/math/muldiv.htm The code is written for 24-bit X 24-bit. I will refer to it as the R-K algorithm .

Both methods were simulated on MPLab 8.92 using MPLab SIM at 4 MHz. The pseudo-Kenyan method uses the instruction set from an enhanced mid-range PIC 16F1829, namely the subwfb instruction. My code is still pretty crude. The common subwf(b) method destroys the number being subtracted from (minuend) and my current code uses a shadow register to reconstruct those registers when needed. I don't like doing that, if it can be avoided. It can be done other ways, but the added loop time for doing that may be self-defeating. For example, this destroys the minuend, which ends up with the result:
Code:
     movf      T2L,w     ;
subwf     T1L,f     ;
movf      T2H,w     ;
subwfb    T1H,f     ;
nop
This preserves the minuend, but takes more time:
Code (Microchip Assembler):
Code:
     movf      T2L,w     ;
subwf     T1L,w
movwf     T3L     ;
movf      T2H,w     ;
subwfb    T1H,w   ;
movwf     T3H
nop
Cutting to the chase, here are my timing results:

The Hemsley method uses rotations and has a relatively fixed time. The alternative method depends on the number of subtractions needed. It can be quite short when a number is divided by itself or long when the divisor is 1. My code does not check (yet) for a zero divisor.

Here's a complete, hopefully working version (trial subtraction method) that I used for simulation:
Code:
;Author: ©JPANHALT
;Date: May 22,2016

list        p=16f1829      ; list directive to define processor
#include    <p16f1829.inc> ; processor specific variable definitions
errorlevel -302,-305

__CONFIG _CONFIG1, _FOSC_INTOSC & _WDTE_OFF & _PWRTE_OFF & _MCLRE_ON & _CP_OFF & _CPD_OFF & _BOREN_OFF & _CLKOUTEN_ON & _IESO_OFF & _FCMEN_OFF
__CONFIG _CONFIG2, _WRT_OFF & _PLLEN_OFF & _STVREN_OFF & _BORV_19 & _LVP_OFF

cblock    0x20
T1U
T1H
T1L
T2U
T2H
T2L
count
resultL
resultH
endc

org 0x0000

goto      Start
Start
;T1(24bit)/T2(16bit)
movlw     high(881)
movwf     T2H
movlw     low(881)
movwf     T2L
call      Refresh
clrf      count
clrf      T2U
clrf      resultL
clrf      resultH
nop
Kenyan
incf      count,f
movf      T2L,w
subwf     T1L,f
movf      T2H,w
subwfb    T1H,f
movf      T2U,w
subwfb    T1U,f
btfss     STATUS,0
goto      Calculate
lslf      T2L
rlf       T2H
rlf       T2U
goto      Kenyan

Calculate
call      Refresh

Cycle
;T1 registers are destroyed during subtraction.  When the subtrahend is too
;big, which results in a zero being rotated into the result, those registers
;must be recovered.  Back-up registers can be used; however, several steps are
;required to restore them.  Doing a trial subtraction and preserving the
;minuend registers was slightly faster.
movf      T2L,w
subwf     T1L,w
movf      T2H,w
subwfb    T1H,w
movf      T2U,w
subwfb    T1U,w
btfss     STATUS,0       ;subtrahend too big divide by 2 and repeat
goto      Result
movf      T2L,w          ;subtractand OK, repeat subtraction w/saves in f
subwf     T1L,f
movf      T2H,w
subwfb    T1H,f
movf      T2U,w
subwfb    T1U,f          ;carry will be set
Result
rlf       resultL,f      ;carry bit is set or clear correctly from above
rlf       resultH,f
decfsz    count,f
goto      DivideT2
goto      Finish
DivideT2
lsrf      T2U
rrf       T2H
rrf       T2L
goto      Cycle

Finish
nop                      ;answer in resultH and resultL

;Subroutine

Refresh
movlw     upper(55536)   ;(673312)
movwf     T1U
movlw     high(55536)    ;(673312)
movwf     T1H
movlw     low(55536)     ;(673312)
movwf     T1L
return

end
None of the code is really polished. I did it just out of curiosity to get something that worked. "Individual results may vary. "

Regards, John

#### BobW

##### Active Member
I have PIC assembly code for generating the control word for an AD9850/51 to 1 Hz resolution. It uses a combination of math operations and lookup tables. I can either post it or email it to you.

#### jpanhalt

##### Well-Known Member
Please post it to be consistent with ETO guidelines. I keep a file of Assembly routines, since that is all I know.

Admittedly, I enjoy Assembly, but feel like John (the other one) crying into the wilderness.

John

#### BobW

##### Active Member
My only concern was that the source code includes a lot of other stuff making it rather large. It's a complete sweep generator application. But I will post it here (I may clean out the extraneous bits first). Won't be able to get to it until later tonight though.

#### jpanhalt

##### Well-Known Member
Hi Bob,

Either way is OK with me. I am used to stripping what I need for my subroutines file.

I enjoy just diddling with Assembly from time to time with no immediate purpose in mind.

John

#### BobW

##### Active Member
Here is my assembly code:
Code:
; Code for generating Frequency Control Word for AD9850 and AD9851
;
; R. Weaver, 2017-08-01
;
; Assembler directives to allow for different reference frequencies
;    #define RefOsc180MHz ;Comment out this line for 125MHz reference osc (AD9850)
#define LeftShifts 1 ; Set this value according to the following table:
; Ref Osc Freq       LeftShifts
; ------------       ----------
;   19-37 MHz           3
;   38-75 MHz           2
;  76-150 MHz           1
; 151-200 MHz           0
;
;  ****
;  Note that for a reference frequency other than 125 or 180 MHz
;  Changes must be made to the lookup table code.
;  ****

;Macro definitions
bank0 macro
bcf STATUS,RP0
endm
bank1 macro
bsf STATUS,RP0
endm
;More intuitive skip definitions
skpeq macro
btfss status,z
endm
skpneq macro
btfsc status,z
endm
skpneg macro
btfsc status,c
endm
skppos macro
btfss status,c
endm
SaveAccHi macro rName ;Store 4 high bytes of acc to spec'd register set
movlw rName
call SaveAccHiInd
endm
SaveAccLo macro rName ;Store 4 low bytes of acc to spec'd register set
movlw rName
call SaveAccLoInd
endm
LoadAccHi macro rName ;Load 4 high bytes of acc from spec'd register set
movlw rName
endm
LoadAccLo macro rName ;Load 4 low bytes of acc from spec'd register set
movlw rName
endm
LoadRLo macro rName ;Load 4 low bytes of aux reg from spec'd register set
movlw rName
endm
LoadRHi macro rName ;Load 4 high bytes of aux reg from spec'd register set
movlw rName
endm
ClearNregisters macro nReg, rName ;Clear nReg registers starting at rName
movlw nReg
movwf xfrCount
movlw rName
call ClearNregistersInd
endm
ShiftUpNregisters macro nReg, rName ;Shift nReg registers starting at rName
movlw nReg
movwf xfrCount
movlw rName
call ShiftUpNregistersInd
endm

;
;variable memory area
cblock 0x20
;   DDS variables:
freqDigits:8 ;8 digit (decimal) frequency in Hz, F(0) is LSB, F(7) is MSB
;   Math Library variables
acc0,acc1,acc2,acc3,acc4,acc5,acc6,acc7 ;accumulator for 32/64 bit math
R0,R1,R2,R3,R4,R5,R6,R7 ; 32/64 bit math aux. registers
;   misc temporary registers
count,temp
xfrCount ;used to count registers moved by the load, save & clear routines
endc

;Lookup tables for DDS frequency ctrl word
#ifdef RefOsc180MHz ;Lookup tables for 180MHz ref frequency
Mpy_B0 ;Least significant byte
andlw 0x0F
retlw 0x00
retlw 0xE1
retlw 0xC1
retlw 0xA2
retlw 0x82
retlw 0x63
retlw 0x44
retlw 0x24
retlw 0x05
retlw 0xE5
Mpy_B1
andlw 0x0F
retlw 0x00
retlw 0xDE
retlw 0xBD
retlw 0x9C
retlw 0x7B
retlw 0x5A
retlw 0x39
retlw 0x18
retlw 0xF7
retlw 0xD5
Mpy_B2
andlw 0x0F
retlw 0x00
retlw 0x65
retlw 0xCB
retlw 0x31
retlw 0x97
retlw 0xFD
retlw 0x63
retlw 0xC9
retlw 0x2E
retlw 0x94
Mpy_B3
andlw 0x0F
retlw 0x00
retlw 0xDC
retlw 0xB8
retlw 0x95
retlw 0x71
retlw 0x4D
retlw 0x2A
retlw 0x06
retlw 0xE3
retlw 0xBF
Mpy_B4 ;Most significant byte
andlw 0x0F
retlw 0x00
retlw 0x17
retlw 0x2F
retlw 0x47
retlw 0x5F
retlw 0x77
retlw 0x8F
retlw 0xA7
retlw 0xBE
retlw 0xD6
#else ;Lookup tables for 125 MHz reference oscillator
Mpy_B0 ;Least significant byte
andlw 0x0F
retlw 0x00
retlw 0x27
retlw 0x4E
retlw 0x75
retlw 0x9B
retlw 0xC2
retlw 0xE9
retlw 0x10
retlw 0x37
retlw 0x5E
Mpy_B1
andlw 0x0F
retlw 0x00
retlw 0xE8
retlw 0xD0
retlw 0xB8
retlw 0xA0
retlw 0x88
retlw 0x70
retlw 0x59
retlw 0x41
retlw 0x29
Mpy_B2
andlw 0x0F
retlw 0x00
retlw 0x0B
retlw 0x17
retlw 0x23
retlw 0x2F
retlw 0x3B
retlw 0x47
retlw 0x53
retlw 0x5F
retlw 0x6B
Mpy_B3
andlw 0x0F
retlw 0x00
retlw 0x2E
retlw 0x5C
retlw 0x8A
retlw 0xB8
retlw 0xE6
retlw 0x14
retlw 0x42
retlw 0x70
retlw 0x9E
Mpy_B4 ;Most significant byte
andlw 0x0F
retlw 0x00
retlw 0x11
retlw 0x22
retlw 0x33
retlw 0x44
retlw 0x55
retlw 0x67
retlw 0x78
retlw 0x89
retlw 0x9A
#endif
;
;
;
MakeCtlWd   ;Generate DDS frequency control word
call clrAcc64 ; Init the acc to zero
MakeWdEntry
movwf FSR    ; and store in indirect ptr reg
movlw 0x08  ; load count of digits
movwf count
;    call clrAcc64 ; Init the acc to zero
goto LoopMCWenter ;skip the x10 for the 1st digit
LoopMCW
call AccTimes10
LoopMCWenter
movf indf,w
decf fsr,f
decfsz count,f
goto LoopMCW
;
;All digits have been processed.
;If Ref freq is 150MHz or lower then do left shift here
;
if LeftShifts != 0
movlw LeftShifts
call AccLShiftN
endif
btfsc acc3,7    ; Check fractional part to see if round up is req'd
incfsz acc4,f   ;Yes, so increment integer part, and propagate carries.
goto MCWshiftWord
incfsz acc5,f
goto MCWshiftWord
incfsz acc6,f
goto MCWshiftWord
incf acc7,f
MCWshiftWord    ;move the result from acc4..acc7 to acc0..acc3
return  ;Frequency/sweep Control Word is in acc0..acc3
AddDigit    ; Add next frequency digit to acc and mpy by 10
movwf temp
call mpy_B0
movwf R0
movf temp,w
call mpy_B1
movwf R1
movf temp,w
call mpy_B2
movwf R2
movf temp,w
call mpy_B3
movwf R3
movf temp,w
call mpy_B4
movwf R4
clrf R5
clrf R6
clrf R7
return
AccTimes10    ;Multiply acc x 10
;copy acc to aux reg.
call AccToR ;move acc contents to aux reg
movlw 0x02
call accLShiftN ;Shift left twice to get x4
call Add64    ; add the saved value back to the acc to get x5
movlw 0x01
call accLShiftN ;shift left once more to get x10
return
AccLShiftN ;Shift acc left N times where N is in Wreg
movwf count1
LoopALSN
bcf status,c
rlf acc0,f
rlf acc1,f
rlf acc2,f
rlf acc3,f
rlf acc4,f
rlf acc5,f
rlf acc6,f
rlf acc7,f
decfsz count1,f
goto LoopALSN
return
;
;Start of general purpose multi-byte math routines
;
clrAcc64 ;clear 64 bit accumulator registers acc and r
ClearNregisters d'16',acc0
return

SaveAccHiInd ; Save high 4 bytes of acc to reg indirect
movwf fsr ;start address is in w
movf acc4,w
movwf indf
movf acc5,w
incf fsr,f
movwf indf
movf acc6,w
incf fsr,f
movwf indf
movf acc7,w
incf fsr,f
movwf indf
return
SaveAccLoInd ; Save low 4 bytes of acc to reg indirect
movwf fsr ;start address is in w
movf acc0,w
movwf indf
movf acc1,w
incf fsr,f
movwf indf
movf acc2,w
incf fsr,f
movwf indf
movf acc3,w
incf fsr,f
movwf indf
return
movwf fsr ;start address is in w
movf indf,w
movwf acc4
incf fsr,f
movf indf,w
movwf acc5
incf fsr,f
movf indf,w
movwf acc6
incf fsr,f
movf indf,w
movwf acc7
return
movwf fsr ;start address is in w
movf indf,w
movwf acc0
incf fsr,f
movf indf,w
movwf acc1
incf fsr,f
movf indf,w
movwf acc2
incf fsr,f
movf indf,w
movwf acc3
return
LoadRxLoInd ; Load low 4 bytes of aux reg r0 from reg indirect
movwf fsr ;start address is in w
movf indf,w
movwf r0
incf fsr,f
movf indf,w
movwf r1
incf fsr,f
movf indf,w
movwf r2
incf fsr,f
movf indf,w
movwf r3
return
LoadRxHiInd ; Load low 4 bytes of aux reg r0 from reg indirect
movwf fsr ;start address is in w
movf indf,w
movwf r4
incf fsr,f
movf indf,w
movwf r5
incf fsr,f
movf indf,w
movwf r6
incf fsr,f
movf indf,w
movwf r7
return
ClearNregistersInd
;Clear specified number of contiguous registers where number is in xfrCount
;and start register address is in w
movwf fsr
LoopCNRI
clrf indf
incf fsr
decfsz xfrCount,f
goto LoopCNRI
return
ShiftUpNregistersInd
;Shift bytes up across N registers
;highest register is in w
movwf fsr
LoopSUNRI
decf fsr,f
movf indf,w
incf fsr,f
movwf indf
decf fsr,f
decfsz xfrCount,f
goto LoopSUNRI
return

AccToR ;move acc contents to aux reg.
;No indirects here, so that calling code can use FSR without it
;getting corrupted here.
movf acc0,w
movwf R0
movf acc1,w
movwf R1
movf acc2,w
movwf R2
movf acc3,w
movwf R3
movf acc4,w
movwf R4
movf acc5,w
movwf R5
movf acc6,w
movwf R6
movf acc7,w
movwf R7
return

movf R0,w
btfsc status,c
call carry1
movf R1,w
btfsc status,c
call carry2
movf R2,w
btfsc status,c
incf acc3,f
movf R3,w
return
neg32    ;Negate the accumulator
comf acc0,f
comf acc1,f
comf acc2,f
comf acc3,f
incfsz acc0,f ;increment and check for carry
return ;no carry, so return
carry1
incfsz acc1,f
return
carry2
incfsz acc2,f
return
incf acc3,f
return
movf R0,w
btfsc status,c
call carry64_1
movf R1,w
btfsc status,c
call carry64_2
movf R2,w
btfsc status,c
call carry64_3
movf R3,w
btfsc status,c
call carry64_4
movf R4,w
btfsc status,c
call carry64_5
movf R5,w
btfsc status,c
call carry64_6
movf R6,w
btfsc status,c
incf acc7,f
movf R7,w
return
carry64_1   ;Chained skips to propagate carries
incfsz acc1,f
return
carry64_2
incfsz acc2,f
return
carry64_3
incfsz acc3,f
return
carry64_4
incfsz acc4,f
return
carry64_5
incfsz acc5,f
return
carry64_6
incfsz acc6,f
return
incf acc7,f
return
NegDecimal ; Negate 8-digit unpacked decimal number in acc0..7
movlw 0x08
movwf count
movlw acc0
movwf fsr
LoopND
movf indf,w
sublw 0x09
movwf indf
incf fsr,f
decfsz count
goto LoopND
;Acc now contains nines complement. Now increment to get tens complement
ClearNregisters 8,r0
incf r0,f
AddDecimal ; Eight digit unpacked BCD add; operands in acc0..7 and r0..7
;fsr indf
movlw 0xF6
movwf r0
movwf r1
movwf r2
movwf r3
movwf r4
movwf r5
movwf r6
movwf r7
;new code
movlw acc0
movwf fsr
movlw 0x08
movwf count
movf indf,w
sublw 0x09
rrf temp,f ; save carry bit
incf fsr,f
decfsz count,f
comf temp,f ;
btfsc status,z
movlw r0
movwf fsr
movlw 0x08
movwf count
rrf temp,f
btfss status,c
clrf indf
incf fsr
decfsz count,f
goto LoopAD2
The desired frequency in decimal Hz should be entered in the register array freqdigits. Then call subroutine MakeCtlWd to generate the frequency control word. The result is in registers acc0..acc3.

The calculation is done digit by digit. The 5 byte constant is looked up for the current digit and added to the accumulated result. This is then multiplied by 10. This is repeated until all 8 digits have been processed. The lookup table has 5 byte precision to prevent round off error.

I trimmed this out of a much larger program, so there are probably one or two small subroutines here that aren't used.

Last edited:

#### jpanhalt

##### Well-Known Member
Thank you, Bob,

I will definitely be playing with it. My method (http://www.electro-tech-online.com/articles/demo-assembly-code-ad9850-dds-signal-generator.755/ ) used brute force and the math part took a little more than 1.3 ms.

Right now, though, I am kept busy dealing with rodent damage to my vehicles* and downed ash trees. Bad weather helps, but we have had a very dry Summer so far. Beautiful weather, but the ground is like cinder block.

Regards, John

*Chewed through ABS in multiple places on my F-150, vacuum control of the front locking axle hubs, and the electronic shift on fly (ESOF) for FWD. About $260 in non-Ford dealer parts (twice that at dealer prices) plus the time to diagnose and fix. I think the rascals are trained at Ford to eat where a repair is not practical, like right where the wire enters a casting on the FWD shift motor. #### BobW ##### Active Member A bit of additional info: 1. The way I arranged this, nothing more complicated than multiply by 10 is required. 2. I realized after writing the code, that it could be done without the lookup tables, by storing just a single constant based on the reference oscillator frequency, but I have all of this (and much more code) programmed into a 16F630 and started to run low on register memory. Eliminating the lookup tables would have required more registers which I didn't have enough of. On the other hand, using a single constant would allow for easier fine tuning in case the reference oscillator is slightly off frequency. 3. The code has provision for both 180 MHz ref oscillator (AD9851) and 125 MHz ref oscillator (AD9850). There is a set of lookup tables for each one, and the correct one is chosen by setting an assembler directive (see code comments). If a different ref oscillator frequency is needed then a new set of lookup table data needs to be generated. I created a spreadsheet that does this. See attached. 4. The calculation principle may be a bit obscure. So, if there is any confusion about how it works, let me know and I'll provide more detail. Sorry to hear about your rodent problem. Ford should learn not to use flavor ingredients in their parts that appeal to rodents. It's been hot and dry here too. #### Attachments • 18.5 KB Views: 28 #### jpanhalt ##### Well-Known Member Most Helpful Member I have heard that peanut oil is used in its PVC insulation. Not sure exactly how that happens (proteins?), but that is what they attack. Peanuts are like sex to rodents. By comparison, have yet to have a rodent attack regular household wiring. What I have learned to do is a good scan (Autel, not HF scanner) of my system, then investigate for broken wires. Ford dealer was clueless on an A/C problem with my Fusion, which usually means replacing parts until the fault is located. A several hundred$ repair of the A/C on my Fusion ended up costing me a piece of wire, some heat-shrink tubing, and the Autel scanner that has more than paid for itself.

Today, the hardest part of working under a car for me is getting on the crawler at my age. My new dog is great company, except he does get a little close at times. His breath is enough to keep you awake.

John

#### throbscottle

##### Well-Known Member
Wow. When I eventually have enough time off work without a load of other things to do and I'm not actually knackered, (which is a big one, I seem to be permanently knackered) I'll study these!
Thanks guys, big help
Bob, is the lookup table working the same way it does when you do trig using the cordic method? Going to take me a while to understand what's gonig on there!
Sounds like you need a pit, John!
Feed the dog raw, he'll thank you for it and his breath won't smell (well, hardly anyway). You just have to make sure it's balanced.

#### BobW

##### Active Member
Bob, is the lookup table working the same way it does when you do trig using the cordic method.
No. Actually in the simplest version of the algorithm, the lookup table is simply the values of 2^32/FREF multiplied by the digits 0-9. If we refer to the value 2^32/FREF as constant k, and the tuning control word as W, then:

W = k * FOUT

Eight decimal digits allow us to represent frequencies from 0 to 99,999,999 Hz with 1 Hz resolution. We can represent FOUT as the sum of its eight decimal digits (D0..D7) multiplied by the appropriate powers of 10. Hence:

FOUT = D7*10^7 + D6*10^6 + D5*10^5 + D4*10^4 + D3*10^3 + D2*10^2 + D1*10^1 + D0*10^0

D0 is the least significant digit and D7 is the most significant.
And then the tuning word will be:

W = k*(D7*10^7 + D6*10^6 + D5*10^5 + D4*10^4 + D3*10^3 + D2*10^2 + D1*10^1 + D0*10^0)

We can rearrange this formula as follows:

W = k*D7*10^7 + k*D6*10^6 + k*D5*10^5 + k*D4*10^4 + k*D3*10^3 + k*D2*10^2 + k*D1*10^1 + k*D0*10^0

Since D0..D7 are decimal digits, they can only have the values 0..9, and so we can create a lookup table for the ten possible values of k*DN.
Now that we've done that, the only remaining operations are multiplication by powers of 10. Again rearranging the formula, we get:

W = ((((((k*D7*10 + k*D6)*10 + k*D5)*10 + k*D4)*10 + k*D3)*10 + k*D2)*10 + k*D1)*10 + k*D0

With this rearrangement, it is only necessary to use the lookup table to get the values k*DN and then to multiply by 10, which is very easy.

In order to minimize roundoff error, I used 5 byte precision rather than 4 bytes. That's the reason for 5 lookup tables. There's one for each byte in the table output value. That's basically how the algorithm works for a 180 MHz ref frequency as used by an AD9851. Since the AD9850 ref frequency is a maximum of 125 MHz, I made a slight change to the algorithm to maximize precision. Instead of making k = 2^32/FREF, I use the smallest power of two that won't overflow 5 bytes. For FREF = 125MHz, the value is 2^31. The algorithm proceeds exactly the same way with the only difference being that the final result needs to be shifted left one bit position.

I find that it's often possible to rearrange a calculation to make it much easier to do in assembly language, and that is the case here.

Last edited:

#### throbscottle

##### Well-Known Member
Ohhhhhhhhh I see now. Elegant.
But isn't all that multiplying by 10 messy?
I'd been wondering about the rounding error. I couldn't get past outright truncation.

#### BobW

##### Active Member
The multiplication by 10 is in the same loop that sums the lookup table value. Seven lines of code (not counting comments).

Code:
AccTimes10    ;Multiply acc x 10
;copy acc to aux reg.
call AccToR ;move acc contents to aux reg
movlw 0x02
call accLShiftN ;Shift left twice to get x4
call Add64    ; add the saved value back to the acc to get x5
movlw 0x01
call accLShiftN ;shift left once more to get x10
return

#### throbscottle

##### Well-Known Member
Oh of course, x8 and add 2 more, only done a better way.
See if I can find time tomorrow for a bit of contrast and compare...
Thanks for this, saved me many hours (and wasted ticks/registers!) I think leaves some space in my brain to think about other functions...

#### throbscottle

##### Well-Known Member
Ok so it's been nearly a year...
I found your website, Bob and followed the explanation there since you'd made it easier to understand, so I'm writing my own implementation now.
I played some more with the "fast multiplication" I found at convict.lu but there's no gain in using it. Although it did get me to making a working "reverse double-dabble", which I'm now not using. So I'm just adding k onto the accumulator as many times as a bcd digit gives a count down for, then the x10, next digit, etc etc.
I'm using a vfd for my display. The controller chip also has outputs for led's and inputs for keys - more blinkenlampen!
I'd been wanting to include a sweep function too but kind of "knew what I wanted without really knowing what I wanted" with that - so you've helped there too!
I had another look at the spreadsheet you used to calculate the original LUT values. Still don't understand the maths there - particularly how you're using MOD to generate the values and LOG to get the bit width.

#### BobW

##### Active Member
You got there just in the nick of time. I think I only had that page up for a couple of weeks before I rewrote the PIC program so that it no longer needed lookup tables. Since the tables were gone, I removed the link to the spreadsheet file. Maybe I should put it back again.

As for the MOD function, it simply returns the remainder of an integer division. To convert a number into individual bytes, you divide by 256. The remainder is the low order byte, the integer quotient is the high byte. If the number takes more than two bytes, the process is repeated with the integer quotient until it's less than 256. So for a number like 1234567, which would require 3 bytes, you divide by 256 giving a quotient of 4822 and a remainder of 135.
MOD(1234567,256)=135 <-- this is byte 0 (least significant)
INT(1234567/256)=4822
Repeat:
MOD(4822,256)=214 <-- this is byte 1
INT(4822/256)=18 <-- this quotient is now less than 256, so this is byte 2 (most significant), and we're done.

As for the LOG function, it will inherently give you the number of digits (minus 1) required to represent the number in whatever base you're working in. In base 10 (decimal), use the base 10 log to tell how many decimal digits a number has. For example, the number 1234567 has 7 digits,
LOG(1234567,10)=6.09151
Take the integer part, 6, and add 1, giving 7, which is the number of digits in the number.
For calculating the number of bits required to represent the number in binary, change the base of the log function to 2:
LOG(1234567,2)=20.23557
Take the integer part, 20, and add 1, giving 21, which is the number of bits required to represent 1234567.

#### throbscottle

##### Well-Known Member
Sticking with the BCD maths, I started to implement the sweep function going on your principle of subtracting half the sweep range from the centre frequency as a start value. Taken me all week to figure out how to do subtraction with more than one byte!
Code:
; take 1234 from 4321

TEST
banksel val1
pagesel TEST
movlw 0x43
movwf val1+1
movlw 0x21
movwf val1
movlw 0x12
movwf val2+1
movlw 0x34
movwf val2

movlw 0x99 ; convert BCD to 9's complement
movwf temp
movfw val1
subwf temp,w
movwf val1
movfw val1+1
subwf temp,w
movwf val1+1
movfw val1
subwf temp,w
movwf val1 ; val1 ends up with result
clrw
btfss STATUS,DC
movlw 0x06 ; low digit adjust
btfsc STATUS,C
goto NEXT1
incf val1+1,f ; borrow from next digit, is add because took from 9
iorlw 0x60 ; high digit adjust
NEXT1
movfw val1+1
subwf temp,w
movwf val1+1
clrw
btfss STATUS,DC
movlw 0x06
btfsc STATUS,C
goto NEXT2
incf val1+2,f
iorlw 0x60
NEXT2
subwf val1+1,f
I'm sure it could be more compact though.