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