Find the closest possible frequency of a pulser that uses limited number of fixed time intervals to generate the signal (CAEN V2718)

Question

Disclaimer

I've created this post not only for the code review, but rather because I think it would be helpful for the users of the CAEN V2718 module.

Problem description

I have a VME board (V2718) which may be used as a square form generator. To set the frequency and duty cycle user must specify three parameters :

1. The period N (in terms of the number of time-steps) : $$ 1 \leq N \leq 255 \tag{1a}\label{1a} $$

2. The duty cycle D (in terms of the number of time-steps) : $$ 1 \leq D \leq N-1 \tag{1b}\label{1b} $$

3. The time-step d (used to express the period and duty cycle) which may be either of the following : $$ d_1 = 25n\text{s},\, d_2 = 1600n\text{s},\, d_3 = 410\mu\text{s},\, d_4 = 104m\text{s} $$

For example, to generate the 100 Hz square (50% d.c.) one can choose the following configuration: $$ d = d_3,\,N = 24,\, D=N/2\quad \Rightarrow\quad T = d N\approx 0.01\text{ s} $$

Usually the primary quantities are frequency and duty cycle --- not period and duty cycle. And, of course, people don't want to specify some strange parameters like N, D, and d --- they want frequencies and duties! So the problem is to write a function that chooses the raw parameters to generate the desirable square waveform.

Algorithm

To find the closest possible frequency to the desirable frequency f_set we need to find the closest possible period corresponding to the f_set. Given f_set the algorithm is the following

1. Find two T- and T+ which can be obtained by the module such that : $$ T_{set} \in [T_{-}, T_{+}) $$

• If either of T's is absent then let the absent one be equal to the present one
• If both of T's are absent then return false

2. Calculate the errors : $$ e_{+} := |f_{\text{set}} - f_{+}|,\quad e_{-} := |f_{\text{set}} - f_{-}| $$

3. Choose the corresponding T based on the errors' values.

The picture shows the idea (although see the Useful notes):

Useful notes:

• N cannot be 1. See(\ref{1b})
• T+ and T- can be expressed by using different time-steps (in other words, the red triangles on the picture may be the ticks of the different time scales)

Code

bool V2718Pulser::SetSquare( uint32_t freq, uint8_t duty )
{
    struct { double expo; double num; CVTimeUnits unit; } ss[4] = { { 1000000000., 25., cvUnit25ns },
                                                                    { 10000000.,   16., cvUnit1600ns },
                                                                    { 100000.,     41., cvUnit410us },
                                                                    { 1000.,      104., cvUnit104ms } };

    const uint32_t MAX_PERIOD = 0xff; // MAX_N
    const uint32_t MIN_PERIOD = 0x02; // MIN_N

    struct { uint32_t n; int u; } sPlus, sMinus; // T+, T-

    uint32_t n0 = 0;// The reference point

    if( freq > 0 )
    {
        for( int i = 0; i < 4; ++i )
        {
            n0 = std::floor( ss[i].expo / ss[i].num / freq );

            if( (n0 >= MIN_PERIOD) && (n0 < MAX_PERIOD) )
            {
                // T_set falls between the ticks of the same time-scale
                sMinus = { n0,     i };
                sPlus  = { n0 + 1, i };
                break;
            }
            else if( n0 == 1 )
            {
                // Bad value because of the duty cycle --- should be changed
                if( i > 0 )
                {
                    sMinus = { MAX_PERIOD, i - 1 }; // the last tick of the finer time-scale
                    sPlus  = { n0 + 1,     i     }; // use the current time-scale, but the next tick
                }
                else
                {
                    // There is no finer time-scale
                    sPlus = sMinus = { n0 + 1, i };
                }
                break;
            }
            else if( n0 == MAX_PERIOD )
            {
                if( i < 3 )
                {
                    sMinus = { n0,         i     }; 
                    sPlus  = { MIN_PERIOD, i + 1 }; // use the 2-th tick of the coarser time-scale
                }
                else
                {
                    // The is no coarser time-scale
                    sPlus = sMinus = { n0, i };
                }
                break;
            }
        }

        if( n0 != 0 )
        {
            double errorPlus    = std::fabs( freq - ss[sPlus.u].expo / ss[sPlus.u].num / sPlus.n );
            double errorMinus   = std::fabs( freq - ss[sMinus.u].expo / ss[sMinus.u].num / sMinus.n );
            if( errorPlus < errorMinus )
            {
                n0 = sPlus.n;
                fTimeUnit = ss[sPlus.u].unit; // d
            }
            else
            {
                n0 = sMinus.n;
                fTimeUnit = ss[sMinus.u].unit; // d
            }

            duty = ((duty > 0) ? ((duty < 100) ? duty : 99) : 1);
            uint32_t width = n0 * duty / 100;// < MAX_PERIOD
            width = (width > 0) ? width : 1;

            if( width < n0 )
            {
                fPeriod = n0; // N
                fWidth = width; // D

                return true;
            }
        }
    }

    return false;
}

