party.py

#!/usr/bin/python3
# The MIT License (MIT)
# 
# Copyright (c) 2020 Jeff Epler
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
import argparse
import fractions
import functools
import math

# Find all the ways to partition 'seq' into two subsequences a and b.
# Arbitrarily, the first item is always in the 'a' sequence, so e.g.,
# the partitions of 2 items [p, q] are just 2: [[p], [q]] and [[p, q], []].
# the partitions [[q], [p]] and [[], [p, q]] will not be returned.
#
# The empty sequence results in the single partition [[], []]
def partitions(seq):
    seq = list(seq)

    # Special case: empty sequence
    if len(seq) == 0:
        yield [], []

    else:
        for j in range(1, 2**len(seq), 2):
            yield (
                [si for i, si in enumerate(seq) if j & (1<<i)],
                [si for i, si in enumerate(seq) if (~j) & (1<<i)]
            )

# Convert various representations to fractions
# besides what fractions.Fraction will parse, you can write
#  - underscores within numbers, as place separators
#        1_000_000
#  - Mixed fractions
#        1+3/4  or 1 3/4 (both equal to 7/4)
#        or 1-3/4 (equal to 1/4)
#  - Exponentials in any part
#        1e3 5_000e3/1e9
#  - Common suffixes: M, G, K, m, u
def F(n):
    if isinstance(n, str):
        n = n.replace("_", "")
        n = n.replace("+", " ")
        n = n.replace("-", " -")
        n = n.replace("m", "e-3")
        n = n.replace("u", "e-6")
        n = n.replace("k", "e3")
        n = n.replace("M", "e6")
        n = n.replace("G", "e9")
        if ' ' in n:  # Accomodate 1 1/3, 1+1/3
            w, n = n.rsplit(None, 1)
        else:
            w = 0
        if '/' in n:
            n, d = n.split('/')
        else:
            d = 1
        return F(w) + fractions.Fraction(n) / fractions.Fraction(d)
    return fractions.Fraction(n)

def flcm(a, b):
    a = F(a)
    b = F(b)
    p = a.numerator * b.denominator
    q = b.numerator * a.denominator
    r = math.gcd(p, q)
    d = a.denominator * b.denominator
    return fractions.Fraction(p * q, r * d)

def err_str(x):
    if x == 0: return "0"
    if x < 1e-12: return "%g" % x
    if x < 1e-9: return "%.3fppt" % (x * 1e12)
    if x < 1e-6: return "%.3fppb" % (x * 1e9)
    return "%.3fppm" % (x * 1e6)

def place_in_range(freq, f_low):
    while freq < f_low / 2048:
        freq *= 2
    return freq

def ilog2(x):
    j = 0
    while x > 1:
        j += 1
        x /= 2
    return j

def calculate_freq(clocks, f_low, f_high):
    freq = functools.reduce(flcm, clocks, 1)
    n = (f_low + freq - 1) // freq
    r = freq * n
    return r

(
    INTEGER_MULTIPLIER_DOUBLE_INTEGER_DIVISOR,
    FRACTIONAL_MULTIPLIER_DOUBLE_INTEGER_DIVISOR,
    INTEGER_MULTIPLIER_INTEGER_DIVISOR,
    FRACTIONAL_MULTIPLIER_INTEGER_DIVISOR,
    INTEGER_MULTIPLIER_FRACTIONAL_DIVISOR,
    FRACTIONAL_MULTIPLIER_FRACTIONAL_DIVISOR,
) = range(6)

plan_names = {
    INTEGER_MULTIPLIER_DOUBLE_INTEGER_DIVISOR: "Integer multiplier, double integer divisor",
    FRACTIONAL_MULTIPLIER_DOUBLE_INTEGER_DIVISOR: "Fractional multiplier, double integer divisor",
    INTEGER_MULTIPLIER_INTEGER_DIVISOR: "Integer multiplier, integer divisior",
    FRACTIONAL_MULTIPLIER_INTEGER_DIVISOR: "Fractional multiplier, integer divisor",
    INTEGER_MULTIPLIER_FRACTIONAL_DIVISOR: "Integer multiplier, fractional divisor",
    FRACTIONAL_MULTIPLIER_FRACTIONAL_DIVISOR: "Fractional multiplier, fractional divisor",
}

MAX_DENOM = 1048575

