How to calculate the parameter values for a fractional frequency divider with post-scaler

Question

I need to calculate the parameter for a fractional frequency divider on a micro-controller:

fout = fin / (i + n/d) / p

for a given fin and fout. All parameters are integers with the following allowed ranges:

i: integer part 2..256
d: divisor 1..512
n: dividend 0..511
p: post-scaler 1..64

The calculation has to be done on the controller itself and should be reasonably fast (around 1 ms)

A simple approach using brute-force searching through all 32768 possible d/p combinations (with different rounding for n and i to check for above/below) will be too slow. Using only the lowest possible p to keep i below 256 would miss better fitting parameters.

Is there a suitable algorithm to calculate the parameters? The closest question i could find is this but this does not really give a solutions and does not include the additional complexity of the post-scaler.

If there is no clever algorithm, what optimizations for the brute-force approach are possible? One easy one is to only check for d 256..512 because for d 1..255 there is a combination with identical value in the upper half. Is it somehow possible to eliminate some d/p combination that are also redundant?

I need different parameter sets depending on the situation:

1. the parameters that produce the closest possible fout above a given frequency.
2. the parameters that produce the closest possible fout below a given frequency.
3. the parameters that produce the closest possible or exact fout for a given frequency.

1 and 2 will be used to increase/decrease the fout by the smallest possible steps.

Answer 1

NOTE: I'm leaving the original answer, even though I figured out how to make it provably correct. I'll insert that at the bottom.

---

Here is an approach that will find a good answer very quickly. I don't have a proof, but I think it will always be optimal. If not always, then often.

If we rearrange:

fout = fin / (i + n/d) / p

we get

fin / fout = (i + n/d) * p = (d*i*p + n*p) / d

So we're looking for a small fraction that comes close to a known number, that you can create from your parameters.

Now let's talk about the Stern-Brocot Tree. (Incidentally Brocot was a clock maker trying to come up with good gearing systems, which isn't that different from your problem.)

Suppose that x is a positive number. We start the tree with:

0 / 1 < x  < 1 / 0 (aka infinity)

(It is usual to start with 1 instead of infinity on the right, but this works and lets you handle numbers bigger than 1.)

At each step we have:

a / b < x < c / d

Our next guess is (a + c) / (b + d). We compare that, replace one side, and continue.

The numbers we generate turn out to include all of the best possible rational approximations. Meaning that if we have a / b in this sequence, then there is no fraction with denominator smaller than b that is as good an approximation. And any fraction with that property WILL occur in this sequence.

Let's demonstrate it on finding approximations to e = 2.718281828459045 with a denominator at most 512. This is doable by hand.

   0 /   1 < e <    1 /   0, guess    1 /   1 = 1                 < e
   1 /   1 < e <    1 /   0, guess    2 /   1 = 2                 < e
   2 /   1 < e <    1 /   0, guess    3 /   1 = 3                 > e
   2 /   1 < e <    3 /   1, guess    5 /   2 = 2.5               < e
   5 /   2 < e <    3 /   1, guess    8 /   3 = 2.666666666666666 < e
   8 /   3 < e <    3 /   1, guess   11 /   4 = 2.75              > e
   8 /   3 < e <   11 /   4, guess   19 /   7 = 2.714285714285714 < e
  19 /   7 < e <   11 /   4, guess   30 /  11 = 2.727272727272727 > e
  19 /   7 < e <   30 /  11, guess   49 /  18 = 2.722222222222222 > e
  19 /   7 < e <   49 /  18, guess   68 /  25 = 2.72              > e
  19 /   7 < e <   68 /  25, guess   87 /  32 = 2.7185            > e
  19 /   7 < e <   87 /  32, guess  106 /  39 = 2.717948717948718 < e
 106 /  39 < e <   87 /  32, guess  193 /  71 = 2.718309859154929 > e
 106 /  39 < e <  193 /  71, guess  299 / 110 = 2.718181818181818 < e
 299 / 110 < e <  193 /  71, guess  492 / 181 = 2.718232044198895 < e
 492 / 181 < e <  193 /  71, guess  685 / 252 = 2.718253968253968 < e
 685 / 252 < e <  193 /  71, guess  878 / 323 = 2.718266253869969 < e
 878 / 323 < e <  193 /  71, guess 1071 / 394 = 2.718274111675127 < e
