Trying spline interpolation and simple gauss filter.

config.h

@@ -61,9 +61,9 @@
 #define REDUCE_LOUDNESS 0
 #define USE_XFADE_DATA 1
 // CHORUS parameters
-#define INTERPOLATION_WINDOW_SIZE 13 // use only odd numbers!!!
-//#define INTERPOLATE QUADRATIC // LINEAR QUADRATIC COSINE CUBIC LAGRANGE 
-#define CHORUS_WAVEFORM WAVEFORM_SINE // WAVEFORM_SINE WAVEFORM_SAWTOOTH WAVEFORM_SAWTOOTH_REVERSE WAVEFORM_SQUARE WAVEFORM_TRIANGLE
+#define INTERPOLATION_WINDOW_SIZE 7 // use only odd numbers!!!
+#define INTERPOLATE Catmull
+#define CHORUS_WAVEFORM WAVEFORM_SINE // WAVEFORM_SINE WAVEFORM_TRIANGLE WAVEFORM_SAWTOOTH WAVEFORM_SAWTOOTH_REVERSE
 #define CHORUS_DELAY_LENGTH_SAMPLES (14*AUDIO_BLOCK_SAMPLES) // one AUDIO_BLOCK_SAMPLES = 2.902ms; you need doubled length, e.g. delay point is 20ms, so you need up to 40ms delay!
 
 //*************************************************************************************************

effect_modulated_delay.cpp

@@ -25,7 +25,7 @@
 #include <Audio.h>
 #include "limits.h"
 #include "effect_modulated_delay.h"
-#include "interpolation.h"
+#include "spline.h"
 #include "config.h"
 
 /******************************************************************/
@@ -57,12 +57,8 @@ boolean AudioEffectModulatedDelay::begin(short *delayline, int d_length)
 
   _delayline = delayline;
   _delay_length = _max_delay_length = d_length;
-  _delay_length_half = _delay_length / 2;
 
   memset(_delayline, 0, sizeof(int16_t)*_delay_length);
-#ifdef INTERPOLATE
-  modulation_interpolate = new interpolation();
-#endif
 
   return (true);
 }
@@ -86,9 +82,7 @@ void AudioEffectModulatedDelay::update(void)
   int8_t j;
   float x[INTERPOLATION_WINDOW_SIZE];
   float y[INTERPOLATION_WINDOW_SIZE];
-  modulation_interpolate->valuelenXY(INTERPOLATION_WINDOW_SIZE);
-  modulation_interpolate->valueX(x);
-  modulation_interpolate->valueY(y);
+  Spline s(x, y, INTERPOLATION_WINDOW_SIZE, INTERPOLATE);
 #endif
 
   bp = block->data;
@@ -104,46 +98,30 @@ void AudioEffectModulatedDelay::update(void)
       _delayline[_circ_idx] = *bp;
 
       // calculate modulation index
-      mod_idx = float(*mp) / SHRT_MAX * _delay_length_half; // calculate an index with modulation as a float(!!!)
+      mod_idx = float(*mp) / SHRT_MAX * float(_delay_length / 2); // calculate an index with modulation as a float(!!!)
 
 #ifdef INTERPOLATE
       // get x/y values around mod_idx
       uint8_t c = 0;
       int16_t c_mod_idx = int(mod_idx + 0.5) + _circ_idx;
+      int32_t avg=0;
       for (j = INTERPOLATION_WINDOW_SIZE / -2; j <= INTERPOLATION_WINDOW_SIZE / 2; j++)
       {
-        int16_t jc_mod_idx = (c_mod_idx + j) % _delay_length;
+        int16_t jc_mod_idx = (c_mod_idx + j) % _delay_length - 1;
         if (jc_mod_idx < 0)
-          y[c] = float(_delayline[_delay_length + jc_mod_idx]);
+          y[c] = float(_delayline[_delay_length - 1 + jc_mod_idx]);
         else
           y[c] = float(_delayline[jc_mod_idx]);
         x[c] = float(c);
+        avg += y[c];
         c++; // because 42 is the answer! ;-)
       }
 
