Why NaN from FOS here?

Questions
23

I could see uses for that – not every synth needs higher precision, but some synths would benefit from it.

Btw I think it may be wrong for FOS/SOS to call _next_1() at the end of the Ctor. Feedforward and feedback delay samples should be zero at first, shouldn’t they? It looks to me like we’re calling the function because “the function is supposed to be called during init” but I personally have found more than my fair share of edge cases where weird things happen because of the reflexive application of that code pattern.

hjh

24

I didn’t put that much thinking into it, you may be right. The current implementation just checks buffer length and input rates. I was just suggesting supplementing this aspect. Yes, maybe some unpredictable behavior comes from what you said.

25

The next call in the Ctor isn’t your issue btw – that’s in the sources now, and probably dates back to SC2. I’m just wondering if JMc did that just because of following a design (“ctor should call next_1”) or if it was really thought through. I suspect the former (since there are a lot of UGens).

hjh

26

I didn’t know it was a common practice. (For me, it kind of makes sense) The rationale is to avoid undefined variables. But how it creates instability would require a bit more study and time, I don’t have now.

27

Don’t suppose either of you would mind writing a one or two sentence summary of the cause of the issue? I’m not very DSP literate and this is a long nuanced thread dealing with a pretty big problem - nans in the audio stream are infectious.

28

Jordan,

I listed 3 main reason that those things can happen. That’s not much to it. It’s probably an initialization problem of a filter. The problem is that the equation requires preceding values, which do not exist. So, one can either just not set them, or one can set them to zero, or one can in the initialization function already do a first calculation. The real-world implementation and workaround is different than just the theory. James is suspicious of a common solution. I propose other things to prevent it, and I don’t know if his suspicious is the problem or not.

29

There are a couple of intertwined issues here, which are hard to separate.

NaN in this case is coming from Sweep outputting 0 as its pre-sample, and then reciprocal produces NaN (unlike 1/0 = inf).

The comments about safe initialization of filters are useful, but I think perhaps a bit tangential.

Now… why is Sweep’s pre-sample 0, when its first real output value is nonzero?

Checking the source code, I think it’s just a bit of sloppiness (which is perhaps inevitable, since JMc was just one person writing code for hundreds of UGens). Sweep’s documentation says “Starts a linear raise by rate/sec from zero” – I overlooked this at first. But the calculation functions increment first before outputting – so the zero will get swallowed. Pre-incrementing already makes the actual behavior deviate from the documented intention.

But then, in the constructor,

    ZOUT0(0) = unit->mLevel = 0.f;

Huh…? But we already know from the next functions that the output will be pre-incremented. So setting ZOUT(0) to 0 is… baffling? A mistake? Or should the next functions instead read like LOOP1(inNumSamples, ZXP(out) = level; level += rate;); ? So that the initial 0 will be output as a normal sample.

Because if Sweep’s behavior were consistent with its documentation, then I would have known to write (Sweep.ar(0, SampleRate.ir) + 1).reciprocal, and done, no problem.

Unfortunately, I suspect there are a lot of such cases that we haven’t even found yet. As far as I can see, there was never a clear decision about the meaning of the Ctor pre-sample. Some UGens (like Sweep) initialize to y[-1] and then output y[0] in the first cycle; other UGens initialize to y[0]… which basically reads like, coding shortcuts were taken under time pressure.

So a short summary, I guess, would be:

  1. FOS is fine.
  2. Sweep should not mislead users into thinking that the output will always be > 0.
  3. UGen initialization is a bit sloppy in places.

hjh

2 Likes
31

Ok, sorry, I think our eyes were focusing different points to prevent it. The original code was buggy because of Sweep, and I was trying to “protect” the filter to go out of control because of nan/inf/anything bad.

Since the check for nan and inf would happen just once, it could be useful as part of initialization. But I may have expanded the topic a bit, it’s true. Food for future thought when discussing filter implementations.

32

Last word, just for reference. I’ve implemented some of the checks in c++, for my tinysynth engine. Just a prototype, checks on initialization for now, but theoretically it could check more times, I just need to figure out how to do it outside the audio thread. But may be useful in general, since it’s c++ not Haskell. (the design is minimalist, but works just fine, and with template to use any type; it also implements SIMD operations [that means, at least in my laptop, 8 samples at a time] ). Of course, needs + work when I have time.


