1 //  ---------------------------------------------------------------------------
2 //  This file is part of reSID, a MOS6581 SID emulator engine.
3 //  Copyright (C) 2010  Dag Lem <resid@nimrod.no>
4 //
5 //  This program is free software; you can redistribute it and/or modify
6 //  it under the terms of the GNU General Public License as published by
7 //  the Free Software Foundation; either version 2 of the License, or
8 //  (at your option) any later version.
9 //
10 //  This program is distributed in the hope that it will be useful,
11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 //  GNU General Public License for more details.
14 //
15 //  You should have received a copy of the GNU General Public License
16 //  along with this program; if not, write to the Free Software
17 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18 //  ---------------------------------------------------------------------------
19 
20 #ifndef RESID_FILTER_H
21 #define RESID_FILTER_H
22 
23 #include "resid-config.h"
24 
25 namespace reSID
26 {
27 
28 // ----------------------------------------------------------------------------
29 // The SID filter is modeled with a two-integrator-loop biquadratic filter,
30 // which has been confirmed by Bob Yannes to be the actual circuit used in
31 // the SID chip.
32 //
33 // Measurements show that excellent emulation of the SID filter is achieved,
34 // except when high resonance is combined with high sustain levels.
35 // In this case the SID op-amps are performing less than ideally and are
36 // causing some peculiar behavior of the SID filter. This however seems to
37 // have more effect on the overall amplitude than on the color of the sound.
38 //
39 // The theory for the filter circuit can be found in "Microelectric Circuits"
40 // by Adel S. Sedra and Kenneth C. Smith.
41 // The circuit is modeled based on the explanation found there except that
42 // an additional inverter is used in the feedback from the bandpass output,
43 // allowing the summer op-amp to operate in single-ended mode. This yields
44 // filter outputs with levels independent of Q, which corresponds with the
45 // results obtained from a real SID.
46 //
47 // We have been able to model the summer and the two integrators of the circuit
48 // to form components of an IIR filter.
49 // Vhp is the output of the summer, Vbp is the output of the first integrator,
50 // and Vlp is the output of the second integrator in the filter circuit.
51 //
52 // According to Bob Yannes, the active stages of the SID filter are not really
53 // op-amps. Rather, simple NMOS inverters are used. By biasing an inverter
54 // into its region of quasi-linear operation using a feedback resistor from
55 // input to output, a MOS inverter can be made to act like an op-amp for
56 // small signals centered around the switching threshold.
57 //
58 // In 2008, Michael Huth facilitated closer investigation of the SID 6581
59 // filter circuit by publishing high quality microscope photographs of the die.
60 // Tommi Lempinen has done an impressive work on re-vectorizing and annotating
61 // the die photographs, substantially simplifying further analysis of the
62 // filter circuit.
63 //
64 // The filter schematics below are reverse engineered from these re-vectorized
65 // and annotated die photographs. While the filter first depicted in reSID 0.9
66 // is a correct model of the basic filter, the schematics are now completed
67 // with the audio mixer and output stage, including details on intended
68 // relative resistor values. Also included are schematics for the NMOS FET
69 // voltage controlled resistors (VCRs) used to control cutoff frequency, the
70 // DAC which controls the VCRs, the NMOS op-amps, and the output buffer.
71 //
72 //
73 // SID filter / mixer / output
74 // ---------------------------
75 //
76 //                ---------------------------------------------------
77 //               |                                                   |
78 //               |                         --1R1-- \--  D7           |
79 //               |              ---R1--   |           |              |
80 //               |             |       |  |--2R1-- \--| D6           |
81 //               |    ------------<A]-----|           |     $17      |
82 //               |   |                    |--4R1-- \--| D5  1=open   | (3.5R1)
83 //               |   |                    |           |              |
84 //               |   |                     --8R1-- \--| D4           | (7.0R1)
85 //               |   |                                |              |
86 // $17           |   |                    (CAP2B)     |  (CAP1B)     |
87 // 0=to mixer    |    --R8--    ---R8--        ---C---|       ---C---|
88 // 1=to filter   |          |  |       |      |       |      |       |
89 //                ------R8--|-----[A>--|--Rw-----[A>--|--Rw-----[A>--|
90 //     ve (EXT IN)          |          |              |              |
91 // D3  \ ---------------R8--|          |              | (CAP2A)      | (CAP1A)
92 //     |   v3               |          | vhp          | vbp          | vlp
93 // D2  |   \ -----------R8--|     -----               |              |
94 //     |   |   v2           |    |                    |              |
95 // D1  |   |   \ -------R8--|    |    ----------------               |
96 //     |   |   |   v1       |    |   |                               |
97 // D0  |   |   |   \ ---R8--     |   |    ---------------------------
98 //     |   |   |   |             |   |   |
99 //     R6  R6  R6  R6            R6  R6  R6
100 //     |   |   |   | $18         |   |   |  $18
101 //     |    \  |   | D7: 1=open   \   \   \ D6 - D4: 0=open
102 //     |   |   |   |             |   |   |
103 //      ---------------------------------                          12V
104 //                 |
105 //                 |               D3  --/ --1R2--                  |
106 //                 |    ---R8--       |           |   ---R2--       |
107 //                 |   |       |   D2 |--/ --2R2--|  |       |  ||--
108 //                  ------[A>---------|           |-----[A>-----||
109 //                                 D1 |--/ --4R2--| (4.25R2)    ||--
110 //                        $18         |           |                 |
111 //                        0=open   D0  --/ --8R2--  (8.75R2)        |
112 //
113 //                                                                  vo (AUDIO
114 //                                                                      OUT)
115 //
116 //
117 // v1  - voice 1
118 // v2  - voice 2
119 // v3  - voice 3
120 // ve  - ext in
121 // vhp - highpass output
122 // vbp - bandpass output
123 // vlp - lowpass output
124 // vo  - audio out
125 // [A> - single ended inverting op-amp (self-biased NMOS inverter)
126 // Rn  - "resistors", implemented with custom NMOS FETs
127 // Rw  - cutoff frequency resistor (VCR)
128 // C   - capacitor
129 //
130 // Notes:
131 //
132 // R2  ~  2.0*R1
133 // R6  ~  6.0*R1
134 // R8  ~  8.0*R1
135 // R24 ~ 24.0*R1
136 //
137 // The Rn "resistors" in the circuit are implemented with custom NMOS FETs,
138 // probably because of space constraints on the SID die. The silicon substrate
139 // is laid out in a narrow strip or "snake", with a strip length proportional
140 // to the intended resistance. The polysilicon gate electrode covers the entire
141 // silicon substrate and is fixed at 12V in order for the NMOS FET to operate
142 // in triode mode (a.k.a. linear mode or ohmic mode).
143 //
144 // Even in "linear mode", an NMOS FET is only an approximation of a resistor,
145 // as the apparant resistance increases with increasing drain-to-source
146 // voltage. If the drain-to-source voltage should approach the gate voltage
147 // of 12V, the NMOS FET will enter saturation mode (a.k.a. active mode), and
148 // the NMOS FET will not operate anywhere like a resistor.
149 //
150 //
151 //
152 // NMOS FET voltage controlled resistor (VCR)
153 // ------------------------------------------
154 //
155 //                Vw
156 //
157 //                |
158 //                |
159 //                R1
160 //                |
161 //          --R1--|
162 //         |    __|__
163 //         |    -----
164 //         |    |   |
165 // vi ----------     -------- vo
166 //         |           |
167 //          ----R24----
168 //
169 //
170 // vi  - input
171 // vo  - output
172 // Rn  - "resistors", implemented with custom NMOS FETs
173 // Vw  - voltage from 11-bit DAC (frequency cutoff control)
174 //
175 // Notes:
176 //
177 // An approximate value for R24 can be found by using the formula for the
178 // filter cutoff frequency:
179 //
180 // FCmin = 1/(2*pi*Rmax*C)
181 //
182 // Assuming that a the setting for minimum cutoff frequency in combination with
183 // a low level input signal ensures that only negligible current will flow
184 // through the transistor in the schematics above, values for FCmin and C can
185 // be substituted in this formula to find Rmax.
186 // Using C = 470pF and FCmin = 220Hz (measured value), we get:
187 //
188 // FCmin = 1/(2*pi*Rmax*C)
189 // Rmax = 1/(2*pi*FCmin*C) = 1/(2*pi*220*470e-12) ~ 1.5MOhm
190 //
191 // From this it follows that:
192 // R24 =  Rmax   ~ 1.5MOhm
193 // R1  ~  R24/24 ~  64kOhm
194 // R2  ~  2.0*R1 ~ 128kOhm
195 // R6  ~  6.0*R1 ~ 384kOhm
196 // R8  ~  8.0*R1 ~ 512kOhm
197 //
198 // Note that these are only approximate values for one particular SID chip,
199 // due to process variations the values can be substantially different in
200 // other chips.
201 //
202 //
203 //
204 // Filter frequency cutoff DAC
205 // ---------------------------
206 //
207 //
208 //        12V  10   9   8   7   6   5   4   3   2   1   0   VGND
209 //          |   |   |   |   |   |   |   |   |   |   |   |     |   Missing
210 //         2R  2R  2R  2R  2R  2R  2R  2R  2R  2R  2R  2R    2R   termination
211 //          |   |   |   |   |   |   |   |   |   |   |   |     |
212 //     Vw ----R---R---R---R---R---R---R---R---R---R---R--   ---
213 //
214 // Bit on:  12V
215 // Bit off:  5V (VGND)
216 //
217 // As is the case with all MOS 6581 DACs, the termination to (virtual) ground
218 // at bit 0 is missing.
219 //
220 // Furthermore, the control of the two VCRs imposes a load on the DAC output
221 // which varies with the input signals to the VCRs. This can be seen from the
222 // VCR figure above.
223 //
224 //
225 //
226 // "Op-amp" (self-biased NMOS inverter)
227 // ------------------------------------
228 //
229 //
230 //                        12V
231 //
232 //                         |
233 //              -----------|
234 //             |           |
235 //             |     ------|
236 //             |    |      |
237 //             |    |  ||--
238 //             |     --||
239 //             |       ||--
240 //         ||--            |
241 // vi -----||              |--------- vo
242 //         ||--            |   |
243 //             |       ||--    |
244 //             |-------||      |
245 //             |       ||--    |
246 //         ||--            |   |
247 //       --||              |   |
248 //      |  ||--            |   |
249 //      |      |           |   |
250 //      |       -----------|   |
251 //      |                  |   |
252 //      |                      |
253 //      |                 GND  |
254 //      |                      |
255 //       ----------------------
256 //
257 //
258 // vi  - input
259 // vo  - output
260 //
261 // Notes:
262 //
263 // The schematics above are laid out to show that the "op-amp" logically
264 // consists of two building blocks; a saturated load NMOS inverter (on the
265 // right hand side of the schematics) with a buffer / bias input stage
266 // consisting of a variable saturated load NMOS inverter (on the left hand
267 // side of the schematics).
268 //
269 // Provided a reasonably high input impedance and a reasonably low output
270 // impedance, the "op-amp" can be modeled as a voltage transfer function
271 // mapping input voltage to output voltage.
272 //
273 //
274 //
275 // Output buffer (NMOS voltage follower)
276 // -------------------------------------
277 //
278 //
279 //            12V
280 //
281 //             |
282 //             |
283 //         ||--
284 // vi -----||
285 //         ||--
286 //             |
287 //             |------ vo
288 //             |     (AUDIO
289 //            Rext    OUT)
290 //             |
291 //             |
292 //
293 //            GND
294 //
295 // vi   - input
296 // vo   - output
297 // Rext - external resistor, 1kOhm
298 //
299 // Notes:
300 //
301 // The external resistor Rext is needed to complete the NMOS voltage follower,
302 // this resistor has a recommended value of 1kOhm.
303 //
304 // Die photographs show that actually, two NMOS transistors are used in the
305 // voltage follower. However the two transistors are coupled in parallel (all
306 // terminals are pairwise common), which implies that we can model the two
307 // transistors as one.
308 //
309 // ----------------------------------------------------------------------------
310 
311 // Compile-time computation of op-amp summer and mixer table offsets.
312 
313 // The highpass summer has 2 - 6 inputs (bandpass, lowpass, and 0 - 4 voices).
314 template<int i>
315 struct summer_offset
316 {
317   enum { value = summer_offset<i - 1>::value + ((2 + i - 1) << 16) };
318 };
319 
320 template<>
321 struct summer_offset<0>
322 {
323   enum { value = 0 };
324 };
325 
326 // The mixer has 0 - 7 inputs (0 - 4 voices and 0 - 3 filter outputs).
327 template<int i>
328 struct mixer_offset
329 {
330   enum { value = mixer_offset<i - 1>::value + ((i - 1) << 16) };
331 };
332 
333 template<>
334 struct mixer_offset<1>
335 {
336   enum { value = 1 };
337 };
338 
339 template<>
340 struct mixer_offset<0>
341 {
342   enum { value = 0 };
343 };
344 
345 
346 class Filter
347 {
348 public:
349   Filter();
350 
351   void enable_filter(bool enable);
352   void adjust_filter_bias(double dac_bias);
353   void set_chip_model(chip_model model);
354   void set_voice_mask(reg4 mask);
355 
356   void clock(int voice1, int voice2, int voice3);
357   void clock(cycle_count delta_t, int voice1, int voice2, int voice3);
358   void reset();
359 
360   // Write registers.
361   void writeFC_LO(reg8);
362   void writeFC_HI(reg8);
363   void writeRES_FILT(reg8);
364   void writeMODE_VOL(reg8);
365 
366   // SID audio input (16 bits).
367   void input(short sample);
368 
369   // SID audio output (16 bits).
370   short output();
371 
372 protected:
373   void set_sum_mix();
374   void set_w0();
375   void set_Q();
376 
377   // Filter enabled.
378   bool enabled;
379 
380   // Filter cutoff frequency.
381   reg12 fc;
382 
383   // Filter resonance.
384   reg8 res;
385 
386   // Selects which voices to route through the filter.
387   reg8 filt;
388 
389   // Selects which filter outputs to route into the mixer.
390   reg4 mode;
391 
392   // Output master volume.
393   reg4 vol;
394 
395   // Used to mask out EXT IN if not connected, and for test purposes
396   // (voice muting).
397   reg8 voice_mask;
398 
399   // Select which inputs to route into the summer / mixer.
400   // These are derived from filt, mode, and voice_mask.
401   reg8 sum;
402   reg8 mix;
403 
404   // State of filter.
405   int Vhp; // highpass
406   int Vbp; // bandpass
407   int Vbp_x, Vbp_vc;
408   int Vlp; // lowpass
409   int Vlp_x, Vlp_vc;
410   // Filter / mixer inputs.
411   int ve;
412   int v3;
413   int v2;
414   int v1;
415 
416   // Cutoff frequency DAC voltage, resonance.
417   int Vddt_Vw_2, Vw_bias;
418   int _8_div_Q;
419   // FIXME: Temporarily used for MOS 8580 emulation.
420   int w0;
421   int _1024_div_Q;
422 
423   chip_model sid_model;
424 
425    typedef struct {
426     unsigned short vx;
427     short dvx;
428   } opamp_t;
429 
430   typedef struct {
431     int vo_N16;  // Fixed point scaling for 16 bit op-amp output.
432     int kVddt;   // K*(Vdd - Vth)
433     int n_snake;
434     int voice_scale_s14;
435     int voice_DC;
436     int ak;
437     int bk;
438     int vc_min;
439     int vc_max;
440 
441     // Reverse op-amp transfer function.
442     unsigned short opamp_rev[1 << 16];
443     // Lookup tables for gain and summer op-amps in output stage / filter.
444     unsigned short summer[summer_offset<5>::value];
445     unsigned short gain[16][1 << 16];
446     unsigned short mixer[mixer_offset<8>::value];
447     // Cutoff frequency DAC output voltage table. FC is an 11 bit register.
448     unsigned short f0_dac[1 << 11];
449   } model_filter_t;
450 
451   int solve_gain(opamp_t* opamp, int n, int vi_t, int& x, model_filter_t& mf);
452   int solve_integrate_6581(int dt, int vi_t, int& x, int& vc, model_filter_t& mf);
453 
454   // VCR - 6581 only.
455   static unsigned short vcr_kVg[1 << 16];
456   static unsigned short vcr_n_Ids_term[1 << 16];
457   // Common parameters.
458   static model_filter_t model_filter[2];
459 
460 friend class SID;
461 };
462 
463 
464 // ----------------------------------------------------------------------------
465 // Inline functions.
466 // The following functions are defined inline because they are called every
467 // time a sample is calculated.
468 // ----------------------------------------------------------------------------
469 
470 #if RESID_INLINING || defined(RESID_FILTER_CC)
471 
472 // ----------------------------------------------------------------------------
473 // SID clocking - 1 cycle.
474 // ----------------------------------------------------------------------------
475 RESID_INLINE
476 void Filter::clock(int voice1, int voice2, int voice3)
477 {
478   model_filter_t& f = model_filter[sid_model];
479 
480   v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
481   v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
482   v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
483 
484   // Sum inputs routed into the filter.
485   int Vi = 0;
486   int offset = 0;
487 
488   switch (sum & 0xf) {
489   case 0x0:
490     Vi = 0;
491     offset = summer_offset<0>::value;
492     break;
493   case 0x1:
494     Vi = v1;
495     offset = summer_offset<1>::value;
496     break;
497   case 0x2:
498     Vi = v2;
499     offset = summer_offset<1>::value;
500     break;
501   case 0x3:
502     Vi = v2 + v1;
503     offset = summer_offset<2>::value;
504     break;
505   case 0x4:
506     Vi = v3;
507     offset = summer_offset<1>::value;
508     break;
509   case 0x5:
510     Vi = v3 + v1;
511     offset = summer_offset<2>::value;
512     break;
513   case 0x6:
514     Vi = v3 + v2;
515     offset = summer_offset<2>::value;
516     break;
517   case 0x7:
518     Vi = v3 + v2 + v1;
519     offset = summer_offset<3>::value;
520     break;
521   case 0x8:
522     Vi = ve;
523     offset = summer_offset<1>::value;
524     break;
525   case 0x9:
526     Vi = ve + v1;
527     offset = summer_offset<2>::value;
528     break;
529   case 0xa:
530     Vi = ve + v2;
531     offset = summer_offset<2>::value;
532     break;
533   case 0xb:
534     Vi = ve + v2 + v1;
535     offset = summer_offset<3>::value;
536     break;
537   case 0xc:
538     Vi = ve + v3;
539     offset = summer_offset<2>::value;
540     break;
541   case 0xd:
542     Vi = ve + v3 + v1;
543     offset = summer_offset<3>::value;
544     break;
545   case 0xe:
546     Vi = ve + v3 + v2;
547     offset = summer_offset<3>::value;
548     break;
549   case 0xf:
550     Vi = ve + v3 + v2 + v1;
551     offset = summer_offset<4>::value;
552     break;
553   }
554 
555   // Calculate filter outputs.
556   if (sid_model == 0) {
557     // MOS 6581.
558     Vlp = solve_integrate_6581(1, Vbp, Vlp_x, Vlp_vc, f);
559     Vbp = solve_integrate_6581(1, Vhp, Vbp_x, Vbp_vc, f);
560     Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
561   }
562   else {
563     // MOS 8580. FIXME: Not yet using op-amp model.
564 
565     // delta_t = 1 is converted to seconds given a 1MHz clock by dividing
566     // with 1 000 000.
567 
568     int dVbp = w0*(Vhp >> 4) >> 16;
569     int dVlp = w0*(Vbp >> 4) >> 16;
570     Vbp -= dVbp;
571     Vlp -= dVlp;
572     Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
573   }
574 }
575 
576 // ----------------------------------------------------------------------------
577 // SID clocking - delta_t cycles.
578 // ----------------------------------------------------------------------------
579 RESID_INLINE
580 void Filter::clock(cycle_count delta_t, int voice1, int voice2, int voice3)
581 {
582   model_filter_t& f = model_filter[sid_model];
583 
584   v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
585   v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
586   v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
587 
588   // Enable filter on/off.
589   // This is not really part of SID, but is useful for testing.
590   // On slow CPUs it may be necessary to bypass the filter to lower the CPU
591   // load.
592   if (unlikely(!enabled)) {
593     return;
594   }
595 
596   // Sum inputs routed into the filter.
597   int Vi = 0;
598   int offset = 0;
599 
600   switch (sum & 0xf) {
601   case 0x0:
602     Vi = 0;
603     offset = summer_offset<0>::value;
604     break;
605   case 0x1:
606     Vi = v1;
607     offset = summer_offset<1>::value;
608     break;
609   case 0x2:
610     Vi = v2;
611     offset = summer_offset<1>::value;
612     break;
613   case 0x3:
614     Vi = v2 + v1;
615     offset = summer_offset<2>::value;
616     break;
617   case 0x4:
618     Vi = v3;
619     offset = summer_offset<1>::value;
620     break;
621   case 0x5:
622     Vi = v3 + v1;
623     offset = summer_offset<2>::value;
624     break;
625   case 0x6:
626     Vi = v3 + v2;
627     offset = summer_offset<2>::value;
628     break;
629   case 0x7:
630     Vi = v3 + v2 + v1;
631     offset = summer_offset<3>::value;
632     break;
633   case 0x8:
634     Vi = ve;
635     offset = summer_offset<1>::value;
636     break;
637   case 0x9:
638     Vi = ve + v1;
639     offset = summer_offset<2>::value;
640     break;
641   case 0xa:
642     Vi = ve + v2;
643     offset = summer_offset<2>::value;
644     break;
645   case 0xb:
646     Vi = ve + v2 + v1;
647     offset = summer_offset<3>::value;
648     break;
649   case 0xc:
650     Vi = ve + v3;
651     offset = summer_offset<2>::value;
652     break;
653   case 0xd:
654     Vi = ve + v3 + v1;
655     offset = summer_offset<3>::value;
656     break;
657   case 0xe:
658     Vi = ve + v3 + v2;
659     offset = summer_offset<3>::value;
660     break;
661   case 0xf:
662     Vi = ve + v3 + v2 + v1;
663     offset = summer_offset<4>::value;
664     break;
665   }
666 
667   // Maximum delta cycles for filter fixpoint iteration to converge
668   // is approximately 3.
669   cycle_count delta_t_flt = 3;
670 
671   if (sid_model == 0) {
672     // MOS 6581.
673     while (delta_t) {
674       if (unlikely(delta_t < delta_t_flt)) {
675         delta_t_flt = delta_t;
676       }
677 
678       // Calculate filter outputs.
679       Vlp = solve_integrate_6581(delta_t_flt, Vbp, Vlp_x, Vlp_vc, f);
680       Vbp = solve_integrate_6581(delta_t_flt, Vhp, Vbp_x, Vbp_vc, f);
681       Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
682 
683       delta_t -= delta_t_flt;
684     }
685   }
686   else {
687     // MOS 8580. FIXME: Not yet using op-amp model.
688     while (delta_t) {
689       if (delta_t < delta_t_flt) {
690         delta_t_flt = delta_t;
691       }
692 
693       // delta_t is converted to seconds given a 1MHz clock by dividing
694       // with 1 000 000. This is done in two operations to avoid integer
695       // multiplication overflow.
696 
697       // Calculate filter outputs.
698       int w0_delta_t = w0*delta_t_flt >> 2;
699 
700       int dVbp = w0_delta_t*(Vhp >> 4) >> 14;
701       int dVlp = w0_delta_t*(Vbp >> 4) >> 14;
702       Vbp -= dVbp;
703       Vlp -= dVlp;
704       Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
705 
706       delta_t -= delta_t_flt;
707     }
708   }
709 }
710 
711 
712 // ----------------------------------------------------------------------------
713 // SID audio input (16 bits).
714 // ----------------------------------------------------------------------------
715 RESID_INLINE
716 void Filter::input(short sample)
717 {
718   // Scale to three times the peak-to-peak for one voice and add the op-amp
719   // "zero" DC level.
720   // NB! Adding the op-amp "zero" DC level is a (wildly inaccurate)
721   // approximation of feeding the input through an AC coupling capacitor.
722   // This could be implemented as a separate filter circuit, however the
723   // primary use of the emulator is not to process external signals.
724   // The upside is that the MOS8580 "digi boost" works without a separate (DC)
725   // input interface.
726   // Note that the input is 16 bits, compared to the 20 bit voice output.
727   model_filter_t& f = model_filter[sid_model];
728   ve = (sample*f.voice_scale_s14*3 >> 14) + f.mixer[0];
729 }
730 
731 
732 // ----------------------------------------------------------------------------
733 // SID audio output (16 bits).
734 // ----------------------------------------------------------------------------
735 RESID_INLINE
736 short Filter::output()
737 {
738   model_filter_t& f = model_filter[sid_model];
739 
740   // Writing the switch below manually would be tedious and error-prone;
741   // it is rather generated by the following Perl program:
742 
743   /*
744 my @i = qw(v1 v2 v3 ve Vlp Vbp Vhp);
745 for my $mix (0..2**@i-1) {
746     print sprintf("  case 0x%02x:\n", $mix);
747     my @sum;
748     for (@i) {
749         unshift(@sum, $_) if $mix & 0x01;
750         $mix >>= 1;
751     }
752     my $sum = join(" + ", @sum) || "0";
753     print "    Vi = $sum;\n";
754     print "    offset = mixer_offset<" . @sum . ">::value;\n";
755     print "    break;\n";
756 }
757   */
758 
759   // Sum inputs routed into the mixer.
760   int Vi = 0;
761   int offset = 0;
762 
763   switch (mix & 0x7f) {
764   case 0x00:
765     Vi = 0;
766     offset = mixer_offset<0>::value;
767     break;
768   case 0x01:
769     Vi = v1;
770     offset = mixer_offset<1>::value;
771     break;
772   case 0x02:
773     Vi = v2;
774     offset = mixer_offset<1>::value;
775     break;
776   case 0x03:
777     Vi = v2 + v1;
778     offset = mixer_offset<2>::value;
779     break;
780   case 0x04:
781     Vi = v3;
782     offset = mixer_offset<1>::value;
783     break;
784   case 0x05:
785     Vi = v3 + v1;
786     offset = mixer_offset<2>::value;
787     break;
788   case 0x06:
789     Vi = v3 + v2;
790     offset = mixer_offset<2>::value;
791     break;
792   case 0x07:
793     Vi = v3 + v2 + v1;
794     offset = mixer_offset<3>::value;
795     break;
796   case 0x08:
797     Vi = ve;
798     offset = mixer_offset<1>::value;
799     break;
800   case 0x09:
801     Vi = ve + v1;
802     offset = mixer_offset<2>::value;
803     break;
804   case 0x0a:
805     Vi = ve + v2;
806     offset = mixer_offset<2>::value;
807     break;
808   case 0x0b:
809     Vi = ve + v2 + v1;
810     offset = mixer_offset<3>::value;
811     break;
812   case 0x0c:
813     Vi = ve + v3;
814     offset = mixer_offset<2>::value;
815     break;
816   case 0x0d:
817     Vi = ve + v3 + v1;
818     offset = mixer_offset<3>::value;
819     break;
820   case 0x0e:
821     Vi = ve + v3 + v2;
822     offset = mixer_offset<3>::value;
823     break;
824   case 0x0f:
825     Vi = ve + v3 + v2 + v1;
826     offset = mixer_offset<4>::value;
827     break;
828   case 0x10:
829     Vi = Vlp;
830     offset = mixer_offset<1>::value;
831     break;
832   case 0x11:
833     Vi = Vlp + v1;
834     offset = mixer_offset<2>::value;
835     break;
836   case 0x12:
837     Vi = Vlp + v2;
838     offset = mixer_offset<2>::value;
839     break;
840   case 0x13:
841     Vi = Vlp + v2 + v1;
842     offset = mixer_offset<3>::value;
843     break;
844   case 0x14:
845     Vi = Vlp + v3;
846     offset = mixer_offset<2>::value;
847     break;
848   case 0x15:
849     Vi = Vlp + v3 + v1;
850     offset = mixer_offset<3>::value;
851     break;
852   case 0x16:
853     Vi = Vlp + v3 + v2;
854     offset = mixer_offset<3>::value;
855     break;
856   case 0x17:
857     Vi = Vlp + v3 + v2 + v1;
858     offset = mixer_offset<4>::value;
859     break;
860   case 0x18:
861     Vi = Vlp + ve;
862     offset = mixer_offset<2>::value;
863     break;
864   case 0x19:
865     Vi = Vlp + ve + v1;
866     offset = mixer_offset<3>::value;
867     break;
868   case 0x1a:
869     Vi = Vlp + ve + v2;
870     offset = mixer_offset<3>::value;
871     break;
872   case 0x1b:
873     Vi = Vlp + ve + v2 + v1;
874     offset = mixer_offset<4>::value;
875     break;
876   case 0x1c:
877     Vi = Vlp + ve + v3;
878     offset = mixer_offset<3>::value;
879     break;
880   case 0x1d:
881     Vi = Vlp + ve + v3 + v1;
882     offset = mixer_offset<4>::value;
883     break;
884   case 0x1e:
885     Vi = Vlp + ve + v3 + v2;
886     offset = mixer_offset<4>::value;
887     break;
888   case 0x1f:
889     Vi = Vlp + ve + v3 + v2 + v1;
890     offset = mixer_offset<5>::value;
891     break;
892   case 0x20:
893     Vi = Vbp;
894     offset = mixer_offset<1>::value;
895     break;
896   case 0x21:
897     Vi = Vbp + v1;
898     offset = mixer_offset<2>::value;
899     break;
900   case 0x22:
901     Vi = Vbp + v2;
902     offset = mixer_offset<2>::value;
903     break;
904   case 0x23:
905     Vi = Vbp + v2 + v1;
906     offset = mixer_offset<3>::value;
907     break;
908   case 0x24:
909     Vi = Vbp + v3;
910     offset = mixer_offset<2>::value;
911     break;
912   case 0x25:
913     Vi = Vbp + v3 + v1;
914     offset = mixer_offset<3>::value;
915     break;
916   case 0x26:
917     Vi = Vbp + v3 + v2;
918     offset = mixer_offset<3>::value;
919     break;
920   case 0x27:
921     Vi = Vbp + v3 + v2 + v1;
922     offset = mixer_offset<4>::value;
923     break;
924   case 0x28:
925     Vi = Vbp + ve;
926     offset = mixer_offset<2>::value;
927     break;
928   case 0x29:
929     Vi = Vbp + ve + v1;
930     offset = mixer_offset<3>::value;
931     break;
932   case 0x2a:
933     Vi = Vbp + ve + v2;
934     offset = mixer_offset<3>::value;
935     break;
936   case 0x2b:
937     Vi = Vbp + ve + v2 + v1;
938     offset = mixer_offset<4>::value;
939     break;
940   case 0x2c:
941     Vi = Vbp + ve + v3;
942     offset = mixer_offset<3>::value;
943     break;
944   case 0x2d:
945     Vi = Vbp + ve + v3 + v1;
946     offset = mixer_offset<4>::value;
947     break;
948   case 0x2e:
949     Vi = Vbp + ve + v3 + v2;
950     offset = mixer_offset<4>::value;
951     break;
952   case 0x2f:
953     Vi = Vbp + ve + v3 + v2 + v1;
954     offset = mixer_offset<5>::value;
955     break;
956   case 0x30:
957     Vi = Vbp + Vlp;
958     offset = mixer_offset<2>::value;
959     break;
960   case 0x31:
961     Vi = Vbp + Vlp + v1;
962     offset = mixer_offset<3>::value;
963     break;
964   case 0x32:
965     Vi = Vbp + Vlp + v2;
966     offset = mixer_offset<3>::value;
967     break;
968   case 0x33:
969     Vi = Vbp + Vlp + v2 + v1;
970     offset = mixer_offset<4>::value;
971     break;
972   case 0x34:
973     Vi = Vbp + Vlp + v3;
974     offset = mixer_offset<3>::value;
975     break;
976   case 0x35:
977     Vi = Vbp + Vlp + v3 + v1;
978     offset = mixer_offset<4>::value;
979     break;
980   case 0x36:
981     Vi = Vbp + Vlp + v3 + v2;
982     offset = mixer_offset<4>::value;
983     break;
984   case 0x37:
985     Vi = Vbp + Vlp + v3 + v2 + v1;
986     offset = mixer_offset<5>::value;
987     break;
988   case 0x38:
989     Vi = Vbp + Vlp + ve;
990     offset = mixer_offset<3>::value;
991     break;
992   case 0x39:
993     Vi = Vbp + Vlp + ve + v1;
994     offset = mixer_offset<4>::value;
995     break;
996   case 0x3a:
997     Vi = Vbp + Vlp + ve + v2;
998     offset = mixer_offset<4>::value;
999     break;
1000   case 0x3b:
1001     Vi = Vbp + Vlp + ve + v2 + v1;
1002     offset = mixer_offset<5>::value;
1003     break;
1004   case 0x3c:
1005     Vi = Vbp + Vlp + ve + v3;
1006     offset = mixer_offset<4>::value;
1007     break;
1008   case 0x3d:
1009     Vi = Vbp + Vlp + ve + v3 + v1;
1010     offset = mixer_offset<5>::value;
1011     break;
1012   case 0x3e:
1013     Vi = Vbp + Vlp + ve + v3 + v2;
1014     offset = mixer_offset<5>::value;
1015     break;
1016   case 0x3f:
1017     Vi = Vbp + Vlp + ve + v3 + v2 + v1;
1018     offset = mixer_offset<6>::value;
1019     break;
1020   case 0x40:
1021     Vi = Vhp;
1022     offset = mixer_offset<1>::value;
1023     break;
1024   case 0x41:
1025     Vi = Vhp + v1;
1026     offset = mixer_offset<2>::value;
1027     break;
1028   case 0x42:
1029     Vi = Vhp + v2;
1030     offset = mixer_offset<2>::value;
1031     break;
1032   case 0x43:
1033     Vi = Vhp + v2 + v1;
1034     offset = mixer_offset<3>::value;
1035     break;
1036   case 0x44:
1037     Vi = Vhp + v3;
1038     offset = mixer_offset<2>::value;
1039     break;
1040   case 0x45:
1041     Vi = Vhp + v3 + v1;
1042     offset = mixer_offset<3>::value;
1043     break;
1044   case 0x46:
1045     Vi = Vhp + v3 + v2;
1046     offset = mixer_offset<3>::value;
1047     break;
1048   case 0x47:
1049     Vi = Vhp + v3 + v2 + v1;
1050     offset = mixer_offset<4>::value;
1051     break;
1052   case 0x48:
1053     Vi = Vhp + ve;
1054     offset = mixer_offset<2>::value;
1055     break;
1056   case 0x49:
1057     Vi = Vhp + ve + v1;
1058     offset = mixer_offset<3>::value;
1059     break;
1060   case 0x4a:
1061     Vi = Vhp + ve + v2;
1062     offset = mixer_offset<3>::value;
1063     break;
1064   case 0x4b:
1065     Vi = Vhp + ve + v2 + v1;
1066     offset = mixer_offset<4>::value;
1067     break;
1068   case 0x4c:
1069     Vi = Vhp + ve + v3;
1070     offset = mixer_offset<3>::value;
1071     break;
1072   case 0x4d:
1073     Vi = Vhp + ve + v3 + v1;
1074     offset = mixer_offset<4>::value;
1075     break;
1076   case 0x4e:
1077     Vi = Vhp + ve + v3 + v2;
1078     offset = mixer_offset<4>::value;
1079     break;
1080   case 0x4f:
1081     Vi = Vhp + ve + v3 + v2 + v1;
1082     offset = mixer_offset<5>::value;
1083     break;
1084   case 0x50:
1085     Vi = Vhp + Vlp;
1086     offset = mixer_offset<2>::value;
1087     break;
1088   case 0x51:
1089     Vi = Vhp + Vlp + v1;
1090     offset = mixer_offset<3>::value;
1091     break;
1092   case 0x52:
1093     Vi = Vhp + Vlp + v2;
1094     offset = mixer_offset<3>::value;
1095     break;
1096   case 0x53:
1097     Vi = Vhp + Vlp + v2 + v1;
1098     offset = mixer_offset<4>::value;
1099     break;
1100   case 0x54:
1101     Vi = Vhp + Vlp + v3;
1102     offset = mixer_offset<3>::value;
1103     break;
1104   case 0x55:
1105     Vi = Vhp + Vlp + v3 + v1;
1106     offset = mixer_offset<4>::value;
1107     break;
1108   case 0x56:
1109     Vi = Vhp + Vlp + v3 + v2;
1110     offset = mixer_offset<4>::value;
1111     break;
1112   case 0x57:
1113     Vi = Vhp + Vlp + v3 + v2 + v1;
1114     offset = mixer_offset<5>::value;
1115     break;
1116   case 0x58:
1117     Vi = Vhp + Vlp + ve;
1118     offset = mixer_offset<3>::value;
1119     break;
1120   case 0x59:
1121     Vi = Vhp + Vlp + ve + v1;
1122     offset = mixer_offset<4>::value;
1123     break;
1124   case 0x5a:
1125     Vi = Vhp + Vlp + ve + v2;
1126     offset = mixer_offset<4>::value;
1127     break;
1128   case 0x5b:
1129     Vi = Vhp + Vlp + ve + v2 + v1;
1130     offset = mixer_offset<5>::value;
1131     break;
1132   case 0x5c:
1133     Vi = Vhp + Vlp + ve + v3;
1134     offset = mixer_offset<4>::value;
1135     break;
1136   case 0x5d:
1137     Vi = Vhp + Vlp + ve + v3 + v1;
1138     offset = mixer_offset<5>::value;
1139     break;
1140   case 0x5e:
1141     Vi = Vhp + Vlp + ve + v3 + v2;
1142     offset = mixer_offset<5>::value;
1143     break;
1144   case 0x5f:
1145     Vi = Vhp + Vlp + ve + v3 + v2 + v1;
1146     offset = mixer_offset<6>::value;
1147     break;
1148   case 0x60:
1149     Vi = Vhp + Vbp;
1150     offset = mixer_offset<2>::value;
1151     break;
1152   case 0x61:
1153     Vi = Vhp + Vbp + v1;
1154     offset = mixer_offset<3>::value;
1155     break;
1156   case 0x62:
1157     Vi = Vhp + Vbp + v2;
1158     offset = mixer_offset<3>::value;
1159     break;
1160   case 0x63:
1161     Vi = Vhp + Vbp + v2 + v1;
1162     offset = mixer_offset<4>::value;
1163     break;
1164   case 0x64:
1165     Vi = Vhp + Vbp + v3;
1166     offset = mixer_offset<3>::value;
1167     break;
1168   case 0x65:
1169     Vi = Vhp + Vbp + v3 + v1;
1170     offset = mixer_offset<4>::value;
1171     break;
1172   case 0x66:
1173     Vi = Vhp + Vbp + v3 + v2;
1174     offset = mixer_offset<4>::value;
1175     break;
1176   case 0x67:
1177     Vi = Vhp + Vbp + v3 + v2 + v1;
1178     offset = mixer_offset<5>::value;
1179     break;
1180   case 0x68:
1181     Vi = Vhp + Vbp + ve;
1182     offset = mixer_offset<3>::value;
1183     break;
1184   case 0x69:
1185     Vi = Vhp + Vbp + ve + v1;
1186     offset = mixer_offset<4>::value;
1187     break;
1188   case 0x6a:
1189     Vi = Vhp + Vbp + ve + v2;
1190     offset = mixer_offset<4>::value;
1191     break;
1192   case 0x6b:
1193     Vi = Vhp + Vbp + ve + v2 + v1;
1194     offset = mixer_offset<5>::value;
1195     break;
1196   case 0x6c:
1197     Vi = Vhp + Vbp + ve + v3;
1198     offset = mixer_offset<4>::value;
1199     break;
1200   case 0x6d:
1201     Vi = Vhp + Vbp + ve + v3 + v1;
1202     offset = mixer_offset<5>::value;
1203     break;
1204   case 0x6e:
1205     Vi = Vhp + Vbp + ve + v3 + v2;
1206     offset = mixer_offset<5>::value;
1207     break;
1208   case 0x6f:
1209     Vi = Vhp + Vbp + ve + v3 + v2 + v1;
1210     offset = mixer_offset<6>::value;
1211     break;
1212   case 0x70:
1213     Vi = Vhp + Vbp + Vlp;
1214     offset = mixer_offset<3>::value;
1215     break;
1216   case 0x71:
1217     Vi = Vhp + Vbp + Vlp + v1;
1218     offset = mixer_offset<4>::value;
1219     break;
1220   case 0x72:
1221     Vi = Vhp + Vbp + Vlp + v2;
1222     offset = mixer_offset<4>::value;
1223     break;
1224   case 0x73:
1225     Vi = Vhp + Vbp + Vlp + v2 + v1;
1226     offset = mixer_offset<5>::value;
1227     break;
1228   case 0x74:
1229     Vi = Vhp + Vbp + Vlp + v3;
1230     offset = mixer_offset<4>::value;
1231     break;
1232   case 0x75:
1233     Vi = Vhp + Vbp + Vlp + v3 + v1;
1234     offset = mixer_offset<5>::value;
1235     break;
1236   case 0x76:
1237     Vi = Vhp + Vbp + Vlp + v3 + v2;
1238     offset = mixer_offset<5>::value;
1239     break;
1240   case 0x77:
1241     Vi = Vhp + Vbp + Vlp + v3 + v2 + v1;
1242     offset = mixer_offset<6>::value;
1243     break;
1244   case 0x78:
1245     Vi = Vhp + Vbp + Vlp + ve;
1246     offset = mixer_offset<4>::value;
1247     break;
1248   case 0x79:
1249     Vi = Vhp + Vbp + Vlp + ve + v1;
1250     offset = mixer_offset<5>::value;
1251     break;
1252   case 0x7a:
1253     Vi = Vhp + Vbp + Vlp + ve + v2;
1254     offset = mixer_offset<5>::value;
1255     break;
1256   case 0x7b:
1257     Vi = Vhp + Vbp + Vlp + ve + v2 + v1;
1258     offset = mixer_offset<6>::value;
1259     break;
1260   case 0x7c:
1261     Vi = Vhp + Vbp + Vlp + ve + v3;
1262     offset = mixer_offset<5>::value;
1263     break;
1264   case 0x7d:
1265     Vi = Vhp + Vbp + Vlp + ve + v3 + v1;
1266     offset = mixer_offset<6>::value;
1267     break;
1268   case 0x7e:
1269     Vi = Vhp + Vbp + Vlp + ve + v3 + v2;
1270     offset = mixer_offset<6>::value;
1271     break;
1272   case 0x7f:
1273     Vi = Vhp + Vbp + Vlp + ve + v3 + v2 + v1;
1274     offset = mixer_offset<7>::value;
1275     break;
1276   }
1277 
1278   // Sum the inputs in the mixer and run the mixer output through the gain.
1279   if (sid_model == 0) {
1280     return (short)(f.gain[vol][f.mixer[offset + Vi]] - (1 << 15));
1281   }
1282   else {
1283     // FIXME: Temporary code for MOS 8580, should use code above.
1284     /* do hard clipping here, else some tunes manage to overflow this
1285        (eg /MUSICIANS/L/Linus/64_Forever.sid, starting at 0:44) */
1286     int tmp = Vi*(int)vol >> 4;
1287     if (tmp < -32768) tmp = -32768;
1288     if (tmp > 32767) tmp = 32767;
1289     return (short)tmp;  }
1290 }
1291 
1292 
1293 /*
1294 Find output voltage in inverting gain and inverting summer SID op-amp
1295 circuits, using a combination of Newton-Raphson and bisection.
1296 
1297              ---R2--
1298             |       |
1299   vi ---R1-----[A>----- vo
1300             vx
1301 
1302 From Kirchoff's current law it follows that
1303 
1304   IR1f + IR2r = 0
1305 
1306 Substituting the triode mode transistor model K*W/L*(Vgst^2 - Vgdt^2)
1307 for the currents, we get:
1308 
1309   n*((Vddt - vx)^2 - (Vddt - vi)^2) + (Vddt - vx)^2 - (Vddt - vo)^2 = 0
1310 
1311 Our root function f can thus be written as:
1312 
1313   f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - vo)^2 = 0
1314 
1315 We are using the mapping function x = vo - vx -> vx. We thus substitute
1316 for vo = vx + x and get:
1317 
1318   f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - (vx + x))^2 = 0
1319 
1320 Using substitution constants
1321 
1322   a = n + 1
1323   b = Vddt
1324   c = n*(Vddt - vi)^2
1325 
1326 the equations for the root function and its derivative can be written as:
1327 
1328   f = a*(b - vx)^2 - c - (b - (vx + x))^2
1329   df = 2*((b - (vx + x))*(dvx + 1) - a*(b - vx)*dvx)
1330 */
1331 RESID_INLINE
1332 int Filter::solve_gain(opamp_t* opamp, int n, int vi, int& x, model_filter_t& mf)
1333 {
1334   // Note that all variables are translated and scaled in order to fit
1335   // in 16 bits. It is not necessary to explicitly translate the variables here,
1336   // since they are all used in subtractions which cancel out the translation:
1337   // (a - t) - (b - t) = a - b
1338 
1339   // Start off with an estimate of x and a root bracket [ak, bk].
1340   // f is increasing, so that f(ak) < 0 and f(bk) > 0.
1341   int ak = mf.ak, bk = mf.bk;
1342 
1343   int a = n + (1 << 7);              // Scaled by 2^7
1344   int b = mf.kVddt;                  // Scaled by m*2^16
1345   int b_vi = b - vi;                 // Scaled by m*2^16
1346   if (b_vi < 0) b_vi = 0;
1347   int c = n*int(unsigned(b_vi)*unsigned(b_vi) >> 12);    // Scaled by m^2*2^27
1348 
1349   for (;;) {
1350     int xk = x;
1351 
1352     // Calculate f and df.
1353     int vx = opamp[x].vx;      // Scaled by m*2^16
1354     int dvx = opamp[x].dvx;    // Scaled by 2^11
1355 
1356     // f = a*(b - vx)^2 - c - (b - vo)^2
1357     // df = 2*((b - vo)*(dvx + 1) - a*(b - vx)*dvx)
1358     //
1359     int vo = vx + (x << 1) - (1 << 16);
1360     if (vo >= (1 << 16)) {
1361       vo = (1 << 16) - 1;
1362     }
1363     else if (vo < 0) {
1364       vo = 0;
1365     }
1366     int b_vx = b - vx;
1367     if (b_vx < 0) b_vx = 0;
1368     int b_vo = b - vo;
1369     if (b_vo < 0) b_vo = 0;
1370     // The dividend is scaled by m^2*2^27.
1371     int f = a*int(unsigned(b_vx)*unsigned(b_vx) >> 12) - c - int(unsigned(b_vo)*unsigned(b_vo) >> 5);
1372     // The divisor is scaled by m*2^11.
1373     int df = ((b_vo*(dvx + (1 << 11)) >> 1) - (a*(b_vx*dvx >> 8))) >> 14;
1374     // The resulting quotient is thus scaled by m*2^16.
1375 
1376     // Newton-Raphson step: xk1 = xk - f(xk)/f'(xk)
1377     // If f(xk) or f'(xk) are zero then we can't improve further.
1378     if (df) {
1379         x -= f/df;
1380     }
1381     if (unlikely(x == xk)) {
1382       // No further root improvement possible.
1383       return vo;
1384     }
1385 
1386     // Narrow down root bracket.
1387     if (f < 0) {
1388       // f(xk) < 0
1389       ak = xk;
1390     }
1391     else {
1392       // f(xk) > 0
1393       bk = xk;
1394     }
1395 
1396     if (unlikely(x <= ak) || unlikely(x >= bk)) {
1397       // Bisection step (ala Dekker's method).
1398       x = (ak + bk) >> 1;
1399       if (unlikely(x == ak)) {
1400         // No further bisection possible.
1401         return vo;
1402       }
1403     }
1404   }
1405 }
1406 
1407 
1408 /*
1409 Find output voltage in inverting integrator SID op-amp circuits, using a
1410 single fixpoint iteration step.
1411 
1412 A circuit diagram of a MOS 6581 integrator is shown below.
1413 
1414                  ---C---
1415                 |       |
1416   vi -----Rw-------[A>----- vo
1417        |      | vx
1418         --Rs--
1419 
1420 From Kirchoff's current law it follows that
1421 
1422   IRw + IRs + ICr = 0
1423 
1424 Using the formula for current through a capacitor, i = C*dv/dt, we get
1425 
1426   IRw + IRs + C*(vc - vc0)/dt = 0
1427   dt/C*(IRw + IRs) + vc - vc0 = 0
1428   vc = vc0 - n*(IRw(vi,vx) + IRs(vi,vx))
1429 
1430 which may be rewritten as the following iterative fixpoint function:
1431 
1432   vc = vc0 - n*(IRw(vi,g(vc)) + IRs(vi,g(vc)))
1433 
1434 To accurately calculate the currents through Rs and Rw, we need to use
1435 transistor models. Rs has a gate voltage of Vdd = 12V, and can be
1436 assumed to always be in triode mode. For Rw, the situation is rather
1437 more complex, as it turns out that this transistor will operate in
1438 both subthreshold, triode, and saturation modes.
1439 
1440 The Shichman-Hodges transistor model routinely used in textbooks may
1441 be written as follows:
1442 
1443   Ids = 0                          , Vgst < 0               (subthreshold mode)
1444   Ids = K/2*W/L*(2*Vgst - Vds)*Vds , Vgst >= 0, Vds < Vgst  (triode mode)
1445   Ids = K/2*W/L*Vgst^2             , Vgst >= 0, Vds >= Vgst (saturation mode)
1446 
1447   where
1448   K   = u*Cox (conductance)
1449   W/L = ratio between substrate width and length
1450   Vgst = Vg - Vs - Vt (overdrive voltage)
1451 
1452 This transistor model is also called the quadratic model.
1453 
1454 Note that the equation for the triode mode can be reformulated as
1455 independent terms depending on Vgs and Vgd, respectively, by the
1456 following substitution:
1457 
1458   Vds = Vgst - (Vgst - Vds) = Vgst - Vgdt
1459 
1460   Ids = K*W/L*(2*Vgst - Vds)*Vds
1461   = K*W/L*(2*Vgst - (Vgst - Vgdt)*(Vgst - Vgdt)
1462   = K*W/L*(Vgst + Vgdt)*(Vgst - Vgdt)
1463   = K*W/L*(Vgst^2 - Vgdt^2)
1464 
1465 This turns out to be a general equation which covers both the triode
1466 and saturation modes (where the second term is 0 in saturation mode).
1467 The equation is also symmetrical, i.e. it can calculate negative
1468 currents without any change of parameters (since the terms for drain
1469 and source are identical except for the sign).
1470 
1471 FIXME: Subthreshold as function of Vgs, Vgd.
1472   Ids = I0*e^(Vgst/(n*VT))       , Vgst < 0               (subthreshold mode)
1473 
1474 The remaining problem with the textbook model is that the transition
1475 from subthreshold the triode/saturation is not continuous.
1476 
1477 Realizing that the subthreshold and triode/saturation modes may both
1478 be defined by independent (and equal) terms of Vgs and Vds,
1479 respectively, the corresponding terms can be blended into (equal)
1480 continuous functions suitable for table lookup.
1481 
1482 The EKV model (Enz, Krummenacher and Vittoz) essentially performs this
1483 blending using an elegant mathematical formulation:
1484 
1485   Ids = Is*(if - ir)
1486   Is = 2*u*Cox*Ut^2/k*W/L
1487   if = ln^2(1 + e^((k*(Vg - Vt) - Vs)/(2*Ut))
1488   ir = ln^2(1 + e^((k*(Vg - Vt) - Vd)/(2*Ut))
1489 
1490 For our purposes, the EKV model preserves two important properties
1491 discussed above:
1492 
1493 - It consists of two independent terms, which can be represented by
1494   the same lookup table.
1495 - It is symmetrical, i.e. it calculates current in both directions,
1496   facilitating a branch-free implementation.
1497 
1498 Rw in the circuit diagram above is a VCR (voltage controlled resistor),
1499 as shown in the circuit diagram below.
1500 
1501                    Vw
1502 
1503                    |
1504            Vdd     |
1505               |---|
1506              _|_   |
1507            --    --| Vg
1508           |      __|__
1509           |      -----  Rw
1510           |      |   |
1511   vi ------------     -------- vo
1512 
1513 
1514 In order to calculalate the current through the VCR, its gate voltage
1515 must be determined.
1516 
1517 Assuming triode mode and applying Kirchoff's current law, we get the
1518 following equation for Vg:
1519 
1520 u*Cox/2*W/L*((Vddt - Vg)^2 - (Vddt - vi)^2 + (Vddt - Vg)^2 - (Vddt - Vw)^2) = 0
1521 2*(Vddt - Vg)^2 - (Vddt - vi)^2 - (Vddt - Vw)^2 = 0
1522 (Vddt - Vg) = sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
1523 
1524 Vg = Vddt - sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
1525 
1526 */
1527 RESID_INLINE
1528 int Filter::solve_integrate_6581(int dt, int vi, int& vx, int& vc, model_filter_t& mf)
1529 {
1530   // Note that all variables are translated and scaled in order to fit
1531   // in 16 bits. It is not necessary to explicitly translate the variables here,
1532   // since they are all used in subtractions which cancel out the translation:
1533   // (a - t) - (b - t) = a - b
1534 
1535   int kVddt = mf.kVddt;      // Scaled by m*2^16
1536 
1537   // "Snake" voltages for triode mode calculation.
1538   unsigned int Vgst = kVddt - vx;
1539   unsigned int Vgdt = kVddt - vi;
1540   unsigned int Vgdt_2 = Vgdt*Vgdt;
1541 
1542   // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
1543   int n_I_snake = mf.n_snake*(int(Vgst*Vgst - Vgdt_2) >> 15);
1544 
1545   // VCR gate voltage.       // Scaled by m*2^16
1546   // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
1547   int kVg = vcr_kVg[(Vddt_Vw_2 + (Vgdt_2 >> 1)) >> 16];
1548 
1549   // VCR voltages for EKV model table lookup.
1550   int Vgs = kVg - vx;
1551   if (Vgs < 0) Vgs = 0;
1552   int Vgd = kVg - vi;
1553   if (Vgd < 0) Vgd = 0;
1554 
1555   // VCR current, scaled by m*2^15*2^15 = m*2^30
1556   int n_I_vcr = int(unsigned(vcr_n_Ids_term[Vgs] - vcr_n_Ids_term[Vgd]) << 15);
1557 
1558   // Change in capacitor charge.
1559   vc -= (n_I_snake + n_I_vcr)*dt;
1560 
1561 /*
1562   // FIXME: Determine whether this check is necessary.
1563   if (vc < mf.vc_min) {
1564     vc = mf.vc_min;
1565   }
1566   else if (vc > mf.vc_max) {
1567     vc = mf.vc_max;
1568   }
1569 */
1570 
1571   // vx = g(vc)
1572   vx = mf.opamp_rev[(vc >> 15) + (1 << 15)];
1573 
1574   // Return vo.
1575   return vx + (vc >> 14);
1576 }
1577 
1578 #endif // RESID_INLINING || defined(RESID_FILTER_CC)
1579 
1580 } // namespace reSID
1581 
1582 #endif // not RESID_FILTER_H
1583