Efficient sine oscillator

[xoxos]

i'd love to chat about this sometime when i'm not coursing thru kvr intent on accomplishing something unrelated, and when the subject is fresh in my mind so i don't embarass myself.

really though, iirc, my best results generating arbitrary saws with iFFT use a sinc function to set the magnitude of each bin when rendering an arbitrary frequency. i don't want to think about what mystran wrote right now, but i had a terribly difficult time for months trying to evince any method that came anywhere near as close as that or that anyone talked about anywhere on the internet. even though it defies what one "understands" about bins.

i know there's "a" reference out there somewhere which iirc has some information on fft pitch shifting which was somehow impermeable to me. what i didn't do, but would be next, if i had continued to apply myself to the task, would have been using the sinc function to weight the phase of each frequency i am rendering to the bin (remembering that i am accomplishing this with like a 128 bin sinc function, so i would be averaging bin results for hundreds of frequencies per bin).

you can see why i stopped (iirc each bin's phase was set for the closest partial). i should dig up those pleading threads about fft pitch shifting, because this is the same matter as arbitrary iFFT synthesis. i am sure that translating phase to a given bin was trivial considering my crude yet effective grasp of algebra, but haha i could just go on talking and i'm sure nobody believes me or is in the slightest bit concerned as to how i arrived at such extents.

or.. maybe that's all it takes.. is just averaging out all those sinc phase values.. and i'll be there..

[andy-cytomic]

Code: Select all

init:
g = tan(pi*cutoff/samplerate);
g0 = 2/(1 + g*g);
g1 = g*g0;
loop:
const v1 = g0*ic1eq - g1*ic2eq;
const v2 = g1*ic1eq + g0*ic2eq;
ic1eq = v1 - ic1eq;
ic2eq = v2 - ic2eq;

Actually, I've just gone through the math and there is a more efficient way, you can roll the final addition into the coefficients:

Code: Select all

init:
const g = tan(pi*cutoff/samplerate);
const gg = 2/(1 + g*g)
g0 = gg-1
g1 = g*gg
c = cos(2*pi*startphase)
s = sin(2*pi*startphase)
loop:
const t0 = g0*c - g1*s
const t1 = g1*c + g0*s
c = t0
s = t1

I had a hunch that mystrans version where the coefficients start from zero and rise up with frequency would have slightly better numerical performance since you end up have state += g*val type increments, but I've tested out both methods and they are pretty much identical, the error after a week of running is:

0.0001 hz:
error svf mystran=2.17963e-07
error svf andy=2.17963e-07

1 hz:
error svf mystran=3.26555e-07
error svf andy=3.26547e-07

So there is isn't really any practical difference between the two as far as I can tell.

[andy-cytomic]

I may be wrong, but if you're going to use continuous multiplies to create rotations like this, it's probably best to normalize your results once in a while, even when using doubles. The slight errors are compounded through millions of iterations.
It reminds me of a story where someone left a hardware reverb on for days on end, and it started to make annoying humming sounds that slowly got louder and louder - a floating point rounding error where an oscillator got larger and larger over time perhaps?

Doubles are pretty precise if you use the right algorithm. Here are the results after running the oscillator for a week (26671680000 samples) at different frequencies, the "tr" version is the one both mystran and I posted, the "fe" one is a standard forward euler svf like itoa posted:
(updated, I was using the wrong metric before, I now calculate the norm, divide by 1, and subtract 1, and take the abs)

freq=200 samplerate=44,100
error svf tr=4.8e-7
error svf fe=0.0045

freq=2,000 samplerate=44,100
error svf tr=2.8e-7
error svf fe=0.058

freq=20,000 samplerate=44,100
error svf tr=4e-9
error svf fe=8.66

I also tried some square wave modulation between 200hz and 20000hz and the fe would blow up sometimes, or have an error of around 8 at others, the tr stays at an error of around 2.5e-7 or lower no matter what you do.

(ps: so that means that for the trapezoidal version after a week of modulating it has hard as you want the amplitude is around 2e-6 of a dB away from 1.0)

[mystran]

Also.. even if you do normalize once a block or something, the less accurate formulas tend to accumulate a lot more phase-drift; even where it's not fatal, it can be quite annoying if you're trying to do something like additive synthesis.

[DaveHoskins]

I may be wrong, but if you're going to use continuous multiplies to create rotations like this, it's probably best to normalize your results once in a while, even when using doubles. The slight errors are compounded through millions of iterations.
It reminds me of a story where someone left a hardware reverb on for days on end, and it started to make annoying humming sounds that slowly got louder and louder - a floating point rounding error where an oscillator got larger and larger over time perhaps?

Doubles are pretty precise if you use the right algorithm. Here are the results after running the oscillator for a week (26671680000 samples) at different frequencies, the "tr" version is the one both mystran and I posted, the "fe" one is a standard forward euler svf like itoa posted:
(updated, I was using the wrong metric before, I now calculate the norm, divide by 1, and subtract 1, and take the abs)

freq=200 samplerate=44,100
error svf tr=4.8e-7
error svf fe=0.0045

freq=2,000 samplerate=44,100
error svf tr=2.8e-7
error svf fe=0.058

freq=20,000 samplerate=44,100
error svf tr=4e-9
error svf fe=8.66

I also tried some square wave modulation between 200hz and 20000hz and the fe would blow up sometimes, or have an error of around 8 at others, the tr stays at an error of around 2.5e-7 or lower no matter what you do.

(ps: so that means that for the trapezoidal version after a week of modulating it has hard as you want the amplitude is around 2e-6 of a dB away from 1.0)

Thanks for doing an actual test, I guess there are other factors like the complexity of each loop, and multiplies and divides involved, and how small or large the rotation is perhaps?
I had a rotation using floats that only lasted a minute at certain frequencies, using floats.

Anyway, what happens at week 3! Is there a sudden exponential explosion?
I'm still going to normalize occasionally, I think.

[andy-cytomic]

Anyway, what happens at week 3! Is there a sudden exponential explosion?
I'm still going to normalize occasionally, I think.

The error looks to be roughly linear with time:

Code: Select all

algorithm:
            const double t0 = q0*c0 - q1*s0;
            const double t1 = q1*c0 + q0*s0;
            c0 = t0 - c0;
            s0 = t1 - s0;

samplerate = 44100

20000 hz:
min=1 err=4.84834e-12
min=2 err=9.73088e-12
min=4 err=1.95433e-11
min=8 err=3.91331e-11
min=16 err=7.80123e-11
min=32 err=1.56644e-10
min=64 err=3.12756e-10
min=128 err=6.26296e-10
min=256 err=1.25338e-09
min=512 err=2.50805e-09
min=1024 err=5.01869e-09
min=2048 err=1.00376e-08
min=4096 err=2.00766e-08

2 hz:
min=1 err=1.93492e-10
min=2 err=3.8556e-10
min=4 err=7.74586e-10
min=8 err=1.55843e-09
min=16 err=3.12185e-09
min=32 err=6.25071e-09
min=64 err=1.25084e-08
min=128 err=2.50072e-08
min=256 err=5.00369e-08
min=512 err=1.00086e-07
min=1024 err=2.00169e-07
min=2048 err=4.00321e-07
min=4096 err=8.00701e-07

So taking the 2hz case as the largest error it would take 6e7 minutes to get to 0.1 dB away from 1, which is around 100 years. So yes, re-normalising every month or so would be a good idea if you plan on running your oscillator for years at a time.

[Nalodx]

Looks very good. So basically its a 2 pole state variable filter where the state values correspond
to real and imaginary values. Output is directly feeded back to the input using trapezoidal integration or zero delay feedback. This is a nice backdoor to complex math (multiple angle formulars, etc).
Doing it this way you are bypassing the problem like the chamberlin based one which makes problems at around nyquist/4 or make even more problems if you try to modulate fast (fm).
I am always excited if i hear filter makes problems up to nyquist. If feedback gets solved in a correct way this vanishes.
2 pole state variable filter always offers some surprises.

Mystran and your code seems to be more compact than mine.
Do you face the same issue that if you do it with floats (sse2) the error adds up to quickly. With FPU/or doubles(SSE2) is works like a charm.
(Maybe i should buy a faster machine!!! SS3 then 4 doubles could be handled).

[andy-cytomic]

Looks very good. So basically its a 2 pole state variable filter where the state values correspond
to real and imaginary values. Output is directly feeded back to the input using trapezoidal integration or zero delay feedback. This is a nice backdoor to complex math (multiple angle formulars, etc).
Doing it this way you are bypassing the problem like the chamberlin based one which makes problems at around nyquist/4 or make even more problems if you try to modulate fast (fm).
I am always excited if i hear filter makes problems up to nyquist. If feedback gets solved in a correct way this vanishes.
2 pole state variable filter always offers some surprises.

Mystran and your code seems to be more compact than mine.
Do you face the same issue that if you do it with floats (sse2) the error adds up to quickly. With FPU/or doubles(SSE2) is works like a charm.
(Maybe i should buy a faster machine!!! SS3 then 4 doubles could be handled).

The state variable filter in its continuous form is a direct state space representation of a biquad, the values at each capacitor correspond to the state, so this is directly a sin / cos generator by design.

There is no zero delay here, I think you are getting confused. The phase delay of the filter is the cutoff frequency, delay is required to make it work. The equations result from using trapezoidal integration, which is an implicit numerical integration method, so it that sense it is not an explicit method (like the chamberlin svf) if that is what you mean by "delay", but such making such distinctions is not useful in general.

For correct operation at low frequencies you will need double precision numbers to have enough precision so that the accumulation of state can still be resolved. This is easier to see in mystrans += form, so if the += is operating on a large state and trying to increment it by a very small number problems will arrise:

eg:
float v = 1.0;
v += 1e-20f; // this won't work!

[mystran]

I'd like to point out that while Andy pointed out the similarity with my formula (well, not really "my formula" in the sense that I'd come up with it; it's a standard well-known recurrence, found in many FFTs among other things) and TR, doing any kind of fancy integration isn't really necessary at all, you can arrive at the same results just optimizing numerical performance of standard incremental cos/sin rotation.

Basically, all you need to do is to observe that the cosine term is losing a lot of precision, because small angles result in cosines very close to 1, which biases the coeffs with a fixed large exponent and results in rounding. So then you consult your tables of trigonometric identities, and find cos(2*q) = 1-2*sin(q)^2 (as it's usually listed, but one can then substitute p = q/2 to get what we want), which splits the coefficient into a difference that has a fixed unit value (so multiply free and no rounding) and a variable term that approaches zero for small angles (so gets stored with small exponent for good precision).

So from (x+i*y)*(cos(p)+i*sin(p)) we can rewrite it as:
(x+i*y)*((1 - 2*sin(p/2)^2) + i*sin(p)) = (x+i*y)+(x+i*y)*(-2*sin(p/2)^2 + i*sin(p))

There's no fancy integration or filtering required here. It's simply the usual complex rotation formula, but with a substitution to keep numbers in a better range and reorder the calculations in order to better preserve precision. It's also worth noting that when implemented, if you actually try to optimize (or your compiler does that for you; fortunately they usually won't, even with fast-math) by calculating (1-2*sin(p/2)^2) first (to save the final additition) then you lose all the advantages and get numerical performance identical to just using cos/sin.

[Nalodx]

thanks for info andy. I think i have just seen your divide by (1.0f+g*g) which
is what stays when solving the feedback of 2 cascaded DF2 Integrators where each
feed back to input. (Sample Input related code of SVF removed). However, same add problems with float SSE2), double,float(FPU) works of course.
I got this after looking at Neotecs stuff. Its okay if you want to modulate freq very fast (6 Khz). using only one coefficient.
If freq is constant i would prefer the standard rotation.

