How to make this infinite decay envelope (y=1/(x+1)^c curve) more CPU efficient?

[mikejm]

I've built a basic ADSR envelope adapted from both of:
http://www.martin-finke.de/blog/article ... envelopes/
http://www.earlevel.com/main/2013/06/03 ... adsr-code/

My design requirement was for the decay to be infinite (based on y= 1/(x+1)^c) where "decay time" is time to an amplitude of 1/e relative to the peak value. That's how envelope time works in Reaktor and I've grown accustomed to it.

The envelope works perfectly well, but I'm doing additive synthesis where I have hundreds of these things running together, and I am finding them quite costly in terms of CPU. I can reduce the CPU with a sample rate reduction for the envelope, so it's only outputting a calculation every x number of samples, but this can introduce clicks on the attack phase. I'd rather figure out how to make it more internally efficient if possible.

Here's the graphing function if you are curious what I mean by y=1/(x+1)^c:
https://www.desmos.com/calculator/sby91g32yu

My code is that when I enter the Decay case:

Code: Select all

        case ENVELOPE_STAGE_DECAY:
                        currentSampleIndex = 0;			
                        sustainLevel = sustainPercent * maxLevel;
			decayExponent = 1 / log(decayTimeSeconds+1);
			break; 

"decayExponent" is meant to represent the "c" in 1/(x+1)^c, so it's only calculated once at the case initiation to save calculations per output...

Then my nextSample() function for decay is this:

