1 //  This file is part of reSID, a MOS6581 SID emulator engine.
2 //  Copyright (C) 2004  Dag Lem <resid@nimrod.no>
3 //
4 //  This program is free software; you can redistribute it and/or modify
5 //  it under the terms of the GNU General Public License as published by
6 //  the Free Software Foundation; either version 2 of the License, or
7 //  (at your option) any later version.
8 //
9 //  This program is distributed in the hope that it will be useful,
10 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
11 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 //  GNU General Public License for more details.
13 //
14 //  You should have received a copy of the GNU General Public License
15 //  along with this program; if not, write to the Free Software
16 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
17 //  ---------------------------------------------------------------------------
18 //  Filter distortion code written by Antti S. Lankila 2007 - 2009.
19 
20 #ifndef __FILTERFP_H__
21 #define __FILTERFP_H__
22 
23 #include <math.h>
24 #include "siddefsfp.h"
25 
26 // ----------------------------------------------------------------------------
27 // The SID filter is modeled with a two-integrator-loop biquadratic filter,
28 // which has been confirmed by Bob Yannes to be the actual circuit used in
29 // the SID chip.
30 //
31 // Measurements show that excellent emulation of the SID filter is achieved,
32 // except when high resonance is combined with high sustain levels.
33 // In this case the SID op-amps are performing less than ideally and are
34 // causing some peculiar behavior of the SID filter. This however seems to
35 // have more effect on the overall amplitude than on the color of the sound.
36 //
37 // The theory for the filter circuit can be found in "Microelectric Circuits"
38 // by Adel S. Sedra and Kenneth C. Smith.
39 // The circuit is modeled based on the explanation found there except that
40 // an additional inverter is used in the feedback from the bandpass output,
41 // allowing the summer op-amp to operate in single-ended mode. This yields
42 // inverted filter outputs with levels independent of Q, which corresponds with
43 // the results obtained from a real SID.
44 //
45 // We have been able to model the summer and the two integrators of the circuit
46 // to form components of an IIR filter.
47 // Vhp is the output of the summer, Vbp is the output of the first integrator,
48 // and Vlp is the output of the second integrator in the filter circuit.
49 //
50 // According to Bob Yannes, the active stages of the SID filter are not really
51 // op-amps. Rather, simple NMOS inverters are used. By biasing an inverter
52 // into its region of quasi-linear operation using a feedback resistor from
53 // input to output, a MOS inverter can be made to act like an op-amp for
54 // small signals centered around the switching threshold.
55 //
56 // Qualified guesses at SID filter schematics are depicted below.
57 //
58 // SID filter
59 // ----------
60 //
61 //     -----------------------------------------------
62 //    |                                               |
63 //    |            ---Rq--                            |
64 //    |           |       |                           |
65 //    |  ------------<A]-----R1---------              |
66 //    | |                               |             |
67 //    | |                        ---C---|      ---C---|
68 //    | |                       |       |     |       |
69 //    |  --R1--    ---R1--      |---Rs--|     |---Rs--|
70 //    |        |  |       |     |       |     |       |
71 //     ----R1--|-----[A>--|--R-----[A>--|--R-----[A>--|
72 //             |          |             |             |
73 // vi -----R1--           |             |             |
74 //
75 //                       vhp           vbp           vlp
76 //
77 //
78 // vi  - input voltage
79 // vhp - highpass output
80 // vbp - bandpass output
81 // vlp - lowpass output
82 // [A> - op-amp
83 // R1  - summer resistor
84 // Rq  - resistor array controlling resonance (4 resistors)
85 // R   - NMOS FET voltage controlled resistor controlling cutoff frequency
86 // Rs  - shunt resitor
87 // C   - capacitor
88 //
89 //
90 //
91 // SID integrator
92 // --------------
93 //
94 //                                   V+
95 //
96 //                                   |
97 //                                   |
98 //                              -----|
99 //                             |     |
100 //                             | ||--
101 //                              -||
102 //                   ---C---     ||->
103 //                  |       |        |
104 //                  |---Rs-----------|---- vo
105 //                  |                |
106 //                  |            ||--
107 // vi ----     -----|------------||
108 //        |   ^     |            ||->
109 //        |___|     |                |
110 //        -----     |                |
111 //          |       |                |
112 //          |---R2--                 |
113 //          |
114 //          R1                       V-
115 //          |
116 //          |
117 //
118 //          Vw
119 //
120 // ----------------------------------------------------------------------------
121 class FilterFP
122 {
123 public:
124   FilterFP();
125 
126   void enable_filter(bool enable);
127   void set_chip_model(chip_model model);
128   void set_distortion_properties(float, float, float);
129   void set_type3_properties(float, float, float, float);
130   void set_type4_properties(float, float);
131   void set_clock_frequency(float);
132   void set_nonlinearity(float);
133 
134   inline
135   float clock(float voice1, float voice2, float voice3,
136 	      float ext_in);
137   void reset();
138 
139   // Write registers.
140   void writeFC_LO(reg8);
141   void writeFC_HI(reg8);
142   void writeRES_FILT(reg8);
143   void writeMODE_VOL(reg8);
144 
145 private:
146   void set_Q();
147   void set_w0();
148   float type3_w0(const float dist);
149   float type4_w0();
150   void calculate_helpers();
151   void nuke_denormals();
152   float waveshaper1(float value);
153   float waveshaper2(float value);
154 
155   // Filter enabled.
156   bool enabled;
157 
158   // 6581/8580 filter model (XXX: we should specialize in separate classes)
159   chip_model model;
160 
161   // Filter cutoff frequency.
162   reg12 fc;
163 
164   // Filter resonance.
165   reg8 res;
166 
167   // Selects which inputs to route through filter.
168   reg8 filt;
169 
170   // Switch voice 3 off.
171   reg8 voice3off;
172 
173   // Highpass, bandpass, and lowpass filter modes.
174   reg8 hp_bp_lp;
175 
176   // Output master volume.
177   reg4 vol;
178   float volf; /* avoid integer-to-float conversion at output */
179 
180   // clock
181   float clock_frequency;
182 
183   /* Distortion params for Type3 */
184   float attenuation, distortion_nonlinearity, intermixing_leaks;
185 
186   /* Type3 params. */
187   float type3_baseresistance, type3_offset, type3_steepness, type3_minimumfetresistance;
188 
189   /* Type4 params */
190   float type4_k, type4_b;
191 
192   // State of filter.
193   float Vhp, Vbp, Vlp;
194 
195   /* Resonance/Distortion/Type3/Type4 helpers. */
196   float type4_w0_cache, _1_div_Q, type3_fc_kink_exp, distortion_CT;
197 
198   float nonlinearity;
199 
200 friend class SIDFP;
201 };
202 
203 // ----------------------------------------------------------------------------
204 // Inline functions.
205 // The following functions are defined inline because they are called every
206 // time a sample is calculated.
207 // ----------------------------------------------------------------------------
208 
209 const float sidcaps_6581 = 470e-12f;
210 const float FC_TO_OSC = 512.f;
211 
212 inline
fastexp(float val)213 static float fastexp(float val) {
214     typedef union {
215         int i;
216         float f;
217     } conv;
218 
219     conv tmp;
220 
221     /* single precision fp has 1 + 8 + 23 bits, exponent bias is 127.
222      * It therefore follows that we need to shift left by 23 bits, and to
223      * calculate exp(x) instead of pow(2, x) we divide the power by ln(2). */
224     const float a = static_cast<float>((1 << 23) / M_LN2);
225     /* The other factor corrects for the exponent bias so that 2^0 = 1. */
226     const float b = static_cast<float>((1 << 23) * 127);
227     /* According to "A Fast, Compact Approximation of the Exponential Function"
228      * by Nicol N. Schraudolph, 60801.48 yields the minimum RMS error for the
229      * piecewise-linear approximation when using doubles (20 bits residual).
230      * We have 23 bits, so we scale this value by 8. */
231     const float c = 60801.48f * 8.f;
232 
233     /* Parenthesis are important: C standard disallows folding subtraction.
234      * Unfortunately GCC appears to generate a write to memory rather than
235      * handle this conversion entirely in registers. */
236     tmp.i = static_cast<int>(a * val + (b - c + 0.5f));
237     return tmp.f;
238 }
239 
240 inline
type3_w0(const float dist)241 float FilterFP::type3_w0(const float dist)
242 {
243     float fetresistance = type3_fc_kink_exp;
244     if (dist > 0) {
245         fetresistance *= fastexp(dist * type3_steepness);
246     }
247     const float dynamic_resistance = type3_minimumfetresistance + fetresistance;
248 
249     /* 2 parallel resistors */
250     const float _1_div_resistance = (type3_baseresistance + dynamic_resistance) / (type3_baseresistance * dynamic_resistance);
251     /* 1.f / (clock * caps * resistance) */
252     return distortion_CT * _1_div_resistance;
253 }
254 
255 inline
type4_w0()256 float FilterFP::type4_w0()
257 {
258     const float freq = type4_k * fc + type4_b;
259     return 2.f * static_cast<float>(M_PI) * freq / clock_frequency;
260 }
261 
waveshaper1(float value)262 inline float FilterFP::waveshaper1(float value)
263 {
264     if (value > distortion_nonlinearity) {
265         value -= (value - distortion_nonlinearity) * 0.5f;
266     }
267     return value;
268 }
269 
270 // ----------------------------------------------------------------------------
271 // SID clocking - 1 cycle.
272 // ----------------------------------------------------------------------------
273 inline
clock(float voice1,float voice2,float voice3,float ext_in)274 float FilterFP::clock(float voice1,
275 		   float voice2,
276 		   float voice3,
277 		   float ext_in)
278 {
279     float Vi = 0.f, Vf = 0.f;
280 
281     // Route voices into or around filter.
282     ((filt & 1) ? Vi : Vf) += voice1;
283     ((filt & 2) ? Vi : Vf) += voice2;
284     // NB! Voice 3 is not silenced by voice3off if it is routed through
285     // the filter.
286     if (filt & 4) {
287 	Vi += voice3;
288     } else if (! voice3off) {
289 	Vf += voice3;
290     }
291     ((filt & 8) ? Vi : Vf) += ext_in;
292 
293     if (hp_bp_lp & 1) {
294 	Vf += Vlp;
295     }
296     if (hp_bp_lp & 2) {
297 	Vf += Vbp;
298     }
299     if (hp_bp_lp & 4) {
300 	Vf += Vhp;
301     }
302 
303     if (model == MOS6581) {
304         Vlp -= Vbp * type3_w0(Vbp);
305         Vbp -= Vhp * type3_w0(Vhp);
306         Vhp = (Vbp * _1_div_Q - Vlp - Vi * 0.85f) * attenuation;
307 
308         /* output strip mixing to filter state */
309         if (hp_bp_lp & 1) {
310             Vlp += Vf * intermixing_leaks;
311         }
312         if (hp_bp_lp & 2) {
313             Vbp += Vf * intermixing_leaks;
314         }
315         if (hp_bp_lp & 4) {
316             Vhp += Vf * intermixing_leaks;
317         }
318 
319         /* saturate. This is likely the output inverter saturation. */
320         Vf *= volf;
321         Vf = waveshaper1(Vf);
322     } else {
323         /* On the 8580, BP appears mixed in phase with the rest. */
324         Vlp += Vbp * type4_w0_cache;
325         Vbp += Vhp * type4_w0_cache;
326         Vhp = -Vbp * _1_div_Q - Vlp - Vi;
327         Vf *= volf;
328     }
329 
330     return Vf;
331 }
332 
333 inline
nuke_denormals()334 void FilterFP::nuke_denormals()
335 {
336     /* We only need this for systems that don't do -msse and -mfpmath=sse */
337     if (Vbp > -1e-12f && Vbp < 1e-12f)
338         Vbp = 0;
339     if (Vlp > -1e-12f && Vlp < 1e-12f)
340         Vlp = 0;
341 }
342 
343 #endif // not __FILTER_H__
344