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 <math.h>
#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
Links
|