Hammond style vibrato/chorus via circuit simulation

[MarkT]

I wanted to improve the accuracy of my vibrato/chorus for my Hammond organ MIDI synth code, as its
a distinctive part of the sound. (Note this is separate from the Leslie rotating speaker effect also part
of the sound).

So I found an old Hammond vibrato line-unit on eBay:


These are basically multi-stage LC delay units that connect to a wonderfully old-school rotating scanner
unit (which I won't go into here...) - the various delay taps are effectively interpolated between as the
scanner rotates.

The line-unit taps used can be selected by a complex switch array, as the iconic Hamond B3 vibrato/chorus
unit has 19 taps of which 9 are used at a time.

Injecting square wave signal into the line-unit and looking at a variety of taps, the first few:

and every 5 taps or so:

Shows the progressing delay and the ripples in the response due to the analog filtering of the LCR array.
(Note the line unit I bought isn't an exact match for the B3 unit I went on to emulate, but those are
much harder to source)

So I wrote code to (crudely) emulate the LC array line-unit as part of the chorus/vibrato effect. To get
reasonable convincing results I found up-sampling to 88200 sample rate was needed, and each sample
triggers 19 stage LCR circuit simulation via the differential equations for the C's and L's. Thus there are
probably close to 1000 operations per original sample, and floating point too. This is relying on the
grunt of the T4.0/T4.1 to handle.

Taking several traces from various taps (with the code modified to output 30% of the zeroth tap plus a
specific other tap):

Shows that the simulation is producing similar time-domain signals, and delays. The zeroth tap signal produced an
initial rising edge for triggering the 'scope, this was purely for aligning the traces - this image is actually a composite
of 5 separate traces overlapped by post-processing the images. The code doesn't do this, it passes the various taps
to a simulation of the rotating scanner.

This convinced me the code was doing a plausible job of simulation of the line unit, as well as the sound of it...

Anyway here's the library for those that want to try this out - be warned it takes about 17% of CPU on a T4 at 600MHz
or something like that, so its not for a slow board or one without hardware FP.

HammondVibrato.h:

#ifndef __HAMMONDVIBRATO_H__
#define __HAMMONDVIBRATO_H__
/* Hammond-style vibrato/scanner unit simulator
 * Copyright (c) 2020, Mark Tillotson
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include "Arduino.h"
#include "AudioStream.h"

class AudioEffectHammondVibrato : public AudioStream
{
public:
  AudioEffectHammondVibrato(void) : AudioStream(1, inputQueueArray)
  {
    enabled = true ;
    vibrato = true ;
    depth = 1 ;

    scanner_ph = 0 ;
    scanner_seg = 0 ;
    tap = 0 ;
    ntap = 1 ;
  }

  void enable (bool on) ;        // on/off
  void vibrato_mode () ;         // set vibrato mode
  void chorus_mode () ;          // set chorus mode
  void set_depth (int _depth) ;  // 1, 2 or 3 
  
  virtual void update(void);

  
protected:
  float butterworth (float s) ;
  void simulate_lineunit (float s) ;
  void set_taps ();

  
private:
  int scanner_ph ;  // scanner segment tracking
  int scanner_seg ;
  int tap, ntap ;

  bool enabled ;    // configuration
  bool vibrato ;
  int depth ;

  float filtvals[2] ;  // butterworth coeffs

  float V[19] ;  // circuit state, node voltages/currents
  float I[18] ;
  float Vg, deltaVg ;
  float out[19] ;  // tap outputs

  audio_block_t * inputQueueArray [1] ;
};

#endif

HammondVibrato.cpp:

#include "HammondVibrato.h"
/* Hammond-style vibrato/scanner unit simulator
 * Copyright (c) 2020, Mark Tillotson
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

void AudioEffectHammondVibrato::enable (bool on)        { enabled = on ; }
void AudioEffectHammondVibrato::vibrato_mode ()         { vibrato = true ; }
void AudioEffectHammondVibrato::chorus_mode ()          { vibrato = false ; }

void AudioEffectHammondVibrato::set_depth (int _depth)
{
  if (_depth < 1) _depth = 1 ;
  if (_depth > 3) _depth = 3 ;
  depth = _depth ;
  set_taps () ;
}

#define SAMPS_PER_SEG 403    // int (round (44100.0 / (6.83*16))), samples per scanner segment 


// somewhat smoothed trianglur cross-fade envelope
static const float crossfade[SAMPS_PER_SEG+1] =
{
  0.00000, 0.00005, 0.00021, 0.00047, 0.00084, 0.00131, 0.00189, 0.00257, 0.00336, 0.00426, 0.00525, 0.00636, 0.00757, 0.00888, 0.01030, 0.01182,
  0.01345, 0.01518, 0.01702, 0.01897, 0.02102, 0.02317, 0.02543, 0.02779, 0.03026, 0.03284, 0.03548, 0.03813, 0.04078, 0.04342, 0.04607, 0.04872,
  0.05136, 0.05401, 0.05666, 0.05931, 0.06195, 0.06460, 0.06725, 0.06989, 0.07254, 0.07519, 0.07783, 0.08048, 0.08313, 0.08577, 0.08842, 0.09107,
  0.09371, 0.09636, 0.09901, 0.10165, 0.10430, 0.10695, 0.10959, 0.11224, 0.11489, 0.11754, 0.12018, 0.12283, 0.12548, 0.12812, 0.13077, 0.13342,
  0.13606, 0.13871, 0.14136, 0.14400, 0.14665, 0.14930, 0.15194, 0.15459, 0.15724, 0.15988, 0.16253, 0.16518, 0.16782, 0.17047, 0.17312, 0.17577,
  0.17841, 0.18106, 0.18371, 0.18635, 0.18900, 0.19165, 0.19429, 0.19694, 0.19959, 0.20223, 0.20488, 0.20753, 0.21017, 0.21282, 0.21547, 0.21811,
  0.22076, 0.22341, 0.22605, 0.22870, 0.23135, 0.23400, 0.23664, 0.23929, 0.24194, 0.24458, 0.24723, 0.24988, 0.25252, 0.25517, 0.25782, 0.26046,
  0.26311, 0.26576, 0.26840, 0.27105, 0.27370, 0.27634, 0.27899, 0.28164, 0.28428, 0.28693, 0.28958, 0.29222, 0.29487, 0.29752, 0.30017, 0.30281,
  0.30546, 0.30811, 0.31075, 0.31340, 0.31605, 0.31869, 0.32134, 0.32399, 0.32663, 0.32928, 0.33193, 0.33457, 0.33722, 0.33987, 0.34251, 0.34516,
  0.34781, 0.35045, 0.35310, 0.35575, 0.35840, 0.36104, 0.36369, 0.36634, 0.36898, 0.37163, 0.37428, 0.37692, 0.37957, 0.38222, 0.38486, 0.38751,
  0.39016, 0.39280, 0.39545, 0.39810, 0.40074, 0.40339, 0.40604, 0.40868, 0.41133, 0.41398, 0.41663, 0.41927, 0.42192, 0.42457, 0.42721, 0.42986,
  0.43251, 0.43515, 0.43780, 0.44045, 0.44309, 0.44574, 0.44839, 0.45103, 0.45368, 0.45633, 0.45897, 0.46162, 0.46427, 0.46691, 0.46956, 0.47221,
  0.47486, 0.47750, 0.48015, 0.48280, 0.48544, 0.48809, 0.49074, 0.49338, 0.49603, 0.49868, 0.50132, 0.50397, 0.50662, 0.50926, 0.51191, 0.51456,
  0.51720, 0.51985, 0.52250, 0.52514, 0.52779, 0.53044, 0.53309, 0.53573, 0.53838, 0.54103, 0.54367, 0.54632, 0.54897, 0.55161, 0.55426, 0.55691,
  0.55955, 0.56220, 0.56485, 0.56749, 0.57014, 0.57279, 0.57543, 0.57808, 0.58073, 0.58337, 0.58602, 0.58867, 0.59132, 0.59396, 0.59661, 0.59926,
  0.60190, 0.60455, 0.60720, 0.60984, 0.61249, 0.61514, 0.61778, 0.62043, 0.62308, 0.62572, 0.62837, 0.63102, 0.63366, 0.63631, 0.63896, 0.64160,
  0.64425, 0.64690, 0.64955, 0.65219, 0.65484, 0.65749, 0.66013, 0.66278, 0.66543, 0.66807, 0.67072, 0.67337, 0.67601, 0.67866, 0.68131, 0.68395,
  0.68660, 0.68925, 0.69189, 0.69454, 0.69719, 0.69983, 0.70248, 0.70513, 0.70778, 0.71042, 0.71307, 0.71572, 0.71836, 0.72101, 0.72366, 0.72630,
  0.72895, 0.73160, 0.73424, 0.73689, 0.73954, 0.74218, 0.74483, 0.74748, 0.75012, 0.75277, 0.75542, 0.75806, 0.76071, 0.76336, 0.76600, 0.76865,
  0.77130, 0.77395, 0.77659, 0.77924, 0.78189, 0.78453, 0.78718, 0.78983, 0.79247, 0.79512, 0.79777, 0.80041, 0.80306, 0.80571, 0.80835, 0.81100,
  0.81365, 0.81629, 0.81894, 0.82159, 0.82423, 0.82688, 0.82953, 0.83218, 0.83482, 0.83747, 0.84012, 0.84276, 0.84541, 0.84806, 0.85070, 0.85335,
  0.85600, 0.85864, 0.86129, 0.86394, 0.86658, 0.86923, 0.87188, 0.87452, 0.87717, 0.87982, 0.88246, 0.88511, 0.88776, 0.89041, 0.89305, 0.89570,
  0.89835, 0.90099, 0.90364, 0.90629, 0.90893, 0.91158, 0.91423, 0.91687, 0.91952, 0.92217, 0.92481, 0.92746, 0.93011, 0.93275, 0.93540, 0.93805,
  0.94069, 0.94334, 0.94599, 0.94864, 0.95128, 0.95393, 0.95658, 0.95922, 0.96187, 0.96452, 0.96716, 0.96974, 0.97221, 0.97457, 0.97683, 0.97898,
  0.98103, 0.98298, 0.98482, 0.98655, 0.98818, 0.98970, 0.99112, 0.99243, 0.99364, 0.99475, 0.99574, 0.99664, 0.99743, 0.99811, 0.99869, 0.99916,
  0.99953, 0.99979, 0.99995, 1.00000,
};


void AudioEffectHammondVibrato::update ()
{
  if (!enabled)  // if effect not enabled simply pass-through
  {
    audio_block_t * block = receiveReadOnly (0) ;
    if (block != NULL)
    {
      transmit (block) ;
      release (block) ;
      return ;
    }
  }
  
  audio_block_t * output = allocate () ;
  if (output == NULL) // can't do anything without a buffer
  {
    //Serial.print (".") ;
    return ;
  }

  audio_block_t * input = receiveReadOnly (0) ;
  int16_t * data = input == NULL ? NULL : input->data ;

  for (int i = 0 ; i < AUDIO_BLOCK_SAMPLES ; i++)
  {
    float s = data == NULL ? 0.0 : data[i] ;

    simulate_lineunit (butterworth (2*s)) ;  // interpolate to 88200 SPS, pass to circuit simulator
    simulate_lineunit (butterworth (0)) ;

    // simulate scanner unit with smoothed triangular cross-fade
    float interpolated = crossfade [scanner_ph]                 * out[ntap] +
    			 crossfade [SAMPS_PER_SEG - scanner_ph] * out[tap] ;
    
    output->data[i] = (int16_t) interpolated ;
    
    scanner_ph += 1 ;  // update the scanner segment state
    if (scanner_ph >= SAMPS_PER_SEG)
    {
      scanner_seg += 1 ;
      scanner_seg &= 0xF ;
      scanner_ph = 0 ;
      set_taps () ;
    }
  }
  
  if (input != NULL)
    release (input) ;
  transmit (output) ;
  release (output) ;
}

//////////////////////////////////// line unit circuit simulation //////////////////////////////////

static const float dt = 1.0 / (2 * 44100) ;
static const float Rind = 200.0 ;  // perhaps 500 is less resonant?

// inductances in henries
static const float ind [] =
{
  0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5, 0.5, 0.5,
  0.5, 0.5, 0.5
};

#define nF4 (4.0e-9)
#define nF1 (1.0e-9)

// capacitances in farads
static const float cap[] =
{
  0.0, nF4, nF4, nF4, nF4,
  nF4, nF4, nF4, nF4, nF4,
  nF4, nF4, nF4, nF4, nF4,
  nF4, nF4, nF4, nF1
};

// conductances of the resistor dividers
static const float cond[] =
{
  1/(27e3+68e3),
  1/(56e3+150e3),
  1/(39e3+150e3),
  1/(33e3+180e3),
  1/(18e3+180e3),
  1/(12e3+180e3),
  0, 0, 0, 0, 0,
  0, 0, 0, 0, 0,
  0, 0, 1/15e3
};

// gain factors due to resistive dividers
static const float gain[] =
{
  0.716, 0.728, 0.794, 0.845, 0.909, 0.937, 1, 1, 1, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 1
};

// crude differential equation simulation
void AudioEffectHammondVibrato::simulate_lineunit (float val)
{
  // first stage
  V[0] = val ;
  I[0] += (V[0] - V[1]) * dt / ind[0] ;
  float Iout = (V[0] - Vg) * cond[0] ;
  out[0] = V[0] * gain[0] ;

  for (int node = 1 ; node < 18 ; node++)
  {
    // middle stages
    float vdiff = V[node] - V[node+1] ;
    float vrind = I[node] * Rind ;
    vdiff -= vrind ;
    I[node] += vdiff * dt / ind[node] ;
    float Io = I[node-1] - I[node] ;
    float iR = (V[node] - Vg) * cond[node] ;
    float iC = Io - iR ;
    V[node] += deltaVg ;
    V[node] += iC * dt / cap[node] ;  // need to compensate for delta Vg
    out[node] = V[node] * gain[node] ;
    Iout += Io ;
  }
  // final stage
  float Io = I[17] ;
  float iR = V[18] * cond[18] ;
  float iC = Io - iR ;
  V[18] += iC * dt / cap[18] ;
  out[18] = V[18] * gain[18] ;
  Iout += Io ;
  if (!vibrato)  // for chorus the "ground" side of LC network is floated through a resistor
  {
    deltaVg = Iout * 22e3 - Vg ;
    Vg = Iout * 22e3 ;
  }
  else
  {
    deltaVg = 0 ;
    Vg = 0 ;
  }
}


//////////////////////////////////// tap point setup //////////////////////////////////


static byte tap_points_1[] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 4, 3, 2, 1, 0 };
static byte tap_points_2[] = { 0, 1, 2, 4, 6, 8,10,11,12,11,10, 8, 6, 4, 2, 1, 0 };
static byte tap_points_3[] = { 0, 1, 3, 6, 9,12,15,17,18,17,15,12, 9, 6, 3, 1, 0 };

void AudioEffectHammondVibrato::set_taps ()
{
  switch (depth)
  {
  case 1:
    tap = tap_points_1 [scanner_seg] ;
    ntap = tap_points_1 [scanner_seg+1] ;
    break ;
  case 2:
    tap = tap_points_2 [scanner_seg] ;
    ntap = tap_points_2 [scanner_seg+1] ;
    break ;
  case 3:
    tap = tap_points_3 [scanner_seg] ;
    ntap = tap_points_3 [scanner_seg+1] ;
    break ;
  }
}

//////////////////////////////////// oversampling interpolation filter //////////////////////////////////

#define a0 0.0572
#define a1 -1.21888
#define a2  0.44768

// interpolation filter for up-sampling before circuit sim.
float AudioEffectHammondVibrato::butterworth (float t)
{
  t -= a1 * filtvals[0] ;
  t -= a2 * filtvals[1] ;
  float u = t ;
  t += 2 * filtvals[0] ;
  t += filtvals[1] ;
  filtvals[1] = filtvals[0] ;
  filtvals[0] = u ;
  
  return a0 * t ;
}