libsidplayfp  1.0.3
filter.h
00001 //  ---------------------------------------------------------------------------
00002 //  This file is part of reSID, a MOS6581 SID emulator engine.
00003 //  Copyright (C) 2010  Dag Lem <resid@nimrod.no>
00004 //
00005 //  This program is free software; you can redistribute it and/or modify
00006 //  it under the terms of the GNU General Public License as published by
00007 //  the Free Software Foundation; either version 2 of the License, or
00008 //  (at your option) any later version.
00009 //
00010 //  This program is distributed in the hope that it will be useful,
00011 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00012 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013 //  GNU General Public License for more details.
00014 //
00015 //  You should have received a copy of the GNU General Public License
00016 //  along with this program; if not, write to the Free Software
00017 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00018 //  ---------------------------------------------------------------------------
00019 
00020 #ifndef RESID_FILTER_H
00021 #define RESID_FILTER_H
00022 
00023 #include "resid-config.h"
00024 
00025 namespace reSID
00026 {
00027 
00028 // ----------------------------------------------------------------------------
00029 // The SID filter is modeled with a two-integrator-loop biquadratic filter,
00030 // which has been confirmed by Bob Yannes to be the actual circuit used in
00031 // the SID chip.
00032 //
00033 // Measurements show that excellent emulation of the SID filter is achieved,
00034 // except when high resonance is combined with high sustain levels.
00035 // In this case the SID op-amps are performing less than ideally and are
00036 // causing some peculiar behavior of the SID filter. This however seems to
00037 // have more effect on the overall amplitude than on the color of the sound.
00038 //
00039 // The theory for the filter circuit can be found in "Microelectric Circuits"
00040 // by Adel S. Sedra and Kenneth C. Smith.
00041 // The circuit is modeled based on the explanation found there except that
00042 // an additional inverter is used in the feedback from the bandpass output,
00043 // allowing the summer op-amp to operate in single-ended mode. This yields
00044 // filter outputs with levels independent of Q, which corresponds with the
00045 // results obtained from a real SID.
00046 //
00047 // We have been able to model the summer and the two integrators of the circuit
00048 // to form components of an IIR filter.
00049 // Vhp is the output of the summer, Vbp is the output of the first integrator,
00050 // and Vlp is the output of the second integrator in the filter circuit.
00051 //
00052 // According to Bob Yannes, the active stages of the SID filter are not really
00053 // op-amps. Rather, simple NMOS inverters are used. By biasing an inverter
00054 // into its region of quasi-linear operation using a feedback resistor from
00055 // input to output, a MOS inverter can be made to act like an op-amp for
00056 // small signals centered around the switching threshold.
00057 //
00058 // In 2008, Michael Huth facilitated closer investigation of the SID 6581
00059 // filter circuit by publishing high quality microscope photographs of the die.
00060 // Tommi Lempinen has done an impressive work on re-vectorizing and annotating
00061 // the die photographs, substantially simplifying further analysis of the
00062 // filter circuit.
00063 // 
00064 // The filter schematics below are reverse engineered from these re-vectorized
00065 // and annotated die photographs. While the filter first depicted in reSID 0.9
00066 // is a correct model of the basic filter, the schematics are now completed
00067 // with the audio mixer and output stage, including details on intended
00068 // relative resistor values. Also included are schematics for the NMOS FET
00069 // voltage controlled resistors (VCRs) used to control cutoff frequency, the
00070 // DAC which controls the VCRs, the NMOS op-amps, and the output buffer.
00071 //
00072 //
00073 // SID filter / mixer / output
00074 // ---------------------------
00075 // 
00076 //                ---------------------------------------------------
00077 //               |                                                   |
00078 //               |                         --1R1-- \--  D7           |
00079 //               |              ---R1--   |           |              |
00080 //               |             |       |  |--2R1-- \--| D6           |
00081 //               |    ------------<A]-----|           |     $17      |
00082 //               |   |                    |--4R1-- \--| D5  1=open   | (3.5R1)
00083 //               |   |                    |           |              |
00084 //               |   |                     --8R1-- \--| D4           | (7.0R1)
00085 //               |   |                                |              |
00086 // $17           |   |                    (CAP2B)     |  (CAP1B)     |
00087 // 0=to mixer    |    --R8--    ---R8--        ---C---|       ---C---| 
00088 // 1=to filter   |          |  |       |      |       |      |       |
00089 //                ------R8--|-----[A>--|--Rw-----[A>--|--Rw-----[A>--|
00090 //     ve (EXT IN)          |          |              |              |
00091 // D3  \ ---------------R8--|          |              | (CAP2A)      | (CAP1A)
00092 //     |   v3               |          | vhp          | vbp          | vlp
00093 // D2  |   \ -----------R8--|     -----               |              |
00094 //     |   |   v2           |    |                    |              |
00095 // D1  |   |   \ -------R8--|    |    ----------------               |
00096 //     |   |   |   v1       |    |   |                               |
00097 // D0  |   |   |   \ ---R8--     |   |    ---------------------------
00098 //     |   |   |   |             |   |   |
00099 //     R6  R6  R6  R6            R6  R6  R6
00100 //     |   |   |   | $18         |   |   |  $18
00101 //     |    \  |   | D7: 1=open   \   \   \ D6 - D4: 0=open
00102 //     |   |   |   |             |   |   |
00103 //      ---------------------------------                          12V
00104 //                 |
00105 //                 |               D3  --/ --1R2--                  |
00106 //                 |    ---R8--       |           |   ---R2--       |
00107 //                 |   |       |   D2 |--/ --2R2--|  |       |  ||--
00108 //                  ------[A>---------|           |-----[A>-----||
00109 //                                 D1 |--/ --4R2--| (4.25R2)    ||--
00110 //                        $18         |           |                 |
00111 //                        0=open   D0  --/ --8R2--  (8.75R2)        |
00112 //
00113 //                                                                  vo (AUDIO
00114 //                                                                      OUT)
00115 //
00116 //
00117 // v1  - voice 1
00118 // v2  - voice 2
00119 // v3  - voice 3
00120 // ve  - ext in
00121 // vhp - highpass output
00122 // vbp - bandpass output
00123 // vlp - lowpass output
00124 // vo  - audio out
00125 // [A> - single ended inverting op-amp (self-biased NMOS inverter)
00126 // Rn  - "resistors", implemented with custom NMOS FETs
00127 // Rw  - cutoff frequency resistor (VCR)
00128 // C   - capacitor
00129 //
00130 // Notes:
00131 //
00132 // R2  ~  2.0*R1
00133 // R6  ~  6.0*R1
00134 // R8  ~  8.0*R1
00135 // R24 ~ 24.0*R1
00136 //
00137 // The Rn "resistors" in the circuit are implemented with custom NMOS FETs,
00138 // probably because of space constraints on the SID die. The silicon substrate
00139 // is laid out in a narrow strip or "snake", with a strip length proportional
00140 // to the intended resistance. The polysilicon gate electrode covers the entire
00141 // silicon substrate and is fixed at 12V in order for the NMOS FET to operate
00142 // in triode mode (a.k.a. linear mode or ohmic mode).
00143 //
00144 // Even in "linear mode", an NMOS FET is only an approximation of a resistor,
00145 // as the apparant resistance increases with increasing drain-to-source
00146 // voltage. If the drain-to-source voltage should approach the gate voltage
00147 // of 12V, the NMOS FET will enter saturation mode (a.k.a. active mode), and
00148 // the NMOS FET will not operate anywhere like a resistor.
00149 //
00150 // 
00151 // 
00152 // NMOS FET voltage controlled resistor (VCR)
00153 // ------------------------------------------
00154 //
00155 //                Vw
00156 //
00157 //                |
00158 //                |
00159 //                R1
00160 //                |
00161 //          --R1--|
00162 //         |    __|__
00163 //         |    -----
00164 //         |    |   |
00165 // vi ----------     -------- vo
00166 //         |           |
00167 //          ----R24----
00168 //
00169 //
00170 // vi  - input
00171 // vo  - output
00172 // Rn  - "resistors", implemented with custom NMOS FETs
00173 // Vw  - voltage from 11-bit DAC (frequency cutoff control)
00174 // 
00175 // Notes:
00176 //
00177 // An approximate value for R24 can be found by using the formula for the
00178 // filter cutoff frequency:
00179 //
00180 // FCmin = 1/(2*pi*Rmax*C)
00181 //
00182 // Assuming that a the setting for minimum cutoff frequency in combination with
00183 // a low level input signal ensures that only negligible current will flow
00184 // through the transistor in the schematics above, values for FCmin and C can
00185 // be substituted in this formula to find Rmax.
00186 // Using C = 470pF and FCmin = 220Hz (measured value), we get:
00187 //
00188 // FCmin = 1/(2*pi*Rmax*C)
00189 // Rmax = 1/(2*pi*FCmin*C) = 1/(2*pi*220*470e-12) ~ 1.5MOhm
00190 //
00191 // From this it follows that:
00192 // R24 =  Rmax   ~ 1.5MOhm
00193 // R1  ~  R24/24 ~  64kOhm
00194 // R2  ~  2.0*R1 ~ 128kOhm
00195 // R6  ~  6.0*R1 ~ 384kOhm
00196 // R8  ~  8.0*R1 ~ 512kOhm
00197 //
00198 // Note that these are only approximate values for one particular SID chip,
00199 // due to process variations the values can be substantially different in
00200 // other chips.
00201 // 
00202 // 
00203 // 
00204 // Filter frequency cutoff DAC
00205 // ---------------------------
00206 //
00207 //
00208 //        12V  10   9   8   7   6   5   4   3   2   1   0   VGND
00209 //          |   |   |   |   |   |   |   |   |   |   |   |     |   Missing
00210 //         2R  2R  2R  2R  2R  2R  2R  2R  2R  2R  2R  2R    2R   termination
00211 //          |   |   |   |   |   |   |   |   |   |   |   |     |
00212 //     Vw ----R---R---R---R---R---R---R---R---R---R---R--   ---
00213 //
00214 // Bit on:  12V
00215 // Bit off:  5V (VGND)
00216 //
00217 // As is the case with all MOS 6581 DACs, the termination to (virtual) ground
00218 // at bit 0 is missing.
00219 //
00220 // Furthermore, the control of the two VCRs imposes a load on the DAC output
00221 // which varies with the input signals to the VCRs. This can be seen from the
00222 // VCR figure above.
00223 //
00224 // 
00225 // 
00226 // "Op-amp" (self-biased NMOS inverter)
00227 // ------------------------------------
00228 //                  
00229 //                  
00230 //                        12V
00231 //
00232 //                         |
00233 //              -----------|
00234 //             |           |
00235 //             |     ------|
00236 //             |    |      |
00237 //             |    |  ||--
00238 //             |     --||
00239 //             |       ||--
00240 //         ||--            |
00241 // vi -----||              |--------- vo
00242 //         ||--            |   |
00243 //             |       ||--    |
00244 //             |-------||      |
00245 //             |       ||--    |
00246 //         ||--            |   |
00247 //       --||              |   |
00248 //      |  ||--            |   |
00249 //      |      |           |   |
00250 //      |       -----------|   |
00251 //      |                  |   |
00252 //      |                      |
00253 //      |                 GND  |
00254 //      |                      |
00255 //       ----------------------
00256 //
00257 //
00258 // vi  - input
00259 // vo  - output
00260 //
00261 // Notes:
00262 //
00263 // The schematics above are laid out to show that the "op-amp" logically
00264 // consists of two building blocks; a saturated load NMOS inverter (on the
00265 // right hand side of the schematics) with a buffer / bias input stage
00266 // consisting of a variable saturated load NMOS inverter (on the left hand
00267 // side of the schematics).
00268 //
00269 // Provided a reasonably high input impedance and a reasonably low output
00270 // impedance, the "op-amp" can be modeled as a voltage transfer function
00271 // mapping input voltage to output voltage.
00272 //
00273 //
00274 //
00275 // Output buffer (NMOS voltage follower)
00276 // -------------------------------------
00277 //
00278 //
00279 //            12V
00280 //
00281 //             |
00282 //             |
00283 //         ||--
00284 // vi -----||
00285 //         ||--
00286 //             |
00287 //             |------ vo
00288 //             |     (AUDIO
00289 //            Rext    OUT)
00290 //             |
00291 //             |
00292 //
00293 //            GND
00294 //
00295 // vi   - input
00296 // vo   - output
00297 // Rext - external resistor, 1kOhm
00298 //
00299 // Notes:
00300 //
00301 // The external resistor Rext is needed to complete the NMOS voltage follower,
00302 // this resistor has a recommended value of 1kOhm.
00303 //
00304 // Die photographs show that actually, two NMOS transistors are used in the
00305 // voltage follower. However the two transistors are coupled in parallel (all
00306 // terminals are pairwise common), which implies that we can model the two
00307 // transistors as one.
00308 //
00309 // ----------------------------------------------------------------------------
00310 
00311 // Compile-time computation of op-amp summer and mixer table offsets.
00312 
00313 // The highpass summer has 2 - 6 inputs (bandpass, lowpass, and 0 - 4 voices).
00314 template<int i>
00315 struct summer_offset
00316 {
00317   enum { value = summer_offset<i - 1>::value + ((2 + i - 1) << 16) };
00318 };
00319 
00320 template<>
00321 struct summer_offset<0>
00322 {
00323   enum { value = 0 };
00324 };
00325 
00326 // The mixer has 0 - 7 inputs (0 - 4 voices and 0 - 3 filter outputs).
00327 template<int i>
00328 struct mixer_offset
00329 {
00330   enum { value = mixer_offset<i - 1>::value + ((i - 1) << 16) };
00331 };
00332               
00333 template<>
00334 struct mixer_offset<1>
00335 {
00336   enum { value = 1 };
00337 };
00338 
00339 template<>
00340 struct mixer_offset<0>
00341 {
00342   enum { value = 0 };
00343 };
00344 
00345 
00346 class Filter
00347 {
00348 public:
00349   Filter();
00350 
00351   void enable_filter(bool enable);
00352   void adjust_filter_bias(double dac_bias);
00353   void set_chip_model(chip_model model);
00354   void set_voice_mask(reg4 mask);
00355 
00356   void clock(int voice1, int voice2, int voice3);
00357   void clock(cycle_count delta_t, int voice1, int voice2, int voice3);
00358   void reset();
00359 
00360   // Write registers.
00361   void writeFC_LO(reg8);
00362   void writeFC_HI(reg8);
00363   void writeRES_FILT(reg8);
00364   void writeMODE_VOL(reg8);
00365 
00366   // SID audio input (16 bits).
00367   void input(short sample);
00368 
00369   // SID audio output (16 bits).
00370   short output();
00371 
00372 protected:
00373   void set_sum_mix();
00374   void set_w0();
00375   void set_Q();
00376 
00377   // Filter enabled.
00378   bool enabled;
00379 
00380   // Filter cutoff frequency.
00381   reg12 fc;
00382 
00383   // Filter resonance.
00384   reg8 res;
00385 
00386   // Selects which voices to route through the filter.
00387   reg8 filt;
00388 
00389   // Selects which filter outputs to route into the mixer.
00390   reg4 mode;
00391 
00392   // Output master volume.
00393   reg4 vol;
00394 
00395   // Used to mask out EXT IN if not connected, and for test purposes
00396   // (voice muting).
00397   reg8 voice_mask;
00398 
00399   // Select which inputs to route into the summer / mixer.
00400   // These are derived from filt, mode, and voice_mask.
00401   reg8 sum;
00402   reg8 mix;
00403 
00404   // State of filter.
00405   int Vhp; // highpass
00406   int Vbp; // bandpass
00407   int Vbp_x, Vbp_vc;
00408   int Vlp; // lowpass
00409   int Vlp_x, Vlp_vc;
00410   // Filter / mixer inputs.
00411   int ve;
00412   int v3;
00413   int v2;
00414   int v1;
00415 
00416   // Cutoff frequency DAC voltage, resonance.
00417   int Vddt_Vw_2, Vw_bias;
00418   int _8_div_Q;
00419   // FIXME: Temporarily used for MOS 8580 emulation.
00420   int w0;
00421   int _1024_div_Q;
00422 
00423   chip_model sid_model;
00424 
00425   typedef struct {
00426     int vo_N16;  // Fixed point scaling for 16 bit op-amp output.
00427     int kVddt;   // K*(Vdd - Vth)
00428     int n_snake;
00429     int voice_scale_s14;
00430     int voice_DC;
00431     int ak;
00432     int bk;
00433     int vc_min;
00434     int vc_max;
00435 
00436     // Reverse op-amp transfer function.
00437     unsigned short opamp_rev[1 << 16];
00438     // Lookup tables for gain and summer op-amps in output stage / filter.
00439     unsigned short summer[summer_offset<5>::value];
00440     unsigned short gain[16][1 << 16];
00441     unsigned short mixer[mixer_offset<8>::value];
00442     // Cutoff frequency DAC output voltage table. FC is an 11 bit register.
00443     unsigned short f0_dac[1 << 11];
00444   } model_filter_t;
00445 
00446   int solve_gain(int* opamp, int n, int vi_t, int& x, model_filter_t& mf);
00447   int solve_integrate_6581(int dt, int vi_t, int& x, int& vc, model_filter_t& mf);
00448 
00449   // VCR - 6581 only.
00450   static unsigned short vcr_kVg[1 << 16];
00451   static unsigned short vcr_n_Ids_term[1 << 16];
00452   // Common parameters.
00453   static model_filter_t model_filter[2];
00454 
00455 friend class SID;
00456 };
00457 
00458 
00459 // ----------------------------------------------------------------------------
00460 // Inline functions.
00461 // The following functions are defined inline because they are called every
00462 // time a sample is calculated.
00463 // ----------------------------------------------------------------------------
00464 
00465 #if RESID_INLINING || defined(RESID_FILTER_CC)
00466 
00467 // ----------------------------------------------------------------------------
00468 // SID clocking - 1 cycle.
00469 // ----------------------------------------------------------------------------
00470 RESID_INLINE
00471 void Filter::clock(int voice1, int voice2, int voice3)
00472 {
00473   model_filter_t& f = model_filter[sid_model];
00474 
00475   v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
00476   v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
00477   v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
00478 
00479   // Sum inputs routed into the filter.
00480   int Vi = 0;
00481   int offset = 0;
00482 
00483   switch (sum & 0xf) {
00484   case 0x0:
00485     Vi = 0;
00486     offset = summer_offset<0>::value;
00487     break;
00488   case 0x1:
00489     Vi = v1;
00490     offset = summer_offset<1>::value;
00491     break;
00492   case 0x2:
00493     Vi = v2;
00494     offset = summer_offset<1>::value;
00495     break;
00496   case 0x3:
00497     Vi = v2 + v1;
00498     offset = summer_offset<2>::value;
00499     break;
00500   case 0x4:
00501     Vi = v3;
00502     offset = summer_offset<1>::value;
00503     break;
00504   case 0x5:
00505     Vi = v3 + v1;
00506     offset = summer_offset<2>::value;
00507     break;
00508   case 0x6:
00509     Vi = v3 + v2;
00510     offset = summer_offset<2>::value;
00511     break;
00512   case 0x7:
00513     Vi = v3 + v2 + v1;
00514     offset = summer_offset<3>::value;
00515     break;
00516   case 0x8:
00517     Vi = ve;
00518     offset = summer_offset<1>::value;
00519     break;
00520   case 0x9:
00521     Vi = ve + v1;
00522     offset = summer_offset<2>::value;
00523     break;
00524   case 0xa:
00525     Vi = ve + v2;
00526     offset = summer_offset<2>::value;
00527     break;
00528   case 0xb:
00529     Vi = ve + v2 + v1;
00530     offset = summer_offset<3>::value;
00531     break;
00532   case 0xc:
00533     Vi = ve + v3;
00534     offset = summer_offset<2>::value;
00535     break;
00536   case 0xd:
00537     Vi = ve + v3 + v1;
00538     offset = summer_offset<3>::value;
00539     break;
00540   case 0xe:
00541     Vi = ve + v3 + v2;
00542     offset = summer_offset<3>::value;
00543     break;
00544   case 0xf:
00545     Vi = ve + v3 + v2 + v1;
00546     offset = summer_offset<4>::value;
00547     break;
00548   }
00549 
00550   // Calculate filter outputs.
00551   if (sid_model == 0) {
00552     // MOS 6581.
00553     Vlp = solve_integrate_6581(1, Vbp, Vlp_x, Vlp_vc, f);
00554     Vbp = solve_integrate_6581(1, Vhp, Vbp_x, Vbp_vc, f);
00555     Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
00556   }
00557   else {
00558     // MOS 8580. FIXME: Not yet using op-amp model.
00559 
00560     // delta_t = 1 is converted to seconds given a 1MHz clock by dividing
00561     // with 1 000 000.
00562 
00563     int dVbp = w0*(Vhp >> 4) >> 16;
00564     int dVlp = w0*(Vbp >> 4) >> 16;
00565     Vbp -= dVbp;
00566     Vlp -= dVlp;
00567     Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
00568   }
00569 }
00570 
00571 // ----------------------------------------------------------------------------
00572 // SID clocking - delta_t cycles.
00573 // ----------------------------------------------------------------------------
00574 RESID_INLINE
00575 void Filter::clock(cycle_count delta_t, int voice1, int voice2, int voice3)
00576 {
00577   model_filter_t& f = model_filter[sid_model];
00578 
00579   v1 = (voice1*f.voice_scale_s14 >> 18) + f.voice_DC;
00580   v2 = (voice2*f.voice_scale_s14 >> 18) + f.voice_DC;
00581   v3 = (voice3*f.voice_scale_s14 >> 18) + f.voice_DC;
00582 
00583   // Enable filter on/off.
00584   // This is not really part of SID, but is useful for testing.
00585   // On slow CPUs it may be necessary to bypass the filter to lower the CPU
00586   // load.
00587   if (unlikely(!enabled)) {
00588     return;
00589   }
00590 
00591   // Sum inputs routed into the filter.
00592   int Vi = 0;
00593   int offset = 0;
00594 
00595   switch (sum & 0xf) {
00596   case 0x0:
00597     Vi = 0;
00598     offset = summer_offset<0>::value;
00599     break;
00600   case 0x1:
00601     Vi = v1;
00602     offset = summer_offset<1>::value;
00603     break;
00604   case 0x2:
00605     Vi = v2;
00606     offset = summer_offset<1>::value;
00607     break;
00608   case 0x3:
00609     Vi = v2 + v1;
00610     offset = summer_offset<2>::value;
00611     break;
00612   case 0x4:
00613     Vi = v3;
00614     offset = summer_offset<1>::value;
00615     break;
00616   case 0x5:
00617     Vi = v3 + v1;
00618     offset = summer_offset<2>::value;
00619     break;
00620   case 0x6:
00621     Vi = v3 + v2;
00622     offset = summer_offset<2>::value;
00623     break;
00624   case 0x7:
00625     Vi = v3 + v2 + v1;
00626     offset = summer_offset<3>::value;
00627     break;
00628   case 0x8:
00629     Vi = ve;
00630     offset = summer_offset<1>::value;
00631     break;
00632   case 0x9:
00633     Vi = ve + v1;
00634     offset = summer_offset<2>::value;
00635     break;
00636   case 0xa:
00637     Vi = ve + v2;
00638     offset = summer_offset<2>::value;
00639     break;
00640   case 0xb:
00641     Vi = ve + v2 + v1;
00642     offset = summer_offset<3>::value;
00643     break;
00644   case 0xc:
00645     Vi = ve + v3;
00646     offset = summer_offset<2>::value;
00647     break;
00648   case 0xd:
00649     Vi = ve + v3 + v1;
00650     offset = summer_offset<3>::value;
00651     break;
00652   case 0xe:
00653     Vi = ve + v3 + v2;
00654     offset = summer_offset<3>::value;
00655     break;
00656   case 0xf:
00657     Vi = ve + v3 + v2 + v1;
00658     offset = summer_offset<4>::value;
00659     break;
00660   }
00661 
00662   // Maximum delta cycles for filter fixpoint iteration to converge
00663   // is approximately 3.
00664   cycle_count delta_t_flt = 3;
00665 
00666   if (sid_model == 0) {
00667     // MOS 6581.
00668     while (delta_t) {
00669       if (unlikely(delta_t < delta_t_flt)) {
00670         delta_t_flt = delta_t;
00671       }
00672 
00673       // Calculate filter outputs.
00674       Vlp = solve_integrate_6581(delta_t_flt, Vbp, Vlp_x, Vlp_vc, f);
00675       Vbp = solve_integrate_6581(delta_t_flt, Vhp, Vbp_x, Vbp_vc, f);
00676       Vhp = f.summer[offset + f.gain[_8_div_Q][Vbp] + Vlp + Vi];
00677 
00678       delta_t -= delta_t_flt;
00679     }
00680   }
00681   else {
00682     // MOS 8580. FIXME: Not yet using op-amp model.
00683     while (delta_t) {
00684       if (delta_t < delta_t_flt) {
00685         delta_t_flt = delta_t;
00686       }
00687 
00688       // delta_t is converted to seconds given a 1MHz clock by dividing
00689       // with 1 000 000. This is done in two operations to avoid integer
00690       // multiplication overflow.
00691 
00692       // Calculate filter outputs.
00693       int w0_delta_t = w0*delta_t_flt >> 2;
00694 
00695       int dVbp = w0_delta_t*(Vhp >> 4) >> 14;
00696       int dVlp = w0_delta_t*(Vbp >> 4) >> 14;
00697       Vbp -= dVbp;
00698       Vlp -= dVlp;
00699       Vhp = (Vbp*_1024_div_Q >> 10) - Vlp - Vi;
00700 
00701       delta_t -= delta_t_flt;
00702     }
00703   }
00704 }
00705 
00706 
00707 // ----------------------------------------------------------------------------
00708 // SID audio input (16 bits).
00709 // ----------------------------------------------------------------------------
00710 RESID_INLINE
00711 void Filter::input(short sample)
00712 {
00713   // Scale to three times the peak-to-peak for one voice and add the op-amp
00714   // "zero" DC level.
00715   // NB! Adding the op-amp "zero" DC level is a (wildly inaccurate)
00716   // approximation of feeding the input through an AC coupling capacitor.
00717   // This could be implemented as a separate filter circuit, however the
00718   // primary use of the emulator is not to process external signals.
00719   // The upside is that the MOS8580 "digi boost" works without a separate (DC)
00720   // input interface.
00721   // Note that the input is 16 bits, compared to the 20 bit voice output.
00722   model_filter_t& f = model_filter[sid_model];
00723   ve = (sample*f.voice_scale_s14*3 >> 14) + f.mixer[0];
00724 }
00725 
00726 
00727 // ----------------------------------------------------------------------------
00728 // SID audio output (16 bits).
00729 // ----------------------------------------------------------------------------
00730 RESID_INLINE
00731 short Filter::output()
00732 {
00733   model_filter_t& f = model_filter[sid_model];
00734 
00735   // Writing the switch below manually would be tedious and error-prone;
00736   // it is rather generated by the following Perl program:
00737 
00738   /*
00739 my @i = qw(v1 v2 v3 ve Vlp Vbp Vhp);
00740 for my $mix (0..2**@i-1) {
00741     print sprintf("  case 0x%02x:\n", $mix);
00742     my @sum;
00743     for (@i) {
00744         unshift(@sum, $_) if $mix & 0x01;
00745         $mix >>= 1;
00746     }
00747     my $sum = join(" + ", @sum) || "0";
00748     print "    Vi = $sum;\n";
00749     print "    offset = mixer_offset<" . @sum . ">::value;\n";
00750     print "    break;\n";
00751 }
00752   */
00753 
00754   // Sum inputs routed into the mixer.
00755   int Vi = 0;
00756   int offset = 0;
00757 
00758   switch (mix & 0x7f) {
00759   case 0x00:
00760     Vi = 0;
00761     offset = mixer_offset<0>::value;
00762     break;
00763   case 0x01:
00764     Vi = v1;
00765     offset = mixer_offset<1>::value;
00766     break;
00767   case 0x02:
00768     Vi = v2;
00769     offset = mixer_offset<1>::value;
00770     break;
00771   case 0x03:
00772     Vi = v2 + v1;
00773     offset = mixer_offset<2>::value;
00774     break;
00775   case 0x04:
00776     Vi = v3;
00777     offset = mixer_offset<1>::value;
00778     break;
00779   case 0x05:
00780     Vi = v3 + v1;
00781     offset = mixer_offset<2>::value;
00782     break;
00783   case 0x06:
00784     Vi = v3 + v2;
00785     offset = mixer_offset<2>::value;
00786     break;
00787   case 0x07:
00788     Vi = v3 + v2 + v1;
00789     offset = mixer_offset<3>::value;
00790     break;
00791   case 0x08:
00792     Vi = ve;
00793     offset = mixer_offset<1>::value;
00794     break;
00795   case 0x09:
00796     Vi = ve + v1;
00797     offset = mixer_offset<2>::value;
00798     break;
00799   case 0x0a:
00800     Vi = ve + v2;
00801     offset = mixer_offset<2>::value;
00802     break;
00803   case 0x0b:
00804     Vi = ve + v2 + v1;
00805     offset = mixer_offset<3>::value;
00806     break;
00807   case 0x0c:
00808     Vi = ve + v3;
00809     offset = mixer_offset<2>::value;
00810     break;
00811   case 0x0d:
00812     Vi = ve + v3 + v1;
00813     offset = mixer_offset<3>::value;
00814     break;
00815   case 0x0e:
00816     Vi = ve + v3 + v2;
00817     offset = mixer_offset<3>::value;
00818     break;
00819   case 0x0f:
00820     Vi = ve + v3 + v2 + v1;
00821     offset = mixer_offset<4>::value;
00822     break;
00823   case 0x10:
00824     Vi = Vlp;
00825     offset = mixer_offset<1>::value;
00826     break;
00827   case 0x11:
00828     Vi = Vlp + v1;
00829     offset = mixer_offset<2>::value;
00830     break;
00831   case 0x12:
00832     Vi = Vlp + v2;
00833     offset = mixer_offset<2>::value;
00834     break;
00835   case 0x13:
00836     Vi = Vlp + v2 + v1;
00837     offset = mixer_offset<3>::value;
00838     break;
00839   case 0x14:
00840     Vi = Vlp + v3;
00841     offset = mixer_offset<2>::value;
00842     break;
00843   case 0x15:
00844     Vi = Vlp + v3 + v1;
00845     offset = mixer_offset<3>::value;
00846     break;
00847   case 0x16:
00848     Vi = Vlp + v3 + v2;
00849     offset = mixer_offset<3>::value;
00850     break;
00851   case 0x17:
00852     Vi = Vlp + v3 + v2 + v1;
00853     offset = mixer_offset<4>::value;
00854     break;
00855   case 0x18:
00856     Vi = Vlp + ve;
00857     offset = mixer_offset<2>::value;
00858     break;
00859   case 0x19:
00860     Vi = Vlp + ve + v1;
00861     offset = mixer_offset<3>::value;
00862     break;
00863   case 0x1a:
00864     Vi = Vlp + ve + v2;
00865     offset = mixer_offset<3>::value;
00866     break;
00867   case 0x1b:
00868     Vi = Vlp + ve + v2 + v1;
00869     offset = mixer_offset<4>::value;
00870     break;
00871   case 0x1c:
00872     Vi = Vlp + ve + v3;
00873     offset = mixer_offset<3>::value;
00874     break;
00875   case 0x1d:
00876     Vi = Vlp + ve + v3 + v1;
00877     offset = mixer_offset<4>::value;
00878     break;
00879   case 0x1e:
00880     Vi = Vlp + ve + v3 + v2;
00881     offset = mixer_offset<4>::value;
00882     break;
00883   case 0x1f:
00884     Vi = Vlp + ve + v3 + v2 + v1;
00885     offset = mixer_offset<5>::value;
00886     break;
00887   case 0x20:
00888     Vi = Vbp;
00889     offset = mixer_offset<1>::value;
00890     break;
00891   case 0x21:
00892     Vi = Vbp + v1;
00893     offset = mixer_offset<2>::value;
00894     break;
00895   case 0x22:
00896     Vi = Vbp + v2;
00897     offset = mixer_offset<2>::value;
00898     break;
00899   case 0x23:
00900     Vi = Vbp + v2 + v1;
00901     offset = mixer_offset<3>::value;
00902     break;
00903   case 0x24:
00904     Vi = Vbp + v3;
00905     offset = mixer_offset<2>::value;
00906     break;
00907   case 0x25:
00908     Vi = Vbp + v3 + v1;
00909     offset = mixer_offset<3>::value;
00910     break;
00911   case 0x26:
00912     Vi = Vbp + v3 + v2;
00913     offset = mixer_offset<3>::value;
00914     break;
00915   case 0x27:
00916     Vi = Vbp + v3 + v2 + v1;
00917     offset = mixer_offset<4>::value;
00918     break;
00919   case 0x28:
00920     Vi = Vbp + ve;
00921     offset = mixer_offset<2>::value;
00922     break;
00923   case 0x29:
00924     Vi = Vbp + ve + v1;
00925     offset = mixer_offset<3>::value;
00926     break;
00927   case 0x2a:
00928     Vi = Vbp + ve + v2;
00929     offset = mixer_offset<3>::value;
00930     break;
00931   case 0x2b:
00932     Vi = Vbp + ve + v2 + v1;
00933     offset = mixer_offset<4>::value;
00934     break;
00935   case 0x2c:
00936     Vi = Vbp + ve + v3;
00937     offset = mixer_offset<3>::value;
00938     break;
00939   case 0x2d:
00940     Vi = Vbp + ve + v3 + v1;
00941     offset = mixer_offset<4>::value;
00942     break;
00943   case 0x2e:
00944     Vi = Vbp + ve + v3 + v2;
00945     offset = mixer_offset<4>::value;
00946     break;
00947   case 0x2f:
00948     Vi = Vbp + ve + v3 + v2 + v1;
00949     offset = mixer_offset<5>::value;
00950     break;
00951   case 0x30:
00952     Vi = Vbp + Vlp;
00953     offset = mixer_offset<2>::value;
00954     break;
00955   case 0x31:
00956     Vi = Vbp + Vlp + v1;
00957     offset = mixer_offset<3>::value;
00958     break;
00959   case 0x32:
00960     Vi = Vbp + Vlp + v2;
00961     offset = mixer_offset<3>::value;
00962     break;
00963   case 0x33:
00964     Vi = Vbp + Vlp + v2 + v1;
00965     offset = mixer_offset<4>::value;
00966     break;
00967   case 0x34:
00968     Vi = Vbp + Vlp + v3;
00969     offset = mixer_offset<3>::value;
00970     break;
00971   case 0x35:
00972     Vi = Vbp + Vlp + v3 + v1;
00973     offset = mixer_offset<4>::value;
00974     break;
00975   case 0x36:
00976     Vi = Vbp + Vlp + v3 + v2;
00977     offset = mixer_offset<4>::value;
00978     break;
00979   case 0x37:
00980     Vi = Vbp + Vlp + v3 + v2 + v1;
00981     offset = mixer_offset<5>::value;
00982     break;
00983   case 0x38:
00984     Vi = Vbp + Vlp + ve;
00985     offset = mixer_offset<3>::value;
00986     break;
00987   case 0x39:
00988     Vi = Vbp + Vlp + ve + v1;
00989     offset = mixer_offset<4>::value;
00990     break;
00991   case 0x3a:
00992     Vi = Vbp + Vlp + ve + v2;
00993     offset = mixer_offset<4>::value;
00994     break;
00995   case 0x3b:
00996     Vi = Vbp + Vlp + ve + v2 + v1;
00997     offset = mixer_offset<5>::value;
00998     break;
00999   case 0x3c:
01000     Vi = Vbp + Vlp + ve + v3;
01001     offset = mixer_offset<4>::value;
01002     break;
01003   case 0x3d:
01004     Vi = Vbp + Vlp + ve + v3 + v1;
01005     offset = mixer_offset<5>::value;
01006     break;
01007   case 0x3e:
01008     Vi = Vbp + Vlp + ve + v3 + v2;
01009     offset = mixer_offset<5>::value;
01010     break;
01011   case 0x3f:
01012     Vi = Vbp + Vlp + ve + v3 + v2 + v1;
01013     offset = mixer_offset<6>::value;
01014     break;
01015   case 0x40:
01016     Vi = Vhp;
01017     offset = mixer_offset<1>::value;
01018     break;
01019   case 0x41:
01020     Vi = Vhp + v1;
01021     offset = mixer_offset<2>::value;
01022     break;
01023   case 0x42:
01024     Vi = Vhp + v2;
01025     offset = mixer_offset<2>::value;
01026     break;
01027   case 0x43:
01028     Vi = Vhp + v2 + v1;
01029     offset = mixer_offset<3>::value;
01030     break;
01031   case 0x44:
01032     Vi = Vhp + v3;
01033     offset = mixer_offset<2>::value;
01034     break;
01035   case 0x45:
01036     Vi = Vhp + v3 + v1;
01037     offset = mixer_offset<3>::value;
01038     break;
01039   case 0x46:
01040     Vi = Vhp + v3 + v2;
01041     offset = mixer_offset<3>::value;
01042     break;
01043   case 0x47:
01044     Vi = Vhp + v3 + v2 + v1;
01045     offset = mixer_offset<4>::value;
01046     break;
01047   case 0x48:
01048     Vi = Vhp + ve;
01049     offset = mixer_offset<2>::value;
01050     break;
01051   case 0x49:
01052     Vi = Vhp + ve + v1;
01053     offset = mixer_offset<3>::value;
01054     break;
01055   case 0x4a:
01056     Vi = Vhp + ve + v2;
01057     offset = mixer_offset<3>::value;
01058     break;
01059   case 0x4b:
01060     Vi = Vhp + ve + v2 + v1;
01061     offset = mixer_offset<4>::value;
01062     break;
01063   case 0x4c:
01064     Vi = Vhp + ve + v3;
01065     offset = mixer_offset<3>::value;
01066     break;
01067   case 0x4d:
01068     Vi = Vhp + ve + v3 + v1;
01069     offset = mixer_offset<4>::value;
01070     break;
01071   case 0x4e:
01072     Vi = Vhp + ve + v3 + v2;
01073     offset = mixer_offset<4>::value;
01074     break;
01075   case 0x4f:
01076     Vi = Vhp + ve + v3 + v2 + v1;
01077     offset = mixer_offset<5>::value;
01078     break;
01079   case 0x50:
01080     Vi = Vhp + Vlp;
01081     offset = mixer_offset<2>::value;
01082     break;
01083   case 0x51:
01084     Vi = Vhp + Vlp + v1;
01085     offset = mixer_offset<3>::value;
01086     break;
01087   case 0x52:
01088     Vi = Vhp + Vlp + v2;
01089     offset = mixer_offset<3>::value;
01090     break;
01091   case 0x53:
01092     Vi = Vhp + Vlp + v2 + v1;
01093     offset = mixer_offset<4>::value;
01094     break;
01095   case 0x54:
01096     Vi = Vhp + Vlp + v3;
01097     offset = mixer_offset<3>::value;
01098     break;
01099   case 0x55:
01100     Vi = Vhp + Vlp + v3 + v1;
01101     offset = mixer_offset<4>::value;
01102     break;
01103   case 0x56:
01104     Vi = Vhp + Vlp + v3 + v2;
01105     offset = mixer_offset<4>::value;
01106     break;
01107   case 0x57:
01108     Vi = Vhp + Vlp + v3 + v2 + v1;
01109     offset = mixer_offset<5>::value;
01110     break;
01111   case 0x58:
01112     Vi = Vhp + Vlp + ve;
01113     offset = mixer_offset<3>::value;
01114     break;
01115   case 0x59:
01116     Vi = Vhp + Vlp + ve + v1;
01117     offset = mixer_offset<4>::value;
01118     break;
01119   case 0x5a:
01120     Vi = Vhp + Vlp + ve + v2;
01121     offset = mixer_offset<4>::value;
01122     break;
01123   case 0x5b:
01124     Vi = Vhp + Vlp + ve + v2 + v1;
01125     offset = mixer_offset<5>::value;
01126     break;
01127   case 0x5c:
01128     Vi = Vhp + Vlp + ve + v3;
01129     offset = mixer_offset<4>::value;
01130     break;
01131   case 0x5d:
01132     Vi = Vhp + Vlp + ve + v3 + v1;
01133     offset = mixer_offset<5>::value;
01134     break;
01135   case 0x5e:
01136     Vi = Vhp + Vlp + ve + v3 + v2;
01137     offset = mixer_offset<5>::value;
01138     break;
01139   case 0x5f:
01140     Vi = Vhp + Vlp + ve + v3 + v2 + v1;
01141     offset = mixer_offset<6>::value;
01142     break;
01143   case 0x60:
01144     Vi = Vhp + Vbp;
01145     offset = mixer_offset<2>::value;
01146     break;
01147   case 0x61:
01148     Vi = Vhp + Vbp + v1;
01149     offset = mixer_offset<3>::value;
01150     break;
01151   case 0x62:
01152     Vi = Vhp + Vbp + v2;
01153     offset = mixer_offset<3>::value;
01154     break;
01155   case 0x63:
01156     Vi = Vhp + Vbp + v2 + v1;
01157     offset = mixer_offset<4>::value;
01158     break;
01159   case 0x64:
01160     Vi = Vhp + Vbp + v3;
01161     offset = mixer_offset<3>::value;
01162     break;
01163   case 0x65:
01164     Vi = Vhp + Vbp + v3 + v1;
01165     offset = mixer_offset<4>::value;
01166     break;
01167   case 0x66:
01168     Vi = Vhp + Vbp + v3 + v2;
01169     offset = mixer_offset<4>::value;
01170     break;
01171   case 0x67:
01172     Vi = Vhp + Vbp + v3 + v2 + v1;
01173     offset = mixer_offset<5>::value;
01174     break;
01175   case 0x68:
01176     Vi = Vhp + Vbp + ve;
01177     offset = mixer_offset<3>::value;
01178     break;
01179   case 0x69:
01180     Vi = Vhp + Vbp + ve + v1;
01181     offset = mixer_offset<4>::value;
01182     break;
01183   case 0x6a:
01184     Vi = Vhp + Vbp + ve + v2;
01185     offset = mixer_offset<4>::value;
01186     break;
01187   case 0x6b:
01188     Vi = Vhp + Vbp + ve + v2 + v1;
01189     offset = mixer_offset<5>::value;
01190     break;
01191   case 0x6c:
01192     Vi = Vhp + Vbp + ve + v3;
01193     offset = mixer_offset<4>::value;
01194     break;
01195   case 0x6d:
01196     Vi = Vhp + Vbp + ve + v3 + v1;
01197     offset = mixer_offset<5>::value;
01198     break;
01199   case 0x6e:
01200     Vi = Vhp + Vbp + ve + v3 + v2;
01201     offset = mixer_offset<5>::value;
01202     break;
01203   case 0x6f:
01204     Vi = Vhp + Vbp + ve + v3 + v2 + v1;
01205     offset = mixer_offset<6>::value;
01206     break;
01207   case 0x70:
01208     Vi = Vhp + Vbp + Vlp;
01209     offset = mixer_offset<3>::value;
01210     break;
01211   case 0x71:
01212     Vi = Vhp + Vbp + Vlp + v1;
01213     offset = mixer_offset<4>::value;
01214     break;
01215   case 0x72:
01216     Vi = Vhp + Vbp + Vlp + v2;
01217     offset = mixer_offset<4>::value;
01218     break;
01219   case 0x73:
01220     Vi = Vhp + Vbp + Vlp + v2 + v1;
01221     offset = mixer_offset<5>::value;
01222     break;
01223   case 0x74:
01224     Vi = Vhp + Vbp + Vlp + v3;
01225     offset = mixer_offset<4>::value;
01226     break;
01227   case 0x75:
01228     Vi = Vhp + Vbp + Vlp + v3 + v1;
01229     offset = mixer_offset<5>::value;
01230     break;
01231   case 0x76:
01232     Vi = Vhp + Vbp + Vlp + v3 + v2;
01233     offset = mixer_offset<5>::value;
01234     break;
01235   case 0x77:
01236     Vi = Vhp + Vbp + Vlp + v3 + v2 + v1;
01237     offset = mixer_offset<6>::value;
01238     break;
01239   case 0x78:
01240     Vi = Vhp + Vbp + Vlp + ve;
01241     offset = mixer_offset<4>::value;
01242     break;
01243   case 0x79:
01244     Vi = Vhp + Vbp + Vlp + ve + v1;
01245     offset = mixer_offset<5>::value;
01246     break;
01247   case 0x7a:
01248     Vi = Vhp + Vbp + Vlp + ve + v2;
01249     offset = mixer_offset<5>::value;
01250     break;
01251   case 0x7b:
01252     Vi = Vhp + Vbp + Vlp + ve + v2 + v1;
01253     offset = mixer_offset<6>::value;
01254     break;
01255   case 0x7c:
01256     Vi = Vhp + Vbp + Vlp + ve + v3;
01257     offset = mixer_offset<5>::value;
01258     break;
01259   case 0x7d:
01260     Vi = Vhp + Vbp + Vlp + ve + v3 + v1;
01261     offset = mixer_offset<6>::value;
01262     break;
01263   case 0x7e:
01264     Vi = Vhp + Vbp + Vlp + ve + v3 + v2;
01265     offset = mixer_offset<6>::value;
01266     break;
01267   case 0x7f:
01268     Vi = Vhp + Vbp + Vlp + ve + v3 + v2 + v1;
01269     offset = mixer_offset<7>::value;
01270     break;
01271   }
01272 
01273   // Sum the inputs in the mixer and run the mixer output through the gain.
01274   if (sid_model == 0) {
01275     return (short)(f.gain[vol][f.mixer[offset + Vi]] - (1 << 15));
01276   }
01277   else {
01278     // FIXME: Temporary code for MOS 8580, should use code above.
01279     return Vi*vol >> 4;
01280   }
01281 }
01282 
01283 
01284 /*
01285 Find output voltage in inverting gain and inverting summer SID op-amp
01286 circuits, using a combination of Newton-Raphson and bisection.
01287 
01288              ---R2--
01289             |       |
01290   vi ---R1-----[A>----- vo
01291             vx
01292 
01293 From Kirchoff's current law it follows that
01294 
01295   IR1f + IR2r = 0
01296 
01297 Substituting the triode mode transistor model K*W/L*(Vgst^2 - Vgdt^2)
01298 for the currents, we get:
01299 
01300   n*((Vddt - vx)^2 - (Vddt - vi)^2) + (Vddt - vx)^2 - (Vddt - vo)^2 = 0
01301 
01302 Our root function f can thus be written as:
01303 
01304   f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - vo)^2 = 0
01305 
01306 We are using the mapping function x = vo - vx -> vx. We thus substitute
01307 for vo = vx + x and get:
01308 
01309   f = (n + 1)*(Vddt - vx)^2 - n*(Vddt - vi)^2 - (Vddt - (vx + x))^2 = 0
01310 
01311 Using substitution constants
01312 
01313   a = n + 1
01314   b = Vddt
01315   c = n*(Vddt - vi)^2
01316 
01317 the equations for the root function and its derivative can be written as:
01318 
01319   f = a*(b - vx)^2 - c - (b - (vx + x))^2
01320   df = 2*((b - (vx + x))*(dvx + 1) - a*(b - vx)*dvx)
01321 */
01322 RESID_INLINE
01323 int Filter::solve_gain(int* opamp, int n, int vi, int& x, model_filter_t& mf)
01324 {
01325   // Note that all variables are translated and scaled in order to fit
01326   // in 16 bits. It is not necessary to explicitly translate the variables here,
01327   // since they are all used in subtractions which cancel out the translation:
01328   // (a - t) - (b - t) = a - b
01329 
01330   // Start off with an estimate of x and a root bracket [ak, bk].
01331   // f is increasing, so that f(ak) < 0 and f(bk) > 0.
01332   int ak = mf.ak, bk = mf.bk;
01333 
01334   int a = n + (1 << 7);              // Scaled by 2^7
01335   int b = mf.kVddt;                  // Scaled by m*2^16
01336   int b_vi = b - vi;                 // Scaled by m*2^16
01337   if (b_vi < 0) b_vi = 0;
01338   int c = n*int(unsigned(b_vi)*unsigned(b_vi) >> 12);    // Scaled by m^2*2^27
01339 
01340   for (;;) {
01341     int xk = x;
01342 
01343     // Calculate f and df.
01344     int vx_dvx = opamp[x];
01345     int vx = vx_dvx & 0xffff;  // Scaled by m*2^16
01346     int dvx = vx_dvx >> 16;    // Scaled by 2^11
01347 
01348     // f = a*(b - vx)^2 - c - (b - vo)^2
01349     // df = 2*((b - vo)*(dvx + 1) - a*(b - vx)*dvx)
01350     //
01351     int vo = vx + (x << 1) - (1 << 16);
01352     if (vo >= (1 << 16)) {
01353       vo = (1 << 16) - 1;
01354     }
01355     else if (vo < 0) {
01356       vo = 0;
01357     }
01358     int b_vx = b - vx;
01359     if (b_vx < 0) b_vx = 0;
01360     int b_vo = b - vo;
01361     if (b_vo < 0) b_vo = 0;
01362     // The dividend is scaled by m^2*2^27.
01363     int f = a*int(unsigned(b_vx)*unsigned(b_vx) >> 12) - c - int(unsigned(b_vo)*unsigned(b_vo) >> 5);
01364     // The divisor is scaled by m*2^11.
01365     int df = (b_vo*(dvx + (1 << 11)) - a*(b_vx*dvx >> 7)) >> 15;
01366     // The resulting quotient is thus scaled by m*2^16.
01367 
01368     // Newton-Raphson step: xk1 = xk - f(xk)/f'(xk)
01369     x -= f/df;
01370     if (unlikely(x == xk)) {
01371       // No further root improvement possible.
01372       return vo;
01373     }
01374 
01375     // Narrow down root bracket.
01376     if (f < 0) {
01377       // f(xk) < 0
01378       ak = xk;
01379     }
01380     else {
01381       // f(xk) > 0
01382       bk = xk;
01383     }
01384 
01385     if (unlikely(x <= ak) || unlikely(x >= bk)) {
01386       // Bisection step (ala Dekker's method).
01387       x = (ak + bk) >> 1;
01388       if (unlikely(x == ak)) {
01389         // No further bisection possible.
01390         return vo;
01391       }
01392     }
01393   }
01394 }
01395 
01396 
01397 /*
01398 Find output voltage in inverting integrator SID op-amp circuits, using a
01399 single fixpoint iteration step.
01400 
01401 A circuit diagram of a MOS 6581 integrator is shown below.
01402 
01403                  ---C---
01404                 |       |
01405   vi -----Rw-------[A>----- vo
01406        |      | vx
01407         --Rs--
01408 
01409 From Kirchoff's current law it follows that
01410 
01411   IRw + IRs + ICr = 0
01412 
01413 Using the formula for current through a capacitor, i = C*dv/dt, we get
01414 
01415   IRw + IRs + C*(vc - vc0)/dt = 0
01416   dt/C*(IRw + IRs) + vc - vc0 = 0
01417   vc = vc0 - n*(IRw(vi,vx) + IRs(vi,vx))
01418 
01419 which may be rewritten as the following iterative fixpoint function:
01420 
01421   vc = vc0 - n*(IRw(vi,g(vc)) + IRs(vi,g(vc)))
01422 
01423 To accurately calculate the currents through Rs and Rw, we need to use
01424 transistor models. Rs has a gate voltage of Vdd = 12V, and can be
01425 assumed to always be in triode mode. For Rw, the situation is rather
01426 more complex, as it turns out that this transistor will operate in
01427 both subthreshold, triode, and saturation modes.
01428 
01429 The Shichman-Hodges transistor model routinely used in textbooks may
01430 be written as follows:
01431 
01432   Ids = 0                          , Vgst < 0               (subthreshold mode) 
01433   Ids = K/2*W/L*(2*Vgst - Vds)*Vds , Vgst >= 0, Vds < Vgst  (triode mode)
01434   Ids = K/2*W/L*Vgst^2             , Vgst >= 0, Vds >= Vgst (saturation mode)
01435 
01436   where
01437   K   = u*Cox (conductance)
01438   W/L = ratio between substrate width and length
01439   Vgst = Vg - Vs - Vt (overdrive voltage)
01440 
01441 This transistor model is also called the quadratic model.
01442 
01443 Note that the equation for the triode mode can be reformulated as
01444 independent terms depending on Vgs and Vgd, respectively, by the
01445 following substitution:
01446 
01447   Vds = Vgst - (Vgst - Vds) = Vgst - Vgdt
01448 
01449   Ids = K*W/L*(2*Vgst - Vds)*Vds
01450   = K*W/L*(2*Vgst - (Vgst - Vgdt)*(Vgst - Vgdt)
01451   = K*W/L*(Vgst + Vgdt)*(Vgst - Vgdt)
01452   = K*W/L*(Vgst^2 - Vgdt^2)
01453 
01454 This turns out to be a general equation which covers both the triode
01455 and saturation modes (where the second term is 0 in saturation mode).
01456 The equation is also symmetrical, i.e. it can calculate negative
01457 currents without any change of parameters (since the terms for drain
01458 and source are identical except for the sign).
01459 
01460 FIXME: Subthreshold as function of Vgs, Vgd.
01461   Ids = I0*e^(Vgst/(n*VT))       , Vgst < 0               (subthreshold mode) 
01462 
01463 The remaining problem with the textbook model is that the transition
01464 from subthreshold the triode/saturation is not continuous.
01465 
01466 Realizing that the subthreshold and triode/saturation modes may both
01467 be defined by independent (and equal) terms of Vgs and Vds,
01468 respectively, the corresponding terms can be blended into (equal)
01469 continuous functions suitable for table lookup.
01470 
01471 The EKV model (Enz, Krummenacher and Vittoz) essentially performs this
01472 blending using an elegant mathematical formulation:
01473 
01474   Ids = Is*(if - ir)
01475   Is = 2*u*Cox*Ut^2/k*W/L
01476   if = ln^2(1 + e^((k*(Vg - Vt) - Vs)/(2*Ut))
01477   ir = ln^2(1 + e^((k*(Vg - Vt) - Vd)/(2*Ut))
01478 
01479 For our purposes, the EKV model preserves two important properties
01480 discussed above:
01481 
01482 - It consists of two independent terms, which can be represented by
01483   the same lookup table.
01484 - It is symmetrical, i.e. it calculates current in both directions,
01485   facilitating a branch-free implementation.
01486 
01487 Rw in the circuit diagram above is a VCR (voltage controlled resistor),
01488 as shown in the circuit diagram below.
01489 
01490                    Vw
01491                    
01492                    |
01493            Vdd     |
01494               |---|  
01495              _|_   |
01496            --    --| Vg
01497           |      __|__
01498           |      -----  Rw
01499           |      |   |
01500   vi ------------     -------- vo
01501 
01502 
01503 In order to calculalate the current through the VCR, its gate voltage
01504 must be determined.
01505 
01506 Assuming triode mode and applying Kirchoff's current law, we get the
01507 following equation for Vg:
01508 
01509 u*Cox/2*W/L*((Vddt - Vg)^2 - (Vddt - vi)^2 + (Vddt - Vg)^2 - (Vddt - Vw)^2) = 0
01510 2*(Vddt - Vg)^2 - (Vddt - vi)^2 - (Vddt - Vw)^2 = 0
01511 (Vddt - Vg) = sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
01512 
01513 Vg = Vddt - sqrt(((Vddt - vi)^2 + (Vddt - Vw)^2)/2)
01514 
01515 */
01516 RESID_INLINE
01517 int Filter::solve_integrate_6581(int dt, int vi, int& vx, int& vc,
01518                                  model_filter_t& mf)
01519 {
01520   // Note that all variables are translated and scaled in order to fit
01521   // in 16 bits. It is not necessary to explicitly translate the variables here,
01522   // since they are all used in subtractions which cancel out the translation:
01523   // (a - t) - (b - t) = a - b
01524 
01525   int kVddt = mf.kVddt;      // Scaled by m*2^16
01526 
01527   // "Snake" voltages for triode mode calculation.
01528   unsigned int Vgst = kVddt - vx;
01529   unsigned int Vgdt = kVddt - vi;
01530   unsigned int Vgdt_2 = Vgdt*Vgdt;
01531 
01532   // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
01533   int n_I_snake = mf.n_snake*(int(Vgst*Vgst - Vgdt_2) >> 15);
01534 
01535   // VCR gate voltage.       // Scaled by m*2^16
01536   // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
01537   int kVg = vcr_kVg[(Vddt_Vw_2 + (Vgdt_2 >> 1)) >> 16];
01538 
01539   // VCR voltages for EKV model table lookup.
01540   int Vgs = kVg - vx;
01541   if (Vgs < 0) Vgs = 0;
01542   int Vgd = kVg - vi;
01543   if (Vgd < 0) Vgd = 0;
01544 
01545   // VCR current, scaled by m*2^15*2^15 = m*2^30
01546   int n_I_vcr = (vcr_n_Ids_term[Vgs] - vcr_n_Ids_term[Vgd]) << 15;
01547 
01548   // Change in capacitor charge.
01549   vc -= (n_I_snake + n_I_vcr)*dt;
01550 
01551 /*
01552   // FIXME: Determine whether this check is necessary.
01553   if (vc < mf.vc_min) {
01554     vc = mf.vc_min;
01555   }
01556   else if (vc > mf.vc_max) {
01557     vc = mf.vc_max;
01558   }
01559 */
01560 
01561   // vx = g(vc)
01562   vx = mf.opamp_rev[(vc >> 15) + (1 << 15)];
01563 
01564   // Return vo.
01565   return vx + (vc >> 14);
01566 }
01567 
01568 #endif // RESID_INLINING || defined(RESID_FILTER_CC)
01569 
01570 } // namespace reSID
01571 
01572 #endif // not RESID_FILTER_H