-      modulation_interpolate->valueI(mod_idx - int(mod_idx + 0.5));
-
-#if INTERPOLATE == LINEAR
-      *bp = int(modulation_interpolate->LinearInterpolate() + 0.5);
-#elif INTERPOLATE == QUDRATIC
-      *bp = int(modulation_interpolate->QuadraticInterpolate() + 0.5);
-#elif INTERPOLATE == COSINE
-      *bp = int(modulation_interpolate->CosineInterpolate() + 0.5);
-#elif INTERPOLATE == CUBIC
-      *bp = int(modulation_interpolate->CubicInterpolate() + 0.5);
-#elif INTERPOLATE == LAGRANGE
-      *bp = int(modulation_interpolate->LagrangeInterpolate() + 0.5);
-#else
-      // No interpolation - should sound really bad...
-      int16_t c_mod_idx = (int(mod_idx + 0.5) + _circ_idx) % _delay_length;
-      if (c_mod_idx < 0)
-        *bp = _delayline[_delay_length - 1 + c_mod_idx];
-      else
-        *bp = _delayline[c_mod_idx];
-#endif
+      //*bp = int(s.value(mod_idx - int(mod_idx + 0.5)) + 0.5);
+      *bp = avg / INTERPOLATION_WINDOW_SIZE;
 #else
       // No interpolation - should sound really bad...
-      int16_t c_mod_idx = (int(mod_idx + 0.5) + _circ_idx) % _delay_length;
+      int16_t c_mod_idx = (int(mod_idx + 0.5) + _circ_idx) % _delay_length - 1;
       if (c_mod_idx < 0)
         *bp = _delayline[_delay_length - 1 + c_mod_idx];
       else
@@ -169,5 +147,4 @@ void AudioEffectModulatedDelay::update(void)
 void AudioEffectModulatedDelay::setDelay(float milliseconds)
 {
   _delay_length = min(AUDIO_SAMPLE_RATE * milliseconds / 500, _max_delay_length);
-  _delay_length_half = _delay_length / 2;
 }

effect_modulated_delay.h

@@ -50,17 +50,13 @@ class AudioEffectModulatedDelay :
     virtual void setDelay(float milliseconds);
 
   private:
+    void _spline_interpolation(float* x, float* a, uint8_t n);
+    int16_t interpolate(int16_t x1, int16_t y1, int16_t x2, int16_t y2);
     audio_block_t *inputQueueArray[2];
-
-#ifdef INTERPOLATE
-    class interpolation *modulation_interpolate;
-#endif
-
     int16_t *_delayline;
     int16_t _circ_idx;
     uint16_t _max_delay_length;
     uint16_t _delay_length;
-    uint16_t _delay_length_half;
 };
 
 #endif

spline.cpp

