|
|
|
|
 |
Entries in this category
VST SDK GUI Switch without
(Allmost) Ready-to-use oscillators
16-to-8-bit first-order dither
18dB/oct resonant 3 pole LPF with tanh() dist
1st and 2nd order pink noise filters
2 Wave shaping things
303 type filter with saturation
3rd order Spline interpollation
5-point spline interpollation
Alias-free waveform generation with analog filtering
Alien Wah
All-Pass Filters, a good explanation
Allocating aligned memory
Another 4-pole lowpass...
Another cheap sinusoidal LFO
another LFO class
Antialiased Lines
Arbitary shaped band-limited waveform generation (using oversampling and low-pass filtering)
Audiable alias free waveform gen using width sine
Bandlimited sawtooth synthesis
Bandlimited waveform generation
Bandlimited waveform generation with hard sync
Bandlimited waveforms synopsis.
Bandlimited waveforms...
Base-2 exp
Biquad C code
Bit quantization/reduction effect
Bit-Reversed Counting
Block/Loop Benchmarking
C++ class implementation of RBJ Filters
C++ gaussian noise generation
Calculate notes (java)
Cascaded resonant lp/hp filter
Center separation in a stereo mixdown
Center separation in a stereo mixdown
Cheap pseudo-sinusoidal lfo
Class for waveguide/delay effects
Clipping without branching
Coefficients for Daubechies wavelets 1-38
Compressor
Constant-time exponent of 2 detector
Conversions on a PowerPC
Cool Sounding Lowpass With Decibel Measured Resonance
Copy-protection schemes
Cubic interpollation
Cubic polynomial envelopes
DC filter
Decimator
Delay time calculation for reverberation
Denormal DOUBLE variables, macro
Denormal numbers
Denormal numbers, the meta-text
DFT
Digital RIAA equalization filter coefficients
Direct form II
Discrete Summation Formula (DSF)
Dither code
Dithering
Double to Int
Drift generator
DSF (super-set of BLIT)
Early echo's with image-mirror technique
ECE320 project: Reverberation w/ parameter control from PC
Envelope detector
Envelope Follower
Envelope follower with different attack and release
Exponential parameter mapping
fast abs/neg/sign for 32bit floats
Fast binary log approximations
Fast exp2 approximation
Fast in-place Walsh-Hadamard Transform
Fast LFO in Delphi...
Fast log2
fast power and root estimates for 32bit floats
Fast sine and cosine calculation
Fast sine wave calculation
Fast square wave generator
FFT
FFT classes in C++ and Object Pascal
Float to int
Float to int (more intel asm)
Float-to-int, coverting an array of floats
Formant filter
frequency warped FIR lattice
Gaussian dithering
Gaussian random numbers
Gaussian White noise
Gaussian White Noise
Guitar feedback
Hermite Interpolator (x86 ASM)
Hermite interpollation
Inverted parabolic envelope
Java FFT
Karlsen
Lo-Fi Crusher
Lock free fifo
Look ahead limiting
Lowpass filter for parameter edge filtering
LP and HP filter
LPC analysis (autocorrelation + Levinson-Durbin recursion)
Magnitude and phase plot of arbitrary IIR function, up to 5th order
MATLAB-Tools for SNDAN
Measuring interpollation noise
MIDI note/frequency conversion
Millimeter to DB (faders...)
Moog VCF
Moog VCF, variation 1
Moog VCF, variation 2
Most simple and smooth feedback delay
Most simple static delay
Noise Shaping Class
Nonblocking multiprocessor/multithread algorithms in C++
Notch filter
One pole LP and HP
One pole, one zero LP/HP
One zero, LP/HP
Parallel combs delay calculation
Peak/Notch filter
Phase equalization
Phase modulation Vs. Frequency modulation
Phase modulation Vs. Frequency modulation II
Phaser code
Pink noise filter
Polyphase Filters
pow(x,4) approximation
Prewarping
Pseudo-Random generator
Pulsewidth modulation
QFT and DQFT (double precision) classes
quick and dirty sine generator
RBJ-Audio-EQ-Cookbook
Reading the compressed WA! parts in gigasampler files
Real basic DSP with Matlab (+ GUI) ...
real value vs display value
Really fast x86 floating point sin/cos
Reasonably accurate/fastish tanh approximation
resampling
Resonant filter
Resonant IIR lowpass (12dB/oct)
Resonant low pass filter
Reverb Filter Generator
Reverberation Algorithms in Matlab
Reverberation techniques
Saturation
SawSin
Simple peak follower
Sin(x) Aproximation (with SSE code)
Sin, Cos, Tan approximation
Sine calculation
smsPitchScale Source Code
Soft saturation
Square Waves
State variable
State Variable Filter (Chamberlin version)
State Variable Filter (Double Sampled, Stable)
Stereo Enhancer
Stilson's Moog filter code
Time compression-expansion using standard phase vocoder
Time domain convolution with O(n^log2(3))
Time domain convolution with O(n^log2(3))
tone detection with Goertzel
Tone detection with Goertzel (x86 ASM)
transistor differential amplifier simulation
Variable-hardness clipping function
Various Biquad filters
Waveform generator using MinBLEPS
WaveShaper
Waveshaper
Waveshaper
Waveshaper (simple description)
Waveshaper :: Gloubi-boulga
Wavetable Synthesis
Weird synthesis
Zoelzer biquad filters
|
|
|
|
|
|
|
|
|
VST SDK GUI Switch without
References : Posted by quintosardo[AT]yahoo[DOT]it
Notes :
In VST GUI an on-vaue is represented by 1.0 and off by 0.0.
Code :
Say you have two signals you want to switch between when the user changes a switch.
You could do:
if(fSwitch == 0.f) //fSwitch is either 0.0 or 1.0
output = input1
else
output = input2
However, you can avoid the branch by doing:
output = input1 * (1.f - fSwitch) + input2 * fSwitch
Which would be like a quick mix. You could make the change clickless by adding a simple one-pole filter:
smooth = filter(fSwitch)
output = input1 * (1.f - smooth) + input2 * smooth
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
(Allmost) Ready-to-use oscillators
Type : waveform generation
References : Ross Bencina, Olli Niemitalo, ...
Notes :
Ross Bencina: original source code poster
Olli Niemitalo: UpdateWithCubicInterpolation
Code :
//this code is meant as an EXAMPLE
//uncomment if you need an FM oscillator
//define FM_OSCILLATOR
/*
members are:
float phase;
int TableSize;
float sampleRate;
float *table, dtable0, dtable1, dtable2, dtable3;
->these should be filled as folows... (remember to wrap around!!!)
table[i] = the wave-shape
dtable0[i] = table[i+1] - table[i];
dtable1[i] = (3.f*(table[i]-table[i+1])-table[i-1]+table[i+2])/2.f
dtable2[i] = 2.f*table[i+1]+table[i-1]-(5.f*table[i]+table[i+2])/2.f
dtable3[i] = (table[i+1]-table[i-1])/2.f
*/
float Oscillator::UpdateWithoutInterpolation(float frequency)
{
int i = (int) phase;
phase += (sampleRate/(float TableSize)/frequency;
if(phase >= (float)TableSize)
phase -= (float)TableSize;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
return table[i] ;
}
float Oscillator::UpdateWithLinearInterpolation(float frequency)
{
int i = (int) phase;
float alpha = phase - (float) i;
phase += (sampleRate/(float)TableSize)/frequency;
if(phase >= (float)TableSize)
phase -= (float)TableSize;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
/*
dtable0[i] = table[i+1] - table[i]; //remember to wrap around!!!
*/
return table[i] + dtable0[i]*alpha;
}
float Oscillator::UpdateWithCubicInterpolation( float frequency )
{
int i = (int) phase;
float alpha = phase - (float) i;
phase += (sampleRate/(float)TableSize)/frequency;
if(phase >= (float)TableSize)
phase -= (float)TableSize;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
/* //remember to wrap around!!!
dtable1[i] = (3.f*(table[i]-table[i+1])-table[i-1]+table[i+2])/2.f
dtable2[i] = 2.f*table[i+1]+table[i-1]-(5.f*table[i]+table[i+2])/2.f
dtable3[i] = (table[i+1]-table[i-1])/2.f
*/
return ((dtable1[i]*alpha + dtable2[i])*alpha + dtable3[i])*alpha+table[i];
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
16-to-8-bit first-order dither
Type : First order error feedforward dithering code
References : Posted by Jon Watte
Notes :
This
is about as simple a dithering algorithm as you can implement, but it's
likely to sound better than just truncating to N bits.
Note that you might not want to carry forward the full difference for
infinity. It's probably likely that the worst performance hit comes
from the saturation conditionals, which can be avoided with appropriate
instructions on many DSPs and integer SIMD type instructions, or CMOV.
Last, if sound quality is paramount (such as when going from > 16
bits to 16 bits) you probably want to use a higher-order dither
function found elsewhere on this site.
Code :
// This code will down-convert and dither a 16-bit signed short
// mono signal into an 8-bit unsigned char signal, using a first
// order forward-feeding error term dither.
#define uchar unsigned char
void dither_one_channel_16_to_8( short * input, uchar * output, int count, int * memory )
{
int m = *memory;
while( count-- > 0 ) {
int i = *input++;
i += m;
int j = i + 32768 - 128;
uchar o;
if( j < 0 ) {
o = 0;
}
else if( j > 65535 ) {
o = 255;
}
else {
o = (uchar)((j>>8)&0xff);
}
m = ((j-32768+128)-i);
*output++ = o;
}
*memory = m;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
1st and 2nd order pink noise filters
Type : Pink noise
References : Posted by umminger[AT]umminger[DOT]com
Notes :
Here are some new lower-order pink noise filter coefficients.
These have approximately equiripple error in decibels from 20hz to 20khz at a 44.1khz sampling rate.
1st order, ~ +/- 3 dB error (not recommended!)
num = [0.05338071119116 -0.03752455712906]
den = [1.00000000000000 -0.97712493947102]
2nd order, ~ +/- 0.9 dB error
num = [ 0.04957526213389 -0.06305581334498 0.01483220320740 ]
den = [ 1.00000000000000 -1.80116083982126 0.80257737639225 ]
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
2 Wave shaping things
References : Posted by Frederic Petrot
Notes :
Makes nice saturations effects that can be easilly computed using cordic
First using a atan function:
y1 using k=16
max is the max value you can reach (32767 would be a good guess)
Harmonics scale down linealy and not that fast
Second using the hyperbolic tangent function:
y2 using k=2
Harmonics scale down linealy very fast
Code :
y1 = (max>>1) * atan(k * x/max)
y2 = max * th(x/max)
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
303 type filter with saturation
Type : Runge-Kutta Filters
References : Posted by Hans Mikelson
Linked file : filters001.txt
Notes :
I
posted a filter to the Csound mailing list a couple of weeks ago that
has a 303 flavor to it. It basically does wacky distortions to the
sound. I used Runge-Kutta for the diff eq. simulation though which
makes it somewhat sluggish.
This is a CSound score!!
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
3rd order Spline interpollation
References : Posted by Dave from Muon Software, originally from Josh Scholar
Notes :
(from Joshua Scholar about Spline interpollation in general...)
According to sampling theory, a perfect interpolation could be found by
replacing each sample with a sinc function centered on that sample,
ringing at your target nyquest frequency, and at each target point you
just sum all of contributions from the sinc functions of every single
point in source.
The sinc function has ringing that dies away very slowly, so each
target sample will have to have contributions from a large neighborhood
of source samples. Luckily, by definition the sinc function is
bandwidth limited, so once we have a source that is prefilitered for
our target nyquest frequency and reasonably oversampled relative to our
nyquest frequency, ordinary interpolation techniques are quite fruitful
even though they would be pretty useless if we hadn't oversampled.
We want an interpolation routine that at very least has the following characteristics:
1. Obviously it's continuous. But since finite differencing a signal (I
don't really know about true differentiation) is equivalent to a low
frequency attenuator that drops only about 6 dB per octave, continuity
at the higher derivatives is important too.
2. It has to be stiff enough to find peaks when our oversampling missed
them. This is where what I said about the combination the sinc
function's limited bandwidth and oversampling making interpolation
possible comes into play.
I've read some papers on splines, but most stuff on splines relates to
graphics and uses a control point descriptions that is completely
irrelevant to our sort of interpolation. In reading this stuff I
quickly came to the conclusion that splines:
1. Are just piecewise functions made of polynomials designed to have some higher order continuity at the transition points.
2. Splines are highly arbitrary, because you can choose arbitrary
derivatives (to any order) at each transition. Of course the more you
specify the higher order the polynomials will be.
3. I already know enough about polynomials to construct any sort of
spline. A polynomial through 'n' points with a derivative specified at
'm[1]' points and second derivatives specified at 'm[2]' points etc.
will be a polynomial of the order n-1+m[1]+m[2]...
A way to construct third order splines (that admittedly doesn't help
you construct higher order splines), is to linear interpolate between
two parabolas. At each point (they are called knots) you have a
parabola going through that point, the previous and the next point.
Between each point you linearly interpolate between the polynomials for
each point. This may help you imagine splines.
As a starting point I used a polynomial through 5 points for each knot
and used MuPad (a free Mathematica like program) to derive a polynomial
going through two points (knots) where at each point it has the same
first two derivatives as a 4th order polynomial through the surrounding
5 points. My intuition was that basing it on polynomials through 3
points wouldn't be enough of a neighborhood to get good continuity.
When I tested it, I found that not only did basing it on 5 point
polynomials do much better than basing it on 3 point ones, but that 7
point ones did nearly as badly as 3 point ones. 5 points seems to be a
sweet spot.
However, I could have set the derivatives to a nearly arbitrary values
- basing the values on those of polynomials through the surrounding
points was just a guess.
I've read that the math of sampling theory has different interpretation
to the sinc function one where you could upsample by making a
polynomial through every point at the same order as the number of
points and this would give you the same answer as sinc function
interpolation (but this only converges perfectly when there are an
infinite number of points). Your head is probably spinning right now -
the only point of mentioning that is to point out that perfect
interpolation is exactly as stiff as a polynomial through the target
points of the same order as the number of target points.
Code :
//interpolates between L0 and H0 taking the previous (L1) and next (H1)
points into account
inline float ThirdInterp(const float x,const float L1,const float L0,const
float H0,const float H1)
{
return
L0 +
.5f*
x*(H0-L1 +
x*(H0 + L0*(-2) + L1 +
x*( (H0 - L0)*9 + (L1 - H1)*3 +
x*((L0 - H0)*15 + (H1 - L1)*5 +
x*((H0 - L0)*6 + (L1 - H1)*2 )))));
}
5 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
5-point spline interpollation
Type : interpollation
References : Joshua Scholar, posted by David Waugh
Code :
//nMask = sizeofwavetable-1 where sizeofwavetable is a power of two.
double interpolate(double* wavetable, int nMask, double location)
{
/* 5-point spline*/
int nearest_sample = (int) location;
double x = location - (double) nearest_sample;
double p0=wavetable[(nearest_sample-2)&nMask];
double p1=wavetable[(nearest_sample-1)&nMask];
double p2=wavetable[nearest_sample];
double p3=wavetable[(nearest_sample+1)&nMask];
double p4=wavetable[(nearest_sample+2)&nMask];
double p5=wavetable[(nearest_sample+3)&nMask];
return p2 + 0.04166666666*x*((p3-p1)*16.0+(p0-p4)*2.0
+ x *((p3+p1)*16.0-p0-p2*30.0- p4
+ x *(p3*66.0-p2*70.0-p4*33.0+p1*39.0+ p5*7.0- p0*9.0
+ x *( p2*126.0-p3*124.0+p4*61.0-p1*64.0- p5*12.0+p0*13.0
+ x *((p3-p2)*50.0+(p1-p4)*25.0+(p5-p0)*5.0)))));
};
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Allocating aligned memory
Type : memory allocation
References : Posted by Benno Senoner
Notes :
we waste up to align_size + sizeof(int) bytes when we alloc a memory area.
We store the aligned_ptr - unaligned_ptr delta in an int located before the aligned area.
This is needed for the free() routine since we need to free all the memory not only the aligned area.
You have to use aligned_free() to free the memory allocated with aligned_malloc() !
Code :
/* align_size has to be a power of two !! */
void *aligned_malloc(size_t size, size_t align_size) {
char *ptr,*ptr2,*aligned_ptr;
int align_mask = align_size - 1;
ptr=(char *)malloc(size + align_size + sizeof(int));
if(ptr==NULL) return(NULL);
ptr2 = ptr + sizeof(int);
aligned_ptr = ptr2 + (align_size - ((size_t)ptr2 & align_mask));
ptr2 = aligned_ptr - sizeof(int);
*((int *)ptr2)=(int)(aligned_ptr - ptr);
return(aligned_ptr);
}
void aligned_free(void *ptr) {
int *ptr2=(int *)ptr - 1;
ptr -= *ptr2;
free(ptr);
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Another 4-pole lowpass...
Type : 4-pole LP/HP
References : Posted by fuzzpilz [AT] gmx [DOT] net
Notes :
Vaguely based on the Stilson/Smith Moog paper, but going in a rather different direction from others I've seen here.
The parameters are peak frequency and peak magnitude (g below); both
are reasonably accurate for magnitudes above 1. DC gain is 1.
The filter has some undesirable properties - e.g. it's unstable for low
peak freqs if implemented in single precision (haven't been able to
cleanly separate it into biquads or onepoles to see if that helps), and
it responds so strongly to parameter changes that it's not advisable to
update the coefficients much more rarely than, say, every eight samples
during sweeps, which makes it somewhat expensive.
I like the sound, however, and the accuracy is nice to have, since many filters are not very strong in that respect.
I haven't looked at the HP again for a while, but IIRC it had approximately the same good and bad sides.
Code :
double coef[9];
double d[4];
double omega; //peak freq
double g; //peak mag
// calculating coefficients:
double k,p,q,a;
double a0,a1,a2,a3,a4;
k=(4.0*g-3.0)/(g+1.0);
p=1.0-0.25*k;p*=p;
// LP:
a=1.0/(tan(0.5*omega)*(1.0+p));
p=1.0+a;
q=1.0-a;
a0=1.0/(k+p*p*p*p);
a1=4.0*(k+p*p*p*q);
a2=6.0*(k+p*p*q*q);
a3=4.0*(k+p*q*q*q);
a4= (k+q*q*q*q);
p=a0*(k+1.0);
coef[0]=p;
coef[1]=4.0*p;
coef[2]=6.0*p;
coef[3]=4.0*p;
coef[4]=p;
coef[5]=-a1*a0;
coef[6]=-a2*a0;
coef[7]=-a3*a0;
coef[8]=-a4*a0;
// or HP:
a=tan(0.5*omega)/(1.0+p);
p=a+1.0;
q=a-1.0;
a0=1.0/(p*p*p*p+k);
a1=4.0*(p*p*p*q-k);
a2=6.0*(p*p*q*q+k);
a3=4.0*(p*q*q*q-k);
a4= (q*q*q*q+k);
p=a0*(k+1.0);
coef[0]=p;
coef[1]=-4.0*p;
coef[2]=6.0*p;
coef[3]=-4.0*p;
coef[4]=p;
coef[5]=-a1*a0;
coef[6]=-a2*a0;
coef[7]=-a3*a0;
coef[8]=-a4*a0;
// per sample:
out=coef[0]*in+d[0];
d[0]=coef[1]*in+coef[5]*out+d[1];
d[1]=coef[2]*in+coef[6]*out+d[2];
d[2]=coef[3]*in+coef[7]*out+d[3];
d[3]=coef[4]*in+coef[8]*out;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Another cheap sinusoidal LFO
References : Posted by info[at]e-phonic[dot]com
Notes :
Some pseudo code for a easy to calculate LFO.
You can even make a rough triangle wave out of this by substracting the output of 2 of these with different phases.
PJ
Code :
r = the rate 0..1
--------------
p += r
if(p > 1) p -= 2;
out = p*(1-abs(p));
--------------
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
another LFO class
References : Posted by mdsp
Linked file : LFO.zip
Notes :
This
LFO uses an unsigned 32-bit phase and increment whose 8 Most
Significant Bits adress a Look-up table while the 24 Least Significant
Bits are used as the fractionnal part.
Note: As the phase overflow automatically, the index is always in the range 0-255.
It performs linear interpolation, but it is easy to add other types of interpolation.
Don't know how good it could be as an oscillator, but I found it good enough for a LFO.
BTW there is also different kind of waveforms.
Modifications:
We could use phase on 64-bit or change the proportion of bits used by the index and the fractionnal part.
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Antialiased Lines
Type : A slow, ugly, and unoptimized but short method to perform antialiased lines in a framebuffer
References : Posted by arguru[AT]smartelectronix[DOT]com
Notes :
Simple code to perform antialiased lines in a 32-bit RGBA (1 byte/component) framebuffer.
pframebuffer <- unsigned char* to framebuffer bytes (important: Y
flipped line order! [like in the way Win32 CreateDIBSection works...])
client_height=framebuffer height in lines
client_width=framebuffer width in pixels (not in bytes)
This doesnt perform any clip checl so it fails if coordinates are set out of bounds.
sorry for the engrish
Code :
//
// By Arguru
//
void PutTransPixel(int const x,int const y,UCHAR const r,UCHAR const g,UCHAR const b,UCHAR const a)
{
unsigned char* ppix=pframebuffer+(x+(client_height-(y+1))*client_width)*4;
ppix[0]=((a*b)+(255-a)*ppix[0])/256;
ppix[1]=((a*g)+(255-a)*ppix[1])/256;
ppix[2]=((a*r)+(255-a)*ppix[2])/256;
}
void LineAntialiased(int const x1,int const y1,int const x2,int const y2,UCHAR const r,UCHAR const g,UCHAR const b)
{
// some useful constants first
double const dw=x2-x1;
double const dh=y2-y1;
double const slx=dh/dw;
double const sly=dw/dh;
// determine wichever raster scanning behaviour to use
if(fabs(slx)<1.0)
{
// x scan
int tx1=x1;
int tx2=x2;
double raster=y1;
if(x1>x2)
{
tx1=x2;
tx2=x1;
raster=y2;
}
for(int x=tx1;x<=tx2;x++)
{
int const ri=int(raster);
double const in_y0=1.0-(raster-ri);
double const in_y1=1.0-(ri+1-raster);
PutTransPixel(x,ri+0,r,g,b,in_y0*255.0);
PutTransPixel(x,ri+1,r,g,b,in_y1*255.0);
raster+=slx;
}
}
else
{
// y scan
int ty1=y1;
int ty2=y2;
double raster=x1;
if(y1>y2)
{
ty1=y2;
ty2=y1;
raster=x2;
}
for(int y=ty1;y<=ty2;y++)
{
int const ri=int(raster);
double const in_x0=1.0-(raster-ri);
double const in_x1=1.0-(ri+1-raster);
PutTransPixel(ri+0,y,r,g,b,in_x0*255.0);
PutTransPixel(ri+1,y,r,g,b,in_x1*255.0);
raster+=sly;
}
}
}
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Arbitary shaped band-limited waveform generation (using oversampling and low-pass filtering)
References : Posted by remage[AT]kac[DOT]poliod[DOT]hu
Code :
Arbitary shaped band-limited waveform generation
(using oversampling and low-pass filtering)
There are many articles about band-limited waveform synthesis
techniques, that provide correct and fast methods for generating
classic analogue waveforms, such as saw, pulse, and triangle
wave. However, generating arbitary shaped band-limited waveforms, such
as the "sawsin" shape (found in this source-code archive), seems to be
quite hard using these techniques.
My analogue waveforms are generated in a _very_ high sampling rate
(actually it's 1.4112 GHz for 44.1 kHz waveforms, using 32x
oversampling). Using this sample-rate, the amplitude of the aliasing
harmonics are negligible (the base analogue waveforms has exponentially
decreasing harmonics amplitudes).
Using a 511-tap windowed sync FIR filter (with Blackman-Harris window,
and 12 kHz cutoff frequency) the harmonics above 20 kHz are killed, the
higher harmonics (that cause the sharp overshoot at step response) are
dampened.
The filtered signal downsampled to 44.1 kHz contains the audible (non-aliased) harmonics only.
This waveform synthesis is performed for wavetables of 4096, 2048,
1024, ... 8, 4, 2 samples. The real-time signal is interpolated from
these waveform-tables, using Hermite-(cubic-)interpolation for the
waveforms, and linear interpolation between the two wavetables near the
required note.
This procedure is quite time-consuming, but the whole waveform (or, in
my implementation, the whole waveform-set) can be precalculated (or
saved at first launch of the synth) and reloaded at synth
initialization.
I don't know if this is a theoretically correct solution, but the
waveforms sound good (no audible aliasing). Please let me know if I'm
wrong...
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Audiable alias free waveform gen using width sine
Type : Very simple
References : Posted by joakim[DOT]dahlstrom[AT]ongame[DOT]com
Notes :
Warning, my english abilities is terribly limited.
How ever, the other day when finally understanding what bandlimited
wave creation is (i am a noobie, been doing DSP stuf on and off for a
half/year) it hit me i can implement one little part in my synths. It's
all about the freq (that i knew), very simple you can reduce alias (the
alias that you can hear that is) extremely by keeping track of your
frequence, the way i solved it is using a factor, afact = 1 -
sin(f*2PI). This means you can do audiable alias free synthesis without
very complex algorithms or very huge tables, even though the sound
becomes kind of low-filtered.
Propably something like this is mentioned b4, but incase it hasn't this is worth looking up
The psuedo code describes it more.
// Druttis
Code :
f := freq factor, 0 - 0.5 (0 to half samplingrate)
afact(f) = 1 - sin(f*2PI)
t := time (0 to ...)
ph := phase shift (0 to 1)
fm := freq mod (0 to 1)
sine(t,f,ph,fm) = sin((t*f+ph)*2PI + 0.5PI*fm*afact(f))
fb := feedback (0 to 1) (1 max saw)
saw(t,f,ph,fm,fb) = sine(t,f,ph,fb*sine(t-1,f,ph,fm))
pm := pulse mod (0 to 1) (1 max pulse)
pw := pulse width (0 to 1) (1 square)
pulse(t,f,ph,fm,fb,pm,pw) = saw(t,f,ph,fm,fb) - (t,f,ph+0.5*pw,fm,fb) * pm
I am not completely sure about fm for saw & pulse since i cant test
that atm. but it should work :) otherwise just make sure fm are 0 for
saw & pulse.
As you can see the saw & pulse wave are very variable.
// Druttis
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Bandlimited sawtooth synthesis
Type : DSF BLIT
References : Posted by emanuel.landeholm [AT] telia.com
Linked file : synthesis002.txt
Notes :
This
is working code for synthesizing a bandlimited sawtooth waveform. The
algorithm is DSF BLIT + leaky integrator. Includes driver code.
There are two parameters you may tweak:
1) Desired attenuation at nyquist. A low value yields a duller sawtooth
but gets rid of those annoying CLICKS when sweeping the frequency up
real high. Must be strictly less than 1.0!
2) Integrator leakiness/cut off. Affects the shape of the waveform to
some extent, esp. at the low end. Ideally you would want to set this
low, but too low a setting will give you problems with DC.
Have fun!
/Emanuel Landeholm
(see linked file)
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Bandlimited waveforms...
References : Posted by Paul Kellet
Notes :
(Quoted from Paul's mail)
Below is another waveform generation method based on a train of sinc
functions (actually an alternating loop along a sinc between t=0 and
t=period/2).
The code integrates the pulse train with a dc offset to get a sawtooth,
but other shapes can be made in the usual ways... Note that 'dc' and
'leak' may need to be adjusted for very high or low frequencies.
I don't know how original it is (I ought to read more) but it is of
usable quality, particularly at low frequencies. There's some scope for
optimisation by using a table for sinc, or maybe a a truncated/windowed
sinc?
I think it should be possible to minimise the aliasing by fine tuning
'dp' to slightly less than 1 so the sincs join together neatly, but I
haven't found the best way to do it. Any comments gratefully received.
Code :
float p=0.0f; //current position
float dp=1.0f; //change in postion per sample
float pmax; //maximum position
float x; //position in sinc function
float leak=0.995f; //leaky integrator
float dc; //dc offset
float saw; //output
//set frequency...
pmax = 0.5f * getSampleRate() / freqHz;
dc = -0.498f/pmax;
//for each sample...
p += dp;
if(p < 0.0f)
{
p = -p;
dp = -dp;
}
else if(p > pmax)
{
p = pmax + pmax - p;
dp = -dp;
}
x= pi * p;
if(x < 0.00001f)
x=0.00001f; //don't divide by 0
saw = leak*saw + dc + (float)sin(x)/(x);
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Base-2 exp
References : Posted by Laurent de Soras
Notes :
Linear approx. between 2 integer values of val. Uses 32-bit integers. Not very efficient but fastest than exp()
This code was designed for x86 (little endian), but could be adapted for big endian processors.
Laurent thinks you just have to change the (*(1 + (int *) &ret))
expressions and replace it by (*(int *) &ret). However, He didn't
test it.
Code :
inline double fast_exp2 (const double val)
{
int e;
double ret;
if (val >= 0)
{
e = int (val);
ret = val - (e - 1);
((*(1 + (int *) &ret)) &= ~(2047 << 20)) += (e + 1023) << 20;
}
else
{
e = int (val + 1023);
ret = val - (e - 1024);
((*(1 + (int *) &ret)) &= ~(2047 << 20)) += e << 20;
}
return (ret);
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Bit quantization/reduction effect
Type : Bit-level noise-generating effect
References : Posted by Jon Watte
Notes :
This
function, run on each sample, will emulate half the effect of running
your signal through a Speak-N-Spell or similar low-bit-depth circuitry.
The other half would come from downsampling with no aliasing control, i
e replicating every N-th sample N times in the output signal.
Code :
short keep_bits_from_16( short input, int keepBits ) {
return (input & (-1 << (16-keepBits)));
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Bit-Reversed Counting
References : Posted by mailbjl[AT]yahoo[DOT]com
Notes :
Bit-reversed
ordering comes up frequently in FFT implementations. Here is a
non-branching algorithm (given in C) that increments the variable "s"
bit-reversedly from 0 to N-1, where N is a power of 2.
Code :
int r = 0; // counter
int s = 0; // bit-reversal of r/2
int N = 256; // N can be any power of 2
int N2 = N << 1; // N<<1 == N*2
do {
printf("%u ", s);
r += 2;
s ^= N - (N / (r&-r));
}
while (r < N2);
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Block/Loop Benchmarking
Type : Benchmarking Tool
References : Posted by arguru[AT]smartelectronix[DOT]com
Notes :
Requires CPU with RDTSC support
Code :
// Block-Process Benchmarking Code using rdtsc
// useful for measure DSP block stuff
// (based on Intel papers)
// 64-bit precission
// VeryUglyCode(tm) by Arguru
// globals
UINT time,time_low,time_high;
// call this just before enter your loop or whatever
void bpb_start()
{
// read time stamp to EAX
__asm rdtsc;
__asm mov time_low,eax;
__asm mov time_high,edx;
}
// call the following function just after your loop
// returns average cycles wasted per sample
UINT bpb_finish(UINT const num_samples)
{
__asm rdtsc
__asm sub eax,time_low;
__asm sub edx,time_high;
__asm div num_samples;
__asm mov time,eax;
return time;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
C++ gaussian noise generation
Type : gaussian noise generation
References : Posted by paul[at]expdigital[dot]co[dot]uk
Notes :
References :
Tobybears delphi noise generator was the basis. Simply converted it to C++.
Link for original is:
http://www.musicdsp.org/archive.php?classid=0#129
The output is in noise.
Code :
/* Include requisits */
#include <cstdlib>
#include <ctime>
/* Generate a new random seed from system time - do this once in your constructor */
srand(time(0));
/* Setup constants */
const static int q = 15;
const static float c1 = (1 << q) - 1;
const static float c2 = ((int)(c1 / 3)) + 1;
const static float c3 = 1.f / c1;
/* random number in range 0 - 1 not including 1 */
float random = 0.f;
/* the white noise */
float noise = 0.f;
for (int i = 0; i < numSamples; i++)
{
random = ((float)rand() / (float)(RAND_MAX + 1));
noise = (2.f * ((random * c2) + (random * c2) + (random * c2)) - 3.f * (c2 - 1.f)) * c3;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Calculate notes (java)
Type : Java class for calculating notes with different in params
References : Posted by larsby[AT]elak[DOT]org
Linked file : Frequency.java
Notes :
Converts
between string notes and frequencies and back. I vaguely remember
writing bits of it, and I got it off the net somwhere so dont ask me
- Larsby
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cascaded resonant lp/hp filter
Type : lp+hp
References : Posted by tobybear[AT]web[DOT]de
Notes :
// Cascaded resonant lowpass/hipass combi-filter
// The original source for this filter is from Paul Kellet from
// the archive. This is a cascaded version in Delphi where the
// output of the lowpass is fed into the highpass filter.
// Cutoff frequencies are in the range of 0<=x<1 which maps to
// 0..nyquist frequency
// input variables are:
// cut_lp: cutoff frequency of the lowpass (0..1)
// cut_hp: cutoff frequency of the hipass (0..1)
// res_lp: resonance of the lowpass (0..1)
// res_hp: resonance of the hipass (0..1)
Code :
var n1,n2,n3,n4:single; // filter delay, init these with 0!
fb_lp,fb_hp:single; // storage for calculated feedback
const p4=1.0e-24; // Pentium 4 denormal problem elimination
function dofilter(inp,cut_lp,res_lp,cut_hp,res_hp:single):single;
begin
fb_lp:=res_lp+res_lp/(1-cut_lp);
fb_hp:=res_hp+res_hp/(1-cut_lp);
n1:=n1+cut_lp*(inp-n1+fb_lp*(n1-n2))+p4;
n2:=n2+cut_lp*(n1-n2);
n3:=n3+cut_hp*(n2-n3+fb_hp*(n3-n4))+p4;
n4:=n4+cut_hp*(n3-n4);
result:=i-n4;
end;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Center separation in a stereo mixdown
References : Posted by Thiburce BELAVENTURE
Notes :
One year ago, i found a little trick to isolate or remove the center in a stereo mixdown.
My method use the time-frequency representation (FFT). I use a min
fuction between left and right channels (for each bin) to create the
pseudo center. I apply a phase correction, and i substract this signal
to the left and right signals.
Then, we can remix them after treatments (or without) to produce a stereo signal in output.
This algorithm (I called it "TBIsolator") is not perfect, but the
result is very nice, better than the phase technic (L substract R...).
I know that it is not mathematically correct, but as an estimation of
the center, the exact match is very hard to obtain. So, it is not so
bad (just listen the result and see).
My implementation use a 4096 FFT size, with overlap-add method (factor
2). With a lower FFT size, the sound will be more dirty, and with a
16384 FFT size, the center will have too much high frequency (I don't
explore why this thing appears).
I just post the TBIsolator code (see FFTReal in this site for implement the FFT engine).
pIns and pOuts buffers use the representation of the FFTReal class (0 to N/2-1: real parts, N/2 to N-1: imaginary parts).
Have fun with the TBIsolator algorithm ! I hope you enjoy it and if you enhance it, contact me (it's my baby...).
P.S.: the following function is not optimized.
Code :
/* ============================================================= */
/* nFFTSize must be a power of 2 */
/* ============================================================= */
/* Usage examples: */
/* - suppress the center: fAmpL = 1.f, fAmpC = 0.f, fAmpR = 1.f */
/* - keep only the center: fAmpL = 0.f, fAmpC = 1.f, fAmpR = 0.f */
/* ============================================================= */
void processTBIsolator(float *pIns[2], float *pOuts[2], long nFFTSize, float fAmpL, float fAmpC, float fAmpR)
{
float fModL, fModR;
float fRealL, fRealC, fRealR;
float fImagL, fImagC, fImagR;
double u;
for ( long i = 0, j = nFFTSize / 2; i < nFFTSize / 2; i++ )
{
fModL = pIns[0][i] * pIns[0][i] + pIns[0][j] * pIns[0][j];
fModR = pIns[1][i] * pIns[1][i] + pIns[1][j] * pIns[1][j];
// min on complex numbers
if ( fModL > fModR )
{
fRealC = fRealR;
fImagC = fImagR;
}
else
{
fRealC = fRealL;
fImagC = fImagL;
}
// phase correction...
u = fabs(atan2(pIns[0][j], pIns[0][i]) - atan2(pIns[1][j], pIns[1][i])) / 3.141592653589;
if ( u >= 1 ) u -= 1.;
u = pow(1 - u*u*u, 24);
fRealC *= (float) u;
fImagC *= (float) u;
// center extraction...
fRealL = pIns[0][i] - fRealC;
fImagL = pIns[0][j] - fImagC;
fRealR = pIns[1][i] - fRealC;
fImagR = pIns[1][j] - fImagC;
// You can do some treatments here...
pOuts[0][i] = fRealL * fAmpL + fRealC * fAmpC;
pOuts[0][j] = fImagL * fAmpL + fImagC * fAmpC;
pOuts[1][i] = fRealR * fAmpR + fRealC * fAmpC;
pOuts[1][j] = fImagR * fAmpR + fImagC * fAmpC;
}
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Center separation in a stereo mixdown
References : Posted by Thiburce BELAVENTURE
Notes :
One year ago, i found a little trick to isolate or remove the center in a stereo mixdown.
My method use the time-frequency representation (FFT). I use a min
fuction between left and right channels (for each bin) to create the
pseudo center. I apply a phase correction, and i substract this signal
to the left and right signals.
Then, we can remix them after treatments (or without) to produce a stereo signal in output.
This algorithm (I called it "TBIsolator") is not perfect, but the
result is very nice, better than the phase technic (L substract R...).
I know that it is not mathematically correct, but as an estimation of
the center, the exact match is very hard to obtain. So, it is not so
bad (just listen the result and see).
My implementation use a 4096 FFT size, with overlap-add method (factor
2). With a lower FFT size, the sound will be more dirty, and with a
16384 FFT size, the center will have too much high frequency (I don't
explore why this thing appears).
I just post the TBIsolator code (see FFTReal in this site for implement the FFT engine).
pIns and pOuts buffers use the representation of the FFTReal class (0 to N/2-1: real parts, N/2 to N-1: imaginary parts).
Have fun with the TBIsolator algorithm ! I hope you enjoy it and if you enhance it, contact me (it's my baby...).
P.S.: the following function is not optimized.
Code :
/* ============================================================= */
/* nFFTSize must be a power of 2 */
/* ============================================================= */
/* Usage examples: */
/* - suppress the center: fAmpL = 1.f, fAmpC = 0.f, fAmpR = 1.f */
/* - keep only the center: fAmpL = 0.f, fAmpC = 1.f, fAmpR = 0.f */
/* ============================================================= */
void processTBIsolator(float *pIns[2], float *pOuts[2], long nFFTSize, float fAmpL, float fAmpC, float fAmpR)
{
float fModL, fModR;
float fRealL, fRealC, fRealR;
float fImagL, fImagC, fImagR;
double u;
for ( long i = 0, j = nFFTSize / 2; i < nFFTSize / 2; i++ )
{
fModL = pIns[0][i] * pIns[0][i] + pIns[0][j] * pIns[0][j];
fModR = pIns[1][i] * pIns[1][i] + pIns[1][j] * pIns[1][j];
// min on complex numbers
if ( fModL > fModR )
{
fRealC = fRealR;
fImagC = fImagR;
}
else
{
fRealC = fRealL;
fImagC = fImagL;
}
// phase correction...
u = fabs(atan2(pIns[0][j], pIns[0][i]) - atan2(pIns[1][j], pIns[1][i])) / 3.141592653589;
if ( u >= 1 ) u -= 1.;
u = pow(1 - u*u*u, 24);
fRealC *= (float) u;
fImagC *= (float) u;
// center extraction...
fRealL = pIns[0][i] - fRealC;
fImagL = pIns[0][j] - fImagC;
fRealR = pIns[1][i] - fRealC;
fImagR = pIns[1][j] - fImagC;
// You can do some treatments here...
pOuts[0][i] = fRealL * fAmpL + fRealC * fAmpC;
pOuts[0][j] = fImagL * fAmpL + fImagC * fAmpC;
pOuts[1][i] = fRealR * fAmpR + fRealC * fAmpC;
pOuts[1][j] = fImagR * fAmpR + fImagC * fAmpC;
}
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cheap pseudo-sinusoidal lfo
References : Posted by fumminger[AT]umminger[DOT]com
Notes :
Although
the code is written in standard C++, this algorithm is really better
suited for dsps where one can take advantage of multiply-accumulate
instructions and where the required phase accumulator can be easily
implemented by masking a counter.
It provides a pretty cheap roughly sinusoidal waveform that is good enough for an lfo.
Code :
// x should be between -1.0 and 1.0
inline
double pseudo_sine(double x)
{
// Compute 2*(x^2-1.0)^2-1.0
x *= x;
x -= 1.0;
x *= x;
// The following lines modify the range to lie between -1.0 and 1.0.
// If a range of between 0.0 and 1.0 is acceptable or preferable
// (as in a modulated delay line) then you can save some cycles.
x *= 2.0;
x -= 1.0;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Class for waveguide/delay effects
Type : IIR filter
References : Posted by arguru[AT]smartelectronix.com
Notes :
Flexible-time, non-sample quantized delay , can be used for stuff like waveguide synthesis or time-based (chorus/flanger) fx.
MAX_WG_DELAY is a constant determining MAX buffer size (in samples)
Code :
class cwaveguide
{
public:
cwaveguide(){clear();}
virtual ~cwaveguide(){};
void clear()
{
counter=0;
for(int s=0;s<MAX_WG_DELAY;s++)
buffer[s]=0;
}
inline float feed(float const in,float const feedback,double const delay)
{
// calculate delay offset
double back=(double)counter-delay;
// clip lookback buffer-bound
if(back<0.0)
back=MAX_WG_DELAY+back;
// compute interpolation left-floor
int const index0=floor_int(back);
// compute interpolation right-floor
int index_1=index0-1;
int index1=index0+1;
int index2=index0+2;
// clip interp. buffer-bound
if(index_1<0)index_1=MAX_WG_DELAY-1;
if(index1>=MAX_WG_DELAY)index1=0;
if(index2>=MAX_WG_DELAY)index2=0;
// get neighbourgh samples
float const y_1= buffer [index_1];
float const y0 = buffer [index0];
float const y1 = buffer [index1];
float const y2 = buffer [index2];
// compute interpolation x
float const x=(float)back-(float)index0;
// calculate
float const c0 = y0;
float const c1 = 0.5f*(y1-y_1);
float const c2 = y_1 - 2.5f*y0 + 2.0f*y1 - 0.5f*y2;
float const c3 = 0.5f*(y2-y_1) + 1.5f*(y0-y1);
float const output=((c3*x+c2)*x+c1)*x+c0;
// add to delay buffer
buffer[counter]=in+output*feedback;
// increment delay counter
counter++;
// clip delay counter
if(counter>=MAX_WG_DELAY)
counter=0;
// return output
return output;
}
float buffer[MAX_WG_DELAY];
int counter;
};
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Clipping without branching
Type : Min, max and clip
References : Posted by Laurent de Soras
Notes :
It
may reduce accuracy for small numbers. I.e. if you clip to [-1; 1],
fractional part of the result will be quantized to 23 bits (or more,
depending on the bit depth of the temporary results). Thus, 1e-20 will
be rounded to 0. The other (positive) side effect is the denormal
number elimination.
Code :
float max (float x, float a)
{
x -= a;
x += fabs (x);
x *= 0.5;
x += a;
return (x);
}
float min (float x, float b)
{
x = b - x;
x += fabs (x)
x *= 0.5;
x = b - x;
return (x);
}
float clip (float x, float a, float b)
{
x1 = fabs (x-a);
x2 = fabs (x-b);
x = x1 + (a+b);
x -= x2;
x *= 0.5;
return (x);
}
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Compressor
Type : Hardknee compressor with RMS look-ahead envelope calculation and adjustable attack/decay
References : Posted by flashinc[AT]mail[DOT]ru
Notes :
RMS is a true way to estimate _musical_ signal energy,
our ears behaves in a same way.
to making all it work,
try this values (as is, routine accepts percents and milliseconds) for first time:
threshold = 50%
slope = 50%
RMS window width = 1 ms
lookahead = 3 ms
attack time = 0.1 ms
release time = 300 ms
This code can be significantly improved in speed by
changing RMS calculation loop to 'running summ'
(keeping the summ in 'window' -
adding next newest sample and subtracting oldest on each step)
Code :
void compress
(
float* wav_in, // signal
int n, // N samples
double threshold, // threshold (percents)
double slope, // slope angle (percents)
int sr, // sample rate (smp/sec)
double tla, // lookahead (ms)
double twnd, // window time (ms)
double tatt, // attack time (ms)
double trel // release time (ms)
)
{
typedef float stereodata[2];
stereodata* wav = (stereodata*) wav_in; // our stereo signal
threshold *= 0.01; // threshold to unity (0...1)
slope *= 0.01; // slope to unity
tla *= 1e-3; // lookahead time to seconds
twnd *= 1e-3; // window time to seconds
tatt *= 1e-3; // attack time to seconds
trel *= 1e-3; // release time to seconds
// attack and release "per sample decay"
double att = (tatt == 0.0) ? (0.0) : exp (-1.0 / (sr * tatt));
double rel = (trel == 0.0) ? (0.0) : exp (-1.0 / (sr * trel));
// envelope
double env = 0.0;
// sample offset to lookahead wnd start
int lhsmp = (int) (sr * tla);
// samples count in lookahead window
int nrms = (int) (sr * twnd);
// for each sample...
for (int i = 0; i < n; ++i)
{
// now compute RMS
double summ = 0;
// for each sample in window
for (int j = 0; j < nrms; ++j)
{
int lki = i + j + lhsmp;
double smp;
// if we in bounds of signal?
// if so, convert to mono
if (lki < n)
smp = 0.5 * wav[lki][0] + 0.5 * wav[lki][1];
else
smp = 0.0; // if we out of bounds we just get zero in smp
summ += smp * smp; // square em..
}
double rms = sqrt (summ / nrms); // root-mean-square
// dynamic selection: attack or release?
double theta = rms > env ? att : rel;
// smoothing with capacitor, envelope extraction...
// here be aware of pIV denormal numbers glitch
env = (1.0 - theta) * rms + theta * env;
// the very easy hard knee 1:N compressor
double gain = 1.0;
if (env > threshold)
gain = gain - (env - threshold) * slope;
// result - two hard kneed compressed channels...
float leftchannel = wav[i][0] * gain;
float rightchannel = wav[i][1] * gain;
}
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Constant-time exponent of 2 detector
References : Posted by Brent Lehman (mailbjl[AT]yahoo.com)
Notes :
In
your common FFT program, you want to make sure that the frame you're
working with has a size that is a power of 2. This tells you in just a
few operations. Granted, you won't be using this algorithm inside a
loop, so the savings aren't that great, but every little hack helps ;)
Code :
// Quit if size isn't a power of 2
if ((-size ^ size) & size) return;
// If size is an unsigned int, the above might not compile.
// You'd want to use this instead:
if (((~size + 1) ^ size) & size) return;
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Conversions on a PowerPC
Type : motorola ASM conversions
References : Posted by James McCartney
Code :
double ftod(float x) { return (double)x;
00000000: 4E800020 blr
// blr == return from subroutine, i.e. this function is a noop
float dtof(double x) { return (float)x;
00000000: FC200818 frsp fp1,fp1
00000004: 4E800020 blr
int ftoi(float x) { return (int)x;
00000000: FC00081E fctiwz fp0,fp1
00000004: D801FFF0 stfd fp0,-16(SP)
00000008: 8061FFF4 lwz r3,-12(SP)
0000000C: 4E800020 blr
int dtoi(double x) { return (int)x;
00000000: FC00081E fctiwz fp0,fp1
00000004: D801FFF0 stfd fp0,-16(SP)
00000008: 8061FFF4 lwz r3,-12(SP)
0000000C: 4E800020 blr
double itod(int x) { return (double)x;
00000000: C8220000 lfd fp1,@1558(RTOC)
00000004: 6C608000 xoris r0,r3,$8000
00000008: 9001FFF4 stw r0,-12(SP)
0000000C: 3C004330 lis r0,17200
00000010: 9001FFF0 stw r0,-16(SP)
00000014: C801FFF0 lfd fp0,-16(SP)
00000018: FC200828 fsub fp1,fp0,fp1
0000001C: 4E800020 blr
float itof(int x) { return (float)x;
00000000: C8220000 lfd fp1,@1558(RTOC)
00000004: 6C608000 xoris r0,r3,$8000
00000008: 9001FFF4 stw r0,-12(SP)
0000000C: 3C004330 lis r0,17200
00000010: 9001FFF0 stw r0,-16(SP)
00000014: C801FFF0 lfd fp0,-16(SP)
00000018: EC200828 fsubs fp1,fp0,fp1
0000001C: 4E800020 blr
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cool Sounding Lowpass With Decibel Measured Resonance
Type : LP 2-pole resonant tweaked butterworth
References : Posted by daniel_jacob_werner [AT] yaho [DOT] com [DOT] au
Notes :
This
algorithm is a modified version of the tweaked butterworth lowpass
filter by Patrice Tarrabia posted on musicdsp.org's archives. It
calculates the coefficients for a second order IIR filter. The
resonance is specified in decibels above the DC gain. It can be made
suitable to use as a SoundFont 2.0 filter by scaling the output so the
overall gain matches the specification (i.e. if resonance is 6dB then
you should scale the output by -3dB). Note that you can replace the
sqrt(2) values in the standard butterworth highpass algorithm with my
"q =" line of code to get a highpass also. How it works: normally q is
the constant sqrt(2), and this value controls resonance. At sqrt(2)
resonance is 0dB, smaller values increase resonance. By multiplying
sqrt(2) by a power ratio we can specify the resonant gain at the cutoff
frequency. The resonance power ratio is calculated with a standard
formula to convert between decibels and power ratios (the powf
statement...).
Good Luck,
Daniel Werner
http://experimentalscene.com/
Code :
float c, csq, resonance, q, a0, a1, a2, b1, b2;
c = 1.0f / (tanf(pi * (cutoff / samplerate)));
csq = c * c;
resonance = powf(10.0f, -(resonancedB * 0.1f));
q = sqrt(2.0f) * resonance;
a0 = 1.0f / (1.0f + (q * c) + (csq));
a1 = 2.0f * a0;
a2 = a0;
b1 = (2.0f * a0) * (1.0f - csq);
b2 = a0 * (1.0f - (q * c) + csq);
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Copy-protection schemes
References : Posted by Moyer, Andy
Notes :
This post of Andy sums up everything there is to know about copy-protection schemes:
"Build a great product and release improvements regularly so that people will
be willing to spend the money on it, thus causing anything that is cracked
to be outdated quickly. Build a strong relationship with your customers,
because if they've already paid for one of your products, and were
satisfied, chances are, they will be more likely to buy another one of your
products. Make your copy protection good enough so that somebody can't just
do a search in Google and enter in a published serial number, but don't make
registered users jump through flaming hoops to be able to use the product.
Also use various approaches to copy protection within a release, and vary
those approaches over multiple releases so that a hacker that cracked your
app's version 1.0 can't just run a recorded macro in a text editor to crack
your version 2.0 software [this being simplified]."
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cubic interpollation
Type : interpollation
References : Posted by Olli Niemitalo
Linked file : other001.gif
Notes :
(see linkfile)
finpos is the fractional, inpos the integer part.
Code :
xm1 = x [inpos - 1];
x0 = x [inpos + 0];
x1 = x [inpos + 1];
x2 = x [inpos + 2];
a = (3 * (x0-x1) - xm1 + x2) / 2;
b = 2*x1 + xm1 - (5*x0 + x2) / 2;
c = (x1 - xm1) / 2;
y [outpos] = (((a * finpos) + b) * finpos + c) * finpos + x0;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cubic polynomial envelopes
Type : envellope generation
References : Posted by Andy Mucho
Notes :
This function runs from:
startlevel at Time=0
midlevel at Time/2
endlevel at Time
At moments of extreme change over small time, the function can generate out
of range (of the 3 input level) numbers, but isn't really a problem in
actual use with real numbers, and sensible/real times..
Code :
time = 32
startlevel = 0
midlevel = 100
endlevel = 120
k = startlevel + endlevel - (midlevel * 2)
r = startlevel
s = (endlevel - startlevel - (2 * k)) / time
t = (2 * k) / (time * time)
bigr = r
bigs = s + t
bigt = 2 * t
for(int i=0;i<time;i++)
{
bigr = bigr + bigs
bigs = bigs + bigt
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
DC filter
Type : 1-pole/1-zero DC filter
References : Posted by andy[DOT]rossol[AT]bluewin[DOT]ch
Notes :
This is based on code found in the document:
"Introduction to Digital Filters (DRAFT)"
Julius O. Smith III (jos@ccrma.stanford.edu)
(http://www-ccrma.stanford.edu/~jos/filters/)
---
Some audio algorithms (asymmetric waveshaping, cascaded filters, ...)
can produce DC offset. This offset can accumulate and reduce the
signal/noise ratio.
So, how to fix it? The example code from Julius O. Smith's document is:
...
y(n) = x(n) - x(n-1) + R * y(n-1)
// "R" between 0.9 .. 1
// n=current (n-1)=previous in/out value
...
"R" depends on sampling rate and the low frequency point. Do not set
"R" to a fixed value (e.g. 0.99) if you don't know the sample rate.
Instead set R to:
(-3dB @ 40Hz): R = 1-(250/samplerate)
(-3dB @ 30Hz): R = 1-(190/samplerate)
(-3dB @ 20Hz): R = 1-(126/samplerate)
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Decimator
Type : Bit-reducer and sample&hold unit
References : Posted by tobyear[AT]web[DOT]de
Notes :
This
is a simple bit and sample rate reduction code, maybe some of you can
use it. The parameters are bits (1..32) and rate (0..1, 1 is the
original samplerate).
Call the function like this:
y=decimate(x);
A VST plugin implementing this algorithm (with full Delphi source code included) can be downloaded from here:
http://tobybear.phreque.com/decimator.zip
Comments/suggestions/improvements are welcome, send them to: tobybear@web.de
Code :
// bits: 1..32
// rate: 0..1 (1 is original samplerate)
********** Pascal source **********
var m:longint;
y,cnt,rate:single;
// call this at least once before calling
// decimate() the first time
procedure setparams(bits:integer;shrate:single);
begin
m:=1 shl (bits-1);
cnt:=1;
rate:=shrate;
end;
function decimate(i:single):single;
begin
cnt:=cnt+rate;
if (cnt>1) then
begin
cnt:=cnt-1;
y:=round(i*m)/m;
end;
result:=y;
end;
********** C source **********
int bits=16;
float rate=0.5;
long int m=1<<(bits-1);
float y=0,cnt=0;
float decimate(float i)
{
cnt+=rate;
if (cnt>=1)
{
cnt-=1;
y=(long int)(i*m)/(float)m;
}
return y;
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Delay time calculation for reverberation
References : Posted by Andy Mucho
Notes :
This is from some notes I had scribbled down from a while back on
automatically calculating diffuse delays. Given an intial delay line gain
and time, calculate the times and feedback gain for numlines delay lines..
Code :
int numlines = 8;
float t1 = 50.0; // d0 time
float g1 = 0.75; // d0 gain
float rev = -3*t1 / log10 (g1);
for (int n = 0; n < numlines; ++n)
{
float dt = t1 / pow (2, (float (n) / numlines));
float g = pow (10, -((3*dt) / rev));
printf ("d%d t=%.3f g=%.3f\n", n, dt, g);
}
The above with t1=50.0 and g1=0.75 yields:
d0 t=50.000 g=0.750
d1 t=45.850 g=0.768
d2 t=42.045 g=0.785
d3 t=38.555 g=0.801
d4 t=35.355 g=0.816
d5 t=32.421 g=0.830
d6 t=29.730 g=0.843
d7 t=27.263 g=0.855
To go more diffuse, chuck in dual feedback paths with a one cycle delay
effectively creating a phase-shifter in the feedback path, then things get
more exciting.. Though what the optimum phase shifts would be I couldn't
tell you right now..
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Denormal DOUBLE variables, macro
References : Posted by Jon Watte
Notes :
Use this macro if you want to find denormal numbers and you're using doubles...
Code :
#if PLATFORM_IS_BIG_ENDIAN
#define INDEX 0
#else
#define INDEX 1
#endif
inline bool is_denormal( double const & d ) {
assert( sizeof( d ) == 2*sizeof( int ) );
int l = ((int *)&d)[INDEX];
return (l&0x7fe00000) != 0;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Denormal numbers, the meta-text
References : Laurent de Soras
Linked file : denormal.pdf
Notes :
This
very interesting paper, written by Laurent de Soras (www.ohmforce.com)
has everything you ever wanted to know about denormal numbers! And it
obviously descibes how you can get rid of them too!
(see linked file)
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
DFT
Type : fourier transform
References : Posted by Andy Mucho
Code :
AnalyseWaveform(float *waveform, int framesize)
{
float aa[MaxPartials];
float bb[MaxPartials];
for(int i=0;i<partials;i++)
{
aa[i]=0;
bb[i]=0;
}
int hfs=framesize/2;
float pd=pi/hfs;
for (i=0;i<framesize;i++)
{
float w=waveform[i];
int im = i-hfs;
for(int h=0;h<partials;h++)
{
float th=(pd*(h+1))*im;
aa[h]+=w*cos(th);
bb[h]+=w*sin(th);
}
}
for (int h=0;h<partials;h++)
amp[h]= sqrt(aa[h]*aa[h]+bb[h]*bb[h])/hfs;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Digital RIAA equalization filter coefficients
Type : RIAA
References : Posted by Frederick Umminger
Notes :
Use at your own risk. Confirm correctness before using. Don't assume I didn't goof something up.
-Frederick Umminger
Code :
The
"turntable-input software" thread inspired me to generate some
coefficients for a digital RIAA equalization filter. These coefficients
were found by matching the magnitude response of the s-domain transfer
function using some proprietary Matlab scripts. The phase response may
or may not be totally whacked.
The s-domain transfer function is
R3(1+R1*C1*s)(1+R2*C2*s)/(R1(1+R2*C2*s) + R2(1+R1*C1*s) + R3(1+R1*C1*s)(1+R2*C2*s))
where
R1 = 883.3k
R2 = 75k
R3 = 604
C1 = 3.6n
C2 = 1n
This is based on the reference circuit found in http://www.hagtech.com/pdf/riaa.pdf
The coefficients of the digital transfer function b(z^-1)/a(z^-1) in descending powers of z, are:
44.1kHz
b = [ 0.02675918611906 -0.04592084787595 0.01921229297239]
a = [ 1.00000000000000 -0.73845850035973 -0.17951755477430]
error +/- 0.25dB
48kHz
b = [ 0.02675918611906 -0.04592084787595 0.01921229297239]
a = [ 1.00000000000000 -0.73845850035973 -0.17951755477430]
error +/- 0.15dB
88.2kHz
b = [ 0.04872204977233 -0.09076930609195 0.04202280710877]
a = [ 1.00000000000000 -0.85197860443215 -0.10921171201431]
error +/- 0.01dB
96kHz
b = [ 0.05265477122714 -0.09864197097385 0.04596474352090 ]
a = [ 1.00000000000000 -0.85835597216218 -0.10600020417219 ]
error +/- 0.006dB
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Direct form II
Type : generic
References : Posted by Fuzzpilz
Notes :
I've
noticed there's no code for direct form II filters in general here,
though probably many of the filter examples use it. I haven't looked at
them all to verify that, but there certainly doesn't seem to be a
snippet describing this.
This is a simple direct form II implementation of a k-pole, k-zero
filter. It's a little faster than (a naive, real-time implementation
of) direct form I, as well as more numerically accurate.
Code :
Direct form I pseudocode:
y[n] = a[0]*x[n] + a[1]*x[n-1] + .. + a[k]*x[n-k]
- b[1]*y[n-1] - .. - b[k]*y[n-k];
Simple equivalent direct form II pseudocode:
y[n] = a[0]*x[n] + d[0];
d[0] = a[1]*x[n] - b[1]*y[n] + d[1];
d[1] = a[2]*x[n] - b[2]*y[n] + d[2];
.
.
d[k-2] = a[k-1]*x[n] - b[k-1]*y[n] + d[k-1];
d[k-1] = a[k]*x[n] - b[k]*y[n];
For example, a biquad:
out = a0*in + a1*h0 + a2*h1 - b1*h2 - b2*h3;
h1 = h0;
h0 = in;
h3 = h2;
h2 = out;
becomes
out = a0*in + d0;
d0 = a1*in - b1*out + d1;
d1 = a2*in - b2*out;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Discrete Summation Formula (DSF)
References : Stylson, Smith and others... (posted by Alexander Kritov)
Notes :
Buzz uses this type of synth.
For cool sounds try to use variable,
for example a=exp(-x/12000)*0.8 // x- num.samples
Code :
double DSF (double x, // input
double a, // a<1.0
double N, // N<SmplFQ/2,
double fi) // phase
{
double s1 = pow(a,N-1.0)*sin((N-1.0)*x+fi);
double s2 = pow(a,N)*sin(N*x+fi);
double s3 = a*sin(x+fi);
double s4 =1.0 - (2*a*cos(x)) +(a*a);
if (s4==0)
return 0;
else
return (sin(fi) - s3 - s2 +s1)/s4;
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Dither code
Type : Dither with noise-shaping
References : Posted by Paul Kellett
Notes :
This
is a simple implementation of highpass triangular-PDF dither (a good
general-purpose dither) with optional 2nd-order noise shaping (which
lowers the noise floor by 11dB below 0.1 Fs).
The code assumes input data is in the range +1 to -1 and doesn't check for overloads!
To save time when generating dither for multiple channels you can
re-use lower bits of a previous random number instead of calling rand()
again. e.g. r3=(r1 & 0x7F)<<8;
Code :
int r1, r2; //rectangular-PDF random numbers
float s1, s2; //error feedback buffers
float s = 0.5f; //set to 0.0f for no noise shaping
float w = pow(2.0,bits-1); //word length (usually bits=16)
float wi= 1.0f/w;
float d = wi / RAND_MAX; //dither amplitude (2 lsb)
float o = wi * 0.5f; //remove dc offset
float in, tmp;
int out;
//for each sample...
r2=r1; //can make HP-TRI dither by
r1=rand(); //subtracting previous rand()
in += s * (s1 + s1 - s2); //error feedback
tmp = in + o + d * (float)(r1 - r2); //dc offset and dither
out = (int)(w * tmp); //truncate downwards
if(tmp<0.0f) out--; //this is faster than floor()
s2 = s1;
s1 = in - wi * (float)out; //error
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Double to Int
Type : pointer cast (round to zero, or 'trunctate')
References : Posted by many people, implementation by Andy M00cho
Notes :
-Platform
independant, literally. You have IEEE FP numbers, this will work, as
long as your not expecting a signed integer back larger than 16bits :)
-Will only work correctly for FP numbers within the range of [-32768.0,32767.0]
-The FPU must be in Double-Precision mode
Code :
typedef double lreal;
typedef float real;
typedef unsigned long uint32;
typedef long int32;
//2^36 * 1.5, (52-_shiftamt=36) uses limited precision to floor
//16.16 fixed point representation
const lreal _double2fixmagic = 68719476736.0*1.5;
const int32 _shiftamt = 16;
#if BigEndian_
#define iexp_ 0
#define iman_ 1
#else
#define iexp_ 1
#define iman_ 0
#endif //BigEndian_
// Real2Int
inline int32 Real2Int(lreal val)
{
val= val + _double2fixmagic;
return ((int32*)&val)[iman_] >> _shiftamt;
}
// Real2Int
inline int32 Real2Int(real val)
{
return Real2Int ((lreal)val);
}
For the x86 assembler freaks here's the assembler equivalent:
__double2fixmagic dd 000000000h,042380000h
fld AFloatingPoint Number
fadd QWORD PTR __double2fixmagic
fstp TEMP
movsx eax,TEMP+2
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Drift generator
Type : Random
References : Posted by quintosardo[AT]yahoo[DOT]it
Notes :
I use this drift to modulate any sound parameter of my synth.
It is very effective if it slightly modulates amplitude or frequency of an FM modulator.
It is based on an incremental random variable, sine-warped.
I like it because it is "continuous" (as opposite to "sample and hold"), and I can set variation rate and max variation.
It can go to upper or lower constraint (+/- max drift) but it gradually
decreases rate of variation when approaching to the limit.
I use it exactly as an LFO (-1.f .. +1.f)
I use a table for sin instead of sin() function because this way I can
change random distribution, by selecting a different curve (different
table) from sine...
I hope that it is clear ... (sigh... :-)
Bye!!!
P.S. Thank you for help in previous submission ;-)
Code :
const int kSamples //Number of samples in fSinTable below
float fSinTable[kSamples] // Tabulated sin() [0 - 2pi[ amplitude [-1.f .. 1.f]
float fWhere// Index
float fRate // Max rate of variation
float fLimit //max or min value
float fDrift // Output
//I assume that random() is a number from 0.f to 1.f, otherwise scale it
fWhere += fRate * random()
//I update this drift in a long-term cycle, so I don't care of branches
if (fWhere >= 1.f) fWhere -= 1.f
else if (fWhere < 0.f) sWhere += 1.f
fDrift = fLimit * fSinTable[(long) (fWhere * kSamples)]
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
DSF (super-set of BLIT)
Type : matlab code
References : Posted by David Lowenfels
Notes :
Discrete Summation Formula ala Moorer
computes equivalent to sum{k=0:N-1}(a^k * sin(beta + k*theta))
modified from Emanuel Landeholm's C code
output should never clip past [-1,1]
If using for BLIT synthesis for virtual analog:
N = maxN;
a = attn_at_Nyquist ^ (1/maxN); %hide top harmonic popping in and out when sweeping frequency
beta = pi/2;
num = 1 - a^N * cos(N*theta) - a*( cos(theta) - a^N * cos(N*theta - theta) ); %don't waste time on beta
You can also get growing harmonics if a > 1, but the min statement
in the code must be removed, and the scaling will be weird.
Code :
function output = dsf( freq, a, H, samples, beta)
%a = rolloff coeffecient
%H = number of harmonic overtones (fundamental not included)
%beta = harmonic phase shift
samplerate = 44.1e3;
freq = freq/samplerate; %normalize frequency
bandlimit = samplerate / 2; %Nyquist
maxN = 1 + floor( bandlimit / freq ); %prevent aliasing
N = min(H+2,maxN);
theta = 2*pi * phasor(freq, samples);
epsilon = 1e-6;
a = min(a, 1-epsilon); %prevent divide by zero
num = sin(beta) - a*sin(beta-theta) - a^N*sin(beta + N*theta) + a^(N+1)*sin(beta+(N-1)*theta);
den = (1 + a * ( a - 2*cos(theta) ));
output = 2*(num ./ den - 1) * freq; %subtract by one to remove DC, scale by freq to normalize
output = output * maxN/N; %OPTIONAL: rescale to give louder output as rolloff increases
function out = phasor(normfreq, samples);
out = mod( (0:samples-1)*normfreq , 1);
out = out * 2 - 1; %make bipolar
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Envelope detector
References : Posted by Bram
Notes :
Basicaly
a one-pole LP filter with different coefficients for attack and release
fed by the abs() of the signal. If you don't need different attack and
decay settings, just use in->abs()->LP
Code :
//attack and release in milliseconds
float ga = (float) exp(-1/(SampleRate*attack));
float gr = (float) exp(-1/(SampleRate*release));
float envelope=0;
for(...)
{
//get your data into 'input'
EnvIn = abs(input);
if(envelope < EnvIn)
{
envelope *= ga;
envelope += (1-ga)*EnvIn;
}
else
{
envelope *= gr;
envelope += (1-gr)*EnvIn;
}
//envelope now contains.........the envelope ;)
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Envelope Follower
References : Posted by ers
Code :
#define V_ENVELOPE_FOLLOWER_NUM_POINTS 2000
class vEnvelopeFollower :
{
public:
vEnvelopeFollower();
virtual ~vEnvelopeFollower();
inline void Calculate(float *b)
{
envelopeVal -= *buff;
if (*b < 0)
envelopeVal += *buff = -*b;
else
envelopeVal += *buff = *b;
if (buff++ == bufferEnd)
buff = buffer;
}
void SetBufferSize(float value);
void GetControlValue(){return envelopeVal / (float)bufferSize;}
private:
float buffer[V_ENVELOPE_FOLLOWER_NUM_POINTS];
float *bufferEnd, *buff, envelopeVal;
int bufferSize;
float val;
};
vEnvelopeFollower::vEnvelopeFollower()
{
bufferEnd = buffer + V_ENVELOPE_FOLLOWER_NUM_POINTS-1;
buff = buffer;
val = 0;
float *b = buffer;
do
{
*b++ = 0;
}while (b <= bufferEnd);
bufferSize = V_ENVELOPE_FOLLOWER_NUM_POINTS;
envelopeVal= 0;
}
vEnvelopeFollower::~vEnvelopeFollower()
{
}
void vEnvelopeFollower::SetBufferSize(float value)
{
bufferEnd = buffer + (bufferSize = 100 + (int)(value * ((float)V_ENVELOPE_FOLLOWER_NUM_POINTS-102)));
buff = buffer;
float val = envelopeVal / bufferSize;
do
{
*buff++ = val;
}while (buff <= bufferEnd);
buff = buffer;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Envelope follower with different attack and release
References : Posted by Bram
Notes :
xxxx_in_ms is xxxx in milliseconds ;-)
Code :
init::
attack_coef = exp(log(0.01)/( attack_in_ms * samplerate * 0.001));
release_coef = exp(log(0.01)/( release_in_ms * samplerate * 0.001));
envelope = 0.0;
loop::
tmp = fabs(in);
if(tmp > envelope)
envelope = attack_coef * (envelope - tmp) + tmp;
else
envelope = release_coef * (envelope - tmp) + tmp;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Exponential parameter mapping
References : Posted by Russell Borogove
Notes :
Use this if you want to do an exponential map of a parameter (mParam) to a range (mMin - mMax).
Output is in mData...
Code :
float logmax = log10f( mMax );
float logmin = log10f( mMin );
float logdata = (mParam * (logmax-logmin)) + logmin;
mData = powf( 10.0f, logdata );
if (mData < mMin)
{
mData = mMin;
}
if (mData > mMax)
{
mData = mMax;
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
fast abs/neg/sign for 32bit floats
Type : floating point functions
References : Posted by tobybear[AT]web[DOT]de
Notes :
Haven't
seen this elsewhere, probably because it is too obvious? Anyway, these
functions are intended for 32-bit floating point numbers only and
should work a bit faster than the regular ones.
fastabs() gives you the absolute value of a float
fastneg() gives you the negative number (faster than multiplying with -1)
fastsgn() gives back +1 for 0 or positive numbers, -1 for negative numbers
Comments are welcome (tobybear[AT]web[DOT]de)
Cheers
Toby (www.tobybear.de)
Code :
// C/C++ code:
float fastabs(float f)
{int i=((*(int*)&f)&0x7fffffff);return (*(float*)&i);}
float fastneg(float f)
{int i=((*(int*)&f)^0x80000000);return (*(float*)&i);}
int fastsgn(float f)
{return 1+(((*(int*)&f)>>31)<<1);}
//Delphi/Pascal code:
function fastabs(f:single):single;
begin i:=longint((@f)^) and $7FFFFFFF;result:=single((@i)^) end;
function fastneg(f:single):single;
begin i:=longint((@f)^) xor $80000000;result:=single((@i)^) end;
function fastsgn(f:single):longint;
begin result:=1+((longint((@f)^) shr 31)shl 1) end;
5 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast binary log approximations
Type : C code
References : Posted by musicdsp.org[AT]mindcontrol.org
Notes :
This
code uses IEEE 32-bit floating point representation knowledge to
quickly compute approximations to the log2 of a value. Both functions
return under-estimates of the actual value, although the second flavour
is less of an under-estimate than the first (and might be sufficient
for using in, say, a dBV/FS level meter).
Running the test program, here's the output:
0.1: -4 -3.400000
1: 0 0.000000
2: 1 1.000000
5: 2 2.250000
100: 6 6.562500
Code :
// Fast logarithm (2-based) approximation
// by Jon Watte
#include <assert.h>
int floorOfLn2( float f ) {
assert( f > 0. );
assert( sizeof(f) == sizeof(int) );
assert( sizeof(f) == 4 );
return (((*(int *)&f)&0x7f800000)>>23)-0x7f;
}
float approxLn2( float f ) {
assert( f > 0. );
assert( sizeof(f) == sizeof(int) );
assert( sizeof(f) == 4 );
int i = (*(int *)&f);
return (((i&0x7f800000)>>23)-0x7f)+(i&0x007fffff)/(float)0x800000;
}
// Here's a test program:
#include <stdio.h>
// insert code from above here
int
main()
{
printf( "0.1: %d %f\n", floorOfLn2( 0.1 ), approxLn2( 0.1 ) );
printf( "1: %d %f\n", floorOfLn2( 1. ), approxLn2( 1. ) );
printf( "2: %d %f\n", floorOfLn2( 2. ), approxLn2( 2. ) );
printf( "5: %d %f\n", floorOfLn2( 5. ), approxLn2( 5. ) );
printf( "100: %d %f\n", floorOfLn2( 100. ), approxLn2( 100. ) );
return 0;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast exp2 approximation
References : Posted by Laurent de Soras
Notes :
Partial approximation of exp2 in fixed-point arithmetic. It is exactly :
[0 ; 1[ -> [0.5 ; 1[
f : x |-> 2^(x-1)
To get the full exp2 function, you have to separate the integer and
fractionnal part of the number. The integer part may be processed by
bitshifting. Process the fractionnal part with the function, and
multiply the two results.
Maximum error is only 0.3 % which is pretty good for two mul ! You get also the continuity of the first derivate.
-- Laurent
Code :
// val is a 16-bit fixed-point value in 0x0 - 0xFFFF ([0 ; 1[)
// Returns a 32-bit fixed-point value in 0x80000000 - 0xFFFFFFFF ([0.5 ; 1[)
unsigned int fast_partial_exp2 (int val)
{
unsigned int result;
__asm
{
mov eax, val
shl eax, 15 ; eax = input [31/31 bits]
or eax, 080000000h ; eax = input + 1 [32/31 bits]
mul eax
mov eax, edx ; eax = (input + 1) ^ 2 [32/30 bits]
mov edx, 2863311531 ; 2/3 [32/32 bits], rounded to +oo
mul edx ; eax = 2/3 (input + 1) ^ 2 [32/30 bits]
add edx, 1431655766 ; + 4/3 [32/30 bits] + 1
mov result, edx
}
return (result);
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast in-place Walsh-Hadamard Transform
Type : wavelet transform
References : Posted by Timo H Tossavainen
Notes :
IIRC, They're also called walsh-hadamard transforms.
Basically like Fourier, but the basis functions are squarewaves with different sequencies.
I did this for a transform data compression study a while back.
Here's some code to do a walsh hadamard transform on long ints in-place
(you need to divide by n to get transform) the order is bit-reversed at
output, IIRC.
The inverse transform is the same as the forward transform (expects
bit-reversed input). i.e. x = 1/n * FWHT(FWHT(x)) (x is a vector)
Code :
void inline wht_bfly (long& a, long& b)
{
long tmp = a;
a += b;
b = tmp - b;
}
// just a integer log2
int inline l2 (long x)
{
int l2;
for (l2 = 0; x > 0; x >>=1)
{
++ l2;
}
return (l2);
}
////////////////////////////////////////////
// Fast in-place Walsh-Hadamard Transform //
////////////////////////////////////////////
void FWHT (std::vector& data)
{
const int log2 = l2 (data.size()) - 1;
for (int i = 0; i < log2; ++i)
{
for (int j = 0; j < (1 << log2); j += 1 << (i+1))
{
for (int k = 0; k < (1<<i); ++k)
{
wht_bfly (data [j + k], data [j + k + (1<<i)]);
}
}
}
}
5 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast LFO in Delphi...
References : Posted by Dambrin Didier ( gol [AT] e-officedirect [DOT] com )
Linked file : LFOGenerator.zip
Notes :
[from Didier's mail...]
[see attached zip file too!]
I was working on a flanger, & needed an LFO for it. I first used a
Sin(), but it was too slow, then tried a big wavetable, but it wasn't
accurate enough.
I then checked the alternate sine generators from your web site, &
while they're good, they all can drift, so you're also wasting too much
CPU in branching for the drift checks.
So I made a quick & easy linear LFO, then a sine-like version of it. Can be useful for LFO's, not to output as sound.
If has no branching & is rather simple. 2 Abs() but apparently they're fast. In all cases faster than a Sin()
It's in delphi, but if you understand it you can translate it if you want.
It uses a 32bit integer counter that overflows, & a power for the sine output.
If you don't know delphi, $ is for hex (h at the end in c++?), Single
is 32bit float, integer is 32bit integer (signed, normally).
Code :
unit Unit1;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, ExtCtrls, ComCtrls;
type
TForm1 = class(TForm)
PaintBox1: TPaintBox;
Bevel1: TBevel;
procedure PaintBox1Paint(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;
var
Form1: TForm1;
implementation
{$R *.DFM}
procedure TForm1.PaintBox1Paint(Sender: TObject);
var n,Pos,Speed:Integer;
Output,Scale,HalfScale,PosMul:Single;
OurSpeed,OurScale:Single;
begin
OurSpeed:=100; // 100 samples per cycle
OurScale:=100; // output in -100..100
Pos:=0; // position in our linear LFO
Speed:=Round($100000000/OurSpeed);
// --- triangle LFO ---
Scale:=OurScale*2;
PosMul:=Scale/$80000000;
// loop
for n:=0 to 299 do
Begin
// inc our 32bit integer LFO pos & let it overflow. It will be seen as signed when read by the math unit
Pos:=Pos+Speed;
Output:=Abs(Pos*PosMul)-OurScale;
// visual
Paintbox1.Canvas.Pixels[n,Round(100+Output)]:=clRed;
End;
// --- sine-like LFO ---
Scale:=Sqrt(OurScale*4);
PosMul:=Scale/$80000000;
HalfScale:=Scale/2;
// loop
for n:=0 to 299 do
Begin
// inc our 32bit integer LFO pos & let it overflow. It will be seen as signed when read by the math unit
Pos:=Pos+Speed;
Output:=Abs(Pos*PosMul)-HalfScale;
Output:=Output*(Scale-Abs(Output));
// visual
Paintbox1.Canvas.Pixels[n,Round(100+Output)]:=clBlue;
End;
end;
end.
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast log2
References : Posted by Laurent de Soras
Code :
inline float fast_log2 (float val)
{
assert (val > 0);
int * const exp_ptr = reinterpret_cast <int *> (&val);
int x = *exp_ptr;
const int log_2 = ((x >> 23) & 255) - 128;
x &= ~(255 << 23);
x += 127 << 23;
*exp_ptr = x;
return (val + log_2);
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
fast power and root estimates for 32bit floats
Type : floating point functions
References : Posted by tobybear[AT]web[DOT]de
Notes :
Original code by Stefan Stenzel (also in this archive, see "pow(x,4) approximation") - extended for more flexibility.
fastpow(f,n) gives a rather *rough* estimate of a float number f to the
power of an integer number n (y=f^n). It is fast but result can be
quite a bit off, since we directly mess with the floating point
exponent.-> use it only for getting rough estimates of the values
and where precision is not that important.
fastroot(f,n) gives the n-th root of f. Same thing concerning precision applies here.
Cheers
Toby (www.tobybear.de)
Code :
//C/C++ source code:
float fastpower(float f,int n)
{
long *lp,l;
lp=(long*)(&f);
l=*lp;l-=0x3F800000l;l<<=(n-1);l+=0x3F800000l;
*lp=l;
return f;
}
float fastroot(float f,int n)
{
long *lp,l;
lp=(long*)(&f);
l=*lp;l-=0x3F800000l;l>>=(n-1);l+=0x3F800000l;
*lp=l;
return f;
}
//Delphi/Pascal source code:
function fastpower(i:single;n:integer):single;
var l:longint;
begin
l:=longint((@i)^);
l:=l-$3F800000;l:=l shl (n-1);l:=l+$3F800000;
result:=single((@l)^);
end;
function fastroot(i:single;n:integer):single;
var l:longint;
begin
l:=longint((@i)^);
l:=l-$3F800000;l:=l shr (n-1);l:=l+$3F800000;
result:=single((@l)^);
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast sine wave calculation
Type : waveform generation
References : James McCartney in Computer Music Journal, also the Julius O. Smith paper
Notes :
(posted by Niels Gorisse)
If you change the frequency, the amplitude rises (pitch lower) or
lowers (pitch rise) a LOT I fixed the first problem by thinking about
what actually goes wrong. The answer was to recalculate the phase for
that frequency and the last value, and then continue normally.
Code :
Variables:
ip = phase of the first output sample in radians
w = freq*pi / samplerate
b1 = 2.0 * cos(w)
Init:
y1=sin(ip-w)
y2=sin(ip-2*w)
Loop:
y0 = b1*y1 - y2
y2 = y1
y1 = y0
output is in y0 (y0 = sin(ip + n*freq*pi / samplerate), n= 0, 1, 2, ... I *think*)
Later note by James McCartney:
if you unroll such a loop by 3 you can even eliminate the assigns!!
y0 = b1*y1 - y2
y2 = b1*y0 - y1
y1 = b1*y2 - y0
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Fast square wave generator
Type : NON-bandlimited osc...
References : Posted by Wolfgang (wschneider[AT]nexoft.de)
Notes :
Produces a square wave -1.0f .. +1.0f.
The resulting waveform is NOT band-limited, so it's propably of not
much use for syntheis. It's rather useful for LFOs and the like, though.
Code :
Idea: use integer overflow to avoid conditional jumps.
// init:
typedef unsigned long ui32;
float sampleRate = 44100.0f; // whatever
float freq = 440.0f; // 440 Hz
float one = 1.0f;
ui32 intOver = 0L;
ui32 intIncr = (ui32)(4294967296.0 / hostSampleRate / freq));
// loop:
(*((ui32 *)&one)) &= 0x7FFFFFFF; // mask out sign bit
(*((ui32 *)&one)) |= (intOver & 0x80000000);
intOver += intIncr;
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
FFT
References : Toth Laszlo
Linked file : rvfft.ps
Linked file : rvfft.cpp
Notes :
A
paper (postscript) and some C++ source for 4 different fft algorithms,
compiled by Toth Laszlo from the Hungarian Academy of Sciences Research
Group on Artificial Intelligence.
Toth says: "I've found that Sorensen's split-radix algorithm was the
fastest, so I use this since then (this means that you may as well
delete the other routines in my source - if you believe my results)."
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Float to int (more intel asm)
References : Posted by Laurent de Soras (via flipcode)
Notes :
[Found this on flipcode, seemed worth posting here too, hopefully Laurent will approve :) -- Ross]
Here is the code I often use. It is not a _ftol replacement, but it
provides useful functions. The current processor rounding mode should
be "to nearest" ; it is the default setting for most of the compilers.
The [fadd st, st (0) / sar i,1] trick fixes the "round to nearest even
number" behaviour. Thus, round_int (N+0.5) always returns N+1 and
floor_int function is appropriate to convert floating point numbers
into fixed point numbers.
---------------------
Laurent de Soras
Audio DSP engineer & software designer
Ohm Force - Digital audio software
http://www.ohmforce.com
Code :
inline int round_int (double x)
{
int i;
static const float round_to_nearest = 0.5f;
__asm
{
fld x
fadd st, st (0)
fadd round_to_nearest
fistp i
sar i, 1
}
return (i);
}
inline int floor_int (double x)
{
int i;
static const float round_toward_m_i = -0.5f;
__asm
{
fld x
fadd st, st (0)
fadd round_toward_m_i
fistp i
sar i, 1
}
return (i);
}
inline int ceil_int (double x)
{
int i;
static const float round_toward_p_i = -0.5f;
__asm
{
fld x
fadd st, st (0)
fsubr round_toward_p_i
fistp i
sar i, 1
}
return (-i);
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Float-to-int, coverting an array of floats
References : Posted by Stefan Stenzel
Notes :
intel only
Code :
void f2short(float *fptr,short *iptr,int n)
{
_asm {
mov ebx,n
mov esi,fptr
mov edi,iptr
lea ebx,[esi+ebx*4] ; ptr after last
mov edx,0x80008000 ; damn endianess confuses...
mov ecx,0x4b004b00 ; damn endianess confuses...
mov eax,[ebx] ; get last value
push eax
mov eax,0x4b014B01
mov [ebx],eax ; mark end
mov ax,[esi+2]
jmp startf2slp
; Pad with nops to make loop start at address divisible
; by 16 + 2, e.g. 0x01408062, don't ask why, but this
; gives best performance. Unfortumately "align 16" does
; not seem to work with my VC.
; below I noted the measured execution times for different
; nop-paddings on my Pentium Pro, 100 conversions.
; saturation: off pos neg
nop ;355 546 563 <- seems to be best
; nop ;951 547 544
; nop ;444 646 643
; nop ;444 646 643
; nop ;944 951 950
; nop ;358 447 644
; nop ;358 447 643
; nop ;358 544 643
; nop ;543 447 643
; nop ;643 447 643
; nop ;1047 546 746
; nop ;545 954 1253
; nop ;545 547 661
; nop ;544 547 746
; nop ;444 947 1147
; nop ;444 548 545
in_range:
mov eax,[esi]
xor eax,edx
saturate:
lea esi,[esi+4]
mov [edi],ax
mov ax,[esi+2]
add edi,2
startf2slp:
cmp ax,cx
je in_range
mov eax,edx
js saturate ; saturate neg -> 0x8000
dec eax ; saturate pos -> 0x7FFF
cmp esi,ebx ; end reached ?
jb saturate
pop eax
mov [ebx],eax ; restore end flag
}
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Formant filter
References : Posted by Alex
Code :
/*
Public source code by alex@smartelectronix.com
Simple example of implementation of formant filter
Vowelnum can be 0,1,2,3,4 <=> A,E,I,O,U
Good for spectral rich input like saw or square
*/
//-------------------------------------------------------------VOWEL COEFFICIENTS
const double coeff[5][11]= {
{ 8.11044e-06,
8.943665402, -36.83889529, 92.01697887, -154.337906, 181.6233289,
-151.8651235, 89.09614114, -35.10298511, 8.388101016, -0.923313471 ///A
},
{4.36215e-06,
8.90438318, -36.55179099, 91.05750846, -152.422234, 179.1170248, ///E
-149.6496211,87.78352223, -34.60687431, 8.282228154, -0.914150747
},
{ 3.33819e-06,
8.893102966, -36.49532826, 90.96543286, -152.4545478, 179.4835618,
-150.315433, 88.43409371, -34.98612086, 8.407803364, -0.932568035 ///I
},
{1.13572e-06,
8.994734087, -37.2084849, 93.22900521, -156.6929844, 184.596544, ///O
-154.3755513, 90.49663749, -35.58964535, 8.478996281, -0.929252233
},
{4.09431e-07,
8.997322763, -37.20218544, 93.11385476, -156.2530937, 183.7080141, ///U
-153.2631681, 89.59539726, -35.12454591, 8.338655623, -0.910251753
}
};
//---------------------------------------------------------------------------------
static double memory[10]={0,0,0,0,0,0,0,0,0,0};
//---------------------------------------------------------------------------------
float formant_filter(float *in, int vowelnum)
{
res= (float) ( coeff[vowelnum][0] *in +
coeff[vowelnum][1] *memory[0] +
coeff[vowelnum][2] *memory[1] +
coeff[vowelnum][3] *memory[2] +
coeff[vowelnum][4] *memory[3] +
coeff[vowelnum][5] *memory[4] +
coeff[vowelnum][6] *memory[5] +
coeff[vowelnum][7] *memory[6] +
coeff[vowelnum][8] *memory[7] +
coeff[vowelnum][9] *memory[8] +
coeff[vowelnum][10] *memory[9] );
memory[9]= memory[8];
memory[8]= memory[7];
memory[7]= memory[6];
memory[6]= memory[5];
memory[5]= memory[4];
memory[4]= memory[3];
memory[3]= memory[2];
memory[2]= memory[1];
memory[1]= memory[0];
memory[0]=(double) res;
return res;
}
8 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
frequency warped FIR lattice
Type : FIR using allpass chain
References : Posted by mail[AT]mutagene[DOT]net
Notes :
Not
at all optimized and pretty hungry in terms of arrays and overhead
(function requires two arrays containing lattice filter's internal
state and ouputs to another two arrays with their next states). In this
implementation I think you'll have to juggle taps1/newtaps in your
processing loop, alternating between one set of arrays and the other
for which to send to wfirlattice).
A frequency-warped lattice filter is just a lattice filter where every
delay has been replaced with an allpass filter. By adjusting the
allpass filters, the frequency response of the filter can be adjusted
(e.g., design an FIR that approximates some filter. Play with with
warping coefficient to "sweep" the FIR up and down without changing any
other coefficients). Much more on warped filters can be found on Aki
Harma's website ( http://www.acoustics.hut.fi/~aqi/ )
Code :
float
wfirlattice(float input, float *taps1, float *taps2, float *reflcof,
float lambda, float *newtaps1, float *newtaps2, int P)
// input is filter input
// taps1,taps2 are previous filter states (init to 0)
// reflcof are reflection coefficients. abs(reflcof) < 1 for stable filter
// lamba is warping (0 = no warping, 0.75 is close to bark scale at 44.1 kHz)
// newtaps1, newtaps2 are new filter states
// P is the order of the filter
{
float forward;
float topline;
forward = input;
topline = forward;
for (int i=0;i<P;i++)
{
newtaps2[i] = topline;
newtaps1[i] = float(lambda)*(-topline + taps1[i]) + taps2[i];
topline = newtaps1[i]+forward*(reflcof[i]);
forward += newtaps1[i]*(reflcof[i]);
taps1[i]=newtaps1[i];
taps2[i]=newtaps2[i];
}
return forward;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Gaussian dithering
Type : Dithering
References : Posted by Aleksey Vaneev (picoder[AT]mail[DOT]ru)
Notes :
It
is a more sophisticated dithering than simple RND. It gives the most
low noise floor for the whole spectrum even without noise-shaping. You
can use as big N as you can afford (it will not hurt), but 4 or 5 is
generally enough.
Code :
Basically, next value is calculated this way (for RND going from -0.5 to 0.5):
dither = (RND+RND+...+RND) / N.
\ /
\_________/
N times
If your RND goes from 0 to 1, then this code is applicable:
dither = (RND+RND+...+RND - 0.5*N) / N.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Gaussian random numbers
Type : random number generation
References : Posted by tobybear[AT]web[DOT]de
Notes :
// Gaussian random numbers
// This algorithm (adapted from "Natur als fraktale Grafik" by
// Reinhard Scholl) implements a generation method for gaussian
// distributed random numbers with mean=0 and variance=1
// (standard gaussian distribution) mapped to the range of
// -1 to +1 with the maximum at 0.
// For only positive results you might abs() the return value.
// The q variable defines the precision, with q=15 the smallest
// distance between two numbers will be 1/(2^q div 3)=1/10922
// which usually gives good results.
// Note: the random() function used is the standard random
// function from Delphi/Pascal that produces *linear*
// distributed numbers from 0 to parameter-1, the equivalent
// C function is probably rand().
Code :
const q=15;
c1=(1 shl q)-1;
c2=(c1 div 3)+1;
c3=1/c1;
function GRandom:single;
begin
result:=(2*(random(c2)+random(c2)+random(c2))-3*(c2-1))*c3;
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Gaussian White noise
References : Posted by Alexey Menshikov
Notes :
Code I use sometimes, but don't remember where I ripped it from.
- Alexey Menshikov
Code :
#define ranf() ((float) rand() / (float) RAND_MAX)
float ranfGauss (int m, float s)
{
static int pass = 0;
static float y2;
float x1, x2, w, y1;
if (pass)
{
y1 = y2;
} else {
do {
x1 = 2.0f * ranf () - 1.0f;
x2 = 2.0f * ranf () - 1.0f;
w = x1 * x1 + x2 * x2;
} while (w >= 1.0f);
w = (float)sqrt (-2.0 * log (w) / w);
y1 = x1 * w;
y2 = x2 * w;
}
pass = !pass;
return ( (y1 * s + (float) m));
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Gaussian White Noise
References : Posted by remage[AT]netposta.hu
Notes :
SOURCE:
Steven W. Smith:
The Scientist and Engineer's Guide to Digital Signal Processing
http://www.dspguide.com
Code :
#define PI 3.1415926536f
float R1 = (float) rand() / (float) RAND_MAX;
float R2 = (float) rand() / (float) RAND_MAX;
float X = (float) sqrt( -2.0f * log( R1 )) * cos( 2.0f * PI * R2 );
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Guitar feedback
References : Posted by Sean Costello
Notes :
It
is fairly simple to simulate guitar feedback with a simple
Karplus-Strong algorithm (this was described in a CMJ article in the
early 90's):
Code :
Run
the output of the Karplus-Strong delay lines into a nonlinear shaping
function for distortion (i.e. 6 parallel delay lines for 6 strings,
going into 1 nonlinear shaping function that simulates an overdriven
amplifier, fuzzbox, etc.);
Run part of the output into a delay line, to simulate the distance from the amplifier to the "strings";
The delay line feeds back into the Karplus-Strong delay lines. By
controlling the amount of the output fed into the delay line, and the
length of the delay line, you can control the intensity and pitch of
the feedback note.
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Hermite Interpolator (x86 ASM)
Type : Hermite interpolator in x86 assembly (for MS VC++)
References : Posted by robert[DOT]bielik[AT]rbcaudio[DOT]com
Notes :
An
"assemblified" variant of Laurent de Soras hermite interpolator. I
tried to do calculations as parallell as I could muster, but there is
almost certainly room for improvements. Right now, it works about 5.3
times (!) faster, not bad to start with...
Parameter explanation:
frac_pos: fractional value [0.0f - 1.0f] to interpolator
pntr: pointer to float array where:
pntr[0] = previous sample (idx = -1)
pntr[1] = current sample (idx = 0)
pntr[2] = next sample (idx = +1)
pntr[3] = after next sample (idx = +2)
The interpolation takes place between pntr[1] and pntr[2].
Regards,
/Robert Bielik
RBC Audio
Code :
const float c_half = 0.5f;
__declspec(naked) float __hermite(float frac_pos, const float* pntr)
{
__asm
{
push ecx;
mov ecx, dword ptr[esp + 12];
//////////////////////////////////////////////////////////////////////////////////////////////////
add ecx,
0x04; // ST(0) ST(1) ST(2) ST(3) ST(4) ST(5) ST(6) ST(7)
fld dword ptr [ecx+4]; // x1
fsub dword ptr [ecx-4]; // x1-xm1
fld dword ptr [ecx]; // x0 x1-xm1
fsub dword ptr [ecx+4]; // v x1-xm1
fld dword ptr [ecx+8]; // x2 v x1-xm1
fsub dword ptr [ecx]; // x2-x0 v x1-xm1
fxch st(2); // x1-m1 v x2-x0
fmul c_half; // c v x2-x0
fxch st(2); // x2-x0 v c
fmul c_half; // 0.5*(x2-x0) v c
fxch st(2); // c v 0.5*(x2-x0)
fst st(3); // c v 0.5*(x2-x0) c
fadd st(0), st(1); // w v 0.5*(x2-x0) c
fxch st(2); // 0.5*(x2-x0) v w c
faddp st(1), st(0); // v+.5(x2-x0) w c
fadd st(0), st(1); // a w c
fadd st(1), st(0); // a b_neg c
fmul dword ptr [esp+8]; // a*frac b_neg c
fsubp st(1), st(0); // a*f-b c
fmul dword ptr [esp+8]; // (a*f-b)*f c
faddp st(1), st(0); // res-x0/f
fmul dword ptr [esp+8]; // res-x0
fadd dword ptr [ecx]; // res
pop ecx;
ret;
}
}
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Hermite interpollation
References : Posted by various
Notes :
These are all different ways to do the same thing : hermite interpollation. Try'm all and benchmark.
Code :
// original
inline float hermite1(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float c2 = y0 - 2.5f * y1 + 2.f * y2 - 0.5f * y3;
float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
return ((c3 * x + c2) * x + c1) * x + c0;
}
// james mccartney
inline float hermite2(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
float c2 = y0 - y1 + c1 - c3;
return ((c3 * x + c2) * x + c1) * x + c0;
}
// james mccartney
inline float hermite3(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float y0my1 = y0 - y1;
float c3 = (y1 - y2) + 0.5f * (y3 - y0my1 - y2);
float c2 = y0my1 + c1 - c3;
return ((c3 * x + c2) * x + c1) * x + c0;
}
// laurent de soras
inline float hermite4(float frac_pos, float xm1, float x0, float x1, float x2)
{
const float c = (x1 - xm1) * 0.5f;
const float v = x0 - x1;
const float w = c + v;
const float a = w + v + (x2 - x0) * 0.5f;
const float b_neg = w + a;
return ((((a * frac_pos) - b_neg) * frac_pos + c) * frac_pos + x0);
}
7 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Inverted parabolic envelope
Type : envellope generation
References : Posted by James McCartney
Code :
dur = duration in samples
midlevel = amplitude at midpoint
beglevel = beginning and ending level (typically zero)
amp = midlevel - beglevel;
rdur = 1.0 / dur;
rdur2 = rdur * rdur;
level = beglevel;
slope = 4.0 * amp * (rdur - rdur2);
curve = -8.0 * amp * rdur2;
...
for (i=0; i<dur; ++i) {
level += slope;
slope += curve;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Java FFT
Type : FFT Analysis
References : Posted by Loreno Heer
Notes :
May not work correctly ;-)
Code :
// WTest.java
/*
Copyright (C) 2003 Loreno Heer, (helohe at bluewin dot ch)
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
public class WTest{
private static double[] sin(double step, int size){
double f = 0;
double[] ret = new double[size];
for(int i = 0; i < size; i++){
ret[i] = Math.sin(f);
f += step;
}
return ret;
}
private static double[] add(double[] a, double[] b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] + b[i];
}
return c;
}
private static double[] sub(double[] a, double[] b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] - b[i];
}
return c;
}
private static double[] add(double[] a, double b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] + b;
}
return c;
}
private static double[] cp(double[] a, int size){
double[] c = new double[size];
for(int i = 0; i < size; i++){
c[i] = a[i];
}
return c;
}
private static double[] mul(double[] a, double b){
double[] c = new double[a.length];
for(int i = 0; i < a.length; i++){
c[i] = a[i] * b;
}
return c;
}
private static void print(double[] value){
for(int i = 0; i < value.length; i++){
System.out.print(i + "," + value[i] + "\n");
}
System.out.println();
}
private static double abs(double[] a){
double c = 0;
for(int i = 0; i < a.length; i++){
c = ((c * i) + Math.abs(a[i])) / (i + 1);
}
return c;
}
private static double[] fft(double[] a, int min, int max, int step){
double[] ret = new double[(max - min) / step];
int i = 0;
for(int d = min; d < max; d = d + step){
double[] f = sin(fc(d), a.length);
double[] dif = sub(a, f);
ret[i] = 1 - abs(dif);
i++;
}
return ret;
}
private static double[] fft_log(double[] a){
double[] ret = new double[1551];
int i = 0;
for(double d = 0; d < 15.5; d = d + 0.01){
double[] f = sin(fc(Math.pow(2,d)), a.length);
double[] dif = sub(a, f);
ret[i] = Math.abs(1 - abs(dif));
i++;
}
return ret;
}
private static double fc(double d){
return d * Math.PI / res;
}
private static void print_log(double[] value){
for(int i = 0; i < value.length; i++){
System.out.print(Math.pow(2,((double)i/100d)) + "," + value[i] + "\n");
}
System.out.println();
}
public static void main(String[] args){
double[] f_0 = sin(fc(440), sample_length); // res / pi =>14005
//double[] f_1 = sin(.02, sample_length);
double[] f_2 = sin(fc(520), sample_length);
//double[] f_3 = sin(.25, sample_length);
//double[] f = add( add( add(f_0, f_1), f_2), f_3);
double[] f = add(f_0, f_2);
//print(f);
double[] d = cp(f,1000);
print_log(fft_log(d));
}
static double length = .2; // sec
static int res = 44000; // resoultion (pro sec)
static int sample_length = res; // resoultion
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Karlsen
Type : 24-dB (4-pole) lowpass
References : Posted by Best Regards,Ove Karlsen
Notes :
There's really not much voodoo going on in the filter itself, it's a simple as possible:
pole1 = (in * frequency) + (pole1 * (1 - frequency));
Most of you can probably understand that math, it's very similar to how an analog condenser works.
Although, I did have to do some JuJu to add resonance to it.
While studing the other filters, I found that the feedback phase is very important to how the overall
resonance level will be, and so I made a dynamic feedback path, and constant Q approximation by manipulation
of the feedback phase.
A bonus with this filter, is that you can "overdrive" it... Try high input levels..
Code :
// Karlsen 24dB Filter by Ove Karlsen / Synergy-7 in the year 2003.
// b_f = frequency 0..1
// b_q = resonance 0..50
// b_in = input
// to do bandpass, subtract poles from eachother, highpass subtract with input.
float b_inSH = b_in // before the while statement.
while (b_oversample < 2) { //2x oversampling (@44.1khz)
float prevfp;
prevfp = b_fp;
if (prevfp > 1) {prevfp = 1;} // Q-limiter
b_fp = (b_fp * 0.418) + ((b_q * pole4) * 0.582); // dynamic feedback
float intfp;
intfp = (b_fp * 0.36) + (prevfp * 0.64); // feedback phase
b_in = b_inSH - intfp; // inverted feedback
pole1 = (b_in * b_f) + (pole1 * (1 - b_f)); // pole 1
if (pole1 > 1) {pole1 = 1;} else if (pole1 < -1) {pole1 = -1;} // pole 1 clipping
pole2 = (pole1 * b_f) + (pole2 * (1 - b_f)); // pole 2
pole3 = (pole2 * b_f) + (pole3 * (1 - b_f)); // pole 3
pole4 = (pole3 * b_f) + (pole4 * (1 - b_f)); // pole 4
b_oversample++;
}
lowpassout = b_in;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Lo-Fi Crusher
Type : Quantizer / Decimator with smooth control
References : Posted by David Lowenfels
Notes :
Yet another bitcrusher algorithm. But this one has smooth parameter control.
Normfreq goes from 0 to 1.0; (freq/samplerate)
Input is assumed to be between 0 and 1.
Output gain is greater than unity when bits < 1.0;
Code :
function output = crusher( input, normfreq, bits );
step = 1/2^(bits);
phasor = 0;
last = 0;
for i = 1:length(input)
phasor = phasor + normfreq;
if (phasor >= 1.0)
phasor = phasor - 1.0;
last = step * floor( input(i)/step + 0.5 ); %quantize
end
output(i) = last; %sample and hold
end
end
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Look ahead limiting
References : Posted by Wilfried Welti
Notes :
use add_value with all values which enter the look-ahead area,
and remove_value with all value which leave this area. to get
the maximum value in the look-ahead area, use get_max_value.
in the very beginning initialize the table with zeroes.
If you always want to know the maximum amplitude in
your look-ahead area, the thing becomes a sorting
problem. very primitive approach using a look-up table
Code :
void lookup_add(unsigned section, unsigned size, unsigned value)
{
if (section==value)
lookup[section]++;
else
{
size >>= 1;
if (value>section)
{
lookup[section]++;
lookup_add(section+size,size,value);
}
else
lookup_add(section-size,size,value);
}
}
void lookup_remove(unsigned section, unsigned size, unsigned value)
{
if (section==value)
lookup[section]--;
else
{
size >>= 1;
if (value>section)
{
lookup[section]--;
lookup_remove(section+size,size,value);
}
else
lookup_remove(section-size,size,value);
}
}
unsigned lookup_getmax(unsigned section, unsigned size)
{
unsigned max = lookup[section] ? section : 0;
size >>= 1;
if (size)
if (max)
{
max = lookup_getmax((section+size),size);
if (!max) max=section;
}
else
max = lookup_getmax((section-size),size);
return max;
}
void add_value(unsigned value)
{
lookup_add(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value);
}
void remove_value(unsigned value)
{
lookup_remove(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value);
}
unsigned get_max_value()
{
return lookup_getmax(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1);
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Lowpass filter for parameter edge filtering
References : Olli Niemitalo
Linked file : filter001.gif
Notes :
use this filter to smooth sudden parameter changes
(see linkfile!)
Code :
/* - Three one-poles combined in parallel
* - Output stays within input limits
* - 18 dB/oct (approx) frequency response rolloff
* - Quite fast, 2x3 parallel multiplications/sample, no internal buffers
* - Time-scalable, allowing use with different samplerates
* - Impulse and edge responses have continuous differential
* - Requires high internal numerical precision
*/
{
/* Parameters */
// Number of samples from start of edge to halfway to new value
const double scale = 100;
// 0 < Smoothness < 1. High is better, but may cause precision problems
const double smoothness = 0.999;
/* Precalc variables */
double a = 1.0-(2.4/scale); // Could also be set directly
double b = smoothness; // -"-
double acoef = a;
double bcoef = a*b;
double ccoef = a*b*b;
double mastergain = 1.0 / (-1.0/(log(a)+2.0*log(b))+2.0/
(log(a)+log(b))-1.0/log(a));
double again = mastergain;
double bgain = mastergain * (log(a*b*b)*(log(a)-log(a*b)) /
((log(a*b*b)-log(a*b))*log(a*b))
- log(a)/log(a*b));
double cgain = mastergain * (-(log(a)-log(a*b)) /
(log(a*b*b)-log(a*b)));
/* Runtime variables */
long streamofs;
double areg = 0;
double breg = 0;
double creg = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++)
{
/* Update filters */
areg = acoef * areg + fromstream [streamofs];
breg = bcoef * breg + fromstream [streamofs];
creg = ccoef * creg + fromstream [streamofs];
/* Combine filters in parallel */
long temp = again * areg
+ bgain * breg
+ cgain * creg;
/* Check clipping */
if (temp > 32767)
{
temp = 32767;
}
else if (temp < -32768)
{
temp = -32768;
}
/* Store new value */
tostream [streamofs] = temp;
}
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
LP and HP filter
Type : biquad, tweaked butterworth
References : Posted by Patrice Tarrabia
Code :
r = rez amount, from sqrt(2) to ~ 0.1
f = cutoff frequency
(from ~0 Hz to SampleRate/2 - though many
synths seem to filter only up to SampleRate/4)
The filter algo:
out(n) = a1 * in + a2 * in(n-1) + a3 * in(n-2) - b1*out(n-1) - b2*out(n-2)
Lowpass:
c = 1.0 / tan(pi * f / sample_rate);
a1 = 1.0 / ( 1.0 + r * c + c * c);
a2 = 2* a1;
a3 = a1;
b1 = 2.0 * ( 1.0 - c*c) * a1;
b2 = ( 1.0 - r * c + c * c) * a1;
Hipass:
c = tan(pi * f / sample_rate);
a1 = 1.0 / ( 1.0 + r * c + c * c);
a2 = -2*a1;
a3 = a1;
b1 = 2.0 * ( c*c - 1.0) * a1;
b2 = ( 1.0 - r * c + c * c) * a1;
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
LPC analysis (autocorrelation + Levinson-Durbin recursion)
References : Posted by mail[AT]mutagene[DOT]net
Notes :
The
autocorrelation function implements a warped autocorrelation, so that
frequency resolution can be specified by the variable 'lambda'.
Levinson-Durbin recursion calculates autoregression coefficients a and
reflection coefficients (for lattice filter implementation) K. Comments
for Levinson-Durbin function implement matlab version of the same
function.
No optimizations.
Code :
//find the order-P autocorrelation array, R, for the sequence x of length L and warping of lambda
//wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0);
wAutocorrelate(float * x, unsigned int L, float * R, unsigned int P, float lambda)
{
double * dl = new double [L];
double * Rt = new double [L];
double r1,r2,r1t;
R[0]=0;
Rt[0]=0;
r1=0;
r2=0;
r1t=0;
for(unsigned int k=0; k<L;k++)
{
Rt[0]+=double(x[k])*double(x[k]);
dl[k]=r1-double(lambda)*double(x[k]-r2);
r1 = x[k];
r2 = dl[k];
}
for(unsigned int i=1; i<=P; i++)
{
Rt[i]=0;
r1=0;
r2=0;
for(unsigned int k=0; k<L;k++)
{
Rt[i]+=double(dl[k])*double(x[k]);
r1t = dl[k];
dl[k]=r1-double(lambda)*double(r1t-r2);
r1 = r1t;
r2 = dl[k];
}
}
for(i=0; i<=P; i++)
R[i]=float(Rt[i]);
delete[] dl;
delete[] Rt;
}
// Calculate the Levinson-Durbin recursion for the autocorrelation
sequence R of length P+1 and return the autocorrelation coefficients a
and reflection coefficients K
LevinsonRecursion(unsigned int P, float *R, float *A, float *K)
{
double Am1[62];
if(R[0]==0.0) {
for(unsigned int i=1; i<=P; i++)
{
K[i]=0.0;
A[i]=0.0;
}}
else {
double km,Em1,Em;
unsigned int k,s,m;
for (k=0;k<=P;k++){
A[0]=0;
Am1[0]=0; }
A[0]=1;
Am1[0]=1;
km=0;
Em1=R[0];
for (m=1;m<=P;m++) //m=2:N+1
{
double err=0.0f; //err = 0;
for (k=1;k<=m-1;k++) //for k=2:m-1
err += Am1[k]*R[m-k]; // err = err + am1(k)*R(m-k+1);
km = (R[m]-err)/Em1; //km=(R(m)-err)/Em1;
K[m-1] = -float(km);
A[m]=(float)km; //am(m)=km;
for (k=1;k<=m-1;k++) //for k=2:m-1
A[k]=float(Am1[k]-km*Am1[m-k]); // am(k)=am1(k)-km*am1(m-k+1);
Em=(1-km*km)*Em1; //Em=(1-km*km)*Em1;
for(s=0;s<=P;s++) //for s=1:N+1
Am1[s] = A[s]; // am1(s) = am(s)
Em1 = Em; //Em1 = Em;
}
}
return 0;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Magnitude and phase plot of arbitrary IIR function, up to 5th order
Type : magnitude and phase at any frequency
References : Posted by George Yohng
Notes :
Amplitude and phase calculation of IIR equation
run at sample rate "sampleRate" at frequency "F".
AMPLITUDE
-----------
cf_mag(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
PHASE
-----------
cf_phi(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
If you need a frequency diagram, draw a plot for
F=0...sampleRate/2
If you need amplitude in dB, use cf_lin2db(cf_mag(.......))
Set b0=-1 if you have such function:
y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
+ b1*y[n-1] + b2*y[n-2] + b3*y[n-3] + b4*y[n-4] + b5*y[n-5];
Set b0=1 if you have such function:
y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
- b1*y[n-1] - b2*y[n-2] - b3*y[n-3] - b4*y[n-4] - b5*y[n-5];
Do not try to reverse engineer these formulae - they don't give any sense
other than they are derived from transfer function, and they work. :)
Code :
/*
C file can be downloaded from
http://www.yohng.com/dsp/cfsmp.c
*/
#define C_PI 3.14159265358979323846264
double cf_mag(double f,double rate,
double a0,double a1,double a2,double a3,double a4,double a5,
double b0,double b1,double b2,double b3,double b4,double b5)
{
return
sqrt((a0*a0 + a1*a1 + a2*a2 + a3*a3 + a4*a4 + a5*a5 +
2*(a0*a1 + a1*a2 + a2*a3 + a3*a4 + a4*a5)*cos((2*f*C_PI)/rate) +
2*(a0*a2 + a1*a3 + a2*a4 + a3*a5)*cos((4*f*C_PI)/rate) +
2*a0*a3*cos((6*f*C_PI)/rate) + 2*a1*a4*cos((6*f*C_PI)/rate) +
2*a2*a5*cos((6*f*C_PI)/rate) + 2*a0*a4*cos((8*f*C_PI)/rate) +
2*a1*a5*cos((8*f*C_PI)/rate) + 2*a0*a5*cos((10*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) +
2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) +
2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) +
2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) +
2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate)));
}
double cf_phi(double f,double rate,
double a0,double a1,double a2,double a3,double a4,double a5,
double b0,double b1,double b2,double b3,double b4,double b5)
{
atan2((a0*b0 + a1*b1 + a2*b2 + a3*b3 + a4*b4 + a5*b5 +
(a0*b1 + a1*(b0 + b2) + a2*(b1 + b3) + a5*b4 + a3*(b2 + b4) +
a4*(b3 + b5))*cos((2*f*C_PI)/rate) +
((a0 + a4)*b2 + (a1 + a5)*b3 + a2*(b0 + b4) +
a3*(b1 + b5))*cos((4*f*C_PI)/rate) + a3*b0*cos((6*f*C_PI)/rate) +
a4*b1*cos((6*f*C_PI)/rate) + a5*b2*cos((6*f*C_PI)/rate) +
a0*b3*cos((6*f*C_PI)/rate) + a1*b4*cos((6*f*C_PI)/rate) +
a2*b5*cos((6*f*C_PI)/rate) + a4*b0*cos((8*f*C_PI)/rate) +
a5*b1*cos((8*f*C_PI)/rate) + a0*b4*cos((8*f*C_PI)/rate) +
a1*b5*cos((8*f*C_PI)/rate) +
(a5*b0 + a0*b5)*cos((10*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*((b0*b1 + b1*b2 + b3*(b2 + b4) + b4*b5)*cos((2*f*C_PI)/rate) +
(b2*(b0 + b4) + b3*(b1 + b5))*cos((4*f*C_PI)/rate) +
(b0*b3 + b1*b4 + b2*b5)*cos((6*f*C_PI)/rate) +
(b0*b4 + b1*b5)*cos((8*f*C_PI)/rate) +
b0*b5*cos((10*f*C_PI)/rate))),
((a1*b0 + a3*b0 + a5*b0 - a0*b1 + a2*b1 + a4*b1 - a1*b2 +
a3*b2 + a5*b2 - a0*b3 - a2*b3 + a4*b3 -
a1*b4 - a3*b4 + a5*b4 - a0*b5 - a2*b5 - a4*b5 +
2*(a3*b1 + a5*b1 - a0*b2 + a4*(b0 + b2) - a1*b3 + a5*b3 +
a2*(b0 - b4) - a0*b4 - a1*b5 - a3*b5)*cos((2*f*C_PI)/rate) +
2*(a3*b0 + a4*b1 + a5*(b0 + b2) - a0*b3 - a1*b4 - a0*b5 - a2*b5)*
cos((4*f*C_PI)/rate) + 2*a4*b0*cos((6*f*C_PI)/rate) +
2*a5*b1*cos((6*f*C_PI)/rate) - 2*a0*b4*cos((6*f*C_PI)/rate) -
2*a1*b5*cos((6*f*C_PI)/rate) + 2*a5*b0*cos((8*f*C_PI)/rate) -
2*a0*b5*cos((8*f*C_PI)/rate))*sin((2*f*C_PI)/rate))/
(b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 +
2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) +
2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) +
2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) +
2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) +
2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate)));
}
double cf_lin2db(double lin)
{
if (lin<9e-51) return -1000; /* prevent invalid operation */
return 20*log10(lin);
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Measuring interpollation noise
References : Posted by Jon Watte
Notes :
You can easily estimate the error by evaluating the actual function and
evaluating your interpolator at each of the mid-points between your
samples. The absolute difference between these values, over the absolute
value of the "correct" value, is your relative error. log10 of your relative
error times 20 is an estimate of your quantization noise in dB. Example:
You have a table for every 0.5 "index units". The value at index unit 72.0
is 0.995 and the value at index unit 72.5 is 0.999. The interpolated value
at index 72.25 is 0.997. Suppose the actual function value at that point was
0.998; you would have an error of 0.001 which is a relative error of 0.001002004..
log10(error) is about -2.99913, which times 20 is about -59.98. Thus, that's
your quantization noise at that position in the table. Repeat for each pair of
samples in the table.
Note: I said "quantization noise" not "aliasing noise". The aliasing noise will,
as far as I know, only happen when you start up-sampling without band-limiting
and get frequency aliasing (wrap-around), and thus is mostly independent of
what specific interpolation mechanism you're using.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
MIDI note/frequency conversion
Type : -
References : Posted by tobybear[AT]web[DOT]de
Notes :
I
get often asked about simple things like MIDI note/frequency
conversion, so I thought I could as well post some source code about
this.
The following is Pascal/Delphi syntax, but it shouldn't be a problem to convert it to almost any language in no time.
Uses for this code are mainly for initializing oscillators to the right
frequency based upon a given MIDI note, but you might also check what
MIDI note is closest to a given frequency for pitch detection etc.
In realtime applications it might be a good idea to get rid of the
power and log2 calculations and generate a lookup table on
initialization.
A full Pascal/Delphi unit with these functions (including lookup table
generation) and a simple demo application can be downloaded here:
http://tobybear.phreque.com/dsp_conv.zip
If you have any comments/suggestions, please send them to: tobybear@web.de
Code :
// MIDI NOTE/FREQUENCY CONVERSIONS
const notes:array[0..11] of string= ('C ','C#','D ','D#','E ','F ','F#','G ','G#','A ','A#','B ');
const base_a4=440; // set A4=440Hz
// converts from MIDI note number to frequency
// example: NoteToFrequency(12)=32.703
function NoteToFrequency(n:integer):double;
begin
if (n>=0)and(n<=119) then
result:=base_a4*power(2,(n-57)/12)
else result:=-1;
end;
// converts from MIDI note number to string
// example: NoteToName(12)='C 1'
function NoteToName(n:integer):string;
begin
if (n>=0)and(n<=119) then
result:=notes[n mod 12]+inttostr(n div 12)
else result:='---';
end;
// converts from frequency to closest MIDI note
// example: FrequencyToNote(443)=57 (A 4)
function FrequencyToNote(f:double):integer;
begin
result:=round(12*log2(f/base_a4))+57;
end;
// converts from string to MIDI note
// example: NameToNote('A4')=57
function NameToNote(s:string):integer;
var c,i:integer;
begin
if length(s)=2 then s:=s[1]+' '+s[2];
if length(s)<>3 then begin result:=-1;exit end;
s:=uppercase(s);
c:=-1;
for i:=0 to 11 do
if notes[i]=copy(s,1,2) then
begin
c:=i;
break
end;
try
i:=strtoint(s[3]);
result:=i*12+c;
except
result:=-1;
end;
if c<0 then result:=-1;
end;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Millimeter to DB (faders...)
References : Posted by James McCartney
Notes :
These two functions reproduce a traditional professional
mixer fader taper.
MMtoDB converts millimeters of fader travel from the
bottom of the fader for a 100 millimeter fader into
decibels. DBtoMM is the inverse.
The taper is as follows from the top:
The top of the fader is +10 dB
100 mm to 52 mm : -5 dB per 12 mm
52 mm to 16 mm : -10 dB per 12 mm
16 mm to 4 mm : -20 dB per 12 mm
4 mm to 0 mm : fade to zero. (in these functions I go to -200dB
which is effectively zero for up to 32 bit audio.)
Code :
float MMtoDB(float mm)
{
float db;
mm = 100. - mm;
if (mm <= 0.) {
db = 10.;
} else if (mm < 48.) {
db = 10. - 5./12. * mm;
} else if (mm < 84.) {
db = -10. - 10./12. * (mm - 48.);
} else if (mm < 96.) {
db = -40. - 20./12. * (mm - 84.);
} else if (mm < 100.) {
db = -60. - 35. * (mm - 96.);
} else db = -200.;
return db;
}
float DBtoMM(float db)
{
float mm;
if (db >= 10.) {
mm = 0.;
} else if (db > -10.) {
mm = -12./5. * (db - 10.);
} else if (db > -40.) {
mm = 48. - 12./10. * (db + 10.);
} else if (db > -60.) {
mm = 84. - 12./20. * (db + 40.);
} else if (db > -200.) {
mm = 96. - 1./35. * (db + 60.);
} else mm = 100.;
mm = 100. - mm;
return mm;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper.
Notes :
Digital approximation of Moog VCF. Fairly easy to calculate coefficients, fairly easy to process algorithm, good sound.
Code :
//Init
cutoff = cutoff freq in Hz
fs = sampling frequency //(e.g. 44100Hz)
res = resonance [0 - 1] //(minimum - maximum)
f = 2 * cutoff / fs; //[0 - 1]
k = 3.6*f - 1.6*f*f -1; //(Empirical tunning)
p = (k+1)*0.5;
scale = e^((1-p)*1.386249;
r = res*scale;
y4 = output;
y1=y2=y3=y4=oldx=oldy1=oldy2=oldy3=0;
//Loop
//--Inverted feed back for corner peaking
x = input - r*y4;
//Four cascaded onepole filters (bilinear transform)
y1=x*p + oldx*p - k*y1;
y2=y1*p+oldy1*p - k*y2;
y3=y2*p+oldy2*p - k*y3;
y4=y3*p+oldy3*p - k*y4;
//Clipper band limited sigmoid
y4 = y4 - (y4^3)/6;
oldx = x;
oldy1 = y1;
oldy2 = y2;
oldy3 = y3;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF, variation 1
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper., Paul Kellett version
Notes :
The
second "q =" line previously used exp() - I'm not sure if what I've
done is any faster, but this line needs playing with anyway as it
controls which frequencies will self-oscillate. I
think it could be tweaked to sound better than it currently does.
Highpass / Bandpass :
They are only 6dB/oct, but still seem musically useful - the 'fruity' sound of the 24dB/oct lowpass is retained.
Code :
// Moog 24 dB/oct resonant lowpass VCF
// References: CSound source code, Stilson/Smith CCRMA paper.
// Modified by paul.kellett@maxim.abel.co.uk July 2000
float f, p, q; //filter coefficients
float b0, b1, b2, b3, b4; //filter buffers (beware denormals!)
float t1, t2; //temporary buffers
// Set coefficients given frequency & resonance [0.0...1.0]
q = 1.0f - frequency;
p = frequency + 0.8f * frequency * q;
f = p + p - 1.0f;
q = resonance * (1.0f + 0.5f * q * (1.0f - q + 5.6f * q * q));
// Filter (in [-1.0...+1.0])
in -= q * b4; //feedback
t1 = b1; b1 = (in + b0) * p - b1 * f;
t2 = b2; b2 = (b1 + t1) * p - b2 * f;
t1 = b3; b3 = (b2 + t2) * p - b3 * f;
b4 = (b3 + t1) * p - b4 * f;
b4 = b4 - b4 * b4 * b4 * 0.166667f; //clipping
b0 = in;
// Lowpass output: b4
// Highpass output: in - b4;
// Bandpass output: 3.0f * (b3 - b4);
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF, variation 2
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper., Timo Tossavainen (?) version
Notes :
in[x] and out[x] are member variables, init to 0.0 the controls:
fc = cutoff, nearly linear [0,1] -> [0, fs/2]
res = resonance [0, 4] -> [no resonance, self-oscillation]
Code :
Tdouble MoogVCF::run(double input, double fc, double res)
{
double f = fc * 1.16;
double fb = res * (1.0 - 0.15 * f * f);
input -= out4 * fb;
input *= 0.35013 * (f*f)*(f*f);
out1 = input + 0.3 * in1 + (1 - f) * out1; // Pole 1
in1 = input;
out2 = out1 + 0.3 * in2 + (1 - f) * out2; // Pole 2
in2 = out1;
out3 = out2 + 0.3 * in3 + (1 - f) * out3; // Pole 3
in3 = out2;
out4 = out3 + 0.3 * in4 + (1 - f) * out4; // Pole 4
in4 = out3;
return out4;
}
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Most simple and smooth feedback delay
Type : Feedback delay
References : Posted by antiprosynthesis[AT]hotmail[DOT]com
Notes :
fDlyTime = delay time parameter (0-1)
i = input index
j = delay index
Code :
if( i >= SampleRate )
i = 0;
j = i - (fDlyTime * SampleRate);
if( j < 0 )
j += SampleRate;
Output = DlyBuffer[ i++ ] = Input + (DlyBuffer[ j ] * fFeedback);
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Most simple static delay
Type : Static delay
References : Posted by antiprosynthesis[AT]hotmail[DOT]com
Notes :
This
is the most simple static delay (just delays the input sound an amount
of samples). Very useful for newbies also probably very easy to change
in a feedback delay (for comb filters for example).
Note: fDlyTime is the delay time parameter (0 to 1)
i = input index
j = output index
Code :
if( i >= SampleRate )
i = 0;
DlyBuffer[ i ] = Input;
j = i - (fDlyTime * SampleRate);
i++;
if( j < 0 )
j = SampleRate + j;
Output = DlyBuffer[ j ];
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Noise Shaping Class
Type : Dithering with 9th order noise shaping
References : Posted by cshei[AT]indiana.edu
Linked file : NS9dither16.h
Notes :
This
is an implemetation of a 9th order noise shaping & dithering class,
that runs quite fast (it has one function that uses Intel x86 assembly,
but you can replace it with a different rounding function if you are
running on a non-Intel platform). _aligned_malloc and _aligned_free
require the MSVC++ Processor Pack, available from www.microsoft.com.
You can replace them with "new" and "delete," but allocating aligned
memory seems to make it run faster. Also, you can replace ZeroMemory
with a memset that sets the memory to 0 if you aren't using Win32.
Input should be floats from -32768 to 32767 (processS will clip at
these points for you, but clipping is bad when you are trying to
convert floats to shorts). Note to reviewer - it would probably be
better if you put the code in a file such as NSDither.h and have a link
to it - it's rather long.
(see linked file)
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Notch filter
Type : 2 poles 2 zeros IIR
References : Posted by Olli Niemitalo
Notes :
Creates
a muted spot in the spectrum with adjustable steepness. A complex
conjugate pair of zeros on the z- plane unit circle and neutralizing
poles approaching at the same angles from inside the unit circle.
Code :
Parameters:
0 =< freq =< samplerate/2
0 =< q < 1 (The higher, the narrower)
AlgoAlgo=double pi = 3.141592654;
double sqrt2 = sqrt(2.0);
double freq = 2050; // Change! (zero & pole angle)
double q = 0.4; // Change! (pole magnitude)
double z1x = cos(2*pi*freq/samplerate);
double a0a2 = (1-q)*(1-q)/(2*(fabs(z1x)+1)) + q;
double a1 = -2*z1x*a0a2;
double b1 = -2*z1x*q;
double b2 = q*q;
double reg0, reg1, reg2;
unsigned int streamofs;
reg1 = 0;
reg2 = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++)
{
reg0 = a0a2 * ((double)fromstream[streamofs]
+ fromstream[streamofs+2])
+ a1 * fromstream[streamofs+1]
- b1 * reg1
- b2 * reg2;
reg2 = reg1;
reg1 = reg0;
int temp = reg0;
/* Check clipping */
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
/* Store new value */
tostream[streamofs] = temp;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One pole LP and HP
References : Posted by Bram
Code :
LP:
recursion: tmp = (1-p)*in + p*tmp with output = tmp
coefficient: p = (2-cos(x)) - sqrt((2-cos(x))^2 - 1) with x = 2*pi*cutoff/samplerate
coeficient approximation: p = (1 - 2*cutoff/samplerate)^2
HP:
recursion: tmp = (p-1)*in - p*tmp with output = tmp
coefficient: p = (2+cos(x)) - sqrt((2+cos(x))^2 - 1) with x = 2*pi*cutoff/samplerate
coeficient approximation: p = (2*cutoff/samplerate)^2
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One pole, one zero LP/HP
References : Posted by mistert[AT]inwind[DOT]it
Code :
void SetLPF(float fCut, float fSampling)
{
float w = 2.0 * fSampling;
float Norm;
fCut *= 2.0F * PI;
Norm = 1.0 / (fCut + w);
b1 = (w - fCut) * Norm;
a0 = a1 = fCut * Norm;
}
void SetHPF(float fCut, float fSampling)
{
float w = 2.0 * fSampling;
float Norm;
fCut *= 2.0F * PI;
Norm = 1.0 / (fCut + w);
a0 = w * Norm;
a1 = -a0;
b1 = (w - fCut) * Norm;
}
Where
out[n] = in[n]*a0 + in[n-1]*a1 + out[n-1]*b1;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One zero, LP/HP
References : Posted by Bram
Notes :
LP is only 'valid' for cutoffs > samplerate/4
HP is only 'valid' for cutoffs < samplerate/4
Code :
theta = cutoff*2*pi / samplerate
LP:
H(z) = (1+p*z^(-1)) / (1+p)
out[i] = 1/(1+p) * in[i] + p/(1+p) * in[i-1];
p = (1-2*cos(theta)) - sqrt((1-2*cos(theta))^2 - 1)
Pi/2 < theta < Pi
HP:
H(z) = (1-p*z^(-1)) / (1+p)
out[i] = 1/(1+p) * in[i] - p/(1+p) * in[i-1];
p = (1+2*cos(theta)) - sqrt((1+2*cos(theta))^2 - 1)
0 < theta < Pi/2
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Parallel combs delay calculation
References : Posted by Juhana Sadeharju ( kouhia[AT]nic[DOT]funet[DOT]fi )
Notes :
This formula can be found from a patent related to parallel combs
structure. The formula places the first echoes coming out of parallel
combs to uniformly distributed sequence. If T_ ,...,T_n are the delay
lines in increasing order, the formula can be derived by setting
T_(k-1)/T_k = Constant and T_n/(2*T_1) = Constant, where 2*T_1 is the
echo coming just after the echo T_n.
I figured this out myself as
it is not told in the patent. The formula is not the best which one can
come up. I use a search method to find echo sequences which are uniform
enough for long enough time. The formula is uniform for a short time
only.
The formula doesn't work good for series allpass and FDN structures,
for which a similar formula can be derived with the same idea. The
search method works for these structures as well.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Peak/Notch filter
Type : peak/notch
References : Posted by tobybear[AT]web[DOT]de
Notes :
// Peak/Notch filter
// I don't know anymore where this came from, just found it on
// my hard drive :-)
// Seems to be a peak/notch filter with adjustable slope
// steepness, though slope gets rather wide the lower the
// frequency is.
// "cut" and "steep" range is from 0..1
// Try to feed it with white noise, then the peak output does
// rather well eliminate all other frequencies except the given
// frequency in higher frequency ranges.
Code :
var f,r:single;
outp,outp1,outp2:single; // init these with 0!
const p4=1.0e-24; // Pentium 4 denormal problem elimination
function PeakNotch(inp,cut,steep:single;ftype:integer):single;
begin
r:=steep*0.99609375;
f:=cos(pi*cut);
a0:=(1-r)*sqrt(r*(r-4*(f*f)+2)+1);
b1:=2*f*r;
b2:=-(r*r);
outp:=a0*inp+b1*outp1+b2*outp2+p4;
outp2:=outp1;
outp1:=outp;
if ftype=0 then
result:=outp //peak
else
result:=inp-outp; //notch
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Phase equalization
Type : Allpass
References : Posted by Uli the Grasso
Notes :
The
idea is simple: One can equalize the phase response of a system, for
example of a loudspeaker, by approximating its phase response by an FIR
filter and then turn around the coefficients of the filter. At
http://grassomusic.de/english/phaseeq.htm you find more info and an
Octave script.
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Phase modulation Vs. Frequency modulation II
References : Posted by James McCartney
Notes :
The
difference between FM & PM in a digital oscillator is that FM is
added to the frequency before the phase integration, while PM is added
to the phase after the phase integration. Phase integration is when the
old phase for the oscillator is added to the current frequency (in
radians per sample) to get the new phase for the oscillator. The
equivalent PM modulator to obtain the same waveform as FM is the
integral of the FM modulator. Since the integral of sine waves are
inverted cosine waves this is no problem. In modulators with multiple
partials, the equivalent PM modulator will have different relative
partial amplitudes. For example, the integral of a square wave is a
triangle wave; they have the same harmonic content, but the relative
partial amplitudes are different. These differences make no difference
since we are not trying to exactly recreate FM, but real (or nonreal)
instruments.
The reason PM is better is because in PM and FM there can be non-zero
energy produced at 0 Hz, which in FM will produce a shift in pitch if
the FM wave is used again as a modulator, however in PM the DC
component will only produce a phase shift. Another reason PM is better
is that the modulation index (which determines the number of sidebands
produced and which in normal FM is calculated as the modulator
amplitude divided by frequency of modulator) is not dependant on the
frequency of the modulator, it is always equal to the amplitude of the
modulator in radians. The benefit of solving the DC frequency shift
problem, is that cascaded carrier-modulator pairs and feedback
modulation are possible. The simpler calculation of modulation index
makes it easier to have voices keep the same harmonic structure
throughout all pitches.
The basic mathematics of phase modulation are available in any text on electronic communication theory.
Below is some C code for a digital oscillator that implements FM,PM,and
AM. It illustrates the difference in implementation of FM & PM. It
is only meant as an example, and not as an efficient implementation.
Code :
/* Example implementation of digital oscillator with FM, PM, & AM */
#define PI 3.14159265358979
#define RADIANS_TO_INDEX (512.0 / (2.0 * PI))
typedef struct{ /* oscillator data */
double freq; /* oscillator frequency in radians per sample */
double phase; /* accumulated oscillator phase in radians */
double wavetable[512]; /* waveform lookup table */
} OscilRec;
/* oscil - compute 1 sample of oscillator output whose freq. phase and
* wavetable are in the OscilRec structure pointed to by orec.
*/
double oscil(orec, fm, pm, am)
OscilRec *orec; /* pointer to the oscil's data */
double fm; /* frequency modulation input in radians per sample */
double pm; /* phase modulation input in radians */
double am; /* amplitude modulation input in any units you want */
{
long tableindex; /* index into wavetable */
double instantaneous_freq; /* oscillator freq + freq modulation */
double instantaneous_phase; /* oscillator phase + phase modulation */
double output; /* oscillator output */
instantaneous_freq = orec->freq + fm; /* get instantaneous freq */
orec->phase += instantaneous_freq; /* accumulate phase */
instantaneous_phase = orec->phase + pm; /* get instantaneous phase */
/* convert to lookup table index */
tableindex = RADIANS_TO_INDEX * instantaneous_phase;
tableindex &= 511; /* make it mod 512 === eliminate multiples of 2*k*PI */
output = orec->wavetable[tableindex] * am; /* lookup and mult by am input */
return (output); /* return oscillator output */
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
pow(x,4) approximation
References : Posted by Stefan Stenzel
Notes :
Very hacked, but it gives a rough estimate of x**4 by modifying exponent and mantissa.
Code :
float p4fast(float in)
{
long *lp,l;
lp=(long *)(&in);
l=*lp;
l-=0x3F800000l; /* un-bias */
l<<=2; /* **4 */
l+=0x3F800000l; /* bias */
*lp=l;
/* compiler will read this from memory since & operator had been used */
return in;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Prewarping
Type : explanation
References : Posted by robert bristow-johnson (better known as "rbj" )
Notes :
prewarping is simply recognizing the warping that the BLT introduces.
to determine frequency response, we evaluate the digital H(z) at
z=exp(j*w*T) and we evaluate the analog Ha(s) at s=j*W . the following
will confirm the jw to unit circle mapping and will show exactly what the
mapping is (this is the same stuff in the textbooks):
the BLT says: s = (2/T) * (z-1)/(z+1)
substituting: s = j*W = (2/T) * (exp(j*w*T) - 1) / (exp(j*w*T) + 1)
j*W = (2/T) * (exp(j*w*T/2) - exp(-j*w*T/2)) / (exp(j*w*T/2) + exp(-j*w*T/2))
= (2/T) * (j*2*sin(w*T/2)) / (2*cos(w*T/2))
= j * (2/T) * tan(w*T/2)
or
analog W = (2/T) * tan(w*T/2)
so when the real input frequency is w, the digital filter will behave with
the same amplitude gain and phase shift as the analog filter will have at a
hypothetical frequency of W. as w*T approaches pi (Nyquist) the digital
filter behaves as the analog filter does as W -> inf. for each degree of
freedom that you have in your design equations, you can adjust the analog
design frequency to be just right so that when the deterministic BLT
warping does its thing, the resultant warped frequency comes out just
right. for a simple LPF, you have only one degree of freedom, the cutoff
frequency. you can precompensate it so that the true cutoff comes out
right but that is it, above the cutoff, you will see that the LPF dives
down to -inf dB faster than an equivalent analog at the same frequencies.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Pseudo-Random generator
Type : Linear Congruential, 32bit
References : Hal Chamberlain, "Musical Applications of Microprocessors" (Posted by Phil Burk)
Notes :
This can be used to generate random numeric sequences or to synthesise a white noise audio signal.
If you only use some of the bits, use the most significant bits by shifting right.
Do not just mask off the low bits.
Code :
/* Calculate pseudo-random 32 bit number based on linear congruential method. */
unsigned long GenerateRandomNumber( void )
{
/* Change this for different random sequences. */
static unsigned long randSeed = 22222;
randSeed = (randSeed * 196314165) + 907633515;
return randSeed;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Pulsewidth modulation
Type : waveform generation
References : Steffan Diedrichsen
Notes :
Take
an upramping sawtooth and its inverse, a downramping sawtooth. Adding
these two waves with a well defined delay between 0 and period (1/f)
results in a square wave with a duty cycle ranging from 0 to 100%.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
QFT and DQFT (double precision) classes
References : Posted by Joshua Scholar
Linked file : qft.tar_1.gz
Notes :
Since it's a Visual C++ project (though it has relatively portable C++) I
guess the main audience are PC users. As such I'm including a zip file.
Some PC users wouldn't know what to do with a tgz file.
The QFT and DQFT (double precision) classes supply the following functions:
1. Real valued FFT and inverse FFT functions. Note that separate arraysare used
for real and imaginary component of the resulting spectrum.
2. Decomposition of a spectrum into a separate spectrum of the evensamples
and a spectrum of the odd samples. This can be useful for buildingfilter banks.
3. Reconstituting a spectrum from separate spectrums of the even samples
and odd samples. This can be useful for building filter banks.
4. A discrete Sin transform (a QFT decomposes an FFT into a DST and DCT).
5. A discrete Cos transfrom.
6. Since a QFT does it's last stage calculating from the outside in thelast part
can be left unpacked and only calculated as needed in the case wherethe entire
spectrum isn't needed (I used this for calculating correlations andconvolutions
where I only needed half of the results).
ReverseNoUnpack()
UnpackStep()
and NegUnpackStep()
implement this functionality
NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseNoUnpack() does not.
7. Also if you only want the first half of the results you can call ReverseHalf()
NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseHalf() does not.
8. QFT is less numerically stable than regular FFTs. With singleprecision calculations,
a block length of 2^15 brings the accuracy down to being barelyaccurate enough.
At that size, single precision calculations tested sound files wouldoccasionally have
a sample off by 2, and a couple off by 1 per block. Full volume whitenoise would generate
a few samples off by as much as 6 per block at the end, beginning and middle.
No matter what the inputs the errors are always at the same positions in the block.
There some sort of cancelation that gets more delicate as the block size gets bigger.
For the sake of doing convolutions and the like where the forward transform is
done only once for one of the inputs, I created a AccurateForward() function.
It uses a regular FFT algorithm for blocks larger than 2^12, and decomposes into even and
odd FFTs recursively.
In any case you can always use the double precision routines to get more accuracy.
DQFT even has routines that take floats as inputs and return double precision
spectrum outputs.
As for portability:
1. The files qft.cpp and dqft.cpp start with defines:
#define _USE_ASM
If you comment those define out, then what's left is C++ with no assembly language.
2. There is unnecessary windows specific code in "criticalSection.h"
I used a critical section because objects are not reentrant (each object has
permanent scratch pad memory), but obviously critical sections are operating
system specific. In any case that code can easily be taken out.
If you look at my code and see that there's an a test built in the examples
that makes sure that the results are in the ballpark of being right. It
wasn't that I expected the answers to be far off, it was that I uncommenting
the "no assembly language" versions of some routines and I wanted to make
sure that they weren't broken.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
quick and dirty sine generator
Type : sine generator
References : Posted by couriervst[AT]hotmail[DOT]com
Notes :
this is part of my library, although I've seen a lot of sine generators, I've never seen the simplest one, so I try to do it,
tell me something, I've try it and work so tell me something about it
Code :
PSPsample PSPsin1::doOsc(int numCh)
{
double x=0;
double t=0;
if(m_time[numCh]>m_sampleRate) //re-init cycle
m_time[numCh]=0;
if(m_time[numCh]>0)
{
t =(double)(((double)m_time[numCh])/(double)m_sampleRate);
x=(m_2PI *(double)(t)*m_freq);
}
else
x=0;
PSPsample r=(PSPsample) sin(x+m_phase)*m_amp;
m_time[numCh]++;
return r;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Real basic DSP with Matlab (+ GUI) ...
Type : Like effects racks, made with Matlab !
References : Posted by guillaume[DOT]carniato[AT]meletu[DOT]univ-valenciennes[DOT]fr
Linked file : http://www.xenggeng.fr.st/ici/guitou/Matlab Music.zip
Notes :
You need Matlab v6.0 or more to run this stuff...
Code :
take a look at http://www.xenggeng.fr.st/ici/guitou/Matlab Music.zip
I'm now working on a Matlab - sequencer, which will certainly use
'Matlab Music'. I'm interested in integrating WaveWarp in this project;
it's a toolbox that allow you to make real time DSP with Matlab.
If you're ok to improve this version (add effects, improve effects
quality,anything else...) let's go ! Email me if you're interested in
developing this beginner work...
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
real value vs display value
Type : Macro
References : Posted by emil[AT]arpanet[DOT]no
Notes :
REALVAL converts the vst param at given ranges to a display value.
VSTVAL does the opposite.
a = start
b = end
Code :
#define REALVAL(a, b, vstval) (a + (vstval)*(b-a))
#define VSTVAL(a, b, realval) ((realval-a)/(b-a))
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Really fast x86 floating point sin/cos
References : Posted by rerdavies[AT]msn[DOT]com
Linked file : sincos.zip
Notes :
Frightful code from the Intel Performance optimization front. Not for the squeamish.
The following code calculates sin and cos of a floating point value on
x86 platforms to 20 bits precision with 2 multiplies and two adds. The
basic principle is to use sin(x+y) and cos(x+y) identities to generate
the result from lookup tables. Each lookup table takes care of 10 bits
of precision in the input. The same principle can be used to generate
sin/cos to full (! Really. Full!) 24-bit float precision using two
8-bit tables, and one 10 bit table (to provide guard bits), for a net
speed gain of about 4x over fsin/fcos, and 8x if you want both sin and
cos. Note that microsoft compilers have trouble keeping doubles aligned
properly on the stack (they must be 8-byte aligned in order not to
incur a massive alignment penalty). As a result, this class should NOT
be allocated on the stack. Add it as a member variable to any class
that uses it.
e.g.
class CSomeClass {
CQuickTrig m_QuickTrig;
...
mQuickTrig.QuickSinCos(dAngle,fSin,fCos);
...
}
Code :
(see attached file)
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Reasonably accurate/fastish tanh approximation
References : Posted by Fuzzpilz
Notes :
Fairly obvious, but maybe not obvious enough, since I've seen calls to tanh() in code snippets here.
It's this, basically:
tanh(x) = sinh(x)/cosh(x)
= (exp(x) - exp(-x))/(exp(x) + exp(-x))
= (exp(2x) - 1)/(exp(2x) + 1)
Combine this with a somewhat less accurate approximation for exp than
usual (I use a third-order Taylor approximation below), and you're set.
Useful for waveshaping.
Notes on the exp approximation:
It only works properly for input values above x, but since tanh is odd, that isn't a problem.
exp(x) = 1 + x + x^2/(2!) + x^3/(3!) + ...
Breaking the Taylor series off after the third term, I get
1 + x + x^2/2 + x^3/6.
I can save some multiplications by using
6 + x * (6 + x * (3 + x))
instead; a division by 6 becomes necessary, but is lumped into the additions in the tanh part:
(a/6 - 1)/(a/6 + 1) = (a - 6)/(a + 6).
Accuracy:
I haven't looked at this in very great detail, but it's always <=
the real tanh (>= for x<0), and the greatest deviation occurs at
about +/- 1.46, where the result is ca. .95 times the correct value.
This is still faster than tanh if you use a better approximation for the exponential, even if you simply call exp.
There are probably additional ways of improving parts of this, and
naturally if you're going to use it you'll want to figure out whether
your particular application offers additional ways of simplifying it,
but it's a good start.
Code :
/*
single precision absolute value, a lot faster than fabsf() (if you use
MSVC++ 6 Standard - others' implementations might be less slow) */
float sabs(float a)
{
int b=(*((int *)(&a)))&0x7FFFFFFF;
return *((float *)(&b));
}
/* approximates tanh(x/2) rather than tanh(x) - depending on how you're
using this, fixing that could well be wasting a multiplication (though
that isn't much, and it could be done with an integer addition in sabs
instead) */
float tanh2(float x)
{
float a=sabs(x);
a=6+a*(6+a*(3+a));
return ((x<0)?-1:1)*(a-6)/(a+6); /* instead of using <, you can
also check directly whether the sign bit is set ((*((int
*)(&x)))&0x80000000), but it's not really worth it */
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
resampling
Type : linear interpolated aliased resampling of a wave file
References : Posted by mail[AT]mroc[DOT]de
Notes :
som resampling stuff. the code is heavily used in MSynth, but do not lough about ;-)
perhaps, prefiltering would reduce aliasing.
Code :
signed short* pSample = ...;
unsigned int sampleLength = ...;
// stretch sample to length of one bar...
float playPosDelta = sampleLength / ( ( 240.0f / bpm ) * samplingRate );
// requires for position calculation...
float playpos1 = 0.0f;
unsigned int iter = 0;
// required for interpolation...
unsigned int i1, i2;
float* pDest = ....;
float* pDestStop = pDest + len;
for( float *pt=pDest;pt<pDestStop;++pt )
{
// linear interpolation...
i1 = (unsigned int)playpos;
i2 = i1 + 1;
(*pt) = ((pSample[i2]-pSample[i1]) * (playpos - i1) + pSample[i1]);
// position calculation preventing float sumation error...
playpos1 = (++iter) * playposIncrement;
}
...
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Resonant filter
References : Posted by Paul Kellett
Notes :
This filter consists of two first order low-pass filters in
series, with some of the difference between the two filter
outputs fed back to give a resonant peak.
You can use more filter stages for a steeper cutoff but the
stability criteria get more complicated if the extra stages
are within the feedback loop.
Code :
//set feedback amount given f and q between 0 and 1
fb = q + q/(1.0 - f);
//for each sample...
buf0 = buf0 + f * (in - buf0 + fb * (buf0 - buf1));
buf1 = buf1 + f * (buf0 - buf1);
out = buf1;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Resonant IIR lowpass (12dB/oct)
Type : Resonant IIR lowpass (12dB/oct)
References : Posted by Olli Niemitalo
Notes :
Hard to calculate coefficients, easy to process algorithm
Code :
resofreq = pole frequency
amp = magnitude at pole frequency (approx)
double pi = 3.141592654;
/* Parameters. Change these! */
double resofreq = 5000;
double amp = 1.0;
DOUBLEWORD streamofs;
double w = 2.0*pi*resofreq/samplerate; // Pole angle
double q = 1.0-w/(2.0*(amp+0.5/(1.0+w))+w-2.0); // Pole magnitude
double r = q*q;
double c = r+1.0-2.0*cos(w)*q;
double vibrapos = 0;
double vibraspeed = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++) {
/* Accelerate vibra by signal-vibra, multiplied by lowpasscutoff */
vibraspeed += (fromstream[streamofs] - vibrapos) * c;
/* Add velocity to vibra's position */
vibrapos += vibraspeed;
/* Attenuate/amplify vibra's velocity by resonance */
vibraspeed *= r;
/* Check clipping */
temp = vibrapos;
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
/* Store new value */
tostream[streamofs] = temp;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Reverb Filter Generator
Type : FIR
References : Posted by Stephen McGovern
Notes :
This
is a MATLAB function that makes a rough calculation of a room's impulse
response. The output can then be convolved with an audio clip to
produce good and realistic sounding reverb. I have written a paper
discussing the theory used by this algorithm. It is available at
http://stevem.us/rir.html.
NOTES:
1) Large values of N will use large amounts of memory.
2) The output is normalized to the largest value of the
output.
Code :
function [h]=rir(fs, mic, n, r, rm, src);
%RIR Room Impulse Response.
% [h] = RIR(FS, MIC, N, R, RM, SRC) calculates the impulse % response
% of a room.
%
% FS = sample rate.
% MIC = row vector, MIC=[X Y Z], giving the x,y,z % % coordinates of
% the microphone.
% N = The program will account for (2*N+1)^3 virtual % sources
% R = reflection coefficient for the walls, 0<R<1.
% RM = row vector, RM=[X Y Z], giving the dimensions % of the room.
% SRC = row vector, SRC=[X Y Z], giving the x,y,z % coordinates of
% the sound source.
%
% NOTES:
%
% 1) To implement this filter, you will need to do a fast
% convolution. The program FCONV.m will do this. It is % available
% at: http://stevem.us/code/fconv.m
% 2) All distances are in meters.
% 3) If this is your first time running this program, set % N equal
% to 10 or less, and R=0.9.
% 4) I've written an article discussing the theory behind % this
% algorithm. It can be found at
% http://stevem.us/rir.html.
%
%
%Version 1.0
%Coded by: Stephen G. McGovern, 2003.
%The comments below refer to equations in my paper.
nn=[-n:1:n]; % Index for the
% sequence
rms= nn+0.5-0.5*(-1).^nn; % Part of equations
% 2,3,& 4
srcs=(-1).^(nn); % part of equations
% 2,3,& 4
xi=[srcs*src(1)+rms*rm(1)-mic(1)]; % Equation 2
yj=[srcs*src(2)+rms*rm(2)-mic(2)]; % Equation 3
zk=[srcs*src(3)+rms*rm(3)-mic(3)]; % Equation 4
[i,j,k]=meshgrid(xi,yj,zk); % convert vectors to
% 3D matrices
d=sqrt(i.^2+j.^2+k.^2); % Equation 5
time=d./343; % Similar to equation
% 6
time=round(time*fs); % Quantized delay time
b=1./(4*pi*((d).^(2))); % Equation 8
[e,f,g] = meshgrid(nn, nn, nn); % convert vectors to
% 3D matrices
c=r.^(abs(e)+abs(f)+abs(g)); % Equation 9
e= b.*c; % Equation 10
h=full(sparse(time(:),1,e(:))); % Equivalent to
% equation 11
h=h/max(abs(h)); % Normalize h
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Reverberation Algorithms in Matlab
References : Posted by Gautham J. Mysore (gauthamjm [AT] yahoo [DOT] com)
Linked file : MATLABReverb.zip
Notes :
These
M-files implement a few reverberation algorithms (based on Schroeder's
and Moorer's algorithms). Each of the M-files include a short
description.
There are 5 M-files that implement reverberation. They are:
- schroeder1.m
- schroeder2.m
- schroeder3.m
- moorer.m
- stereoverb.m
The remaining 8 M-files implement filters, delay lines etc. Most of
these are used in the above M-files. They can also be used as building
blocks for other reverberation algorithms.
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Reverberation techniques
References : Posted by Sean Costello
Notes :
*
Parallel comb filters, followed by series allpass filters. This was the
original design by Schroeder, and was extended by Moorer. Has a VERY
metallic sound for sharp transients.
* Several allpass filters in serie (also proposed by Schroeder). Also suffers from metallic sound.
* 2nd-order comb and allpass filters (described by Moorer). Not supposed to give much of an advantage over first order sections.
* Nested allpass filters, where an allpass filter will replace the
delay line in another allpass filter. Pioneered by Gardner. Haven't
heard the results.
* Strange allpass amp delay line based structure in Jon Dattorro
article (JAES). Four allpass filters are used as an input to a cool
"figure-8" feedback loop, where four allpass reverberators are used in
series with
a few delay lines. Outputs derived from various taps in structure.
Supposedly based on a Lexicon reverb design. Modulating delay lines are
used in some of the allpass structures to "spread out" the eigentones.
* Feedback Delay Networks. Pioneered by Puckette/Stautner, with Jot
conducting extensive recent research. Sound VERY good, based on initial
experiments. Modulating delay lines and feedback matrixes used to
spread out eigentones.
* Waveguide-based reverbs, where the reverb structure is based upon the
junction of many waveguides. Julius Smith developed these. Recently,
these have been shown to be essentially equivalent to the feedback
delay network reverbs. Also sound very nice. Modulating delay lines and
scattering values used to spread out eigentones.
* Convolution-based reverbs, where the sound to be reverbed is
convolved with the impulse response of a room, or with
exponentially-decaying white noise. Supposedly the best sound, but very
computationally expensive, and not very flexible.
* FIR-based reverbs. Essentially the same as convolution. Probably not
used, but shorter FIR filters are probably used in combination with
many of the above techniques, to provide early reflections.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Saturation
Type : Waveshaper
References : Posted by Bram
Notes :
when
the input is below a certain threshold (t) these functions return the
input, if it goes over that threshold, they return a soft shaped
saturation.
Neigther claims to be fast ;-)
Code :
float saturate(float x, float t)
{
if(fabs(x)<t)
return x
else
{
if(x > 0.f);
return t + (1.f-t)*tanh((x-t)/(1-t));
else
return -(t + (1.f-t)*tanh((-x-t)/(1-t)));
}
}
or
float sigmoid(x)
{
if(fabs(x)<1)
return x*(1.5f - 0.5f*x*x);
else
return x > 0.f ? 1.f : -1.f;
}
float saturate(float x, float t)
{
if(abs(x)<t)
return x
else
{
if(x > 0.f);
return t + (1.f-t)*sigmoid((x-t)/((1-t)*1.5f));
else
return -(t + (1.f-t)*sigmoid((-x-t)/((1-t)*1.5f)));
}
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
SawSin
Type : Oscillator shape
References : Posted by Alexander Kritov
Code :
double sawsin(double x)
{
double t = fmod(x/(2*M_PI),(double)1.0);
if (t>0.5)
return -sin(x);
if (t<=0.5)
return (double)2.0*t-1.0;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Simple peak follower
Type : amplitude analysis
References : Posted by Phil Burk
Notes :
This
simple peak follower will give track the peaks of a signal. It will
rise rapidly when the input is rising, and then decay exponentially
when the input drops. It can be used to drive VU meters, or used in an
automatic gain control circuit.
Code :
// halfLife = time in seconds for output to decay to half value after an impulse
static float output = 0.0;
float scalar = pow( 0.5, 1.0/(halfLife * sampleRate)));
if( input < 0.0 )
input = -input; /* Absolute value. */
if ( input >= output )
{
/* When we hit a peak, ride the peak to the top. */
output = input;
}
else
{
/* Exponential decay of output when signal is low. */
output = output * scalar;
/*
** When current gets close to 0.0, set current to 0.0 to prevent FP underflow
** which can cause a severe performance degradation due to a flood
** of interrupts.
*/
if( output < VERY_SMALL_FLOAT ) output = 0.0;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Sin(x) Aproximation (with SSE code)
References : Posted by williamk[AT]wusik[DOT]com
Notes :
Sin Aproximation: sin(x) = x + ( x * (-x * x / 6));
This is very handy and fast, but not precise. Below you will find a simple SSE code.
Remember that all movaps command requires 16 bit aligned variables.
Code :
SSE code for computing only ONE value (scalar)
Replace all "ss" with "ps" if you want to calculate 4 values. And instead of "movps" use "movaps".
movss xmm1, xmm0 ; xmm0 = x
mulss xmm1, Filter_GenVal[k_n1] ; * -1
mulss xmm1, xmm0 ; -x * x
divss xmm1, Filter_GenVal[k_6] ; / 6
mulss xmm1, xmm0
addss xmm0, xmm1
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Sin, Cos, Tan approximation
References : http://www.wild-magic.com
Linked file : approx.h
Notes :
Code for approximation of cos, sin, tan and inv sin, etc.
Surprisingly accurate and very usable.
[edit by bram]
this code is taken literaly from
http://www.wild-magic.com/SourceCode.html
Go have a look at the MgcMath.h and MgcMath.cpp files in their library...
[/edit]
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Sine calculation
Type : waveform generation, Taylor approximation of sin()
References : Posted by Phil Burk
Notes :
Code
from JSyn for a sine wave generator based on a Taylor Expansion. It is
not as efficient as the filter methods, but it has linear frequency
control and is, therefore, suitable for FM or other time varying
applications where accurate frequency is needed. The sine generated is
accurate to at least 16 bits.
Code :
for(i=0; i < nSamples ; i++)
{
//Generate sawtooth phasor to provide phase for sine generation
IncrementWrapPhase(phase, freqPtr[i]);
//Wrap phase back into region where results are more accurate
if(phase > 0.5)
yp = 1.0 - phase;
else
{
if(phase < -0.5)
yp = -1.0 - phase;
else
yp = phase;
}
x = yp * PI;
x2 = x*x;
//Taylor expansion out to x**9/9! factored into multiply-adds
fastsin = x*(x2*(x2*(x2*(x2*(1.0/362880.0)
- (1.0/5040.0))
+ (1.0/120.0))
- (1.0/6.0))
+ 1.0);
outPtr[i] = fastsin * amplPtr[i];
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Soft saturation
Type : waveshaper
References : Posted by Bram de Jong
Notes :
This only works for positive values of x. a should be in the range 0..1
Code :
x < a:
f(x) = x
x > a:
f(x) = a + (x-a)/(1+((x-a)/(1-a))^2)
x > 1:
f(x) = (a+1)/2
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Square Waves
Type : waveform generation
References : Posted by Sean Costello
Notes :
One way to do a square wave:
You need two buzz generators (see Dodge & Jerse, or the Csound
source code, for implementation details). One of the buzz generators
runs at the desired square wave frequency, while the second buzz
generator is exactly one octave above this pitch. Subtract the higher
octave buzz generator's output from the lower buzz generator's output -
the result should be a signal with all odd harmonics, all at equal
amplitude. Filter the resultant signal (maybe integrate it). Voila, a
bandlimited square wave! Well, I think it should work...
The one question I have with the above technique is whether it produces
a waveform that truly resembles a square wave in the time domain. Even
if the number of harmonics, and the relative ratio of the harmonics, is
identical to an "ideal" bandwidth-limited square wave, it may have an
entirely different waveshape. No big deal, unless the signal is
processed by a nonlinearity, in which case the results of the nonlinear
processing will be far different than the processing of a waveform that
has a similar shape to a square wave.
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State variable
Type : 12db resonant low, high or bandpass
References : Effect Deisgn Part 1, Jon Dattorro, J. Audio Eng. Soc., Vol 45, No. 9, 1997 September
Notes :
Digital approximation of Chamberlin two-pole low pass. Easy to calculate coefficients, easy to process algorithm.
Code :
cutoff = cutoff freq in Hz
fs = sampling frequency //(e.g. 44100Hz)
f = 2 sin (pi * cutoff / fs) //[approximately]
q = resonance/bandwidth [0 < q <= 1] most res: q=1, less: q=0
low = lowpass output
high = highpass output
band = bandpass output
notch = notch output
scale = q
low=high=band=0;
//--beginloop
low = low + f * band;
high = scale * input - low - q*band;
band = f * high + band;
notch = high + low;
//--endloop
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State Variable Filter (Chamberlin version)
References : Hal Chamberlin, "Musical Applications of Microprocessors," 2nd Ed, Hayden Book Company 1985. pp 490-492.
Code :
//Input/Output
I - input sample
L - lowpass output sample
B - bandpass output sample
H - highpass output sample
N - notch output sample
F1 - Frequency control parameter
Q1 - Q control parameter
D1 - delay associated with bandpass output
D2 - delay associated with low-pass output
// parameters:
Q1 = 1/Q
// where Q1 goes from 2 to 0, ie Q goes from .5 to infinity
// simple frequency tuning with error towards nyquist
// F is the filter's center frequency, and Fs is the sampling rate
F1 = 2*pi*F/Fs
// ideal tuning:
F1 = 2 * sin( pi * F / Fs )
// algorithm
// loop
L = D2 + F1 * D1
H = I - L - Q1*D1
B = F1 * H + D1
N = H + L
// store delays
D1 = B
D2 = L
// outputs
L,H,B,N
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State Variable Filter (Double Sampled, Stable)
Type : 2 Pole Low, High, Band, Notch and Peaking
References : Posted by Andrew Simper
Notes :
Thanks to Laurent de Soras for the stability limit
and Steffan Diedrichsen for the correct notch output.
Code :
input = input buffer;
output = output buffer;
fs = sampling frequency;
fc = cutoff frequency normally something like:
440.0*pow(2.0, (midi_note - 69.0)/12.0);
res = resonance 0 to 1;
drive = internal distortion 0 to 0.1
freq = 2.0*sin(PI*MIN(0.25, fc/(fs*2))); // the fs*2 is because it's double sampled
damp = MIN(2.0*(1.0 - pow(res, 0.25)), MIN(2.0, 2.0/freq - freq*0.5));
notch = notch output
low = low pass output
high = high pass output
band = band pass output
peak = peaking output = low - high
--
double sampled svf loop:
for (i=0; i<numSamples; i++)
{
in = input[i];
notch = in - damp*band;
low = low + freq*band;
high = notch - low;
band = freq*high + band - drive*band*band*band;
out = 0.5*(notch or low or high or band or peak);
notch = in - damp*band;
low = low + freq*band;
high = notch - low;
band = freq*high + band - drive*band*band*band;
out += 0.5*(same out as above);
output[i] = out;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Stereo Enhancer
References : Posted by kurmisk[at]inbox[DOT]lv
Notes :
Stereo Enhanca
Code :
// WideCoeff 0.0 .... 1.5
#define StereoEnhanca(SamplL,SamplR,MonoSign, \
DeltaLeft,WideCoeff ) \
MonoSign = (SamplL + SamplR)/2.0; \
DeltaLeft = SamplL - MonoSign; \
DeltaLeft = DeltaLeft * WideCoeff; \
SamplL=SamplL + DeltaLeft; \
SamplR=SamplR - DeltaLeft;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Stilson's Moog filter code
Type : 4-pole LP, with fruity BP/HP
References : Posted by DFL
Notes :
Mind your p's and Q's...
This code was borrowed from Tim Stilson, and rewritten by me into a pd extern (moog~) available here:
http://www-ccrma.stanford.edu/~dfl/pd/index.htm
I ripped out the essential code and pasted it here...
Code :
WARNING: messy code follows ;)
// table to fixup Q in order to remain constant for various pole
frequencies, from Tim Stilson's code @ CCRMA (also in CLM distribution)
static float gaintable[199] = { 0.999969, 0.990082, 0.980347, 0.970764,
0.961304, 0.951996, 0.94281, 0.933777, 0.924866, 0.916077, 0.90741,
0.898865, 0.89044
2, 0.882141 , 0.873962, 0.865906, 0.857941, 0.850067, 0.842346,
0.834686, 0.827148, 0.819733, 0.812378, 0.805145, 0.798004, 0.790955,
0.783997, 0.77713, 0.77
0355, 0.763672, 0.75708 , 0.75058, 0.744141, 0.737793, 0.731537,
0.725342, 0.719238, 0.713196, 0.707245, 0.701355, 0.695557, 0.689819,
0.684174, 0.678558, 0.
673035, 0.667572, 0.66217, 0.65686, 0.651581, 0.646393, 0.641235,
0.636169, 0.631134, 0.62619, 0.621277, 0.616425, 0.611633, 0.606903,
0.602234, 0.597626, 0.
593048, 0.588531, 0.584045, 0.579651, 0.575287 , 0.570953, 0.566681,
0.562469, 0.558289, 0.554169, 0.550079, 0.546051, 0.542053, 0.538116,
0.53421, 0.530334,
0.52652, 0.522736, 0.518982, 0.515289, 0.511627, 0.507996 , 0.504425,
0.500885, 0.497375, 0.493896, 0.490448, 0.487061, 0.483704, 0.480377,
0.477081, 0.4738
16, 0.470581, 0.467377, 0.464203, 0.46109, 0.457977, 0.454926,
0.451874, 0.448883, 0.445892, 0.442932, 0.440033, 0.437134, 0.434265,
0.431427, 0.428619, 0.42
5842, 0.423096, 0.42038, 0.417664, 0.415009, 0.412354, 0.409729,
0.407135, 0.404572, 0.402008, 0.399506, 0.397003, 0.394501, 0.392059,
0.389618, 0.387207, 0.
384827, 0.382477, 0.380127, 0.377808, 0.375488, 0.37323, 0.370972,
0.368713, 0.366516, 0.364319, 0.362122, 0.359985, 0.357849, 0.355713,
0.353607, 0.351532,
0.349457, 0.347412, 0.345398, 0.343384, 0.34137, 0.339417, 0.337463,
0.33551, 0.333588, 0.331665, 0.329773, 0.327911, 0.32605, 0.324188,
0.322357, 0.320557,
0.318756, 0.316986, 0.315216, 0.313446, 0.311707, 0.309998, 0.308289,
0.30658, 0.304901, 0.303223, 0.301575, 0.299927, 0.298309, 0.296692,
0.295074, 0.293488
, 0.291931, 0.290375, 0.288818, 0.287262, 0.285736, 0.284241, 0.282715,
0.28125, 0.279755, 0.27829, 0.276825, 0.275391, 0.273956, 0.272552,
0.271118, 0.26974
5, 0.268341, 0.266968, 0.265594, 0.264252, 0.262909, 0.261566, 0.260223, 0.258911, 0.257599, 0.256317, 0.255035, 0.25375 };
static inline float saturate( float input ) { //clamp without branching
#define _limit 0.95
float x1 = fabsf( input + _limit );
float x2 = fabsf( input - _limit );
return 0.5 * (x1 - x2);
}
static inline float crossfade( float amount, float a, float b ) {
return (1-amount)*a + amount*b;
}
//code for setting Q
float ix, ixfrac;
int ixint;
ix = x->p * 99;
ixint = floor( ix );
ixfrac = ix - ixint;
Q = resonance * crossfade( ixfrac, gaintable[ ixint + 99 ], gaintable[ ixint + 100 ] );
//code for setting pole coefficient based on frequency
float fc = 2 * frequency / x->srate;
float x2 = fc*fc;
float x3 = fc*x2;
p = -0.69346 * x3 - 0.59515 * x2 + 3.2937 * fc - 1.0072; //cubic fit by DFL, not 100% accurate but better than nothing...
}
process loop:
float state[4], output; //should be global scope / preserved between calls
int i,pole;
float temp, input;
for ( i=0; i < numSamples; i++ ) {
input = *(in++);
output = 0.25 * ( input - output ); //negative feedback
for( pole = 0; pole < 4; pole++) {
temp = state[pole];
output = saturate( output + p * (output - temp));
state[pole] = output;
output = saturate( output + temp );
}
lowpass = output;
highpass = input - output;
bandpass = 3 * x->state[2] - x->lowpass; //got this one from paul kellet
*out++ = lowpass;
output *= Q; //scale the feedback
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Time domain convolution with O(n^log2(3))
References : Wilfried Welti
Notes :
[Quoted from Wilfrieds mail...]
I found last weekend that it is possible to do convolution in time
domain (no complex numbers, 100% exact result with int) with
O(n^log2(3)) (about O(n^1.58)).
Due to smaller overhead compared to FFT-based convolution, it should be the fastest algorithm for medium sized FIR's.
Though, it's slower as FFT-based convolution for large n.
It's pretty easy:
Let's say we have two finite signals of length 2n, which we want
convolve : A and B. Now we split both signals into parts of size n, so
we get A = A1 + A2, and B = B1 +B2.
Now we can write:
(1) A*B = (A1+A2)*(B1+B2) = A1*B1 + A2*B1 + A1*B2 + A2*B2
where * means convolution.
This we knew already: We can split a convolution into four convolutions of halved size.
Things become interesting when we start shifting blocks in time:
Be z a signal which has the value 1 at x=1 and zero elsewhere.
Convoluting a signal X with z is equivalent to shifting X by one
rightwards. When I define z^n as n-fold convolution of z with itself,
like: z^1 = z, z^2 = z*z, z^0 = z shifted leftwards by 1 = impulse at
x=0, and so on, I can use it to shift signals:
X * z^n means shifting the signal X by the value n rightwards.
X * z^-n means shifting the signal X by the value n leftwards.
Now we look at the following term:
(2) (A1 + A2 * z^-n) * (B1 + B2 * z^-n)
This is a convolution of two blocks of size n: We shift A2 by n leftwards so it completely overlaps A1, then we add them.
We do the same thing with B1 and B2. Then we convolute the two resulting blocks.
now let's transform this term:
(3) (A1 + A2 * z^-n) * (B1 + B2 * z^-n)
= A1*B1 + A1*B2*z^-n + A2*z^-n*B1 + A2*z^ n*B2*z^-n
= A1*B1 + (A1*B2 + A2*B1)*z^-n + A2*B2*z^-2n
(4) (A1 + A2 * z^-n) * (B1 + B2 * z^-n) - A1*B1 - A2*B2*z^-2n
= (A1*B2 + A2*B1)*z^-n
Now we convolute both sides of the equation (4) by z^n:
(5) (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n - A1*B1*z^n - A2*B2*z^-n
= (A1*B2 + A2*B1)
Now we see that the right part of equation (5) appears within equation
(1), so we can replace this appearance by the left part of eq (5).
(6) A*B = (A1+A2)*(B1+B2) = A1*B1 + A2*B1 + A1*B2 + A2*B2
= A1*B1
+ (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n - A1*B1*z^n - A2*B2*z^-n
+ A2*B2
Voila!
We have constructed the convolution of A*B with only three convolutions
of halved size. (Since the convolutions with z^n and z^-n are only
shifts
of blocks with size n, they of course need only n operations for processing :)
This can be used to construct an easy recursive algorithm of Order O(n^log2(3))
Code :
void convolution(value* in1, value* in2, value* out, value* buffer, int size)
{
value* temp1 = buffer;
value* temp2 = buffer + size/2;
int i;
// clear output.
for (i=0; i<size*2; i++) out[i] = 0;
// Break condition for recursion: 1x1 convolution is multiplication.
if (size == 1)
{
out[0] = in1[0] * in2[0];
return;
}
// first calculate (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n
signal_add(in1, in1+size/2, temp1, size/2);
signal_add(in2, in2+size/2, temp2, size/2);
convolution(temp1, temp2, out+size/2, buffer+size, size/2);
// then add A1*B1 and substract A1*B1*z^n
convolution(in1, in2, temp1, buffer+size, size/2);
signal_add_to(out, temp1, size);
signal_sub_from(out+size/2, temp1, size);
// then add A2*B2 and substract A2*B2*z^-n
convolution(in1+size/2, in2+size/2, temp1, buffer+size, size/2);
signal_add_to(out+size, temp1, size);
signal_sub_from(out+size/2, temp1, size);
}
"value" may be a suitable type like int or float.
Parameter "size" is the size of the input signals and must be a power of 2. out and buffer must point to arrays of size 2*n.
Just to be complete, the helper functions:
void signal_add(value* in1, value* in2, value* out, int size)
{
int i;
for (i=0; i<size; i++) out[i] = in1[i] + in2[i];
}
void signal_sub_from(value* out, value* in, int size)
{
int i;
for (i=0; i<size; i++) out[i] -= in[i];
}
void signal_add_to(value* out, value* in, int size)
{
int i;
for (i=0; i<size; i++) out[i] += in[i];
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Time domain convolution with O(n^log2(3))
References : Posted by Magnus Jonsson
Notes :
[see other code by Wilfried Welti too!]
Code :
void mul_brute(float *r, float *a, float *b, int w)
{
for (int i = 0; i < w+w; i++)
r[i] = 0;
for (int i = 0; i < w; i++)
{
float *rr = r+i;
float ai = a[i];
for (int j = 0; j < w; j++)
rr[j] += ai*b[j];
}
}
// tmp must be of length 2*w
void mul_knuth(float *r, float *a, float *b, int w, float *tmp)
{
if (w < 30)
{
mul_brute(r, a, b, w);
}
else
{
int m = w>>1;
for (int i = 0; i < m; i++)
{
r[i ] = a[m+i]-a[i ];
r[i+m] = b[i ]-b[m+i];
}
mul_knuth(tmp, r , r+m, m, tmp+w);
mul_knuth(r , a , b , m, tmp+w);
mul_knuth(r+w, a+m, b+m, m, tmp+w);
for (int i = 0; i < m; i++)
{
float bla = r[m+i]+r[w+i];
r[m+i] = bla+r[i ]+tmp[i ];
r[w+i] = bla+r[w+m+i]+tmp[i+m];
}
}
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
tone detection with Goertzel
Type : Goertzel
References : Posted by espenr[AT]ii[DOT]uib[DOT]no
Linked file : http://www.ii.uib.no/~espenr/tonedetect.zip
Notes :
Goertzel
is basically DFT of parts of a spectrum not the total spectrum as you
normally do with FFT. So if you just want to check out the power for
some frequencies this could be better. Is good for DTFM detection I've
heard.
The WNk isn't calculated 100% correctly, but it seems to work so ;)
Yeah and the code is C++ so you might have to do some small adjustment
to compile it as C.
Code :
/** Tone detect by Goertzel algorithm
*
* This program basically searches for tones (sines) in a sample and reports the different dB it finds for
* different frequencies. Can easily be extended with some thresholding to report true/false on detection.
* I'm far from certain goertzel it implemented 100% correct, but it works :)
*
* Hint, the SAMPLERATE, BUFFERSIZE, FREQUENCY, NOISE and SIGNALVOLUME all affects the outcome of the reported dB. Tweak
* em to find the settings best for your application. Also, seems to be pretty sensitive to noise (whitenoise anyway) which
* is a bit sad. Also I don't know if the goertzel really likes float values for the frequency ... And using 44100 as
* samplerate for detecting 6000 Hz tone is kinda silly I know :)
*
* Written by: Espen Riskedal, espenr@ii.uib.no, july-2002
*/
#include <iostream>
#include <cmath>
#include <cstdlib>
using std::rand;
// math stuff
using std::cos;
using std::abs;
using std::exp;
using std::log10;
// iostream stuff
using std::cout;
using std::endl;
#define PI 3.14159265358979323844
// change the defines if you want to
#define SAMPLERATE 44100
#define BUFFERSIZE 8820
#define FREQUENCY 6000
#define NOISE 0.05
#define SIGNALVOLUME 0.8
/** The Goertzel algorithm computes the k-th DFT coefficient of the input signal using a second-order filter.
* http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html.
* Basiclly it just does a DFT of the frequency we want to check, and none of the others (FFT calculates for all frequencies).
*/
float goertzel(float *x, int N, float frequency, int samplerate) {
float Skn, Skn1, Skn2;
Skn = Skn1 = Skn2 = 0;
for (int i=0; i<N; i++) {
Skn2 = Skn1;
Skn1 = Skn;
Skn = 2*cos(2*PI*frequency/samplerate)*Skn1 - Skn2 + x[i];
}
float WNk = exp(-2*PI*frequency/samplerate); // this one ignores complex stuff
//float WNk = exp(-2*j*PI*k/N);
return (Skn - WNk*Skn1);
}
/** Generates a tone of the specified frequency
* Gotten
from:
http://groups.google.com/groups?hl=en&lr=&ie=UTF-8&oe=UTF-8&safe=off&selm=3c641e%243jn%40uicsl.csl.uiuc.edu
*/
float *makeTone(int samplerate, float frequency, int length, float gain=1.0) {
//y(n) = 2 * cos(A) * y(n-1) - y(n-2)
//A= (frequency of interest) * 2 * PI / (sampling frequency)
//A is in radians.
// frequency of interest MUST be <= 1/2 the sampling frequency.
float *tone = new float[length];
float A = frequency*2*PI/samplerate;
for (int i=0; i<length; i++) {
if (i > 1) tone[i]= 2*cos(A)*tone[i-1] - tone[i-2];
else if (i > 0) tone[i] = 2*cos(A)*tone[i-1] - (cos(A));
else tone[i] = 2*cos(A)*cos(A) - cos(2*A);
}
for (int i=0; i<length; i++) tone[i] = tone[i]*gain;
return tone;
}
/** adds whitenoise to a sample */
void *addNoise(float *sample, int length, float gain=1.0) {
for (int i=0; i<length; i++) sample[i] += (2*(rand()/(float)RAND_MAX)-1)*gain;
}
/** returns the signal power/dB */
float power(float value) {
return 20*log10(abs(value));
}
int main(int argc, const char* argv) {
cout << "Samplerate: " << SAMPLERATE << "Hz\n";
cout << "Buffersize: " << BUFFERSIZE << " samples\n";
cout << "Correct frequency is: " << FREQUENCY << "Hz\n";
cout << " - signal volume: " << SIGNALVOLUME*100 << "%\n";
cout << " - white noise: " << NOISE*100 << "%\n";
float *tone = makeTone(SAMPLERATE, FREQUENCY, BUFFERSIZE, SIGNALVOLUME);
addNoise(tone, BUFFERSIZE,NOISE);
int stepsize = FREQUENCY/5;
for (int i=0; i<10; i++) {
int freq = stepsize*i;
cout << "Trying freq: " << freq <<
"Hz -> dB: " << power(goertzel(tone, BUFFERSIZE, freq,
SAMPLERATE)) << endl;
}
delete tone;
return 0;
}
7 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Tone detection with Goertzel (x86 ASM)
Type : Tone detection with Goertzel in x86 assembly
References : Posted by Christian[AT]savioursofsoul[DOT]de
Notes :
This is an "assemblified" version of the Goertzel Tone Detector. It is about 2 times faster than the original code.
The code has been tested and it works fine.
Hope you can use it. I'm gonna try to build a Tuner (as VST-Plugin). I
hope, that this will work :-\ If anyone is intrested, please let me
know.
Christian
Code :
function Goertzel_x87(Buffer :Psingle; BLength:Integer; frequency: Single; samplerate: Single):Single;
asm
mov ecx,BLength
mov eax,Buffer
fld x2
fldpi
fmulp
fmul frequency
fdiv samplerate
fld st(0)
fcos
fld x2
fmulp
fxch st(1)
fldz
fsub st(0),st(1)
fstp st(1)
fldl2e
fmul
fld st(0)
frndint
fsub st(1),st(0)
fxch st(1)
f2xm1
fld1
fadd
fscale
fstp st(1)
fldz
fldz
fldz
@loopStart:
fxch st(1)
fxch st(2)
fstp st(0)
fld st(3)
fmul st(0),st(1)
fsub st(0),st(2)
fld [eax].Single
faddp
add eax,4
loop @loopStart
@loopEnd:
fxch st(3)
fmulp st(2), st(0)
fsub st(0),st(1)
fstp result
ffree st(2)
ffree st(1)
ffree st(0)
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
transistor differential amplifier simulation
Type : Waveshaper
References : Posted by Christian[at]savioursofsoul[dot]de
Notes :
Writting
an exam about electronic components, i learned several equations about
simulating that stuff. One simplified equation was the tanh(x) formula
for the differential amplifier. It is not exact, but since the
amplifiers are driven with only small amplitudes the behaviour is most
often even advanced linear.
The fact, that the amp is differential, means, that the 2n order is eliminated, so the sound is also similar to a tube.
For a very fast use, this code is in pure assembly language (not
optimized with SSE-Code yet) and performs in VST-Plugins very fast.
The code was written in delphi and if you want to translate the
assembly code, you should know, the the parameters passing is done via
registers. So pinp=EAX pout=EDX sf=ECX.
Code :
procedure Transistor(pinp,pout : PSingle; sf:Integer; Faktor: Single);
asm
fld Faktor
@Start:
fld [eax].single
fmul st(0),st(1)
fldl2e
fmul
fld st(0)
frndint
fsub st(1),st
fxch st(1)
f2xm1
fld1
fadd
fscale { result := z * 2**i }
fstp st(1)
fld st(0)
fmulp
fld st(0)
fld1
faddp
fld1
fsubp st(2),st(0)
fdivp
fstp [edx].single
add eax,4
add edx,4
loop @Start
fstp st(0)
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Variable-hardness clipping function
References : Posted by Laurent de Soras
Linked file : laurent.gif
Notes :
k >= 1 is the "clipping hardness". 1 gives a smooth clipping, and a high value gives hardclipping.
Don't set k too high, because the formula use the pow() function, which
use exp() and would overflow easily. 100 seems to be a reasonable value
for "hardclipping"
Code :
f (x) = sign (x) * pow (atan (pow (abs (x), k)), (1 / k));
4 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Waveform generator using MinBLEPS
References : Posted by locke[AT]rpgfan.demon.co.uk
Linked file : MinBLEPS.zip
Notes :
C code and project file for MSVC6 for a bandwidth-limited saw/square (with PWM) generator using MinBLEPS.
This code is based on Eli's MATLAB MinBLEP code and uses his original minblep.mat file.
Instead of keeping a list of all active MinBLEPS, the output of each
MinBLEP is stored in a buffer, in which all consequent MinBLEPS and the
waveform output are added together. This optimization makes it fast
enough to be used realtime.
Produces slight aliasing when sweeping high frequencies. I don't know
wether Eli's original code does the same, because I don't have MATLAB.
Any help would be appreciated.
The project name is 'hardsync', because it's easy to generate hardsync using MinBLEPS.
Code :
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
WaveShaper
Type : waveshaper
References : Posted by Bram de Jong
Notes :
where x (in [-1..1] will be distorted and a is a distortion parameter that goes from 1 to infinity
The equation is valid for positive and negativ values.
If a is 1, it results in a slight distortion and with bigger a's the signal get's more funky.
A good thing about the shaper is that feeding it with bigger-than-one
values, doesn't create strange fx. The maximum this function will reach is
1.2 for a=1.
Code :
f(x,a) = x*(abs(x) + a)/(x^2 + (a-1)*abs(x) + 1)
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Waveshaper
Type : waveshaper
References : Posted by Jon Watte
Notes :
A favourite of mine is using a sin() function instead.
This will have the "unfortunate" side effect of removing
odd harmonics if you take it to the extreme: a triangle
wave gets mapped to a pure sine wave.
This will work with a going from .1 or so to a= 5 and bigger!
The mathematical limits for a = 0 actually turns it into a linear
function at that point, but unfortunately FPUs aren't that good
with calculus :-) Once a goes above 1, you start getting clipping
in addition to the "soft" wave shaping. It starts getting into
more of an effect and less of a mastering tool, though :-)
Seeing as this is just various forms of wave shaping, you
could do it all with a look-up table, too. In my version, that would
get rid of the somewhat-expensive sin() function.
Code :
(input: a == "overdrive amount")
z = M_PI * a;
s = 1/sin(z)
b = 1/a
if (x > b)
f(x) = 1
else
f(x) = sin(z*x)*s
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Waveshaper
References : Posted by Partice Tarrabia and Bram de Jong
Notes :
amount should be in [-1..1[ Plot it and stand back in astonishment! ;)
Code :
x = input in [-1..1]
y = output
k = 2*amount/(1-amount);
f(x) = (1+k)*x/(1+k*abs(x))
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Waveshaper (simple description)
Type : Polynomial; Distortion
References : Posted by Jon Watte
Notes :
> The other question; what's a 'waveshaper' algorithm. Is it simply another
> word for distortion?
A typical "waveshaper" is some function which takes an input sample value
X and transforms it to an output sample X'. A typical implementation would
be a look-up table of some number of points, and some level of interpolation
between those points (say, cubic). When people talk about a wave shaper,
this is most often what they mean. Note that a wave shaper, as opposed to a
filter, does not have any state. The mapping from X -> X' is stateless.
Some wave shapers are implemented as polynomials, or using other math
functions. Hard clipping is a wave shaper implemented using the min() and
max() functions (or the three-argument clamp() function, which is the same
thing). A very mellow and musical-sounding distortion is implemented using
a third-degree polynomial; something like X' = (3/2)X - (1/2)X^3. The nice
thing with polynomial wave shapers is that you know that the maximum they
will expand bandwidth is their order. Thus, you need to oversample 3x to
make sure that a third-degree polynomial is aliasing free. With a lookup
table based wave shaper, you don't know this (unless you treat an N-point
table as an N-point polynomial :-)
Code :
float waveshape_distort( float in ) {
return 1.5f * in - 0.5f * in *in * in;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Weird synthesis
References : Posted by Andy M00cho
Notes :
(quoted from Andy's mail...)
What I've done in a soft-synth I've been working on is used what I've
termed Fooglers, no reason, just liked the name :) Anyway all I've done
is use a *VERY* short delay line of 256 samples and then use 2
controllable taps into the delay with High Frequency Damping, and a
feedback parameter.
Using a tiny fixed delay size of approx. 4.8ms (really 256 samples/1k
memory with floats) means this costs, in terms of cpu consumption
practically nothing, and the filter is a real simple 1 pole low-pass
filter. Maybe not DSP'litically correct but all I wanted was to avoid
the high frequencies trashing the delay line when high feedbacks
(99%->99.9%) are used (when the fun starts ;).
I've been getting some really sexy sounds out of this idea, and of
course you can have the delay line tuneable if you choose to use
fractional taps, but I'm happy with it as it is.. 1 nice simple, yet
powerful addition to the base oscillators.
In reality you don't need 2 taps, but I found that using 2 added that extra element of funkiness...
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Zoelzer biquad filters
Type : biquad IIR
References : Udo Zoelzer: Digital Audio Signal Processing (John Wiley & Sons, ISBN 0 471 97226 6), Chris Townsend
Notes :
Here's the formulas for the Low Pass, Peaking, and Low Shelf, which should
cover the basics. I tried to convert the formulas so they are little more consistent.
Also, the Zolzer low pass/shelf formulas didn't have adjustable Q, so I added that for
consistency with Roberts formulas as well. I think someone may want to check that I did
it right.
------------ Chris Townsend
I mistranscribed the low shelf cut formulas.
Hopefully this is correct. Thanks to James McCartney for noticing.
------------ Chris Townsend
Code :
omega = 2*PI*frequency/sample_rate
K=tan(omega/2)
Q=Quality Factor
V=gain
LPF: b0 = K^2
b1 = 2*K^2
b2 = K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
peakingEQ:
boost:
b0 = 1 + V*K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - V*K/Q + K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
cut:
b0 = 1 + K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - K/Q + K^2
a0 = 1 + V*K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - V*K/Q + K^2
lowShelf:
boost:
b0 = 1 + sqrt(2*V)*K + V*K^2
b1 = 2*(V*K^2 - 1)
b2 = 1 - sqrt(2*V)*K + V*K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
cut:
b0 = 1 + K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - K/Q + K^2
a0 = 1 + sqrt(2*V)*K + V*K^2
a1 = 2*(v*K^2 - 1)
a2 = 1 - sqrt(2*V)*K + V*K^2
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
