You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

369 lines
13KB

  1. // Mutable Instruments Ripples emulation for VCV Rack
  2. // Copyright (C) 2020 Tyler Coy
  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 3 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, see <https://www.gnu.org/licenses/>.
  16. #include <cmath>
  17. #include <algorithm>
  18. #include <random>
  19. #include "rack.hpp"
  20. #include "aafilter.hpp"
  21. using namespace rack;
  22. namespace ripples
  23. {
  24. // Frequency knob
  25. static const float kFreqKnobMin = 20.f;
  26. static const float kFreqKnobMax = 20000.f;
  27. static const float kFreqKnobVoltage =
  28. std::log2f(kFreqKnobMax / kFreqKnobMin);
  29. // Calculate base and multiplier values to pass to configParam so that the
  30. // knob value is labeled in Hz.
  31. // Model the knob as a generic V/oct input with 100k input impedance.
  32. // Assume the internal knob voltage `v` is on the interval [0, vmax] and
  33. // let `p` be the position of the knob varying linearly along [0, 1]. Then,
  34. // freq = fmin * 2^v
  35. // v = vmax * p
  36. // vmax = log2(fmax / fmin)
  37. // freq = fmin * 2^(log2(fmax / fmin) * p)
  38. // = fmin * (fmax / fmin)^p
  39. static const float kFreqKnobDisplayBase = kFreqKnobMax / kFreqKnobMin;
  40. static const float kFreqKnobDisplayMultiplier = kFreqKnobMin;
  41. // Frequency CV amplifier
  42. // The 2164's gain constant is -33mV/dB. Multiply by 6dB/1V to find the
  43. // nominal gain of the amplifier.
  44. static const float kVCAGainConstant = -33e-3f;
  45. static const float kPlus6dB = 20.f * std::log10(2.f);
  46. static const float kFreqAmpGain = kVCAGainConstant * kPlus6dB;
  47. static const float kFreqInputR = 100e3f;
  48. static const float kFreqAmpR = -kFreqAmpGain * kFreqInputR;
  49. static const float kFreqAmpC = 560e-12f;
  50. // Resonance CV amplifier
  51. static const float kResInputR = 22e3f;
  52. static const float kResKnobV = 12.f;
  53. static const float kResKnobR = 62e3f;
  54. static const float kResAmpR = 47e3f;
  55. static const float kResAmpC = 560e-12f;
  56. // Gain CV amplifier
  57. static const float kGainInputR = 27e3f;
  58. static const float kGainNormalV = 12.f;
  59. static const float kGainNormalR = 15e3f;
  60. static const float kGainAmpR = 47e3f;
  61. static const float kGainAmpC = 560e-12f;
  62. // Filter core
  63. static const float kFilterMaxCutoff = kFreqKnobMax;
  64. static const float kFilterCellR = 33e3f;
  65. static const float kFilterCellRC =
  66. 1.f / (2.f * M_PI * kFilterMaxCutoff);
  67. static const float kFilterCellC = kFilterCellRC / kFilterCellR;
  68. static const float kFilterInputR = 100e3f;
  69. static const float kFilterInputGain = kFilterCellR / kFilterInputR;
  70. static const float kFilterCellSelfModulation = 0.01f;
  71. // Filter core feedback path
  72. static const float kFeedbackRt = 22e3f;
  73. static const float kFeedbackRb = 1e3f;
  74. static const float kFeedbackR = kFeedbackRt + kFeedbackRb;
  75. static const float kFeedbackGain = kFeedbackRb / kFeedbackR;
  76. // Filter core feedforward path
  77. static const float kFeedforwardRt = 300e3f;
  78. static const float kFeedforwardRb = 1e3f;
  79. static const float kFeedforwardR = kFeedforwardRt + kFeedforwardRb;
  80. static const float kFeedforwardGain = kFeedforwardRb / kFeedforwardR;
  81. static const float kFeedforwardC = 220e-9f;
  82. // Filter output amplifiers
  83. static const float kLP2Gain = -100e3f / 39e3f;
  84. static const float kLP4Gain = -100e3f / 33e3f;
  85. static const float kBP2Gain = -100e3f / 39e3f;
  86. // VCA
  87. static const float kVCAInputC = 4.7e-6f;
  88. static const float kVCAInputRt = 100e3f;
  89. static const float kVCAInputRb = 1e3f;
  90. static const float kVCAInputR = kVCAInputRt + kVCAInputRb;
  91. static const float kVCAInputGain = kVCAInputRb / kVCAInputR;
  92. static const float kVCAOutputR = 100e3f;
  93. // Voltage-to-current converters
  94. // Saturation voltage at BJT collector
  95. static const float kVtoICollectorVSat = -10.f;
  96. // Opamp saturation voltage
  97. static const float kOpampSatV = 10.6f;
  98. class RipplesEngine
  99. {
  100. public:
  101. struct Frame
  102. {
  103. // Parameters
  104. float res_knob; // 0 to 1 linear
  105. float freq_knob; // 0 to 1 linear
  106. float fm_knob; // -1 to 1 linear
  107. // Inputs
  108. float res_cv;
  109. float freq_cv;
  110. float fm_cv;
  111. float input;
  112. float gain_cv;
  113. bool gain_cv_present;
  114. // Outputs
  115. float bp2;
  116. float lp2;
  117. float lp4;
  118. float lp4vca;
  119. };
  120. RipplesEngine()
  121. {
  122. setSampleRate(1.f);
  123. }
  124. void setSampleRate(float sample_rate)
  125. {
  126. sample_time_ = 1.f / sample_rate;
  127. cell_voltage_ = 0.f;
  128. aa_filter_.Init(sample_rate);
  129. float oversample_rate =
  130. sample_rate * aa_filter_.GetOversamplingFactor();
  131. float freq_cut = 1.f / (2.f * M_PI * kFreqAmpR * kFreqAmpC);
  132. float res_cut = 1.f / (2.f * M_PI * kResAmpR * kResAmpC);
  133. float gain_cut = 1.f / (2.f * M_PI * kGainAmpR * kGainAmpC);
  134. float ff_cut = 1.f / (2.f * M_PI * kFeedforwardR * kFeedforwardC);
  135. auto cutoffs = simd::float_4(ff_cut, freq_cut, res_cut, gain_cut);
  136. rc_filters_.setCutoffFreq(cutoffs / oversample_rate);
  137. float vca_cut = 1.f / (2.f * M_PI * kVCAInputR * kVCAInputC);
  138. vca_hpf_.setCutoffFreq(vca_cut / oversample_rate);
  139. }
  140. void process(Frame& frame)
  141. {
  142. // Calculate equivalent frequency CV
  143. float v_oct = 0.f;
  144. v_oct += (frame.freq_knob - 1.f) * kFreqKnobVoltage;
  145. v_oct += frame.freq_cv;
  146. v_oct += frame.fm_cv * frame.fm_knob;
  147. v_oct = std::min(v_oct, 0.f);
  148. // Calculate resonance control current
  149. float i_reso = VtoIConverter(kResAmpR, frame.res_cv, kResInputR,
  150. frame.res_knob * kResKnobV, kResKnobR);
  151. // Calculate gain control current
  152. float gain_cv = frame.gain_cv;
  153. float gain_input_r = kGainInputR;
  154. if (!frame.gain_cv_present)
  155. {
  156. gain_cv = kGainNormalV;
  157. gain_input_r += kGainNormalR;
  158. }
  159. float i_vca = VtoIConverter(kGainAmpR, gain_cv, gain_input_r);
  160. // Pack and upsample inputs
  161. int oversampling_factor = aa_filter_.GetOversamplingFactor();
  162. float timestep = sample_time_ / oversampling_factor;
  163. // Add noise to input to bootstrap self-oscillation
  164. float input = frame.input + 1e-6 * (random::uniform() - 0.5f);
  165. auto inputs = simd::float_4(input, v_oct, i_reso, i_vca);
  166. inputs *= oversampling_factor;
  167. simd::float_4 outputs;
  168. for (int i = 0; i < oversampling_factor; i++)
  169. {
  170. inputs = aa_filter_.ProcessUp((i == 0) ? inputs : 0.f);
  171. outputs = CoreProcess(inputs, timestep);
  172. outputs = aa_filter_.ProcessDown(outputs);
  173. }
  174. frame.bp2 = outputs[0];
  175. frame.lp2 = outputs[1];
  176. frame.lp4 = outputs[2];
  177. frame.lp4vca = outputs[3];
  178. }
  179. protected:
  180. float sample_time_;
  181. simd::float_4 cell_voltage_;
  182. ripples::AAFilter<simd::float_4> aa_filter_;
  183. dsp::TRCFilter<simd::float_4> rc_filters_;
  184. dsp::TRCFilter<float> vca_hpf_;
  185. // High-rate processing core
  186. // inputs: vector containing (input, v_oct, i_reso, i_vca)
  187. // returns: vector containing (bp2, lp2, lp4, lp4vca)
  188. simd::float_4 CoreProcess(simd::float_4 inputs, float timestep)
  189. {
  190. rc_filters_.process(inputs);
  191. // Lowpass the control signals
  192. simd::float_4 control = rc_filters_.lowpass();
  193. float v_oct = control[1];
  194. float i_reso = control[2];
  195. float i_vca = control[3];
  196. // Highpass the input signal to generate the resonance feedforward
  197. float feedforward = rc_filters_.highpass()[0];
  198. // The 2164's input terminal is a virtual ground, so we can model the
  199. // vca-integrator cell like so:
  200. // ___
  201. // ┌───────────┤___├───────────┐
  202. // │ R │
  203. // │ ┌───┤├───┤
  204. // │ A*ix │ C │
  205. // ___ │ ┌───┐ │ │╲ │
  206. // vin ──┤___├──┤ ┌─┤ → ├───┴─┤-╲____├── vout
  207. // R │ ↓ix │ └───┘ ┌─┤+╱
  208. // ╧ ╧ ╧ │╱
  209. //
  210. // We can see that:
  211. // ix = (vin + vout) / R
  212. // A*ix = -C * dvout/dt
  213. // Thus,
  214. // dvout/dt = -A/(RC) * (vin + vout)
  215. // Calculate -A / RC
  216. simd::float_4 rad_per_s = -std::exp2f(v_oct) / kFilterCellRC;
  217. // Emulate the filter core
  218. cell_voltage_ = StepRK2(timestep, cell_voltage_, [&](simd::float_4 vout)
  219. {
  220. // vout contains the initial cell voltages (v0, v1 v2, v3)
  221. // Rotate cell voltages. vin will contain (v3, v0, v1, v2)
  222. simd::float_4 vin =
  223. _mm_shuffle_ps(vout.v, vout.v, _MM_SHUFFLE(2, 1, 0, 3));
  224. // The core input is the filter input plus the resonance signal
  225. float vp = feedforward * kFeedforwardGain;
  226. float vn = vout[3] * kFeedbackGain;
  227. float res = kFilterCellR * OTAVCA(vp, vn, i_reso);
  228. simd::float_4 in = inputs[0] * kFilterInputGain + res;
  229. // Replace lowest element of vin with lowest element from in
  230. vin = _mm_move_ss(vin.v, in.v);
  231. // Now, vin contains (in, v0, v1, v2)
  232. // and vout contains (v0, v1, v2, v3)
  233. // Their sum gives us vin + vout for each cell
  234. simd::float_4 vsum = vin + vout;
  235. simd::float_4 dvout = rad_per_s * vsum;
  236. // Generate some even-order harmonics via self-modulation
  237. dvout *= (1.f + vsum * kFilterCellSelfModulation);
  238. return dvout;
  239. });
  240. cell_voltage_ = simd::clamp(cell_voltage_, -kOpampSatV, kOpampSatV);
  241. float lp1 = cell_voltage_[0];
  242. float lp2 = cell_voltage_[1];
  243. float lp4 = cell_voltage_[3];
  244. float bp2 = (lp1 + lp2) * kBP2Gain;
  245. vca_hpf_.process(lp4);
  246. float lp4vca = vca_hpf_.highpass();
  247. lp4vca = -kVCAOutputR * OTAVCA(0.f, lp4vca * kVCAInputGain, i_vca);
  248. lp2 *= kLP2Gain;
  249. lp4 *= kLP4Gain;
  250. return simd::float_4(bp2, lp2, lp4, lp4vca);
  251. }
  252. // Solves an ODE system using the 2nd order Runge-Kutta method
  253. template <typename T, typename F>
  254. T StepRK2(float dt, T y, F f)
  255. {
  256. T k1 = f(y);
  257. T k2 = f(y + k1 * dt / 2.f);
  258. return y + dt * k2;
  259. }
  260. // Model of Ripples nonlinear CV voltage-to-current converters
  261. float VtoIConverter(
  262. float rfb, // Amplifier feedback resistor
  263. float vc, float rc, // CV voltage and input resistor
  264. float vp = 0.f, float rp = 1e12f) // Knob voltage and resistor
  265. {
  266. // Find nominal voltage at the BJT collector, ignoring nonlinearity
  267. float vnom = -(vc * rfb / rc + vp * rfb / rp);
  268. // Apply clipping - naive for now
  269. float vout = std::max(vnom, kVtoICollectorVSat);
  270. // Find voltage at the opamp's negative terminal
  271. float nrc = rp * rfb;
  272. float nrp = rc * rfb;
  273. float nrfb = rc * rp;
  274. float vneg = (vc * nrc + vp * nrp + vout * nrfb) / (nrc + nrp + nrfb);
  275. // Find output current
  276. float iout = (vneg - vout) / rfb;
  277. return std::max(iout, 0.f);
  278. }
  279. // Model of LM13700 OTA VCA, neglecting linearizing diodes
  280. // vp: voltage at positive input terminal
  281. // vn: voltage at negative input terminal
  282. // i_abc: amplifier bias current
  283. // returns: OTA output current
  284. template <typename T>
  285. T OTAVCA(T vp, T vn, T i_abc)
  286. {
  287. // For the derivation of this equation, see this fantastic paper:
  288. // http://www.openmusiclabs.com/files/otadist.pdf
  289. // Thanks guest!
  290. //
  291. // i_out = i_abc * (e^(vi/vt) - 1) / (e^(vi/vt) + 1)
  292. // or equivalently,
  293. // i_out = i_abc * tanh(vi / (2vt))
  294. const float kTemperature = 40.f; // Silicon temperature in Celsius
  295. const float kKoverQ = 8.617333262145e-5;
  296. const float kKelvin = 273.15f; // 0C in K
  297. const float kVt = kKoverQ * (kTemperature + kKelvin);
  298. const float kZlim = 2.f * std::sqrt(3.f);
  299. T vi = vp - vn;
  300. T zlim = kZlim;
  301. T z = math::clamp(vi / (2 * kVt), -zlim, zlim);
  302. // Pade approximant of tanh(z)
  303. T z2 = z * z;
  304. T q = 12.f + z2;
  305. T p = 12.f * z * q / (36.f * z2 + q * q);
  306. return i_abc * p;
  307. }
  308. };
  309. }