#1698 Self-oscillating 8580 filter resonance

[encore64]

While finishing the Sonic the Hedgehog-music I noticed in x64sc/vsid that when I used filter resonance, on 8580 it would self-oscillate and sustain after the tone was released. I found a work-around by turning off the resonance at the end of a note, but obviously that is not needed on a real 8580.

Here's a test-file, it plays 3 tones:
1. A pulsetone with lowpass filter with resonance set to F and the frequency sweeping down and up repeatedly. When gate is set to off, the resonance can be heard self-oscillating as a tone.
2. A pulsetone with lowpass filter with resonance set to 0 and the frequency sweeping down and up repeatedly. When gate is set to off, the resonance does not self-oscillate.
3. A pulsetone with lowpass filter with resonance set to F, frequency sweeping down and then the resonance sweeping down and up repeatedly. When gate is set to off, the self-oscillation can be heard until the resonance is sweeping down towards 0.

I've attached a .sid-file as well as a sample from 8580 on my 64 reloaded, as well as from vsid 3,6,2-dev_rev_41765. I also edited and made an exaggerated version so it's easy to hear the issue itself. Hope this helps.

[Leandro Nini]

I guess you're using a low quality sampling method. Can you please check in the Audio/SID settings page?

[encore64]

Hi Leandro.
If you mean the ReSID sampling method-option, I'm using the default "Resampling" (as the filename also says) for the audio-clip, but it's the same behaviour in "Fast", "Interpolation" and "Fast resampling" too.

[Leandro Nini]

Hmm, you're right, it's not the resampling.
It's not happening with residfp which is supposedly the same code... will check what's going wrong.

[Leandro Nini]

I think I've found the cause of the problem, if you remove the + (1 << 11) term from https://sourceforge.net/p/vice-emu/code/HEAD/tree/trunk/vice/src/resid/filter8580new.h#l1535 the ringing disappears.

Not entirely sure about the math behind this, it seems that the +1 should not be there in the function derivative:

f = a*(b - vx)^2 - c - (b - (vx + x))^2
df = 2*((b - (vx + x))*(dvx + 1) - a*(b - vx)*dvx)

[gpz]

I asked our math guru about the +1 and he says it depends on what kind of derivative this is actually ment to be. for "partial derivative" the +1 is wrong (i hope i am relaying this correctly)

[Leandro Nini]

If my math is not too rusty this might make sense. Still have some doubt on why signs are inverted though:

f = a*(b_vx)^2 - c - (b_vo)^2
  = a*(b_vx)^2 - c - (b_vx - x)^2
  = a*(b_vx)^2 - c - (b_vx^2 + x^2 - 2*x*b_vx)


df/dvx = 2*a*b_vx - 2*b_vx + 2*x
       = 2*(a*b_vx - b_vx + x)
       = 2*(a*b_vx - (b_vo - x) + x)
       = 2*(a*b_vx - b_vo)

df = 2*(a*b_vx - b_vo)*dvx
   = 2*(a*b_vx*dvx - b_vo*dvx)

[Leandro Nini]

Here's the patch, it affects both models so watch out for regressions.
It looks good with the tunes I've tested and UBSan doesn't complain.

[gpz]

applied patch in r41943 - please test!

[encore64]

Hi Leandro, great that you're working on it! I've tested r41944 and only 8580 with the .sid that I provided in the first post.

Test 1 - the ringing has improved, now only a very faint ringing can be heard when gate is set to off, but it's static in frequency - it doesn't follow the filter frequency like it did before.

Test 2 - (resonance set to 0) here there was no ringing before - now there is! So this one is worse.

Test 3 - the ringing have improved.

[Bjoern Froehlich]

Hi there, gpz pointed me to this discussion so I dare to ask some details concerning the math in here.

If I understood the whole thing correctly, f is a root function, and calculating df is needed to apply Newton's method.
- I take it that directly resolving f=0 to vo is no option due to the sqrt-term that will occur.
- If we go the path and solve f=0 via Newton's method, I remember it to need the derivative of the variable the root function depends on. Since the terms have been slightly adjusted by substituting vo = vx + x, this should give the dependency f = f(x).
- Thus f(x) = a*(b - vx)^2 - c - (b - (vx + x))^2
yields df/dx = 2*(b - vx - x)
-> What is the reason for the used df being more complex than this?
-> Why are terms for dvxexplicitely calculated in the code? If any, it should be dx, but I don't see a reason to actually calculate them (after all, these are just basis elements of the cotangential space these df-elements live in)
If we need to consider more than just x as variable (say vx and x), then it might make sense to calculate the total differential of f, but this will turn the whole thing into something more-dimensional (and -while still possible-Newton's method will surely look more evil, driving the motivation to just use a sqrt-function into new heights).

Thanks in advance for any enlightments ;)

[Leandro Nini]

Hi Bjoern, thanks for joining the discussion.

