4
\$\begingroup\$

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
  1. Calculate the errors : $$ e_{+} := |f_{\text{set}} - f_{+}|,\quad e_{-} := |f_{\text{set}} - f_{-}| $$

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

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

enter image description here

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:

enter image description here

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.

\$\endgroup\$
2
  • \$\begingroup\$ If N cannot be 1, then your bounds statement at the beginning is incorrect. \$\endgroup\$
    – Reinderien
    Aug 28, 2021 at 0:46
  • \$\begingroup\$ @Reinderien, You're right. But the word "cannot" here is equal to "useless" because formally it can be 1 but in this case what the D can be : either 0 or 1. In this case the module just don't generate anything :) \$\endgroup\$
    – LRDPRDX
    Aug 28, 2021 at 12:15

2 Answers 2

1
\$\begingroup\$

@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.

\$\endgroup\$
3
  • \$\begingroup\$ First of all, thank you for your answer. 1) I see you define the error as a f^2 + D^2. But I actually don't mind about precision of the duty cycle. More over, it is not obvious that the frequency and duty cycle have to be included "symmetrically" in the error formula. 2) I agree that your approach to run through all ranges followed by picking the closest frequency is better here. 3) Yes, I think you're right about the actual values of ranges -- it makes sense, more over, I see this on my oscilloscope. \$\endgroup\$
    – LRDPRDX
    Aug 28, 2021 at 7:50
  • \$\begingroup\$ BTW, I think that Python is redundant here. You could explain all of this stuff in words and formulas. \$\endgroup\$
    – LRDPRDX
    Aug 28, 2021 at 7:55
  • \$\begingroup\$ I guess I could do that? But it wouldn't be as useful of an answer. Having actual code execute an alternate algorithm is somewhat far from redundant, particularly since I wouldn't be confident in my results otherwise. \$\endgroup\$
    – Reinderien
    Aug 28, 2021 at 8:36
2
\$\begingroup\$

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!

\$\endgroup\$
1
  • 1
    \$\begingroup\$ Thank you for your answer, I'll take it into account. \$\endgroup\$
    – LRDPRDX
    Aug 28, 2021 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.