Why NaN from FOS here?

Questions
9

I’d expect inf here…? 0/0 is NaN; 1/0 is inf, IIRC.

But I tested Sweep’s output and y(0) from Sweep = the increment and not 0. { var trig = Impulse.ar(0); Poll.ar(trig, Sweep.ar(0, SampleRate.ir)); FreeSelf.kr(trig <= 0); Silent.ar(1) }.play; – I’m not at the computer now but I’m pretty sure this prints 1.0.

So if FOS is getting polluted by 1/0, the only way this could happen is if the FOS constructor is receiving y(-1) as the coefficient’s initial sample.

My opinion is that it is always a mistake when a UGen returns y(-1) as its pre-sample.

Some UGens’ pre-sample is y(0); others’ is y(-1). This means that downstream units can’t assume anything about the meaning of the initial value. The only sense in which this is a good thing is that y(-1) may be easier to produce, from a code perspective. But IMO this is not a good rationale. And in fact, there’s a whole slew of weird bugs resulting from this inconsistency in the UGen library.

hjh

1 Like
10

EDIT: isn’t it (1/0) * 0 ? I think that’s what I wanted to write

But the Nan propagates in the first samples so it’s all a mess.

My personal opinion is that filters could avoid the user to ser "bad coefficients’ and "bad initialization’.

On top of that, also apply that function I posted to avoid another cause of troubles.

No, I don’t think no SuperCollider user is always safe and “knows what they’re doing”. And SC interface has not hint about it.

Unless it’s important to let the user to set bad coefficients, optionally.

EDIT: Also, to give an example, if I’m experimenting with a filter knowing the coefficients are limited to the range that the filter is stable, I will explore much more its possibilities.

Don’t you think? For example, this website implements those filters and impose those limits. I think without it, I would not play with it in the same way: Filter Playground

11

Ah… found it.

{ FOS.ar(DC.ar(1), Sweep.ar(0, 1) > 0, [0, DC.ar(0)], [0, DC.ar(0)]) }.plot;

sc-fos-bug

The top line picks up 0 from Sweep’s pre-sample and interpolates from this value to 1.

void FOS_next_k(FOS* unit, int inNumSamples) {
    float* out = ZOUT(0);
    float* in = ZIN(0);
    float next_a0 = ZIN0(1);
    float next_a1 = ZIN0(2);
    float next_b1 = ZIN0(3);

    double y1 = unit->m_y1;
    double a0 = unit->m_a0;
    double a1 = unit->m_a1;
    double b1 = unit->m_b1;
    double a0_slope = CALCSLOPE(next_a0, a0);
    double a1_slope = CALCSLOPE(next_a1, a1);
    double b1_slope = CALCSLOPE(next_b1, b1);

    ...
}

But the bottom line shows that Sweep’s true output is always nonzero: Sweep.ar(0, 1) > 0 is always 1 for sample index >= 0. This means the top line is picking up y(-1).

I haven’t yet heard a rationale why pre-sample = y(-1) is good, except for “it’s too hard to change” (which is just a cop-out).

hjh

1 Like
12

FOS help:

Formula is equivalent to:

out(i) = (a0 * in(i)) + (a1 * in(i-1)) + (b1 * out(i-1))

FOS source:

    double y0 = in + b1 * y1;
    ZOUT0(0) = a0 * y0 + a1 * y1;
    y1 = y0;

These certainly don’t behave the same.

(
r = Routine { |inval|
	var avg = 0.0, n = 0, recip;
	var avg2 = 0.0, a0, a1 = 0, b1, y0, y1 = 0;
	loop {
		n = n + 1;
		a0 = n.reciprocal;
		b1 = 1 - a0;

		avg = (a0 * inval) + (b1 * avg);
		y0 = inval + (b1 * y1);
		avg2 = a0 * y0 + (a1 * y1);
		inval = [avg, avg2].yield;
		y1 = y0;
	}
};
)

10.do { r.next(1).postln };

[1.0, 1.0]
[1.0, 0.75]
[1.0, 0.66666666666667]
[1.0, 0.625]
[1.0, 0.6]
[1.0, 0.58333333333333]
[1.0, 0.57142857142857]
[1.0, 0.5625]
[1.0, 0.55555555555556]
[1.0, 0.55]

So “equivalent to” is flat wrong in the docs.

hjh