f is actually a fuction of vx, and x a "mapping function x = vo - vx -> vx" (from the comment in the code) and we derive with respect to vx. Not sure if we shall consider x as dependent on vx here.

dvx is precalculated so we can get df by multiplying the derivative df/dvx by dvx, if I get it right, then we can obtain the Newton-Raphson term as f/df.

As a side note, in residfp the code is slightly different, vx and vo are independent and the derivative, done with respect to vo, is simpler:

df = 2*((b - vo)*dvo - a*(b - vx))

[Bjoern Froehlich]

Hi Leandro,

thanks for having me ;)

Though I'm sure a lot of the code does some shifting & scaling magic I did not yet understand completely, I do not agree on how the Newton-Raphson-method (if we want to use the full term, hehe;)) is currently applied.

Given an equation we want to resolve to a certain variable, like x^2 = a, instead of acutally resolving it, a root function is introduced. This root function depends on the variable we originally wanted to resolve it to - so the approach shifts the challenge of resolving an equation to finding a zero of the corresponding root function. The example above'd give f = x^2 -a, and for the Newton-Raphson-method it is necessary to handle the root function as f = f(x).

The equation we have to deal with here (that's at least what I understand) is one that brings certain voltage values together. For a given vx, we now need to determine vo. Doing this via Newton-Raphson means: vx is assumed to be given (read: constant), vo is the variable the root function f depends on. By using this substitution vo = vx + x, we make xthe new variable for f - vx will continue to be just a given value!

So, for my understanding, applying Newton-Raphson is just doing
x_{n+1} = x_n - f(x_n) / f'(x_n)
until x_n makes fsmall enough.

Since fand f'can be calculated without any dxas explicit term, I also do stick to my thinking that the whole thing can be done without the need to actually add any differential terms to the calculation, e.g.
f(x_n) = a*(b - v*x_n)^2 - c - (b - (vx + x_n))^2 = gives a value, say f_n

f'(x_n) = 2*(b - vx - x_n) = gives a value, say f'_n

x_{n+1} = x_n - (f_n / f'_n) = new value for x

Ok, a lot of text, somehow the same I already wrote in my first post - but this looks pretty correct for me this way.

[Leandro Nini]

Ok, I understand your point now.
The Newton step is x = x - f(x) / f'(x) where f'(x) = df/dx, while in the code it does x = x - f/df which is not theoretically correct.
However trying to implement the suggested change produces no sound :/

[Leandro Nini]

Just for the record, I've finally been able to make sense of this.

We have the root function of the variable vx:
f(vx) = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - (vx + x))^2
where
vx + x = vo
We can expand the function to:

f(vx) = a*(b - vx)^2 - c - (b - (vx + x))^2
      = a*(b^2 + vx^2 - 2*b*vx) - c - (b^2 + (vx + x)^2 - 2*b*(vx + x))
      = a*b^2 + a*vx^2 - 2*a*b*vx - c - b^2 - (vx + x)^2 + 2*b*(vx + x)
      = a*b^2 + a*vx^2 - 2*a*b*vx - c - b^2 - vx^2 - x^2 - 2*x*vx + 2*b*vx + 2*b*x

Then we calculate and expand the derivative:

f'(vx) = 2*a*vx - 2*a*b  - 2*vx - 2*x + 2*b
       = 2*(a*vx - a*b - vx - x + b)
       = 2*(a*(vx - b) + b - (vx + x))
       = 2*(b - (vx + x) - a*(b - vx))

Then, given f'(x) = df/dx, we have the resulting

df = 2*((b - (vx + x)) - a*(b - vx))*dvx

so the patch is ok even from a theoretical point of view.

[gpz]

So can we close this ticket? Was everything applied and all is fine? I am getting confused again just be scrolling over the posts here :)

Perhaps the above should be added in comments at a strategic place in the code? :)

[Leandro Nini]

Yes, the original problem is fixed, I'll post a patch for the comments in the upcoming days.

[gpz]

cool, please post on the patches tracker - i am closing this ticket

[bozz64]

The patch also fixes the issue with tune: /MUSICIANS/P/Psych858o/90_Seconds.sid

[bozz64]

However, the patch introduces a new problem with playing samples on a 8580 when no voices are played, see issue:

https://sourceforge.net/p/vice-emu/bugs/1696/#4523

[Leandro Nini]

Hmm, it looks more like a rounding issue.
@encore64 @bozz64 can you please try the attached (ugly) patch and see if the problems are gone?

[bozz64]

The digi tune is now played correctly with this patch but this tune plays really weird now:

/MUSICIANS/P/Psych858o/90_Seconds.sid

[Leandro Nini]

What do you mean by weird? It sounds good to me.

[bozz64]

I applied the patch again and now it sounds good. Perhaps I made a mistake in applying the patch. I tested the digi tune again. All is good now. I will keep testing the coming days.

[Leandro Nini]

Sorry, I think I got the scaling wrong in the above patch, will send and update later.