init: buf1 = 1.0
buf2 = 0;

g = tan(pi*cutoff/samplerate); //for low freqs no tan needed

float low = (buf2 + g * (buf1))/ (1.0f+g*g);
float band = (buf1-g*low);

buf1 = band - g*low; //buf1 re/cos
buf2 = low + g*band; //buf2 im/sin

[andy-cytomic]

..., if you actually try to optimize (or your compiler does that for you; fortunately they usually won't, even with fast-math) by calculating (1-2*sin(p/2)^2) first (to save the final additition) then you lose all the advantages and get numerical performance identical to just using cos/sin.

Thanks for pointing that out Mystran. When I tested your version I initialised the formula using double precision tan and a division, which made the coefficients less accurate so my other version looked like it was doing better. I used the sin version and indeed the numerical performance of the += version was best (which is what I expected!)

Just remember that if you want to use a sin approximation to initialise this algorithm you need them to be exact otherwise you will kill how accurate the results are. So you have several tradeoffs:

[*] making two calls to math sin, most accurate coefficients
[*] making a single call to tan / tan approximation and using a division, not as accurate coefficients, but both coefficients will match each other well even if your tan wasn't so accurate (you will just have a less accurate frequency of the oscillator
[*] using either the += version so taking an extra addition but being more accurate
[*] using the version without the += which is fewer operations but less accurate (but still runs for weeks without issue)

[mystran]

Here's another thought.. if one has two coefficient pair as above (with cos-1 and sine), then one can also add the angles directly (which could be quite useful for interpolation, etc) without calculating cosines.

Basically, let a0,b0 be the first angle (well, the coefficient pair), a1,b1 be the second one, then we can calculate:
(1+a0+i*b0)*(1+a1+i*b1) - 1 =
(1 + a0 + i*b0 + a1*a0 + a1 + i*a1*b0 + i*b1 + i*b1*a0 - b0*b1) - 1 =
(a0 + a1 + a0*a1 - b0*b1) + i*(b0 + b1 + a1*b0 + a0*b1)

The last form has 4 extra additions over regular complex multiply, but saves one from having to add 1 (and potentially kill accuracy for small angles/slow frequencies) to any of the coefficients.

Alternatively in complex numbers I believe it becomes:

c = c0 + c1 + c0*c1

I've not actually tested that exact variant numerically, but I use bunch of very similar formulations, so should work.

[xoxos]

you guys..

needed to do one of these yesterday. used the two multiply, no-transcendentals-in-coefficients form from the first post. to inhibit explosion under duress i hardclipped the output state.

unlike andy, i have no idea if that will sustain under all sample rates and modulations thrown at it, but i can complete a thousand shitty implementations a week (i'll just stash these notes for later..)

[andy-cytomic]

Slightly OT: but I've implemented a "sin only" coefficient calculation like mystran has done for the sin osc, but this is for the full linear SVF. I expect that like the sin osc this form with the best numerical properties (especially with low cutoffs), and also the most accurate coefficients (eg the cutoff is where you want it):

www.cytomic.com/files/dsp/SvfLinearTrapezoidalSin.pdf

[mystran]

Slightly OT: but I've implemented a "sin only" coefficient calculation like mystran has done for the sin osc, but this is for the full linear SVF. I expect that like the sin osc this form with the best numerical properties (especially with low cutoffs), and also the most accurate coefficients (eg the cutoff is where you want it):

I'd like to object against your idea that converting from DF coeffs to SVF is necessarily a bad idea.. because in some cases you might not want to use BLT at all and it might be the case that you simply end up with a bunch of "abstract" digital poles/zeroes, yet you'd still prefer the SVF as an implementation structure. Obviously converting cookbook filters is silly (you could just as well prewarp the cutoff and take rest of the coeffs from the analog proto directly) but for other purposes it can be quite useful.

edit: I do admit though, that if you're actually solving/fitting/whatever the digital poles/zeroes directly, it might often make sense to deal with the poles (at least) in polar form, and in that case converting it into a DF polynomial first would be a bit pointless...