1 Like
14

James, I think some design would have to be changed to avoid those things. The SynthDef would have to be checked before initialization, because the calculation is… not quite that simple.

Do you know how that could be done with current model?

Just for illustration, this is a code in Haskell that calculates that. It involves complex numbers and some computations that may not be cheap.

EDIT: Would be nice to benchmark the cost of checking this in an optimized real-time code. Actually, it doesn’t need to check at audio rate, but much lower rates or demand rates

import Data.Complex (Complex((:+)))

data IIRParams = IIRParams
    { b0 :: Float
    , b1 :: Float
    , b2 :: Float
    , a0 :: Float
    , a1 :: Float
    , a2 :: Float
    }
    deriving (Show)

-- Example of an unstable filter (IIR BAD Coefficients)
unstableFilter :: IIRParams
unstableFilter = IIRParams
    { b0 = 0.5
    , b1 = 0.5
    , b2 = 0.5
    , a0 = 1.0
    , a1 = -1.8  -- These coefficients
    , a2 = 0.9   -- produce poles outside 'circle'
    }

-- Calculate poles of IIR filters
calculatePoles :: IIRParams -> (Complex Float, Complex Float)
calculatePoles (IIRParams _ _ _ a0 a1 a2) =
    let discriminant = a1 * a1 - 4 * a0 * a2
        twoA = 2 * a0
        sqrtDisc = sqrt discriminant
        realPart = (-a1) / twoA
        imagPart = sqrtDisc / twoA
        root1 = realPart :+ imagPart
        root2 = realPart :+ (-imagPart)
    in (root1, root2)


-- Magnitude of complex numbers
magnitude :: Complex Float -> Float
magnitude (r :+ i) = sqrt (r * r + i * i)

-- Check IIR filter is stable
isStable :: IIRParams -> Bool
isStable params =
    let (pole1, pole2) = calculatePoles params
    in magnitude pole1 < 1 && magnitude pole2 < 1

main :: IO ()
main = do
    let params = unstableFilter
    print params
    print $ calculatePoles params
    print $ isStable params  -- should be False

The theory comes from here: Audio EQ Cookbook

One needs to calculate two things:

  1. The calculatePoles finds the roots of the characteristic equation; the equation that represents the poles of the filter. The characteristic equation is the filter transfer function denominators setting them equal to zero, then multiplying this by z^2

$a_0 + a_1 z^{ -1} + a_2 z^{-2} = 0$ then multiply this by z^2 → $ a_0 z^2 + a_1 z + a_2 = 0 $
image ----> image

  1. The isStable function ensures both poles have magnitudes less than 1.
16

I appreciate the thorough approach above, but perhaps that’s for a different thread?

This thread is about a very specific case. The idea is to produce a running average – not a moving average (where SC already has the incorrectly-named RunningSum, which you can divide by the window size to get a moving average) – a true running average, which keeps accumulating indefinitely.

I was using FOS to try to achieve this, because a moving average can be implemented as a first-order recursive filter:

moving_avg(x(n)) = (1/n) * x(n) + (1 - 1/n) * y(n-1)

I.e., a0 = 1/n, b1 = 1 - (1/n).

The answer to the title question has two parts:

First, 1/0 and 0.reciprocal return different results in the server (while returning the same result in the language):

[0.reciprocal, 1/0];
-> [inf, inf]

s.boot;

(
f = { |sig|
	{
		var trig = Impulse.ar(0);
		Poll.ar(trig, sig.value);
		Silent.ar(1)
	}.play;
};

f.({ var x = DC.ar(0); [x.reciprocal, 1/x] });
)

UGen(UnaryOpUGen): -nan
UGen(BinaryOpUGen): inf

Second, that 0.reciprocal should never be entering FOS at all. The first sample, y(0), of Sweep.ar(trig, SampleRate.ir) is 1.0. Thus it can be entering FOS only in the Ctor. My opinion is that it’s a remaining bit of sloppiness in the SC plugin interfaces, that the Ctor “pre-sample” currently is either y(0) or y(-1), at the whim of each plug-in, which of course has no idea which value is expected downstream. The only way that you can have a consistent expectation downstream is if that initial value being produced is consistent.

The NaN’s here are not FOS’s fault. They’re Sweep’s fault.

