|
|
|
|
 |
Entries in this category
VST SDK GUI Switch without
(Allmost) Ready-to-use oscillators
1 pole LPF for smooth parameter changes
1-RC and C filter
16-Point Fast Integer Sinc Interpolator.
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
3 Band Equaliser
303 type filter with saturation
3rd order Spline interpollation
4th order Linkwitz-Riley filters
5-point spline interpollation
Alias-free waveform generation with analog filtering
Alien Wah
All-Pass Filters, a good explanation
Allocating aligned memory
AM Formantic Synthesis
Another 4-pole lowpass...
Another cheap sinusoidal LFO
another LFO class
Antialiased Lines
antialiased square generator
Arbitary shaped band-limited waveform generation (using oversampling and low-pass filtering)
Audiable alias free waveform gen using width sine
Automatic PDC system
Band Limited PWM Generator
Band Limited waveforms my way
Bandlimited sawtooth synthesis
Bandlimited waveform generation
Bandlimited waveform generation with hard sync
Bandlimited waveforms synopsis.
Bandlimited waveforms...
Base-2 exp
Bass Booster
Beat Detector Class
Biquad C code
Biquad, Butterworth, Chebyshev N-order, M-channel optimized filters
Bit quantization/reduction effect
Bit-Reversed Counting
Block/Loop Benchmarking
Branchless Clipping
Butterworth
Butterworth Optimized C++ Class
C# Oscilator class
C++ class implementation of RBJ Filters
C++ gaussian noise generation
C-Weighed Filter
Calculate notes (java)
Cascaded resonant lp/hp filter
Center separation in a stereo mixdown
Center separation in a stereo mixdown
Cheap pseudo-sinusoidal lfo
chebyshev waveshaper (using their recursive definition)
Class for waveguide/delay effects
Clipping without branching
Coefficients for Daubechies wavelets 1-38
Compressor
Constant-time exponent of 2 detector
Conversion and normalization of 16-bit sample to a floating point number
Conversions on a PowerPC
Cool Sounding Lowpass With Decibel Measured Resonance
Copy-protection schemes
Cubic interpollation
Cubic polynomial envelopes
Cure for malicious samples
DC filter
Decimator
Delay time calculation for reverberation
Delphi Class implementation of the RBJ filters
Denormal DOUBLE variables, macro
Denormal numbers
Denormal numbers, the meta-text
Denormalization preventer
Denormalization preventer
DFT
Digital RIAA equalization filter coefficients
DIRAC - Free C/C++ Library for Time and Pitch Manipulation of Audio Based on Time-Frequency Transforms
Direct form II
Direct Form II biquad
Direct pink noise synthesis with auto-correlated generator
Discrete Summation Formula (DSF)
Dither code
Dithering
Double to Int
Drift generator
DSF (super-set of BLIT)
dynamic convolution
Early echo's with image-mirror technique
Easy noise generation
ECE320 project: Reverberation w/ parameter control from PC
Envelope detector
Envelope Detector class (C++)
Envelope Follower
Envelope follower with different attack and release
Exponential curve for
Exponential parameter mapping
Fast & small sine generation tutorial
fast abs/neg/sign for 32bit floats
Fast binary log approximations
Fast cube root, square root, and reciprocal for x86/SSE CPUs.
Fast Downsampling With Antialiasing
fast exp() approximations
Fast exp2 approximation
Fast Exponential Envelope Generator
Fast Float Random Numbers
Fast in-place Walsh-Hadamard Transform
Fast LFO in Delphi...
Fast log2
fast power and root estimates for 32bit floats
Fast rounding functions in pascal
Fast sign for 32 bit floats
Fast SIN approximation for usage in e.g. additive synthesizers
Fast sine and cosine calculation
Fast sine wave calculation
Fast square wave generator
Fast Whitenoise Generator
FFT
FFT classes in C++ and Object Pascal
Float to int
Float to int (more intel asm)
Float to integer conversion
Float-to-int, coverting an array of floats
fold back distortion
Formant filter
Frequency response from biquad coefficients
frequency warped FIR lattice
Gaussian dithering
Gaussian random numbers
Gaussian White noise
Gaussian White Noise
Generator
Guitar feedback
Hermite Interpolator (x86 ASM)
Hermite interpollation
Hilbert Filter Coefficient Calculation
Hiqh quality /2 decimators
Inverted parabolic envelope
Java FFT
Karlsen
Karlsen Fast Ladder
Linear interpolation
Lo-Fi Crusher
Lock free fifo
Look ahead limiting
Lookahead Limiter
Lowpass filter for parameter edge filtering
LP and HP filter
LPC analysis (autocorrelation + Levinson-Durbin recursion)
LPF 24dB/Oct
Magnitude and phase plot of arbitrary IIR function, up to 5th order
Matlab Time Domain Impulse Response Inverter/Divider
MATLAB-Tools for SNDAN
matlab/octave code for minblep table generation
MDCT and IMDCT based on FFTW3
Measuring interpollation noise
MIDI note/frequency conversion
Millimeter to DB (faders...)
Moog Filter
Moog VCF
Moog VCF, variation 1
Moog VCF, variation 2
Most simple and smooth feedback delay
Most simple static delay
Motorola 56300 Disassembler
Noise Shaping Class
Nonblocking multiprocessor/multithread algorithms in C++
Notch filter
One pole filter, LP and HP
One pole LP and HP
One pole, one zero LP/HP
One zero, LP/HP
One-Liner IIR Filters (1st order)
Output Limiter using Envelope Follower in C++
PADsynth synthesys method
Parabolic shaper
Parallel combs delay calculation
Peak/Notch filter
Perfect LP4 filter
Phase equalization
Phase modulation Vs. Frequency modulation
Phase modulation Vs. Frequency modulation II
Phaser code
Piecewise quadratic approximate exponential function
Pink noise filter
please add it as a comment to the Denormalization preventer entry (no comments are allowed now) thanks
Plot Filter (Analyze filter characteristics)
Plotting R B-J Equalisers in Excel
Polynominal Waveshaper
Polyphase Filters
Polyphase Filters (Delphi)
Poor Man's FIWIZ
pow(x,4) approximation
Prewarping
PRNG for non-uniformly distributed values from trigonometric identity
Pseudo-Random generator
PulseQuad
Pulsewidth modulation
QFT and DQFT (double precision) classes
Quick & Dirty Sine
quick and dirty sine generator
randja compressor
rational tanh approximation
RBJ Audio-EQ-Cookbook
RBJ Audio-EQ-Cookbook
RBJ Wavetable 101
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
Remez Exchange Algorithm (Parks/McClellan)
Remez Remez (Parks/McClellan)
resampling
Resonant filter
Resonant IIR lowpass (12dB/oct)
Resonant low pass filter
Reverb Filter Generator
Reverberation Algorithms in Matlab
Reverberation techniques
Rossler and Lorenz Oscillators
Saturation
SawSin
Simple biquad filter from apple.com
Simple Compressor class (C++)
Simple peak follower
Simple Tilt equalizer
Simple Time Stretching-Granular Synthesizer
Sin(x) Aproximation (with SSE code)
Sin, Cos, Tan approximation
Sine calculation
Smooth random LFO Generator
smsPitchScale Source Code
Soft saturation
Spuc's open source filters
Square Waves
State variable
State Variable Filter (Chamberlin version)
State Variable Filter (Double Sampled, Stable)
Stereo Enhancer
Stereo Field Rotation Via Transformation Matrix
Stereo Width Control (Obtained Via Transfromation Matrix)
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)
Trammell Pink Noise (C++ class)
transistor differential amplifier simulation
Type : LPF 24dB/Oct
Variable-hardness clipping function
Various Biquad filters
Vintage VU meters tutorial
Waveform generator using MinBLEPS
WaveShaper
Waveshaper
Waveshaper
Waveshaper (simple description)
Waveshaper :: Gloubi-boulga
Wavetable Synthesis
Weird synthesis
Windowed Sinc FIR Generator
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
3 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];
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
1 pole LPF for smooth parameter changes
Type : 1-pole LPF class
References : Posted by zioguido@gmail.com
Notes :
This is a very simple class that I'm using in my plugins for smoothing parameter changes that directly affect audio stream.
It's a 1-pole LPF, very easy on CPU.
Change the value of variable "a" (0~1) for slower or a faster response.
Of course you can also use it as a lowpass filter for audio signals.
Code :
class CParamSmooth
{
public:
CParamSmooth() { a = 0.99f; b = 1.f - a; z = 0; };
~CParamSmooth();
inline float Process(float in) { z = (in * b) + (z * a); return z; }
private:
float a, b, z;
};
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
1-RC and C filter
Type : Simple 2-pole LP
References : Posted by madbrain[AT]videotron[DOT]ca
Notes :
This filter is called 1-RC and C since it uses these two parameters. C and R correspond to raw cutoff and inverted resonance, and have a range from 0 to 1.
Code :
//Parameter calculation
//cutoff and resonance are from 0 to 127
c = pow(0.5, (128-cutoff) / 16.0);
r = pow(0.5, (resonance+24) / 16.0);
//Loop:
v0 = (1-r*c)*v0 - (c)*v1 + (c)*input;
v1 = (1-r*c)*v1 + (c)*v0;
output = v1;
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
16-Point Fast Integer Sinc Interpolator.
References : Posted by mumart[at]gmail[dot]com
Notes :
This is designed for fast upsampling with good quality using only a 32-bit accumulator. Sound quality is very good. Conceptually it resamples the input signal 32768x and performs nearest-neighbour to get the requested sample rate. As a result downsampling will result in aliasing.
The provided Sinc table is Blackman-Harris windowed with a slight lowpass. The table entries are 16-bit and are 16x linear-oversampled. It should be pretty easy to figure out how to make your own table for it.
Code provided is in Java. Converting to C/MMX etc. should be pretty trivial.
Remember the interpolator requires a number of samples before and after the sample to be interpolated, so you can't resample the whole of a passed input buffer in one go.
Have fun,
Martin
Code :
public class SincResampler {
private final int FP_SHIFT = 15;
private final int FP_ONE = 1 << FP_SHIFT;
private final int FP_MASK = FP_ONE - 1;
private final int POINT_SHIFT = 4; // 16 points
private final int OVER_SHIFT = 4; // 16x oversampling
private final short[] table = {
0, -7, 27, -71, 142, -227, 299, 32439, 299, -227, 142, -71, 27, -7, 0, 0,
0, 0, -5, 36, -142, 450, -1439, 32224, 2302, -974, 455, -190, 64, -15, 2, 0,
0, 6, -33, 128, -391, 1042, -2894, 31584, 4540, -1765, 786, -318, 105, -25, 3, 0,
0, 10, -55, 204, -597, 1533, -4056, 30535, 6977, -2573, 1121, -449, 148, -36, 5, 0,
-1, 13, -71, 261, -757, 1916, -4922, 29105, 9568, -3366, 1448, -578, 191, -47, 7, 0,
-1, 15, -81, 300, -870, 2185, -5498, 27328, 12263, -4109, 1749, -698, 232, -58, 9, 0,
-1, 15, -86, 322, -936, 2343, -5800, 25249, 15006, -4765, 2011, -802, 269, -68, 10, 0,
-1, 15, -87, 328, -957, 2394, -5849, 22920, 17738, -5298, 2215, -885, 299, -77, 12, 0,
0, 14, -83, 319, -938, 2347, -5671, 20396, 20396, -5671, 2347, -938, 319, -83, 14, 0,
0, 12, -77, 299, -885, 2215, -5298, 17738, 22920, -5849, 2394, -957, 328, -87, 15, -1,
0, 10, -68, 269, -802, 2011, -4765, 15006, 25249, -5800, 2343, -936, 322, -86, 15, -1,
0, 9, -58, 232, -698, 1749, -4109, 12263, 27328, -5498, 2185, -870, 300, -81, 15, -1,
0, 7, -47, 191, -578, 1448, -3366, 9568, 29105, -4922, 1916, -757, 261, -71, 13, -1,
0, 5, -36, 148, -449, 1121, -2573, 6977, 30535, -4056, 1533, -597, 204, -55, 10, 0,
0, 3, -25, 105, -318, 786, -1765, 4540, 31584, -2894, 1042, -391, 128, -33, 6, 0,
0, 2, -15, 64, -190, 455, -974, 2302, 32224, -1439, 450, -142, 36, -5, 0, 0,
0, 0, -7, 27, -71, 142, -227, 299, 32439, 299, -227, 142, -71, 27, -7, 0
};
/*
private final int POINT_SHIFT = 1; // 2 points
private final int OVER_SHIFT = 0; // 1x oversampling
private final short[] table = {
32767, 0,
0 , 32767
};
*/
private final int POINTS = 1 << POINT_SHIFT;
private final int INTERP_SHIFT = FP_SHIFT - OVER_SHIFT;
private final int INTERP_BITMASK = ( 1 << INTERP_SHIFT ) - 1;
/*
input - array of input samples
inputPos - sample position ( must be at least POINTS/2 + 1, ie. 7 )
inputFrac - fractional sample position ( 0 <= inputFrac < FP_ONE )
step - number of input samples per output sample * FP_ONE
lAmp - left output amplitude ( 1.0 = FP_ONE )
lBuf - left output buffer
rAmp - right output amplitude ( 1.0 = FP_ONE )
rBuf - right output buffer
pos - offset into output buffers
count - number of output samples to produce
*/
public void resample( short[] input, int inputPos, int inputFrac, int step,
int lAmp, int[] lBuf, int rAmp, int[] rBuf, int pos, int count ) {
for( int p = 0; p < count; p++ ) {
int tabidx1 = ( inputFrac >> INTERP_SHIFT ) << POINT_SHIFT;
int tabidx2 = tabidx1 + POINTS;
int bufidx = inputPos - POINTS/2 + 1;
int a1 = 0, a2 = 0;
for( int t = 0; t < POINTS; t++ ) {
a1 += table[ tabidx1++ ] * input[ bufidx ] >> 15;
a2 += table[ tabidx2++ ] * input[ bufidx ] >> 15;
bufidx++;
}
int out = a1 + ( ( a2 - a1 ) * ( inputFrac & INTERP_BITMASK ) >> INTERP_SHIFT );
lBuf[ pos ] += out * lAmp >> FP_SHIFT;
rBuf[ pos ] += out * rAmp >> FP_SHIFT;
pos++;
inputFrac += step;
inputPos += inputFrac >> FP_SHIFT;
inputFrac &= FP_MASK;
}
}
}
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)
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
3 Band Equaliser
References : Posted by Neil C
Notes :
Simple 3 band equaliser with adjustable low and high frequencies ...
Fairly fast algo, good quality output (seems to be accoustically transparent with all gains set to 1.0)
How to use ...
1. First you need to declare a state for your eq
EQSTATE eq;
2. Now initialise the state (we'll assume your output frequency is 48Khz)
set_3band_state(eq,880,5000,480000);
Your EQ bands are now as follows (approximatley!)
low band = 0Hz to 880Hz
mid band = 880Hz to 5000Hz
high band = 5000Hz to 24000Hz
3. Set the gains to some values ...
eq.lg = 1.5; // Boost bass by 50%
eq.mg = 0.75; // Cut mid by 25%
eq.hg = 1.0; // Leave high band alone
4. You can now EQ some samples
out_sample = do_3band(eq,in_sample)
Have fun and mail me if any problems ... etanza at lycos dot co dot uk
Neil C / Etanza Systems, 2006 :)
Code :
First the header file ....
//---------------------------------------------------------------------------
//
// 3 Band EQ :)
//
// EQ.H - Header file for 3 band EQ
//
// (c) Neil C / Etanza Systems / 2K6
//
// Shouts / Loves / Moans = etanza at lycos dot co dot uk
//
// This work is hereby placed in the public domain for all purposes, including
// use in commercial applications.
//
// The author assumes NO RESPONSIBILITY for any problems caused by the use of
// this software.
//
//----------------------------------------------------------------------------
#ifndef __EQ3BAND__
#define __EQ3BAND__
// ------------
//| Structures |
// ------------
typedef struct
{
// Filter #1 (Low band)
double lf; // Frequency
double f1p0; // Poles ...
double f1p1;
double f1p2;
double f1p3;
// Filter #2 (High band)
double hf; // Frequency
double f2p0; // Poles ...
double f2p1;
double f2p2;
double f2p3;
// Sample history buffer
double sdm1; // Sample data minus 1
double sdm2; // 2
double sdm3; // 3
// Gain Controls
double lg; // low gain
double mg; // mid gain
double hg; // high gain
} EQSTATE;
// ---------
//| Exports |
// ---------
extern void init_3band_state(EQSTATE* es, int lowfreq, int highfreq, int mixfreq);
extern double do_3band(EQSTATE* es, double sample);
#endif // #ifndef __EQ3BAND__
//---------------------------------------------------------------------------
Now the source ...
//----------------------------------------------------------------------------
//
// 3 Band EQ :)
//
// EQ.C - Main Source file for 3 band EQ
//
// (c) Neil C / Etanza Systems / 2K6
//
// Shouts / Loves / Moans = etanza at lycos dot co dot uk
//
// This work is hereby placed in the public domain for all purposes, including
// use in commercial applications.
//
// The author assumes NO RESPONSIBILITY for any problems caused by the use of
// this software.
//
//----------------------------------------------------------------------------
// NOTES :
//
// - Original filter code by Paul Kellet (musicdsp.pdf)
//
// - Uses 4 first order filters in series, should give 24dB per octave
//
// - Now with P4 Denormal fix :)
//----------------------------------------------------------------------------
// ----------
//| Includes |
// ----------
#include <math.h>
#include "eq.h"
// -----------
//| Constants |
// -----------
static double vsa = (1.0 / 4294967295.0); // Very small amount (Denormal Fix)
// ---------------
//| Initialise EQ |
// ---------------
// Recommended frequencies are ...
//
// lowfreq = 880 Hz
// highfreq = 5000 Hz
//
// Set mixfreq to whatever rate your system is using (eg 48Khz)
void init_3band_state(EQSTATE* es, int lowfreq, int highfreq, int mixfreq)
{
// Clear state
memset(es,0,sizeof(EQSTATE));
// Set Low/Mid/High gains to unity
es->lg = 1.0;
es->mg = 1.0;
es->hg = 1.0;
// Calculate filter cutoff frequencies
es->lf = 2 * sin(M_PI * ((double)lowfreq / (double)mixfreq));
es->hf = 2 * sin(M_PI * ((double)highfreq / (double)mixfreq));
}
// ---------------
//| EQ one sample |
// ---------------
// - sample can be any range you like :)
//
// Note that the output will depend on the gain settings for each band
// (especially the bass) so may require clipping before output, but you
// knew that anyway :)
double do_3band(EQSTATE* es, double sample)
{
// Locals
double l,m,h; // Low / Mid / High - Sample Values
// Filter #1 (lowpass)
es->f1p0 += (es->lf * (sample - es->f1p0)) + vsa;
es->f1p1 += (es->lf * (es->f1p0 - es->f1p1));
es->f1p2 += (es->lf * (es->f1p1 - es->f1p2));
es->f1p3 += (es->lf * (es->f1p2 - es->f1p3));
l = es->f1p3;
// Filter #2 (highpass)
es->f2p0 += (es->hf * (sample - es->f2p0)) + vsa;
es->f2p1 += (es->hf * (es->f2p0 - es->f2p1));
es->f2p2 += (es->hf * (es->f2p1 - es->f2p2));
es->f2p3 += (es->hf * (es->f2p2 - es->f2p3));
h = es->sdm3 - es->f2p3;
// Calculate midrange (signal - (low + high))
m = es->sdm3 - (h + l);
// Scale, Combine and store
l *= es->lg;
m *= es->mg;
h *= es->hg;
// Shuffle history buffer
es->sdm3 = es->sdm2;
es->sdm2 = es->sdm1;
es->sdm1 = sample;
// Return result
return(l + m + h);
}
//----------------------------------------------------------------------------
13 comment(s) | 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!!
1 comment(s) | 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 |
|
|
|
|
|
|
|
|
4th order Linkwitz-Riley filters
Type : LP/HP - LR4
References : Posted by neolit123 at gmail dot com
Notes :
Original from T. Lossius - ttblue project
Optimized version in pseudo-code.
[! The filter is unstable for fast automation changes in the lower frequency range. Parameter interpolation and/or oversampling should fix this. !]
The sum of the Linkwitz-Riley (Butterworth squared) HP and LP outputs, will result an all-pass filter at Fc and flat magnitude response - close to ideal for crossovers.
Lubomir I. Ivanov
Code :
//-----------------------------------------
// [code]
//-----------------------------------------
//fc -> cutoff frequency
//pi -> 3.14285714285714
//srate -> sample rate
//================================================
// shared for both lp, hp; optimizations here
//================================================
wc=2*pi*fc;
wc2=wc*wc;
wc3=wc2*wc;
wc4=wc2*wc2;
k=wc/tan(pi*fc/srate);
k2=k*k;
k3=k2*k;
k4=k2*k2;
sqrt2=sqrt(2);
sq_tmp1=sqrt2*wc3*k;
sq_tmp2=sqrt2*wc*k3;
a_tmp=4*wc2*k2+2*sq_tmp1+k4+2*sq_tmp2+wc4;
b1=(4*(wc4+sq_tmp1-k4-sq_tmp2))/a_tmp;
b2=(6*wc4-8*wc2*k2+6*k4)/a_tmp;
b3=(4*(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp;
b4=(k4-2*sq_tmp1+wc4-2*sq_tmp2+4*wc2*k2)/a_tmp;
//================================================
// low-pass
//================================================
a0=wc4/a_tmp;
a1=4*wc4/a_tmp;
a2=6*wc4/a_tmp;
a3=a1;
a4=a0;
//=====================================================
// high-pass
//=====================================================
a0=k4/a_tmp;
a1=-4*k4/a_tmp;
a2=6*k4/a_tmp;
a3=a1;
a4=a0;
//=====================================================
// sample loop - same for lp, hp
//=====================================================
tempx=input;
tempy=a0*tempx+a1*xm1+a2*xm2+a3*xm3+a4*xm4-b1*ym1-b2*ym2-b3*ym3-b4*ym4;
xm4=xm3;
xm3=xm2;
xm2=xm1;
xm1=tempx;
ym4=ym3;
ym3=ym2;
ym2=ym1;
ym1=tempy;
output=tempy;
1 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)))));
};
2 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 |
|
|
|
|
|
|
|
|
AM Formantic Synthesis
References : Posted by Paul Sernine
Notes :
Here is another tutorial from Doc Rochebois.
It performs formantic synthesis without filters and without grains. Instead, it uses "double carrier amplitude modulation" to pitch shift formantic waveforms. Just beware the phase relationships to avoid interferences. Some patches of the DX7 used the same trick but phase interferences were a problem. Here, Thierry Rochebois avoids them by using cosine-phased waveforms.
Various formantic waveforms are precalculated and put in tables, they correspond to different formant widths.
The runtime uses many intances (here 4) of these and pitch shifts them with double carriers (to preserve the harmonicity of the signal).
This is a tutorial code, it can be optimized in many ways.
Have Fun
Paul
Code :
// FormantsAM.cpp
// Thierry Rochebois' "Formantic Synthesis by Double Amplitude Modulation"
// Based on a tutorial by Thierry Rochebois.
// Comments by Paul Sernine.
// The spectral content of the signal is obtained by adding amplitude modulated formantic
// waveforms. The amplitude modulations spectraly shift the formantic waveforms.
// Continuous spectral shift, without losing the harmonic structure, is obtained
// by using crossfaded double carriers (multiple of the base frequency).
// To avoid unwanted interference artifacts, phase relationships must be of the
// "cosine type".
// The output is a 44100Hz 16bit stereo PCM file.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//Approximates cos(pi*x) for x in [-1,1].
inline float fast_cos(const float x)
{
float x2=x*x;
return 1+x2*(-4+2*x2);
}
//Length of the table
#define L_TABLE (256+1) //The last entry of the table equals the first (to avoid a modulo)
//Maximal formant width
#define I_MAX 64
//Table of formants
float TF[L_TABLE*I_MAX];
//Formantic function of width I (used to fill the table of formants)
float fonc_formant(float p,const float I)
{
float a=0.5f;
int hmax=int(10*I)>L_TABLE/2?L_TABLE/2:int(10*I);
float phi=0.0f;
for(int h=1;h<hmax;h++)
{
phi+=3.14159265359f*p;
float hann=0.5f+0.5f*fast_cos(h*(1.0f/hmax));
float gaussienne=0.85f*exp(-h*h/(I*I));
float jupe=0.15f;
float harmonique=cosf(phi);
a+=hann*(gaussienne+jupe)*harmonique;
}
return a;
}
//Initialisation of the table TF with the fonction fonc_formant.
void init_formant(void)
{ float coef=2.0f/(L_TABLE-1);
for(int I=0;I<I_MAX;I++)
for(int P=0;P<L_TABLE;P++)
TF[P+I*L_TABLE]=fonc_formant(-1+P*coef,float(I));
}
//This function emulates the function fonc_formant
// thanks to the table TF. A bilinear interpolation is
// performed
float formant(float p,float i)
{
i=i<0?0:i>I_MAX-2?I_MAX-2:i; // width limitation
float P=(L_TABLE-1)*(p+1)*0.5f; // phase normalisation
int P0=(int)P; float fP=P-P0; // Integer and fractional
int I0=(int)i; float fI=i-I0; // parts of the phase (p) and width (i).
int i00=P0+L_TABLE*I0; int i10=i00+L_TABLE;
//bilinear interpolation.
return (1-fI)*(TF[i00] + fP*(TF[i00+1]-TF[i00]))
+ fI*(TF[i10] + fP*(TF[i10+1]-TF[i10]));
}
// Double carrier.
// h : position (float harmonic number)
// p : phase
float porteuse(const float h,const float p)
{
float h0=floor(h); //integer and
float hf=h-h0; //decimal part of harmonic number.
// modulos pour ramener p*h0 et p*(h0+1) dans [-1,1]
float phi0=fmodf(p* h0 +1+1000,2.0f)-1.0f;
float phi1=fmodf(p*(h0+1)+1+1000,2.0f)-1.0f;
// two carriers.
float Porteuse0=fast_cos(phi0); float Porteuse1=fast_cos(phi1);
// crossfade between the two carriers.
return Porteuse0+hf*(Porteuse1-Porteuse0);
}
int main()
{
//Formant table for various french vowels (you can add your own)
float F1[]={ 730, 200, 400, 250, 190, 350, 550, 550, 450};
float A1[]={ 1.0f, 0.5f, 1.0f, 1.0f, 0.7f, 1.0f, 1.0f, 0.3f, 1.0f};
float F2[]={ 1090, 2100, 900, 1700, 800, 1900, 1600, 850, 1100};
float A2[]={ 2.0f, 0.5f, 0.7f, 0.7f,0.35f, 0.3f, 0.5f, 1.0f, 0.7f};
float F3[]={ 2440, 3100, 2300, 2100, 2000, 2500, 2250, 1900, 1500};
float A3[]={ 0.3f,0.15f, 0.2f, 0.4f, 0.1f, 0.3f, 0.7f, 0.2f, 0.2f};
float F4[]={ 3400, 4700, 3000, 3300, 3400, 3700, 3200, 3000, 3000};
float A4[]={ 0.2f, 0.1f, 0.2f, 0.3f, 0.1f, 0.1f, 0.3f, 0.2f, 0.3f};
float f0,dp0,p0=0.0f;
int F=7; //number of the current formant preset
float f1,f2,f3,f4,a1,a2,a3,a4;
f1=f2=f3=f4=100.0f;a1=a2=a3=a4=0.0f;
init_formant();
FILE *f=fopen("sortie.pcm","wb");
for(int ns=0;ns<10*44100;ns++)
{
if(0==(ns%11025)){F++;F%=8;} //formant change
f0=12*powf(2.0f,4-4*ns/(10*44100.0f)); //sweep
f0*=(1.0f+0.01f*sinf(ns*0.0015f)); //vibrato
dp0=f0*(1/22050.0f);
float un_f0=1.0f/f0;
p0+=dp0; //phase increment
p0-=2*(p0>1);
{ //smoothing of the commands.
float r=0.001f;
f1+=r*(F1[F]-f1);f2+=r*(F2[F]-f2);f3+=r*(F3[F]-f3);f4+=r*(F4[F]-f4);
a1+=r*(A1[F]-a1);a2+=r*(A2[F]-a2);a3+=r*(A3[F]-a3);a4+=r*(A4[F]-a4);
}
//The f0/fn coefficients stand for a -3dB/oct spectral enveloppe
float out=
a1*(f0/f1)*formant(p0,100*un_f0)*porteuse(f1*un_f0,p0)
+0.7f*a2*(f0/f2)*formant(p0,120*un_f0)*porteuse(f2*un_f0,p0)
+ a3*(f0/f3)*formant(p0,150*un_f0)*porteuse(f3*un_f0,p0)
+ a4*(f0/f4)*formant(p0,300*un_f0)*porteuse(f4*un_f0,p0);
short s=short(15000.0f*out);
fwrite(&s,2,1,f);fwrite(&s,2,1,f); //fichier raw pcm stereo
}
fclose(f);
return 0;
}
3 comment(s) | 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;
1 comment(s) | 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));
--------------
4 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.
3 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;
}
}
}
7 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
antialiased square generator
Type : 1st April edition
References : Posted by Paul Sernine
Notes :
It is based on a code by Thierry Rochebois, obfuscated by me.
It generates a 16bit MONO raw pcm file. Have Fun.
Code :
//sqrfish.cpp
#include <math.h>
#include <stdio.h>
//obfuscation P.Sernine
int main() {float ccc,cccc=0,CC=0,cc=0,CCCC,
CCC,C,c; FILE *CCCCCCC=fopen("sqrfish.pcm",
"wb" ); int ccccc= 0; float CCCCC=6.89e-6f;
for(int CCCCCC=0;CCCCCC<1764000;CCCCCC++ ){
if(!(CCCCCC%7350)){if(++ccccc>=30){ ccccc =0;
CCCCC*=2;}CCC=1;}ccc=CCCCC*expf(0.057762265f*
"aiakahiafahadfaiakahiahafahadf"[ccccc]);CCCC
=0.75f-1.5f*ccc;cccc+=ccc;CCC*=0.9999f;cccc-=
2*(cccc>1);C=cccc+CCCC*CC; c=cccc+CCCC*cc; C
-=2*(C>1);c-=2*(c>1);C+=2*(C<-1); c+=1+2
*(c<-1);c-=2*(c>1);C=C*C*(2 *C*C-4);
c=c*c*(2*c*c-4); short cccccc=short(15000.0f*
CCC*(C-c )*CCC);CC=0.5f*(1+C+CC);cc=0.5f*(1+
c+cc); fwrite(&cccccc,2,1,CCCCCCC);}
//algo by Thierry Rochebois
fclose(CCCCCCC);
return 0000000;}
8 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...
5 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
5 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Automatic PDC system
Type : the type that actually works, completely
References : Posted by Tebello Thejane
Linked file : pdc.pdf
Notes :
No, people, implementing PDC is actually not as difficult as you might think it is.
This paper presents a solution to the problem of latency inherent in audio effects processors, and the two appendices give examples of the method being applied on Cubase SX (with an example which its native half-baked PDC fails to solve properly) as well as a convoluted example in FL Studio (taking advantage of the flexible routing capabilities introduced in version 6 of the software). All that's necessary to understand it is a grasp of basic algebra and an intermediate understanding of how music production software works (no need to understand the Laplace transform, linear processes, sigma and integral notation... YAY!).
Please do send me any feedback (kudos, errata, flames, job offers, questions, comments) you might have - my email address is included in the paper - or simply use musicdsp.org's own commenting system.
Tebello Thejane.
Code :
(I have sent the PDF to Bram as he suggested)
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Band Limited PWM Generator
Type : PWM generator
References : Posted by paul_sernine75 AT hotmail DOT fr
Notes :
This is a commented and deobfuscated version of my 1st April fish. It is based on a tutorial code by Thierry Rochebois. I just translated and added comments.
Regards,
Paul Sernine.
Code :
// SelfPMpwm.cpp
// Antialised PWM oscillator
// Based on a tutorial code by Thierry Rochebois (98).
// Itself inspired by US patent 4249447 by Norio Tomisawa (81).
// Comments added/translated by P.Sernine (06).
// This program generates a 44100Hz-raw-PCM-mono-wavefile.
// It is based on Tomisawa's self-phase-modulated sinewave generators.
// Rochebois uses a common phase accumulator to feed two half-Tomisawa-
// oscillators. Each half-Tomisawa-oscillator generates a bandlimited
// sawtooth (band limitation depending on the feedback coeff B).
// These half oscillators are phase offseted according to the desired
// pulse width. They are finally combined to obtain the PW signal.
// Note: the anti-"hunting" filter is a critical feature of a good
// implementation of Tomisawa's method.
#include <math.h>
#include <stdio.h>
const float pi=3.14159265359f;
int main()
{
float freq,dphi; //!< frequency (Hz) and phase increment(rad/sample)
float dphif=0; //!< filtered (anti click) phase increment
float phi=-pi; //!< phase
float Y0=0,Y1=0; //!< feedback memories
float PW=pi; //!< pulse width ]0,2pi[
float B=2.3f; //!< feedback coef
FILE *f=fopen("SelfPMpwm.pcm","wb");
// séquence ('a'=mi=E)
// you can edit this if you prefer another melody.
static char seq[]="aiakahiafahadfaiakahiahafahadf"; //!< sequence
int note=sizeof(seq)-2; //!< note number in the sequence
int octave=0; //!< octave number
float env,envf=0; //!< envelopped and filtered envelopped
for(int ns=0;ns<8*(sizeof(seq)-1)*44100/6;ns++)
{
//waveform control --------------------------------------------------
//Frequency
//freq=27.5f*powf(2.0f,8*ns/(8*30*44100.0f/6)); //sweep
freq=27.5f*powf(2.0f,octave+(seq[note]-'a'-5)/12.0f);
//freq*=(1.0f+0.01f*sinf(ns*0.0015f)); //vibrato
dphi=freq*(pi/22050.0f); // phase increment
dphif+=0.1f*(dphi-dphif);
//notes and enveloppe trigger
if((ns%(44100/6))==0)
{
note++;
if(note>=(sizeof(seq)-1))// sequence loop
{
note=0;
octave++;
}
env=1; //env set
//PW=pi*(0.4+0.5f*(rand()%1000)/1000.0f); //random PW
}
env*=0.9998f; // exp enveloppe
envf+=0.1f*(env-envf); // de-clicked enveloppe
B=1.0f; // feedback coefficient
//try this for a nice bass sound:
//B*=envf*envf; // feedback controlled by enveloppe
B*=2.3f*(1-0.0001f*freq); // feedback limitation
if(B<0)
B=0;
//waveform generation -----------------------------------------------
//Common phase
phi+=dphif; // phase increment
if(phi>=pi)
phi-=2*pi; // phase wrapping
// "phase" half Tomisawa generator 0
// B*Y0 -> self phase modulation
float out0=cosf(phi+B*Y0); // half-output 0
Y0=0.5f*(out0+Y0); // anti "hunting" filter
// "phase+PW" half Tomisawa generator 1
// B*Y1 -> self phase modulation
// PW -> phase offset
float out1=cosf(phi+B*Y1+PW); // half-output 1
Y1=0.5f*(out1+Y1); // anti "hunting" filter
// combination, enveloppe and output
short s=short(15000.0f*(out0-out1)*envf);
fwrite(&s,2,1,f); // file output
}
fclose(f);
return 0;
}
7 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Band Limited waveforms my way
Type : classic Sawtooth example
References : Posted by Anton Savov (antto mail bg)
Notes :
This is my <ugly> C++ code for generating a single cycle of a Sawtooth in a table
normaly i create my "fundamental" table big enough to hold on around 20-40Hz in the current Sampling rate
also, i create the table twice as big, i do "mip-maps" then
so the size should be a power of two, say 1024 for 44100Hz = 44100/1024 = ~43.066Hz
then the mip-maps are with decreasing sizes (twice) 512, 256, 128, 64, 32, 16, 8, 4, and 2
if the "gibbs" effect is what i think it is - then i have a simple solution
here is my crappy code:
Code :
int sz = 1024; // the size of the table
int i = 0;
float *table; // pointer to the table
double scale = 1.0;
double pd; // phase
double omega = 1.0 / (double)(sz);
while (i < sz)
{
double amp = scale;
double x = 0.0; // the sample
double h = 1; // harmonic number (starts from 1)
double dd; // fix high frequency "ring"
pd = (double)(i) / (double)(sz); // calc phase
double hpd = pd; // phase of the harmonic
while (true) // start looping for this sample
{
if ((omega * h) < 0.5) // harmonic frequency is in range?
{
dd = cos(omega * h * 2 * pi);
x = x + (amp * dd * sin(hpd * 2 * pi));
h = h + 1;
hpd = pd * h;
amp = 1.0 / h;
}
else { break; }
}
table[i] = x;
++i;
}
the peaks are around +/- 0.8
a square can be generated by just changing h = h+2; the peaks would be +/- 0.4
any bugs/improvements?
3 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 |
|
|
|
|
|
|
|
|
Bass Booster
Type : LP and SUM
References : Posted by Johny Dupej
Notes :
This function adds a low-passed signal to the original signal. The low-pass has a quite wide response.
Params:
selectivity - frequency response of the LP (higher value gives a steeper one) [70.0 to 140.0 sounds good]
ratio - how much of the filtered signal is mixed to the original
gain2 - adjusts the final volume to handle cut-offs (might be good to set dynamically)
Code :
#define saturate(x) __min(__max(-1.0,x),1.0)
float BassBoosta(float sample)
{
static float selectivity, gain1, gain2, ratio, cap;
gain1 = 1.0/(selectivity + 1.0);
cap= (sample + cap*selectivity )*gain1;
sample = saturate((sample + cap*ratio)*gain2);
return sample;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Beat Detector Class
References : Posted by DSPMaster[at]free[dot]fr
Notes :
This class was designed for a VST plugin. Basically, it's just a 2nd order LP filter, followed by an enveloppe detector (thanks Bram), feeding a Schmitt trigger. The rising edge detector provides a 1-sample pulse each time a beat is detected. Code is self documented...
Note : The class uses a fixed comparison level, you may need to change it.
Code :
// ***** BEATDETECTOR.H *****
#ifndef BeatDetectorH
#define BeatDetectorH
class TBeatDetector
{
private:
float KBeatFilter; // Filter coefficient
float Filter1Out, Filter2Out;
float BeatRelease; // Release time coefficient
float PeakEnv; // Peak enveloppe follower
bool BeatTrigger; // Schmitt trigger output
bool PrevBeatPulse; // Rising edge memory
public:
bool BeatPulse; // Beat detector output
TBeatDetector();
~TBeatDetector();
virtual void setSampleRate(float SampleRate);
virtual void AudioProcess (float input);
};
#endif
// ***** BEATDETECTOR.CPP *****
#include "BeatDetector.h"
#include "math.h"
#define FREQ_LP_BEAT 150.0f // Low Pass filter frequency
#define T_FILTER 1.0f/(2.0f*M_PI*FREQ_LP_BEAT) // Low Pass filter time constant
#define BEAT_RTIME 0.02f // Release time of enveloppe detector in second
TBeatDetector::TBeatDetector()
// Beat detector constructor
{
Filter1Out=0.0;
Filter2Out=0.0;
PeakEnv=0.0;
BeatTrigger=false;
PrevBeatPulse=false;
setSampleRate(44100);
}
TBeatDetector::~TBeatDetector()
{
// Nothing specific to do...
}
void TBeatDetector::setSampleRate (float sampleRate)
// Compute all sample frequency related coeffs
{
KBeatFilter=1.0/(sampleRate*T_FILTER);
BeatRelease=(float)exp(-1.0f/(sampleRate*BEAT_RTIME));
}
void TBeatDetector::AudioProcess (float input)
// Process incoming signal
{
float EnvIn;
// Step 1 : 2nd order low pass filter (made of two 1st order RC filter)
Filter1Out=Filter1Out+(KBeatFilter*(input-Filter1Out));
Filter2Out=Filter2Out+(KBeatFilter*(Filter1Out-Filter2Out));
// Step 2 : peak detector
EnvIn=fabs(Filter2Out);
if (EnvIn>PeakEnv) PeakEnv=EnvIn; // Attack time = 0
else
{
PeakEnv*=BeatRelease;
PeakEnv+=(1.0f-BeatRelease)*EnvIn;
}
// Step 3 : Schmitt trigger
if (!BeatTrigger)
{
if (PeakEnv>0.3) BeatTrigger=true;
}
else
{
if (PeakEnv<0.15) BeatTrigger=false;
}
// Step 4 : rising edge detector
BeatPulse=false;
if ((BeatTrigger)&&(!PrevBeatPulse))
BeatPulse=true;
PrevBeatPulse=BeatTrigger;
}
5 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Biquad, Butterworth, Chebyshev N-order, M-channel optimized filters
Type : LP, HP, BP, BS, Shelf, Notch, Boost
References : Posted by thevinn at yahoo dot com
Code :
/*
"A Collection of Useful C++ Classes for Digital Signal Processing"
By Vincent Falco
Please direct all comments to either the music-dsp mailing list or
the DSP and Plug-in Development forum:
http://music.columbia.edu/cmc/music-dsp/
http://www.kvraudio.com/forum/viewforum.php?f=33
http://www.kvraudio.com/forum/
Support is provided for performing N-order Dsp floating point filter
operations on M-channel data with a caller specified floating point type.
The implementation breaks a high order IIR filter down into a series of
cascaded second order stages. Tests conclude that numerical stability is
maintained even at higher orders. For example the Butterworth low pass
filter is stable at up to 53 poles.
Processing functions are provided to use either Direct Form I or Direct
Form II of the filter transfer function. Direct Form II is slightly faster
but can cause discontinuities in the output if filter parameters are changed
during processing. Direct Form I is slightly slower, but maintains fidelity
even when parameters are changed during processing.
To support fast parameter changes, filters provide two functions for
adjusting parameters. A high accuracy Setup() function, and a faster
form called SetupFast() that uses approximations for trigonometric
functions. The approximations work quite well and should be suitable for
most applications.
Channels are stored in an interleaved format with M samples per frame
arranged contiguously. A single class instance can process all M channels
simultaneously in an efficient manner. A 'skip' parameter causes the
processing function to advance by skip additional samples in the destination
buffer in between every frame. Through manipulation of the skip paramter it
is possible to exclude channels from processing (for example, only processing
the left half of stereo interleaved data). For multichannel data which is
not interleaved, it will be necessary to instantiate multiple instance of
the filter and set skip=0.
There are a few other utility classes and functions included that may prove useful.
Classes:
Complex
CascadeStages
Biquad
BiquadLowPass
BiquadHighPass
BiquadBandPass1
BiquadBandPass2
BiquadBandStop
BiquadAllPass
BiquadPeakEq
BiquadLowShelf
BiquadHighShelf
PoleFilter
Butterworth
ButterLowPass
ButterHighPass
ButterBandPass
ButterBandStop
Chebyshev1
Cheby1LowPass
Cheby1HighPass
Cheby1BandPass
Cheby1BandStop
Chebyshev2
Cheby2LowPass
Cheby2HighPass
Cheby2BandPass
Cheby2BandStop
EnvelopeFollower
AutoLimiter
Functions:
zero()
copy()
mix()
scale()
interleave()
deinterleave()
Order for PoleFilter derived classes is specified in the number of poles,
except for band pass and band stop filters, for which the number of pole pairs
is specified.
For some filters there are two versions of Setup(), the one called
SetupFast() uses approximations to trigonometric functions for speed.
This is an option if you are doing frequent parameter changes to the filter.
There is an example function at the bottom that shows how to use the classes.
Filter ideas are based on a java applet (http://www.falstad.com/dfilter/)
developed by Paul Falstad.
All of this code was written by the author Vincent Falco except where marked.
--------------------------------------------------------------------------------
License: MIT License (http://www.opensource.org/licenses/mit-license.php)
Copyright (c) 2009 by Vincent Falco
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*
*/
/*
To Do:
- Shelving, peak, all-pass for Butterworth, Chebyshev, and Elliptic.
The Biquads have these versions and I would like the others to have them
as well. It would also be super awesome if higher order filters could
have a "Q" parameter for resonance but I'm not expecting miracles, it
would require redistributing the poles and zeroes. But if there's
a research paper or code out there...I could incorporate it.
- Denormal testing and fixing
I'd like to know if denormals are a problem. And if so, it would be nice
to have a small function that can reproduce the denormal problem. This
way I can test the fix under all conditions. I will include the function
as a "unit test" object in the header file so anyone can verify its
correctness. But I'm a little lost.
- Optimized template specializations for stages=1, channels={1,2}
There are some pretty obvious optimizations I am saving for "the end".
I don't want to do them until the code is finalized.
- Optimized template specializations for SSE, other special instructions
- Optimized trigonometric functions for fast parameter changes
- Elliptic curve based filter coefficients
- More utility functions for manipulating sample buffers
- Need fast version of pow( 10, x )
*/
#ifndef __DSP_FILTER__
#define __DSP_FILTER__
#include <cmath>
#include <cfloat>
#include <assert.h>
#include <memory.h>
#include <stdlib.h>
//#define DSP_USE_STD_COMPLEX
#ifdef DSP_USE_STD_COMPLEX
#include <complex>
#endif
#define DSP_SSE3_OPTIMIZED
#ifdef DSP_SSE3_OPTIMIZED
//#include <xmmintrin.h>
//#include <emmintrin.h>
#include <pmmintrin.h>
#endif
namespace Dsp
{
//--------------------------------------------------------------------------
// WARNING: Here there be templates
//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
//
// Configuration
//
//--------------------------------------------------------------------------
// Regardless of the type of sample that the filter operates on (e.g.
// float or double), all calculations are performed using double (or
// better) for stability and accuracy. This controls the underlying
// type used for calculations:
typedef double CalcT;
typedef int Int32; // Must be 32 bits
typedef __int64 Int64; // Must be 64 bits
// This is used to prevent denormalization.
const CalcT vsa=1.0 / 4294967295.0; // for CalcT as float
// These constants are so important, I made my own copy. If you improve
// the resolution of CalcT be sure to add more significant digits to these.
const CalcT kPi =3.1415926535897932384626433832795;
const CalcT kPi_2 =1.57079632679489661923;
const CalcT kLn2 =0.693147180559945309417;
const CalcT kLn10 =2.30258509299404568402;
//--------------------------------------------------------------------------
// Some functions missing from <math.h>
template<typename T>
inline T acosh( T x )
{
return log( x+::sqrt(x*x-1) );
}
//--------------------------------------------------------------------------
//
// Fast Trigonometric Functions
//
//--------------------------------------------------------------------------
// Three approximations for both sine and cosine at a given angle.
// The faster the routine, the larger the error.
// From http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/
// Tuned for maximum pole stability. r must be in the range 0..pi
// This one breaks down considerably at the higher angles. It is
// included only for educational purposes.
inline void fastestsincos( CalcT r, CalcT *sn, CalcT *cs )
{
const CalcT c=0.70710678118654752440; // sqrt(2)/2
CalcT v=(2-4*c)*r*r+c;
if(r<kPi_2)
{
*sn=v+r; *cs=v-r;
}
else
{
*sn=r+v; *cs=r-v;
}
}
// Lower precision than ::fastsincos() but still decent
inline void fastersincos( CalcT x, CalcT *sn, CalcT *cs )
{
//always wrap input angle to -PI..PI
if (x < -kPi) x += 2*kPi;
else if (x > kPi) x -= 2*kPi;
//compute sine
if (x < 0) *sn = 1.27323954 * x + 0.405284735 * x * x;
else *sn = 1.27323954 * x - 0.405284735 * x * x;
//compute cosine: sin(x + PI/2) = cos(x)
x += kPi_2;
if (x > kPi ) x -= 2*kPi;
if (x < 0) *cs = 1.27323954 * x + 0.405284735 * x * x;
else *cs = 1.27323954 * x - 0.405284735 * x * x;
}
// Slower than ::fastersincos() but still faster than
// sin(), and with the best accuracy of these routines.
inline void fastsincos( CalcT x, CalcT *sn, CalcT *cs )
{
CalcT s, c;
//always wrap input angle to -PI..PI
if (x < -kPi) x += 2*kPi;
else if (x > kPi) x -= 2*kPi;
//compute sine
if (x < 0)
{
s = 1.27323954 * x + .405284735 * x * x;
if (s < 0) s = .225 * (s * -s - s) + s;
else s = .225 * (s * s - s) + s;
}
else
{
s = 1.27323954 * x - 0.405284735 * x * x;
if (s < 0) s = .225 * (s * -s - s) + s;
else s = .225 * (s * s - s) + s;
}
*sn=s;
//compute cosine: sin(x + PI/2) = cos(x)
x += kPi_2;
if (x > kPi ) x -= 2*kPi;
if (x < 0)
{
c = 1.27323954 * x + 0.405284735 * x * x;
if (c < 0) c = .225 * (c * -c - c) + c;
else c = .225 * (c * c - c) + c;
}
else
{
c = 1.27323954 * x - 0.405284735 * x * x;
if (c < 0) c = .225 * (c * -c - c) + c;
else c = .225 * (c * c - c) + c;
}
*cs=c;
}
// Faster approximations to sqrt()
// From http://ilab.usc.edu/wiki/index.php/Fast_Square_Root
// The faster the routine, the more error in the approximation.
// Log Base 2 Approximation
// 5 times faster than sqrt()
inline float fastsqrt1( float x )
{
union { Int32 i; float x; } u;
u.x = x;
u.i = (Int32(1)<<29) + (u.i >> 1) - (Int32(1)<<22);
return u.x;
}
inline double fastsqrt1( double x )
{
union { Int64 i; double x; } u;
u.x = x;
u.i = (Int64(1)<<61) + (u.i >> 1) - (Int64(1)<<51);
return u.x;
}
// Log Base 2 Approximation with one extra Babylonian Step
// 2 times faster than sqrt()
inline float fastsqrt2( float x )
{
float v=fastsqrt1( x );
v = 0.5f * (v + x/v); // One Babylonian step
return v;
}
inline double fastsqrt2(const double x)
{
double v=fastsqrt1( x );
v = 0.5f * (v + x/v); // One Babylonian step
return v;
}
// Log Base 2 Approximation with two extra Babylonian Steps
// 50% faster than sqrt()
inline float fastsqrt3( float x )
{
float v=fastsqrt1( x );
v = v + x/v;
v = 0.25f* v + x/v; // Two Babylonian steps
return v;
}
inline double fastsqrt3(const double x)
{
double v=fastsqrt1( x );
v = v + x/v;
v = 0.25 * v + x/v; // Two Babylonian steps
return v;
}
//--------------------------------------------------------------------------
//
// Complex
//
//--------------------------------------------------------------------------
#ifdef DSP_USE_STD_COMPLEX
template<typename T>
inline std::complex<T> polar( const T &m, const T &a )
{
return std::polar( m, a );
}
template<typename T>
inline T norm( const std::complex<T> &c )
{
return std::norm( c );
}
template<typename T>
inline T abs( const std::complex<T> &c )
{
return std::abs(c);
}
template<typename T, typename To>
inline std::complex<T> addmul( const std::complex<T> &c, T v, const std::complex<To> &c1 )
{
return std::complex<T>( c.real()+v*c1.real(), c.imag()+v*c1.imag() );
}
template<typename T>
inline T arg( const std::complex<T> &c )
{
return std::arg( c );
}
template<typename T>
inline std::complex<T> recip( const std::complex<T> &c )
{
T n=1.0/Dsp::norm(c);
return std::complex<T>( n*c.real(), n*c.imag() );
}
template<typename T>
inline std::complex<T> sqrt( const std::complex<T> &c )
{
return std::sqrt( c );
}
typedef std::complex<CalcT> Complex;
#else
//--------------------------------------------------------------------------
// "Its always good to have a few extra wheels in case one goes flat."
template<typename T>
struct ComplexT
{
ComplexT();
ComplexT( T r_, T i_=0 );
template<typename To>
ComplexT( const ComplexT<To> &c );
T imag ( void ) const;
T real ( void ) const;
ComplexT & neg ( void );
ComplexT & conj ( void );
template<typename To>
ComplexT & add ( const ComplexT<To> &c );
template<typename To>
ComplexT & sub ( const ComplexT<To> &c );
template<typename To>
ComplexT & mul ( const ComplexT<To> &c );
template<typename To>
ComplexT & div ( const ComplexT<To> &c );
template<typename To>
ComplexT & addmul ( T v, const ComplexT<To> &c );
ComplexT operator- ( void ) const;
ComplexT operator+ ( T v ) const;
ComplexT operator- ( T v ) const;
ComplexT operator* ( T v ) const;
ComplexT operator/ ( T v ) const;
ComplexT & operator+= ( T v );
ComplexT & operator-= ( T v );
ComplexT & operator*= ( T v );
ComplexT & operator/= ( T v );
template<typename To>
ComplexT operator+ ( const ComplexT<To> &c ) const;
template<typename To>
ComplexT operator- ( const ComplexT<To> &c ) const;
template<typename To>
ComplexT operator* ( const ComplexT<To> &c ) const;
template<typename To>
ComplexT operator/ ( const ComplexT<To> &c ) const;
template<typename To>
ComplexT & operator+= ( const ComplexT<To> &c );
template<typename To>
ComplexT & operator-= ( const ComplexT<To> &c );
template<typename To>
ComplexT & operator*= ( const ComplexT<To> &c );
template<typename To>
ComplexT & operator/= ( const ComplexT<To> &c );
private:
ComplexT & add ( T v );
ComplexT & sub ( T v );
ComplexT & mul ( T c, T d );
ComplexT & mul ( T v );
ComplexT & div ( T v );
T r;
T i;
};
//--------------------------------------------------------------------------
template<typename T>
inline ComplexT<T>::ComplexT()
{
}
template<typename T>
inline ComplexT<T>::ComplexT( T r_, T i_ )
{
r=r_;
i=i_;
}
template<typename T>
template<typename To>
inline ComplexT<T>::ComplexT( const ComplexT<To> &c )
{
r=c.r;
i=c.i;
}
template<typename T>
inline T ComplexT<T>::imag( void ) const
{
return i;
}
template<typename T>
inline T ComplexT<T>::real( void ) const
{
return r;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::neg( void )
{
r=-r;
i=-i;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::conj( void )
{
i=-i;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::add( T v )
{
r+=v;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::sub( T v )
{
r-=v;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::mul( T c, T d )
{
T ac=r*c;
T bd=i*d;
// must do i first
i=(r+i)*(c+d)-(ac+bd);
r=ac-bd;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::mul( T v )
{
r*=v;
i*=v;
return *this;
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::div( T v )
{
r/=v;
i/=v;
return *this;
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::add( const ComplexT<To> &c )
{
r+=c.r;
i+=c.i;
return *this;
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::sub( const ComplexT<To> &c )
{
r-=c.r;
i-=c.i;
return *this;
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::mul( const ComplexT<To> &c )
{
return mul( c.r, c.i );
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::div( const ComplexT<To> &c )
{
T s=1.0/norm(c);
return mul( c.r*s, -c.i*s );
}
template<typename T>
inline ComplexT<T> ComplexT<T>::operator-( void ) const
{
return ComplexT<T>(*this).neg();
}
template<typename T>
inline ComplexT<T> ComplexT<T>::operator+( T v ) const
{
return ComplexT<T>(*this).add( v );
}
template<typename T>
inline ComplexT<T> ComplexT<T>::operator-( T v ) const
{
return ComplexT<T>(*this).sub( v );
}
template<typename T>
inline ComplexT<T> ComplexT<T>::operator*( T v ) const
{
return ComplexT<T>(*this).mul( v );
}
template<typename T>
inline ComplexT<T> ComplexT<T>::operator/( T v ) const
{
return ComplexT<T>(*this).div( v );
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::operator+=( T v )
{
return add( v );
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::operator-=( T v )
{
return sub( v );
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::operator*=( T v )
{
return mul( v );
}
template<typename T>
inline ComplexT<T> &ComplexT<T>::operator/=( T v )
{
return div( v );
}
template<typename T>
template<typename To>
inline ComplexT<T> ComplexT<T>::operator+( const ComplexT<To> &c ) const
{
return ComplexT<T>(*this).add(c);
}
template<typename T>
template<typename To>
inline ComplexT<T> ComplexT<T>::operator-( const ComplexT<To> &c ) const
{
return ComplexT<T>(*this).sub(c);
}
template<typename T>
template<typename To>
inline ComplexT<T> ComplexT<T>::operator*( const ComplexT<To> &c ) const
{
return ComplexT<T>(*this).mul(c);
}
template<typename T>
template<typename To>
inline ComplexT<T> ComplexT<T>::operator/( const ComplexT<To> &c ) const
{
return ComplexT<T>(*this).div(c);
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::operator+=( const ComplexT<To> &c )
{
return add( c );
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::operator-=( const ComplexT<To> &c )
{
return sub( c );
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::operator*=( const ComplexT<To> &c )
{
return mul( c );
}
template<typename T>
template<typename To>
inline ComplexT<T> &ComplexT<T>::operator/=( const ComplexT<To> &c )
{
return div( c );
}
//--------------------------------------------------------------------------
template<typename T>
inline ComplexT<T> polar( const T &m, const T &a )
{
return ComplexT<T>( m*cos(a), m*sin(a) );
}
template<typename T>
inline T norm( const ComplexT<T> &c )
{
return c.real()*c.real()+c.imag()*c.imag();
}
template<typename T>
inline T abs( const ComplexT<T> &c )
{
return ::sqrt( c.real()*c.real()+c.imag()*c.imag() );
}
template<typename T, typename To>
inline ComplexT<T> addmul( const ComplexT<T> &c, T v, const ComplexT<To> &c1 )
{
return ComplexT<T>( c.real()+v*c1.real(), c.imag()+v*c1.imag() );
}
template<typename T>
inline T arg( const ComplexT<T> &c )
{
return atan2( c.imag(), c.real() );
}
template<typename T>
inline ComplexT<T> recip( const ComplexT<T> &c )
{
T n=1.0/norm(c);
return ComplexT<T>( n*c.real(), -n*c.imag() );
}
template<typename T>
inline ComplexT<T> sqrt( const ComplexT<T> &c )
{
return polar( ::sqrt(abs(c)), arg(c)*0.5 );
}
//--------------------------------------------------------------------------
typedef ComplexT<CalcT> Complex;
#endif
//--------------------------------------------------------------------------
//
// Numerical Analysis
//
//--------------------------------------------------------------------------
// Implementation of Brent's Method provided by
// John D. Cook (http://www.johndcook.com/)
// The return value of Minimize is the minimum of the function f.
// The location where f takes its minimum is returned in the variable minLoc.
// Notation and implementation based on Chapter 5 of Richard Brent's book
// "Algorithms for Minimization Without Derivatives".
template<class TFunction>
CalcT BrentMinimize
(
TFunction& f, // [in] objective function to minimize
CalcT leftEnd, // [in] smaller value of bracketing interval
CalcT rightEnd, // [in] larger value of bracketing interval
CalcT epsilon, // [in] stopping tolerance
CalcT& minLoc // [out] location of minimum
)
{
CalcT d, e, m, p, q, r, tol, t2, u, v, w, fu, fv, fw, fx;
static const CalcT c = 0.5*(3.0 - ::sqrt(5.0));
static const CalcT SQRT_DBL_EPSILON = ::sqrt(DBL_EPSILON);
CalcT& a = leftEnd; CalcT& b = rightEnd; CalcT& x = minLoc;
v = w = x = a + c*(b - a); d = e = 0.0;
fv = fw = fx = f(x);
int counter = 0;
loop:
counter++;
m = 0.5*(a + b);
tol = SQRT_DBL_EPSILON*::fabs(x) + epsilon; t2 = 2.0*tol;
// Check stopping criteria
if (::fabs(x - m) > t2 - 0.5*(b - a))
{
p = q = r = 0.0;
if (::fabs(e) > tol)
{
// fit parabola
r = (x - w)*(fx - fv);
q = (x - v)*(fx - fw);
p = (x - v)*q - (x - w)*r;
q = 2.0*(q - r);
(q > 0.0) ? p = -p : q = -q;
r = e; e = d;
}
if (::fabs(p) < ::fabs(0.5*q*r) && p < q*(a - x) && p < q*(b - x))
{
// A parabolic interpolation step
d = p/q;
u = x + d;
// f must not be evaluated too close to a or b
if (u - a < t2 || b - u < t2)
d = (x < m) ? tol : -tol;
}
else
{
// A golden section step
e = (x < m) ? b : a;
e -= x;
d = c*e;
}
// f must not be evaluated too close to x
if (::fabs(d) >= tol)
u = x + d;
else if (d > 0.0)
u = x + tol;
else
u = x - tol;
fu = f(u);
// Update a, b, v, w, and x
if (fu <= fx)
{
(u < x) ? b = x : a = x;
v = w; fv = fw;
w = x; fw = fx;
x = u; fx = fu;
}
else
{
(u < x) ? a = u : b = u;
if (fu <= fw || w == x)
{
v = w; fv = fw;
w = u; fw = fu;
}
else if (fu <= fv || v == x || v == w)
{
v = u; fv = fu;
}
}
goto loop; // Yes, the dreaded goto statement. But the code
// here is faithful to Brent's orginal pseudocode.
}
return fx;
}
//--------------------------------------------------------------------------
//
// Infinite Impulse Response Filters
//
//--------------------------------------------------------------------------
// IIR filter implementation using multiple second-order stages.
class CascadeFilter
{
public:
// Process data in place using Direct Form I
// skip is added after each frame.
// Direct Form I is more suitable when the filter parameters
// are changed often. However, it is slightly slower.
template<typename T>
void ProcessI( size_t frames, T *dest, int skip=0 );
// Process data in place using Direct Form II
// skip is added after each frame.
// Direct Form II is slightly faster than Direct Form I,
// but changing filter parameters on stream can result
// in discontinuities in the output. It is best suited
// for a filter whose parameters are set only once.
template<typename T>
void ProcessII( size_t frames, T *dest, int skip=0 );
// Convenience function that just calls ProcessI.
// Feel free to change the implementation.
template<typename T>
void Process( size_t frames, T *dest, int skip=0 );
// Determine response at angular frequency (0<=w<=kPi)
Complex Response( CalcT w );
// Clear the history buffer.
void Clear( void );
protected:
struct Hist;
struct Stage;
// for m_nchan==2
#ifdef DSP_SSE3_OPTIMIZED
template<typename T>
void ProcessISSEStageStereo( size_t frames, T *dest, Stage *stage, Hist *h, int skip );
template<typename T>
void ProcessISSEStereo( size_t frames, T *dest, int skip );
#endif
protected:
void Reset ( void );
void Normalize ( CalcT scale );
void SetAStage ( CalcT x1, CalcT x2 );
void SetBStage ( CalcT x0, CalcT x1, CalcT x2 );
void SetStage ( CalcT a1, CalcT a2, CalcT b0, CalcT b1, CalcT b2 );
protected:
struct Hist
{
CalcT v[4];
};
struct Stage
{
CalcT a[3]; // a[0] unused
CalcT b[3];
void Reset( void );
};
struct ResponseFunctor
{
CascadeFilter *m_filter;
CalcT operator()( CalcT w );
ResponseFunctor( CascadeFilter *filter );
};
int m_nchan;
int m_nstage;
Stage * m_stagep;
Hist * m_histp;
};
//--------------------------------------------------------------------------
template<typename T>
void CascadeFilter::ProcessI( size_t frames, T *dest, int skip )
{
#ifdef DSP_SSE3_OPTIMIZED
if( m_nchan==2 )
ProcessISSEStereo( frames, dest, skip );
else
#endif
while( frames-- )
{
Hist *h=m_histp;
for( int j=m_nchan;j;j-- )
{
CalcT in=CalcT(*dest);
Stage *s=m_stagep;
for( int i=m_nstage;i;i--,h++,s++ )
{
CalcT out;
out=s->b[0]*in + s->b[1]*h->v[0] + s->b[2]*h->v[1] +
s->a[1]*h->v[2] + s->a[2]*h->v[3];
h->v[1]=h->v[0]; h->v[0]=in;
h->v[3]=h->v[2]; h->v[2]=out;
in=out;
}
*dest++=T(in);
}
dest+=skip;
}
}
// A good compiler already produces code that is optimized even for
// the general case. The only way to make it go faster is to
// to implement it in assembler or special instructions. Like this:
#ifdef DSP_SSE3_OPTIMIZED
// ALL SSE OPTIMIZATIONS ASSUME CalcT as double
template<typename T>
inline void CascadeFilter::ProcessISSEStageStereo(
size_t frames, T *dest, Stage *s, Hist *h, int skip )
{
assert( m_nchan==2 );
#if 1
CalcT b0=s->b[0];
__m128d m0=_mm_loadu_pd( &s->a[1] ); // a1 , a2
__m128d m1=_mm_loadu_pd( &s->b[1] ); // b1 , b2
__m128d m2=_mm_loadu_pd( &h[0].v[0] ); // h->v[0] , h->v[1]
__m128d m3=_mm_loadu_pd( &h[0].v[2] ); // h->v[2] , h->v[3]
__m128d m4=_mm_loadu_pd( &h[1].v[0] ); // h->v[0] , h->v[1]
__m128d m5=_mm_loadu_pd( &h[1].v[2] ); // h->v[2] , h->v[3]
while( frames-- )
{
CalcT in, b0in, out;
__m128d m6;
__m128d m7;
in=CalcT(*dest);
b0in=b0*in;
m6=_mm_mul_pd ( m1, m2 ); // b1*h->v[0] , b2*h->v[1]
m7=_mm_mul_pd ( m0, m3 ); // a1*h->v[2] , a2*h->v[3]
m6=_mm_add_pd ( m6, m7 ); // b1*h->v[0] + a1*h->v[2], b2*h->v[1] + a2*h->v[3]
m7=_mm_load_sd( &b0in ); // b0*in , 0
m6=_mm_add_sd ( m6, m7 ); // b1*h->v[0] + a1*h->v[2] + in*b0 , b2*h->v[1] + a2*h->v[3] + 0
m6=_mm_hadd_pd( m6, m7 ); // b1*h->v[0] + a1*h->v[2] + in*b0 + b2*h->v[1] + a2*h->v[3], in*b0
_mm_store_sd( &out, m6 );
m6=_mm_loadh_pd( m6, &in ); // out , in
m2=_mm_shuffle_pd( m6, m2, _MM_SHUFFLE2( 0, 1 ) ); // h->v[0]=in , h->v[1]=h->v[0]
m3=_mm_shuffle_pd( m6, m3, _MM_SHUFFLE2( 0, 0 ) ); // h->v[2]=out, h->v[3]=h->v[2]
*dest++=T(out);
in=CalcT(*dest);
b0in=b0*in;
m6=_mm_mul_pd ( m1, m4 ); // b1*h->v[0] , b2*h->v[1]
m7=_mm_mul_pd ( m0, m5 ); // a1*h->v[2] , a2*h->v[3]
m6=_mm_add_pd ( m6, m7 ); // b1*h->v[0] + a1*h->v[2], b2*h->v[1] + a2*h->v[3]
m7=_mm_load_sd( &b0in ); // b0*in , 0
m6=_mm_add_sd ( m6, m7 ); // b1*h->v[0] + a1*h->v[2] + in*b0 , b2*h->v[1] + a2*h->v[3] + 0
m6=_mm_hadd_pd( m6, m7 ); // b1*h->v[0] + a1*h->v[2] + in*b0 + b2*h->v[1] + a2*h->v[3], in*b0
_mm_store_sd( &out, m6 );
m6=_mm_loadh_pd( m6, &in ); // out , in
m4=_mm_shuffle_pd( m6, m4, _MM_SHUFFLE2( 0, 1 ) ); // h->v[0]=in , h->v[1]=h->v[0]
m5=_mm_shuffle_pd( m6, m5, _MM_SHUFFLE2( 0, 0 ) ); // h->v[2]=out, h->v[3]=h->v[2]
*dest++=T(out);
dest+=skip;
}
// move history from registers back to state
_mm_storeu_pd( &h[0].v[0], m2 );
_mm_storeu_pd( &h[0].v[2], m3 );
_mm_storeu_pd( &h[1].v[0], m4 );
_mm_storeu_pd( &h[1].v[2], m5 );
#else
// Template-specialized version from which the assembly was modeled
CalcT a1=s->a[1];
CalcT a2=s->a[2];
CalcT b0=s->b[0];
CalcT b1=s->b[1];
CalcT b2=s->b[2];
while( frames-- )
{
CalcT in, out;
in=CalcT(*dest);
out=b0*in+b1*h[0].v[0]+b2*h[0].v[1] +a1*h[0].v[2]+a2*h[0].v[3];
h[0].v[1]=h[0].v[0]; h[0].v[0]=in;
h[0].v[3]=h[0].v[2]; h[0].v[2]=out;
in=out;
*dest++=T(in);
in=CalcT(*dest);
out=b0*in+b1*h[1].v[0]+b2*h[1].v[1] +a1*h[1].v[2]+a2*h[1].v[3];
h[1].v[1]=h[1].v[0]; h[1].v[0]=in;
h[1].v[3]=h[1].v[2]; h[1].v[2]=out;
in=out;
*dest++=T(in);
dest+=skip;
}
#endif
}
// Note there could be a loss of accuracy here. Unlike the original version
// of Process...() we are applying each stage to all of the input data.
// Since the underlying type T could be float, the results from this function
// may be different than the unoptimized version. However, it is much faster.
template<typename T>
void CascadeFilter::ProcessISSEStereo( size_t frames, T *dest, int skip )
{
assert( m_nchan==2 );
Stage *s=m_stagep;
Hist *h=m_histp;
for( int i=m_nstage;i;i--,h+=2,s++ )
{
ProcessISSEStageStereo( frames, dest, s, h, skip );
}
}
#endif
template<typename T>
void CascadeFilter::ProcessII( size_t frames, T *dest, int skip )
{
while( frames-- )
{
Hist *h=m_histp;
for( int j=m_nchan;j;j-- )
{
CalcT in=CalcT(*dest);
Stage *s=m_stagep;
for( int i=m_nstage;i;i--,h++,s++ )
{
CalcT d2=h->v[2]=h->v[1];
CalcT d1=h->v[1]=h->v[0];
CalcT d0=h->v[0]=
in+s->a[1]*d1 + s->a[2]*d2;
in=s->b[0]*d0 + s->b[1]*d1 + s->b[2]*d2;
}
*dest++=T(in);
}
dest+=skip;
}
}
template<typename T>
inline void CascadeFilter::Process( size_t frames, T *dest, int skip )
{
ProcessI( frames, dest, skip );
}
inline Complex CascadeFilter::Response( CalcT w )
{
Complex ch( 1 );
Complex cbot( 1 );
Complex czn1=polar( 1., -w );
Complex czn2=polar( 1., -2*w );
Stage *s=m_stagep;
for( int i=m_nstage;i;i-- )
{
Complex ct( s->b[0] );
Complex cb( 1 );
ct=addmul( ct, s->b[1], czn1 );
cb=addmul( cb, -s->a[1], czn1 );
ct=addmul( ct, s->b[2], czn2 );
cb=addmul( cb, -s->a[2], czn2 );
ch*=ct;
cbot*=cb;
s++;
}
return ch/cbot;
}
inline void CascadeFilter::Clear( void )
{
::memset( m_histp, 0, m_nstage*m_nchan*sizeof(m_histp[0]) );
}
inline void CascadeFilter::Stage::Reset( void )
{
a[1]=0; a[2]=0;
b[0]=1; b[1]=0; b[2]=0;
}
inline void CascadeFilter::Reset( void )
{
Stage *s=m_stagep;
for( int i=m_nstage;i;i--,s++ )
s->Reset();
}
// Apply scale factor to stage coefficients.
inline void CascadeFilter::Normalize( CalcT scale )
{
// We are throwing the normalization into the first
// stage. In theory it might be nice to spread it around
// to preserve numerical accuracy.
Stage *s=m_stagep;
s->b[0]*=scale; s->b[1]*=scale; s->b[2]*=scale;
}
inline void CascadeFilter::SetAStage( CalcT x1, CalcT x2 )
{
Stage *s=m_stagep;
for( int i=m_nstage;i;i-- )
{
if( s->a[1]==0 && s->a[2]==0 )
{
s->a[1]=x1;
s->a[2]=x2;
s=0;
break;
}
if( s->a[2]==0 && x2==0 )
{
s->a[2]=-s->a[1]*x1;
s->a[1]+=x1;
s=0;
break;
}
s++;
}
assert( s==0 );
}
inline void CascadeFilter::SetBStage( CalcT x0, CalcT x1, CalcT x2 )
{
Stage *s=m_stagep;
for( int i=m_nstage;i;i-- )
{
if( s->b[1]==0 && s->b[2]==0 )
{
s->b[0]=x0;
s->b[1]=x1;
s->b[2]=x2;
s=0;
break;
}
if( s->b[2]==0 && x2==0 )
{
// (b0 + z b1)(x0 + z x1) = (b0 x0 + (b1 x0+b0 x1) z + b1 x1 z^2)
s->b[2]=s->b[1]*x1;
s->b[1]=s->b[1]*x0+s->b[0]*x1;
s->b[0]*=x0;
s=0;
break;
}
s++;
}
assert( s==0 );
}
// Optimized version for Biquads
inline void CascadeFilter::SetStage(
CalcT a1, CalcT a2, CalcT b0, CalcT b1, CalcT b2 )
{
assert( m_nstage==1 );
Stage *s=&m_stagep[0];
s->a[1]=a1; s->a[2]=a2;
s->b[0]=b0; s->b[1]=b1; s->b[2]=b2;
}
inline CalcT CascadeFilter::ResponseFunctor::operator()( CalcT w )
{
return -Dsp::abs(m_filter->Response( w ));
}
inline CascadeFilter::ResponseFunctor::ResponseFunctor( CascadeFilter *filter )
{
m_filter=filter;
}
//--------------------------------------------------------------------------
template<int stages, int channels>
class CascadeStages : public CascadeFilter
{
public:
CascadeStages();
private:
Hist m_hist [stages*channels];
Stage m_stage [stages];
};
//--------------------------------------------------------------------------
template<int stages, int channels>
CascadeStages<stages, channels>::CascadeStages( void )
{
m_nchan=channels;
m_nstage=stages;
m_stagep=m_stage;
m_histp=m_hist;
Clear();
}
//--------------------------------------------------------------------------
//
// Biquad Second Order IIR Filters
//
//--------------------------------------------------------------------------
// Formulas from http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
template<int channels>
class Biquad : public CascadeStages<1, channels>
{
protected:
void Setup( const CalcT a[3], const CalcT b[3] );
};
//--------------------------------------------------------------------------
template<int channels>
inline void Biquad<channels>::Setup( const CalcT a[3], const CalcT b[3] )
{
Reset();
// transform Biquad coefficients
CalcT ra0=1/a[0];
SetAStage( -a[1]*ra0, -a[2]*ra0 );
SetBStage( b[0]*ra0, b[1]*ra0, b[2]*ra0 );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadLowPass : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadLowPass<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT a0 = 1 / ( 1 + alph );
CalcT b1 = 1 - cs;
CalcT b0 = a0 * b1 * 0.5;
CalcT a1 = 2 * cs;
CalcT a2 = alph - 1;
SetStage( a1*a0, a2*a0, b0, b1*a0, b0 );
}
template<int channels>
void BiquadLowPass<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadLowPass<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadHighPass : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadHighPass<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT a0 = -1 / ( 1 + alph );
CalcT b1 = -( 1 + cs );
CalcT b0 = a0 * b1 * -0.5;
CalcT a1 = -2 * cs;
CalcT a2 = 1 - alph;
SetStage( a1*a0, a2*a0, b0, b1*a0, b0 );
}
template<int channels>
void BiquadHighPass<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadHighPass<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
// Constant skirt gain, peak gain=Q
template<int channels>
class BiquadBandPass1 : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadBandPass1<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT a0 = -1 / ( 1 + alph );
CalcT b0 = a0 * ( sn * -0.5 );
CalcT a1 = -2 * cs;
CalcT a2 = 1 - alph;
SetStage( a1*a0, a2*a0, b0, 0, -b0 );
}
template<int channels>
void BiquadBandPass1<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadBandPass1<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
// Constant 0dB peak gain
template<int channels>
class BiquadBandPass2 : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadBandPass2<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT b0 = -alph;
CalcT b2 = alph;
CalcT a0 = -1 / ( 1 + alph );
CalcT a1 = -2 * cs;
CalcT a2 = 1 - alph;
SetStage( a1*a0, a2*a0, b0*a0, 0, b2*a0 );
}
template<int channels>
void BiquadBandPass2<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadBandPass2<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadBandStop : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadBandStop<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT a0 = 1 / ( 1 + alph );
CalcT b1 = a0 * ( -2 * cs );
CalcT a2 = alph - 1;
SetStage( -b1, a2*a0, a0, b1, a0 );
}
template<int channels>
void BiquadBandStop<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadBandStop<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadAllPass: public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT q );
void SetupFast ( CalcT normFreq, CalcT q );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT q );
};
//--------------------------------------------------------------------------
template<int channels>
void BiquadAllPass<channels>::SetupCommon( CalcT sn, CalcT cs, CalcT q )
{
CalcT alph = sn / ( 2 * q );
CalcT b2 = 1 + alph;
CalcT a0 = 1 / b2;
CalcT b0 =( 1 - alph ) * a0;
CalcT b1 = -2 * cs * a0;
SetStage( -b1, -b0, b0, b1, b2*a0 );
}
template<int channels>
void BiquadAllPass<channels>::Setup( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
SetupCommon( sn, cs, q );
}
template<int channels>
void BiquadAllPass<channels>::SetupFast( CalcT normFreq, CalcT q )
{
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
SetupCommon( sn, cs, q );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadPeakEq: public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT dB, CalcT bandWidth );
void SetupFast ( CalcT normFreq, CalcT dB, CalcT bandWidth );
protected:
void SetupCommon ( CalcT sn, CalcT cs, CalcT alph, CalcT A );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadPeakEq<channels>::SetupCommon(
CalcT sn, CalcT cs, CalcT alph, CalcT A )
{
CalcT t=alph*A;
CalcT b0 = 1 - t;
CalcT b2 = 1 + t;
t=alph/A;
CalcT a0 = 1 / ( 1 + t );
CalcT a2 = t - 1;
CalcT b1 = a0 * ( -2 * cs );
CalcT a1 = -b1;
SetStage( a1, a2*a0, b0*a0, b1, b2*a0 );
}
template<int channels>
void BiquadPeakEq<channels>::Setup( CalcT normFreq, CalcT dB, CalcT bandWidth )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
CalcT alph = sn * sinh( kLn2/2 * bandWidth * w0/sn );
SetupCommon( sn, cs, alph, A );
}
template<int channels>
void BiquadPeakEq<channels>::SetupFast( CalcT normFreq, CalcT dB, CalcT bandWidth )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
CalcT alph = sn * sinh( kLn2/2 * bandWidth * w0/sn );
SetupCommon( sn, cs, alph, A );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadLowShelf : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 );
void SetupFast ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 );
protected:
void SetupCommon ( CalcT cs, CalcT A, CalcT sa );
};
//--------------------------------------------------------------------------
template<int channels>
inline void BiquadLowShelf<channels>::SetupCommon(
CalcT cs, CalcT A, CalcT sa )
{
CalcT An = A-1;
CalcT Ap = A+1;
CalcT Ancs = An*cs;
CalcT Apcs = Ap*cs;
CalcT b0 = A * (Ap - Ancs + sa );
CalcT b2 = A * (Ap - Ancs - sa );
CalcT b1 = 2 * A * (An - Apcs);
CalcT a2 = sa - (Ap + Ancs);
CalcT a0 = 1 / (Ap + Ancs + sa );
CalcT a1 = 2 * (An + Apcs);
SetStage( a1*a0, a2*a0, b0*a0, b1*a0, b2*a0 );
}
template<int channels>
void BiquadLowShelf<channels>::Setup( CalcT normFreq, CalcT dB, CalcT shelfSlope )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
CalcT al = sn / 2 * ::sqrt( (A + 1/A) * (1/shelfSlope - 1) + 2 );
CalcT sa = 2 * ::sqrt( A ) * al;
SetupCommon( cs, A, sa );
}
// This could be optimized further
template<int channels>
void BiquadLowShelf<channels>::SetupFast( CalcT normFreq, CalcT dB, CalcT shelfSlope )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
CalcT al = sn / 2 * fastsqrt1( (A + 1/A) * (1/shelfSlope - 1) + 2 );
CalcT sa = 2 * fastsqrt1( A ) * al;
SetupCommon( cs, A, sa );
}
//--------------------------------------------------------------------------
template<int channels>
class BiquadHighShelf : public Biquad<channels>
{
public:
void Setup ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 );
void SetupFast ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 );
protected:
void SetupCommon ( CalcT cs, CalcT A, CalcT sa );
};
//--------------------------------------------------------------------------
template<int channels>
void BiquadHighShelf<channels>::SetupCommon(
CalcT cs, CalcT A, CalcT sa )
{
CalcT An = A-1;
CalcT Ap = A+1;
CalcT Ancs = An*cs;
CalcT Apcs = Ap*cs;
CalcT b0 = A * (Ap + Ancs + sa );
CalcT b1 = -2 * A * (An + Apcs);
CalcT b2 = A * (Ap + Ancs - sa );
CalcT a0 = (Ap - Ancs + sa );
CalcT a2 = Ancs + sa - Ap;
CalcT a1 = -2 * (An - Apcs);
SetStage( a1/a0, a2/a0, b0/a0, b1/a0, b2/a0 );
}
template<int channels>
void BiquadHighShelf<channels>::Setup( CalcT normFreq, CalcT dB, CalcT shelfSlope )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT cs = cos(w0);
CalcT sn = sin(w0);
CalcT alph = sn / 2 * ::sqrt( (A + 1/A) * (1/shelfSlope - 1) + 2 );
CalcT sa = 2 * ::sqrt( A ) * alph;
SetupCommon( cs, A, sa );
}
template<int channels>
void BiquadHighShelf<channels>::SetupFast( CalcT normFreq, CalcT dB, CalcT shelfSlope )
{
CalcT A = pow( 10, dB/40 );
CalcT w0 = 2 * kPi * normFreq;
CalcT sn, cs;
fastsincos( w0, &sn, &cs );
CalcT alph = sn / 2 * fastsqrt1( (A + 1/A) * (1/shelfSlope - 1) + 2 );
CalcT sa = 2 * fastsqrt1( A ) * alph;
SetupCommon( cs, A, sa );
}
//--------------------------------------------------------------------------
//
// General N-Pole IIR Filter
//
//--------------------------------------------------------------------------
template<int stages, int channels>
class PoleFilter : public CascadeStages<stages, channels>
{
public:
PoleFilter();
virtual int CountPoles ( void )=0;
virtual int CountZeroes ( void )=0;
virtual Complex GetPole ( int i )=0;
virtual Complex GetZero ( int i )=0;
protected:
virtual Complex GetSPole ( int i, CalcT wc )=0;
protected:
// Determines the method of obtaining
// unity gain coefficients in the passband.
enum Hint
{
// No normalizating
hintNone,
// Use Brent's method to find the maximum
hintBrent,
// Use the response at a given frequency
hintPassband
};
Complex BilinearTransform ( const Complex &c );
Complex BandStopTransform ( int i, const Complex &c );
Complex BandPassTransform ( int i, const Complex &c );
Complex GetBandStopPole ( int i );
Complex GetBandStopZero ( int i );
Complex GetBandPassPole ( int i );
Complex GetBandPassZero ( int i );
void Normalize ( void );
void Prepare ( void );
virtual void BrentHint ( CalcT *w0, CalcT *w1 );
virtual CalcT PassbandHint( void );
protected:
Hint m_hint;
int m_n;
CalcT m_wc;
CalcT m_wc2;
};
//--------------------------------------------------------------------------
template<int stages, int channels>
inline PoleFilter<stages, channels>::PoleFilter( void )
{
m_hint=hintNone;
}
template<int stages, int channels>
inline Complex PoleFilter<stages, channels>::BilinearTransform( const Complex &c )
{
return (c+1.)/(-c+1.);
}
template<int stages, int channels>
inline Complex PoleFilter<stages, channels>::BandStopTransform( int i, const Complex &c )
{
CalcT a=cos((m_wc+m_wc2)*.5) /
cos((m_wc-m_wc2)*.5);
CalcT b=tan((m_wc-m_wc2)*.5);
Complex c2(0);
c2=addmul( c2, 4*(b*b+a*a-1), c );
c2+=8*(b*b-a*a+1);
c2*=c;
c2+=4*(a*a+b*b-1);
c2=Dsp::sqrt( c2 );
c2*=((i&1)==0)?.5:-.5;
c2+=a;
c2=addmul( c2, -a, c );
Complex c3( b+1 );
c3=addmul( c3, b-1, c );
return c2/c3;
}
template<int stages, int channels>
inline Complex PoleFilter<stages, channels>::BandPassTransform( int i, const Complex &c )
{
CalcT a= cos((m_wc+m_wc2)*0.5)/
cos((m_wc-m_wc2)*0.5);
CalcT b=1/tan((m_wc-m_wc2)*0.5);
Complex c2(0);
c2=addmul( c2, 4*(b*b*(a*a-1)+1), c );
c2+=8*(b*b*(a*a-1)-1);
c2*=c;
c2+=4*(b*b*(a*a-1)+1);
c2=Dsp::sqrt( c2 );
if ((i & 1) == 0)
c2=-c2;
c2=addmul( c2, 2*a*b, c );
c2+=2*a*b;
Complex c3(0);
c3=addmul( c3, 2*(b-1), c );
c3+=2*(1+b);
return c2/c3;
}
template<int stages, int channels>
Complex PoleFilter<stages, channels>::GetBandStopPole( int i )
{
Complex c=GetSPole( i/2, kPi_2 );
c=BilinearTransform( c );
c=BandStopTransform( i, c );
return c;
}
template<int stages, int channels>
Complex PoleFilter<stages, channels>::GetBandStopZero( int i )
{
return BandStopTransform( i, Complex( -1 ) );
}
template<int stages, int channels>
Complex PoleFilter<stages, channels>::GetBandPassPole( int i )
{
Complex c=GetSPole( i/2, kPi_2 );
c=BilinearTransform( c );
c=BandPassTransform( i, c );
return c;
}
template<int stages, int channels>
Complex PoleFilter<stages, channels>::GetBandPassZero( int i )
{
return Complex( (i>=m_n)?1:-1 );
}
template<int stages, int channels>
void PoleFilter<stages, channels>::Normalize( void )
{
switch( m_hint )
{
default:
case hintNone:
break;
case hintPassband:
{
CalcT w=PassbandHint();
ResponseFunctor f(this);
CalcT mag=-f(w);
CascadeStages::Normalize( 1/mag );
}
break;
case hintBrent:
{
ResponseFunctor f(this);
CalcT w0, w1, wmin, mag;
BrentHint( &w0, &w1 );
mag=-BrentMinimize( f, w0, w1, 1e-4, wmin );
CascadeStages::Normalize( 1/mag );
}
break;
}
}
template<int stages, int channels>
void PoleFilter<stages, channels>::Prepare( void )
{
if( m_wc2<1e-8 )
m_wc2=1e-8;
if( m_wc >kPi-1e-8 )
m_wc =kPi-1e-8;
Reset();
Complex c;
int poles=CountPoles();
for( int i=0;i<poles;i++ )
{
c=GetPole( i );
if( ::abs(c.imag())<1e-6 )
c=Complex( c.real(), 0 );
if( c.imag()==0 )
SetAStage( c.real(), 0 );
else if( c.imag()>0 )
SetAStage( 2*c.real(), -Dsp::norm(c) );
}
int zeroes=CountZeroes();
for( int i=0;i<zeroes;i++ )
{
c=GetZero( i );
if( ::abs(c.imag())<1e-6 )
c=Complex( c.real(), 0 );
if( c.imag()==0 )
SetBStage( -c.real(), 1, 0 );
else if( c.imag()>0 )
SetBStage( Dsp::norm(c), -2*c.real(), 1 );
}
Normalize();
}
template<int stages, int channels>
void PoleFilter<stages, channels>::BrentHint( CalcT *w0, CalcT *w1 )
{
// best that this never executes
*w0=1e-4;
*w1=kPi-1e-4;
}
template<int stages, int channels>
CalcT PoleFilter<stages, channels>::PassbandHint( void )
{
// should never get here
assert( 0 );
return kPi_2;
}
//--------------------------------------------------------------------------
//
// Butterworth Response IIR Filter
//
//--------------------------------------------------------------------------
// Butterworth filter response characteristic.
// Maximally flat magnitude response in the passband at the
// expense of a more shallow rolloff in comparison to other types.
template<int poles, int channels>
class Butterworth : public PoleFilter<int((poles+1)/2), channels>
{
public:
Butterworth();
// cutoffFreq = freq / sampleRate
void Setup ( CalcT cutoffFreq );
virtual int CountPoles ( void );
virtual int CountZeroes ( void );
virtual Complex GetPole ( int i );
protected:
Complex GetSPole ( int i, CalcT wc );
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Butterworth<poles, channels>::Butterworth( void )
{
m_hint=hintPassband;
}
template<int poles, int channels>
void Butterworth<poles, channels>::Setup( CalcT cutoffFreq )
{
m_n=poles;
m_wc=2*kPi*cutoffFreq;
Prepare();
}
template<int poles, int channels>
int Butterworth<poles, channels>::CountPoles( void )
{
return poles;
}
template<int poles, int channels>
int Butterworth<poles, channels>::CountZeroes( void )
{
return poles;
}
template<int poles, int channels>
Complex Butterworth<poles, channels>::GetPole( int i )
{
return BilinearTransform( GetSPole( i, m_wc ) );
}
template<int poles, int channels>
Complex Butterworth<poles, channels>::GetSPole( int i, CalcT wc )
{
return polar( tan(wc*0.5), kPi_2+(2*i+1)*kPi/(2*m_n) );
}
//--------------------------------------------------------------------------
// Low Pass Butterworth filter
// Stable up to 53 poles (frequency min=0.13% of Nyquist)
template<int poles, int channels>
class ButterLowPass : public Butterworth<poles, channels>
{
public:
Complex GetZero ( int i );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Complex ButterLowPass<poles, channels>::GetZero( int i )
{
return Complex( -1 );
}
template<int poles, int channels>
CalcT ButterLowPass<poles, channels>::PassbandHint( void )
{
return 0;
}
//--------------------------------------------------------------------------
// High Pass Butterworth filter
// Maximally flat magnitude response in the passband.
// Stable up to 110 poles (frequency max=97% of Nyquist)
template<int poles, int channels>
class ButterHighPass : public Butterworth<poles, channels>
{
public:
Complex GetZero( int i );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Complex ButterHighPass<poles, channels>::GetZero( int i )
{
return Complex( 1 );
}
template<int poles, int channels>
CalcT ButterHighPass<poles, channels>::PassbandHint( void )
{
return kPi;
}
//--------------------------------------------------------------------------
// Band Pass Butterworth filter
// Stable up to 80 pairs
template<int pairs, int channels>
class ButterBandPass : public Butterworth<pairs*2, channels>
{
public:
// centerFreq = freq / sampleRate
// normWidth = freqWidth / sampleRate
void Setup ( CalcT centerFreq, CalcT normWidth );
virtual int CountPoles ( void );
virtual int CountZeroes ( void );
virtual Complex GetPole ( int i );
virtual Complex GetZero ( int i );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int pairs, int channels>
void ButterBandPass<pairs, channels>::Setup( CalcT centerFreq, CalcT normWidth )
{
m_n=pairs;
CalcT angularWidth=2*kPi*normWidth;
m_wc2=2*kPi*centerFreq-(angularWidth/2);
m_wc =m_wc2+angularWidth;
Prepare();
}
template<int pairs, int channels>
int ButterBandPass<pairs, channels>::CountPoles( void )
{
return pairs*2;
}
template<int pairs, int channels>
int ButterBandPass<pairs, channels>::CountZeroes( void )
{
return pairs*2;
}
template<int pairs, int channels>
Complex ButterBandPass<pairs, channels>::GetPole( int i )
{
return GetBandPassPole( i );
}
template<int pairs, int channels>
Complex ButterBandPass<pairs, channels>::GetZero( int i )
{
return GetBandPassZero( i );
}
template<int poles, int channels>
CalcT ButterBandPass<poles, channels>::PassbandHint( void )
{
return (m_wc+m_wc2)/2;
}
//--------------------------------------------------------------------------
// Band Stop Butterworth filter
// Stable up to 109 pairs
template<int pairs, int channels>
class ButterBandStop : public Butterworth<pairs*2, channels>
{
public:
// centerFreq = freq / sampleRate
// normWidth = freqWidth / sampleRate
void Setup ( CalcT centerFreq, CalcT normWidth );
virtual int CountPoles ( void );
virtual int CountZeroes ( void );
virtual Complex GetPole ( int i );
virtual Complex GetZero ( int i );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int pairs, int channels>
void ButterBandStop<pairs, channels>::Setup( CalcT centerFreq, CalcT normWidth )
{
m_n=pairs;
CalcT angularWidth=2*kPi*normWidth;
m_wc2=2*kPi*centerFreq-(angularWidth/2);
m_wc =m_wc2+angularWidth;
Prepare();
}
template<int pairs, int channels>
int ButterBandStop<pairs, channels>::CountPoles( void )
{
return pairs*2;
}
template<int pairs, int channels>
int ButterBandStop<pairs, channels>::CountZeroes( void )
{
return pairs*2;
}
template<int pairs, int channels>
Complex ButterBandStop<pairs, channels>::GetPole( int i )
{
return GetBandStopPole( i );
}
template<int pairs, int channels>
Complex ButterBandStop<pairs, channels>::GetZero( int i )
{
return GetBandStopZero( i );
}
template<int poles, int channels>
CalcT ButterBandStop<poles, channels>::PassbandHint( void )
{
if( (m_wc+m_wc2)/2<kPi_2 )
return kPi;
else
return 0;
}
//--------------------------------------------------------------------------
//
// Chebyshev Response IIR Filter
//
//--------------------------------------------------------------------------
// Type I Chebyshev filter characteristic.
// Minimum error between actual and ideal response at the expense of
// a user-definable amount of ripple in the passband.
template<int poles, int channels>
class Chebyshev1 : public PoleFilter<int((poles+1)/2), channels>
{
public:
Chebyshev1();
// cutoffFreq = freq / sampleRate
virtual void Setup ( CalcT cutoffFreq, CalcT rippleDb );
virtual int CountPoles ( void );
virtual int CountZeroes ( void );
virtual Complex GetPole ( int i );
virtual Complex GetZero ( int i );
protected:
void SetupCommon ( CalcT rippleDb );
virtual Complex GetSPole ( int i, CalcT wc );
protected:
CalcT m_sgn;
CalcT m_eps;
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Chebyshev1<poles, channels>::Chebyshev1()
{
m_hint=hintBrent;
}
template<int poles, int channels>
void Chebyshev1<poles, channels>::Setup( CalcT cutoffFreq, CalcT rippleDb )
{
m_n=poles;
m_wc=2*kPi*cutoffFreq;
SetupCommon( rippleDb );
}
template<int poles, int channels>
void Chebyshev1<poles, channels>::SetupCommon( CalcT rippleDb )
{
m_eps=::sqrt( 1/::exp( -rippleDb*0.1*kLn10 )-1 );
Prepare();
// This moves the bottom of the ripples to 0dB gain
//CascadeStages::Normalize( pow( 10, rippleDb/20.0 ) );
}
template<int poles, int channels>
int Chebyshev1<poles, channels>::CountPoles( void )
{
return poles;
}
template<int poles, int channels>
int Chebyshev1<poles, channels>::CountZeroes( void )
{
return poles;
}
template<int poles, int channels>
Complex Chebyshev1<poles, channels>::GetPole( int i )
{
return BilinearTransform( GetSPole( i, m_wc ) )*m_sgn;
}
template<int poles, int channels>
Complex Chebyshev1<poles, channels>::GetZero( int i )
{
return Complex( -m_sgn );
}
template<int poles, int channels>
Complex Chebyshev1<poles, channels>::GetSPole( int i, CalcT wc )
{
int n = m_n;
CalcT ni = 1.0/n;
CalcT alpha = 1/m_eps+::sqrt(1+1/(m_eps*m_eps));
CalcT pn = pow( alpha, ni );
CalcT nn = pow( alpha, -ni );
CalcT a = 0.5*( pn - nn );
CalcT b = 0.5*( pn + nn );
CalcT theta = kPi_2 + (2*i+1) * kPi/(2*n);
Complex c = polar( tan( 0.5*(m_sgn==-1?(kPi-wc):wc) ), theta );
return Complex( a*c.real(), b*c.imag() );
}
//--------------------------------------------------------------------------
// Low Pass Chebyshev Type I filter
template<int poles, int channels>
class Cheby1LowPass : public Chebyshev1<poles, channels>
{
public:
Cheby1LowPass();
void Setup ( CalcT cutoffFreq, CalcT rippleDb );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Cheby1LowPass<poles, channels>::Cheby1LowPass()
{
m_sgn=1;
m_hint=hintPassband;
}
template<int poles, int channels>
void Cheby1LowPass<poles, channels>::Setup( CalcT cutoffFreq, CalcT rippleDb )
{
Chebyshev1::Setup( cutoffFreq, rippleDb );
// move peak of ripple down to 0dB
if( !(poles&1) )
CascadeStages::Normalize( pow( 10, -rippleDb/20.0 ) );
}
template<int poles, int channels>
CalcT Cheby1LowPass<poles, channels>::PassbandHint( void )
{
return 0;
}
//--------------------------------------------------------------------------
// High Pass Chebyshev Type I filter
template<int poles, int channels>
class Cheby1HighPass : public Chebyshev1<poles, channels>
{
public:
Cheby1HighPass();
void Setup ( CalcT cutoffFreq, CalcT rippleDb );
protected:
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int poles, int channels>
Cheby1HighPass<poles, channels>::Cheby1HighPass()
{
m_sgn=-1;
m_hint=hintPassband;
}
template<int poles, int channels>
void Cheby1HighPass<poles, channels>::Setup( CalcT cutoffFreq, CalcT rippleDb )
{
Chebyshev1::Setup( cutoffFreq, rippleDb );
// move peak of ripple down to 0dB
if( !(poles&1) )
CascadeStages::Normalize( pow( 10, -rippleDb/20.0 ) );
}
template<int poles, int channels>
CalcT Cheby1HighPass<poles, channels>::PassbandHint( void )
{
return kPi;
}
//--------------------------------------------------------------------------
// Band Pass Chebyshev Type I filter
template<int pairs, int channels>
class Cheby1BandPass : public Chebyshev1<pairs*2, channels>
{
public:
Cheby1BandPass();
void Setup ( CalcT centerFreq, CalcT normWidth, CalcT rippleDb );
int CountPoles ( void );
int CountZeroes ( void );
Complex GetPole ( int i );
Complex GetZero ( int i );
protected:
void BrentHint ( CalcT *w0, CalcT *w1 );
//CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int pairs, int channels>
Cheby1BandPass<pairs, channels>::Cheby1BandPass()
{
m_sgn=1;
m_hint=hintBrent;
}
template<int pairs, int channels>
void Cheby1BandPass<pairs, channels>::Setup( CalcT centerFreq, CalcT normWidth, CalcT rippleDb )
{
m_n=pairs;
CalcT angularWidth=2*kPi*normWidth;
m_wc2=2*kPi*centerFreq-(angularWidth/2);
m_wc =m_wc2+angularWidth;
SetupCommon( rippleDb );
}
template<int pairs, int channels>
int Cheby1BandPass<pairs, channels>::CountPoles( void )
{
return pairs*2;
}
template<int pairs, int channels>
int Cheby1BandPass<pairs, channels>::CountZeroes( void )
{
return pairs*2;
}
template<int pairs, int channels>
Complex Cheby1BandPass<pairs, channels>::GetPole( int i )
{
return GetBandPassPole( i );
}
template<int pairs, int channels>
Complex Cheby1BandPass<pairs, channels>::GetZero( int i )
{
return GetBandPassZero( i );
}
template<int poles, int channels>
void Cheby1BandPass<poles, channels>::BrentHint( CalcT *w0, CalcT *w1 )
{
CalcT d=1e-4*(m_wc-m_wc2)/2;
*w0=m_wc2+d;
*w1=m_wc-d;
}
/*
// Unfortunately, this doesn't work at the frequency extremes
// Maybe we can inverse pre-warp the center point to make sure
// it stays put after bilinear and bandpass transformation.
template<int poles, int channels>
CalcT Cheby1BandPass<poles, channels>::PassbandHint( void )
{
return (m_wc+m_wc2)/2;
}
*/
//--------------------------------------------------------------------------
// Band Stop Chebyshev Type I filter
template<int pairs, int channels>
class Cheby1BandStop : public Chebyshev1<pairs*2, channels>
{
public:
Cheby1BandStop();
void Setup ( CalcT centerFreq, CalcT normWidth, CalcT rippleDb );
int CountPoles ( void );
int CountZeroes ( void );
Complex GetPole ( int i );
Complex GetZero ( int i );
protected:
void BrentHint ( CalcT *w0, CalcT *w1 );
CalcT PassbandHint ( void );
};
//--------------------------------------------------------------------------
template<int pairs, int channels>
Cheby1BandStop<pairs, channels>::Cheby1BandStop()
{
m_sgn=1;
m_hint=hintPassband;
}
template<int pairs, int channels>
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)));
}
2 comment(s) | 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);
2 comment(s) | 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 |
|
|
|
|
|
|
|
|
Branchless Clipping
Type : Clipping at 0dB, with none of the usual 'if..then..'
References : Posted by musicdsp[AT]dsparsons[DOT]co[DOT]uk
Notes :
I was working on something that I wanted to ensure that the signal never went above 0dB, and a branchless solution occurred to me.
It works by playing with the structure of a single type, shifting the sign bit down to make a new mulitplicand.
calling MaxZerodB(mydBSample) will ensure that it will never stray over 0dB.
By playing with signs or adding/removing offsets, this offers a complete branchless limiting solution, no matter whether dB or not (after all, they're all just numbers...).
eg:
Limit to <=0 : sample:=MaxZerodB(sample);
Limit to <=3 : sample:=MaxZerodB(sample-3)+3;
Limit to <=-4 : sample:=MaxZerodB(sample+4)-4;
Limit to >=0 : sample:=-MaxZerodB(-sample);
Limit to >=2 : sample:=-MaxZerodB(-sample+2)+2;
Limit to >=-1.5: sample:=-MaxZerodB(-sample-1.5)-1.5;
Whether it actually saves any CPU cycles remains to be seen, but it was an interesting diversion for half an hour :)
[Translating from pascal to other languages shouldn't be too hard, and for doubles, you'll need to fiddle it abit :)]
Code :
function MaxZerodB(dBin:single):single;
var tmp:longint;
begin
//given that leftmost bit of a longint indicates the negative,
// if we shift that down to bit0, and multiply dBin by that
// it will return dBin, or zero :)
tmp:=(longint((@dBin)^) and $80000000) shr 31;
result:=dBin*tmp;
end;
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Butterworth
Type : LPF 24dB/Oct
References : Posted by Christian[at]savioursofsoul[dot]de
Code :
First calculate the prewarped digital frequency:
K = tan(Pi * Frequency / Samplerate);
Now calc some intermediate variables: (see 'Factors of Polynoms' at http://en.wikipedia.org/wiki/Butterworth_filter, especially if you want a higher order like 48dB/Oct)
a = 0.76536686473 * Q * K;
b = 1.84775906502 * Q * K;
K = K*K; (to optimize it a little bit)
Calculate the first biquad:
A0 = (K+a+1);
A1 = 2*(1-K);
A2 =(a-K-1);
B0 = K;
B1 = 2*B0;
B2 = B0;
Calculate the second biquad:
A3 = (K+b+1);
A4 = 2*(1-K);
A5 = (b-K-1);
B3 = K;
B4 = 2*B3;
B5 = B3;
Then calculate the output as follows:
Stage1 = B0*Input + State0;
State0 = B1*Input + A1/A0*Stage1 + State1;
State1 = B2*Input + A2/A0*Stage1;
Output = B3*Stage1 + State2;
State2 = B4*Stage1 + A4/A3*Output + State2;
State3 = B5*Stage1 + A5/A3*Output;
9 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Butterworth Optimized C++ Class
Type : 24db Resonant Lowpass
References : Posted by neotec
Notes :
This ist exactly the same as posted by "Zxform" (filters004.txt). The only difference is, that this version is an optimized one.
Parameters:
Cutoff [0.f -> Nyquist.f]
Resonance [0.f -> 1.f]
There are some minima and maxima defined, to make ist sound nice in all situations. This class is part of some of my VST Plugins, and works well and executes fast.
Code :
// FilterButterworth24db.h
#pragma once
class CFilterButterworth24db
{
public:
CFilterButterworth24db(void);
~CFilterButterworth24db(void);
void SetSampleRate(float fs);
void Set(float cutoff, float q);
float Run(float input);
private:
float t0, t1, t2, t3;
float coef0, coef1, coef2, coef3;
float history1, history2, history3, history4;
float gain;
float min_cutoff, max_cutoff;
};
// FilterButterworth24db.cpp
#include <math.h>
#define BUDDA_Q_SCALE 6.f
#include "FilterButterworth24db.h"
CFilterButterworth24db::CFilterButterworth24db(void)
{
this->history1 = 0.f;
this->history2 = 0.f;
this->history3 = 0.f;
this->history4 = 0.f;
this->SetSampleRate(44100.f);
this->Set(22050.f, 0.0);
}
CFilterButterworth24db::~CFilterButterworth24db(void)
{
}
void CFilterButterworth24db::SetSampleRate(float fs)
{
float pi = 4.f * atanf(1.f);
this->t0 = 4.f * fs * fs;
this->t1 = 8.f * fs * fs;
this->t2 = 2.f * fs;
this->t3 = pi / fs;
this->min_cutoff = fs * 0.01f;
this->max_cutoff = fs * 0.45f;
}
void CFilterButterworth24db::Set(float cutoff, float q)
{
if (cutoff < this->min_cutoff)
cutoff = this->min_cutoff;
else if(cutoff > this->max_cutoff)
cutoff = this->max_cutoff;
if(q < 0.f)
q = 0.f;
else if(q > 1.f)
q = 1.f;
float wp = this->t2 * tanf(this->t3 * cutoff);
float bd, bd_tmp, b1, b2;
q *= BUDDA_Q_SCALE;
q += 1.f;
b1 = (0.765367f / q) / wp;
b2 = 1.f / (wp * wp);
bd_tmp = this->t0 * b2 + 1.f;
bd = 1.f / (bd_tmp + this->t2 * b1);
this->gain = bd * 0.5f;
this->coef2 = (2.f - this->t1 * b2);
this->coef0 = this->coef2 * bd;
this->coef1 = (bd_tmp - this->t2 * b1) * bd;
b1 = (1.847759f / q) / wp;
bd = 1.f / (bd_tmp + this->t2 * b1);
this->gain *= bd;
this->coef2 *= bd;
this->coef3 = (bd_tmp - this->t2 * b1) * bd;
}
float CFilterButterworth24db::Run(float input)
{
float output = input * this->gain;
float new_hist;
output -= this->history1 * this->coef0;
new_hist = output - this->history2 * this->coef1;
output = new_hist + this->history1 * 2.f;
output += this->history2;
this->history2 = this->history1;
this->history1 = new_hist;
output -= this->history3 * this->coef2;
new_hist = output - this->history4 * this->coef3;
output = new_hist + this->history3 * 2.f;
output += this->history4;
this->history4 = this->history3;
this->history3 = new_hist;
return output;
}
12 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
C# Oscilator class
Type : Sine, Saw, Variable Pulse, Triangle, C64 Noise
References : Posted by neotec
Notes :
Parameters:
Pitch: The Osc's pitch in Cents [0 - 14399] startig at A -> 6.875Hz
Pulsewidth: [0 - 65535] -> 0% to 99.99%
Value: The last Output value, a set to this property 'syncs' the Oscilator
Code :
public class SynthOscilator
{
public enum OscWaveformType
{
SAW, PULSE, TRI, NOISE, SINE
}
public int Pitch
{
get
{
return this._Pitch;
}
set
{
this._Pitch = this.MinMax(0, value, 14399);
this.OscStep = WaveSteps[this._Pitch];
}
}
public int PulseWidth
{
get
{
return this._PulseWidth;
}
set
{
this._PulseWidth = this.MinMax(0, value, 65535);
}
}
public OscWaveformType Waveform
{
get
{
return this._WaveForm;
}
set
{
this._WaveForm = value;
}
}
public int Value
{
get
{
return this._Value;
}
set
{
this._Value = 0;
this.OscNow = 0;
}
}
private int _Pitch;
private int _PulseWidth;
private int _Value;
private OscWaveformType _WaveForm;
private int OscNow;
private int OscStep;
private int ShiftRegister;
public const double BaseFrequence = 6.875;
public const int SampleRate = 44100;
public static int[] WaveSteps = new int[0];
public static int[] SineTable = new int[0];
public SynthOscilator()
{
if (WaveSteps.Length == 0)
this.CalcSteps();
if (SineTable.Length == 0)
this.CalcSine();
this._Pitch = 7200;
this._PulseWidth = 32768;
this._WaveForm = OscWaveformType.SAW;
this.ShiftRegister = 0x7ffff8;
this.OscNow = 0;
this.OscStep = WaveSteps[this._Pitch];
this._Value = 0;
}
private void CalcSteps()
{
WaveSteps = new int[14400];
for (int i = 0; i < 14400; i++)
{
double t0, t1, t2;
t0 = Math.Pow(2.0, (double)i / 1200.0);
t1 = BaseFrequence * t0;
t2 = (t1 * 65536.0) / (double)this.SampleRate;
WaveSteps[i] = (int)Math.Round(t2 * 4096.0);
}
}
private void CalcSine()
{
SineTable = new int[65536];
double s = Math.PI / 32768.0;
for (int i = 0; i < 65536; i++)
{
double v = Math.Sin((double)i * s) * 32768.0;
int t = (int)Math.Round(v) + 32768;
if (t < 0)
t = 0;
else if (t > 65535)
t = 65535;
SineTable[i] = t;
}
}
public override int Run()
{
int ret = 32768;
int osc = this.OscNow >> 12;
switch (this._WaveForm)
{
case OscWaveformType.SAW:
ret = osc;
break;
case OscWaveformType.PULSE:
if (osc < this.PulseWidth)
ret = 65535;
else
ret = 0;
break;
case OscWaveformType.TRI:
if (osc < 32768)
ret = osc << 1;
else
ret = 131071 - (osc << 1);
break;
case OscWaveformType.NOISE:
ret = ((this.ShiftRegister & 0x400000) >> 11) |
((this.ShiftRegister & 0x100000) >> 10) |
((this.ShiftRegister & 0x010000) >> 7) |
((this.ShiftRegister & 0x002000) >> 5) |
((this.ShiftRegister & 0x000800) >> 4) |
((this.ShiftRegister & 0x000080) >> 1) |
((this.ShiftRegister & 0x000010) << 1) |
((this.ShiftRegister & 0x000004) << 2);
ret <<= 4;
break;
case OscWaveformType.SINE:
ret = SynthTools.SineTable[osc];
break;
default:
break;
}
this.OscNow += this.OscStep;
if (this.OscNow > 0xfffffff)
{
int bit0 = ((this.ShiftRegister >> 22) ^ (this.ShiftRegister >> 17)) & 0x1;
this.ShiftRegister <<= 1;
this.ShiftRegister &= 0x7fffff;
this.ShiftRegister |= bit0;
}
this.OscNow &= 0xfffffff;
this._Value = ret - 32768;
return this._Value;
}
public int MinMax(int a, int b, int c)
{
if (b < a)
return a;
else if (b > c)
return c;
else
return b;
}
}
no comments on this item | 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;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
C-Weighed Filter
Type : digital implementation (after bilinear transform)
References : Posted by Christian@savioursofsoul.de
Notes :
unoptimized version!
Code :
First prewarp the frequency of both poles:
K1 = tan(0.5*Pi*20.6 / SampleRate) // for 20.6Hz
K2 = tan(0.5*Pi*12200 / SampleRate) // for 12200Hz
Then calculate the both biquads:
b0 = 1
b1 = 0
b2 =-1
a0 = ((K1+1)*(K1+1)*(K2+1)*(K2+1));
a1 =-4*(K1*K1*K2*K2+K1*K1*K2+K1*K2*K2-K1-K2-1)*t;
a2 =- ((K1-1)*(K1-1)*(K2-1)*(K2-1))*t;
and:
b3 = 1
b4 = 0
b5 =-1
a3 = ((K1+1)*(K1+1)*(K2+1)*(K2+1));
a4 =-4*(K1*K1*K2*K2+K1*K1*K2+K1*K2*K2-K1-K2-1)*t;
a5 =- ((K1-1)*(K1-1)*(K2-1)*(K2-1))*t;
Now use an equation for calculating the biquads like this:
Stage1 = b0*Input + State0;
State0 = + a1/a0*Stage1 + State1;
State1 = b2*Input + a2/a0*Stage1;
Output = b3*Stage1 + State2;
State2 = + a4/a3*Output + State2;
State3 = b5*Stage1 + a5/a3*Output;
1 comment(s) | 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 |
|
|
|
|
|
|
|
|
chebyshev waveshaper (using their recursive definition)
Type : chebyshev
References : Posted by mdsp
Notes :
someone asked for it on kvr-audio.
I use it in an unreleased additive synth.
There's no oversampling needed in my case since I feed it with a pure sinusoid and I control the order to not have frequencies above Fs/2. Otherwise you should oversample by the order you'll use in the function or bandlimit the signal before the waveshaper. unless you really want that aliasing effect... :)
I hope the code is self-explaining, otherwise there's plenty of sites explaining chebyshev polynoms and their applications.
Code :
float chebyshev(float x, float A[], int order)
{
// To = 1
// T1 = x
// Tn = 2.x.Tn-1 - Tn-2
// out = sum(Ai*Ti(x)) , i C {1,..,order}
float Tn_2 = 1.0f;
float Tn_1 = x;
float Tn;
float out = A[0]*Tn_1;
for(int n=2;n<=order;n++)
{
Tn = 2.0f*x*Tn_1 - Tn_2;
out += A[n-1]*Tn;
Tn_2 = Tn_1;
Tn_1 = Tn;
}
return out;
}
4 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;
}
}
6 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);
4 comment(s) | 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
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Cure for malicious samples
Type : Filters Denormals, NaNs, Infinities
References : Posted by urs[AT]u-he[DOT]com
Notes :
A lot of incidents can happen during processing samples. A nasty one is denormalization, which makes cpus consume insanely many cycles for easiest instructions.
But even worse, if you have NaNs or Infinities inside recursive structures, maybe due to division by zero, all subsequent samples that are multiplied with these values will get "infected" and become NaN or Infinity. Your sound makes BLIPPP and that was it, silence from the speakers.
Thus I've written a small function that sets all of these cases to 0.0f.
You'll notice that I treat a buffer of floats as unsigned integers. And I avaoid branches by using comparison results as 0 or 1.
When compiled with GCC, this function should not create any "hidden" branches, but you should verify the assembly code anyway. Sometimes some parenthesis do the trick...
;) Urs
Code :
#ifndef UInt32
#define UInt32 unsigned int
#endif
void erase_All_NaNs_Infinities_And_Denormals( float* inSamples, int& inNumberOfSamples )
{
UInt32* inArrayOfFloats = (UInt32*) inSamples;
for ( int i = 0; i < inNumberOfSamples; i++ )
{
UInt32 sample = *inArrayOfFloats;
UInt32 exponent = sample & 0x7F800000;
// exponent < 0x7F800000 is 0 if NaN or Infinity, otherwise 1
// exponent > 0 is 0 if denormalized, otherwise 1
int aNaN = exponent < 0x7F800000;
int aDen = exponent > 0;
*inArrayOfFloats++ = sample * ( aNaN & aDen );
}
}
4 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 |
|
|
|
|
|
|
|
|
Delphi Class implementation of the RBJ filters
Type : Delphi class implementation of the RBJ filters
References : Posted by veryangrymobster@hotmail.com
Notes :
I haven't tested this code thoroughly as it's pretty much a straight conversion from Arguru c++ implementation.
Code :
{
RBJ Audio EQ Cookbook Filters
A pascal conversion of arguru[AT]smartelectronix[DOT]com's
c++ implementation.
WARNING:This code is not FPU undernormalization safe.
Filter Types
0-LowPass
1-HiPass
2-BandPass CSG
3-BandPass CZPG
4-Notch
5-AllPass
6-Peaking
7-LowShelf
8-HiShelf
}
unit uRbjEqFilters;
interface
uses math;
type
TRbjEqFilter=class
private
b0a0,b1a0,b2a0,a1a0,a2a0:single;
in1,in2,ou1,ou2:single;
fSampleRate:single;
fMaxBlockSize:integer;
fFilterType:integer;
fFreq,fQ,fDBGain:single;
fQIsBandWidth:boolean;
procedure SetQ(NewQ:single);
public
out1:array of single;
constructor create(SampleRate:single;MaxBlockSize:integer);
procedure CalcFilterCoeffs(pFilterType:integer;pFreq,pQ,pDBGain:single;pQIsBandWidth:boolean);overload;
procedure CalcFilterCoeffs;overload;
function Process(input:single):single; overload;
procedure Process(Input:psingle;sampleframes:integer); overload;
property FilterType:integer read fFilterType write fFilterType;
property Freq:single read fFreq write fFreq;
property q:single read fQ write SetQ;
property DBGain:single read fDBGain write fDBGain;
property QIsBandWidth:boolean read fQIsBandWidth write fQIsBandWidth;
end;
implementation
constructor TRbjEqFilter.create(SampleRate:single;MaxBlockSize:integer);
begin
fMaxBlockSize:=MaxBlockSize;
setLength(out1,fMaxBlockSize);
fSampleRate:=SampleRate;
fFilterType:=0;
fFreq:=500;
fQ:=0.3;
fDBGain:=0;
fQIsBandWidth:=true;
in1:=0;
in2:=0;
ou1:=0;
ou2:=0;
end;
procedure TRbjEqFilter.SetQ(NewQ:single);
begin
fQ:=(1-NewQ)*0.98;
end;
procedure TRbjEqFilter.CalcFilterCoeffs(pFilterType:integer;pFreq,pQ,pDBGain:single;pQIsBandWidth:boolean);
begin
FilterType:=pFilterType;
Freq:=pFreq;
Q:=pQ;
DBGain:=pDBGain;
QIsBandWidth:=pQIsBandWidth;
CalcFilterCoeffs;
end;
procedure TRbjEqFilter.CalcFilterCoeffs;
var
alpha,a0,a1,a2,b0,b1,b2:single;
A,beta,omega,tsin,tcos:single;
begin
//peaking, LowShelf or HiShelf
if fFilterType>=6 then
begin
A:=power(10.0,(DBGain/40.0));
omega:=2*pi*fFreq/fSampleRate;
tsin:=sin(omega);
tcos:=cos(omega);
if fQIsBandWidth then
alpha:=tsin*sinh(log2(2.0)/2.0*fQ*omega/tsin)
else
alpha:=tsin/(2.0*fQ);
beta:=sqrt(A)/fQ;
// peaking
if fFilterType=6 then
begin
b0:=1.0+alpha*A;
b1:=-2.0*tcos;
b2:=1.0-alpha*A;
a0:=1.0+alpha/A;
a1:=-2.0*tcos;
a2:=1.0-alpha/A;
end else
// lowshelf
if fFilterType=7 then
begin
b0:=(A*((A+1.0)-(A-1.0)*tcos+beta*tsin));
b1:=(2.0*A*((A-1.0)-(A+1.0)*tcos));
b2:=(A*((A+1.0)-(A-1.0)*tcos-beta*tsin));
a0:=((A+1.0)+(A-1.0)*tcos+beta*tsin);
a1:=(-2.0*((A-1.0)+(A+1.0)*tcos));
a2:=((A+1.0)+(A-1.0)*tcos-beta*tsin);
end;
// hishelf
if fFilterType=8 then
begin
b0:=(A*((A+1.0)+(A-1.0)*tcos+beta*tsin));
b1:=(-2.0*A*((A-1.0)+(A+1.0)*tcos));
b2:=(A*((A+1.0)+(A-1.0)*tcos-beta*tsin));
a0:=((A+1.0)-(A-1.0)*tcos+beta*tsin);
a1:=(2.0*((A-1.0)-(A+1.0)*tcos));
a2:=((A+1.0)-(A-1.0)*tcos-beta*tsin);
end;
end else //other filter types
begin
omega:=2*pi*fFreq/fSampleRate;
tsin:=sin(omega);
tcos:=cos(omega);
if fQIsBandWidth then
alpha:=tsin*sinh(log2(2)/2*fQ*omega/tsin)
else
alpha:=tsin/(2*fQ);
//lowpass
if fFilterType=0 then
begin
b0:=(1-tcos)/2;
b1:=1-tcos;
b2:=(1-tcos)/2;
a0:=1+alpha;
a1:=-2*tcos;
a2:=1-alpha;
end else //hipass
if fFilterType=1 then
begin
b0:=(1+tcos)/2;
b1:=-(1+tcos);
b2:=(1+tcos)/2;
a0:=1+alpha;
a1:=-2*tcos;
a2:=1-alpha;
end else //bandpass CSG
if fFilterType=2 then
begin
b0:=tsin/2;
b1:=0;
b2:=-tsin/2;
a0:=1+alpha;
a1:=-1*tcos;
a2:=1-alpha;
end else //bandpass CZPG
if fFilterType=3 then
begin
b0:=alpha;
b1:=0.0;
b2:=-alpha;
a0:=1.0+alpha;
a1:=-2.0*tcos;
a2:=1.0-alpha;
end else //notch
if fFilterType=4 then
begin
b0:=1.0;
b1:=-2.0*tcos;
b2:=1.0;
a0:=1.0+alpha;
a1:=-2.0*tcos;
a2:=1.0-alpha;
end else //allpass
if fFilterType=5 then
begin
b0:=1.0-alpha;
b1:=-2.0*tcos;
b2:=1.0+alpha;
a0:=1.0+alpha;
a1:=-2.0*tcos;
a2:=1.0-alpha;
end;
end;
b0a0:=b0/a0;
b1a0:=b1/a0;
b2a0:=b2/a0;
a1a0:=a1/a0;
a2a0:=a2/a0;
end;
function TRbjEqFilter.Process(input:single):single;
var
LastOut:single;
begin
// filter
LastOut:= b0a0*input + b1a0*in1 + b2a0*in2 - a1a0*ou1 - a2a0*ou2;
// push in/out buffers
in2:=in1;
in1:=input;
ou2:=ou1;
ou1:=LastOut;
// return output
result:=LastOut;
end;
{
the process method is overloaded.
use Process(input:single):single;
for per sample processing
use Process(Input:psingle;sampleframes:integer);
for block processing. The input is a pointer to
the start of an array of single which contains
the audio data.
i.e.
RBJFilter.Process(@WaveData[0],256);
}
procedure TRbjEqFilter.Process(Input:psingle;sampleframes:integer);
var
i:integer;
LastOut:single;
begin
for i:=0 to SampleFrames-1 do
begin
// filter
LastOut:= b0a0*(input^)+ b1a0*in1 + b2a0*in2 - a1a0*ou1 - a2a0*ou2;
//LastOut:=input^;
// push in/out buffers
in2:=in1;
in1:=input^;
ou2:=ou1;
ou1:=LastOut;
Out1[i]:=LastOut;
inc(input);
end;
end;
end.
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;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