1071 / 394 < e <  193 /  71, guess 1264 / 465 = 2.718279569892473 < e
1264 / 465 < e <  193 /  71, guess 1457 / 536 ...and we're out of range

Optimization note, the first few steps are literally just finding what integers you're between. You can skip them by starting with int(fin/fout) / 1 <= x < (1 + int(fin/fout)) / 1. From there, at most 511 guesses are needed, and usually you need considerably less than that.

And now our candidates for a good fraction below the number we want are:

1264/465, 1071/394, 878/323, ...

And our candidates for a good fraction above the number we want are:

193/71, 87/32, 68/25, ...

Now how do we try to meet the parameters? Well, we have the following rules. If the fraction we want is a/b then:

1. d = b
2. p is at most min(64, a//d//2) where // is integer division. (If p < 1 then no solution.)
3. p divides a and is as large as possible.
4. i = min(a//p//d, 256).
5. n = a//p - i*d and we don't find a solution if 511 < n.

Now let's attempt to get the lower bound 1264/465.

1. d = b = 465
2. p is at most min(64, a//d//2) = 1. Note that a//d is the integer component of a/b which will usually be the integer part of fout/fin. So the upper bound on p in practice is simply min(int(fin/fout/2), 64).
3. p = 1.
4. i = min(a//p//d, 256) = min(1264//1//465, 256) = min(2, 256) = 2
5. n = a//p - i*d = 1264//1 - 2*465 = 334

We don't have to look at the next lower bound because this one worked.

Now let's try to get the upper bound 193/71.

1. d = b = 71
2. p is at most min(64, a//d//2) = 1.
3. p = 1.
4. i = min(a//p//d, 256) = min(193//1//71, 256) = 2.
5. n = a//p - i*d = 193//1//71 - 2*71 = 51

We don't have to look at the next upper bound because this one worked.

Side note: The cases where you have to check multiple bounds to get an answer will be when 1 < p is required. So 256 < fin/fout.

I know this was a lot of writing. But a calculation that can be written out by hand like this should be able to be executed on your microcontroller well within the desired millisecond.

---

OK, the above did not always work at the top end of the range. Here is a fix.

First of all, p only has a maximum of 64 values. So only a finite and generally small list of values for p will leave us with

2 <= fin / fout / p < 257

For each of them we have i = int(fin / fout / p). If we have an integer then n = 0 and d = 1. Else we just need to find n/d by walking the Stern-Brocot tree as above.

So for each of at most 64 values of p get a candidate upper and lower candidate. And then pick the best one as your final answer.

If that's too slow, there are possible speedups with either continued fractions, or a binary search past the long stretches of < x over and over again. But I suspect that the complexity won't pay for itself with such small numbers.

Answer 2

OK I see it like this (C++):

//---------------------------------------------------------------------------
struct divider_cfg
    {
    // fout = (fin / (i + n/d)) / p
    unsigned int i; // integer part 2..256
    unsigned int n; // dividend 0..511
    unsigned int d; // divisor 1..512
    unsigned int p; // post-scaler 1..64
    };
//---------------------------------------------------------------------------
float get_fout(float fin,const divider_cfg &cfg)
    {
    float d;
    d=cfg.n; d/=cfg.d; d+=cfg.i; d/=cfg.p;
    return fin/d;
    }
//---------------------------------------------------------------------------
void get_cfg_f32(float fin,float fout,divider_cfg &cfg)
    {
    // fout = (fin / (i + n/d)) / p
    float d0,d1,d,dd;
    divider_cfg c;
    d0=fin/fout;                            // wanted divider
    dd=512.0;                               // initial error (bigger than possible)
    for (c.p=1;c.p<=64;c.p++)               // brute force p
        {
        d1=d0/c.p; c.i=floor(d1); d1-=c.i;  // compute i from p
        for (c.d=1;c.d<=512;c.d++)          // brute force d
            {
            c.n=floor(d1*c.d);              // compute n from d
            d=c.n; d/=c.d; d+=c.i; d/=c.p;  // compute divider
            d=fabs(d-d0);                   // compute error
            if (d<dd){ dd=d; cfg=c; }       // remember better config
            }
        }
    }
//---------------------------------------------------------------------------
// this enables float computation of the wanted divider for better precision
// the search itself is not using float so it does not affect speed too much

    #define _use_float
//---------------------------------------------------------------------------
#ifdef _use_float
void get_cfg_u32(float fin,float fout,divider_cfg &cfg)
#else
void get_cfg_u32(unsigned int fin,unsigned int fout,divider_cfg &cfg)
#endif
    {
    // fout = (fin / (i + n/d)) / p
    // (i + n/d)*p = fin/fout
    // (i + n/d) = (fin/fout)/p
    unsigned int idn,d0,d,dd,sh;
    divider_cfg c;

    #ifndef _use_float
    // d0 = wanted divider no FPU
    // lower precision on 32 bits need more
    unsigned int di,df;
    sh=32-9;                                // fixed point shift: log2(max_unsigned_int) - 9
    di=fin/fout;
    df=fin%fout;
    for (d=0,dd=1;dd<=df;d++,dd<<=1);       // d = ceil(log2(df))
    d+=sh;                                  // number of bits required for (df<<sh)
    if (d<=32) d=0; else d-=32;             // d number of rounded bits to prevent overflow of (df<<sh)
    d0=(di<<sh) + ((df<<(sh-d))/(fout>>d)); // wanted divider   (fin<<sh)/fout
    #else
    // d0 = wanted divider using FPU
    sh=32-9;                                // fixed point shift: log2(max_unsigned_int) - 9
    float fd=float(fin)/float(fout);
    fd*=1<<sh; d0=fd;                       // wanted divider   (fin<<sh)/fout
    #endif

    dd=0xFFFFFFFF;                          // initial error (biggest possible)
    for (c.p=1;c.p<=64;c.p++)               // brute force p
     for (c.d=1;c.d<=512;c.d++)             // brute force d
        {
        if (c.d==508)
         c.d=c.d;

        d=d0/c.p;                           //  = (i<<sh) + ((n<<sh)/d) = (fin<<sh)/(fout*p)
        c.i=d>>sh;
        d-=c.i<<sh;                         //  = ((n<<sh)/d)
        c.n=(d*c.d)>>sh;

        d=((c.i<<sh)+((c.n<<sh)/c.d))/c.p;  // compute divider << sh
        d=abs(d-d0);                        // compute error
        if (d<dd){ dd=d; cfg=c; }           // remember better config
        }
    }
#undef _use_float
//---------------------------------------------------------------------------

the integer version assumes 32bit unsigned integer type... both versions brute force p,d and recomputes the rest remembering better solution. The integer version uses configuration #define _use_float which enables more precise divider computation, in case your platform does not have float at all or you have more than 32bit integers you can use the fully integer version (just rewrite all the 32 to whatever bitwidth you have)

ON my old PC its:

              fin  = 16000000.00
              fout = 3005000.00
[  10.696 ms] fout (f32) =   3005002.25 err(        2.25) f32 | ((  5)+(134/413))/  1
[   1.305 ms] fout (i32) =   3005002.25 err(        2.25) u32 | ((  5)+(134/413))/  1

however on MCU it will be slower of coarse and without FPU I would definately go for the u32 (fixed point) version (its precision using the _use_float is the same as float version anyway).

The sh is precision (number of bits after decimal point) of fixed point math. I made changes so overflow should not occur anymore.