[sc-users] Re: Sources of NaN in a UGen?

["Dan Stowell" <danstowell@xxxxxxxxx>]

OK -

I think I've spotted the problem. I forgot about (inf - inf) being
NaN, so that's probably what's happening.

Thanks anyway
Dan



2006/5/14, Dan Stowell <danstowell@xxxxxxxxx>:
> Hi -
>
> Writing a UGen to represent the "Chua's Circuit" chaotic oscillator, I
> find that the UGen tends very quickly to outputting NaN, although I
> can't tell why. The meat of the UGen only does multiplication and
> addition (there are a couple of trivial divisions). Can anyone suggest
> why this might be happening?
>
> The UGen code is:
>
> void ChuaL_next(ChuaL *unit, int inNumSamples)
> {
>         float *out = ZOUT(0);
>         float freq = ZIN0(0);
>         double a = ZIN0(1);
>         double b = ZIN0(2);
>         double c = ZIN0(3);
>         double d = ZIN0(4);
>         double rr = ZIN0(5); // Reciprocal of r - to save on a
> division operation
>         double h = ZIN0(6);
>         double x0 = ZIN0(7);
>         double y0 = ZIN0(8);
>         double z0 = ZIN0(9);
>
>         double xn = unit->xn;
>         double yn = unit->yn;
>         double zn = unit->zn;
>         float counter = unit->counter;
>         double xnm1 = unit->xnm1;
>         double ynm1 = unit->ynm1;
>         double znm1 = unit->znm1;
>         double frac = unit->frac;
>
>         float samplesPerCycle;
>         double slope;
>         if(freq < unit->mRate->mSampleRate){
>                 samplesPerCycle = unit->mRate->mSampleRate /
> sc_max(freq, 0.001f);
>                 slope = 1.f / samplesPerCycle;
>         }
>         else {
>                 samplesPerCycle = 1.f;
>                 slope = 1.f;
>         }
>
>         if((unit->x0 != x0) || (unit->y0 != y0) || (unit->z0 != z0)){
>                 xnm1 = xn;
>                 ynm1 = yn;
>                 znm1 = zn;
>                 unit->x0 = xn = x0;
>                 unit->y0 = yn = y0;
>                 unit->z0 = zn = z0;
>         }
>
>         double dx = xn - xnm1;
>
>         for (int i=0; i<inNumSamples; ++i) {
>                 if(counter >= samplesPerCycle){
>                         counter -= samplesPerCycle;
>                         frac = 0.f;
>
>                         xnm1 = xn;
>                         ynm1 = yn;
>                         znm1 = zn;
>
>                         double k1x, k2x, k3x, k4x,
>                                 k1y, k2y, k3y, k4y,
>                                 k1z, k2z, k3z, k4z,
>                                 kxHalf, kyHalf, kzHalf,
>                                 ymx, ymx3, ymx3b, hOverR, xnm1plus;
>
>                         hOverR = h * rr;
>
>                         // 4th order Runge-Kutta approximate solution
> for differential equations
>
>                         ymx = ynm1 - xnm1;
>                         ymx3 = ymx * ymx * ymx;
>                         ymx3b = ymx3 * b;
>                         k1x = hOverR * (d * ymx - c * xnm1 - (a * xnm1
> * xnm1 * xnm1) - ymx3b);
>                         k1y = - (h * (d * ymx - znm1 - ymx3b));
>                         k1z = h * ynm1;
>                         kxHalf = k1x * 0.5;
>                         kyHalf = k1y * 0.5;
>                         kzHalf = k1z * 0.5;
>
>                         ymx = ynm1 + kyHalf - xnm1 - kxHalf;
>                         ymx3 = ymx * ymx * ymx;
>                         ymx3b = ymx3 * b;
>                         xnm1plus = xnm1 + kxHalf;
>                         k2x = hOverR * (d * ymx - c * (xnm1 + kxHalf)
> - (a * xnm1plus *
> xnm1plus * xnm1plus) - ymx3b);
>                         k2y = - (h * (d * ymx - znm1 - kzHalf - ymx3b));
>                         k2z = h * (ynm1 + kyHalf);
>                         kxHalf = k2x * 0.5;
>                         kyHalf = k2y * 0.5;
>                         kzHalf = k2z * 0.5;
>
>                         ymx = ynm1 + kyHalf - xnm1 - kxHalf;
>                         ymx3 = ymx * ymx * ymx;
>                         ymx3b = ymx3 * b;
>                         xnm1plus = xnm1 + kxHalf;
>                         k3x = hOverR * (d * ymx - c * (xnm1 + kxHalf)
> - (a * xnm1plus *
> xnm1plus * xnm1plus) - ymx3b);
>                         k3y = - (h * (d * ymx - znm1 - kzHalf - ymx3b));
>                         k3z = h * (ynm1 + kyHalf);
>
>                         xnm1plus = xnm1 + k3x;
>                         k4x = hOverR * (d * ymx - c * (xnm1 + k3x) -
> (a * xnm1plus *
> xnm1plus * xnm1plus) - ymx3b);
>                         k4y = - (h * (d * ymx - znm1 - k3z - ymx3b));
>                         k4z = h * (ynm1 + k3y);
>
>                         xn = xn + (k1x + 2.0*(k2x + k3x) + k4x) * ONESIXTH;
>                         yn = yn + (k1y + 2.0*(k2y + k3y) + k4y) * ONESIXTH;
>                         zn = zn + (k1z + 2.0*(k2z + k3z) + k4z) * ONESIXTH;
>
>                         dx = xn - xnm1;
>                 }
>                 counter++;
>                 ZXP(out) = (xnm1 + dx * frac) * 0.5f;
>                 frac += slope;
>         }
>
>         unit->xn = xn;
>         unit->yn = yn;
>         unit->zn = zn;
>         unit->counter = counter;
>         unit->xnm1 = xnm1;
>         unit->ynm1 = ynm1;
>         unit->znm1 = znm1;
>         unit->frac = frac;
> }
>
> void ChuaL_Ctor(ChuaL* unit){
>         SETCALC(ChuaL_next);
>
>         unit->x0 = unit->xn = unit->xnm1 = ZIN0(7);
>         unit->y0 = unit->yn = unit->ynm1 = ZIN0(8);
>         unit->z0 = unit->zn = unit->znm1 = ZIN0(9);
>         unit->counter = 0.f;
>         unit->frac = 0.f;
>
>         ChuaL_next(unit, 1);
> }
>
>
> and the SC class:
>
> ChuaL : ChaosGen {
>         *ar { arg freq=22050, a=0.3286, b=0.9336, c= -0.8126, d=0.399,
> rr=0.6.reciprocal, h=0.05, xi=0.1, yi=0, zi=0, mul=1.0, add=0.0;
>                 ^this.multiNew('audio', freq, a, b, c, d, rr, h, xi,
> yi, zi).madd(mul, add)
>         }
> }
>
> Thanks in advance for any advice -
> Dan
>
> P.S. Not sure if this is an "sc-users" or "sc-dev" issue so sorry for
> sending to both lists, for those who read both.
>
> --
> http://www.mcld.co.uk


--
http://www.mcld.co.uk