Result

I created the "frequency error graph" of this module for the frequencies from 1 to 10000000 Hz:

Bad frequencies

You can see that there are frequencies with a very big error (~30%). One of such frequency is 1836 Hz. Indeed,

$$ T_{\text{bad}} = \frac{1}{1836} = 0.000544662 = 544.662 \times 10^{-6}\text{s} = 544.662 \mu\text{s} $$

which means the closest times are $$ T_{-} = 255 \times 1600n\text{s} = 408\mu\text{s},\quad T_{+} = 2 \times 410\mu\text{s} = 820\mu\text{s} $$

so

$$ f_{-} = 2450.98\text{ Hz},\quad f_{+} = 1219.51\text{ Hz} $$

Very bad frequencies

There are 2 frequencies (2440, 2441) that even "don't exist" (the function returns false), i.e. neither T+ nor T- were found. Indeed,

$$ T_{\text{very bad}} = \frac{1}{2440} = 0.000409836 = 409.836\mu\text{s} $$

and

$$ 255 \times 1600n\text{s} < T_{\text{very bad}} < 410\mu\text{s} $$

so

$$ \left\lfloor \frac{T_{\text{very bad}}}{d_{i}} \right\rfloor \text{is either 0 or } >255 $$

and the algorithm fails to find T's.

Answer 1

@JDługosz has you covered in terms of C++ syntax. Algorithmically, a quick bit of Python illustrates where your method can be improved. Basically your iteration through all four timescales is correct, but your plus/minus logic is not, leading to the degenerate cases you've already identified. You don't need any explicit conditions at all, and can run through all of the timescales with defined minima and maxima, accounting for error in each case:

from typing import NamedTuple, Iterable, Tuple

d_opts = (25e-9, 1.6e-6, 410e-6, 104e-3)


class Config(NamedTuple):
    f: float
    f_act: float
    f_err: float
    duty: float
    duty_act: float
    duty_err: float
    N: int
    D: int
    d: float

    @classmethod
    def approximate(cls, f: float, duty: float, d: float) -> 'Config':
        # Forward
        N_float = 1 / f / d
        N = max(2, min(255, round(N_float)))
        D_float = duty * N
        D = max(1, min(N-1, round(D_float)))

        # Backward
        f_act = 1 / N / d
        duty_act = D / N

        return cls(
            f=f, f_act=f_act, f_err=f_act/f - 1,
            duty=duty, duty_act=duty_act, duty_err=duty_act/duty - 1,
            N=N, D=D, d=d,
        )

    @classmethod
    def all_approx(cls, f: float, duty: float) -> Iterable[Tuple[
        float,  # total error
        'Config',
    ]]:
        for d in d_opts:
            config = cls.approximate(f, duty, d)
            err = config.f_err**2 + config.duty_err**2
            yield err, config

    @classmethod
    def closest(cls, f: float, duty: float) -> 'Config':
        err, config = min(cls.all_approx(f, duty))
        return config

    def __str__(self) -> str:
        return (
            f'{self.f:7g} Hz @ {self.duty:.0%} ~ '
            f'{self.f_act:9.5g} Hz @ {self.duty_act:.2%}, '
            f'Δ {self.f_err:+.1e}, {self.duty_err:+.1e}, '
            f'N={self.N:3}, D={self.D:3}, d={self.d:.2e} s'
        )


