|
|
|
|
 |
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 |
|
|
|
|
|
|
|
|
303 type filter with saturation
Type : Runge-Kutta Filters
References : Posted by Hans Mikelson
Linked file : filters001.txt
Notes :
I
posted a filter to the Csound mailing list a couple of weeks ago that
has a 303 flavor to it. It basically does wacky distortions to the
sound. I used Runge-Kutta for the diff eq. simulation though which
makes it somewhat sluggish.
This is a CSound score!!
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Another 4-pole lowpass...
Type : 4-pole LP/HP
References : Posted by fuzzpilz [AT] gmx [DOT] net
Notes :
Vaguely based on the Stilson/Smith Moog paper, but going in a rather different direction from others I've seen here.
The parameters are peak frequency and peak magnitude (g below); both
are reasonably accurate for magnitudes above 1. DC gain is 1.
The filter has some undesirable properties - e.g. it's unstable for low
peak freqs if implemented in single precision (haven't been able to
cleanly separate it into biquads or onepoles to see if that helps), and
it responds so strongly to parameter changes that it's not advisable to
update the coefficients much more rarely than, say, every eight samples
during sweeps, which makes it somewhat expensive.
I like the sound, however, and the accuracy is nice to have, since many filters are not very strong in that respect.
I haven't looked at the HP again for a while, but IIRC it had approximately the same good and bad sides.
Code :
double coef[9];
double d[4];
double omega; //peak freq
double g; //peak mag
// calculating coefficients:
double k,p,q,a;
double a0,a1,a2,a3,a4;
k=(4.0*g-3.0)/(g+1.0);
p=1.0-0.25*k;p*=p;
// LP:
a=1.0/(tan(0.5*omega)*(1.0+p));
p=1.0+a;
q=1.0-a;
a0=1.0/(k+p*p*p*p);
a1=4.0*(k+p*p*p*q);
a2=6.0*(k+p*p*q*q);
a3=4.0*(k+p*q*q*q);
a4= (k+q*q*q*q);
p=a0*(k+1.0);
coef[0]=p;
coef[1]=4.0*p;
coef[2]=6.0*p;
coef[3]=4.0*p;
coef[4]=p;
coef[5]=-a1*a0;
coef[6]=-a2*a0;
coef[7]=-a3*a0;
coef[8]=-a4*a0;
// or HP:
a=tan(0.5*omega)/(1.0+p);
p=a+1.0;
q=a-1.0;
a0=1.0/(p*p*p*p+k);
a1=4.0*(p*p*p*q-k);
a2=6.0*(p*p*q*q+k);
a3=4.0*(p*q*q*q-k);
a4= (q*q*q*q+k);
p=a0*(k+1.0);
coef[0]=p;
coef[1]=-4.0*p;
coef[2]=6.0*p;
coef[3]=-4.0*p;
coef[4]=p;
coef[5]=-a1*a0;
coef[6]=-a2*a0;
coef[7]=-a3*a0;
coef[8]=-a4*a0;
// per sample:
out=coef[0]*in+d[0];
d[0]=coef[1]*in+coef[5]*out+d[1];
d[1]=coef[2]*in+coef[6]*out+d[2];
d[2]=coef[3]*in+coef[7]*out+d[3];
d[3]=coef[4]*in+coef[8]*out;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
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 |
|
|
|
|
|
|
|
|
Cool Sounding Lowpass With Decibel Measured Resonance
Type : LP 2-pole resonant tweaked butterworth
References : Posted by daniel_jacob_werner [AT] yaho [DOT] com [DOT] au
Notes :
This
algorithm is a modified version of the tweaked butterworth lowpass
filter by Patrice Tarrabia posted on musicdsp.org's archives. It
calculates the coefficients for a second order IIR filter. The
resonance is specified in decibels above the DC gain. It can be made
suitable to use as a SoundFont 2.0 filter by scaling the output so the
overall gain matches the specification (i.e. if resonance is 6dB then
you should scale the output by -3dB). Note that you can replace the
sqrt(2) values in the standard butterworth highpass algorithm with my
"q =" line of code to get a highpass also. How it works: normally q is
the constant sqrt(2), and this value controls resonance. At sqrt(2)
resonance is 0dB, smaller values increase resonance. By multiplying
sqrt(2) by a power ratio we can specify the resonant gain at the cutoff
frequency. The resonance power ratio is calculated with a standard
formula to convert between decibels and power ratios (the powf
statement...).
Good Luck,
Daniel Werner
http://experimentalscene.com/
Code :
float c, csq, resonance, q, a0, a1, a2, b1, b2;
c = 1.0f / (tanf(pi * (cutoff / samplerate)));
csq = c * c;
resonance = powf(10.0f, -(resonancedB * 0.1f));
q = sqrt(2.0f) * resonance;
a0 = 1.0f / (1.0f + (q * c) + (csq));
a1 = 2.0f * a0;
a2 = a0;
b1 = (2.0f * a0) * (1.0f - csq);
b2 = a0 * (1.0f - (q * c) + csq);
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
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 |
|
|
|
|
|
|
|
|
Digital RIAA equalization filter coefficients
Type : RIAA
References : Posted by Frederick Umminger
Notes :
Use at your own risk. Confirm correctness before using. Don't assume I didn't goof something up.
-Frederick Umminger
Code :
The
"turntable-input software" thread inspired me to generate some
coefficients for a digital RIAA equalization filter. These coefficients
were found by matching the magnitude response of the s-domain transfer
function using some proprietary Matlab scripts. The phase response may
or may not be totally whacked.
The s-domain transfer function is
R3(1+R1*C1*s)(1+R2*C2*s)/(R1(1+R2*C2*s) + R2(1+R1*C1*s) + R3(1+R1*C1*s)(1+R2*C2*s))
where
R1 = 883.3k
R2 = 75k
R3 = 604
C1 = 3.6n
C2 = 1n
This is based on the reference circuit found in http://www.hagtech.com/pdf/riaa.pdf
The coefficients of the digital transfer function b(z^-1)/a(z^-1) in descending powers of z, are:
44.1kHz
b = [ 0.02675918611906 -0.04592084787595 0.01921229297239]
a = [ 1.00000000000000 -0.73845850035973 -0.17951755477430]
error +/- 0.25dB
48kHz
b = [ 0.02675918611906 -0.04592084787595 0.01921229297239]
a = [ 1.00000000000000 -0.73845850035973 -0.17951755477430]
error +/- 0.15dB
88.2kHz
b = [ 0.04872204977233 -0.09076930609195 0.04202280710877]
a = [ 1.00000000000000 -0.85197860443215 -0.10921171201431]
error +/- 0.01dB
96kHz
b = [ 0.05265477122714 -0.09864197097385 0.04596474352090 ]
a = [ 1.00000000000000 -0.85835597216218 -0.10600020417219 ]
error +/- 0.006dB
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Direct form II
Type : generic
References : Posted by Fuzzpilz
Notes :
I've
noticed there's no code for direct form II filters in general here,
though probably many of the filter examples use it. I haven't looked at
them all to verify that, but there certainly doesn't seem to be a
snippet describing this.
This is a simple direct form II implementation of a k-pole, k-zero
filter. It's a little faster than (a naive, real-time implementation
of) direct form I, as well as more numerically accurate.
Code :
Direct form I pseudocode:
y[n] = a[0]*x[n] + a[1]*x[n-1] + .. + a[k]*x[n-k]
- b[1]*y[n-1] - .. - b[k]*y[n-k];
Simple equivalent direct form II pseudocode:
y[n] = a[0]*x[n] + d[0];
d[0] = a[1]*x[n] - b[1]*y[n] + d[1];
d[1] = a[2]*x[n] - b[2]*y[n] + d[2];
.
.
d[k-2] = a[k-1]*x[n] - b[k-1]*y[n] + d[k-1];
d[k-1] = a[k]*x[n] - b[k]*y[n];
For example, a biquad:
out = a0*in + a1*h0 + a2*h1 - b1*h2 - b2*h3;
h1 = h0;
h0 = in;
h3 = h2;
h2 = out;
becomes
out = a0*in + d0;
d0 = a1*in - b1*out + d1;
d1 = a2*in - b2*out;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Formant filter
References : Posted by Alex
Code :
/*
Public source code by alex@smartelectronix.com
Simple example of implementation of formant filter
Vowelnum can be 0,1,2,3,4 <=> A,E,I,O,U
Good for spectral rich input like saw or square
*/
//-------------------------------------------------------------VOWEL COEFFICIENTS
const double coeff[5][11]= {
{ 8.11044e-06,
8.943665402, -36.83889529, 92.01697887, -154.337906, 181.6233289,
-151.8651235, 89.09614114, -35.10298511, 8.388101016, -0.923313471 ///A
},
{4.36215e-06,
8.90438318, -36.55179099, 91.05750846, -152.422234, 179.1170248, ///E
-149.6496211,87.78352223, -34.60687431, 8.282228154, -0.914150747
},
{ 3.33819e-06,
8.893102966, -36.49532826, 90.96543286, -152.4545478, 179.4835618,
-150.315433, 88.43409371, -34.98612086, 8.407803364, -0.932568035 ///I
},
{1.13572e-06,
8.994734087, -37.2084849, 93.22900521, -156.6929844, 184.596544, ///O
-154.3755513, 90.49663749, -35.58964535, 8.478996281, -0.929252233
},
{4.09431e-07,
8.997322763, -37.20218544, 93.11385476, -156.2530937, 183.7080141, ///U
-153.2631681, 89.59539726, -35.12454591, 8.338655623, -0.910251753
}
};
//---------------------------------------------------------------------------------
static double memory[10]={0,0,0,0,0,0,0,0,0,0};
//---------------------------------------------------------------------------------
float formant_filter(float *in, int vowelnum)
{
res= (float) ( coeff[vowelnum][0] *in +
coeff[vowelnum][1] *memory[0] +
coeff[vowelnum][2] *memory[1] +
coeff[vowelnum][3] *memory[2] +
coeff[vowelnum][4] *memory[3] +
coeff[vowelnum][5] *memory[4] +
coeff[vowelnum][6] *memory[5] +
coeff[vowelnum][7] *memory[6] +
coeff[vowelnum][8] *memory[7] +
coeff[vowelnum][9] *memory[8] +
coeff[vowelnum][10] *memory[9] );
memory[9]= memory[8];
memory[8]= memory[7];
memory[7]= memory[6];
memory[6]= memory[5];
memory[5]= memory[4];
memory[4]= memory[3];
memory[3]= memory[2];
memory[2]= memory[1];
memory[1]= memory[0];
memory[0]=(double) res;
return res;
}
8 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
frequency warped FIR lattice
Type : FIR using allpass chain
References : Posted by mail[AT]mutagene[DOT]net
Notes :
Not
at all optimized and pretty hungry in terms of arrays and overhead
(function requires two arrays containing lattice filter's internal
state and ouputs to another two arrays with their next states). In this
implementation I think you'll have to juggle taps1/newtaps in your
processing loop, alternating between one set of arrays and the other
for which to send to wfirlattice).
A frequency-warped lattice filter is just a lattice filter where every
delay has been replaced with an allpass filter. By adjusting the
allpass filters, the frequency response of the filter can be adjusted
(e.g., design an FIR that approximates some filter. Play with with
warping coefficient to "sweep" the FIR up and down without changing any
other coefficients). Much more on warped filters can be found on Aki
Harma's website ( http://www.acoustics.hut.fi/~aqi/ )
Code :
float
wfirlattice(float input, float *taps1, float *taps2, float *reflcof,
float lambda, float *newtaps1, float *newtaps2, int P)
// input is filter input
// taps1,taps2 are previous filter states (init to 0)
// reflcof are reflection coefficients. abs(reflcof) < 1 for stable filter
// lamba is warping (0 = no warping, 0.75 is close to bark scale at 44.1 kHz)
// newtaps1, newtaps2 are new filter states
// P is the order of the filter
{
float forward;
float topline;
forward = input;
topline = forward;
for (int i=0;i<P;i++)
{
newtaps2[i] = topline;
newtaps1[i] = float(lambda)*(-topline + taps1[i]) + taps2[i];
topline = newtaps1[i]+forward*(reflcof[i]);
forward += newtaps1[i]*(reflcof[i]);
taps1[i]=newtaps1[i];
taps2[i]=newtaps2[i];
}
return forward;
}
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Karlsen
Type : 24-dB (4-pole) lowpass
References : Posted by Best Regards,Ove Karlsen
Notes :
There's really not much voodoo going on in the filter itself, it's a simple as possible:
pole1 = (in * frequency) + (pole1 * (1 - frequency));
Most of you can probably understand that math, it's very similar to how an analog condenser works.
Although, I did have to do some JuJu to add resonance to it.
While studing the other filters, I found that the feedback phase is very important to how the overall
resonance level will be, and so I made a dynamic feedback path, and constant Q approximation by manipulation
of the feedback phase.
A bonus with this filter, is that you can "overdrive" it... Try high input levels..
Code :
// Karlsen 24dB Filter by Ove Karlsen / Synergy-7 in the year 2003.
// b_f = frequency 0..1
// b_q = resonance 0..50
// b_in = input
// to do bandpass, subtract poles from eachother, highpass subtract with input.
float b_inSH = b_in // before the while statement.
while (b_oversample < 2) { //2x oversampling (@44.1khz)
float prevfp;
prevfp = b_fp;
if (prevfp > 1) {prevfp = 1;} // Q-limiter
b_fp = (b_fp * 0.418) + ((b_q * pole4) * 0.582); // dynamic feedback
float intfp;
intfp = (b_fp * 0.36) + (prevfp * 0.64); // feedback phase
b_in = b_inSH - intfp; // inverted feedback
pole1 = (b_in * b_f) + (pole1 * (1 - b_f)); // pole 1
if (pole1 > 1) {pole1 = 1;} else if (pole1 < -1) {pole1 = -1;} // pole 1 clipping
pole2 = (pole1 * b_f) + (pole2 * (1 - b_f)); // pole 2
pole3 = (pole2 * b_f) + (pole3 * (1 - b_f)); // pole 3
pole4 = (pole3 * b_f) + (pole4 * (1 - b_f)); // pole 4
b_oversample++;
}
lowpassout = b_in;
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Lowpass filter for parameter edge filtering
References : Olli Niemitalo
Linked file : filter001.gif
Notes :
use this filter to smooth sudden parameter changes
(see linkfile!)
Code :
/* - Three one-poles combined in parallel
* - Output stays within input limits
* - 18 dB/oct (approx) frequency response rolloff
* - Quite fast, 2x3 parallel multiplications/sample, no internal buffers
* - Time-scalable, allowing use with different samplerates
* - Impulse and edge responses have continuous differential
* - Requires high internal numerical precision
*/
{
/* Parameters */
// Number of samples from start of edge to halfway to new value
const double scale = 100;
// 0 < Smoothness < 1. High is better, but may cause precision problems
const double smoothness = 0.999;
/* Precalc variables */
double a = 1.0-(2.4/scale); // Could also be set directly
double b = smoothness; // -"-
double acoef = a;
double bcoef = a*b;
double ccoef = a*b*b;
double mastergain = 1.0 / (-1.0/(log(a)+2.0*log(b))+2.0/
(log(a)+log(b))-1.0/log(a));
double again = mastergain;
double bgain = mastergain * (log(a*b*b)*(log(a)-log(a*b)) /
((log(a*b*b)-log(a*b))*log(a*b))
- log(a)/log(a*b));
double cgain = mastergain * (-(log(a)-log(a*b)) /
(log(a*b*b)-log(a*b)));
/* Runtime variables */
long streamofs;
double areg = 0;
double breg = 0;
double creg = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++)
{
/* Update filters */
areg = acoef * areg + fromstream [streamofs];
breg = bcoef * breg + fromstream [streamofs];
creg = ccoef * creg + fromstream [streamofs];
/* Combine filters in parallel */
long temp = again * areg
+ bgain * breg
+ cgain * creg;
/* Check clipping */
if (temp > 32767)
{
temp = 32767;
}
else if (temp < -32768)
{
temp = -32768;
}
/* Store new value */
tostream [streamofs] = temp;
}
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
LP and HP filter
Type : biquad, tweaked butterworth
References : Posted by Patrice Tarrabia
Code :
r = rez amount, from sqrt(2) to ~ 0.1
f = cutoff frequency
(from ~0 Hz to SampleRate/2 - though many
synths seem to filter only up to SampleRate/4)
The filter algo:
out(n) = a1 * in + a2 * in(n-1) + a3 * in(n-2) - b1*out(n-1) - b2*out(n-2)
Lowpass:
c = 1.0 / tan(pi * f / sample_rate);
a1 = 1.0 / ( 1.0 + r * c + c * c);
a2 = 2* a1;
a3 = a1;
b1 = 2.0 * ( 1.0 - c*c) * a1;
b2 = ( 1.0 - r * c + c * c) * a1;
Hipass:
c = tan(pi * f / sample_rate);
a1 = 1.0 / ( 1.0 + r * c + c * c);
a2 = -2*a1;
a3 = a1;
b1 = 2.0 * ( c*c - 1.0) * a1;
b2 = ( 1.0 - r * c + c * c) * a1;
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper.
Notes :
Digital approximation of Moog VCF. Fairly easy to calculate coefficients, fairly easy to process algorithm, good sound.
Code :
//Init
cutoff = cutoff freq in Hz
fs = sampling frequency //(e.g. 44100Hz)
res = resonance [0 - 1] //(minimum - maximum)
f = 2 * cutoff / fs; //[0 - 1]
k = 3.6*f - 1.6*f*f -1; //(Empirical tunning)
p = (k+1)*0.5;
scale = e^((1-p)*1.386249;
r = res*scale;
y4 = output;
y1=y2=y3=y4=oldx=oldy1=oldy2=oldy3=0;
//Loop
//--Inverted feed back for corner peaking
x = input - r*y4;
//Four cascaded onepole filters (bilinear transform)
y1=x*p + oldx*p - k*y1;
y2=y1*p+oldy1*p - k*y2;
y3=y2*p+oldy2*p - k*y3;
y4=y3*p+oldy3*p - k*y4;
//Clipper band limited sigmoid
y4 = y4 - (y4^3)/6;
oldx = x;
oldy1 = y1;
oldy2 = y2;
oldy3 = y3;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF, variation 1
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper., Paul Kellett version
Notes :
The
second "q =" line previously used exp() - I'm not sure if what I've
done is any faster, but this line needs playing with anyway as it
controls which frequencies will self-oscillate. I
think it could be tweaked to sound better than it currently does.
Highpass / Bandpass :
They are only 6dB/oct, but still seem musically useful - the 'fruity' sound of the 24dB/oct lowpass is retained.
Code :
// Moog 24 dB/oct resonant lowpass VCF
// References: CSound source code, Stilson/Smith CCRMA paper.
// Modified by paul.kellett@maxim.abel.co.uk July 2000
float f, p, q; //filter coefficients
float b0, b1, b2, b3, b4; //filter buffers (beware denormals!)
float t1, t2; //temporary buffers
// Set coefficients given frequency & resonance [0.0...1.0]
q = 1.0f - frequency;
p = frequency + 0.8f * frequency * q;
f = p + p - 1.0f;
q = resonance * (1.0f + 0.5f * q * (1.0f - q + 5.6f * q * q));
// Filter (in [-1.0...+1.0])
in -= q * b4; //feedback
t1 = b1; b1 = (in + b0) * p - b1 * f;
t2 = b2; b2 = (b1 + t1) * p - b2 * f;
t1 = b3; b3 = (b2 + t2) * p - b3 * f;
b4 = (b3 + t1) * p - b4 * f;
b4 = b4 - b4 * b4 * b4 * 0.166667f; //clipping
b0 = in;
// Lowpass output: b4
// Highpass output: in - b4;
// Bandpass output: 3.0f * (b3 - b4);
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Moog VCF, variation 2
Type : 24db resonant lowpass
References : CSound source code, Stilson/Smith CCRMA paper., Timo Tossavainen (?) version
Notes :
in[x] and out[x] are member variables, init to 0.0 the controls:
fc = cutoff, nearly linear [0,1] -> [0, fs/2]
res = resonance [0, 4] -> [no resonance, self-oscillation]
Code :
Tdouble MoogVCF::run(double input, double fc, double res)
{
double f = fc * 1.16;
double fb = res * (1.0 - 0.15 * f * f);
input -= out4 * fb;
input *= 0.35013 * (f*f)*(f*f);
out1 = input + 0.3 * in1 + (1 - f) * out1; // Pole 1
in1 = input;
out2 = out1 + 0.3 * in2 + (1 - f) * out2; // Pole 2
in2 = out1;
out3 = out2 + 0.3 * in3 + (1 - f) * out3; // Pole 3
in3 = out2;
out4 = out3 + 0.3 * in4 + (1 - f) * out4; // Pole 4
in4 = out3;
return out4;
}
6 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Notch filter
Type : 2 poles 2 zeros IIR
References : Posted by Olli Niemitalo
Notes :
Creates
a muted spot in the spectrum with adjustable steepness. A complex
conjugate pair of zeros on the z- plane unit circle and neutralizing
poles approaching at the same angles from inside the unit circle.
Code :
Parameters:
0 =< freq =< samplerate/2
0 =< q < 1 (The higher, the narrower)
AlgoAlgo=double pi = 3.141592654;
double sqrt2 = sqrt(2.0);
double freq = 2050; // Change! (zero & pole angle)
double q = 0.4; // Change! (pole magnitude)
double z1x = cos(2*pi*freq/samplerate);
double a0a2 = (1-q)*(1-q)/(2*(fabs(z1x)+1)) + q;
double a1 = -2*z1x*a0a2;
double b1 = -2*z1x*q;
double b2 = q*q;
double reg0, reg1, reg2;
unsigned int streamofs;
reg1 = 0;
reg2 = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++)
{
reg0 = a0a2 * ((double)fromstream[streamofs]
+ fromstream[streamofs+2])
+ a1 * fromstream[streamofs+1]
- b1 * reg1
- b2 * reg2;
reg2 = reg1;
reg1 = reg0;
int temp = reg0;
/* Check clipping */
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
/* Store new value */
tostream[streamofs] = temp;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One pole LP and HP
References : Posted by Bram
Code :
LP:
recursion: tmp = (1-p)*in + p*tmp with output = tmp
coefficient: p = (2-cos(x)) - sqrt((2-cos(x))^2 - 1) with x = 2*pi*cutoff/samplerate
coeficient approximation: p = (1 - 2*cutoff/samplerate)^2
HP:
recursion: tmp = (p-1)*in - p*tmp with output = tmp
coefficient: p = (2+cos(x)) - sqrt((2+cos(x))^2 - 1) with x = 2*pi*cutoff/samplerate
coeficient approximation: p = (2*cutoff/samplerate)^2
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One pole, one zero LP/HP
References : Posted by mistert[AT]inwind[DOT]it
Code :
void SetLPF(float fCut, float fSampling)
{
float w = 2.0 * fSampling;
float Norm;
fCut *= 2.0F * PI;
Norm = 1.0 / (fCut + w);
b1 = (w - fCut) * Norm;
a0 = a1 = fCut * Norm;
}
void SetHPF(float fCut, float fSampling)
{
float w = 2.0 * fSampling;
float Norm;
fCut *= 2.0F * PI;
Norm = 1.0 / (fCut + w);
a0 = w * Norm;
a1 = -a0;
b1 = (w - fCut) * Norm;
}
Where
out[n] = in[n]*a0 + in[n-1]*a1 + out[n-1]*b1;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
One zero, LP/HP
References : Posted by Bram
Notes :
LP is only 'valid' for cutoffs > samplerate/4
HP is only 'valid' for cutoffs < samplerate/4
Code :
theta = cutoff*2*pi / samplerate
LP:
H(z) = (1+p*z^(-1)) / (1+p)
out[i] = 1/(1+p) * in[i] + p/(1+p) * in[i-1];
p = (1-2*cos(theta)) - sqrt((1-2*cos(theta))^2 - 1)
Pi/2 < theta < Pi
HP:
H(z) = (1-p*z^(-1)) / (1+p)
out[i] = 1/(1+p) * in[i] - p/(1+p) * in[i-1];
p = (1+2*cos(theta)) - sqrt((1+2*cos(theta))^2 - 1)
0 < theta < Pi/2
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Peak/Notch filter
Type : peak/notch
References : Posted by tobybear[AT]web[DOT]de
Notes :
// Peak/Notch filter
// I don't know anymore where this came from, just found it on
// my hard drive :-)
// Seems to be a peak/notch filter with adjustable slope
// steepness, though slope gets rather wide the lower the
// frequency is.
// "cut" and "steep" range is from 0..1
// Try to feed it with white noise, then the peak output does
// rather well eliminate all other frequencies except the given
// frequency in higher frequency ranges.
Code :
var f,r:single;
outp,outp1,outp2:single; // init these with 0!
const p4=1.0e-24; // Pentium 4 denormal problem elimination
function PeakNotch(inp,cut,steep:single;ftype:integer):single;
begin
r:=steep*0.99609375;
f:=cos(pi*cut);
a0:=(1-r)*sqrt(r*(r-4*(f*f)+2)+1);
b1:=2*f*r;
b2:=-(r*r);
outp:=a0*inp+b1*outp1+b2*outp2+p4;
outp2:=outp1;
outp1:=outp;
if ftype=0 then
result:=outp //peak
else
result:=inp-outp; //notch
end;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Phase equalization
Type : Allpass
References : Posted by Uli the Grasso
Notes :
The
idea is simple: One can equalize the phase response of a system, for
example of a loudspeaker, by approximating its phase response by an FIR
filter and then turn around the coefficients of the filter. At
http://grassomusic.de/english/phaseeq.htm you find more info and an
Octave script.
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Prewarping
Type : explanation
References : Posted by robert bristow-johnson (better known as "rbj" )
Notes :
prewarping is simply recognizing the warping that the BLT introduces.
to determine frequency response, we evaluate the digital H(z) at
z=exp(j*w*T) and we evaluate the analog Ha(s) at s=j*W . the following
will confirm the jw to unit circle mapping and will show exactly what the
mapping is (this is the same stuff in the textbooks):
the BLT says: s = (2/T) * (z-1)/(z+1)
substituting: s = j*W = (2/T) * (exp(j*w*T) - 1) / (exp(j*w*T) + 1)
j*W = (2/T) * (exp(j*w*T/2) - exp(-j*w*T/2)) / (exp(j*w*T/2) + exp(-j*w*T/2))
= (2/T) * (j*2*sin(w*T/2)) / (2*cos(w*T/2))
= j * (2/T) * tan(w*T/2)
or
analog W = (2/T) * tan(w*T/2)
so when the real input frequency is w, the digital filter will behave with
the same amplitude gain and phase shift as the analog filter will have at a
hypothetical frequency of W. as w*T approaches pi (Nyquist) the digital
filter behaves as the analog filter does as W -> inf. for each degree of
freedom that you have in your design equations, you can adjust the analog
design frequency to be just right so that when the deterministic BLT
warping does its thing, the resultant warped frequency comes out just
right. for a simple LPF, you have only one degree of freedom, the cutoff
frequency. you can precompensate it so that the true cutoff comes out
right but that is it, above the cutoff, you will see that the LPF dives
down to -inf dB faster than an equivalent analog at the same frequencies.
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Resonant filter
References : Posted by Paul Kellett
Notes :
This filter consists of two first order low-pass filters in
series, with some of the difference between the two filter
outputs fed back to give a resonant peak.
You can use more filter stages for a steeper cutoff but the
stability criteria get more complicated if the extra stages
are within the feedback loop.
Code :
//set feedback amount given f and q between 0 and 1
fb = q + q/(1.0 - f);
//for each sample...
buf0 = buf0 + f * (in - buf0 + fb * (buf0 - buf1));
buf1 = buf1 + f * (buf0 - buf1);
out = buf1;
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Resonant IIR lowpass (12dB/oct)
Type : Resonant IIR lowpass (12dB/oct)
References : Posted by Olli Niemitalo
Notes :
Hard to calculate coefficients, easy to process algorithm
Code :
resofreq = pole frequency
amp = magnitude at pole frequency (approx)
double pi = 3.141592654;
/* Parameters. Change these! */
double resofreq = 5000;
double amp = 1.0;
DOUBLEWORD streamofs;
double w = 2.0*pi*resofreq/samplerate; // Pole angle
double q = 1.0-w/(2.0*(amp+0.5/(1.0+w))+w-2.0); // Pole magnitude
double r = q*q;
double c = r+1.0-2.0*cos(w)*q;
double vibrapos = 0;
double vibraspeed = 0;
/* Main loop */
for (streamofs = 0; streamofs < streamsize; streamofs++) {
/* Accelerate vibra by signal-vibra, multiplied by lowpasscutoff */
vibraspeed += (fromstream[streamofs] - vibrapos) * c;
/* Add velocity to vibra's position */
vibrapos += vibraspeed;
/* Attenuate/amplify vibra's velocity by resonance */
vibraspeed *= r;
/* Check clipping */
temp = vibrapos;
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
/* Store new value */
tostream[streamofs] = temp;
}
1 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Reverb Filter Generator
Type : FIR
References : Posted by Stephen McGovern
Notes :
This
is a MATLAB function that makes a rough calculation of a room's impulse
response. The output can then be convolved with an audio clip to
produce good and realistic sounding reverb. I have written a paper
discussing the theory used by this algorithm. It is available at
http://stevem.us/rir.html.
NOTES:
1) Large values of N will use large amounts of memory.
2) The output is normalized to the largest value of the
output.
Code :
function [h]=rir(fs, mic, n, r, rm, src);
%RIR Room Impulse Response.
% [h] = RIR(FS, MIC, N, R, RM, SRC) calculates the impulse % response
% of a room.
%
% FS = sample rate.
% MIC = row vector, MIC=[X Y Z], giving the x,y,z % % coordinates of
% the microphone.
% N = The program will account for (2*N+1)^3 virtual % sources
% R = reflection coefficient for the walls, 0<R<1.
% RM = row vector, RM=[X Y Z], giving the dimensions % of the room.
% SRC = row vector, SRC=[X Y Z], giving the x,y,z % coordinates of
% the sound source.
%
% NOTES:
%
% 1) To implement this filter, you will need to do a fast
% convolution. The program FCONV.m will do this. It is % available
% at: http://stevem.us/code/fconv.m
% 2) All distances are in meters.
% 3) If this is your first time running this program, set % N equal
% to 10 or less, and R=0.9.
% 4) I've written an article discussing the theory behind % this
% algorithm. It can be found at
% http://stevem.us/rir.html.
%
%
%Version 1.0
%Coded by: Stephen G. McGovern, 2003.
%The comments below refer to equations in my paper.
nn=[-n:1:n]; % Index for the
% sequence
rms= nn+0.5-0.5*(-1).^nn; % Part of equations
% 2,3,& 4
srcs=(-1).^(nn); % part of equations
% 2,3,& 4
xi=[srcs*src(1)+rms*rm(1)-mic(1)]; % Equation 2
yj=[srcs*src(2)+rms*rm(2)-mic(2)]; % Equation 3
zk=[srcs*src(3)+rms*rm(3)-mic(3)]; % Equation 4
[i,j,k]=meshgrid(xi,yj,zk); % convert vectors to
% 3D matrices
d=sqrt(i.^2+j.^2+k.^2); % Equation 5
time=d./343; % Similar to equation
% 6
time=round(time*fs); % Quantized delay time
b=1./(4*pi*((d).^(2))); % Equation 8
[e,f,g] = meshgrid(nn, nn, nn); % convert vectors to
% 3D matrices
c=r.^(abs(e)+abs(f)+abs(g)); % Equation 9
e= b.*c; % Equation 10
h=full(sparse(time(:),1,e(:))); % Equivalent to
% equation 11
h=h/max(abs(h)); % Normalize h
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State variable
Type : 12db resonant low, high or bandpass
References : Effect Deisgn Part 1, Jon Dattorro, J. Audio Eng. Soc., Vol 45, No. 9, 1997 September
Notes :
Digital approximation of Chamberlin two-pole low pass. Easy to calculate coefficients, easy to process algorithm.
Code :
cutoff = cutoff freq in Hz
fs = sampling frequency //(e.g. 44100Hz)
f = 2 sin (pi * cutoff / fs) //[approximately]
q = resonance/bandwidth [0 < q <= 1] most res: q=1, less: q=0
low = lowpass output
high = highpass output
band = bandpass output
notch = notch output
scale = q
low=high=band=0;
//--beginloop
low = low + f * band;
high = scale * input - low - q*band;
band = f * high + band;
notch = high + low;
//--endloop
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State Variable Filter (Chamberlin version)
References : Hal Chamberlin, "Musical Applications of Microprocessors," 2nd Ed, Hayden Book Company 1985. pp 490-492.
Code :
//Input/Output
I - input sample
L - lowpass output sample
B - bandpass output sample
H - highpass output sample
N - notch output sample
F1 - Frequency control parameter
Q1 - Q control parameter
D1 - delay associated with bandpass output
D2 - delay associated with low-pass output
// parameters:
Q1 = 1/Q
// where Q1 goes from 2 to 0, ie Q goes from .5 to infinity
// simple frequency tuning with error towards nyquist
// F is the filter's center frequency, and Fs is the sampling rate
F1 = 2*pi*F/Fs
// ideal tuning:
F1 = 2 * sin( pi * F / Fs )
// algorithm
// loop
L = D2 + F1 * D1
H = I - L - Q1*D1
B = F1 * H + D1
N = H + L
// store delays
D1 = B
D2 = L
// outputs
L,H,B,N
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
State Variable Filter (Double Sampled, Stable)
Type : 2 Pole Low, High, Band, Notch and Peaking
References : Posted by Andrew Simper
Notes :
Thanks to Laurent de Soras for the stability limit
and Steffan Diedrichsen for the correct notch output.
Code :
input = input buffer;
output = output buffer;
fs = sampling frequency;
fc = cutoff frequency normally something like:
440.0*pow(2.0, (midi_note - 69.0)/12.0);
res = resonance 0 to 1;
drive = internal distortion 0 to 0.1
freq = 2.0*sin(PI*MIN(0.25, fc/(fs*2))); // the fs*2 is because it's double sampled
damp = MIN(2.0*(1.0 - pow(res, 0.25)), MIN(2.0, 2.0/freq - freq*0.5));
notch = notch output
low = low pass output
high = high pass output
band = band pass output
peak = peaking output = low - high
--
double sampled svf loop:
for (i=0; i<numSamples; i++)
{
in = input[i];
notch = in - damp*band;
low = low + freq*band;
high = notch - low;
band = freq*high + band - drive*band*band*band;
out = 0.5*(notch or low or high or band or peak);
notch = in - damp*band;
low = low + freq*band;
high = notch - low;
band = freq*high + band - drive*band*band*band;
out += 0.5*(same out as above);
output[i] = out;
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Stilson's Moog filter code
Type : 4-pole LP, with fruity BP/HP
References : Posted by DFL
Notes :
Mind your p's and Q's...
This code was borrowed from Tim Stilson, and rewritten by me into a pd extern (moog~) available here:
http://www-ccrma.stanford.edu/~dfl/pd/index.htm
I ripped out the essential code and pasted it here...
Code :
WARNING: messy code follows ;)
// table to fixup Q in order to remain constant for various pole
frequencies, from Tim Stilson's code @ CCRMA (also in CLM distribution)
static float gaintable[199] = { 0.999969, 0.990082, 0.980347, 0.970764,
0.961304, 0.951996, 0.94281, 0.933777, 0.924866, 0.916077, 0.90741,
0.898865, 0.89044
2, 0.882141 , 0.873962, 0.865906, 0.857941, 0.850067, 0.842346,
0.834686, 0.827148, 0.819733, 0.812378, 0.805145, 0.798004, 0.790955,
0.783997, 0.77713, 0.77
0355, 0.763672, 0.75708 , 0.75058, 0.744141, 0.737793, 0.731537,
0.725342, 0.719238, 0.713196, 0.707245, 0.701355, 0.695557, 0.689819,
0.684174, 0.678558, 0.
673035, 0.667572, 0.66217, 0.65686, 0.651581, 0.646393, 0.641235,
0.636169, 0.631134, 0.62619, 0.621277, 0.616425, 0.611633, 0.606903,
0.602234, 0.597626, 0.
593048, 0.588531, 0.584045, 0.579651, 0.575287 , 0.570953, 0.566681,
0.562469, 0.558289, 0.554169, 0.550079, 0.546051, 0.542053, 0.538116,
0.53421, 0.530334,
0.52652, 0.522736, 0.518982, 0.515289, 0.511627, 0.507996 , 0.504425,
0.500885, 0.497375, 0.493896, 0.490448, 0.487061, 0.483704, 0.480377,
0.477081, 0.4738
16, 0.470581, 0.467377, 0.464203, 0.46109, 0.457977, 0.454926,
0.451874, 0.448883, 0.445892, 0.442932, 0.440033, 0.437134, 0.434265,
0.431427, 0.428619, 0.42
5842, 0.423096, 0.42038, 0.417664, 0.415009, 0.412354, 0.409729,
0.407135, 0.404572, 0.402008, 0.399506, 0.397003, 0.394501, 0.392059,
0.389618, 0.387207, 0.
384827, 0.382477, 0.380127, 0.377808, 0.375488, 0.37323, 0.370972,
0.368713, 0.366516, 0.364319, 0.362122, 0.359985, 0.357849, 0.355713,
0.353607, 0.351532,
0.349457, 0.347412, 0.345398, 0.343384, 0.34137, 0.339417, 0.337463,
0.33551, 0.333588, 0.331665, 0.329773, 0.327911, 0.32605, 0.324188,
0.322357, 0.320557,
0.318756, 0.316986, 0.315216, 0.313446, 0.311707, 0.309998, 0.308289,
0.30658, 0.304901, 0.303223, 0.301575, 0.299927, 0.298309, 0.296692,
0.295074, 0.293488
, 0.291931, 0.290375, 0.288818, 0.287262, 0.285736, 0.284241, 0.282715,
0.28125, 0.279755, 0.27829, 0.276825, 0.275391, 0.273956, 0.272552,
0.271118, 0.26974
5, 0.268341, 0.266968, 0.265594, 0.264252, 0.262909, 0.261566, 0.260223, 0.258911, 0.257599, 0.256317, 0.255035, 0.25375 };
static inline float saturate( float input ) { //clamp without branching
#define _limit 0.95
float x1 = fabsf( input + _limit );
float x2 = fabsf( input - _limit );
return 0.5 * (x1 - x2);
}
static inline float crossfade( float amount, float a, float b ) {
return (1-amount)*a + amount*b;
}
//code for setting Q
float ix, ixfrac;
int ixint;
ix = x->p * 99;
ixint = floor( ix );
ixfrac = ix - ixint;
Q = resonance * crossfade( ixfrac, gaintable[ ixint + 99 ], gaintable[ ixint + 100 ] );
//code for setting pole coefficient based on frequency
float fc = 2 * frequency / x->srate;
float x2 = fc*fc;
float x3 = fc*x2;
p = -0.69346 * x3 - 0.59515 * x2 + 3.2937 * fc - 1.0072; //cubic fit by DFL, not 100% accurate but better than nothing...
}
process loop:
float state[4], output; //should be global scope / preserved between calls
int i,pole;
float temp, input;
for ( i=0; i < numSamples; i++ ) {
input = *(in++);
output = 0.25 * ( input - output ); //negative feedback
for( pole = 0; pole < 4; pole++) {
temp = state[pole];
output = saturate( output + p * (output - temp));
state[pole] = output;
output = saturate( output + temp );
}
lowpass = output;
highpass = input - output;
bandpass = 3 * x->state[2] - x->lowpass; //got this one from paul kellet
*out++ = lowpass;
output *= Q; //scale the feedback
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Time domain convolution with O(n^log2(3))
References : Wilfried Welti
Notes :
[Quoted from Wilfrieds mail...]
I found last weekend that it is possible to do convolution in time
domain (no complex numbers, 100% exact result with int) with
O(n^log2(3)) (about O(n^1.58)).
Due to smaller overhead compared to FFT-based convolution, it should be the fastest algorithm for medium sized FIR's.
Though, it's slower as FFT-based convolution for large n.
It's pretty easy:
Let's say we have two finite signals of length 2n, which we want
convolve : A and B. Now we split both signals into parts of size n, so
we get A = A1 + A2, and B = B1 +B2.
Now we can write:
(1) A*B = (A1+A2)*(B1+B2) = A1*B1 + A2*B1 + A1*B2 + A2*B2
where * means convolution.
This we knew already: We can split a convolution into four convolutions of halved size.
Things become interesting when we start shifting blocks in time:
Be z a signal which has the value 1 at x=1 and zero elsewhere.
Convoluting a signal X with z is equivalent to shifting X by one
rightwards. When I define z^n as n-fold convolution of z with itself,
like: z^1 = z, z^2 = z*z, z^0 = z shifted leftwards by 1 = impulse at
x=0, and so on, I can use it to shift signals:
X * z^n means shifting the signal X by the value n rightwards.
X * z^-n means shifting the signal X by the value n leftwards.
Now we look at the following term:
(2) (A1 + A2 * z^-n) * (B1 + B2 * z^-n)
This is a convolution of two blocks of size n: We shift A2 by n leftwards so it completely overlaps A1, then we add them.
We do the same thing with B1 and B2. Then we convolute the two resulting blocks.
now let's transform this term:
(3) (A1 + A2 * z^-n) * (B1 + B2 * z^-n)
= A1*B1 + A1*B2*z^-n + A2*z^-n*B1 + A2*z^ n*B2*z^-n
= A1*B1 + (A1*B2 + A2*B1)*z^-n + A2*B2*z^-2n
(4) (A1 + A2 * z^-n) * (B1 + B2 * z^-n) - A1*B1 - A2*B2*z^-2n
= (A1*B2 + A2*B1)*z^-n
Now we convolute both sides of the equation (4) by z^n:
(5) (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n - A1*B1*z^n - A2*B2*z^-n
= (A1*B2 + A2*B1)
Now we see that the right part of equation (5) appears within equation
(1), so we can replace this appearance by the left part of eq (5).
(6) A*B = (A1+A2)*(B1+B2) = A1*B1 + A2*B1 + A1*B2 + A2*B2
= A1*B1
+ (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n - A1*B1*z^n - A2*B2*z^-n
+ A2*B2
Voila!
We have constructed the convolution of A*B with only three convolutions
of halved size. (Since the convolutions with z^n and z^-n are only
shifts
of blocks with size n, they of course need only n operations for processing :)
This can be used to construct an easy recursive algorithm of Order O(n^log2(3))
Code :
void convolution(value* in1, value* in2, value* out, value* buffer, int size)
{
value* temp1 = buffer;
value* temp2 = buffer + size/2;
int i;
// clear output.
for (i=0; i<size*2; i++) out[i] = 0;
// Break condition for recursion: 1x1 convolution is multiplication.
if (size == 1)
{
out[0] = in1[0] * in2[0];
return;
}
// first calculate (A1 + A2 * z^-n)*(B1 + B2 * z^-n)*z^n
signal_add(in1, in1+size/2, temp1, size/2);
signal_add(in2, in2+size/2, temp2, size/2);
convolution(temp1, temp2, out+size/2, buffer+size, size/2);
// then add A1*B1 and substract A1*B1*z^n
convolution(in1, in2, temp1, buffer+size, size/2);
signal_add_to(out, temp1, size);
signal_sub_from(out+size/2, temp1, size);
// then add A2*B2 and substract A2*B2*z^-n
convolution(in1+size/2, in2+size/2, temp1, buffer+size, size/2);
signal_add_to(out+size, temp1, size);
signal_sub_from(out+size/2, temp1, size);
}
"value" may be a suitable type like int or float.
Parameter "size" is the size of the input signals and must be a power of 2. out and buffer must point to arrays of size 2*n.
Just to be complete, the helper functions:
void signal_add(value* in1, value* in2, value* out, int size)
{
int i;
for (i=0; i<size; i++) out[i] = in1[i] + in2[i];
}
void signal_sub_from(value* out, value* in, int size)
{
int i;
for (i=0; i<size; i++) out[i] -= in[i];
}
void signal_add_to(value* out, value* in, int size)
{
int i;
for (i=0; i<size; i++) out[i] += in[i];
}
3 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Time domain convolution with O(n^log2(3))
References : Posted by Magnus Jonsson
Notes :
[see other code by Wilfried Welti too!]
Code :
void mul_brute(float *r, float *a, float *b, int w)
{
for (int i = 0; i < w+w; i++)
r[i] = 0;
for (int i = 0; i < w; i++)
{
float *rr = r+i;
float ai = a[i];
for (int j = 0; j < w; j++)
rr[j] += ai*b[j];
}
}
// tmp must be of length 2*w
void mul_knuth(float *r, float *a, float *b, int w, float *tmp)
{
if (w < 30)
{
mul_brute(r, a, b, w);
}
else
{
int m = w>>1;
for (int i = 0; i < m; i++)
{
r[i ] = a[m+i]-a[i ];
r[i+m] = b[i ]-b[m+i];
}
mul_knuth(tmp, r , r+m, m, tmp+w);
mul_knuth(r , a , b , m, tmp+w);
mul_knuth(r+w, a+m, b+m, m, tmp+w);
for (int i = 0; i < m; i++)
{
float bla = r[m+i]+r[w+i];
r[m+i] = bla+r[i ]+tmp[i ];
r[w+i] = bla+r[w+m+i]+tmp[i+m];
}
}
}
no comments on this item | add a comment | nofrills version |
|
|
|
|
|
|
|
|
Zoelzer biquad filters
Type : biquad IIR
References : Udo Zoelzer: Digital Audio Signal Processing (John Wiley & Sons, ISBN 0 471 97226 6), Chris Townsend
Notes :
Here's the formulas for the Low Pass, Peaking, and Low Shelf, which should
cover the basics. I tried to convert the formulas so they are little more consistent.
Also, the Zolzer low pass/shelf formulas didn't have adjustable Q, so I added that for
consistency with Roberts formulas as well. I think someone may want to check that I did
it right.
------------ Chris Townsend
I mistranscribed the low shelf cut formulas.
Hopefully this is correct. Thanks to James McCartney for noticing.
------------ Chris Townsend
Code :
omega = 2*PI*frequency/sample_rate
K=tan(omega/2)
Q=Quality Factor
V=gain
LPF: b0 = K^2
b1 = 2*K^2
b2 = K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
peakingEQ:
boost:
b0 = 1 + V*K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - V*K/Q + K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
cut:
b0 = 1 + K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - K/Q + K^2
a0 = 1 + V*K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - V*K/Q + K^2
lowShelf:
boost:
b0 = 1 + sqrt(2*V)*K + V*K^2
b1 = 2*(V*K^2 - 1)
b2 = 1 - sqrt(2*V)*K + V*K^2
a0 = 1 + K/Q + K^2
a1 = 2*(K^2 - 1)
a2 = 1 - K/Q + K^2
cut:
b0 = 1 + K/Q + K^2
b1 = 2*(K^2 - 1)
b2 = 1 - K/Q + K^2
a0 = 1 + sqrt(2*V)*K + V*K^2
a1 = 2*(v*K^2 - 1)
a2 = 1 - sqrt(2*V)*K + V*K^2
2 comment(s) | add a comment | nofrills version |
|
|
|
|
|
|
