# MinBLEPs

A MinBLEP is a minimum-phase bandlimited step. What is this and why do I want it? Well when you make an oscillator for your synthesizer and you use a mathematical waveform such as square and saw waves you will get aliasing because these waveforms contain harmonics above the nyquist frequency. There are various ways to solve this problem, the main two being oversampling and bandlimited impulse trains (BLITs).

The BLIT method is quite good, but I found MinBLEPs which are related to be better for my synth, as it is a faster method and supports hard sync. To understand the MinBLEP theory you will need an understanding of BLITs. Basically, you use the sinc function to generate a bandlimited impulse train, you integrate that and voila - you have a bandlimited triangle, saw or square wave.

MinBLEPs take the idea further and take a windowed sinc, perform a minimum phase reconstruction and then integrate the result and store it in a table. Now to make an oscillator you just insert a MinBLEP at each discontinuity in the waveform. So for a square wave you insert a MinBLEP where the waveform inverts, for saw wave you insert a MinBLEP where the value inverts, but you generate the ramp as normal. Note that for a discontinuity that is positive you insert a positive MinBLEP, for a negative discontinuity you insert a negative (inverted) MinBLEP.

To implement a MinBLEP oscillator, a simple method is to have an accumulator and a circular buffer. The accumulator stores a single floating point value which starts at 0.0f. The circular buffer is all zeros at the start, and you write the MinBLEPs that you add at each discontinuity into this buffer. Output from the oscillator is generated by adding the circular buffer to the accumulator value. When you insert a MinBLEP, add or subtract 1.0f from the accumulator value. Insert the MinBLEP offset by -1.0f for a positive, 1.0f for negative, so at it's first value (which is near zero) and offset it will cancel out the effect of the accumulator, but when it reaches the end of the MinBLEP, the accumulator value will no longer be cancelled as the sum of the offset and the final MinBLEP value of 1.0f will equal zero and the signal will be at the new inverted value after the step.

I hope everyone can understand this description and code, it's a confusing topic :).

# The C++ Code

 // MinBLEP Generation Code // By Daniel Werner // This Code Is Public Domain #include #define PI 3.14159265358979f // SINC Function inline float SINC(float x) {   float pix;   if(x == 0.0f)     return 1.0f;   else   {     pix = PI * x;   return sinf(pix) / pix;   } } // Generate Blackman Window inline 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 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 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 inline float cabs(float x, float y) {   return sqrtf((x * x) + (y * y)); } // Complex Exponential inline void cexp(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 void RealCepstrum(int n, float *signal, float *realCepstrum) {   float *realTime, *imagTime, *realFreq, *imagFreq;   int i;   realTime = new float[n];   imagTime = new float[n];   realFreq = new float[n];   imagFreq = new float[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(cabs(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];   delete realTime;   delete imagTime;   delete realFreq;   delete imagFreq; } // Compute Minimum Phase Reconstruction Of Signal void MinimumPhase(int n, float *realCepstrum, float *minimumPhase) {   int i, nd2;   float *realTime, *imagTime, *realFreq, *imagFreq;   nd2 = n / 2;   realTime = new float[n];   imagTime = new float[n];   realFreq = new float[n];   imagFreq = new float[n];   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++)     cexp(realFreq[i], imagFreq[i], &realFreq[i], &imagFreq[i]);   InverseDFT(n, realTime, imagTime, realFreq, imagFreq);   for(i = 0; i < n; i++)     minimumPhase[i] = realTime[i];   delete realTime;   delete imagTime;   delete realFreq;   delete imagFreq; } // Generate MinBLEP And Return It In An Array Of Floating Point Values float *GenerateMinBLEP(int zeroCrossings, int overSampling) {   int i, n;   float r, a, b;   float *buffer1, *buffer2, *minBLEP;   n = (zeroCrossings * 2 * overSampling) + 1;   buffer1 = new float[n];   buffer2 = new float[n];   // Generate Sinc   a = (float)-zeroCrossings;   b = (float)zeroCrossings;   for(i = 0; i < n; i++)   {     r = ((float)i) / ((float)(n - 1));     buffer1[i] = SINC(a + (r * (b - a)));   }   // Window Sinc   BlackmanWindow(n, buffer2);   for(i = 0; i < n; i++)     buffer1[i] *= buffer2[i];   // Minimum Phase Reconstruction   RealCepstrum(n, buffer1, buffer2);   MinimumPhase(n, buffer2, buffer1);   // Integrate Into MinBLEP   minBLEP = new float[n];   a = 0.0f;   for(i = 0; i < n; i++)   {     a += buffer1[i];     minBLEP[i] = a;   }   // Normalize   a = minBLEP[n - 1];   a = 1.0f / a;   for(i = 0; i < n; i++)     minBLEP[i] *= a;   delete buffer1;   delete buffer2;   return minBLEP; }

# References

 Author Title Eli Brandt Hard Sync Without Aliasing Tim Stilson, Julius Smith Alias-Free Digital Synthesis of Classic Analog Waveforms