Code: Select all

	else if (currentStage == ENVELOPE_STAGE_DECAY) {
		
		currentLevel = maxLevel / pow((currentSampleIndex/sampleRate)+1.0, decayExponent);
		currentSampleIndex++;
		
		if (currentLevel <= (sustainLevel + 0.000001)) {
			currentLevel = sustainLevel;
			enterStage(ENVELOPE_STAGE_SUSTAIN);

So the "x" in my 1/x^c equation is given by (currentSampleIndex/sampleRate) + 1.0, because I want x=1 to start at 1. I can save one division by instead multiplying currentSampleIndex by "deltaTime" that is precalculated as 1/sampleRate. But I still don't think that's as optimal as it can get.

Is there another cheaper way of calculating these types of infinite y=1/x^c decay curves? eg. Using recursive math?

I have no formal training in math or coding or DSP so I'm just trying to figure it out as I go and maybe I'm doing things in a really dumb way. I would very much appreciate any pointers or clarification if possible.

Thanks

[syntonica]

You've got three speed issues: the Divide, the Pow and the Compare. They are all dreadful speedwise when you apply them to each and every sample.

Idea one: There's nothing special about that curve that can't be modeled by a faster function, one that doesn't involve dividing. Find something that tends to zero, like (1-tanh(x)). That's of course, maybe not the best function, just a for instance.

Idea two: I essentially use a lookup table of 1024 values created with a cubic (or a couple other functions for different decay curves) that is scaled by the release time. This is enough values that any stepping is not audible and by note end, it's smoothly landed at zero. You can pre-calculate a LUT of any size you prefer. Carefully implemented, you can eliminate two, if not all three of the above.

(As an aside, never, ever recurse if you can help it. They are the slowest of functions. They are the sloths of functions. They will eat all your pizza and drink all your beer and not leave. They... well, you get the idea.)

[chi cubed]

maybe currentLevel = maxLevel * pow((currentSampleIndex/sampleRate)+1.0, -decayExponent); with the division replaced by a subtraction before the exp

[mikejm]

maybe currentLevel = maxLevel * pow((currentSampleIndex/sampleRate)+1.0, -decayExponent); with the division replaced by a subtraction before the exp

Hey good idea. So I started with:

currentLevel = maxLevel / pow((currentSampleIndex/sampleRate)+1.0, decayExponent);

I replaced the 1/sampleRate with "deltaTime":

currentLevel = maxLevel / pow((currentSampleIndex * deltaTime)+1.0, decayExponent);

Finally if I do your suggestion I get:

currentLevel = maxLevel * pow((currentSampleIndex * deltaTime)+1.0, negativeDecayExponent);

Took out two divisions per output. Not bad.

[Miles1981]

You can use vectorized pow calls, and if you have an integer exponent, even better. The power function will do the division inside, so indeed better to use a multiplication.

Have a look at svml and fmath.

[Z1202]

All exponential (and linear) envelope curves can be generated recursively by
x[n] = a*x[n-1]+b
This is also how Reaktor Core envelopes work.
In comparison, pow function is very expensive. The problem with the recursion is that at some settings (esp. with a>1) there can be quickly accumulating precision losses, which can be mitigated to an extent by choosing a "better" origin for x.

There has been recently a discussion about that on music-dsp mailing list, I don't have a link to the thread, but it shouldn't be difficult to find.

[Max M.]

The problem with the recursion is that at some settings (esp. with a>1) there can be quickly accumulating precision losses

And there comes the usual trick of using an integrator in the feedback:
k = 1 - a;
x[n] = x[n-1] - k*x[n-1] + b;
Doesn't it? (Though I assume this is already suggested at the mentioned mailing list thread).

[Z1202]

The problem with the recursion is that at some settings (esp. with a>1) there can be quickly accumulating precision losses

And there comes the usual trick of using an integrator in the feedback:
k = 1 - a;
x[n] = x[n-1] - k*x[n-1] + b;
Doesn't it? (Though I assume this is already suggested at the mentioned mailing list thread).

Hmmm, I'm not sure how the integrator helps with a>1. At any rate with the natural order of operations in your formula we first compute x[n-1] - k*x[n-1], whereas IIRC the major source of precision losses is usually the addition of b (therefore the trick is to choose the origin of x such that b is small, where the sign also can matter, or even zero).

[2DaT]

Replace pow(x,c) with exp(log(x)*c). Use single precision, because single precision transcedentals are cheaper.

[Max M.]

Z1202 At any rate with the natural order of operations in your formula we first compute x[n-1] - k*x[n-1] ...

Well, probably I'd better to initially write it in proper statements instead of abstract oneliner:
x[n] = x[n-1]; // not loosing anything here
x[n] -= x[n]*k; // maximum precision when a->1 (i.e. longer decay) and x->0

It's just exactly how we get better numerical stability and SNR for floating point filter implementations. So I can't see any significant difference in this case (unless you're actually assume weird curves like `b->1 and a->1.5...2`).

[Z1202]

Z1202 At any rate with the natural order of operations in your formula we first compute x[n-1] - k*x[n-1] ...

Well, probably I'd better to initially write it in proper statements instead of abstract oneliner:
x[n] = x[n-1]; // not loosing anything here
x[n] -= x[n]*k; // maximum precision when a->1 (i.e. longer decay) and x->0

It's just exactly how we get better numerical stability and SNR for floating point filter implementations. So I can't see any significant difference in this case (unless you're actually assume weird curves like `b->1 and a->1.5...2`).

IIRC for a>>1 you tend to have a difference of two large values of comparable magnitude when adding b. That's what I meant saying that b is the major source of precision loss. In comparison computing x[n-1] + x[n-1]*(a-1) vs. a*x[n-1] probably doesn't make that much difference.

[mikejm]

Replace pow(x,c) with exp(log(x)*c). Use single precision, because single precision transcedentals are cheaper.

Okay so I'll switch it to floats. Also if that's cheaper, I can switch that too. So I get:

currentLevel = maxLevel * exp(log((currentSampleIndex * deltaTime)+1.0)*negativeDecayExponent);

I can also instead of multiplying and adding (currentSampleIndex * deltaTime +1) each time just have a deltaTimeIndex initialized at the start of the case to 1.0 and write it like this:

currentLevel = maxLevel * exp(log(deltaTimeIndex) * negativeDecayExponent);
deltaTimeIndex = deltaTimeIndex + deltaTime;

Did I get my brackets/concepts there correct? Now I've reduced it to two multiplications and two natural transcendentals...

All exponential (and linear) envelope curves can be generated recursively by
x[n] = a*x[n-1]+b
This is also how Reaktor Core envelopes work.
In comparison, pow function is very expensive. The problem with the recursion is that at some settings (esp. with a>1) there can be quickly accumulating precision losses, which can be mitigated to an extent by choosing a "better" origin for x.

There has been recently a discussion about that on music-dsp mailing list, I don't have a link to the thread, but it shouldn't be difficult to find.

The problem with the recursion is that at some settings (esp. with a>1) there can be quickly accumulating precision losses

And there comes the usual trick of using an integrator in the feedback:
k = 1 - a;
x[n] = x[n-1] - k*x[n-1] + b;
Doesn't it? (Though I assume this is already suggested at the mentioned mailing list thread).

So let's say I have an equation of

y= 1/(x+1)^c, and I have calculated c. I have sampleRate, deltaTime variables.

How would I convert that to a:
k = 1 - a;
x[n] = x[n-1] - k*x[n-1] + b;

Type format? What is the math to calculate k, a, b in this case? How do I graph it on desmos.com to confirm?

[sonigen]

since..

1/a^b = a^-b

you can simplify your eq to

y = (x+1)^-c

but that's all you can do i think. You cant get rid of that pow function and you cant convert it to what vadim is talking about since one is y=x^c and one is y=c^x, the later can be done recursively. Tbh if you want low cpu you need to forget about the equation you're using. It also has the problem that if you change 'c' midway you will get a get a sudden jump in the envelope, you don't get that problem with a recursive form, if you change the coefs it just changes the slope from the current position.

[mikejm]

since..

1/a^b = a^-b

you can simplify your eq to

y = (x+1)^-c

but that's all you can do i think. You cant get rid of that pow function and you cant convert it to what vadim is talking about since one is y=x^c and one is y=c^x, the later can be done recursively. Tbh if you want low cpu you need to forget about the equation you're using. It also has the problem that if you change 'c' midway you will get a get a sudden jump in the envelope, you don't get that problem with a recursive form, if you change the coefs it just changes the slope from the current position.

Yes exactly. That's what I thought about x^c vs c^x and recursion. We did work out a few other optimizations I listed in my post above though.

[Max M.]

Z1202
That's what I meant saying that b is the major source of precision loss.

Now I see what you mean. And I agree.

mikejm
What is the math to calculate k, a, b in this case? How do I graph it on desmos.com to confirm?

Roughly:
For `x[n] = a*x[n-1]` the impulse responce is `y = 1*a^n` (where n is number of iterations starting from 0)
For `x[n] = a*x[n-1] - b` it becomes `y = (1-b*n)*a^n`
and `n` -> `T in seconds` is obviously `n/sampleRate`.
E.g.:
https://www.desmos.com/calculator/cdjre0xvfq
https://www.desmos.com/calculator/2emxae2d6m

The rest depends on how exactly you're going to use it. Specifically if you need a separate "curve shape" parameter. In simple case we can for example choose values of the two curve points:
end point - to represent decay time,
some mid point - for its x or y value to specify the "curveness".
Then it becomes a usual "two-equations -> two-unknowns" solution.