@@ -0,0 +1,110 @@
+#include "Arduino.h"
+#include "spline.h"
+#include <math.h>
+
+Spline::Spline(void) {
+  _prev_point = 0;
+}
+
+Spline::Spline( float x[], float y[], int numPoints, int degree )
+{
+  setPoints(x,y,numPoints);
+  setDegree(degree);
+  _prev_point = 0;
+}
+
+Spline::Spline( float x[], float y[], float m[], int numPoints )
+{
+  setPoints(x,y,m,numPoints);
+  setDegree(Hermite);
+  _prev_point = 0;
+}
+
+void Spline::setPoints( float x[], float y[], int numPoints ) {
+  _x = x;
+  _y = y;
+  _length = numPoints;
+}
+
+void Spline::setPoints( float x[], float y[], float m[], int numPoints ) {
+  _x = x;
+  _y = y;
+  _m = m;
+  _length = numPoints;
+}
+
+void Spline::setDegree( int degree ){
+  _degree = degree;
+}
+
+float Spline::value( float x )
+{
+  if( _x[0] > x ) { 
+    return _y[0]; 
+  }
+  else if ( _x[_length-1] < x ) { 
+    return _y[_length-1]; 
+  }
+  else {
+    for(int i = 0; i < _length; i++ )
+    {
+      int index = ( i + _prev_point ) % _length;
+      
+      if( _x[index] == x ) {
+        _prev_point = index;
+        return _y[index];
+      } else if( (_x[index] < x) && (x < _x[index+1]) ) {
+        _prev_point = index;
+        return calc( x, index );
+      }
+    }    
+  }
+}
+
+float Spline::calc( float x, int i )
+{
+  switch( _degree ) {
+    case 0:
+      return _y[i];
+    case 1:
+      if( _x[i] == _x[i+1] ) {
+        // Avoids division by 0
+        return _y[i];
+      } else {
+        return _y[i] + (_y[i+1] - _y[i]) * ( x - _x[i]) / ( _x[i+1] - _x[i] );      
+      }
+    case Hermite:
+      return hermite( ((x-_x[i]) / (_x[i+1]-_x[i])), _y[i], _y[i+1], _m[i], _m[i+1], _x[i], _x[i+1] );
+    case Catmull:
+      if( i == 0 ) {
+        // x prior to spline start - first point used to determine tangent
+        return _y[1];
+      } else if( i == _length-2 ) {
+        // x after spline end - last point used to determine tangent
+        return _y[_length-2];
+      } else {
+        float t = (x-_x[i]) / (_x[i+1]-_x[i]);
+        float m0 = (i==0 ? 0 : catmull_tangent(i) );
+        float m1 = (i==_length-1 ? 0 : catmull_tangent(i+1) );
+        return hermite( t, _y[i], _y[i+1], m0, m1, _x[i], _x[i+1]);        
+      }
+  }
+}
+
+float Spline::hermite( float t, float p0, float p1, float m0, float m1, float x0, float x1 ) {
+  return (hermite_00(t)*p0) + (hermite_10(t)*(x1-x0)*m0) + (hermite_01(t)*p1) + (hermite_11(t)*(x1-x0)*m1);
+}
+float Spline::hermite_00( float t ) { return (2*pow(t,3)) - (3*pow(t,2)) + 1;}
+float Spline::hermite_10( float t ) { return pow(t,3) - (2*pow(t,2)) + t; }
+float Spline::hermite_01( float t ) { return (3*pow(t,2)) - (2*pow(t,3)); }
+float Spline::hermite_11( float t ) { return pow(t,3) - pow(t,2); }
+
+float Spline::catmull_tangent( int i ) 
+{ 
+  if( _x[i+1] == _x[i-1] ) {
+    // Avoids division by 0
+    return 0;
+  } else {  
+    return (_y[i+1] - _y[i-1]) / (_x[i+1] - _x[i-1]);    
+  } 
+}

spline.h

@@ -0,0 +1,42 @@
+/*
+  Library for 1-d splines
+  Copyright Ryan Michael
+  Licensed under the GPLv3 
+*/
+#ifndef spline_h
+#define spline_h
+
+#include "Arduino.h"
+
+#define Hermite 10
+#define Catmull 11
+
+class Spline
+{
+  public:
+    Spline( void );
+    Spline( float x[], float y[], int numPoints, int degree = 1 );
+    Spline( float x[], float y[], float m[], int numPoints );
+    float value( float x );
+    void setPoints( float x[], float y[], int numPoints );
+    void setPoints( float x[], float y[], float m[], int numPoints );
+    void setDegree( int degree );
+    
+  private:
+    float calc( float, int);
+    float* _x;
+    float* _y;
+    float* _m;
+    int _degree;
+    int _length;
+    int _prev_point;
+    
+    float hermite( float t, float p0, float p1, float m0, float m1, float x0, float x1 );
+    float hermite_00( float t );
+    float hermite_10( float t );
+    float hermite_01( float t );
+    float hermite_11( float t );
+    float catmull_tangent( int i );
+};
+
+#endif