ref: 08efd6ad79ec698103bcfc7567e036b4d1d9ef10
parent: 6e63ff953faee4cf27d626eb0d7dcc5b1410fb91
author: Matthew Wang <mjw7@princeton.edu>
date: Mon Apr 6 12:03:51 EDT 2020
a start on minblep oscs
--- a/LEAF/Inc/leaf-oscillators.h
+++ b/LEAF/Inc/leaf-oscillators.h
@@ -15,7 +15,7 @@
#include "leaf-math.h"
#include "leaf-mempool.h"
-#include "leaf-tables.h"
+#include "leaf-tables.h"
#include "leaf-filters.h"
/*!
@@ -302,6 +302,7 @@
{
float phase;
float inc,freq;
+ uint8_t phaseDidReset;
} _tPhasor;
@@ -599,6 +600,53 @@
/*! @} */
//==============================================================================
+
+ typedef struct _tMinBLEP
+ {
+ float coeffs[6];
+
+ // FilterState
+ float* x1, x2, y1, y2;
+
+ float ratio, lastRatio;
+
+ float overSamplingRatio;
+ int zeroCrossings;
+
+ float lastValue;
+ float lastDelta; // previous derivative ...
+
+ // Tweaking the Blep F
+ float proportionalBlepFreq;
+ uint8_t returnDerivative; // set this to return the FIRST DERIVATIVE of the blep (for first der. discontinuities)
+
+ int blepIndex;
+ int numActiveBleps;
+ //currentActiveBlepOffsets;
+ float* offset;
+ float* freqMultiple;
+ float* pos_change_magnitude;
+ float* vel_change_magnitude;
+
+ int minBlepSize;
+ float* minBlepArray;
+ float* minBlepDerivArray;
+ } _tMinBLEP;
+
+ typedef _tMinBLEP* tMinBLEP;
+
+ void tMinBLEP_init (tMinBLEP* const minblep);
+ void tMinBLEP_free (tMinBLEP* const minblep);
+ void tMinBLEP_initToPool (tMinBLEP* const minblep, tMempool* const pool);
+ void tMinBLEP_freeFromPool (tMinBLEP* const minblep, tMempool* const pool);
+
+ // pass in audio, identify discontinuities and return the audio with minbleps inserted
+ float tMinBLEP_tick (tMinBLEP* const minblep, float input);
+ void tMinBLEP_buildBLEP (tMinBLEP* const minblep);
+ void tMinBLEP_addBLEP (tMinBLEP* const minblep, float offset, float posChange, float velChange);
+
+ void tMinBLEP_setOversamplingRatio (tMinBLEP* const minblep, float ratio);
+ void tMinBLEP_setNumZeroCrossings (tMinBLEP* const minblep, int numCrossings);
#ifdef __cplusplus
}
--- a/LEAF/Src/leaf-oscillators.c
+++ b/LEAF/Src/leaf-oscillators.c
@@ -343,6 +343,7 @@
p->phase = 0.0f;
p->inc = 0.0f;
+ p->phaseDidReset = 0;
}
void tPhasor_freeFromPool(tPhasor* const ph, tMempool* const mp)
@@ -370,8 +371,14 @@
p->phase += p->inc;
- if (p->phase >= 1.0f) p->phase -= 1.0f;
+ p->phaseDidReset = 0;
+ if (p->phase >= 1.0f)
+ {
+ p->phaseDidReset = 1;
+ p->phase -= 1.0f;
+ }
+
return p->phase;
}
@@ -649,4 +656,419 @@
{
_tNeuron* n = *nr;
n->current = current;
+}
+
+
+
+void tMinBLEP_init (tMinBLEP* const minblep)
+{
+ tMinBLEP_initToPool(minblep, &leaf.mempool);
+}
+
+void tMinBLEP_free (tMinBLEP* const minblep)
+{
+ tMinBLEP_freeFromPool(minblep, &leaf.mempool);
+}
+
+void tMinBLEP_initToPool (tMinBLEP* const minblep, tMempool* const mp)
+{
+ _tMempool* m = *mp;
+ _tMinBLEP* mb = *minblep = (_tMinBLEP*) mpool_alloc(sizeof(_tMinBLEP), m);
+
+ mb->overSamplingRatio = 16;
+ mb->zeroCrossings = 16;
+ mb->returnDerivative = 0;
+ mb->proportionalBlepFreq = 0.5; // defaults to NyQuist ....
+
+ mb->lastValue = 0;
+ mb->lastDelta = 0;
+
+// float* x1, x2, y1, y2;
+ // AA FILTER
+// zeromem (coefficients, sizeof (coefficients));
+
+ mb->minBlepSize = (mb->zeroCrossings * 2 * mb->overSamplingRatio) + 1;
+
+ mb->ratio = 1;
+ mb->lastRatio = 1;
+
+// createLowPass (ratio);
+// resetFilters();
+
+ mb->blepIndex = 0;
+ mb->numActiveBleps = 0;
+ //currentActiveBlepOffsets;
+ mb->offset = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+ mb->freqMultiple = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+ mb->pos_change_magnitude = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+ mb->vel_change_magnitude = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+
+ mb->minBlepArray = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+ mb->minBlepDerivArray = (float*) mpool_alloc(sizeof(float) * mb->minBlepSize, m);
+
+ tMinBLEP_buildBLEP(minblep);
+}
+
+void tMinBLEP_freeFromPool (tMinBLEP* const minblep, tMempool* const mp)
+{
+ _tMempool* m = *mp;
+ _tMinBLEP* mb = *minblep;
+
+
+ mpool_free(mb, m);
+}
+
+// SINC Function
+static float SINC(float x)
+{
+ float pix;
+
+ if(x == 0.0f)
+ return 1.0f;
+ else
+ {
+ pix = PI * x;
+ return sinf(pix) / pix;
+ }
+}
+
+// Generate Blackman Window
+static void BlackmanWindow(int n, float *w)
+{
+ int m = n - 1;
+ int i;
+ float f1, f2, fm;
+
+ fm = (float)m;
+ for(i = 0; i <= m; i++)
+ {
+ f1 = (2.0f * PI * (float)i) / fm;
+ f2 = 2.0f * f1;
+ w[i] = 0.42f - (0.5f * cosf(f1)) + (0.08f * cosf(f2));
+ }
+}
+
+// Discrete Fourier Transform
+static void DFT(int n, float *realTime, float *imagTime, float *realFreq, float *imagFreq)
+{
+ int k, i;
+ float sr, si, p;
+
+ for(k = 0; k < n; k++)
+ {
+ realFreq[k] = 0.0f;
+ imagFreq[k] = 0.0f;
+ }
+
+ for(k = 0; k < n; k++)
+ for(i = 0; i < n; i++)
+ {
+ p = (2.0f * PI * (float)(k * i)) / n;
+ sr = cosf(p);
+ si = -sinf(p);
+ realFreq[k] += (realTime[i] * sr) - (imagTime[i] * si);
+ imagFreq[k] += (realTime[i] * si) + (imagTime[i] * sr);
+ }
+}
+
+// Inverse Discrete Fourier Transform
+static void InverseDFT(int n, float *realTime, float *imagTime, float *realFreq, float *imagFreq)
+{
+ int k, i;
+ float sr, si, p;
+
+ for(k = 0; k < n; k++)
+ {
+ realTime[k] = 0.0f;
+ imagTime[k] = 0.0f;
+ }
+
+ for(k = 0; k < n; k++)
+ {
+ for(i = 0; i < n; i++)
+ {
+ p = (2.0f * PI * (float)(k * i)) / n;
+ sr = cosf(p);
+ si = -sinf(p);
+ realTime[k] += (realFreq[i] * sr) + (imagFreq[i] * si);
+ imagTime[k] += (realFreq[i] * si) - (imagFreq[i] * sr);
+ }
+ realTime[k] /= n;
+ imagTime[k] /= n;
+ }
+}
+
+// Complex Absolute Value
+static float complexabs(float x, float y)
+{
+ return sqrtf((x * x) + (y * y));
+}
+
+// Complex Exponential
+static void complexexp(float x, float y, float *zx, float *zy)
+{
+ float expx;
+
+ expx = expf(x);
+ *zx = expx * cosf(y);
+ *zy = expx * sinf(y);
+}
+
+// Compute Real Cepstrum Of Signal
+static void RealCepstrum(int n, float *signal, float *realCepstrum)
+{
+ int i;
+
+ float realTime[n];
+ float imagTime[n];
+ float realFreq[n];
+ float imagFreq[n];
+
+ // Compose Complex FFT Input
+
+ for(i = 0; i < n; i++)
+ {
+ realTime[i] = signal[i];
+ imagTime[i] = 0.0f;
+ }
+
+ // Perform DFT
+
+ DFT(n, realTime, imagTime, realFreq, imagFreq);
+
+ // Calculate Log Of Absolute Value
+
+ for(i = 0; i < n; i++)
+ {
+ realFreq[i] = logf(complexabs(realFreq[i], imagFreq[i]));
+ imagFreq[i] = 0.0f;
+ }
+
+ // Perform Inverse FFT
+
+ InverseDFT(n, realTime, imagTime, realFreq, imagFreq);
+
+ // Output Real Part Of FFT
+ for(i = 0; i < n; i++)
+ realCepstrum[i] = realTime[i];
+}
+
+// Compute Minimum Phase Reconstruction Of Signal
+static void MinimumPhase(int n, float *realCepstrum, float *minimumPhase)
+{
+ int i, nd2;
+ float realTime[n];
+ float imagTime[n];
+ float realFreq[n];
+ float imagFreq[n];
+
+ nd2 = n / 2;
+
+ if((n % 2) == 1)
+ {
+ realTime[0] = realCepstrum[0];
+ for(i = 1; i < nd2; i++)
+ realTime[i] = 2.0f * realCepstrum[i];
+ for(i = nd2; i < n; i++)
+ realTime[i] = 0.0f;
+ }
+ else
+ {
+ realTime[0] = realCepstrum[0];
+ for(i = 1; i < nd2; i++)
+ realTime[i] = 2.0f * realCepstrum[i];
+ realTime[nd2] = realCepstrum[nd2];
+ for(i = nd2 + 1; i < n; i++)
+ realTime[i] = 0.0f;
+ }
+
+ for(i = 0; i < n; i++)
+ imagTime[i] = 0.0f;
+
+ DFT(n, realTime, imagTime, realFreq, imagFreq);
+
+ for(i = 0; i < n; i++)
+ complexexp(realFreq[i], imagFreq[i], &realFreq[i], &imagFreq[i]);
+
+ InverseDFT(n, realTime, imagTime, realFreq, imagFreq);
+
+ for(i = 0; i < n; i++)
+ minimumPhase[i] = realTime[i];
+}
+
+
+// Generate the MinBLEP
+void tMinBLEP_buildBLEP(tMinBLEP* const minblep)
+{
+ _tMinBLEP* m = *minblep;
+
+ int i, n;
+ float r, a, b;
+
+ n = m->minBlepSize;
+
+ // Generate Sinc
+
+ a = -m->overSamplingRatio;
+ b = m->overSamplingRatio;
+ for(i = 0; i < n; i++)
+ {
+ r = ((float)i) / ((float)(n - 1));
+ m->minBlepArray[i] = SINC(a + (r * (b - a)));
+ m->minBlepDerivArray[i] = 0;
+ }
+
+ // Window Sinc
+
+ BlackmanWindow(n, m->minBlepDerivArray);
+ for(i = 0; i < n; i++)
+ m->minBlepArray[i] *= m->minBlepDerivArray[i];
+
+ // Minimum Phase Reconstruction
+
+ RealCepstrum(n, m->minBlepArray, m->minBlepDerivArray);
+ MinimumPhase(n, m->minBlepDerivArray, m->minBlepArray);
+
+ // Integrate Into MinBLEP
+ a = 0.0f;
+ float secondInt = 0.0f;
+ for(i = 0; i < n; i++)
+ {
+ a += m->minBlepArray[i];
+ m->minBlepArray[i] = a;
+
+ secondInt += a;
+ m->minBlepDerivArray[i] = secondInt;
+ }
+
+ // Normalize
+ a = m->minBlepArray[n - 1];
+ a = 1.0f / a;
+ b = 0.0f;
+ for(i = 0; i < n; i++)
+ {
+ m->minBlepArray[i] *= a;
+ b = fmaxf(b, m->minBlepDerivArray[i]);
+ }
+
+ // Normalize ...
+ b = 1.0f/b;
+ for (i = 0; i < n; i++)
+ {
+ m->minBlepDerivArray[i] *= b;
+ m->minBlepDerivArray[i] -= ((float)i) / ((float) n-1);
+
+ // SUBTRACT 1 and invert so the signal (so it goes 1->0)
+ m->minBlepArray[i] -= 1.0f;
+ m->minBlepArray[i] = -m->minBlepArray[i];
+ }
+}
+
+void tMinBLEP_addBLEP (tMinBLEP* const minblep, float offset, float posChange, float velChange)
+{
+ _tMinBLEP* m = *minblep;
+
+ int n = m->minBlepSize;
+
+ m->offset[m->blepIndex] = offset;
+ m->freqMultiple[m->blepIndex] = m->overSamplingRatio*m->proportionalBlepFreq;
+ m->pos_change_magnitude[m->blepIndex] = posChange;
+ m->vel_change_magnitude[m->blepIndex] = velChange;
+
+ m->blepIndex++;
+ if (m->blepIndex >= n) m->blepIndex = 0;
+
+ m->numActiveBleps++;
+}
+
+float tMinBLEP_tick (tMinBLEP* const minblep, float input)
+{
+ _tMinBLEP* m = *minblep;
+ // PROCESS ALL BLEPS -
+ /// for each offset, copy a portion of the blep array to the output ....
+
+ float sample = input;
+
+ int n = m->minBlepSize;
+
+ for (int blep = 1; blep <= m->numActiveBleps; blep++)
+ {
+ int i = (m->blepIndex - blep + n) % n;
+ float adjusted_Freq = m->freqMultiple[i];
+ float exactPosition = m->offset[i];
+
+ double blepPosExact = adjusted_Freq*(exactPosition + 1); // +1 because this needs to trigger on the LOW SAMPLE
+ double blepPosSample = 0;
+ double fraction = modf(blepPosExact, &blepPosSample);
+
+ // LIMIT the scaling on the derivative array
+ // otherwise, it can get TOO large
+ double depthLimited = m->proportionalBlepFreq; //jlimit<double>(.1, 1, proportionalBlepFreq);
+ double blepDeriv_PosExact = depthLimited*m->overSamplingRatio*(exactPosition + 1);
+ double blepDeriv_Sample = 0;
+ double fraction_Deriv = modf(blepDeriv_PosExact, &blepDeriv_Sample);
+
+
+ // DONE ... we reached the end ...
+ if (((int) blepPosExact > n) && ((int) blepDeriv_PosExact > n))
+ break;
+
+ // BLEP has not yet occurred ...
+ if (blepPosExact < 0)
+ continue;
+
+
+ // 0TH ORDER COMPENSATION ::::
+ /// add the BLEP to compensate for discontinuties in the POSITION
+ if ( fabs(m->pos_change_magnitude[i]) > 0 && blepPosSample < n)
+ {
+ // LINEAR INTERPOLATION ::::
+ float lowValue = m->minBlepArray[(int) blepPosSample];
+ float hiValue = lowValue;
+
+ if ((int) blepPosSample + 1 < n)
+ hiValue = m->minBlepArray[(int) blepPosSample + 1];
+
+ float delta = hiValue - lowValue;
+ float exactValue = lowValue + fraction*delta;
+
+ // SCALE by the discontinuity magnitude
+ exactValue *= m->pos_change_magnitude[i];
+
+ // ADD to the thruput
+ sample += exactValue;
+ }
+
+
+ // 1ST ORDER COMPENSATION ::::
+ /// add the BLEP DERIVATIVE to compensate for discontinuties in the VELOCITY
+ if ( fabs(m->vel_change_magnitude[i]) > 0 && blepDeriv_PosExact < n)
+ {
+
+ // LINEAR INTERPOLATION ::::
+ double lowValue = m->minBlepDerivArray[(int) blepDeriv_PosExact];
+ double hiValue = lowValue;
+
+ if ((int) blepDeriv_PosExact + 1 < n)
+ hiValue = m->minBlepDerivArray[(int) blepDeriv_PosExact + 1];
+
+ double delta = hiValue - lowValue;
+ double exactValue = lowValue + fraction_Deriv*delta;
+
+ // SCALE by the discontinuity magnitude`
+ exactValue *= m->vel_change_magnitude[i];
+
+ // ADD to the thruput
+ sample += exactValue;
+
+ }
+
+ // UPDATE ::::
+ m->offset[i] += 1.0f;
+ if (m->offset[i] * adjusted_Freq > n)
+ {
+ m->numActiveBleps--;
+ }
+ }
+ return sample;
}
--- a/LEAF_JUCEPlugin/Source/MyTest.cpp
+++ b/LEAF_JUCEPlugin/Source/MyTest.cpp
@@ -29,6 +29,10 @@
tTriangle tri;
+tMinBLEP minblep;
+
+tPhasor phasor;
+
float gain;
float freq;
float dtime;
@@ -47,6 +51,9 @@
tTriangle_init(&tri);
+ tMinBLEP_init(&minblep);
+
+ tPhasor_init(&phasor);
// tNoise_init(&noise, WhiteNoise);
//
// tAutotune_init(&at, 1, 1024, 512);
@@ -66,6 +73,15 @@
// tSampler_setRate(&samp, 1.763f); // Rate of 1.0
}
+inline double getSawFall(double angle) {
+
+ angle = fmod(angle + double_Pi, 2*double_Pi); // shift x
+ double sample = angle/double_Pi - double(1); // computer as remainder
+
+ return sample;
+
+}
+
float LEAFTest_tick (float input)
{
// float sample = tNoise_tick(&noise);
@@ -83,8 +99,17 @@
// return tAutotune_tick(&at, input)[0];
- tTriangle_setFreq(&tri, x);
- return tTriangle_tick(&tri);
+ tPhasor_setFreq(&phasor, y);
+
+
+ float sample = tPhasor_tick(&phasor) * 2.0f - 1.0f;
+
+ if (phasor->phaseDidReset)
+ {
+ tMinBLEP_addBLEP(&minblep, 0.0f, 2, 0.0f);
+ }
+
+ return tMinBLEP_tick(&minblep, sample);
}
int firstFrame = 1;
@@ -100,7 +125,7 @@
float val = getSliderValue("mod freq");
- x = val * -40000.0f + 20000;
+ x = val ;
// a = val * tBuffer_getBufferLength(&buff);
@@ -108,8 +133,7 @@
val = getSliderValue("mod depth");
- y = val * 49.0f + 1.0f;
- b = val * 20.0f - 5.0f;
+ y = val * 20000.0f + 20.0f;
// b = val * tBuffer_getBufferLength(&buff);
DBG("rate: " + String(b));