Diode clipper simulation instability

[Z1202]

Hello Z1202

So, what you said implies that itoa did not implement the TPT the right way, since stiffness should never be a problem when the trapezoidal integrator is involved.

I'm not sure we're talking about the same thing. Trapezoidal integration can have problems with stiffness. It's just that in my experience they are likely to occur with positive rather than negative feedback. This is when you move close to or past the pole of the BLT (this is described in 3.13 of VAFilterDesign). So, my guess would be that the stiffness problem is due to the convergence problems of N/R but "never" is a little bit too strong word

[mystran]

"Stiffness" refers to the problem where your system has very different time constants with the short ones force a tiny integration time step with a conditionally stable integrator, while contributing next to nothing to the actual solution (since they decay out very quickly).

If you use an A-stable implicit integration scheme, then you can more or less ignore the whole stiffness thing in the traditional sense, since the integration will be stable whenever the system itself is stable (there is a slight gotcha with trapezoidal though, namely it has no damping at the limit to the infinity).

Now, I'm not convinced the "positive feedback" issue is really within the realm of stiffness in the traditional sense, since if I'm not mistaken you are referring to the problem where you are trying to solve an unstable system (or unstable linearization of a forced non-linear system), which can fail to have a finite solution (at least temporarily during Newton) and more importantly can make a mess out of simple LU scheme (please raise your hand if you're excited about pivoting at run-time). That's definitely an interesting problem on it's own right, but I wouldn't mix it with "the stiffness problem" since that one happens with stable linear systems.

[andy-cytomic]

I am trying to simulate attached diode limiter circuit.

Still learning electronics, fortunately found ODE for this model:

dVo/dt = (Vi-Vo) /RC - 2 Is/C * sinh(Vo/Vt)

- Found typical values for reverse current and thermal voltage
- Integrated using trapezoidal integrator and Newton-Rhapson

It works and clips nicely (avg 4 iterations) until ...
certain input level get exceeded, then it explodes. In particular it depends on input level / RC ratio .. The bigger input signal is, the lower cutoff freq is needed to make it stable.

I suppose I should limit the voltage somewhere.. as in real circuits.

Could someone help me and explain what is going on here?

Welcome to the wonderful world of circuit solving! What you have is a stiff equations to solve, so directly applying unmodified Newton-Raphson and Trapezoidal integration won't work. There are loads of things to try to fix the issue, but working out a method that converges fast and with low cpu isn't easy. Any methods you do work out can later be applied to something like the K-Method / DK-Method as the MNA solver stage, so any time you spend here is well invested.

Part of the problem can be overcome by not using Trapezoidal integration. Try using Backwards Euler for starters, if you can't solve the circuit with that then you know the integration scheme isn't the issue and you need to look into how Newton Raphson is failing instead. Backwards Euler is A-Stable and L-Stable (http://en.wikipedia.org/wiki/Backward_Euler_method , Trapezoidal is A-Stable but not L-Stable http://en.wikipedia.org/wiki/Trapezoida ... equations) , which is where the problem lies: http://en.wikipedia.org/wiki/Stiff_equation

I hope that helps! (and sorry if this is a repeat to what other people have already posted, I haven't read the thread in depth)

PS: I would also ditch one of the diodes to start with so you are back to a single exp function instead of an sinh, that will simplify things somewhat and allow you to focus on the core issues.

[andy-cytomic]

... But since you use a trapezoidal integrator, stiffness shouldn't be a problem, so I guess it's Newton Raphson. Or a standard algorithm bug ...

Trapezoidal integration isn't L-Stable, so stiffness is a problem: http://en.wikipedia.org/wiki/Stiff_equation

[aciddose]

My suggestion to use a table was based upon:

btw. Just realised that these RC values means 7MHz cutoff SIC?!

Although the values seem to be 2.2k, 10n = 7234hz, which would make the filter significant in the audible range and although the effect due to the diode is minimal in this case, it could in theory be audible especially for higher frequencies.

Code: Select all

*

R1 in lp 2200
C1 lp gnd 10n
B1 buf gnd V = V(lp)
R2 buf rcbd 2200
X1 rcbd gnd clamp

R3 in rcd 2200
C2 rcd gnd 10n
X2 rcd gnd clamp

V1 in 0 DC 0 SIN(0 5 1000hz)

.tran 1n 20m 0 100n

.subckt clamp a b
D1 a b DX
D2 b a DX
.ends

.model DX D()

.end

I highly recommend looking at applying a spice simulation tool before committing significant resources to writing an accurate model. In many cases you can get 99% there without being perfectly accurate, and the extra 0.999% isn't always worth it.

For example at 7.5k the first harmonic of a sine, the difference between the rc -> buffer -> diode and rc+diode is 1.25db. This isn't exactly audible.

At 1k the picture looks like this:
{img}http://xhip.net/temp/1krcd.png{/img}

Doesn't look very audible...

I'm not 100% certain of these results (exactly how ltspice manages fft and how much influence various configurations have on results) but it doesn't look very convincing from the "is this worth modelling?" standpoint.

[Ivan_C]

Well, for me the fact that trapezoidal integration scheme has not L-stability isn't an issue when you deal with things which have an influence on the audio frequencies range right ? I won't say the same about ringing, Newton-Raphson's method convergence, or the stuff about how to handle unstable systems.

However, since I was able to simulate this stuff without having to deal with any of these issues (thanks to the use of a lookup table for the convergence), I suggested that itoa may have been wrong on something else But, the use of the word "never" for stiffness and trapezoidal was indeed too strong

I highly recommend looking at applying a spice simulation tool before committing significant resources to writing an accurate model. In many cases you can get 99% there without being perfectly accurate, and the extra 0.999% isn't always worth it.

+1

[andy-cytomic]

I agree that using simulators like spice and qucs are very useful to give you an idea of how much difference it will make to include different levels of modelling, but to make any final choices you really need to listen to both and decide for yourself. If standard DAW stomp box models are anything to go by then the answer is no, don't bother, but then again if you are going to compete with standard DAW plugins you need to offer something extra.

This simple diode shaper if modelled "the easy way" by buffering the low pass from the diode clipper will sound very similar to the full model and be a huge amount easier to code, since it involves a single one pole low pass filter and then a waveshaper.

The main difference in tone will happen to the high frequencies of your input signal, 5 khz and up. The full model will have a slightly brighter tone on the high frequencies that dynamically depends on the volume of the input signal, which is exactly the sort of thing guitarists like. But it is a lot more effort and a lot more cpu than the trivial model.

[aciddose]

In this case however you can clearly see what happens and what frequency the filter sits at. So for any model with this particular configuration you're talking at most 1db difference in the audible range.

I'm not certain if the sine test input produces the most radical difference, but if someone was worried about weird wigglywave and what spectra it might generate they can go right ahead and try that with the netlist I posted.

It just seems very difficult to me to justify this sort of model until you get at least a 3db difference, and even then additional filter tweaks can make up that difference in a far less expensive and audibly identical way vs. an iterative model.

[andy-cytomic]

Yes, for this very basic case there are some obvious hacks that will get you very close to the full model. I recommend letting people work it out for themselves by looking at what is going on with the equations, this way people will learn a lot more and then be better positioned to tackle more difficult problems. The hint here is to realise that the diode is just a voltage controlled resistor, when there is a small voltage across them they are a very large resistor, when the voltage over them is large they are a small valued resistor. So try replacing the diode with various values of resistor and plot what happens to the frequency response of the system.

But there is no need to resort to a hack since you can already get an exact solution in one step. You can use a sheared multi-dimensional non-linearity to solve this system directly. This kind of thing is generally called the K-Method or DK-Method and is perfect for this, the most basic of cases, since the dimensionality is low and the resultant manifold easy to fit.

[aciddose]

One step will still be more expensive than a lerp table lookup + native lossy integrator.

In fact it is possible to replace the lossy integrator completely with a FIR filter and compute the table loopup and FIR using SSE. This is critical if you want to combine it with oversampling to handle the aliasing.

[andy-cytomic]

One step will still be more expensive than a lerp table lookup + native lossy integrator...

The method I have in mind is a lookup and an increment of a single state variable for the capacitor, so it is slightly more cpu for the lookup since it is not 1D, but otherwise uses the same cpu as the trivial method.

[aciddose]

Well then you end up with a trade-off between various properties such as hf oscillation vs. damping that becomes fairly complex once the whole system including oversampling and later stages is taken into account.

It becomes very hard to say what would be ideal.

For example take this system followed by another identical system. We would want to minimize hf oscillation because any oscillation near nyquist is going to increase the level of aliasing. We would prefer significant damping instead on the first stage and then possibly some makeup on the second stage. Accuracy within a couple dB at nyquist vs. significant aliasing.

[mystran]

One step will still be more expensive than a lerp table lookup + native lossy integrator...

The method I have in mind is a lookup and an increment of a single state variable for the capacitor, so it is slightly more cpu for the lookup since it is not 1D, but otherwise uses the same cpu as the trivial method.

Actually, since you can separate the non-linear part from the dynamics and since the system only has a single one-dimensional non-linearity (if you lump the diodes together) you can use a simple 1D lookup if you want.

edit: I'm not sure how that helps with the Newton problem though, since as far as I can see, you still need a working iterative scheme to actually populate the table.

[andy-cytomic]

One step will still be more expensive than a lerp table lookup + native lossy integrator...

The method I have in mind is a lookup and an increment of a single state variable for the capacitor, so it is slightly more cpu for the lookup since it is not 1D, but otherwise uses the same cpu as the trivial method.

Actually, since you can separate the non-linear part from the dynamics and since the system only has a single one-dimensional non-linearity (if you lump the diodes together) you can use a simple 1D lookup if you want.

edit: I'm not sure how that helps with the Newton problem though, since as far as I can see, you still need a working iterative scheme to actually populate the table.

Yes, I was thinking of the cap being on the other end of the diodes (like an envelope follower), this is only a 1D lookup since the contribution of current from the input resistor and capacitor can be summed.

Forming this table requires a regular full solution, which is why I previously stated that time is not wasted by working on getting the full version to solve properly.

[JCJR]

The hint here is to realise that the diode is just a voltage controlled resistor, when there is a small voltage across them they are a very large resistor, when the voltage over them is large they are a small valued resistor. So try replacing the diode with various values of resistor and plot what happens to the frequency response.

Apologies belaboring a "minor" point-- Perhaps plotting with different-valued resistors would get close enough, but this wild "effective resistance" swing happens during every cycle, and so the simple plotting of different linear frequency responses for different static resistance substitutions, might not be very close at all to reality. Dunno how dramatic would be the variance, having never tried solving it.

The discharge rate of the circuit will constantly change across the waveform, according to the input and historical state. Which might be significantly different amplitude-dependent frequency response compared to a table of responses calculated from a variety of different-value static linear systems.

Signed Captain Obvious