THEN… having resolved that by forcing audio rate coefficients, which causes the 0.reciprocal to be discarded before it affects FOS results, I found that FOS’s behavior does not match the documented “equivalent” formula.

I guess the UGen is trying to do a DR-II or TDR-II formula, to save one delay register, but the formula must be written incorrectly.

    double y0 = in + b1 * y1;
    ZOUT0(0) = a0 * y0 + a1 * y1;
    y1 = y0;

a0 * y0 + (a1 * y1) expands to (a0 * in) + (a0 * b1 * y1) + (a1 * y1) which reduces to (a0 * in) + ((a0 * b1 + a1) * y1) – now here, I’ll admit that my math isn’t good enough – but in the specific running-average case, first a0 = 1, first b1 = 0. So, first cycle, y0 = in(0), output = y0 * 1 + 0 * 0 = in(0), and second cycle starts with y1 = in(0). Then new a0 = 0.5, new b1 = 0.5, new y0 = in(1) + 0.5 * in(0), new out = 0.5 * (in + 0.5 * in(0)) so the feedback coefficient here is really 0.25, NOT 0.5 as supplied.

This can’t be right.

So I’d like to suggest that the purpose of this thread be to identify what is wrong with the formula that is causing this non-equivalent behavior. General filter safety is a different problem.

hjh

17

Some brainstorming ideas… (sketchy) just workarounds? rs

// In FOS_Ctor
unit->m_a0 = ZIN0(1) == 0 ? 1 : ZIN0(1); // a0 is NEVER zero



if (std::isnan(next_a0) || std::isnan(next_a1) || std::isnan(next_b1) || std::isinf(next_a0) || std::isinf(next_a1) || std::isinf(next_b1)) {
    next_a0 = 1; next_a1 = 0; next_b1 = 0;
}

EDIT
The irony of the SC3 ugen source structure is that it has a special place for initialization, and then it does not do the job to initialize the filters properly )))

18

How does this fix the wrong calculation?

hjh

19

Nothing, just aiming at safe initialization. ex:

unit->m_a0 = ZIN0(1);
unit->m_a1 = ZIN0(2);
unit->m_b1 = ZIN0(3);

// to

unit->m_a0 = std::isnan(ZIN0(1)) || std::isinf(ZIN0(1)) ? 1 : ZIN0(1);
unit->m_a1 = std::isnan(ZIN0(2)) || std::isinf(ZIN0(2)) ? 0 : ZIN0(2);
unit->m_b1 = std::isnan(ZIN0(3)) || std::isinf(ZIN0(3)) ? 0 : ZIN0(3);

The calculation is an equation. If it’s not what the helpfile says… I don’t know what to do. “Correct it”?

20

To create a running average that accumulates indefinitely , you need to keep track of the running sum of the signal and divide it by the number of samples processed so far.

Is that what you said? (maybe I just don’t know the concept)

{   var sig, runningAvg, trig, sr, count;
    sig = SinOsc.ar(440, 0, 0.1);  
    sr = SampleRate.ir;
    count = Sweep.ar(Impulse.ar(0), 1);  
        runningAvg = Integrator.ar(sig) / count;
   [sig, runningAvg]
}.plot(0.3)

image
image758×468 13.2 KB

21

Yes, that’s what a running average is.

(Sigma [i = 1 … n] x(i)) / n

But this is equivalent to the recursive formula. Consider the case of n = 2, average = (a + b) / 2. Then the average for 3 samples = ((a+b)/2) * 2/3 + (c/3) = (a+b)/3 + (c/3) = (a+b+c)/3 and inductively it can be seen that this will extend indefinitely up to any n.

The potential advantage is that after a large number of values, the integrator may reach such a high magnitude that floating point addition may no longer produce a correct result – but the feedback term will remain in the same range (though dividing the new value by n might not be much better, I guess).

FOS’s documentation says you can render this formula using the coefficients that I tried (and my Routine test confirms this). But the behavior is discrepant. I think we should fix the discrepancy, but first I’d like to understand how the formula in FOS’s C++ source is related to the “equivalent” formula in the docs. Or perhaps this is a question for dsp stackexchange…?

hjh

22

Those edge cases make me think it would be useful to instantiate Graphs with another Float type, such as double, or long double (pushing the limit of the compiler), or even higher sampling rates in other cases. Not just a curiosity.

It might be useful to do a Faust version and compile with other flags and types

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.