print('A selection of frequencies:')
for f, duty in (
    (100, 0.5),
    (1836, 0.5),
    (2440, 0.5),
    (2441, 0.5),
    (96e6, 0.37),
):
    print(Config.closest(f, duty))
print()

print('One frequency, best to worst config:')
for err, config in sorted(Config.all_approx(f=2441, duty=0.5)):
    print(config)

outputs much happier results for your degenerate cases:

A selection of frequencies:
    100 Hz @ 50% ~    101.63 Hz @ 50.00%, Δ +1.6e-02, +0.0e+00, N= 24, D= 12, d=4.10e-04 s
   1836 Hz @ 50% ~      2451 Hz @ 50.20%, Δ +3.3e-01, +3.9e-03, N=255, D=128, d=1.60e-06 s
   2440 Hz @ 50% ~      2451 Hz @ 50.20%, Δ +4.5e-03, +3.9e-03, N=255, D=128, d=1.60e-06 s
   2441 Hz @ 50% ~      2451 Hz @ 50.20%, Δ +4.1e-03, +3.9e-03, N=255, D=128, d=1.60e-06 s
9.6e+07 Hz @ 37% ~     2e+07 Hz @ 50.00%, Δ -7.9e-01, +3.5e-01, N=  2, D=  1, d=2.50e-08 s

One frequency, best to worst config:
   2441 Hz @ 50% ~      2451 Hz @ 50.20%, Δ +4.1e-03, +3.9e-03, N=255, D=128, d=1.60e-06 s
   2441 Hz @ 50% ~    1219.5 Hz @ 50.00%, Δ -5.0e-01, +0.0e+00, N=  2, D=  1, d=4.10e-04 s
   2441 Hz @ 50% ~    4.8077 Hz @ 50.00%, Δ -1.0e+00, +0.0e+00, N=  2, D=  1, d=1.04e-01 s
   2441 Hz @ 50% ~ 1.5686e+05 Hz @ 50.20%, Δ +6.3e+01, +3.9e-03, N=255, D=128, d=2.50e-08 s

Hardware deficiencies

A more stable algorithm will fix the cases where no configuration was found; but for some worst-case values like 1836 Hz nothing can be done because the hardware doesn't have very good resolution. 1836 is in fact not a coincidence; it's the average of the gap between the second and third ranges. The two worst errors will thus be found at these frequencies (assuming that the original range values are lies and are actually exact multiples of 256 that the vendor has rounded in their documentation):

$$ \frac 1 2 \left( \frac 1 {255 \times 1.6 \times 10^{-6}} + \frac 1 {2 \times 409.6 \times 10^{-6}} \right) \approx 1835.842 $$

$$ \frac 1 2 \left( \frac 1 {255 \times 409.6 \times 10^{-6}} + \frac 1 {2 \times 0.1048576} \right) \approx 7.171257 $$

In each case you're forced to choose between -33% error or +33% error based on the higher or lower range.

Answer 2

ss is neither static nor const, and a stack-based variable will be initialized each time the function is called. It appears that you don't modify this array at all but only refer to it. If that's correct, make it constexpr.

const uint32_t MAX_PERIOD = 0xff; // MAX_N
const uint32_t MIN_PERIOD = 0x02; // MIN_N

With a modern compiler, use constexpr generally for this. const variables might be compile-time constants if the initializer is a constant expression, but there's no warning if that's not the case. The newer constexpr is more demanding and means that it will be a compile-time constant.

sPlus, sMinus; // T+, T-

The name is confusing since it doesn't match the comment: S+ vs T+.

I see they are uninitialized, and the following rather large and complexly-nested block seems to have the sole purpose of setting these values. Break that out into its own function, and then as well as removing a "meander" from the code, you can write them as initialized constants:

const auto [n0, sPlus, sMinus] = find_T();

That will make the function much easier to understand too; it's not easy to see where the "real code" is after the nested blocks that set up these values.

if( n0 != 0 )
{
   // entire meat of the function goes here
   return true;
}
return false;

Both the test for n0 and the earlier test for freq should be written as preconditions. Put reverse the sense of the test and return false if it fails the condition. The body of the function is then the "main line" and not under a conditional. Since this is nested, you'll "write to the left" quite a bit more in this case!