// Base class for some filters
template <std::floating_point FloatType, std::size_t NUM_OF_FRAMES>
class filter_base {
public:
  virtual ~filter_base() = default;

  using batch_type = xsimd::batch<FloatType>;
  static constexpr std::size_t batch_size = batch_type::size;

  virtual void set_parameters(FloatType freq, FloatType Q) noexcept = 0;

  bool is_stable_simple() const noexcept {
    return (xsimd::abs(m_a1) < FloatType(1) + m_a2) &&
           (xsimd::abs(m_a2) < FloatType(1));
  }

  bool are_poles_stable() const noexcept {
    std::complex<FloatType> pole1, pole2;
    calculate_poles(pole1, pole2);
    return (std::abs(pole1) < FloatType(1)) && (std::abs(pole2) < FloatType(1));
  }

  std::complex<FloatType> frequency_response(FloatType freq) const noexcept {
    FloatType omega =
        FloatType(math_constants::two_pi_c<FloatType>) * freq / m_sample_rate;
    std::complex<FloatType> z(xsimd::cos(omega), -xsimd::sin(omega));
    std::complex<FloatType> numerator = m_b0 + m_b1 * z + m_b2 * z * z;
    std::complex<FloatType> denominator =
        FloatType(1) + m_a1 * z + m_a2 * z * z;
    return numerator / denominator;
  }

protected:
  alignas(32) FloatType m_a0 = FloatType(1), m_a1 = FloatType(0),
                        m_a2 = FloatType(0);
  alignas(32) FloatType m_b0 = FloatType(0), m_b1 = FloatType(0),
                        m_b2 = FloatType(0);
  FloatType m_sample_rate = FloatType(44100);

private:
  void calculate_poles(std::complex<FloatType> &pole1,
                       std::complex<FloatType> &pole2) const noexcept {
    FloatType discriminant = m_a1 * m_a1 - FloatType(4) * m_a0 * m_a2;
    FloatType two_a = FloatType(2) * m_a0;
    FloatType sqrt_disc = std::sqrt(std::abs(discriminant));
    FloatType real_part = -m_a1 / two_a;
    FloatType imag_part = sqrt_disc / two_a;

    if (discriminant >= FloatType(0)) {
      pole1 = std::complex<FloatType>(real_part + imag_part, FloatType(0));
      pole2 = std::complex<FloatType>(real_part - imag_part, FloatType(0));
    } else {
      pole1 = std::complex<FloatType>(real_part, imag_part);
      pole2 = std::complex<FloatType>(real_part, -imag_part);
    }
  }
};
33

I’m just going to highlight this for the attention of @mike, who has spent quite a bit of time on resolving initial sample issues across a number of UGens.

1 Like
34

@joslloand @mike I’d like to know more about those enhancements [maybe another topic]. Maybe some ideas coincide

35

I’ve been been chipping away at fixing UGen initialization in fits and starts, feel free to read over the RFC for background, some clarification of terms, and the approach to fixing initialization. I also track the progress here.

WRT to this topic specifically…

according to my notes I’ve reviewed FOS and SOS during a first pass, but haven’t finalized it in a PR. A quick look at that project branch, here was my solution:

void FOS_Ctor(FOS* unit) {
    if (unit->mBufLength == 1) {
        SETCALC(FOS_next_1);
    } else {
        if (INRATE(1) == calc_FullRate && INRATE(2) == calc_FullRate && INRATE(3) == calc_FullRate) {
            SETCALC(FOS_next_a);
        } else {
            SETCALC(FOS_next_k);
        }
    };
    unit->m_y1 = 0.f;
    unit->m_a0 = ZIN0(1);
    unit->m_a1 = ZIN0(2);
    unit->m_b1 = ZIN0(3);

    FOS_next_1(unit, 1);
    unit->m_y1 = 0.f; // the fix
}

I could bump this and SOS to the top of the list for the next round of PRs if folks would find it useful to prioritize this. I’m currently taking some time off, but I can swing back around to it in the near future…