libsidplayfp  1.8.3
filter.h
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>
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>
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  int vo_N16; // Fixed point scaling for 16 bit op-amp output.
427  int kVddt; // K*(Vdd - Vth)
428  int n_snake;
429  int voice_scale_s14;
430  int voice_DC;
431  int ak;
432  int bk;
433  int vc_min;
434  int vc_max;
435 
436  // Reverse op-amp transfer function.
437  unsigned short opamp_rev[1 << 16];
438  // Lookup tables for gain and summer op-amps in output stage / filter.
439  unsigned short summer[summer_offset<5>::value];
440  unsigned short gain[16][1 << 16];
441  unsigned short mixer[mixer_offset<8>::value];
442  // Cutoff frequency DAC output voltage table. FC is an 11 bit register.
443  unsigned short f0_dac[1 << 11];
444  } model_filter_t;
445 
446  int solve_gain(int* opamp, int n, int vi_t, int& x, model_filter_t& mf);
447  int solve_integrate_6581(int dt, int vi_t, int& x, int& vc, model_filter_t& mf);
448 
449  // VCR - 6581 only.
450  static unsigned short vcr_kVg[1 << 16];
451  static unsigned short vcr_n_Ids_term[1 << 16];
452  // Common parameters.
453  static model_filter_t model_filter[2];
454 
455 friend class SID;
456 };
457 
458 
459 // ----------------------------------------------------------------------------
460 // Inline functions.
461 // The following functions are defined inline because they are called every
462 // time a sample is calculated.
463 // ----------------------------------------------------------------------------
464 
465 #if RESID_INLINING || defined(RESID_FILTER_CC)
466 
467 // ----------------------------------------------------------------------------
468 // SID clocking - 1 cycle.
469 // ----------------------------------------------------------------------------
470 RESID_INLINE
471 void Filter::clock(int voice1, int voice2, int voice3)
472 {
473  model_filter_t& f = model_filter[sid_model];
474 
475  v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
476  v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
477  v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
478 
479  // Sum inputs routed into the filter.
480  int Vi = 0;
481  int offset = 0;
482 
483  switch (sum & 0xf) {
484  case 0x0:
485  Vi = 0;
486  offset = summer_offset<0>::value;
487  break;
488  case 0x1:
489  Vi = v1;
490  offset = summer_offset<1>::value;
491  break;
492  case 0x2:
493  Vi = v2;
494  offset = summer_offset<1>::value;
495  break;
496  case 0x3:
497  Vi = v2 + v1;
498  offset = summer_offset<2>::value;
499  break;
500  case 0x4:
501  Vi = v3;
502  offset = summer_offset<1>::value;
503  break;
504  case 0x5:
505  Vi = v3 + v1;
506  offset = summer_offset<2>::value;
507  break;
508  case 0x6:
509  Vi = v3 + v2;
510  offset = summer_offset<2>::value;
511  break;
512  case 0x7:
513  Vi = v3 + v2 + v1;
514  offset = summer_offset<3>::value;
515  break;
516  case 0x8:
517  Vi = ve;
518  offset = summer_offset<1>::value;
519  break;
520  case 0x9:
521  Vi = ve + v1;
522  offset = summer_offset<2>::value;
523  break;
524  case 0xa:
525  Vi = ve + v2;
526  offset = summer_offset<2>::value;
527  break;
528  case 0xb:
529  Vi = ve + v2 + v1;
530  offset = summer_offset<3>::value;
531  break;
532  case 0xc:
533  Vi = ve + v3;
534  offset = summer_offset<2>::value;
535  break;
536  case 0xd:
537  Vi = ve + v3 + v1;
538  offset = summer_offset<3>::value;
539  break;
540  case 0xe:
541  Vi = ve + v3 + v2;
542  offset = summer_offset<3>::value;
543  break;
544  case 0xf:
545  Vi = ve + v3 + v2 + v1;
546  offset = summer_offset<4>::value;
547  break;
548  }
549 
550  // Calculate filter outputs.
551  if (sid_model == 0) {
552  // MOS 6581.
553  Vlp = solve_integrate_6581(1, Vbp, Vlp_x, Vlp_vc, f);
554  Vbp = solve_integrate_6581(1, Vhp, Vbp_x, Vbp_vc, f);
555  Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
556  }
557  else {
558  // MOS 8580. FIXME: Not yet using op-amp model.
559 
560  // delta_t = 1 is converted to seconds given a 1MHz clock by dividing
561  // with 1 000 000.
562 
563  int dVbp = w0*(Vhp >> 4) >> 16;
564  int dVlp = w0*(Vbp >> 4) >> 16;
565  Vbp -= dVbp;
566  Vlp -= dVlp;
567  Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
568  }
569 }
570 
571 // ----------------------------------------------------------------------------
572 // SID clocking - delta_t cycles.
573 // ----------------------------------------------------------------------------
574 RESID_INLINE
575 void Filter::clock(cycle_count delta_t, int voice1, int voice2, int voice3)
576 {
577  model_filter_t& f = model_filter[sid_model];
578 
579  v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
580  v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
581  v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
582 
583  // Enable filter on/off.
584  // This is not really part of SID, but is useful for testing.
585  // On slow CPUs it may be necessary to bypass the filter to lower the CPU
586  // load.
587  if (unlikely(!enabled)) {
588  return;
589  }
590 
591  // Sum inputs routed into the filter.
592  int Vi = 0;
593  int offset = 0;
594 
595  switch (sum & 0xf) {
596  case 0x0:
597  Vi = 0;
598  offset = summer_offset<0>::value;
599  break;
600  case 0x1:
601  Vi = v1;
602  offset = summer_offset<1>::value;
603  break;
604  case 0x2:
605  Vi = v2;
606  offset = summer_offset<1>::value;
607  break;
608  case 0x3:
609  Vi = v2 + v1;
610  offset = summer_offset<2>::value;
611  break;
612  case 0x4:
613  Vi = v3;
614  offset = summer_offset<1>::value;
615  break;
616  case 0x5:
617  Vi = v3 + v1;
618  offset = summer_offset<2>::value;
619  break;
620  case 0x6:
621  Vi = v3 + v2;
622  offset = summer_offset<2>::value;
623  break;
624  case 0x7:
625  Vi = v3 + v2 + v1;
626  offset = summer_offset<3>::value;
627  break;
628  case 0x8:
629  Vi = ve;
630  offset = summer_offset<1>::value;
631  break;
632  case 0x9:
633  Vi = ve + v1;
634  offset = summer_offset<2>::value;
635  break;
636  case 0xa:
637  Vi = ve + v2;
638  offset = summer_offset<2>::value;
639  break;
640  case 0xb:
641  Vi = ve + v2 + v1;
642  offset = summer_offset<3>::value;
643  break;
644  case 0xc:
645  Vi = ve + v3;
646  offset = summer_offset<2>::value;
647  break;
648  case 0xd:
649  Vi = ve + v3 + v1;
650  offset = summer_offset<3>::value;
651  break;
652  case 0xe:
653  Vi = ve + v3 + v2;
654  offset = summer_offset<3>::value;
655  break;
656  case 0xf:
657  Vi = ve + v3 + v2 + v1;
658  offset = summer_offset<4>::value;
659  break;
660  }
661 
662  // Maximum delta cycles for filter fixpoint iteration to converge
663  // is approximately 3.
664  cycle_count delta_t_flt = 3;
665 
666  if (sid_model == 0) {
667  // MOS 6581.
668  while (delta_t) {
669  if (unlikely(delta_t < delta_t_flt)) {
670  delta_t_flt = delta_t;
671  }
672 
673  // Calculate filter outputs.
674  Vlp = solve_integrate_6581(delta_t_flt, Vbp, Vlp_x, Vlp_vc, f);
675  Vbp = solve_integrate_6581(delta_t_flt, Vhp, Vbp_x, Vbp_vc, f);
676  Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
677 
678  delta_t -= delta_t_flt;
679  }
680  }
681  else {
682  // MOS 8580. FIXME: Not yet using op-amp model.
683  while (delta_t) {
684  if (delta_t < delta_t_flt) {
685  delta_t_flt = delta_t;
686  }
687 
688  // delta_t is converted to seconds given a 1MHz clock by dividing
689  // with 1 000 000. This is done in two operations to avoid integer
690  // multiplication overflow.
691 
692  // Calculate filter outputs.
693  int w0_delta_t = w0*delta_t_flt >> 2;
694 
695  int dVbp = w0_delta_t*(Vhp >> 4) >> 14;
696  int dVlp = w0_delta_t*(Vbp >> 4) >> 14;
697  Vbp -= dVbp;
698  Vlp -= dVlp;
699  Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
700 
701  delta_t -= delta_t_flt;
702  }
703  }
704 }
705 
706 
707 // ----------------------------------------------------------------------------
708 // SID audio input (16 bits).
709 // ----------------------------------------------------------------------------
710 RESID_INLINE
711 void Filter::input(short sample)
712 {
713  // Scale to three times the peak-to-peak for one voice and add the op-amp
714  // "zero" DC level.
715  // NB! Adding the op-amp "zero" DC level is a (wildly inaccurate)
716  // approximation of feeding the input through an AC coupling capacitor.
717  // This could be implemented as a separate filter circuit, however the
718  // primary use of the emulator is not to process external signals.
719  // The upside is that the MOS8580 "digi boost" works without a separate (DC)
720  // input interface.
721  // Note that the input is 16 bits, compared to the 20 bit voice output.
722  model_filter_t& f = model_filter[sid_model];
723  ve = (sample*f.voice_scale_s14*3 >> 14) + f.mixer[0];
724 }
725 
726 
727 // ----------------------------------------------------------------------------
728 // SID audio output (16 bits).
729 // ----------------------------------------------------------------------------
730 RESID_INLINE
731 short Filter::output()
732 {
733  model_filter_t& f = model_filter[sid_model];
734 
735  // Writing the switch below manually would be tedious and error-prone;
736  // it is rather generated by the following Perl program:
737 
738  /*
739 my @i = qw(v1 v2 v3 ve Vlp Vbp Vhp);
740 for my $mix (0..2**@i-1) {
741  print sprintf(" case 0x%02x:\n", $mix);
742  my @sum;
743  for (@i) {
744  unshift(@sum, $_) if $mix & 0x01;
745  $mix >>= 1;
746  }
747  my $sum = join(" + ", @sum) || "0";
748  print " Vi = $sum;\n";
749  print " offset = mixer_offset<" . @sum . ">::value;\n";
750  print " break;\n";
751 }
752  */
753 
754  // Sum inputs routed into the mixer.
755  int Vi = 0;
756  int offset = 0;
757 
758  switch (mix & 0x7f) {
759  case 0x00:
760  Vi = 0;
761  offset = mixer_offset<0>::value;
762  break;
763  case 0x01:
764  Vi = v1;
765  offset = mixer_offset<1>::value;
766  break;
767  case 0x02:
768  Vi = v2;
769  offset = mixer_offset<1>::value;
770  break;
771  case 0x03:
772  Vi = v2 + v1;
773  offset = mixer_offset<2>::value;
774  break;
775  case 0x04:
776  Vi = v3;
777  offset = mixer_offset<1>::value;
778  break;
779  case 0x05:
780  Vi = v3 + v1;
781  offset = mixer_offset<2>::value;
782  break;
783  case 0x06:
784  Vi = v3 + v2;
785  offset = mixer_offset<2>::value;
786  break;
787  case 0x07:
788  Vi = v3 + v2 + v1;
789  offset = mixer_offset<3>::value;
790  break;
791  case 0x08:
792  Vi = ve;
793  offset = mixer_offset<1>::value;
794  break;
795  case 0x09:
796  Vi = ve + v1;
797  offset = mixer_offset<2>::value;
798  break;
799  case 0x0a:
800  Vi = ve + v2;
801  offset = mixer_offset<2>::value;
802  break;
803  case 0x0b:
804  Vi = ve + v2 + v1;
805  offset = mixer_offset<3>::value;
806  break;
807  case 0x0c:
808  Vi = ve + v3;
809  offset = mixer_offset<2>::value;
810  break;
811  case 0x0d:
812  Vi = ve + v3 + v1;
813  offset = mixer_offset<3>::value;
814  break;
815  case 0x0e:
816  Vi = ve + v3 + v2;
817  offset = mixer_offset<3>::value;
818  break;
819  case 0x0f:
820  Vi = ve + v3 + v2 + v1;
821  offset = mixer_offset<4>::value;
822  break;
823  case 0x10:
824  Vi = Vlp;
825  offset = mixer_offset<1>::value;
826  break;
827  case 0x11:
828  Vi = Vlp + v1;
829  offset = mixer_offset<2>::value;
830  break;
831  case 0x12:
832  Vi = Vlp + v2;
833  offset = mixer_offset<2>::value;
834  break;
835  case 0x13:
836  Vi = Vlp + v2 + v1;
837  offset = mixer_offset<3>::value;
838  break;
839  case 0x14:
840  Vi = Vlp + v3;
841  offset = mixer_offset<2>::value;
842  break;
843  case 0x15:
844  Vi = Vlp + v3 + v1;
845  offset = mixer_offset<3>::value;
846  break;
847  case 0x16:
848  Vi = Vlp + v3 + v2;
849  offset = mixer_offset<3>::value;
850  break;
851  case 0x17:
852  Vi = Vlp + v3 + v2 + v1;
853  offset = mixer_offset<4>::value;
854  break;
855  case 0x18:
856  Vi = Vlp + ve;
857  offset = mixer_offset<2>::value;
858  break;
859  case 0x19:
860  Vi = Vlp + ve + v1;
861  offset = mixer_offset<3>::value;
862  break;
863  case 0x1a:
864  Vi = Vlp + ve + v2;
865  offset = mixer_offset<3>::value;
866  break;
867  case 0x1b:
868  Vi = Vlp + ve + v2 + v1;
869  offset = mixer_offset<4>::value;
870  break;
871  case 0x1c:
872  Vi = Vlp + ve + v3;
873  offset = mixer_offset<3>::value;
874  break;
875  case 0x1d:
876  Vi = Vlp + ve + v3 + v1;
877  offset = mixer_offset<4>::value;
878  break;
879  case 0x1e:
880  Vi = Vlp + ve + v3 + v2;
881  offset = mixer_offset<4>::value;
882  break;
883  case 0x1f:
884  Vi = Vlp + ve + v3 + v2 + v1;
885  offset = mixer_offset<5>::value;
886  break;
887  case 0x20:
888  Vi = Vbp;
889  offset = mixer_offset<1>::value;
890  break;
891  case 0x21:
892  Vi = Vbp + v1;
893  offset = mixer_offset<2>::value;
894  break;
895  case 0x22:
896  Vi = Vbp + v2;
897  offset = mixer_offset<2>::value;
898  break;
899  case 0x23:
900  Vi = Vbp + v2 + v1;
901  offset = mixer_offset<3>::value;
902  break;
903  case 0x24:
904  Vi = Vbp + v3;
905  offset = mixer_offset<2>::value;
906  break;
907  case 0x25:
908  Vi = Vbp + v3 + v1;
909  offset = mixer_offset<3>::value;
910  break;
911  case 0x26:
912  Vi = Vbp + v3 + v2;
913  offset = mixer_offset<3>::value;
914  break;
915  case 0x27:
916  Vi = Vbp + v3 + v2 + v1;
917  offset = mixer_offset<4>::value;
918  break;
919  case 0x28:
920  Vi = Vbp + ve;
921  offset = mixer_offset<2>::value;
922  break;
923  case 0x29:
924  Vi = Vbp + ve + v1;
925  offset = mixer_offset<3>::value;
926  break;
927  case 0x2a:
928  Vi = Vbp + ve + v2;
929  offset = mixer_offset<3>::value;
930  break;
931  case 0x2b:
932  Vi = Vbp + ve + v2 + v1;
933  offset = mixer_offset<4>::value;
934  break;
935  case 0x2c:
936  Vi = Vbp + ve + v3;
937  offset = mixer_offset<3>::value;
938  break;
939  case 0x2d:
940  Vi = Vbp + ve + v3 + v1;
941  offset = mixer_offset<4>::value;
942  break;
943  case 0x2e:
944  Vi = Vbp + ve + v3 + v2;
945  offset = mixer_offset<4>::value;
946  break;
947  case 0x2f:
948  Vi = Vbp + ve + v3 + v2 + v1;
949  offset = mixer_offset<5>::value;
950  break;
951  case 0x30:
952  Vi = Vbp + Vlp;
953  offset = mixer_offset<2>::value;
954  break;
955  case 0x31:
956  Vi = Vbp + Vlp + v1;
957  offset = mixer_offset<3>::value;
958  break;
959  case 0x32:
960  Vi = Vbp + Vlp + v2;
961  offset = mixer_offset<3>::value;
962  break;
963  case 0x33:
964  Vi = Vbp + Vlp + v2 + v1;
965  offset = mixer_offset<4>::value;
966  break;
967  case 0x34:
968  Vi = Vbp + Vlp + v3;
969  offset = mixer_offset<3>::value;
970  break;
971  case 0x35:
972  Vi = Vbp + Vlp + v3 + v1;
973  offset = mixer_offset<4>::value;
974  break;
975  case 0x36:
976  Vi = Vbp + Vlp + v3 + v2;
977  offset = mixer_offset<4>::value;
978  break;
979  case 0x37:
980  Vi = Vbp + Vlp + v3 + v2 + v1;
981  offset = mixer_offset<5>::value;
982  break;
983  case 0x38:
984  Vi = Vbp + Vlp + ve;
985  offset = mixer_offset<3>::value;
986  break;
987  case 0x39:
988  Vi = Vbp + Vlp + ve + v1;
989  offset = mixer_offset<4>::value;
990  break;
991  case 0x3a:
992  Vi = Vbp + Vlp + ve + v2;
993  offset = mixer_offset<4>::value;
994  break;
995  case 0x3b:
996  Vi = Vbp + Vlp + ve + v2 + v1;
997  offset = mixer_offset<5>::value;
998  break;
999  case 0x3c:
1000  Vi = Vbp + Vlp + ve + v3;
1001  offset = mixer_offset<4>::value;
1002  break;
1003  case 0x3d:
1004  Vi = Vbp + Vlp + ve + v3 + v1;
1005  offset = mixer_offset<5>::value;
1006  break;
1007  case 0x3e:
1008  Vi = Vbp + Vlp + ve + v3 + v2;
1009  offset = mixer_offset<5>::value;
1010  break;
1011  case 0x3f:
1012  Vi = Vbp + Vlp + ve + v3 + v2 + v1;
1013  offset = mixer_offset<6>::value;
1014  break;
1015  case 0x40:
1016  Vi = Vhp;
1017  offset = mixer_offset<1>::value;
1018  break;
1019  case 0x41:
1020  Vi = Vhp + v1;
1021  offset = mixer_offset<2>::value;
1022  break;
1023  case 0x42:
1024  Vi = Vhp + v2;
1025  offset = mixer_offset<2>::value;
1026  break;
1027  case 0x43:
1028  Vi = Vhp + v2 + v1;
1029  offset = mixer_offset<3>::value;
1030  break;
1031  case 0x44:
1032  Vi = Vhp + v3;
1033  offset = mixer_offset<2>::value;
1034  break;
1035  case 0x45:
1036  Vi = Vhp + v3 + v1;
1037  offset = mixer_offset<3>::value;
1038  break;
1039  case 0x46:
1040  Vi = Vhp + v3 + v2;
1041  offset = mixer_offset<3>::value;
1042  break;
1043  case 0x47:
1044  Vi = Vhp + v3 + v2 + v1;
1045  offset = mixer_offset<4>::value;
1046  break;
1047  case 0x48:
1048  Vi = Vhp + ve;
1049  offset = mixer_offset<2>::value;
1050  break;
1051  case 0x49:
1052  Vi = Vhp + ve + v1;
1053  offset = mixer_offset<3>::value;
1054  break;
1055  case 0x4a:
1056  Vi = Vhp + ve + v2;
1057  offset = mixer_offset<3>::value;
1058  break;
1059  case 0x4b:
1060  Vi = Vhp + ve + v2 + v1;
1061  offset = mixer_offset<4>::value;
1062  break;
1063  case 0x4c:
1064  Vi = Vhp + ve + v3;
1065  offset = mixer_offset<3>::value;
1066  break;
1067  case 0x4d:
1068  Vi = Vhp + ve + v3 + v1;
1069  offset = mixer_offset<4>::value;
1070  break;
1071  case 0x4e:
1072  Vi = Vhp + ve + v3 + v2;
1073  offset = mixer_offset<4>::value;
1074  break;
1075  case 0x4f:
1076  Vi = Vhp + ve + v3 + v2 + v1;
1077  offset = mixer_offset<5>::value;
1078  break;
1079  case 0x50:
1080  Vi = Vhp + Vlp;
1081  offset = mixer_offset<2>::value;
1082  break;
1083  case 0x51:
1084  Vi = Vhp + Vlp + v1;
1085  offset = mixer_offset<3>::value;
1086  break;
1087  case 0x52:
1088  Vi = Vhp + Vlp + v2;
1089  offset = mixer_offset<3>::value;
1090  break;
1091  case 0x53:
1092  Vi = Vhp + Vlp + v2 + v1;
1093  offset = mixer_offset<4>::value;
1094  break;
1095  case 0x54:
1096  Vi = Vhp + Vlp + v3;
1097  offset = mixer_offset<3>::value;
1098  break;
1099  case 0x55:
1100  Vi = Vhp + Vlp + v3 + v1;
1101  offset = mixer_offset<4>::value;
1102  break;
1103  case 0x56:
1104  Vi = Vhp + Vlp + v3 + v2;
1105  offset = mixer_offset<4>::value;
1106  break;
1107  case 0x57:
1108  Vi = Vhp + Vlp + v3 + v2 + v1;
1109  offset = mixer_offset<5>::value;
1110  break;
1111  case 0x58:
1112  Vi = Vhp + Vlp + ve;
1113  offset = mixer_offset<3>::value;
1114  break;
1115  case 0x59:
1116  Vi = Vhp + Vlp + ve + v1;
1117  offset = mixer_offset<4>::value;
1118  break;
1119  case 0x5a:
1120  Vi = Vhp + Vlp + ve + v2;
1121  offset = mixer_offset<4>::value;
1122  break;
1123  case 0x5b:
1124  Vi = Vhp + Vlp + ve + v2 + v1;
1125  offset = mixer_offset<5>::value;
1126  break;
1127  case 0x5c:
1128  Vi = Vhp + Vlp + ve + v3;
1129  offset = mixer_offset<4>::value;
1130  break;
1131  case 0x5d:
1132  Vi = Vhp + Vlp + ve + v3 + v1;
1133  offset = mixer_offset<5>::value;
1134  break;
1135  case 0x5e:
1136  Vi = Vhp + Vlp + ve + v3 + v2;
1137  offset = mixer_offset<5>::value;
1138  break;
1139  case 0x5f:
1140  Vi = Vhp + Vlp + ve + v3 + v2 + v1;
1141  offset = mixer_offset<6>::value;
1142  break;
1143  case 0x60:
1144  Vi = Vhp + Vbp;
1145  offset = mixer_offset<2>::value;
1146  break;
1147  case 0x61:
1148  Vi = Vhp + Vbp + v1;
1149  offset = mixer_offset<3>::value;
1150  break;
1151  case 0x62:
1152  Vi = Vhp + Vbp + v2;
1153  offset = mixer_offset<3>::value;
1154  break;
1155  case 0x63:
1156  Vi = Vhp + Vbp + v2 + v1;
1157  offset = mixer_offset<4>::value;
1158  break;
1159  case 0x64:
1160  Vi = Vhp + Vbp + v3;
1161  offset = mixer_offset<3>::value;
1162  break;
1163  case 0x65:
1164  Vi = Vhp + Vbp + v3 + v1;
1165  offset = mixer_offset<4>::value;
1166  break;
1167  case 0x66:
1168  Vi = Vhp + Vbp + v3 + v2;
1169  offset = mixer_offset<4>::value;
1170  break;
1171  case 0x67:
1172  Vi = Vhp + Vbp + v3 + v2 + v1;
1173  offset = mixer_offset<5>::value;
1174  break;
1175  case 0x68:
1176  Vi = Vhp + Vbp + ve;
1177  offset = mixer_offset<3>::value;
1178  break;
1179  case 0x69:
1180  Vi = Vhp + Vbp + ve + v1;
1181  offset = mixer_offset<4>::value;
1182  break;
1183  case 0x6a:
1184  Vi = Vhp + Vbp + ve + v2;
1185  offset = mixer_offset<4>::value;
1186  break;
1187  case 0x6b:
1188  Vi = Vhp + Vbp + ve + v2 + v1;
1189  offset = mixer_offset<5>::value;
1190  break;
1191  case 0x6c:
1192  Vi = Vhp + Vbp + ve + v3;
1193  offset = mixer_offset<4>::value;
1194  break;
1195  case 0x6d:
1196  Vi = Vhp + Vbp + ve + v3 + v1;
1197  offset = mixer_offset<5>::value;
1198  break;
1199  case 0x6e:
1200  Vi = Vhp + Vbp + ve + v3 + v2;
1201  offset = mixer_offset<5>::value;
1202  break;
1203  case 0x6f:
1204  Vi = Vhp + Vbp + ve + v3 + v2 + v1;
1205  offset = mixer_offset<6>::value;
1206  break;
1207  case 0x70:
1208  Vi = Vhp + Vbp + Vlp;
1209  offset = mixer_offset<3>::value;
1210  break;
1211  case 0x71:
1212  Vi = Vhp + Vbp + Vlp + v1;
1213  offset = mixer_offset<4>::value;
1214  break;
1215  case 0x72:
1216  Vi = Vhp + Vbp + Vlp + v2;
1217  offset = mixer_offset<4>::value;
1218  break;
1219  case 0x73:
1220  Vi = Vhp + Vbp + Vlp + v2 + v1;
1221  offset = mixer_offset<5>::value;
1222  break;
1223  case 0x74:
1224  Vi = Vhp + Vbp + Vlp + v3;
1225  offset = mixer_offset<4>::value;
1226  break;
1227  case 0x75:
1228  Vi = Vhp + Vbp + Vlp + v3 + v1;
1229  offset = mixer_offset<5>::value;
1230  break;
1231  case 0x76:
1232  Vi = Vhp + Vbp + Vlp + v3 + v2;
1233  offset = mixer_offset<5>::value;
1234  break;
1235  case 0x77:
1236  Vi = Vhp + Vbp + Vlp + v3 + v2 + v1;
1237  offset = mixer_offset<6>::value;
1238  break;
1239  case 0x78:
1240  Vi = Vhp + Vbp + Vlp + ve;
1241  offset = mixer_offset<4>::value;
1242  break;
1243  case 0x79:
1244  Vi = Vhp + Vbp + Vlp + ve + v1;
1245  offset = mixer_offset<5>::value;
1246  break;
1247  case 0x7a:
1248  Vi = Vhp + Vbp + Vlp + ve + v2;
1249  offset = mixer_offset<5>::value;
1250  break;
1251  case 0x7b:
1252  Vi = Vhp + Vbp + Vlp + ve + v2 + v1;
1253  offset = mixer_offset<6>::value;
1254  break;
1255  case 0x7c:
1256  Vi = Vhp + Vbp + Vlp + ve + v3;
1257  offset = mixer_offset<5>::value;
1258  break;
1259  case 0x7d:
1260  Vi = Vhp + Vbp + Vlp + ve + v3 + v1;
1261  offset = mixer_offset<6>::value;
1262  break;
1263  case 0x7e:
1264  Vi = Vhp + Vbp + Vlp + ve + v3 + v2;
1265  offset = mixer_offset<6>::value;
1266  break;
1267  case 0x7f:
1268  Vi = Vhp + Vbp + Vlp + ve + v3 + v2 + v1;
1269  offset = mixer_offset<7>::value;
1270  break;
1271  }
1272 
1273  // Sum the inputs in the mixer and run the mixer output through the gain.
1274  if (sid_model == 0) {
1275  return (short)(f.gain[vol][f.mixer[offset + Vi]] - (1 << 15));
1276  }
1277  else {
1278  // FIXME: Temporary code for MOS 8580, should use code above.
1279  /* do hard clipping here, else some tunes manage to overflow this
1280  (eg /MUSICIANS/L/Linus/64_Forever.sid, starting at 0:44) */
1281  int tmp = Vi*(int)vol >> 4;
1282  if (tmp < -32768) tmp = -32768;
1283  if (tmp > 32767) tmp = 32767;
1284  return (short)tmp; }
1285 }
1286 
1287 
1288 /*
1289 Find output voltage in inverting gain and inverting summer SID op-amp
1290 circuits, using a combination of Newton-Raphson and bisection.
1291 
1292  ---R2--
1293  | |
1294  vi ---R1-----[A>----- vo
1295  vx
1296 
1297 From Kirchoff's current law it follows that
1298 
1299  IR1f + IR2r = 0
1300 
1301 Substituting the triode mode transistor model K*W/L*(Vgst^2 - Vgdt^2)
1302 for the currents, we get:
1303 
1304  n*((Vddt - vx)^2 - (Vddt - vi)^2) + (Vddt - vx)^2 - (Vddt - vo)^2 = 0
1305 
1306 Our root function f can thus be written as:
1307 
1308  f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - vo)^2 = 0
1309 
1310 We are using the mapping function x = vo - vx -> vx. We thus substitute
1311 for vo = vx + x and get:
1312 
1313  f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - (vx + x))^2 = 0
1314 
1315 Using substitution constants
1316 
1317  a = n + 1
1318  b = Vddt
1319  c = n*(Vddt - vi)^2
1320 
1321 the equations for the root function and its derivative can be written as:
1322 
1323  f = a*(b - vx)^2 - c - (b - (vx + x))^2
1324  df = 2*((b - (vx + x))*(dvx + 1) - a*(b - vx)*dvx)
1325 */
1326 RESID_INLINE
1327 int Filter::solve_gain(int* opamp, int n, int vi, int& x, model_filter_t& mf)
1328 {
1329  // Note that all variables are translated and scaled in order to fit
1330  // in 16 bits. It is not necessary to explicitly translate the variables here,
1331  // since they are all used in subtractions which cancel out the translation:
1332  // (a - t) - (b - t) = a - b
1333 
1334  // Start off with an estimate of x and a root bracket [ak, bk].
1335  // f is increasing, so that f(ak) < 0 and f(bk) > 0.
1336  int ak = mf.ak, bk = mf.bk;
1337 
1338  int a = n + (1 << 7); // Scaled by 2^7
1339  int b = mf.kVddt; // Scaled by m*2^16
1340  int b_vi = b - vi; // Scaled by m*2^16
1341  if (b_vi < 0) b_vi = 0;
1342  int c = n*int(unsigned(b_vi)*unsigned(b_vi) >> 12); // Scaled by m^2*2^27
1343 
1344  for (;;) {
1345  int xk = x;
1346 
1347  // Calculate f and df.
1348  int vx_dvx = opamp[x];
1349  int vx = vx_dvx & 0xffff; // Scaled by m*2^16
1350  int dvx = vx_dvx >> 16; // Scaled by 2^11
1351 
1352  // f = a*(b - vx)^2 - c - (b - vo)^2
1353  // df = 2*((b - vo)*(dvx + 1) - a*(b - vx)*dvx)
1354  //
1355  int vo = vx + (x << 1) - (1 << 16);
1356  if (vo >= (1 << 16)) {
1357  vo = (1 << 16) - 1;
1358  }
1359  else if (vo < 0) {
1360  vo = 0;
1361  }
1362  int b_vx = b - vx;
1363  if (b_vx < 0) b_vx = 0;
1364  int b_vo = b - vo;
1365  if (b_vo < 0) b_vo = 0;
1366  // The dividend is scaled by m^2*2^27.
1367  int f = a*int(unsigned(b_vx)*unsigned(b_vx) >> 12) - c - int(unsigned(b_vo)*unsigned(b_vo) >> 5);
1368  // The divisor is scaled by m*2^11.
1369  int df = (b_vo*(dvx + (1 << 11)) - a*(b_vx*dvx >> 7)) >> 15;
1370  // The resulting quotient is thus scaled by m*2^16.
1371 
1372  // Newton-Raphson step: xk1 = xk - f(xk)/f'(xk)
1373  x -= f/df;
1374  if (unlikely(x == xk)) {
1375  // No further root improvement possible.
1376  return vo;
1377  }
1378 
1379  // Narrow down root bracket.
1380  if (f < 0) {
1381  // f(xk) < 0
1382  ak = xk;
1383  }
1384  else {
1385  // f(xk) > 0
1386  bk = xk;
1387  }
1388 
1389  if (unlikely(x <= ak) || unlikely(x >= bk)) {
1390  // Bisection step (ala Dekker's method).
1391  x = (ak + bk) >> 1;
1392  if (unlikely(x == ak)) {
1393  // No further bisection possible.
1394  return vo;
1395  }
1396  }
1397  }
1398 }
1399 
1400 
1401 /*
1402 Find output voltage in inverting integrator SID op-amp circuits, using a
1403 single fixpoint iteration step.
1404 
1405 A circuit diagram of a MOS 6581 integrator is shown below.
1406 
1407  ---C---
1408  | |
1409  vi -----Rw-------[A>----- vo
1410  | | vx
1411  --Rs--
1412 
1413 From Kirchoff's current law it follows that
1414 
1415  IRw + IRs + ICr = 0
1416 
1417 Using the formula for current through a capacitor, i = C*dv/dt, we get
1418 
1419  IRw + IRs + C*(vc - vc0)/dt = 0
1420  dt/C*(IRw + IRs) + vc - vc0 = 0
1421  vc = vc0 - n*(IRw(vi,vx) + IRs(vi,vx))
1422 
1423 which may be rewritten as the following iterative fixpoint function:
1424 
1425  vc = vc0 - n*(IRw(vi,g(vc)) + IRs(vi,g(vc)))
1426 
1427 To accurately calculate the currents through Rs and Rw, we need to use
1428 transistor models. Rs has a gate voltage of Vdd = 12V, and can be
1429 assumed to always be in triode mode. For Rw, the situation is rather
1430 more complex, as it turns out that this transistor will operate in
1431 both subthreshold, triode, and saturation modes.
1432 
1433 The Shichman-Hodges transistor model routinely used in textbooks may
1434 be written as follows:
1435 
1436  Ids = 0 , Vgst < 0 (subthreshold mode)
1437  Ids = K/2*W/L*(2*Vgst - Vds)*Vds , Vgst >= 0, Vds < Vgst (triode mode)
1438  Ids = K/2*W/L*Vgst^2 , Vgst >= 0, Vds >= Vgst (saturation mode)
1439 
1440  where
1441  K = u*Cox (conductance)
1442  W/L = ratio between substrate width and length
1443  Vgst = Vg - Vs - Vt (overdrive voltage)
1444 
1445 This transistor model is also called the quadratic model.
1446 
1447 Note that the equation for the triode mode can be reformulated as
1448 independent terms depending on Vgs and Vgd, respectively, by the
1449 following substitution:
1450 
1451  Vds = Vgst - (Vgst - Vds) = Vgst - Vgdt
1452 
1453  Ids = K*W/L*(2*Vgst - Vds)*Vds
1454  = K*W/L*(2*Vgst - (Vgst - Vgdt)*(Vgst - Vgdt)
1455  = K*W/L*(Vgst + Vgdt)*(Vgst - Vgdt)
1456  = K*W/L*(Vgst^2 - Vgdt^2)
1457 
1458 This turns out to be a general equation which covers both the triode
1459 and saturation modes (where the second term is 0 in saturation mode).
1460 The equation is also symmetrical, i.e. it can calculate negative
1461 currents without any change of parameters (since the terms for drain
1462 and source are identical except for the sign).
1463 
1464 FIXME: Subthreshold as function of Vgs, Vgd.
1465  Ids = I0*e^(Vgst/(n*VT)) , Vgst < 0 (subthreshold mode)
1466 
1467 The remaining problem with the textbook model is that the transition
1468 from subthreshold the triode/saturation is not continuous.
1469 
1470 Realizing that the subthreshold and triode/saturation modes may both
1471 be defined by independent (and equal) terms of Vgs and Vds,
1472 respectively, the corresponding terms can be blended into (equal)
1473 continuous functions suitable for table lookup.
1474 
1475 The EKV model (Enz, Krummenacher and Vittoz) essentially performs this
1476 blending using an elegant mathematical formulation:
1477 
1478  Ids = Is*(if - ir)
1479  Is = 2*u*Cox*Ut^2/k*W/L
1480  if = ln^2(1 + e^((k*(Vg - Vt) - Vs)/(2*Ut))
1481  ir = ln^2(1 + e^((k*(Vg - Vt) - Vd)/(2*Ut))
1482 
1483 For our purposes, the EKV model preserves two important properties
1484 discussed above:
1485 
1486 - It consists of two independent terms, which can be represented by
1487  the same lookup table.
1488 - It is symmetrical, i.e. it calculates current in both directions,
1489  facilitating a branch-free implementation.
1490 
1491 Rw in the circuit diagram above is a VCR (voltage controlled resistor),
1492 as shown in the circuit diagram below.
1493 
1494  Vw
1495 
1496  |
1497  Vdd |
1498  |---|
1499  _|_ |
1500  -- --| Vg
1501  | __|__
1502  | ----- Rw
1503  | | |
1504  vi ------------ -------- vo
1505 
1506 
1507 In order to calculalate the current through the VCR, its gate voltage
1508 must be determined.
1509 
1510 Assuming triode mode and applying Kirchoff's current law, we get the
1511 following equation for Vg:
1512 
1513 u*Cox/2*W/L*((Vddt - Vg)^2 - (Vddt - vi)^2 + (Vddt - Vg)^2 - (Vddt - Vw)^2) = 0
1514 2*(Vddt - Vg)^2 - (Vddt - vi)^2 - (Vddt - Vw)^2 = 0
1515 (Vddt - Vg) = sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
1516 
1517 Vg = Vddt - sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
1518 
1519 */
1520 RESID_INLINE
1521 int Filter::solve_integrate_6581(int dt, int vi, int& vx, int& vc,
1522  model_filter_t& mf)
1523 {
1524  // Note that all variables are translated and scaled in order to fit
1525  // in 16 bits. It is not necessary to explicitly translate the variables here,
1526  // since they are all used in subtractions which cancel out the translation:
1527  // (a - t) - (b - t) = a - b
1528 
1529  int kVddt = mf.kVddt; // Scaled by m*2^16
1530 
1531  // "Snake" voltages for triode mode calculation.
1532  unsigned int Vgst = kVddt - vx;
1533  unsigned int Vgdt = kVddt - vi;
1534  unsigned int Vgdt_2 = Vgdt*Vgdt;
1535 
1536  // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
1537  int n_I_snake = mf.n_snake*(int(Vgst*Vgst - Vgdt_2) >> 15);
1538 
1539  // VCR gate voltage. // Scaled by m*2^16
1540  // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
1541  int kVg = vcr_kVg[(Vddt_Vw_2 + (Vgdt_2 >> 1)) >> 16];
1542 
1543  // VCR voltages for EKV model table lookup.
1544  int Vgs = kVg - vx;
1545  if (Vgs < 0) Vgs = 0;
1546  int Vgd = kVg - vi;
1547  if (Vgd < 0) Vgd = 0;
1548 
1549  // VCR current, scaled by m*2^15*2^15 = m*2^30
1550  int n_I_vcr = (vcr_n_Ids_term[Vgs] - vcr_n_Ids_term[Vgd]) << 15;
1551 
1552  // Change in capacitor charge.
1553  vc -= (n_I_snake + n_I_vcr)*dt;
1554 
1555 /*
1556  // FIXME: Determine whether this check is necessary.
1557  if (vc < mf.vc_min) {
1558  vc = mf.vc_min;
1559  }
1560  else if (vc > mf.vc_max) {
1561  vc = mf.vc_max;
1562  }
1563 */
1564 
1565  // vx = g(vc)
1566  vx = mf.opamp_rev[(vc >> 15) + (1 << 15)];
1567 
1568  // Return vo.
1569  return vx + (vc >> 14);
1570 }
1571 
1572 #endif // RESID_INLINING || defined(RESID_FILTER_CC)
1573 
1574 } // namespace reSID
1575 
1576 #endif // not RESID_FILTER_H
Definition: dac.cc:45
Definition: filter.h:328
Definition: sid.h:34
Definition: filter.h:315
Definition: filter.h:346
Definition: filter.h:425