class Pll:
    def __init__(self, f_in, clocks, *, f_low=600_000_000, f_high=900_000_000):
        self.f_in = F(f_in)
        self.f_low = F(f_low)
        self.f_high = F(f_high)

        self.clocks = [F(c) for c in clocks]

        self.setting_type, self.multiplier = self._calculate()

    def _calculate(self):
        f_in = self.f_in
        clocks = [place_in_range(c, self.f_low) for c in self.clocks]
        clocks2 = [2*c for c in clocks]

        f = calculate_freq(clocks2 + [f_in], self.f_low, self.f_high)
        if f < self.f_high:
            return INTEGER_MULTIPLIER_DOUBLE_INTEGER_DIVISOR, f / f_in

        clocks = [2*c for c in clocks]
        f = calculate_freq(clocks2, self.f_low, self.f_high)
        if f < self.f_high and (f / f_in).denominator <= MAX_DENOM:
            return FRACTIONAL_MULTIPLIER_DOUBLE_INTEGER_DIVISOR, f / f_in

        f = calculate_freq(clocks + [f_in], self.f_low, self.f_high)
        if f < self.f_high:
            return INTEGER_MULTIPLIER_INTEGER_DIVISOR, f / f_in

        while True:
            clocks.pop()
            if not clocks:
                break

            f = calculate_freq(clocks + [f_in], self.f_low, self.f_high)
            if f < self.f_high:
                return INTEGER_MULTIPLIER_INTEGER_DIVISOR, f / f_in

            f = calculate_freq(clocks, self.f_low, self.f_high)
            if f < self.f_high and (f / f_in).denominator <= MAX_DENOM:
                return FRACTIONAL_MULTIPLIER_FRACTIONAL_DIVISOR, f / f_in

        ratio = (self.f_low / f_in).limit_denominator(MAX_DENOM)
        
        return FRACTIONAL_MULTIPLIER_FRACTIONAL_DIVISOR, ratio

    def exact_divider(self, clock):
        clock2 = place_in_range(clock, self.f_low)
        return ((self.f_out / clock2), ilog2(clock2 / clock))

    def divider(self, clock):
        de, r = self.exact_divider(clock)
        d = de.limit_denominator(MAX_DENOM)
        return d, r

    def error(self, clock):
        de, r = self.exact_divider(clock)
        d = de.limit_denominator(MAX_DENOM)
        return (de - d) / de
        
    @property
    def f_out(self):
        return self.f_in * self.multiplier

    @property
    def dividers(self):
        return [self.divider(c) for c in self.clocks]

    @property
    def errors(self):
        return [self.error(c) for c in self.clocks]

    @property
    def score(self):
        if not self.clocks: return 0
        e = len([e for e in self.errors if e == 0])
        did = len([d for d, r in self.dividers if d.denominator == 1 and d.numerator % 2 == 0])
        sid = len([d for d, r in self.dividers if d.denominator == 1])
        lc = len(self.clocks)

        return 6-self.setting_type + fractions.Fraction(e, lc) + fractions.Fraction(did, lc*lc) + fractions.Fraction(did, lc**3)

    def print(self):
        print(f"Frequency plan type: {plan_names[self.setting_type]} ({self.setting_type})")
        print(f"Multiplier = {self.multiplier} ({float(self.multiplier)})")
        print(f"Intermediate frequency = {self.f_out} ({float(self.f_out)})")
        print()
        for c in self.clocks:
            d, r = self.divider(c)
            e = self.error(c)
            print(f"Desired output frequency: {c} ({float(c)})")
            print(f"Divider = {d} ({float(d)})")
            if e == 0:
                print("Exact")
            else:
                c_actual = self.f_out / d / (2**r)
                print(f"Actual output frequency: {c_actual} ({float(c_actual)})")
                print(f"Relative Error: {err_str(e)}")
                print(f"Absolute Error: {float(c - c_actual):.3g}Hz")
            print(f"r_divider = {r} (/ {2**r})")
            print()
 
def plan(f_in, c1, c2):
    p = Pll(f_in, c1)
    q = Pll(f_in, c2)
    return (p.score + q.score, p, q)

parser = argparse.ArgumentParser(
    description='Create a frequency plan for Si5351')
parser.add_argument('frequencies', metavar='freq', nargs='+',
    type=F, help='Integer, ratio, or decimal frequencies')
parser.add_argument('--input-frequency', '-i', metavar='freq',
     default=25000000, type=F,
    help='Input frequency')
args = parser.parse_args()
f_in = args.input_frequency

score, p, q = max((plan(f_in, p, q) for p, q in partitions(args.frequencies)), key=lambda s: s[0])

print(f"Input frequency: {f_in} ({float(f_in)})")
print(f"Frequency plan score: {float(score):.2f}")
print()
print("PLL A:")
p.print()

if q.clocks:
    print()
    print("PLL B:")
    q.print()