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