This blog post will discuss how to find the best rational approximation with bounded numerator or denominator to a a preferred frequency ratio. This will be a bit mathematical but actually nothing more than standard addition, subtraction, multiplication and division will be used. Some example code is given at the end.
Fractional baudrate generation
Some microprocessors uses fractional multipliers and divisors to generate accurate baud rates from system clocks that are not an integer multiple of the bitslice clock frequency. The bitslice frequency is the baudrate times an oversampling count.
fbitslice=oversample⋅baudrate=ratio⋅fsysclock
An example is the USIC in the Infineon XMC1100 and XMC4500 family microprocessors. Here the clock for the bitslices is derived from the system clock, slightly simplified, as
fbitslice=(DCQT+1)baudrate=1(PDIV+1)⋅STEP1024⋅fsysclock
DCQT gives the oversampling count, typically set to 15. PDIV and STEP are register values controlling the frequency division. On the XMC1100 PDIV and STEP are 10 bit wide so they have values from 0 to 1023. Rearranging the above we see that we must find a good approximation
STEP(PDIV+1)≈(DCQT+1)baudrate⋅1024fsysclock
with values of PDIV and STEP between 0 and 1023.
Example: 38400 baud from a 32MHz clock and 16 times oversampling
We enter the values in (2) to find
STEP(PDIV+1)=16⋅38400⋅102432000000=62914560032000000=19.6608
Working by hand we would now factor out common factors from numerator and denominator, but for computer implementations that is unnecessary extra work. As we shall see the algorithm works as well without that step.
Now we calculate the signed integer division with remainder using the rounded value 20 as quotient.
16⋅38400⋅102432000000=20−1085440032000000=20−13200000010854400
Rewrite the denominator in the last fraction by carrying out the division with remainder
3200000010854400=3−56320010854400≈2.948≈3
We enter this value 3 in the previous formula (3) to get the approximation
16⋅38400⋅102432000000≈20−13=593=19.6667
This is good, after one iteration the error is less than 0.008 or about 300ppm.
We can improve this by using a better using a more precise value
3200000010854400=3−56320010854400=3−110854400563200=3−119+153600563200≈3−119=5619
Inserted into (3) this gives
16⋅38400⋅102432000000≈20−15619=20−1956=110156=19.6607
This approximation is better, error is about 3ppm, but the STEP value 1101 is to large for the 10 bit range, so we keep the previous approximation of 59/3.
Best rational approximation with continued fractions
In order to find a general method to calculate good rational approximations as in the previous example we use the theory of
partial quotients, or convergents, of continued fraction expansions and
best rational approximation theory. [TO DO add links ] This is a well developed general theory for finding the best rational apprimations to a number with a given size of the denominator. There are recursive algoritms for calculating succusively better approximations.
So we are given a rational number p/q and we want to find succesive approximations
an/bn to p/q with smaller values for
an and
bn than p and q. For the baud rate generations we want the best approximation where
an and
bn fits in the respective fractional multiplier and divisor registers.
The formula
pq=c0+1c1+1c2+⋱1cn+1rn+1≈c0+1c1+1c2+⋱1cn=an+1bn+1
For the first few values of n we get:
a−1=0,b−1=1
a0=1,b0=0
a1b1=c0=0+1⋅c01+0⋅c0=a−1+a0⋅c0b−1+b0⋅c0
a1=c0,b1=1
a2b2=c0+1c1=c0⋅c1+1c1=a0+a1⋅c1b0+b1⋅c1
a3b3=a0+a1(c1+1c2)b0+b1(c1+1c2)=(a0+a1c1)c2+a1(b0+b1c1)c2+b1=a1+a2c2b1+b2c2
We now can see the recursion formula develop
a4b4=a1+a2(c2+1c3)b1+b2(c2+1c3)=(a1+a2c2)c3+a2(b1+b2c2)c3+b2=a2+a3c3b2+b3c3
an+1bn+1=an−2+an−1(cn−1+1cn)bn−2+bn−1(cn−1+1cn)=(an−2+an−1cn−1)cn+an−1(bn−2+bn−1cn−1)cn+bn−1=an−1+ancnbn−1+bncn
The algorithm
We start with a rational number r = p/q. We will generate a sequence of five values
an,bn,cn,pn and
qn. The quotients
an/bn will be our succesive approximations to r. The values
rn are defined as
pn/qn but they need not be explicitly calculated, the formula for
rn is used as the template for how
pn and
qn are updated.
Startup:
r0=r
p0=p,q0=q
a−1=0,b−1=1
a0=1,b0=0
Loop until q_n+1 is 0 or an or bn are to large
cn=round(rn)=round(pn/qn)
rn+1=1/(rn−cn) This is done calculated terms of
pn and
qn
pn+1=qn
qn+1=pn−cn⋅qn
an+1=an−1+an⋅cn
bn+1=bn−1+bn⋅cn
The rounding operations can be done downwards, keeping all values positive, or towards the nearest integer improving precision but also introducing negative numers and extra complexity. Note that the
cn value is calculated before updating
pn and
qn and is not remembered and used in the next iteration but recalculated.
A skeleton C implementation
Following code lacks error handling and only checks limit for numerator, but it illustrates the algorithm. It does lack a few comments, but follows the algortithm structure closely.
void cfractr(int32_t p, int32_t q,int32_t alim, uint32_t * ares, uint32_t * bres) {
int ap = 0;
int a1 = 1;
int bp = 1;
int b1 = 0;
int cn = a1, anext, bnext, pnext;
while (1) {
/* Signed rounded rational division */
if ((q>0)&&(p>0)||(q<0)&&(p<0))
cn = (p+q/2)/q;
else
cn = (p-q/2)/q;
/* Next value for partial quotients and remainder */
anext = ap + cn*a1;
bnext = bp + cn*b1;
pnext = p-cn*q;
/* Exact value, remainder is 0, break */
if (pnext == 0) {
a1 = anext;
b1 = bnext;
break;
}
/* Numerator too large, break */
if ((anext<-alim)||(anext>alim)) {
break;
}
/* Shift one step before next iteration */
ap = a1;
bp = b1;
a1 = anext;
b1 = bnext;
p = q;
q = pnext;
}
if (a1<0){a1=-a1;b1=-b1;}
*ares = a1;
*bres = b1;
}
