Please Notice That This Version Does NOT Contain The Comments You Can Find in The Fancy
Please Notice That This Version Does NOT Contain The Comments You Can Find in The Fancy
Please Notice That This Version Does NOT Contain The Comments You Can Find in The Fancy
Contents
Notes :
Ross Bencina: original source code poster
Olli Niemitalo: UpdateWithCubicInterpolation
Code :
//this code is meant as an EXAMPLE
/*
members are:
float phase;
int TableSize;
float sampleRate;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
return table[i] ;
}
phase += (sampleRate/(float)TableSize)/frequency;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
/*
dtable0[i] = table[i+1] - table[i]; //remember to wrap around!!!
*/
phase += (sampleRate/(float)TableSize)/frequency;
#ifdef FM_OSCILLATOR
if(phase < 0.f)
phase += (float)TableSize;
#endif
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.
Notes :
Implementation in CSound and Sync Modular...
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.
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
Code :
y1 = (max>>1) * atan(k * x/max)
y2 = max * th(x/max)
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.
Linked files
; ORCHESTRA
;---------------------------------------------------------
; Runge-Kutta Filters
; Coded by Hans Mikelson June, 2000
;---------------------------------------------------------
sr = 44100 ; Sample rate
kr = 44100 ; Kontrol rate
ksmps = 1 ; Samples/Kontrol period
nchnls = 2 ; Normal stereo
zakinit 50, 50
;------------------------------------------------------------
; Envelope (Knob twisting simulation)
;------------------------------------------------------------
instr 1
idur = p3 ; Duration
iamp = p4 ; Amplitude
ilps = p5 ; Loops
iofst = p6 ; Offset
itabl = p7 ; Table
ioutch = p8 ; Output channel
iphase = p9 ; Phase
endin
;---------------------------------------------------------
; Runge-Kutta Freaky Filter
;---------------------------------------------------------
instr 7
idur = p3 ; Duration
iamp = p4 ; Amplitude
kfco2 zkr p5 ; Filter cutoff
kq1 zkr p6 ; Q
ih = .001 ; Diff eq step size
ipanl = sqrt(p8) ; Pan left
ipanr = sqrt(1-p8) ; Pan right
ifqc = cpspch(p9) ; Pitch to frequency
kpb zkr p10 ; Pentic bounce frequency
kpa zkr p11 ; Pentic bounce amount
kasym zkr p12 ; Q assymmetry amount
kasep zkr p13 ; Q asymmetry separation
ay init 0
ay1 init 0
ay2 init 0
ay3 init 0
axs init 0
avxs init 0
ax vco 1, ifqc, 2, 1, 1, 1 ; Square wave
; R-K Section 1
afdbk = kq*ay/(1+exp(-ay*3*kasep)*kasym) ; Only oscillate in one
direction
ak11 = ih*((ax-ay1)*kfc-afdbk)
ak21 = ih*((ax-(ay1+.5*ak11))*kfc-afdbk)
ak31 = ih*((ax-(ay1+.5*ak21))*kfc-afdbk)
ak41 = ih*((ax-(ay1+ak31))*kfc-afdbk)
ay1 = ay1+(ak11+2*ak21+2*ak31+ak41)/6
; R-K Section 2
ak12 = ih*((ay1-ay2)*kfc)
ak22 = ih*((ay1-(ay2+.5*ak12))*kfc)
ak32 = ih*((ay1-(ay2+.5*ak22))*kfc)
ak42 = ih*((ay1-(ay2+ak32))*kfc)
ay2 = ay2+(ak12+2*ak22+2*ak32+ak42)/6
; R-K Section 3
ak13 = ih*((ay2-ay3)*kfc+aaxs)
ak23 = ih*((ay2-(ay3+.5*ak13))*kfc+aaxs)
ak33 = ih*((ay2-(ay3+.5*ak23))*kfc+aaxs)
ak43 = ih*((ay2-(ay3+ak33))*kfc+aaxs)
ay3 = ay3+(ak13+2*ak23+2*ak33+ak43)/6
; R-K Section 4
ak14 = ih*((ay3-ay)*kfc)
ak24 = ih*((ay3-(ay+.5*ak14))*kfc)
ak34 = ih*((ay3-(ay+.5*ak24))*kfc)
ak44 = ih*((ay3-(ay+ak34))*kfc)
ay = ay+(ak14+2*ak24+2*ak34+ak44)/6
endin
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 )))));
}
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*/
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)))));
};
Notes :
(see linkfile)
Linked files
Ok, here is how I did it. I'm not 100% sure that everything is correct.
I'll demonstrate it on a square wave instead, although i havent tried it.
You all know what the waveform looks like. It's a constant -1
followed by constant 1 followed by .....
We need to convolve the impulse response with the waveform and sample it at discrete
intervals.
What if we assume that we already know the result of the last sample?
Then we can "move" the previous result backwards by multiplying it with exp(p) since
r(t+1) = r(t)*exp(p)
Now the problem is reduced to integrate only between the last sample and the current
sample, and add that to the result.
some pseudo-code:
while forever
{
result *= exp(pole);
phase += freq;
result += integrate(waveform(phase-freq*t), exp(t*pole), t=0..1);
}
integrate(k*exp(t*pole)) = k*exp(t*pole)/pole
k = square(phase-freq*t)
k*exp(pole)/pole-k/pole = (exp(pole)-1)*k/pole
the only problem left to solve now is the jumps from +1 to -1 and vice versa.
then you detect if phase goes beyond 0.5 or 1. If so, find out exactly
where between the samples this change takes place.
Since I am terribly bad at explaining things, this is probably just a mess to you :)
float getsample()
{
sum *= exp(pole);
phase += freq;
sum += (exp(pole)*((phase-freq)*pole+freq)-(phase*pole+freq))/pole;
sum -= exp(pole*x)-1.0f;
phase -= 1;
}
return sum.real();
}
Alien Wah
Notes :
"I found this algoritm by "playing around" with complex numbers. Please email me your opinions about it.
Paul."
Linked files
/*
Alien-Wah by Nasca Octavian Paul from Tg. Mures, Romania
e-mail: <[email protected]> or <[email protected]>.
*/
/*
The algorithm was found by me by mistake(I was looking for something else);
I called this effect "Alien Wah" because sounds a bit like wahwah, but more strange.
The ideea of this effect is very simple: It is a feedback delay who uses complex
numbers.
If x[] represents the input and y[] is the output, so a simple feedback delay looks
like this:
y[n]=y[n-delay]*fb+x[n]*(1-fb)
alpha is the phase of the number and is controlled by the LFO(Low Frequency
Oscillator).
If the 'delay' parameter is low, the effect sounds more like wah-wah,
but if it is big, the effect will sound very interesting.
The input x[n] has the real part of the samples from the wavefile and the imaginary
part is zero.
The output of this effect is the real part of y[n].
/*
Alien Wah Parameters
#include <complex.h>
#include <fcntl.h>
#include <sys\stat.h>
#include <io.h>
#include <stdio.h>
#include <math.h>
/*
.raw files are raw files (without header), signed 16 bit,mono
*/
#define infile "a.raw" //input file
#define outfile "b.raw" //input file
#define samplerate 44100
#define lfoskipsamples 25 // How many samples are processed before compute the lfo
value again
struct params
{
float freq,startphase,fb;
int delay;
} awparams;
//alien wah internal parameters
struct alienwahinternals
{
complex *delaybuf;
float lfoskip;
long int t;
complex c;
int k;
} awint;
//effect initialisation
void init(float freq,float startphase,float fb,int delay){
awparams.freq=freq;
awparams.startphase=startphase;
awparams.fb=fb/4+0.74;
awparams.delay=(int)(delay/44100.0*samplerate);
if (delay<1) delay=1;
awint.delaybuf=new complex[awparams.delay];
int i;
for (i=0;i<delay;i++) awint.delaybuf[i]=complex(0,0);
awint.lfoskip=freq*2*3.141592653589/samplerate;
awint.t=0;
}
//process buffer
void process()
{
int i;
float lfo,out;
complex outc;
for(i=0;i<bufsize;i++)
{
if (awint.t++%lfoskipsamples==0)
{
lfo=(1+cos(awint.t*awint.lfoskip+awparams.startphase));
awint.c=complex(cos(lfo)*awparams.fb,sin(lfo)*awparams.fb);
};
outc=awint.c*awint.delaybuf[awint.k]+(1-awparams.fb)*buf1[i];
awint.delaybuf[awint.k]=outc;
if ((++awint.k)>=awparams.delay)
awint.k=0;
out=real(outc)*3; //take real part of outc
if (out<-32768) out=-32768;
else if (out>32767) out=32767; //Prevents clipping
buf2[i]=out;
};
}
int main()
{
char f1,f2;
int readed;
long int filereaded=0;
printf("\n");
f1=open(infile,O_RDONLY|O_BINARY);
remove(outfile);
f2=open(outfile,O_BINARY|O_CREAT,S_IWRITE);
long int i;
do
{
readed=read(f1,buf1,bufsize*2);
process();
write(f2,buf2,readed);
printf("%ld bytes \r",filereaded);
filereaded+=readed;
}while (readed==bufsize*2);
delete(awint.delaybuf);
close(f1);
close(f2);
printf("\n\n");
return(0);
}
Type : information
References : Posted by Olli Niemitalo
Linked file : filters002.txt (this linked file is included below)
Linked files
All-Pass Filters
----------------
|
|
| |
| | .
| | | | .
|____|____|____|____|____|____|____.
Instead, the first peak is negative and has a carefully chosen height.
|
| | .
____|____|____|____|____|____.____.
|
|
|
|
|
This can be achieved by adding scaled in(T) to the scaled output of:
B) The sum of odd peaks minus the sum of even peaks (including the
negative peak which we give number zero) must be 1.
Let's get into business now that we know what we want. We call the
amplitude of the negative peak "a" and the amplitude of the first postive
peak "b". The ratio between adjacent positive peaks, the feedback gain, is
denoted by "g".
a + b/(1-g) = 1.
B) Using similar math... The sum of odd peaks is b/(1-g^2). The sum of
even peaks is a + b*g/(1-g^2). So the second constraint is formed:
b/(1-g^2) - (a + b*g/(1-g^2)) = 1.
a = -g
b = 1-g^2
That can be thought of as two filters put in serial and a third one in
parallel with them. The first of the two serial ones is the delay with
feedback. The other "filters" are just different gains.
For the first code line, the frequency response is the usual set peaks:
1
--------------------
1 - g e^(-i w delay)
Adding the mixing on the second code line gives the flat-magnitude
frequency response of the whole allpass system:
| 1/g - g |
| -------------------- - 1/g | = 1
| 1 - g e^(-i w delay) |
(g^2-1) sin(delay w)
atan -------------------------
(g^2+1) cos(delay w) - 2g
I hope this cleared things out - and that there aren't fatal errors! :)
1/g - g 1
H = --------------- - -
1 - g*exp(-iwT) g
1 - g^2 - (1 - g*exp(-iwT))
H = ----------------------------
(1 - g*exp(-iwT)) * g
-g + exp(-iwT)
H = ----------------
1 - g*exp(-iwT)
1 - g*exp(iwT)
H = exp(-iwT) * ---------------
1 - g*exp(-iwT)
|H| = 1
et
g * sin(wT)
arg (H) = -wT - 2 * arctan ---------------
1 - g * cos(wT)
Ce dephasage est le meme que celui trouve par Olli mais a l'avantage de
separer le retard global et le dephasage.
More generally (in fact, maximally generally) you will get an all-pass
response with any transfer function of the form :
where |b| = 1 and of course n should be large enough that your filter is
causal.
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;
return(aligned_ptr);
}
Notes :
This is working code for synthesizing a bandlimited sawtooth waveform. The algorithm is DSF BLIT + leaky
integrator. Includes driver code.
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
Linked files
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double pi = 3.1415926535897932384626433832795029L;
double twopi = 6.2831853071795864769252867665590058L;
typedef struct
{
double x1, y1;
double a, b;
} lowpass_t;
lp->x1 = x;
lp->y1 = y;
return y;
}
typedef struct
{
double phase; /* phase accumulator */
double aNQ; /* attenuation at nyquist */
double curcps; /* current frequency, updated once per cycle */
double curper; /* current period, updated once per cycle */
lowpass_t leaky; /* leaky integrator */
double N; /* # partials */
double a; /* dsf parameter which controls roll-off */
double aN; /* former to the N */
} blit_t;
P2 = b->curper / 2.0;
b->N = 1.0 + floor(P2); /* # of partials incl. dc */
n = 1.0 -
b->aN * cos(Nbeta) -
b->a * (cosbeta - b->aN * cos(Nbeta - beta));
d = b->curper * (1.0 + b->a * (-2.0 * cosbeta + b->a));
return saw;
}
s = 32767.0 * x;
{
int i, L;
double x, cps, cm;
blit_t b;
L = 1000000;
cps *= cm;
}
return 0;
}
Notes :
(see linkfile)
Linked files
///////////////////////////
// Wavetable reverse lookup
// Given a playback rate of the wavetable, what is wavetables index?
//
// rate = f.wavesize/fs e.g. 4096f/44100
// max partials = nyquist/f = wavesize/2rate e.g. 2048/rate
//
// using max partials we could then do a lookup to find the wavetables index
// in a pre-calculated table
//
// however, we could skip max partials, and lookup a table based on a
// function of f (or rate)
//
// the first few midi notes (0 - 9) differ by < 1 so there are duplicates
// values of (int) f.
// therefore, to get an index to our table (that indexes the wavetables)
// we need 2f
//
// to get 2f from rate we multiply by the constant
// 2f = 2.fs/wavesize e.g. 88200/4096
//
// our lookup table will have a length>25087 to cover the midi range
// we'll make it 32768 in length for easy processing
int a,b,n;
float* data;
float* sinetable=new float[4096];
float* datap;
for(b=0;b<4096;b++)
sinetable[b]=sin(TWOPI*(float)b/4096.0f);
int partials;
int partial;
int partialindex,reverseindex,lastnumpartials;
float max,m;
int* reverse;
// sawtooth
data=new float[128*4096];
reverse=new int[32768];
reverseindex=0;
partialindex=0;
lastnumpartials=-1;
for(n=0;n<128;n++)
{
partials=(int)((SPEED*0.5f)/float(TUNING*(float)pow(2,(float)
(n-69)/12.0f))); //(int) NYQUIST/f
if(partials!=lastnumpartials)
{
datap=&data[partialindex*4096];
for(b=0;b<4096;b++)
datap[b]=0.0f; //blank wavetable
for(a=0;a<partials;a++)
{
partial=a+1;
m=cos((float)a*HALFPI/(float)partials);
//gibbs
m*=m; //gibbs
m/=(float)partial;
for(b=0;b<4096;b++)
datap[b]+=m*sinetable[(b*partial)%4096];
}
lastnumpartials=partials;
a=int(2.0f*TUNING*(float)pow(2,(float) (n-
69)/12.0f)); //2f
for(b=reverseindex;b<=a;b++)
reverse[b]=partialindex;
reverseindex=a+1;
partialindex++;
}
}
for(b=reverseindex;b<32768;b++)
reverse[b]=partialindex-1;
max=0.0;
for(b=0;b<4096;b++)
{
if(fabs(*(data+b))>max) //normalise to richest waveform (0)
max=(float)fabs(*(data+b));
}
for(b=0;b<4096*partialindex;b++)
{
*(data+b)/=max;
}
//ar.Write(data,4096*partialindex*sizeof(float));
//ar.Write(reverse,32768*sizeof(int));
delete [] data;
delete [] reverse;
}
// end sawtooth
// parabola
data=new float[128*4096];
reverse=new int[32768];
reverseindex=0;
partialindex=0;
lastnumpartials=-1;
float sign;
for(n=0;n<128;n++)
{
partials=(int)((SPEED*0.5f)/float(TUNING*(float)pow(2,(float)
(n-69)/12.0f)));
if(partials!=lastnumpartials)
{
datap=&data[partialindex*4096];
for(b=0;b<4096;b++)
datap[b]=PI*PI/3.0f;
sign=-1.0f;
for(a=0;a<partials;a++)
{
partial=a+1;
m=cos((float)a*HALFPI/(float)partials);
//gibbs
m*=m; //gibbs
m/=(float)(partial*partial);
m*=4.0f*sign;
for(b=0;b<4096;b++)
for(b=reverseindex;b<32768;b++)
reverse[b]=partialindex-1;
max=0.0;
for(b=0;b<4096;b++)
{
if(fabs(*(data+b))>max) //normalise to richest waveform (0)
max=(float)fabs(*(data+b));
}
max*=0.5;
for(b=0;b<4096*partialindex;b++)
{
*(data+b)/=max;
*(data+b)-=1.0f;
}
//ar.Write(data,4096*partialindex*sizeof(float));
//ar.Write(reverse,32768*sizeof(int));
delete [] data;
delete [] reverse;
}
// end parabola
/////////////////////////////////////////////////////////////////////////
// An example of playback of a sawtooth wave
// This is not optimised for easy reading
// When optimising you'll need to get this in assembly (especially those
// float to int conversions)
/////////////////////////////////////////////////////////
float index;
float rate;
int wavetableindex;
float ratetofloatfactor;
float* wavetable;
float m;
float temprate;
while(length--)
{
currentsample=(int) index;
nextsample=(currentsample+1) & WAVETABLE_MASK;
m=index-(float) currentsample; //fractional part
*buffer++=(1.0f-m)*wavetable[currentsample]+m*wavetable[nextsample];
//linear interpolation
rate*=slide; //slide coeffecient if required
temprate=rate*fm; //frequency modulation if required
index+=temprate;
if(index>WAVETABLE_SIZEF)
{
//new cycle, respecify wavetable for sliding
wavetableindex=reverse[(int)(ratetofloatfactor*temprate)];
wavetable=&sawtoothdata[wavetableindex*WAVETABLE_SIZE];
index-=WAVETABLE_SIZEF;
}
}
}
Notes :
(see linkfile)
Linked files
After help from this group and reading various literature I have finished my
waveform engine. As requested, I am now going to share some of the things
are have learnt from a practical viewpoint.
Problem:
The other condition of this problem is that the waveforms are generatable in
real-time. Additionally, bonuses are given for solutions which allow
pulse-width modulation and triangle asymmetry.
Wavetables:
You can generate a wavetable for the waveform by summing sine waves up to
the Nyquist frequency. For example, the sawtooth waveform can be generated
by:
The wavetable can then be played back, pitched up or down subject that
pitched f*n < Ts/2. Anything lower will not alias but it may lack some
higher harmonics if pitched too low.
You may need to compromise between number of wavetables and accuracy because
of memory considerations (especially if over-sampling). This means some
wavetables will have to cover larger ranges than they should. As long as
the range is too far in the lower direction rather than higher, you will not
alias (you will just miss some higher harmonics).
BLIT intergration
BLIT stands for bandlimited impluse train. I'm not going to go into the
theory, you'll have to read the paper for that. However, simply put, a
pulse train of the form 10000100001000 etc... is not bandlimited.
[note] in the paper they describe m as the largest odd interger not
exceeding the period when in fact their formular for m (which is the one
that works in practive) makes it (int)p or (int)p +1
Now for the clever bit. Lets start with the square/rectangle wave. Through
intergration we get:
C6 = k0/period
Lets look at some rough and ready code. I find the best thing to do is to
generate one period at a time and then reset everything. The numerical
errors over one period are negligable (based on 32bit float) but if you keep
Highly un-optimised code but you get the point. At each sample(t) the
output values are square and triangle respectively.
Sawtooth:
Now, C2 is a bit tricky. Its the average of BLIT for that period (i.e. 1/n
* Sum (k=0,n) (BLIT(k)). [Note] Blit(0) = 0.
I found the best way to deal with this is to have a lookup table which you
have generated and saved to disk as a file which contains a value of C2 for
every period you are interested in. This is because I know of no easy way
to generate C2 in real-time.
Bandlimited waveforms...
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...
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
Base-2 exp
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);
}
Biquad C code
Notes :
Implementation of the RBJ cookbook, in C.
Linked files
http://www.smartelectronix.com/musicdsp/text/filters005.txt
* Enjoy.
*
* This work is hereby placed in the public domain for all purposes, whether
* commercial, free [as in speech] or educational, etc. Use the code and please
* give me credit if you wish.
*
* Tom St Denis -- http://tomstdenis.home.dhs.org
*/
#ifndef M_LN2
#define M_LN2 0.69314718055994530942
#endif
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
/* filter types */
enum {
LPF, /* low pass filter */
HPF, /* High pass filter */
BPF, /* band pass filter */
NOTCH, /* Notch Filter */
PEQ, /* Peaking band EQ filter */
LSH, /* Low shelf filter */
HSH /* High shelf filter */
};
/* compute result */
result = b->a0 * sample + b->a1 * b->x1 + b->a2 * b->x2 -
b->a3 * b->y1 - b->a4 * b->y2;
return result;
}
b = malloc(sizeof(biquad));
if (b == NULL)
return NULL;
/* setup variables */
A = pow(10, dbGain /40);
switch (type) {
case LPF:
b0 = (1 - cs) /2;
b1 = 1 - cs;
b2 = (1 - cs) /2;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case HPF:
b0 = (1 + cs) /2;
b1 = -(1 + cs);
b2 = (1 + cs) /2;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case BPF:
b0 = alpha;
b1 = 0;
b2 = -alpha;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case NOTCH:
b0 = 1;
b1 = -2 * cs;
b2 = 1;
a0 = 1 + alpha;
a1 = -2 * cs;
a2 = 1 - alpha;
break;
case PEQ:
b0 = 1 + (alpha * A);
b1 = -2 * cs;
b2 = 1 - (alpha * A);
a0 = 1 + (alpha /A);
a1 = -2 * cs;
a2 = 1 - (alpha /A);
break;
case LSH:
b0 = A * ((A + 1) - (A - 1) * cs + beta * sn);
b1 = 2 * A * ((A - 1) - (A + 1) * cs);
b2 = A * ((A + 1) - (A - 1) * cs - beta * sn);
a0 = (A + 1) + (A - 1) * cs + beta * sn;
a1 = -2 * ((A - 1) + (A + 1) * cs);
a2 = (A + 1) + (A - 1) * cs - beta * sn;
break;
case HSH:
b0 = A * ((A + 1) + (A - 1) * cs + beta * sn);
b1 = -2 * A * ((A - 1) + (A + 1) * cs);
b2 = A * ((A + 1) + (A - 1) * cs - beta * sn);
a0 = (A + 1) - (A - 1) * cs + beta * sn;
a1 = 2 * ((A - 1) - (A + 1) * cs);
a2 = (A + 1) - (A - 1) * cs - beta * sn;
break;
default:
free(b);
return NULL;
}
return b;
}
/* crc==3062280887, version==4, Sat Jul 7 00:03:23 2001 */
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)));
}
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
Linked files
textcomponent.setText(Integer.toString((int)(getFrequency(parseNoteSymbol(s)) +
0.5F)));
return;
}
catch(IllegalArgumentException _ex)
{
flag1 = true;
}
return;
}
try
{
textcomponent.setText(makeNoteSymbol(getPitch((new
Float(s)).floatValue())));
return;
}
catch(NumberFormatException _ex)
{
flag1 = true;
}
}
public Frequency(float f)
{
frequency = 1.0F;
frequency = f;
}
public Frequency(String s)
throws IllegalArgumentException
{
frequency = 1.0F;
frequency = valueOf(s);
}
Notes :
Flexible-time, non-sample quantized delay , can be used for stuff like waveguide synthesis or
time-based (chorus/flanger) fx.
Code :
class cwaveguide
{
public:
cwaveguide(){clear();}
virtual ~cwaveguide(){};
void clear()
{
counter=0;
for(int s=0;s<MAX_WG_DELAY;s++)
buffer[s]=0;
}
// 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);
// return output
return output;
}
float buffer[MAX_WG_DELAY];
int counter;
};
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);
}
Linked files
-1.801640704047490915268262912739550962585651469641090625323864e-03,
3.537137999745202484462958363064254310959060059520040012524275e-04};
1.273693403357932600826772332014009770786177480422245995563097e-01,
9.305736460357235116035228983545273226942917998946925868063974e-02,
-7.139414716639708714533609307605064767292611983702150917523756e-02,
-2.945753682187581285828323760141839199388200516064948779769654e-02,
3.321267405934100173976365318215912897978337413267096043323351e-02,
3.606553566956169655423291417133403299517350518618994762730612e-03,
-1.073317548333057504431811410651364448111548781143923213370333e-02,
1.395351747052901165789318447957707567660542855688552426721117e-03,
1.992405295185056117158742242640643211762555365514105280067936e-03,
-6.858566949597116265613709819265714196625043336786920516211903e-04,
-1.164668551292854509514809710258991891527461854347597362819235e-04,
9.358867032006959133405013034222854399688456215297276443521873e-05,
-1.326420289452124481243667531226683305749240960605829756400674e-05};
-9.643212009650708202650320534322484127430880143045220514346402e-02,
1.084913025582218438089010237748152188661630567603334659322512e-02,
4.154627749508444073927094681906574864513532221388374861287078e-02,
-1.221864906974828071998798266471567712982466093116558175344811e-02,
-1.284082519830068329466034471894728496206109832314097633275225e-02,
6.711499008795509177767027068215672450648112185856456740379455e-03,
2.248607240995237599950865211267234018343199786146177099262010e-03,
-2.179503618627760471598903379584171187840075291860571264980942e-03,
6.545128212509595566500430399327110729111770568897356630714552e-06,
3.886530628209314435897288837795981791917488573420177523436096e-04,
-8.850410920820432420821645961553726598738322151471932808015443e-05,
-2.424154575703078402978915320531719580423778362664282239377532e-05,
1.277695221937976658714046362616620887375960941439428756055353e-05,
-1.529071758068510902712239164522901223197615439660340672602696e-06};
2.186706877589065214917475918217517051765774321270432059030273e-01,
-2.716885522787480414142192476181171094604882465683330814311896e-01,
-2.180335299932760447555558812702311911975240669470604752747127e-01,
1.383952138648065910739939690021573713989900463229686119059119e-01,
1.399890165844607012492943162271163440328221555614326181333683e-01,
-8.674841156816968904560822066727795382979149539517503657492964e-02,
-7.154895550404613073584145115173807990958069673129538099990913e-02,
5.523712625921604411618834060533403397913833632511672157671107e-02,
2.698140830791291697399031403215193343375766595807274233284349e-02,
-3.018535154039063518714822623489137573781575406658652624883756e-02,
-5.615049530356959133218371367691498637457297203925810387698680e-03,
1.278949326633340896157330705784079299374903861572058313481534e-02,
-7.462189892683849371817160739181780971958187988813302900435487e-04,
-3.849638868022187445786349316095551774096818508285700493058915e-03,
1.061691085606761843032566749388411173033941582147830863893939e-03,
7.080211542355278586442977697617128983471863464181595371670094e-04,
-3.868319473129544821076663398057314427328902107842165379901468e-04,
-4.177724577037259735267979539839258928389726590132730131054323e-05,
6.875504252697509603873437021628031601890370687651875279882727e-05,
-1.033720918457077394661407342594814586269272509490744850691443e-05,
-4.389704901781394115254042561367169829323085360800825718151049e-06,
1.724994675367812769885712692741798523587894709867356576910717e-06,
-1.787139968311359076334192938470839343882990309976959446994022e-07};
2.579269915531893680925862417616855912944042368767340709160119e-05,
-2.813329626604781364755324777078478665791443876293788904267255e-05,
3.362987181737579803124845210420177472134846655864078187186304e-06,
1.811270407940577083768510912285841160577085925337507850590290e-06,
-6.316882325881664421201597299517657654166137915121195510416641e-07,
6.133359913305752029056299460289788601989190450885396512173845e-08};
-3.283207483639617360909665340725061767581597698151558024679130e-01,
-1.265997522158827028744679110933825505053966260104086162103728e-01,
1.973105895650109927854047044781930142551422414135646917122284e-01,
1.011354891774702721509699856433434802196622545499664876109437e-01,
-1.268156917782863110948571128662331680384792185915017065732137e-01,
-5.709141963167692728911239478651382324161160869845347053990144e-02,
8.110598665416088507965885748555429201024364190954499194020678e-02,
2.231233617810379595339136059534813756232242114093689244020869e-02,
-4.692243838926973733300897059211400507138768125498030602878439e-02,
-3.270955535819293781655360222177494452069525958061609392809275e-03,
2.273367658394627031845616244788448969906713741338339498024864e-02,
-3.042989981354637068592482637907206078633395457225096588287881e-03,
-8.602921520322854831713706413243659917926736284271730611920986e-03,
2.967996691526094872806485060008038269959463846548378995044195e-03,
2.301205242153545624302059869038423604241976680189447476064764e-03,
-1.436845304802976126222890402980384903503674530729935809561434e-03,
-3.281325194098379713954444017520115075812402442728749700195651e-04,
4.394654277686436778385677527317841632289249319738892179465910e-04,
-2.561010956654845882729891210949920221664082061531909655178413e-05,
-8.204803202453391839095482576282189866136273049636764338689593e-05,
2.318681379874595084482068205706277572106695174091895338530734e-05,
6.990600985076751273204549700855378627762758585902057964027481e-06,
-4.505942477222988194102268206378312129713572600716499944918416e-06,
3.016549609994557415605207594879939763476168705217646897702706e-07,
2.957700933316856754979905258816151367870345628924317307354639e-07,
-8.423948446002680178787071296922877068410310942222799622593133e-08,
7.267492968561608110879767441409035034158581719789791088892046e-09};
1.186300338581174657301741592161819084544899417452317405185615e-04,
4.943343605466738130665529516802974834299638313366477765295203e-03,
-1.118732666992497072800658855238650182318060482584970145512687e-03,
-1.340596298336106629517567228251583609823044524685986640323942e-03,
6.284656829651457125619449885420838217551022796301582874349652e-04,
2.135815619103406884039052814341926025873200325996466522543440e-04,
-1.986485523117479485798245416362489554927797880264017876139605e-04,
-1.535917123534724675069770335876717193700472427021513236587288e-07,
3.741237880740038181092208138035393952304292615793985030731363e-05,
-8.520602537446695203919254911655523022437596956226376512305917e-06,
-3.332634478885821888782452033341036827311505907796498439829337e-06,
1.768712983627615455876328730755375176412501359114058815453100e-06,
-7.691632689885176146000152878539598405817397588156525116769908e-08,
-1.176098767028231698450982356561292561347579777695396953528141e-07,
3.068835863045174800935478294933975372450179787894574492930570e-08,
-2.507934454948598267195173183147126731806317144868275819941403e-09};
3.010964316296526339695334454725943632645798938162427168851382e-06,
1.531931476691193063931832381086636031203123032723477463624141e-06,
-6.862755657769142701883554613486732854452740752771392411758418e-07,
1.447088298797844542078219863291615420551673574071367834316167e-08,
4.636937775782604223430857728210948898871748291085962296649320e-08,
-1.116402067035825816390504769142472586464975799284473682246076e-08,
8.666848838997619350323013540782124627289742190273059319122840e-10};
1.483675408901114285014404448710249837385836373168215616427030e-01,
3.677286834460374788614690818452372827430535649696462720334897e-01,
5.784327310095244271421181831735444106385099957908657145590104e-01,
5.079010906221639018391523325390716836568713192498711562711282e-01,
7.372450118363015165570139016530653113725172412104955350368114e-02,
-3.127265804282961918033226222621788537078452535993545440716988e-01,
-2.005684061048870939324361244042200174132905844868237447130382e-01,
1.640931881067664818606223226286885712554385317412228836705888e-01,
1.799731879928913037252154295313083168387840791424988422757762e-01,
-9.711079840911470969274209179691733251456735137994201552926799e-02,
-1.317681376866834107513648518146838345477875022352088357523838e-01,
6.807631439273221556739202147004580559367442550641388181886023e-02,
8.455737636682607503362813659356786494357635805197410905877078e-02,
-5.136425429744413245727949984018884707909441768477091944584584e-02,
-4.653081182750671347875833607846979997825771277976548080904423e-02,
3.697084662069802057615318892988581825637896696876361343354380e-02,
2.058670762756536044060249710676656807281671451609632981487139e-02,
-2.348000134449318868560142854519364987363882333754753819791381e-02,
-6.213782849364658499069336123807608293122900450508440420104462e-03,
1.256472521834337406887017835495604463815382993214296088172221e-02,
3.001373985076435951229129255588255746904937042979316054485183e-04,
-5.455691986156717076595353163071679107868762395367234726592273e-03,
1.044260739186025323350755659184734060807432172611689413745029e-03,
1.827010495657279080112597436850157110235336772062961041154607e-03,
-7.706909881231196232880372722955519781655769913634565757339739e-04,
-4.237873998391800799531947768003976978197438302533528661825758e-04,
3.286094142136787341983758471405935405823323072829619248523697e-04,
4.345899904532003379046992625575076092823809665933575578710696e-05,
-9.405223634815760421845190098352673647881298980040512091599943e-05,
1.137434966212593172736144274866639210339820203135670505287250e-05,
1.737375695756189356163565074505405906859746605867772002320509e-05,
-6.166729316467578372152251668422979152169587307212708981768966e-06,
-1.565179131995160159307426993578204733378112742579926503832095e-06,
1.295182057318877573889711232345068147800395721925682566394936e-06,
-8.779879873361286276888117046153049053917243760475816789226764e-08,
-1.283336228751754417819693932114064887075096030264748079976736e-07,
3.761228749337362366156711648187743399164239397803629022612862e-08,
1.680171404922988885554331183691280245962290247654438114807112e-09,
-2.729623146632976083449327361739104754443221903317745768938846e-09,
5.335938821667489905169783227036804533253011117886586305435615e-10,
-3.602113484339554703794807810939301847299106970237814334104274e-11};
5.510185172419193913452724227212507720514144116478727269717859e-01,
1.813926253638400136259098302138614937264260737638175539416540e-01,
-2.613921480306441118856795735210118413900307577511142987337375e-01,
-2.714020986078430556604069575184718123763697177381058877113471e-01,
9.212540708241805260646030910734894258577648089100630012130261e-02,
2.235736582420402317149513960822561717689875252792817094811874e-01,
-3.303744709428937875006612792463031409461636228731285046551636e-02,
-1.640113215318759250156057837165276039181451149292112929401186e-01,
2.028307457564929974897286607551313323418860610791382310375731e-02,
1.122970436181072886950734465075645977754665593869789965874572e-01,
-2.112621235622724100704783293549467048999443844657058425212982e-02,
-7.020739157490110946204219011957565343899895499962369353294028e-02,
2.176585683449997560776882472168730165799461445156766923497545e-02,
3.849533252256919901057154320407596073180564628069920893870768e-02,
-1.852351365015615979794689960740674782817814176166333519597796e-02,
-1.753710100303584537915846117408613551147985251726558719415169e-02,
1.275194393152828646243157404474947115052750581861997731041018e-02,
6.031840650024162816289878206037841640814102314209075233751820e-03,
-7.075319273706152814194039481466556204493276773483821748740018e-03,
-1.134865473356251691289337120013286756337393784110786907825400e-03,
3.122876449818144997419144765125750522437659393621577492535411e-03,
-2.465014005163512031940473100375377210862560761576109755841161e-04,
-1.061231228886651321139357625683805642193648671030425010215075e-03,
3.194204927099011503676530359692366990929679170022583007683112e-04,
2.567624520078737205563856675376636092314813400664190770435450e-04,
-1.500218503490340967673163290447832236259277810659068637402668e-04,
-3.378894834120903434270962452674534330903724108906662510305045e-05,
4.426071203109246077621875303440935335701832843654692827539837e-05,
-2.635207889249186237209225933170897825432335273771458456888097e-06,
-8.347875567854625544366043748844183086765894974439245409223337e-06,
2.397569546840240057403739507525641239509517148980849889986407e-06,
8.147574834779447778085443041422881439860288287528356019216814e-07,
-5.339005405209421154584783682848780965053642859373536945701365e-07,
1.853091785633965019353699857864654181728710556702529908304185e-08,
5.417549179539278736503176166323685597634496102979977037271945e-08,
-1.399935495437998845130909687361847103274208993447892120341999e-08,
-9.472885901812050535221582074673490573092096712822067564903012e-10,
1.050446453696543404071105111096438573423068913105255997908040e-09,
-1.932405111313417542192651899622541612314066389643607507706887e-10,
1.250203302351040941433216718217504240541423430995137507404787e-11};
2.809855532337118833442626085115402941842959475929278883281409e-01,
-1.872714068851562376981887159775791469060265778441667840307934e-01,
-3.179430789993627375453948489797707550898087789160025182664299e-01,
4.776613684344728187950198323031360866349104994035553200788631e-03,
2.392373887803108551973268291945824822214858134512317715815616e-01,
4.252872964148383258147364472170645232684343235486951540533893e-02,
-1.711753513703468896897638515080572393949165942335556397917666e-01,
-3.877717357792001620177594726199572688446488033750771020190283e-02,
1.210163034692242362312637311149062286659377039046006801523826e-01,
2.098011370914481534980883827326017063121637262728447783605518e-02,
-8.216165420800166702291466006164189460916816748629968198028898e-02,
-4.578436241819221637997516339765068825260159169893967894877272e-03,
5.130162003998087915555334881398688958843078494595140394873884e-02,
-4.944709428125628299815920032649550811877887219282751174798211e-03,
-2.821310709490189098113895361900699228886900995412759197674058e-02,
7.661721881646585897329899904308764405384658404613669817843430e-03,
1.304997087108573583052494067883717533043101857128653233783396e-02,
-6.291435370018187780721843581169343900864298634085743861509767e-03,
-4.746568786323113800477796959513558401732252800905982385017245e-03,
3.736046178282523345179052160810332868725126356493155728625572e-03,
1.153764936839481504858282495202271984454410046682805375157566e-03,
-1.696456818974824394274534636412116243080312601322325642741589e-03,
-4.416184856141520063365958900079406737636243682138363561877750e-05,
5.861270593183109933716735450272894035425792347806515678695765e-04,
-1.181233237969554740613021227756568966806892308457221016257961e-04,
-1.460079817762616838924301818082729036314539476811023255670666e-04,
6.559388639305634085303738560455061974369354538271316071502698e-05,
2.183241460466558363365044032984257709791187640963509380549307e-05,
-2.022888292612697682860859987200455702614855595412267510558659e-05,
1.341157750809114719319937553186023660581084151828593222893663e-08,
3.901100338597702610409014129024223853127911530009766793352492e-06,
-8.980253143938407724149926669980791166378388013293887718404796e-07,
-4.032507756879971624098983247358983425236092110387724315244646e-07,
2.166339653278574639176393978510246335478946697396400359281412e-07,
-5.057645419792500308492508924343248979317507866520688417567606e-10,
-2.255740388176086107368821674947175804005323153443170526520277e-08,
5.157776789671999638950774266313208715015419699643333784626363e-09,
4.748375824256231118094453549799175824526559994333227456737433e-10,
-4.024658644584379774251499574468195118601698713554294941756559e-10,
6.991801157638230974132696433509625934021677793453732225542951e-11,
-4.342782503803710247259037552886749457951053124203814185811297e-12};
5.816368967460577833534892038757085635755639698734580573323031e-01,
3.678850748029466984371319740855532278670733841012809062966976e-01,
-9.717464096463814276130048169040892607068486428294030952842447e-02,
-3.364730796417461309562110148848845218930261030262170601615289e-01,
-8.758761458765466140226687673880006154266689569025041229545538e-02,
2.245378197451017129525176510409543155930843160711989062118482e-01,
1.181552867199598604563067876819931882639429216001523151773895e-01,
-1.505602137505796309518094206831433270850173484773520730186277e-01,
-9.850861528996022153725952822686729410420350758543226219234795e-02,
1.066338050184779528831274540522414711301747903916268438037723e-01,
6.675216449401860666895983072443984697329752470942906490126865e-02,
-7.708411105657419356208567671699032054872853174701595359329826e-02,
-3.717396286112250887598137324046870459877639250821705817221557e-02,
5.361790939877949960629041419546536897037332284703545849594129e-02,
1.554260592910229163981295854603203625062268043511894295387375e-02,
-3.404232046065334099320628584033729153497903561399447916116575e-02,
-3.079836794847036661636693963570288706232460663070983852354326e-03,
1.892280447662762841086581178691039363674755753459524525886478e-02,
-1.989425782202736494289461896386235348901617760816745484282494e-03,
-8.860702618046368399013064252456556969199612331833605310278698e-03,
2.726936258738495739871469244610042793734119359765762028996059e-03,
3.322707773973191780118197357194829286271392998979276105842863e-03,
-1.842484290203331280837780430014195744813667655929909114672154e-03,
-8.999774237462950491085382524008429604309720852269895692000702e-04,
8.772581936748274843488806190175921376284150686011179612908221e-04,
1.153212440466300456460181455345639872216326644527860903202733e-04,
-3.098800990984697989530544245356271119416614147098459162436317e-04,
3.543714523276059005284289830559259809540337561365927850248007e-05,
7.904640003965528255137496303166001735463107762646364003487560e-05,
-2.733048119960041746353244004225286857636045649642652816856524e-05,
-1.277195293199783804144903848434605690990373526086311486716394e-05,
8.990661393062588905369930197413951232059323587543226269327396e-06,
5.232827708153076417963912065899772684403904504491727061662335e-07,
-1.779201332653634562565948556039009149458987774189389221295909e-06,
3.212037518862519094895005816661093988294166712919881121802831e-07,
1.922806790142371601278104244711267420759978799176017569693322e-07,
-8.656941732278507163388031517930974947984281611717187862530250e-08,
-2.611598556111770864259843089151782206922842627174274274741722e-09,
9.279224480081372372250073354726511359667401736947170444723772e-09,
-1.880415755062155537197782595740975189878162661203102565611681e-09,
-2.228474910228168899314793352064795957306403503495743572518755e-10,
1.535901570162657197021927739530721955859277615795931442682785e-10,
-2.527625163465644811048864286169758128142169484216932624854015e-11,
1.509692082823910867903367712096001664979004526477422347957324e-12};
6.227474402514960484193581705107415937690538641013309745983962e-02,
1.950394387167700994245891508369324694703820522489789125908612e-01,
4.132929622783563686116108686666547082846741228042232731476147e-01,
5.736690430342222603195557147853022060758392664086633396520345e-01,
4.391583117891662321931477565794105633815363384084590559889493e-01,
1.774076780986685727823533562031556893226571319881417676492595e-03,
-3.263845936917800216385340830055349953447745005769416287177497e-01,
-1.748399612893925042664835683606584215248582345438816346170042e-01,
1.812918323111226960705459766025430918716233584167982942044424e-01,
1.827554095896723746537533832033286839689931924709760567945595e-01,
-1.043239002859270439148009137202747658420968144330108510179290e-01,
-1.479771932752544935782314546369458188243947772922980064071205e-01,
6.982318611329236513756591683950208955110603212379412334701145e-02,
1.064824052498086303236593797715344405836015002929319291715777e-01,
-5.344856168148319149493577269390074213960237013099439431132086e-02,
-6.865475960403591525454725258715351280947435823354011140858001e-02,
4.223218579637203541206570902753288247790857760067894456114927e-02,
3.853571597111186425832144567362328142994885395255438867968781e-02,
-3.137811036306775484244644776337594435094096964336402798072360e-02,
-1.776090356835818354094298625884058170354129044259951019182732e-02,
2.073492017996382475887790073068984224515077665517103399898854e-02,
5.829580555318887971939315747596613038479561943085291072787359e-03,
-1.178549790619302893728624468402138072504226527540325463847390e-02,
-5.287383992626814439198630765217969804966319971038003993984480e-04,
5.601947239423804853206514239940474788977188460452053462770324e-03,
-9.390582504738289646165698675070641765810790863514339205205998e-04,
-2.145530281567620980305401403432221668847980295600748913748902e-03,
8.383488056543616046381924054554052104937784379435436426690560e-04,
6.161382204574344193703789012696411561214682388271673214197731e-04,
-4.319557074261807466712901913481943478521991611607433971794602e-04,
-1.060574748283803889966150803551837402553866816191659959347053e-04,
1.574795238607493590547765666590811258087715699737771458390360e-04,
-5.277795493037868976293566636015627609248847457646525246271036e-06,
-4.109673996391477816326502438997466532822639385119090230965252e-05,
1.074221540872195031273584409245060623104931330938273936484593e-05,
7.000078682964986734859102495210684809643657474253921074934684e-06,
-3.887400161856795187587790410706550576033603097954065074023128e-06,
-4.650463220640262639231145944536092973446596027469833860001618e-07,
7.939210633709952088373459255067360793370284788682979065122810e-07,
-1.079004237578671411922961583845716126060658213943840375162654e-07,
-8.904466370168590769052983362721567202750591914741016835071257e-08,
3.407795621290730008673832107214820587991557116806912418558069e-08,
2.169328259850323106986222296525930099935873861026310788086221e-09,
-3.776010478532324328184043667556576385639846460337894963138621e-09,
6.780047245828636668305808192607091517605349478677442468580825e-10,
1.002303191046526913509281844136258004034177309673269533418644e-10,
-5.840408185341171468465492447799819262905317576847426970757700e-11,
9.130510016371796243923232926650252570239054815939483900056681e-12,
-5.251871224244435037810503452564279828539007071678724285717464e-13};
-1.309465606856955151282041809232358209226373823424148862843577e-09,
1.521614984778521740775073159445241799352681846880808663329946e-09,
-2.415526928011130660506395791946234018673860470542996426005750e-10,
-4.374986224293654395069947682013996351823060759948583134078918e-11,
2.213662088067662485181472969374945928903854605356443772873438e-11,
-3.295790122476585807069953975043096139541415768606924980926275e-12,
1.828188352882424933624530026056448539377272017834175009418822e-13};
-1.004326041333422601781848560432120920634648692782357855473103e-05,
1.247900317574834146052381692752796047052443265982232422642017e-06,
1.840363734517769191684379309039277810350620305330900536404818e-06,
-6.670215479954892588747450458085225880096882699397256774967304e-07,
-1.757461173209842779903676264971918635870906983281392939812547e-07,
1.490660013535362170989340065033061951960933954388633507264360e-07,
-8.262387315626556965966429243600984899650039704831080988658278e-09,
-1.784138690875710077191713941441263246560738410213624546116655e-08,
5.044047056383436444631252840057862002264087720676808580373667e-09,
6.944540328946226952976704718677697525410051405055662575530111e-10,
-6.077041247229010224760245305596307803830053533836849384680534e-10,
8.492220011056382105461206077240377024404404638947591299761197e-11,
1.867367263783390418963879146175452376940453585791428841004699e-11,
-8.365490471258800799349289794397908900767054085216008197372193e-12,
1.188850533405901520842321749021089497203940688882364518455403e-12,
-6.367772354714857335632692092267254266368934590973693820942617e-14};
1.877120925723650133179338154344873477230567340668548016358682e-03,
1.087053942226062966738944397844498417945523630053411148182206e-03,
-1.000778327085680541055696707760062870925897014530348262794137e-03,
-2.000711363076779808296301110796026470163110202848894744316755e-04,
4.111283454742767033424740543004041500054889660665367490129376e-04,
-2.292018041214499897382298271438084577065170236103859181134525e-05,
-1.293044840080720609161466939678226852440475312744714379499074e-04,
3.645026068562774967665464216602750761690984830805534178557146e-05,
2.913344750169041218495787251929571015775436967652945386217480e-05,
-1.657328395306616289863396387854880512976861409870690029695161e-05,
-3.593644804025187638066915189731950450034629392522542962477168e-06,
4.750609246452552850197117564759363194953518317428400241629683e-06,
-3.029054592052818286474228294307141792053791695855058563299597e-07,
-8.975701750636280734511651941681818767895052287332471537510510e-07,
2.633898386997696553900967704111473475368019612368922599394214e-07,
9.387197411095863026484410601284876812292554863800653292318725e-08,
-6.286156922010786166768503252870590953166867739448102804392389e-08,
1.076591906619196137385201975028785139607670319821266803566785e-09,
7.768978854770062238895964639391324551611701293594055935346266e-09,
-1.893995386171984147774611076618946011337498790609031626697228e-09,
-3.426800863263089001811012278889864200550342566386405676893537e-10,
2.407099453509342962399811991929330725186626582891090462239366e-10,
-2.940589250764532582888473974638273664244682541297835986306504e-11,
-7.832509733627817032356556582819494794884131433810848844709881e-12,
3.152762413370310423797539876893861621418382024668704492620948e-12,
-4.285654870068344101898185073376307686875386259541180967347399e-13,
2.219191311588302960934661700068023727737812918006011019184982e-14};
4.380166467141773250305407710250135373016604593736480428415303e-02,
5.671236574473569492590636983030617493807140224924978946302257e-02,
-3.567339749675960965780819743176056734137251336781389369397564e-02,
-3.226375891935220815954913483392725682165778426411705216010280e-02,
2.707861959529418272206848318420006522973840949600186710327776e-02,
1.528796076985739546052896626042375110302102640936712142026221e-02,
-1.839974386811734118728169880549148389603890445324127330811811e-02,
-5.296859666131086629169938675330494864053932988161015674773617e-03,
1.091563165830488927536881480211929049886878831313700460017968e-02,
6.196717564977244383592534999284255315694546230739551683085460e-04,
-5.530730148192003288871383856487027893918513053091795443517653e-03,
8.433845866620933982126003584365932145598126087481400294999080e-04,
2.324520094060099304385756339638431339131122661576649123053845e-03,
-8.609276968110423879660725173525347077801305237644122054954659e-04,
-7.678782504380918697963922441514742758516706160788123977340073e-04,
5.050948239033467796256544554086554367969638627715114003635557e-04,
1.724825842351709725545759714374272164367933578194910678479473e-04,
-2.161718301169633804271038862087964094429005266172702380483361e-04,
-8.548305467584070994787824796256108217987765582429940610377190e-06,
6.982008370808327851082027193100914402221658444151889697045071e-05,
-1.339716863293971629296314599448901465078920406443516550195793e-05,
-1.636152478725426488654528710478856195004608401773950511915162e-05,
7.252145535890469015723401169934327900622894130695550273452916e-06,
2.327549098493686509557358103785598216688723737824121617676858e-06,
-2.187267676996166416699555236143059249832615777542412142603694e-06,
1.099474338526203304286307383463498542376432972308342428764576e-08,
4.261662326011572446469849114416378817419458434583398455985144e-07,
-1.000414682354500898864979332965559934104686157639553850670490e-07,
-4.764379965139453357729154748688006975561934425368712852985388e-08,
2.605442754977625431940885841950955928085338672381046225838880e-08,
5.553397861397053982967618072672572206490972606026556946910028e-10,
-3.331105680467578245901976412732595596538702049437802824373020e-09,
6.984862691832182584221096665570313611280449991512869846064780e-10,
1.613622978270904360610418704685783656905979134344922647926295e-10,
-9.461387997276802120884525814092001871993910062127702293573920e-11,
1.000105131393171192746337860330428369495110180346654025287492e-11,
3.239428638532286114355931428908079297696045600279108835760520e-12,
-1.185237592101582328254231496310584611948560976394420324137742e-12,
1.543997570847620046003616417646988780670333040868954794039905e-13,
-7.737942630954405708679963277418806436871098329050829841696327e-15};
5.511398409142754983590484577074663132074992263886810324421617e-01,
4.294688082061372955430413148799008354573408538414331312236645e-01,
2.716921249736946422305354732634261873401679092095992827198308e-02,
-3.109551183195075186926560285811004715398678229333522634202008e-01,
-2.179784855235633521693544507220105631639547435903112747133934e-01,
1.401782887652732681656253206993073895422881511380152633441096e-01,
2.249667114737370933697297905066886078307490136415302624018330e-01,
-4.992634916046823977000579399730138693074543903234092797936484e-02,
-1.869623608957154494374577196258383009208655076187653847079167e-01,
1.543698842948893409652995335281236231845293548571166883219023e-02,
1.450895009319931981518942907854879059128872873116921504156674e-01,
-8.139832273469236863527708715566588550006680549152344840146851e-03,
-1.076127733234956326668605511648013952380301953590447106075614e-01,
1.094129745236496925725237900637802669504835743555466811796369e-02,
7.535361174328140695528289751109133941376701984419452638686226e-02,
-1.488002661810482202699555987503429289100801979910046913257306e-02,
-4.861907546485433003537603385831190109391263542044516048871113e-02,
1.615417156598591113619453864586701665635869166193865651960591e-02,
2.804761936675616906861927211659154977049392281479113764697785e-02,
-1.427627527776351943309800140756746087215016194775579070599004e-02,
-1.390055293926652880755898888934447671732373519028670201124816e-02,
1.051763948737184089128633441244991643331033825102031908858652e-02,
5.516163573310992566561289762241160214476622662764637181816550e-03,
-6.520852375874612553325469682628530079210293774541131381751695e-03,
-1.428264223218909891400516038687842292177211292295049238921068e-03,
3.393066776715931928419358796960612411097347419792355896915546e-03,
-6.397901106014600492881202314307290077992972755016494062875201e-05,
-1.459041741985160943114515221598080223845239255190055621901681e-03,
3.431398296904734438118401084929505912208229684629857530009147e-04,
4.998816175637222614896912406679513231966722440032799024979502e-04,
-2.396583469402949615285646688069476140260781708006174912535660e-04,
-1.243411617250228669409179807383399199879641177993453588807726e-04,
1.089584350416766882738651833752634206358441308880869184416670e-04,
1.501335727444532997071651937630983442758297688087711521441229e-05,
-3.631255157860086164261313773172162991107348698083164489165837e-05,
4.034520235184278839752741499546098778993926344831736074409765e-06,
8.795301342692987765440618030678349427367022581211855857458220e-06,
-3.035142365891509630069007852947057220760887215249503512783023e-06,
-1.369060230942940782050489751987123955074404782177163471279285e-06,
9.810015422044371573950976088058064384946146188110905321673802e-07,
5.327250656974915426977440959783080593776012130063170688309127e-08,
-1.975925129170206248152121156696590501303803187231928513867046e-07,
3.616826517331004805247567218405798591329788122337274956172315e-08,
2.328309713821409644308538888589329921141948539678106680777082e-08,
-1.061529602150252306500404266150823962402673780484965538270541e-08,
-6.474311687959861398702581539341954438747926255671605657095807e-10,
1.408568151025177427076547804944585301332087108125727813194374e-09,
-2.524043954153353306183643702933218308617979467184848456565837e-10,
-7.348930032486263904766913919653624379586487437915175106407348e-11,
3.692108808871129411604189196259677640440919369478263728899602e-11,
-3.327008967125979929910636246337150851642079794871116041187279e-12,
-1.324334917243963163878274345609465717294426628053460151843705e-12,
4.445467096291932163298411852093011459626037560439178917611592e-13,
-5.559442050579014337641375730083534521513818164827556763756543e-14,
2.699382879762665647295493928801387173921314576598505507855504e-15};
-5.259809282684322782648914338377962890245975842272425408122506e-05,
-1.294045779405512723950480259110995722517019870286295908085366e-05,
1.824268401980691220603850117995712615809177092802967489081228e-05,
-6.361781532260254953363913076575914206506177493714496098327288e-07,
-4.558309576264423135123964145585288808181431652781253437738445e-06,
1.202889036321620990296134494079846952404216422923750605507047e-06,
7.560047625595947819392627283726711361273296630256477108501994e-07,
-4.285970693151457255418342315045357407199066350632593899896712e-07,
-5.003361868748230293692887222336390314786090450819216035110269e-08,
8.965966311957728376981484572655177545054433542721057470726361e-08,
-1.219924359483373093110396748985081720383992859961285213840740e-08,
-1.104383021722648979552131128575075255513372249283096583736746e-08,
4.250422311980592983740943309197245384991941251563471671065543e-09,
4.384387799940474369553236949848427579687147486892033587998023e-10,
-5.881091462634605628881794361152305108432139465417759716875076e-10,
8.904723796221605490455387579189371137903330749397374037644960e-11,
3.263270741332907875981844980104948375955551273115386408552080e-11,
-1.430918765169202320188022211739750594608742928641485026836608e-11,
1.075610653501062115165734990153347111902874668945095034791947e-12,
5.361482229611801638107331379599434078296259332654994508124989e-13,
-1.663800489433402369889818192962259823988673359967722467427927e-13,
2.000715303810524954375796020597627467104635766752154321244151e-14,
-9.421019139535078421314655362291088223782497046057523323473331e-16};
-2.524858297747649929258392207837724793937727346177294684700378e-02,
3.868706076024496481748675031852528047303323816250150793091832e-02,
1.070326582001954942654534968137727769698168853186071888736311e-02,
-2.572876175473297336123211392278301875687760837710204579628265e-02,
-2.167758617353607324783298657172830203896433848418061622436727e-03,
1.531695411585766548347442266431874060229304787191589430967538e-02,
-1.594288782414604768637856446111392724059836934455189837500244e-03,
-7.953540387057939240459305406538116220678495240302592677582773e-03,
2.389062408165908575935815973439728988151836094753689966108405e-03,
3.480800953405711999411461002429227385937942254778524257436278e-03,
-1.860718214455795912074482150710567824317228203897000129729967e-03,
-1.204309257604658876916644980097327372892008586047095719636829e-03,
1.074380696351291355073899234941719080473877020595209197706651e-03,
2.727305847336937211749282358350196461733595290569540045817329e-04,
-4.908329007590351474487792254066540683724948757382104652497458e-04,
4.393166251766185755059005296958129844094063524324718175254673e-06,
1.780431898251245351831728023200069586928513661382622116969992e-04,
-4.160438516273709306234368807933932360567787692918883118883736e-05,
-4.929564423417301834310231482621574127409950921583062559483686e-05,
2.423335398816890365621188379922041046073808819182024026589770e-05,
9.070805757828453800203677464921508178468256685438211818575040e-06,
-8.866121366757736169176034432364298134186929098274651022820760e-06,
-3.607516102879771631230351118595069330196155459105589342866625e-07,
2.288371276141527305481395545993763010565968667577768164201792e-06,
-4.426923407952870147984002129341809185622768353983550670755106e-07,
-3.985791291985944076942626511739220753169387460984290019185514e-07,
1.822443332571053437467128998002798233969112236553215291639303e-07,
3.377972703730854377516206663481869099376154259897212784144779e-08,
-3.987838198518880722819502850814936369197384392561970319349663e-08,
3.672863576838181340505563759379169099717712645283448779390320e-09,
5.111211857347453839549366593998758891130921028374576213256027e-09,
-1.671392677251932495173219614104411841891545601521784559793012e-09,
-2.496402105246193648073519269370197331176405371538404298745013e-10,
2.426833102305682309891302883361232297664099485514601790344279e-10,
-3.049574453945863430361296931455141500128170151643206937547928e-11,
-1.420236859889936792437077844940412749343225644487770840543290e-11,
5.509414720765524548752673631197714447818740985929081064907524e-12,
-3.343481218953278765982532722689984725170758193566174566492199e-13,
-2.152488386833302618520603545685994753329478275805993737095214e-13,
6.214740247174398315576214699577230693021307854673557214652751e-14,
-7.196510545363322414033654470779070592316600780697558361083151e-15,
3.289373678416306368625564108782095644036415401902518812978798e-16};
1.241524821113768081954449898210969172708199672428635378051285e-01,
2.877650592337145629334256618087718872558560120999651277991839e-01,
4.784787462793710621468610706120519466268010329031345843336104e-01,
5.305550996564631773133260223990794445605699030503652382795600e-01,
2.903663295072749510455945186199530115755664977934564128822650e-01,
-1.282468421744371672912377747048558427612774932943748628650824e-01,
-3.315253015083869417715548463087537345035828886426345397256876e-01,
-1.038919155156404718287260506925867970596448618647006698388596e-01,
2.169072201874275950610018667099322465619408030256534197819784e-01,
1.666017504122074437311574334509261366682993700573488534577890e-01,
-1.273373582238011562843862636988693890108793629966541695807247e-01,
-1.609249271778668063014799490429649196614628857267382976958607e-01,
7.799184693794810738265349531832015087096882277333968473726399e-02,
1.341259602711361284802399913977387999358280900708582462625539e-01,
-5.448296806413904636632671383140642554265865948686157271017286e-02,
-1.029475969928140852342073823689090498245496056845473569066667e-01,
4.357609464963129726428486610925800727137724136370669421246609e-02,
7.318523543679560555546221335452045680757998947493883124934567e-02,
-3.701283841786244960356402125554190040750079009127461655784927e-02,
-4.743855964527776247220681410983851377889756018716427358008296e-02,
3.073974657395934459931226513844134346305562928466993208164603e-02,
2.722835075635419610095839895805858855202745897718117731496534e-02,
-2.367173792282636485046786438094940427456079528043555566867110e-02,
-1.314398001665716086105827506126287041342680578404007359439612e-02,
1.640937419986519252112261495537409592363156309874473310057471e-02,
4.713649260999809905918876125437488856235874027077755004539205e-03,
-1.004550670836151917439146861146431000364858401181337134891421e-02,
-6.194748845153872839014356621835501857322345445234809347431098e-04,
5.334950768759936032170270195983921511565539100791906952901398e-03,
-7.692127975067836975989490900561029844887285335804349474993607e-04,
-2.399453943537055863933124827688081952701780599883067560501870e-03,
8.589959874363661955444898475746536583497522107459291718900058e-04,
8.751999064078688732610570055224339733760304773327228476255647e-04,
-5.527355762144197975516415296735124460550632283763688359649888e-04,
-2.326732140233531635428863212833942245597361085708567528230733e-04,
2.650772397558057819755811309071002543822145660933016957735937e-04,
2.660050018453441903046828468025589086403126180798464347801678e-05,
-9.914697770780134603580350758869378471802751837608461971022567e-05,
1.353117227249649581251887376414486225127346352042209141315562e-05,
2.844951419697807376503080001943765930601242225183893658540032e-05,
-1.057657494257950623848316304755218120233253479317574337409622e-05,
-5.710826510998303938275050074333400305512451419983646591762318e-06,
4.169871758547028398316761659984928804362023643629741358799744e-06,
4.979718101421307748081857636471761057429219265531618602960147e-07,
-1.116306534817008428597995070751765080383261658112656948526954e-06,
1.448195708333185127061180618150009526758658641231104901703561e-07,
2.025990666667859216690536885693725545344933235432307649205497e-07,
-7.526701740412589411177481797841044281662555785969415398369019e-08,
-1.990346501531736915866180448337614967570744211158241514589121e-08,
1.740423332936068076497051274445147160190783847854409836489662e-08,
-8.665744261368722215864741166245385888818567571145958531936939e-10,
-2.316501946995482751582294240136010067415084499025753117941001e-09,
6.446378210323402313101214894500231181606520211579581132442548e-10,
1.300410318609415248880403259300467720631189120978928377152233e-10,
-9.904774537632409015479530333979124540183199174591377762845227e-11,
1.004208735461769864836516428998306778031143650101842361622330e-11,
6.080125354000167254059025929915591291115751734288584563131636e-12,
-2.107879108915301546285370395443778864676275235126044599683271e-12,
9.799451158211597727901178520526388692140586041163624252991805e-14,
8.579194051799733179793112298652600511486581216528683482143106e-14,
-2.317083703906408481078257081903089523234020423092175261925515e-14,
2.587338381935699555813538163144986688834142571207152879144731e-15,
-1.148944754480590128244815794312606245147888158018823490936280e-16};
6.137754586740521089596801883631921221145712545042519987641234e-03,
1.428088794070762107355585870669842132609159040625895090070111e-03,
-3.357644380922383229567732565298665639037348585961127075507937e-03,
7.615969435172736546769649923895317451534703066016116257300160e-06,
1.549637469702362975561719246539787717204438637997824935787688e-03,
-3.346692164250854961608526121524596908041109918361306282201310e-04,
-5.864810318991817532175809224131456738367101035694188223408841e-04,
2.648328819961289039302810122699710966048565368047575218693134e-04,
1.700012283661249043584690194716767771204207742625746308522935e-04,
-1.365883072261161602559926714744746422567509177443594045709653e-04,
-2.976995962848509743944225866488519668585242655980656646544319e-05,
5.304143122913310222538317980686374696005605533475685587486683e-05,
-2.437001526827789860990429478540556752694389693432668831073769e-06,
-1.572442077270281693663288966405861215692805972737981986121447e-05,
4.308047861716731191350493437937513220737450410132878032163179e-06,
3.353345862871309889390877168046133657377105681618708355266688e-06,
-1.895929617693153288493891051875444439753318548105998166574535e-06,
-3.903931733287306166657519468494511920760767388397825775326745e-07,
5.302368616904760917074352633915743250769600635829229600812520e-07,
-3.700308378205124537986402644918879149894035910106489082512364e-08,
-9.990396944534900755781728477561240762191443422318249128866740e-08,
3.008188650719066928230268918661718274504955045022550217051301e-08,
1.084902733789934825266560240100449884702749303326571747323086e-08,
-7.458116552893037631192407611262788593505988638365840409367117e-09,
5.897951310384361575470355861162022501172491937837712969865619e-11,
1.030823345485433383811700481488557422005210168069163779730908e-09,
-2.433545573751672936168877250405940817227367937230289801251648e-10,
-6.407938256501889018430608323235974406219193176918284664973727e-11,
4.000536627253744510742788201354093006471710416671002244302586e-11,
-3.125639357108557540598098228678150768528121565391376265627294e-12,
-2.567065476155081449204643852428401530283519685638256074752850e-12,
8.015088533687900921948605418789324826115616416343391081288979e-13,
-2.597954328893848084315198205094389145706680129208998638802995e-14,
-3.397720856796267431956783825659069596940335130100871912329556e-14,
8.624037434720089202680337663692777682810714650060805832406135e-15,
-9.298012529324185420921555664719863501848315099116725184370339e-16,
4.014628712333488654318569164614220308046021091178184654250982e-17};
4.397519752934862993862182898358763783110745559238982179690132e-02,
-2.944210395891145711100715969898758940722458887377844633443675e-01,
-2.468070369781255270524798278622698446566520718230313889086016e-01,
9.811420416311477050518401371401568038943437322299913514049728e-02,
2.465372776089742110529709111809595434656418762898152706621356e-01,
7.278515095792229009687682299460382878643139026668958884429641e-03,
-1.993372056086496198603363400094784142714162256792182570541036e-01,
-4.586140074639271639145126228774831743002971373998329604574394e-02,
1.541062366276428841776316300420654875883842819413623395358262e-01,
5.027618007353842862036816972809884096761706036019748316890913e-02,
-1.188037543101356316801816931383547446073152951044444224449501e-01,
-3.988085357551317584091699967924044034100374257075864260934102e-02,
9.115678225801654406336059281306715151058903055370522031843771e-02,
2.503872144956848989919484296709846860569180993040383621980546e-02,
-6.820901663681751124880436344265538690580358108714540763125119e-02,
-1.131910031681742794381808082173695022123056280821611354577883e-02,
4.851308354780908538616267662315735632292989749013261207046367e-02,
1.424972661765391603147802607378542396323429657660009755652404e-03,
-3.198072067763969654470293513742344601172739688274251641873778e-02,
3.984040198717004857397179486790082321314291366656151213429068e-03,
1.906359478062535932877576164368198274858108513696832728889209e-02,
-5.657813245058818380424016973516714570499161434975761798379020e-03,
-9.990263473281372348001743806489172665465685056975652497503772e-03,
5.022989106665829004699819220796538830393945994687289792465541e-03,
4.413484835350575251918616780287775585471012556848037301025999e-03,
-3.484541445404883311209541395428535732697661971818727286003028e-03,
-1.503074066296643749549363655363411879858070202740814054964603e-03,
1.990793771851737270404293245701878186600899439513475823305914e-03,
2.776812795712026068152384207605140383490242756921936501940389e-04,
-9.463403823261101964604918059447913047725482130063492242779878e-04,
8.614565758992702032613879159402330909634737204578606399403107e-05,
3.693507284967510502620040341882236687749563414433432842567511e-04,
-1.155118895843527096848376999413102395191976350936666573818799e-04,
-1.131899468084665671727391922924411467938450743565106978099456e-04,
6.694741196930590257104231749283786251555566773398199990337698e-05,
2.375106683660860777161950832380341362257503761490580896617678e-05,
-2.731390824654337912922346414722045404779935825834384250023192e-05,
-1.183471059985615942783182762352360917304348034947412986608322e-06,
8.372218198160788432628056043217491552198857358432112275253310e-06,
-1.586145782434577495502614631566211839722879492827911790709498e-06,
-1.870811602859180713762972281154953528056257451900381097476968e-06,
8.311421279707778528163597405935375886855029592150424544500718e-07,
2.548423522556577831218519052844387478819866531902854523544709e-07,
-2.455377658434232699135878286794578515387138194247693201846263e-07,
2.753249073339512254085076456700241929492720457889076058451072e-09,
4.799043465450992009934526867650497683545716858606119786327559e-08,
-1.156093688817008406756913949175208452083765368825442482226093e-08,
-5.612784343327791397474114357094368557982413895802980814813369e-09,
3.138841695782424018351567952158415003571380699236147752239001e-09,
1.090815553713751810964713058800448676068475673611349566405716e-10,
-4.512545778563249634425200856088490195004077806062978067796020e-10,
8.962418203859611987065968320295929679774693465791367610044773e-11,
3.037429098112535221800013609576297196061786927734556635696416e-11,
-1.599716689261357143200396922409448515398648489795044468046420e-11,
8.876846287217374213524399682895564055949886050748321818411161e-13,
1.070969357114017002424433471621197579059927261727846375968378e-12,
-3.029285026974877268896134589769473854669758797446795757329862e-13,
5.542263182639804235231685861028995158694397223907295269180336e-15,
1.338071386299105896025578761458472955294763310766371178363783e-14,
-3.204628543401749860439316638848579711789176444320134355253750e-15,
3.339971984818693213132578777712503670014459411167839211495237e-16,
-1.403274175373190617489823209168013922564353495443487431242610e-17};
7.387757452855583640107787619408806919082115520707105052944171e-03,
-2.248053187003824706127276829147166466869908326245810952521710e-03,
-3.394523276408398601988475786247462646314228994098320665709345e-03,
1.816871343801423525477184531347879515909226877688306010517914e-03,
1.263934258117477182626760951047019242187910977671449470318766e-03,
-1.111484865318630197259018233162929628309920117691177260742614e-03,
-3.280788470880198419407186455190899535706232295554613820907245e-04,
5.490532773373631230219769273898345809368332716288071475378651e-04,
1.534439023195503211083338679106161291342621676983096723309776e-05,
-2.208944032455493852493630802748509781675182699536797043565515e-04,
4.336726125945695214852398433524024058216834313839357806404424e-05,
7.055138782065465075838703109997365141906130284669094131032488e-05,
-3.098662927619930052417611453170793938796310141219293329658062e-05,
-1.639162496160583099236044020495877311072716199713679670940295e-05,
1.354327718416781810683349121150634031343717637827354228989989e-05,
1.849945003115590390789683032647334516600314304175482456338006e-06,
-4.309941556597092389020622638271988877959028012481278949268461e-06,
4.854731396996411681769911684430785681028852413859386141424939e-07,
1.002121399297177629772998172241869405763288457224082581829033e-06,
-3.494948603445727645895194867933547164628229076947330682199174e-07,
-1.509885388671583553484927666148474078148724554849968758642331e-07,
1.109031232216439389999036327867142640916239658806376290861690e-07,
5.350657515461434290618742656970344024396382191417247602674540e-09,
-2.252193836724805775389816424695618411834716065179297102428180e-08,
4.224485706362419268050011630338101126995607958955688879525896e-09,
2.793974465953982659829387370821677112004867350709951380622807e-09,
-1.297205001469435139867686007585972538983682739297235604327668e-09,
-1.031411129096974965677950646498153071722880698222864687038596e-10,
1.946164894082315021308714557636277980079559327508927751052218e-10,
-3.203398244123241367987902201268363088933939831689591684670080e-11,
-1.398415715537641487959551682557483348661602836709278513081908e-11,
6.334955440973913249611879065201632922100533284261000819747915e-12,
-2.096363194234800541614775742755555713279549381264881030843258e-13,
-4.421612409872105367333572734854401373201808896976552663098518e-13,
1.138052830921439682522395208295427884729893377395129205716662e-13,
-4.518889607463726394454509623712773172513778367070839294449849e-16,
-5.243025691884205832260354503748325334301994904062750850180233e-15,
1.189012387508252879928637969242590755033933791160383262132698e-15,
-1.199280335852879554967035114674445327319437557227036460257649e-16,
4.906615064935203694857690087429901193139905690549533773201453e-18};
3.307757814110146511493637534404611754800768677041577030757306e-01,
4.965911753117180976599171147718708939352414838951726087564419e-01,
4.933560785171007975728485346997317064969513623594359091115804e-01,
2.130505713555785138286743353458562451255624665951160445122307e-01,
-1.828676677083358907975548507946239135218223185041410632924815e-01,
-3.216756378089978628483471725406916361929841940528189059002548e-01,
-6.226650604782432226643360160478765847565862101045597180310490e-02,
2.321259638353531085028708104285994998671615563662858079262996e-01,
1.499851196187170199586403453788927307298226028262603028635758e-01,
-1.417956859730596216710053144522330276392591055375830654519080e-01,
-1.599125651582443618288533214523534937804208844386102639177693e-01,
8.563812155615105741612217814369165313487129645536001850276987e-02,
1.414147340733826800884683119379170594092606174915755283496153e-01,
-5.658645863072738145681787657843320646815509410635114234947902e-02,
-1.147311707107443752394144019458942779715665489230169950201022e-01,
4.309589543304764288137871223616030624246568683595408792078602e-02,
8.720439826203975011910714164154456762073786124233088471855868e-02,
-3.660510340287429567372071039506772372567938710943432838908247e-02,
-6.176620870841315993604736705613246241897497782373337911398117e-02,
3.198987753153780630818381136366859026137035450576631134176875e-02,
4.005498110511594820952087086241114309038577379366732959648548e-02,
-2.689149388089451438550851767715967313417890393287236700072071e-02,
-2.311413402054931680856913553585621248925303865540203357180768e-02,
2.090464525565524340215982365351342094670261491526831672682244e-02,
1.129049727868596484270081487761544232851115891449843967151657e-02,
-1.470188206539868213708986402816605045648481224662435114088245e-02,
-4.131306656031089274123231103326745723188134548520938157995702e-03,
9.214785032197180512031534870181734003522861645903894504302286e-03,
5.625715748403532005741565594881148757066703437214522101740941e-04,
-5.071314509218348093935061417505663002006821323958752649640329e-03,
7.169821821064019257784165364894915621888541496773370435889585e-04,
2.400697781890973183892306914082592143984140550210130139535193e-03,
-8.448626665537775009068937851465856973251363010924003314643612e-04,
-9.424614077227377964015942271780098283910230639908018778588910e-04,
5.810759750532863662020321063678196633409555706981476723988312e-04,
2.817639250380670746018048967535608190123523180612961062603672e-04,
-3.031020460726611993600629020329784682496477106470427787747855e-04,
-4.555682696668420274688683005987764360677217149927938344795290e-05,
1.262043350166170705382346537131817701361522387904917335958705e-04,
-1.155409103833717192628479047983460953381959342642374175822863e-05,
-4.175141648540397797296325065775711309197411926289412468280801e-05,
1.334176149921350382547503457286060922218070031330137601427324e-05,
1.037359184045599795632258335010065103524959844966094870217687e-05,
-6.456730428469619160379910439617575420986972394137121953806236e-06,
-1.550844350118602575853380148525912999401292473185534395740371e-06,
2.149960269939665207789548199790770596890252405076394885606038e-06,
-8.487087586072593071869805266089426629606479876982221840833098e-08,
-5.187733738874144426008474683378542368066310000602823096009187e-07,
1.396377545508355481227961581059961184519872502493462010264633e-07,
8.400351046895965526933587176781279507953080669259318722910523e-08,
-4.884757937459286762082185411608763964041010392101914854918157e-08,
-5.424274800287298511126684174854414928447521710664476410973981e-09,
1.034704539274858480924046490952803937328239537222908159451039e-08,
-1.436329487795135706854539856979275911183628476521636251660849e-09,
-1.349197753983448821850381770889786301246741304307934955997111e-09,
5.261132557357598494535766638772624572100332209198979659077082e-10,
6.732336490189308685740626964182623159759767536724844030164551e-11,
-8.278256522538134727330692938158991115335384611795874767521731e-11,
1.101692934599454551150832622160224231280195362919498540913658e-11,
6.291537317039508581580913620859140835852886308989584198166174e-12,
-2.484789237563642857043361214502760723611468591833262675852242e-12,
2.626496504065252070488282876470525379851429538389481576454618e-14,
1.808661236274530582267084846343959377085922019067808145635263e-13,
-4.249817819571463006966616371554206572863122562744916796556474e-14,
-4.563397162127373109101691643047923747796563449194075621854491e-16,
2.045099676788988907802272564402310095398641092819367167252952e-15,
-4.405307042483461342449027139838301611006835285455050155842865e-16,
4.304596839558790016251867477122791508849697688058169053134463e-17,
-1.716152451088744188732404281737964277713026087224248235541071e-18};
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;
Conversions on a PowerPC
Cubic interpollation
Type : interpollation
References : Posted by Olli Niemitalo
Linked file : other001.gif (this linked file is included below)
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;
Linked files
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
}
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);
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..
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;
}
Denormal numbers
Notes :
this text describes some ways to avoid denormalisation. Denormalisation happens when FPU's go mad processing very
small numbers
Linked files
Denormal numbers
----------------
Here's a small recap on all proposed solutions to prevent the FPU from
denormalizing:
When you feed the FPU really small values (what's the exact 'threshold'
value?) the CPU will go into denormal mode to maintain precision; as such
precision isn't required for audio applications and denormal operations are
MUCH slower than normal operations, we want to avoid them.
All methods have been proposed by people other than me, and things I say
here may be inaccurate or completely wrong :). 'Algorithms' have not been
tested and can be implemented more efficiently no doubt.
If I made some horrible mistakes, or left some stuff out, please let me/the
list know.
To detect if a denormal number has occured, just trace your code and
look up on the STAT FPU register ... if 0x0002 flag is set then a denormal
operation occured (this flag stays fixed until next FINIT)
NOTES:
This is the least computationally efficient method of the lot, but has the
advantage of being inaudible.
Please note that in every feedback loop, you should also check for
denormals (rendering it useless on algorithms with loads of filters,
feedback loops, etc).
CODE:
// inner-loop:
NOTES:
Less overhead than the first solution, but still 2 mem accesses. Because a
number of the values of denormalBuf will be denormals themselves, there
will always be *some* denormal overhead. However, a small percentage
denormals probably isn't a problem.
Use this eq. to calculate the appropriate value of id (presuming rand()
generates evenly distrubuted values):
id = 1/percentageOfDenormalsAllowed * denormalThreshold
(I do not know the exact value of the denormalThreshold, the value at which
the FPU starts to denormal).
CODE:
// on construction:
float **denormalBuf = new float[getBlockSize()];
float id = 1.0E-20;
// inner-loop:
NOTES:
CODE:
// on construction:
float id = 1.0E-25;
// inner-loop:
is1 = *(++in1) + id;
is2 = *(++in2) + id;
id = -id;
..
NOTES:
This is the most efficient method of the lot and is also inaudible.
However, some filters will probably remove the added DC offset, thus
rendering it useless.
CODE:
// inner-loop:
** Fix-up solution
You can also walk through your filter and clamp any numbers which are close
enough to being denormal at the end of each block, as long as your blocks are
not too large. For instance, if you implement EQ using a bi-quad, you can check
the delay slots at the end of each process() call, and if any slot has a magnitude
of 10^-15 or smaller, you just clamp it to 0. This will ensure that your filter
doesn't run for long times with denormal numbers; ideally (depending on the
coefficients) it won't reach 10^-35 from the 10^-15 initial state within the time
of one block of samples.
That solution uses the least cycles, and also has the nice property of generating
absolute-0 output values for long stretches of absolute-0 input values; the others
don't.
Notes :
This very interesting paper, written by Laurent de Soras (www.ohmforce.com) has everything
you ever wanted to know about denormal numbers! And it obviously descibes how you can get
rid of them too!
DFT
int hfs=framesize/2;
float pd=pi/hfs;
for (i=0;i<framesize;i++)
{
float w=waveform[i];
int im = i-hfs;
for(int h=0;h<partials;h++)
{
float th=(pd*(h+1))*im;
aa[h]+=w*cos(th);
bb[h]+=w*sin(th);
}
}
for (int h=0;h<partials;h++)
amp[h]= sqrt(aa[h]*aa[h]+bb[h]*bb[h])/hfs;
}
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.
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
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
96kHz
b = [ 0.05265477122714 -0.09864197097385 0.04596474352090 ]
a = [ 1.00000000000000 -0.85835597216218 -0.10600020417219 ]
error +/- 0.006dB
Notes :
Buzz uses this type of synth.
For cool sounds try to use variable,
for example a=exp(-x/12000)*0.8 // x- num.samples
Code :
double DSF (double x, // input
double a, // a<1.0
double N, // N<SmplFQ/2,
double fi) // phase
{
double s1 = pow(a,N-1.0)*sin((N-1.0)*x+fi);
double s2 = pow(a,N)*sin(N*x+fi);
double s3 = a*sin(x+fi);
double s4 =1.0 - (2*a*cos(x)) +(a*a);
if (s4==0)
return 0;
else
return (sin(fi) - s3 - s2 +s1)/s4;
}
Dither code
Notes :
This is a simple implementation of highpass triangular-PDF dither (a good general-purpose
dither) with optional 2nd-order noise shaping (which lowers the noise floor by 11dB below 0.1
Fs).
The code assumes input data is in the range +1 to -1 and doesn't check for overloads!
To save time when generating dither for multiple channels you can re-use lower bits of a
previous random number instead of calling rand() again. e.g. r3=(r1 & 0x7F)<<8;
Code :
int r1, r2; //rectangular-PDF random numbers
float s1, s2; //error feedback buffers
float s = 0.5f; //set to 0.0f for no noise shaping
float w = pow(2.0,bits-1); //word length (usually bits=16)
float wi= 1.0f/w;
float d = wi / RAND_MAX; //dither amplitude (2 lsb)
float o = wi * 0.5f; //remove dc offset
float in, tmp;
int out;
s2 = s1;
s1 = in - wi * (float)out; //error
Dithering
Notes :
(see linked file)
Linked files
The noise shaping lowers the noise floor by 11dB below 5kHz (@ 44100Hz
sample rate) compared to triangular-PDF dither. The code below assumes
input data is in the range +1 to -1 and doesn't check for overloads!
To save time when generating dither for multiple channels you can do
things like this: r3=(r1 & 0x7F)<<8; instead of calling rand() again.
s2 = s1;
s1 = in - wi * (float)out; //error
--
[email protected]
http://www.maxim.abel.co.uk
Double to Int
Notes :
-Platform independant, literally. You have IEEE FP numbers, this will work, as long as your not
expecting a signed integer back larger than 16bits :)
-Will only work correctly for FP numbers within the range of [-32768.0,32767.0]
-The FPU must be in Double-Precision mode
Code :
typedef double lreal;
typedef float real;
typedef unsigned long uint32;
typedef long int32;
#if BigEndian_
#define iexp_ 0
#define iman_ 1
#else
#define iexp_ 1
#define iman_ 0
#endif //BigEndian_
// Real2Int
inline int32 Real2Int(lreal val)
{
val= val + _double2fixmagic;
return ((int32*)&val)[iman_] >> _shiftamt;
}
// Real2Int
inline int32 Real2Int(real val)
{
return Real2Int ((lreal)val);
}
__double2fixmagic dd 000000000h,042380000h
Notes :
(see linked files)
Donald Schulz's code for computing early echoes using the image-mirror method. There's an english and a german version.
Linked files
/*
From: "Donald Schulz" <[email protected]>
To: <[email protected]>
Subject: [music-dsp] Image-mirror method source code
Date: Sun, 11 Jun 2000 15:01:51 +0200
Have phun,
Donald.
*/
/****************************************************************************
*
* Early Echo Computation using image-mirror method
*
* Position of listener, 2 sound-sources, room-size may be set.
* Four early echo responses are calculated (from left sound source and
* right sound source to left and right ear). Angle with which the sound
* meets the ears is taken into account.
* The early echo response from left sound source to left ear is printed
* to screen for demonstration, the first table contains the delay times
* and the second one the weights.
*
* Program is released into the public domain.
*
* Sorry for german comments :-(
* Some frequently used german words:
* hoerpos : listening position
* breite : width
* laenge : length
* hoehe : height
* daempfung : damping
* links : left
* rechts : right
* Ohr : ear
* Winkel : angle
* Wichtung : weight
* Normierung : normalization
*
*
* If someone does some improvements on this, I (Donald, [email protected])
* would be happy to get the improved code.
*
***************************************************************************/
#include <math.h>
#include <stdio.h>
int early_pos=0;
Raummodell:
H - Hoererposition
L - Spiegelschallquellen des linken Kanales
U - Koordinatenursprung
Raum sei 11 meter breit und 5 meter lang (1 Zeichen = 1 meter)
Linker Kanal stehe bei x=2 y=4 z=?
Hoerer stehe bei x=5 y=1 z=?
___________________________________________________________________
| | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
__L___L________________L___L________________L___L________________L_
L | L | L | L | L | L | L|
| | | | | | |
| | | | | | |
| | | H | | | |
_________________________X_________________________________________
| | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
__L___L________________L___L_________________L__L________________L_
*/
main()
{
int i,j,select;
float dist_max;
float x,y,z,xref,yref,zref;
float x_pos,y_pos,z_pos;
float distance,winkel;
float wichtung;
float normierungr,normierungl;
early_l2l=(float *)malloc(early_length*sizeof(float));
early_l2r=(float *)malloc(early_length*sizeof(float));;
early_r2l=(float *)malloc(early_length*sizeof(float));;
early_r2r=(float *)malloc(early_length*sizeof(float));;
for (x=-ceil(dist_max/(2*laenge));x<=ceil(dist_max/(2*laenge));x++)
for (y=-ceil(dist_max/(2*breite));y<=ceil(dist_max/(2*breite));y++)
for (z=-ceil(dist_max/(2*hoehe));z<=ceil(dist_max/(2*hoehe));z++)
{
xref=2*x*breite;
yref=2*y*laenge;
zref=2*z*hoehe;
for (select=0;select<8;select++) /* vollstaendige Permutation */
{
if (select&1) x_pos=xref+leftposx;
else x_pos=xref-leftposx;
if (select&2) y_pos=yref+leftposy;
else y_pos=yref-leftposy;
if (select&4) z_pos=zref+leftposz;
else z_pos=zref-leftposz;
/* Jetzt steht die absolute Position der Quelle in ?_pos */
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147)<early_length)
{
/* Einfallswinkel im Bereich -Pi/2 bis Pi/2 ermitteln: */
winkel=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=links=>Verstaerkung=1 Pi=rechts=>Verstaerkung=0.22 (?) */
winkel+=3.1415926/2; wichtung = 1 - winkel/4;
}
else /* Klang kommt von hinten: */
{
winkel-=3.1415926/2; wichtung= 1 + winkel/4;
}
/* Early-echo gemaess Winkel und Entfernung gewichten: */
early_l2l[(int) (distance*147.)]+=wichtung/(pow(distance,3.1));
}
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147)<early_length)
{
/* Einfallswinkel im Bereich -Pi/2 bis Pi/2 ermitteln: */
winkel=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=links=>Verstaerkung=1 Pi=rechts=>Verstaerkung=0.22 (?) */
winkel-=3.1415926/2; wichtung = 1 + winkel/4;
}
else /* Klang kommt von hinten: */
{
winkel+=3.1415926/2; wichtung= 1 - winkel/4;
}
/* Early-echo gemaess Winkel und Entfernung gewichten: */
early_l2r[(int) (distance*147.)]+=wichtung/(pow(distance,3.1));
}
}
}
for (x=-ceil(dist_max/(2*laenge));x<=ceil(dist_max/(2*laenge));x++)
for (y=-ceil(dist_max/(2*breite));y<=ceil(dist_max/(2*breite));y++)
for (z=-ceil(dist_max/(2*hoehe));z<=ceil(dist_max/(2*hoehe));z++)
{
xref=2*x*breite;
yref=2*y*laenge;
zref=2*z*hoehe;
for (select=0;select<8;select++) /* vollstaendige Permutation */
{
if (select&1) x_pos=xref+rightposx;
else x_pos=xref-rightposx;
if (select&2) y_pos=yref+rightposy;
else y_pos=yref-rightposy;
if (select&4) z_pos=zref+rightposz;
else z_pos=zref-rightposz;
/* Jetzt steht die absolute Position der Quelle in ?_pos */
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147.)<early_length)
{
/* Einfallswinkel im Bereich -Pi/2 bis Pi/2 ermitteln: */
winkel=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=links=>Verstaerkung=1 Pi=rechts=>Verstaerkung=0.22 (?) */
winkel+=3.1415926/2; wichtung = 1 - winkel/4;
}
else /* Klang kommt von hinten: */
{
winkel-=3.1415926/2; wichtung= 1 + winkel/4;
}
/* Early-echo gemaess Winkel und Entfernung gewichten: */
early_r2l[(int) (distance*147.)]+=wichtung/(pow(distance,3.1));
}
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* Und jetzt aus berechnetem Echo die ersten early_tap_num Werte holen */
/* Erst mal e's zuruecksetzen: */
for (i=0;i<early_tap_num;i++)
{
e_values_l2l[i]=e_values_l2r[i]=0.0;
e_delays_l2l[i]=e_delays_l2r[i]=0; /* Unangenehme Speicherzugriffe vermeiden */
e_values_r2l[i]=e_values_r2r[i]=0.0;
e_delays_r2l[i]=e_delays_r2r[i]=0;
}
if (j==early_tap_num) break;
}
j=0;
for(i=0;i<early_length;i++)
{
if (early_r2l[i]!=0)
{
e_delays_r2l[j]=i;
e_values_r2l[j]=early_r2l[i];
normierungl+=early_r2l[i];
j++;
}
if (j==early_tap_num) break;
}
j=0;
for(i=0;i<early_length;i++)
{
if (early_r2r[i]!=0)
{
e_delays_r2r[j]=i;
e_values_r2r[j]=early_r2r[i];
normierungr+=early_r2r[i];
j++;
}
if (j==early_tap_num) break;
}
for (j=0;j<early_tap_num;j++)
{
e_values_l2l[j]/=normierungr;
e_values_l2r[j]/=normierungr;
e_values_r2l[j]/=normierungr;
e_values_r2r[j]/=normierungr;
}
/* Ausgeben nur der l2l-Werte fuer schnelles Reverb */
printf("int e_delays[%d]={",early_tap_num);
for (j=0;j<(early_tap_num-1);j++)
printf("%d, ",e_delays_l2l[j]);
printf("%d};\n\n",e_delays_l2l[j]);
printf("float e_values[%d]={",early_tap_num);
for (j=0;j<(early_tap_num-1);j++)
printf("%0.4f, ",e_values_l2l[j]);
printf("%0.4f};\n\n",e_values_l2l[j]);
/*
From: "Donald Schulz" <[email protected]>
To: <[email protected]>
Subject: [music-dsp] Image-mirror method source code
Date: Sun, 11 Jun 2000 15:01:51 +0200
Have phun,
Donald.
*/
/*
douglas
*/
/****************************************************************************
*
* Early Echo Computation using image-mirror method
*
* Position of listener, 2 sound-sources, room-size may be set.
* Four early echo responses are calculated (from left sound source and
* right sound source to left and right ear). Angle with which the sound
* meets the ears is taken into account.
* The early echo response from left sound source to left ear is printed
* to screen for demonstration, the first table contains the delay times
* and the second one the weights.
*
* Program is released into the public domain.
*
* Sorry for german comments :-(
* Some frequently used german words:
* hoerpos : listening position
* breite : width
* laenge : length
* hoehe : height
* daempfung : damping
* links : left
* rechts : right
* Ohr : ear
* Winkel : angle
* Wichtung : weight
* Normierung : normalization
*
*
* If someone does some improvements on this, I (Donald, [email protected])
* would be happy to get the improved code.
*
***************************************************************************/
#include <math.h>
#include <stdio.h>
int early_pos=0;
Raummodell:
H - Hoererposition
L - Spiegelschallquellen des linken Kanales
U - Koordinatenursprung
Raum sei 11 meter breit und 5 meter lang (1 Zeichen = 1 meter)
Linker Kanal stehe bei x=2 y=4 z=?
Hoerer stehe bei x=5 y=1 z=?
___________________________________________________________________
| | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
__L___L________________L___L________________L___L________________L_
L | L | L | L | L | L | L|
| | | | | | |
| | | | | | |
| | | H | | | |
_________________________X_________________________________________
| | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
__L___L________________L___L_________________L__L________________L_
*/
main()
{
int i,j,select;
float dist_max;
float x,y,z,xref,yref,zref;
float x_pos,y_pos,z_pos;
float distance,angle;
float weight;
float normalizationr,normalizationl;
early_l2l=(float *)malloc(early_length*sizeof(float));
early_l2r=(float *)malloc(early_length*sizeof(float));;
early_r2l=(float *)malloc(early_length*sizeof(float));;
early_r2r=(float *)malloc(early_length*sizeof(float));;
for (x=-ceil(dist_max/(2*length));x<=ceil(dist_max/(2*length));x++)
for (y=-ceil(dist_max/(2*width));y<=ceil(dist_max/(2*width));y++)
for (z=-ceil(dist_max/(2*height));z<=ceil(dist_max/(2*height));z++)
{
xref=2*x*width;
yref=2*y*length;
zref=2*z*height;
for (select=0;select<8;select++) /* vollstaendige Permutation */
{
if (select&1) x_pos=xref+leftposx;
else x_pos=xref-leftposx;
if (select&2) y_pos=yref+leftposy;
else y_pos=yref-leftposy;
if (select&4) z_pos=zref+leftposz;
else z_pos=zref-leftposz;
/* Jetzt steht die absolute Position der Quelle in ?_pos */
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147)<early_length)
{
/* Einfallsangle im Bereich -Pi/2 bis Pi/2 ermitteln: */
angle=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=left=>Verstaerkung=1 Pi=right=>Verstaerkung=0.22 (?) */
angle+=3.1415926/2; weight = 1 - angle/4;
}
else /* Klang kommt von hinten: */
{
angle-=3.1415926/2; weight= 1 + angle/4;
}
/* Early-echo gemaess angle und Entfernung gewichten: */
early_l2l[(int) (distance*147.)]+=weight/(pow(distance,3.1));
}
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147)<early_length)
{
/* Einfallsangle im Bereich -Pi/2 bis Pi/2 ermitteln: */
angle=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=left=>Verstaerkung=1 Pi=right=>Verstaerkung=0.22 (?) */
angle-=3.1415926/2; weight = 1 + angle/4;
}
else /* Klang kommt von hinten: */
{
angle+=3.1415926/2; weight= 1 - angle/4;
}
/* Early-echo gemaess angle und Entfernung gewichten: */
early_l2r[(int) (distance*147.)]+=weight/(pow(distance,3.1));
}
}
}
for (x=-ceil(dist_max/(2*length));x<=ceil(dist_max/(2*length));x++)
for (y=-ceil(dist_max/(2*width));y<=ceil(dist_max/(2*width));y++)
for (z=-ceil(dist_max/(2*height));z<=ceil(dist_max/(2*height));z++)
{
xref=2*x*width;
yref=2*y*length;
zref=2*z*height;
for (select=0;select<8;select++) /* vollstaendige Permutation */
{
if (select&1) x_pos=xref+rightposx;
else x_pos=xref-rightposx;
if (select&2) y_pos=yref+rightposy;
else y_pos=yref-rightposy;
if (select&4) z_pos=zref+rightposz;
else z_pos=zref-rightposz;
/* Jetzt steht die absolute Position der Quelle in ?_pos */
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* distance*147 = distance/300[m/s]*44100[ATW/s]; bei 32kHZ-> *106 */
if ((distance*147.)<early_length)
{
/* Einfallsangle im Bereich -Pi/2 bis Pi/2 ermitteln: */
angle=atan(y_pos/x_pos);
if (y_pos>0) /* Klang kommt aus vorderem Bereich: */
{ /* 0=left=>Verstaerkung=1 Pi=right=>Verstaerkung=0.22 (?) */
angle+=3.1415926/2; weight = 1 - angle/4;
}
else /* Klang kommt von hinten: */
{
angle-=3.1415926/2; weight= 1 + angle/4;
}
/* Early-echo gemaess angle und Entfernung gewichten: */
early_r2l[(int) (distance*147.)]+=weight/(pow(distance,3.1));
}
distance=sqrt(pow(x_pos,2)+pow(y_pos,2)+pow(z_pos,2));
/* Und jetzt aus berechnetem Echo die ersten early_tap_num Werte holen */
/* Erst mal e's zuruecksetzen: */
for (i=0;i<early_tap_num;i++)
{
e_values_l2l[i]=e_values_l2r[i]=0.0;
e_delays_l2l[i]=e_delays_l2r[i]=0; /* Unangenehme Speicherzugriffe vermeiden */
e_values_r2l[i]=e_values_r2r[i]=0.0;
e_delays_r2l[i]=e_delays_r2r[i]=0;
}
if (j==early_tap_num) break;
}
j=0;
for(i=0;i<early_length;i++)
{
if (early_r2l[i]!=0)
{
e_delays_r2l[j]=i;
e_values_r2l[j]=early_r2l[i];
normalizationl+=early_r2l[i];
j++;
}
if (j==early_tap_num) break;
}
j=0;
for(i=0;i<early_length;i++)
{
if (early_r2r[i]!=0)
{
e_delays_r2r[j]=i;
e_values_r2r[j]=early_r2r[i];
normalizationr+=early_r2r[i];
j++;
}
if (j==early_tap_num) break;
}
for (j=0;j<early_tap_num;j++)
{
e_values_l2l[j]/=normalizationr;
e_values_l2r[j]/=normalizationr;
e_values_r2l[j]/=normalizationr;
e_values_r2r[j]/=normalizationr;
}
/* Ausgeben nur der l2l-Werte fuer schnelles Reverb */
printf("int e_delays[%d]={",early_tap_num);
for (j=0;j<(early_tap_num-1);j++)
printf("%d, ",e_delays_l2l[j]);
printf("%d};\n\n",e_delays_l2l[j]);
printf("float e_values[%d]={",early_tap_num);
for (j=0;j<(early_tap_num-1);j++)
printf("%0.4f, ",e_values_l2l[j]);
printf("%0.4f};\n\n",e_values_l2l[j]);
References : Posted by Brahim Hamadicharef (project by Hua Zheng and Shobhit Jain)
Linked file : rev.txt (this linked file is included below)
Notes :
rev.asm
ECE320 project: Reverberation w/ parameter control from PC
Hua Zheng and Shobhit Jain
12/02/98 ~ 12/11/98
(se linked file)
Linked files
; rev.asm
;
; ECE320 project: Reverberation w/ parameter control from PC
;
; Hua Zheng and Shobhit Jain
; 12/02/98 ~ 12/11/98
;
; bugs fixed
; wrong co coefficients
; using current out point to calculate new in point
; r6 changed in set_dl (now changed to r4)
; initialize er delaylines to be 0 causes no output -- program memory
; periodic pops: getting garbage because external memory is configured to 16k
;================================================
; Initialization
;================================================
SAMPLERATE equ 48
nolist
include 'core302.asm'
list
;---------------------------
; Variables and storage setup
;---------------------------
org x:$0
ser_data ds n_para ; location of serial data chunk
delays ; default rev parameters: delay length
dc 967 ; er1 $3c7
org p:$1000
dl_er1 dsm dl_er_l ; er1 delayline
dl_er2 dsm dl_er_l ; er2 delayline
dl_er3 dsm dl_er_l ; er3 delayline
org y:$8000
dl_co1 dsm dl_co_l ; co1 delayline
dl_co2 dsm dl_co_l ; co2 delayline
dl_co3 dsm dl_co_l ; co3 delayline
dl_co4 dsm dl_co_l ; co4 delayline
dl_co5 dsm dl_co_l ; co5 delayline
dl_co6 dsm dl_co_l ; co6 delayline
org y:$F000
dl_ap1 dsm dl_ap_l ; ap1 delayline
;---------------------------
; Memory usage
;---------------------------
; P:$0000 -- $0200 core file
; P:$0200 -- $0300 progam
;
; X:$0000 -- $1BFF data 7168=149.3ms serial data, parameters,
pointers
; Y:$0000 -- $1BFF data 7168=149.3ms not used
; P:$1000 -- $4FFF data 16384=341.3ms er(85*3=255ms)
; Y:$8000 -- $FFFF data 32768=682.67ms co(80*6=480ms) and ap(10*3=30ms)
;
; X,Y:$1C00 -- $1BFF reserved for system
;
;================================================
; Main program
;================================================
org p:$200
main
;---------------------------
; Initialization
;---------------------------
move #0,x0
; move #dl_er1,r0
; move #dl_er_l,y0
; do #er_number,clear_dl_er_loop
; rep y0
; movem x0,p:(r0)+
; nop
;clear_dl_er_loop
move #dl_co1,r0
move #dl_co_l,y0
do #co_number,clear_dl_co_loop
rep y0
move x0,y:(r0)+
nop
clear_dl_co_loop
move #dl_ap1,r0
move #dl_ap_l,y0
do #ap_number,clear_dl_ap_loop
rep y0
move x0,y:(r0)+
nop
clear_dl_ap_loop
move #comb,r0
rep #co_number
move x0,x:(r0)+ ; init comb filter states
move #lpf,r0
move x0,x:(r0) ; init LPF state
; initialize SCI
movep #$0302,x:M_SCR ; R/T enable, INT disable
movep #$006a,x:M_SCCR ; SCP=0, CD=106 (9636 bps)
movep #7,x:M_PCRE
;---------------------------
; Main loop
; Register usage
; r0: delayline pointer pointer
; r1: coefficients pointer
; r2: comb filter internal state pointer
; r4,r5: used in delayline subroutine
; r6: incoming buffer pointer
; a: output of each segment
;---------------------------
main_loop:
move #dl_p,r0
move #coeffs,r1
waitdata r3,buflen,1
move x:(r3)+,a
move a,x:<in ; save input sample for mix
;---------------------------
; Early reflection
;---------------------------
; temp = in;
; for (int i=0; i<earlyRefNum; i++)
; {
; in = delayLineEarlyRef[c][i]->tick(in);
; temp += earlyRefCoeff[i] * in;
; }
; return temp;
move b,a
;---------------------------
; Comb filter
;---------------------------
; float temp1 = 0., temp2;
; for (int i=0; i<combNum; i++)
; {
; temp2 = delayLineComb[c][i]->tick
; (in + combCoeff[i] * combLastOut[c][i]);
; combLastOut[c][i] = temp2+combDamp[i]*combLastOut[c][i];
; temp1 += temp2;
; }
; return temp1 / float(combNum);
move #comb,r2
move a,y1
clr b
move #(dl_co_l-1),n6
do #co_number,co_loop
move y1,a ; a=in
move x:(r1)+,x0 ; x0=coeff
move x:(r2),y0 ; y0=lastout
mac x0,y0,a x:(r1)+,x0 ; x0=damp
jsr use_dl
move a,x1 ; a=x1=temp2
mac x0,y0,a ; a=lastout
move a,x:(r2)+
add x1,b ; b=temp1
co_loop
; asr #2,b,a
move b,a
;---------------------------
; All-pass filter
;---------------------------
; float temp1, temp2;
; for (int i=0; i<allPassNum; i++)
; {
; temp1 = delayLineAllPass[c][i]->lastOut();
; temp2 = allPassCoeff[i] * (in - temp1);
; delayLineAllPass[c][i]->tick(in + temp2);
; in = temp1 + temp2;
; }
; return in;
move #1,n0
move #0.7,x1
move #(dl_ap_l-1),n6
do #ap_number,ap_loop
move y:(r0+n0),x0 ; x0=temp1
sub x0,a ; a=in-temp1
move a,y0
mpy x1,y0,b ; b=temp2
add b,a ; a=in+temp2
jsr use_dl
add x0,b ; b=temp1+temp2
move b,a ; a=in
ap_loop
;---------------------------
; Brightness
;---------------------------
; lastout = lastout + BW * (in - lastout);
move x:<lpf,b
sub b,a x:(r1)+,x0 ; a=in-lastout, x0=bri
move a,y0
mpy x0,y0,a
add b,a
move a,x:<lpf
;---------------------------
; Mix
;---------------------------
; out = (1-mix)*in + mix*out = in + mix * (out - in);
move x:<in,y0 ; y0=in
sub y0,a x:(r1)+,x0 ; a=out-in, x0=mix
move a,y1 ; y1=out-in
mpy x0,y1,b y0,a ; b=mix*(out-in), a=in
add b,a
;---------------------------
; Spit out
;---------------------------
move a,y:(r3)+
move a,y:(r3)+
move (r3)+
;---------------------------
; Get parameters and reformat them
;---------------------------
jclr #2,x:M_SSR,main_loop ; do nothing if no new data arrived
cmp #RESET,a
jeq reformat_data ; if it's RESET, then reformat data
jmp main_loop
reformat_data:
; order of parameters:
; er1 delay, er1 coeff, er2 ..., er3 ...
; co1 delay, coeff_c, coeff_d, co2 ... , ... , co6
; ap1 delay, ap2, ap3
; brightness
; mix
move #ser_data,r0
move #delays,r1
move #coeffs,r2
do #3,format_er_loop
move x:(r0)+,a ; er delay
asr #20,a,a
; max delay 4096=2^12, max value 256=2^8, scale 256/4096=2^-4
move a0,x:(r1)+
move x:(r0)+,a ; er coeff
asr #9,a,a
move a0,x:(r2)+
format_er_loop
move #>$000001,x0
do #6,format_co_loop
move x:(r0)+,a ; co delay
asr #20,a,a ; max delay 4096=2^12
move a0,a1
;try this: asl #4,a,a
or x0,a
move a1,x:(r1)+
move x:(r0)+,a ; co coeff
asr #9,a,a
move a0,x:(r2)+
move x:(r0)+,a ; co damping
asr #9,a,a
move a0,x:(r2)+
format_co_loop
do #3,format_ap_loop
move x:(r0)+,a ; ap delay
asr #23,a,a ; max delay 528=2^9
move a0,a1
or x0,a
move a1,x:(r1)+
format_ap_loop
jsr set_dl
jmp main_loop
;================================================
; Set all delayline length subroutine
; IN: nothing
rts
;================================================
; Access delayline subroutine
; IN: in and out pointers in r4,r5
; modulo (delayline length-1) in n6
; input sample in a
; OUT: in and out pointers modulo incremented
; output sample in a
;================================================
; inputs[inPoint++] = sample;
; inPoint &= lengthm1;
; lastOutput = inputs[outPoint++];
; outPoint &= lengthm1;
; return lastOutput;
use_dl:
move n6,m4
move n6,m5
move x:(r0)+,r4 ; in point
move x:(r0)-,r5 ; out point
move a,y:(r4)+ ; queue in
move y:(r5)+,a ; queue out
move r4,x:(r0)+ ; modulo incremented in point
move r5,x:(r0)+ ; modulo incremented out point
rts
Envelope detector
Notes :
Basicaly a one-pole LP filter with different coefficients for attack and release fed by the abs() of
the signal. If you don't need different attack and decay settings, just use in->abs()->LP
Code :
//attack and release in milliseconds
float ga = (float) exp(-1/(SampleRate*attack));
float gr = (float) exp(-1/(SampleRate*release));
float envelope=0;
for(...)
{
//get your data into 'input'
EnvIn = abs(input);
Notes :
Use this if you want to do an exponential map of a parameter (mParam) to a range (mMin -
mMax).
Output is in mData...
Code :
float logmax = log10f( mMax );
float logmin = log10f( mMin );
float logdata = (mParam * (logmax-logmin)) + logmin;
Type : C code
References : Posted by musicdsp.org[AT]mindcontrol.org
Notes :
This code uses IEEE 32-bit floating point representation knowledge to quickly compute
approximations to the log2 of a value. Both functions return under-estimates of the actual value,
although the second flavour is less of an under-estimate than the first (and might be sufficient
for using in, say, a dBV/FS level meter).
0.1: -4 -3.400000
1: 0 0.000000
2: 1 1.000000
5: 2 2.250000
100: 6 6.562500
Code :
// Fast logarithm (2-based) approximation
// by Jon Watte
#include <assert.h>
#include <stdio.h>
int
main()
{
printf( "0.1: %d %f\n", floorOfLn2( 0.1 ), approxLn2( 0.1 ) );
printf( "1: %d %f\n", floorOfLn2( 1. ), approxLn2( 1. ) );
printf( "2: %d %f\n", floorOfLn2( 2. ), approxLn2( 2. ) );
printf( "5: %d %f\n", floorOfLn2( 5. ), approxLn2( 5. ) );
printf( "100: %d %f\n", floorOfLn2( 100. ), approxLn2( 100. ) );
return 0;
}
Notes :
Partial approximation of exp2 in fixed-point arithmetic. It is exactly :
[0 ; 1[ -> [0.5 ; 1[
f : x |-> 2^(x-1)
To get the full exp2 function, you have to separate the integer and fractionnal part of the
number. The integer part may be processed by bitshifting. Process the fractionnal part with the
function, and multiply the two results.
Maximum error is only 0.3 % which is pretty good for two mul ! You get also the continuity of
the first derivate.
-- Laurent
Code :
// val is a 16-bit fixed-point value in 0x0 - 0xFFFF ([0 ; 1[)
// Returns a 32-bit fixed-point value in 0x80000000 - 0xFFFFFFFF ([0.5 ;
1[)
unsigned int fast_partial_exp2 (int val)
{
unsigned int result;
__asm
{
mov eax, val
shl eax, 15 ; eax = input [31/31 bits]
or eax, 080000000h ; eax = input + 1 [32/31 bits]
mul eax
mov eax, edx ; eax = (input + 1) ^ 2 [32/30 bits]
mov edx, 2863311531 ; 2/3 [32/32 bits], rounded to +oo
mul edx ; eax = 2/3 (input + 1) ^ 2 [32/30 bits]
add edx, 1431655766 ; + 4/3 [32/30 bits] + 1
mov result, edx
}
return (result);
}
Notes :
IIRC, They're also called walsh-hadamard transforms.
Basically like Fourier, but the basis functions are squarewaves with different sequencies.
I did this for a transform data compression study a while back.
Here's some code to do a walsh hadamard transform on long ints in-place (you need to divide
by n to get transform) the order is bit-reversed at output, IIRC.
The inverse transform is the same as the forward transform (expects bit-reversed input). i.e. x =
1/n * FWHT(FWHT(x)) (x is a vector)
Code :
void inline wht_bfly (long& a, long& b)
{
long tmp = a;
a += b;
b = tmp - b;
}
return (l2);
}
////////////////////////////////////////////
// Fast in-place Walsh-Hadamard Transform //
////////////////////////////////////////////
Fast log2
float s[2];
s[0] = 0.5f;
s[1] = 0.f;
loop:
s[0] = s[0] - a*s[1];
s[1] = s[1] + a*s[0];
output_sine = s[0];
output_cosine = s[1]
Notes :
(posted by Niels Gorisse)
If you change the frequency, the amplitude rises (pitch lower) or lowers (pitch rise) a LOT I
fixed the first problem by thinking about what actually goes wrong. The answer was to
recalculate the phase for that frequency and the last value, and then continue normally.
Code :
Variables:
ip = phase of the first output sample in radians
w = freq*pi / samplerate
b1 = 2.0 * cos(w)
Init:
y1=sin(ip-w)
y2=sin(ip-2*w)
Loop:
y0 = b1*y1 - y2
y2 = y1
y1 = y0
y0 = b1*y1 - y2
y2 = b1*y0 - y1
y1 = b1*y2 - y0
Notes :
Produces a square wave -1.0f .. +1.0f.
The resulting waveform is NOT band-limited, so it's propably of not much use for syntheis. It's
rather useful for LFOs and the like, though.
Code :
Idea: use integer overflow to avoid conditional jumps.
// init:
typedef unsigned long ui32;
// loop:
(*((ui32 *)&one)) &= 0x7FFFFFFF; // mask out sign bit
(*((ui32 *)&one)) |= (intOver & 0x80000000);
intOver += intIncr;
FFT
Notes :
A paper (postscript) and some C++ source for 4 different fft algorithms, compiled by Toth Laszlo from the
Hungarian Academy of Sciences Research Group on Artificial Intelligence.
Toth says: "I've found that Sorensen's split-radix algorithm was the fastest, so I use this since then (this means that
you may as well delete the other routines in my source - if you believe my results)."
Linked files
//
// FFT library
//
// (one-dimensional complex and real FFTs for array
// lengths of 2^n)
//
// Author: Toth Laszlo ([email protected])
//
// Research Group on Artificial Intelligence
// H-6720 Szeged, Aradi vertanuk tere 1, Hungary
//
// Last modified: 97.05.29
/////////////////////////////////////////////////////////
#include <math.h>
#include <stdlib.h>
#include "pi.h"
/////////////////////////////////////////////////////////
//Gives back "i" bit-reversed where "size" is the array
//length
//currently none of the routines call it
long result,mask;
result=0;
for(;size>1;size>>=1){
mask=i&1;
i>>=1;
result<<=1;
result|=mask;
}
////////////////////////////////////////////////////////
//Bit-reverser for the Bruun FFT
//Parameters as of "bitreverse()"
result=0;
add=sizze;
while(1){
if(i!=0) {
while((i&1)==0) { i>>=1; add>>=1;}
i>>=1; add>>=1;
result+=add;
}
else {result<<=1;result+=add; return result;}
if(i!=0) {
while((i&1)==0) { i>>=1; add>>=1;}
i>>=1; add>>=1;
result-=add;
}
else {result<<=1;result-=add; return result;}
}
}
/*assembly version
long bruun_reverse(long i, long sizze){
long result;
result=0;
__asm{
mov edx,sizze
mov eax,i
mov ebx,0
l: bsf cx,eax
jz kesz1
inc cx
shr edx,cl
add ebx,edx
shr eax,cl
bsf cx,eax
jz kesz2
inc cx
shr edx,cl
sub ebx,edx
shr eax,cl
jmp l
kesz1:
shl ebx,1
add ebx,edx
jmp vege
kesz2:
shl ebx,1
sub ebx,edx
/////////////////////////////////////////////////////////
// Decimation-in-freq radix-2 in-place butterfly
// data: array of doubles:
// re(0),im(0),re(1),im(1),...,re(size-1),im(size-1)
// it means that size=array_length/2 !!
//
// suggested use:
// intput in normal order
// output in bit-reversed order
//
// Source: Rabiner-Gold: Theory and Application of DSP,
// Prentice Hall,1978
long angle,astep,dl;
double xr,yr,xi,yi,wr,wi,dr,di,ang;
double *l1, *l2, *end, *ol2;
astep=1;
end=data+size+size;
for(dl=size;dl>1;dl>>=1,astep+=astep){
l1=data;
l2=data+dl;
for(;l2<end;l1=l2,l2=l2+dl){
ol2=l2;
for(angle=0;l1<ol2;l1+=2,l2+=2){
ang=2*pi*angle/size;
wr=cos(ang);
wi=-sin(ang);
xr=*l1+*l2;
xi=*(l1+1)+*(l2+1);
dr=*l1-*l2;
di=*(l1+1)-*(l2+1);
yr=dr*wr-di*wi;
yi=dr*wi+di*wr;
*(l1)=xr;*(l1+1)=xi;
*(l2)=yr;*(l2+1)=yi;
angle+=astep;
}
}
}
}
/////////////////////////////////////////////////////////
// Decimation-in-time radix-2 in-place inverse butterfly
// data: array of doubles:
// re(0),im(0),re(1),im(1),...,re(size-1),im(size-1)
// it means that size=array_length/2 !!
//
// suggested use:
// intput in bit-reversed order
// output in normal order
//
// Source: Rabiner-Gold: Theory and Application of DSP,
// Prentice Hall,1978
long angle,astep,dl;
double xr,yr,xi,yi,wr,wi,dr,di,ang;
double *l1, *l2, *end, *ol2;
astep=size>>1;
end=data+size+size;
for(dl=2;astep>0;dl+=dl,astep>>=1){
l1=data;
l2=data+dl;
for(;l2<end;l1=l2,l2=l2+dl){
ol2=l2;
for(angle=0;l1<ol2;l1+=2,l2+=2){
ang=2*pi*angle/size;
wr=cos(ang);
wi=sin(ang);
xr=*l1;
xi=*(l1+1);
yr=*l2;
yi=*(l2+1);
dr=yr*wr-yi*wi;
di=yr*wi+yi*wr;
*(l1)=xr+dr;*(l1+1)=xi+di;
*(l2)=xr-dr;*(l2+1)=xi-di;
angle+=astep;
}
}
}
}
/////////////////////////////////////////////////////////
// data shuffling into bit-reversed order
// data: array of doubles:
// re(0),im(0),re(1),im(1),...,re(size-1),im(size-1)
// it means that size=array_length/2 !!
//
// Source: Rabiner-Gold: Theory and Application of DSP,
// Prentice Hall,1978
long i,j,k,l,m;
double re,im;
l=size-1;
m=size>>1;
for (i=0,j=0; i<l ; i++){
if (i<j){
re=data[j+j]; im=data[j+j+1];
data[j+j]=data[i+i]; data[j+j+1]=data[i+i+1];
data[i+i]=re; data[i+i+1]=im;
}
k=m;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
}
/////////////////////////////////////////////////////////
// used by realfft
// parameters as above
//
// Source: Brigham: The Fast Fourier Transform
// Prentice Hall, ?
double xr,yr,xi,yi,wr,wi,dr,di,ang,astep;
double *l1, *l2;
l1=data;
l2=data+size+size-2;
xr=*l1;
xi=*(l1+1);
*l1=xr+xi;
*(l1+1)=xr-xi;
l1+=2;
astep=pi/size;
for(ang=astep;l1<=l2;l1+=2,l2-=2,ang+=astep){
xr=(*l1+*l2)/2;
yi=(-(*l1)+(*l2))/2;
yr=(*(l1+1)+*(l2+1))/2;
xi=(*(l1+1)-*(l2+1))/2;
wr=cos(ang);
wi=-sin(ang);
dr=yr*wr-yi*wi;
di=yr*wi+yi*wr;
*l1=xr+dr;
*(l1+1)=xi+di;
*l2=xr-dr;
*(l2+1)=-xi+di;
}
}
/////////////////////////////////////////////////////////
// used by inverse realfft
// parameters as above
//
// Source: Brigham: The Fast Fourier Transform
// Prentice Hall, ?
double xr,yr,xi,yi,wr,wi,dr,di,ang,astep;
double *l1, *l2;
l1=data;
l2=data+size+size-2;
xr=(*l1)/2;
xi=(*(l1+1))/2;
*l1=xr+xi;
*(l1+1)=xr-xi;
l1+=2;
astep=pi/size;
for(ang=astep;l1<=l2;l1+=2,l2-=2,ang+=astep){
xr=(*l1+*l2)/2;
yi=-(-(*l1)+(*l2))/2;
yr=(*(l1+1)+*(l2+1))/2;
xi=(*(l1+1)-*(l2+1))/2;
wr=cos(ang);
wi=-sin(ang);
dr=yr*wr-yi*wi;
di=yr*wi+yi*wr;
*l2=xr+dr;
*(l1+1)=xi+di;
*l1=xr-dr;
*(l2+1)=-xi+di;
}
}
/////////////////////////////////////////////////////////
// in-place Radix-2 FFT for complex values
// data: array of doubles:
// re(0),im(0),re(1),im(1),...,re(size-1),im(size-1)
// it means that size=array_length/2 !!
//
// output is in similar order, normalized by array length
//
// Source: see the routines it calls ...
dif_butterfly(data,size);
unshuffle(data,size);
end=data+size+size;
for(l=data;l<end;l++){*l=*l/size;};
}
/////////////////////////////////////////////////////////
// in-place Radix-2 inverse FFT for complex values
unshuffle(data,size);
inverse_dit_butterfly(data,size);
}
/////////////////////////////////////////////////////////
// in-place Radix-2 FFT for real values
// (by the so-called "packing method")
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),re(1),im(1),re(2),im(2),...,re(size/2-1),im(size/2-1)
// normalized by array length
//
// Source: see the routines it calls ...
size>>=1;
dif_butterfly(data,size);
unshuffle(data,size);
realize(data,size);
end=data+size+size;
div=size+size;
for(l=data;l<end;l++){*l=*l/div;};
}
/////////////////////////////////////////////////////////
// in-place Radix-2 inverse FFT for real values
// (by the so-called "packing method")
// data: array of doubles:
// re(0),re(size/2),re(1),im(1),re(2),im(2),...,re(size/2-1),im(size/2-1)
//
// output:
// re(0),re(1),re(2),...,re(size-1)
size>>=1;
unrealize(data,size);
unshuffle(data,size);
inverse_dit_butterfly(data,size);
end=data+size+size;
for(l=data;l<end;l++){*l=(*l)*2;};
}
/////////////////////////////////////////////////////////
// Bruun FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),...,re(i),im(i)... pairs in
// "bruun-reversed" order
// normalized by array length
//
// Source:
// Bruun: z-Transform DFT Filters and FFT's
// IEEE Trans. ASSP, ASSP-26, No. 1, February 1978
//
// Comments:
// (this version is implemented in a manner that every
// cosine is calculated only once;
// faster than the other version (see next)
end=data+size;
//first filterings, when there're only two taps
size>>=1;
dl=size;
dl2=dl/2;
for(;dl>1;dl>>=1,dl2>>=1){
l0=data;
l3=data+dl;
for(i=0;i<dl;i++){
d0=*l0;
d2=*l3;
*l0=d0+d2;
*l3=d0-d2;
l0++;
l3++;
}
}
l0=data;l1=data+1;
d0=*l0;d1=*l1;
*l0=d0+d1;*l1=d0-d1;
l1+=2;
*l1=-(*l1);
p4=pi/(2*size);
j=1;
kk=1;
dl_o=size/2;
dl2_o=size/4;
while(dl_o>1){
for(k=0;k<kk;k++){
c2=p4*bruun_reverse(j,size);
c=2*cos(c2);
c2=2*sin(c2);
dl=dl_o;
dl2=dl2_o;
for(;dl>1;dl>>=1,dl2>>=1){
l0=data+((dl*j)<<1);
l1=l0+dl2;l2=l0+dl;l3=l1+dl;
for(i=0;i<dl2;i++){
d1=(*l1)*c;
d2=(*l2)*c;
d3=*l3+(*l1);
d0=*l0+(*l2);
*l0=d0+d1;
*l1=d3+d2;
*l2=d0-d1;
*l3=d3-d2;
l0++;
l1++;
l2++;
l3++;
}
}
//real conversion
l3-=4;
*l3=*l3-c*(*l0)/2;
*l0=-c2*(*l0)/2;
*l1=*l1+c*(*l2)/2;
*l2=-c2*(*l2)/2;
j++;
}
dl_o>>=1;
dl2_o>>=1;
kk<<=1;
}
/////////////////////////////////////////////////////////
// Bruun FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),re(1),im(1),re(2),im(2),...,re(size/2-1),im(size/2-1)
// normalized by array length
//
// Source: see the routines it calls ...
void realfft_bruun_unshuffled(double *data, long size){
double *data2;
long i,j,k;
realfft_bruun(data,size);
//unshuffling - cannot be done in-place (?)
data2=(double *)malloc(size*sizeof(double));
for(i=1,k=size>>1;i<k;i++){
j=bruun_reverse(i,k);
data2[j+j]=data[i+i];
data2[j+j+1]=data[i+i+1];
}
for(i=2;i<size;i++) data[i]=data2[i];
free(data2);
}
/////////////////////////////////////////////////////////
// Bruun FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),...,re(i),im(i)... pairs in
// "bruun-reversed" order
end=data+size;
p4=pi/(size);
size>>=1;
dl=size;
dl2=dl/2;
//first filtering, when there're only two taps
for(;dl>1;dl>>=1,dl2>>=1){
l0=data;
l3=data+dl;
for(i=0;i<dl;i++){
d0=*l0;
d2=*l3;
*l0=d0+d2;
*l3=d0-d2;
l0++;
l3++;
}
//the remaining filterings
j=1;
while(l3<end){
l0=l3;l1=l0+dl2;l2=l0+dl;l3=l1+dl;
c=2*cos(p4*bruun_reverse(j,size));
for(i=0;i<dl2;i++){
d0=*l0;
d1=*l1;
d2=*l2;
d3=*l3;
*l0=d0+c*d1+d2;
*l2=d0-c*d1+d2;
*l1=d1+c*d2+d3;
*l3=d1-c*d2+d3;
l0++;
l1++;
l2++;
l3++;
}
j++;
}
}
double *data2;
long i,j,k;
realfft_bruun2(data,size);
//unshuffling - cannot be done in-place (?)
data2=(double *)malloc(size*sizeof(double));
for(i=1,k=size>>1;i<k;i++){
j=bruun_reverse(i,k);
data2[j+j]=data[i+i];
data2[j+j+1]=data[i+i+1];
}
for(i=2;i<size;i++) data[i]=data2[i];
free(data2);
}*/
/////////////////////////////////////////////////////////
// Sorensen in-place split-radix FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
// normalized by array length
//
// Source:
// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987
long i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8;
double t1,t2,t3,t4,t5,t6,a3,ss1,ss3,cc1,cc3,a,e,sqrt2;
sqrt2=sqrt(2.0);
n4=n-1;
//data shuffling
for (i=0,j=0,n2=n/2; i<n4 ; i++){
if (i<j){
t1=data[j];
data[j]=data[i];
data[i]=t1;
}
k=n2;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
/*----------------------*/
data[i0]=t1+data[i1];
data[i1]=t1-data[i1];
}
id<<=1;
i0=id-2;
id<<=1;
} while ( i0<n4 );
/*----------------------*/
//L shaped butterflies
n2=2;
for(k=n;k>2;k>>=1){
n2<<=1;
n4=n2>>2;
n8=n2>>3;
e = 2*pi/(n2);
i1=0;
id=n2<<1;
do{
for (; i1<n; i1+=id){
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
t1=data[i4]+data[i3];
data[i4]-=data[i3];
data[i3]=data[i1]-t1;
data[i1]+=t1;
if (n4!=1){
i0=i1+n8;
i2+=n8;
i3+=n8;
i4+=n8;
t1=(data[i3]+data[i4])/sqrt2;
t2=(data[i3]-data[i4])/sqrt2;
data[i4]=data[i2]-t1;
data[i3]=-data[i2]-t1;
data[i2]=data[i0]-t2;
data[i0]+=t2;
}
}
id<<=1;
i1=id-n2;
id<<=1;
} while ( i1<n );
a=e;
for (j=2; j<=n8; j++){
a3=3*a;
cc1=cos(a);
ss1=sin(a);
cc3=cos(a3);
ss3=sin(a3);
a=j*e;
i=0;
id=n2<<1;
do{
for (; i<n; i+=id){
i1=i+j-1;
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
i5=i+n4-j+1;
i6=i5+n4;
i7=i6+n4;
i8=i7+n4;
t1=data[i3]*cc1+data[i7]*ss1;
t2=data[i7]*cc1-data[i3]*ss1;
t3=data[i4]*cc3+data[i8]*ss3;
t4=data[i8]*cc3-data[i4]*ss3;
t5=t1+t3;
t6=t2+t4;
t3=t1-t3;
t4=t2-t4;
t2=data[i6]+t6;
data[i3]=t6-data[i6];
data[i8]=t2;
t2=data[i2]-t3;
data[i7]=-data[i2]-t3;
data[i4]=t2;
t1=data[i1]+t5;
data[i6]=data[i1]-t5;
data[i1]=t1;
t1=data[i5]+t4;
data[i5]-=t4;
data[i2]=t1;
}
id<<=1;
i=id-n2;
id<<=1;
} while(i<n);
}
}
/////////////////////////////////////////////////////////
// Sorensen in-place split-radix FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),re(1),im(1),re(2),im(2),...,re(size/2-1),im(size/2-1)
// normalized by array length
//
// Source:
// Source: see the routines it calls ...
double *data2;
long i,j;
realfft_split(data,n);
//unshuffling - not in-place
data2=(double *)malloc(n*sizeof(double));
j=n/2;
data2[0]=data[0];
data2[1]=data[j];
for(i=1;i<j;i++) {data2[i+i]=data[i];data2[i+i+1]=data[n-i];}
for(i=0;i<n;i++) data[i]=data2[i];
free(data2);
}
/////////////////////////////////////////////////////////
// Sorensen in-place inverse split-radix FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
//
// output:
// re(0),re(1),re(2),...,re(size-1)
// NOT normalized by array length
//
// Source:
// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987
long i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8,n1;
double t1,t2,t3,t4,t5,a3,ss1,ss3,cc1,cc3,a,e,sqrt2;
sqrt2=sqrt(2.0);
n1=n-1;
n2=n<<1;
for(k=n;k>2;k>>=1){
id=n2;
n2>>=1;
n4=n2>>2;
n8=n2>>3;
e = 2*pi/(n2);
i1=0;
do{
for (; i1<n; i1+=id){
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
t1=data[i1]-data[i3];
data[i1]+=data[i3];
data[i2]*=2;
data[i3]=t1-2*data[i4];
data[i4]=t1+2*data[i4];
if (n4!=1){
i0=i1+n8;
i2+=n8;
i3+=n8;
i4+=n8;
t1=(data[i2]-data[i0])/sqrt2;
t2=(data[i4]+data[i3])/sqrt2;
data[i0]+=data[i2];
data[i2]=data[i4]-data[i3];
data[i3]=2*(-t2-t1);
data[i4]=2*(-t2+t1);
}
}
id<<=1;
i1=id-n2;
id<<=1;
} while ( i1<n1 );
a=e;
for (j=2; j<=n8; j++){
a3=3*a;
cc1=cos(a);
ss1=sin(a);
cc3=cos(a3);
ss3=sin(a3);
a=j*e;
i=0;
id=n2<<1;
do{
for (; i<n; i+=id){
i1=i+j-1;
i2=i1+n4;
i3=i2+n4;
i4=i3+n4;
i5=i+n4-j+1;
i6=i5+n4;
i7=i6+n4;
i8=i7+n4;
t1=data[i1]-data[i6];
data[i1]+=data[i6];
t2=data[i5]-data[i2];
data[i5]+=data[i2];
t3=data[i8]+data[i3];
data[i6]=data[i8]-data[i3];
t4=data[i4]+data[i7];
data[i2]=data[i4]-data[i7];
t5=t1-t4;
t1+=t4;
t4=t2-t3;
t2+=t3;
data[i3]=t5*cc1+t4*ss1;
data[i7]=-t4*cc1+t5*ss1;
data[i4]=t1*cc3-t2*ss3;
data[i8]=t2*cc3+t1*ss3;
}
id<<=1;
i=id-n2;
id<<=1;
} while(i<n1);
}
}
/*----------------------*/
i0=0;
id=4;
do{
for (; i0<n1; i0+=id){
i1=i0+1;
t1=data[i0];
data[i0]=t1+data[i1];
data[i1]=t1-data[i1];
}
id<<=1;
i0=id-2;
id<<=1;
} while ( i0<n1 );
/*----------------------*/
//data shuffling
for (i=0,j=0,n2=n/2; i<n1 ; i++){
if (i<j){
t1=data[j];
data[j]=data[i];
data[i]=t1;
}
k=n2;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
}
/////////////////////////////////////////////////////////
// Sorensen in-place radix-2 FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
// normalized by array length
//
// Source:
// Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
// IEEE Trans. ASSP, ASSP-35, No. 6, June 1987
n4=n-1;
//data shuffling
for (i=0,j=0,n2=n/2; i<n4 ; i++){
if (i<j){
xt=data[j];
data[j]=data[i];
data[i]=xt;
}
k=n2;
while (k<=j){
j-=k;
k>>=1;
}
j+=k;
}
/* -------------------- */
for (i=0; i<n; i += 2)
{
xt = data[i];
data[i] = xt + data[i+1];
data[i+1] = xt - data[i+1];
}
/* ------------------------ */
n2 = 1;
for (k=n;k>2;k>>=1){
n4 = n2;
n2 = n4 << 1;
n1 = n2 << 1;
e = 2*pi/(n1);
for (i=0; i<n; i+=n1){
xt = data[i];
data[i] = xt + data[i+n2];
data[i+n2] = xt-data[i+n2];
data[i+n4+n2] = -data[i+n4+n2];
a = e;
n3=n4-1;
for (j = 1; j <=n3; j++){
i1 = i+j;
i2 = i - j + n2;
i3 = i1 + n2;
i4 = i - j + n1;
cc = cos(a);
ss = sin(a);
a += e;
t1 = data[i3] * cc + data[i4] * ss;
t2 = data[i3] * ss - data[i4] * cc;
data[i4] = data[i2] - t2;
data[i3] = -data[i2] - t2;
data[i2] = data[i1] - t1;
data[i1] += t1;
}
}
}
/////////////////////////////////////////////////////////
// Sorensen in-place split-radix FFT for real values
// data: array of doubles:
// re(0),re(1),re(2),...,re(size-1)
//
// output:
// re(0),re(size/2),re(1),im(1),re(2),im(2),...,re(size/2-1),im(size/2-1)
// normalized by array length
//
// Source:
// Source: see the routines it calls ...
double *data2;
long i,j;
Notes :
(see linkfile)
Float to int
Notes :
intel only
Code :
int truncate(float flt)
{
int i;
static const double half = 0.5f;
_asm
{
fld flt
fsub half
fistp i
}
return i
}
Notes :
intel only
Code :
void f2short(float *fptr,short *iptr,int n)
{
_asm {
mov ebx,n
mov esi,fptr
mov edi,iptr
lea ebx,[esi+ebx*4] ; ptr after last
mov edx,0x80008000 ; damn endianess confuses...
mov ecx,0x4b004b00 ; damn endianess confuses...
mov eax,[ebx] ; get last value
push eax
mov eax,0x4b014B01
mov [ebx],eax ; mark end
mov ax,[esi+2]
jmp startf2slp
Formant filter
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;
}
Gaussian dithering
Type : Dithering
References : Posted by Aleksey Vaneev (picoder[AT]mail[DOT]ru)
Notes :
It is a more sophisticated dithering than simple RND. It gives the most low noise floor for the
whole spectrum even without noise-shaping. You can use as big N as you can afford (it will not
hurt), but 4 or 5 is generally enough.
Code :
Basically, next value is calculated this way (for RND going from -0.5 to
0.5):
dither = (RND+RND+...+RND) / N.
\ /
\_________/
N times
Notes :
Code I use sometimes, but don't remember where I ripped it from.
- Alexey Menshikov
Code :
#define ranf() ((float) rand() / (float) RAND_MAX)
if (pass)
{
y1 = y2;
} else {
do {
x1 = 2.0f * ranf () - 1.0f;
x2 = 2.0f * ranf () - 1.0f;
w = x1 * x1 + x2 * x2;
} while (w >= 1.0f);
Notes :
SOURCE:
Steven W. Smith:
The Scientist and Engineer's Guide to Digital Signal Processing
http://www.dspguide.com
Code :
#define PI 3.1415926536f
Guitar feedback
Notes :
It is fairly simple to simulate guitar feedback with a simple Karplus-Strong algorithm (this was
described in a CMJ article in the early 90's):
Code :
Run the output of the Karplus-Strong delay lines into a nonlinear shaping
function for distortion (i.e. 6 parallel delay lines for 6 strings, going
into 1 nonlinear shaping function that simulates an overdriven amplifier,
fuzzbox, etc.);
Run part of the output into a delay line, to simulate the distance from the
amplifier to the "strings";
The delay line feeds back into the Karplus-Strong delay lines. By
controlling the amount of the output fed into the delay line, and the
length of the delay line, you can control the intensity and pitch of the
feedback note.
Hermite interpollation
Notes :
These are all different ways to do the same thing : hermite interpollation. Try'm all and
benchmark.
Code :
// original
inline float hermite1(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float c2 = y0 - 2.5f * y1 + 2.f * y2 - 0.5f * y3;
float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
// james mccartney
inline float hermite2(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
float c2 = y0 - y1 + c1 - c3;
// james mccartney
inline float hermite3(float x, float y0, float y1, float y2, float y3)
{
// 4-point, 3rd-order Hermite (x-form)
float c0 = y1;
float c1 = 0.5f * (y2 - y0);
float y0my1 = y0 - y1;
float c3 = (y1 - y2) + 0.5f * (y3 - y0my1 - y2);
float c2 = y0my1 + c1 - c3;
// laurent de soras
inline float hermite4(float frac_pos, float xm1, float x0, float x1, float
x2)
{
const float c = (x1 - xm1) * 0.5f;
const float v = x0 - x1;
const float w = c + v;
const float a = w + v + (x2 - x0) * 0.5f;
const float b_neg = w + a;
level = beglevel;
slope = 4.0 * amp * (rdur - rdur2);
curve = -8.0 * amp * rdur2;
...
Notes :
Simple implementation of a lock free FIFO.
Linked files
#include <vector>
#include <exception>
using std::vector;
using std::exception;
T get (void)
{
if (readidx == writeidx)
throw runtime_error ("underrun");
T result = buffer[readidx];
return result;
}
if (newidx == readidx)
throw runtime_error ("overrun");
buffer[writeidx] = datum;
writeidx = newidx;
}
private:
Notes :
use add_value with all values which enter the look-ahead area,
and remove_value with all value which leave this area. to get
the maximum value in the look-ahead area, use get_max_value.
in the very beginning initialize the table with zeroes.
Code :
void lookup_add(unsigned section, unsigned size, unsigned value)
{
if (section==value)
lookup[section]++;
else
{
size >>= 1;
if (value>section)
{
lookup[section]++;
lookup_add(section+size,size,value);
}
else
lookup_add(section-size,size,value);
}
}
}
}
unsigned get_max_value()
{
return lookup_getmax(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1);
}
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];
/* Check clipping */
if (temp > 32767)
{
temp = 32767;
}
else if (temp < -32768)
{
temp = -32768;
}
Linked files
LP and HP filter
Lowpass:
c = 1.0 / tan(pi * f / sample_rate);
Hipass:
c = tan(pi * f / sample_rate);
Notes :
Amplitude and phase calculation of IIR equation
run at sample rate "sampleRate" at frequency "F".
AMPLITUDE
-----------
cf_mag(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
PHASE
-----------
cf_phi(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------
Do not try to reverse engineer these formulae - they don't give any sense
other than they are derived from transfer function, and they work. :)
Code :
/*
C file can be downloaded from
http://www.yohng.com/dsp/cfsmp.c
*/
Notes :
(see linkfile)
Notes :
You can easily estimate the error by evaluating the actual function and
evaluating your interpolator at each of the mid-points between your
samples. The absolute difference between these values, over the absolute
value of the "correct" value, is your relative error. log10 of your relative
error times 20 is an estimate of your quantization noise in dB. Example:
You have a table for every 0.5 "index units". The value at index unit 72.0
is 0.995 and the value at index unit 72.5 is 0.999. The interpolated value
at index 72.25 is 0.997. Suppose the actual function value at that point was
0.998; you would have an error of 0.001 which is a relative error of 0.001002004..
log10(error) is about -2.99913, which times 20 is about -59.98. Thus, that's
your quantization noise at that position in the table. Repeat for each pair of
samples in the table.
Note: I said "quantization noise" not "aliasing noise". The aliasing noise will,
as far as I know, only happen when you start up-sampling without band-limiting
and get frequency aliasing (wrap-around), and thus is mostly independent of
what specific interpolation mechanism you're using.
Millimeter to DB (faders...)
Notes :
These two functions reproduce a traditional professional
mixer fader taper.
Code :
float MMtoDB(float mm)
{
float db;
mm = 100. - mm;
mm = 100. - mm;
return mm;
}
Moog VCF
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)
y1=y2=y3=y4=oldx=oldy1=oldy2=oldy3=0;
//Loop
//--Inverted feed back for corner peaking
x = input - r*y4;
oldx = x;
oldy1 = y1;
oldy2 = y2;
oldy3 = y3;
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 [email protected] July 2000
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));
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);
Notes :
in[x] and out[x] are member variables, init to 0.0 the controls:
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;
}
Notes :
This is an implemetation of a 9th order noise shaping & dithering class, that runs quite fast (it has one function that uses Intel
x86 assembly, but you can replace it with a different rounding function if you are running on a non-Intel platform).
_aligned_malloc and _aligned_free require the MSVC++ Processor Pack, available from www.microsoft.com. You can replace
them with "new" and "delete," but allocating aligned memory seems to make it run faster. Also, you can replace ZeroMemory
with a memset that sets the memory to 0 if you aren't using Win32.
Input should be floats from -32768 to 32767 (processS will clip at these points for you, but clipping is bad when you are trying
to convert floats to shorts). Note to reviewer - it would probably be better if you put the code in a file such as NSDither.h and
have a link to it - it's rather long.
Linked files
#pragma once
#include <malloc.h>
#include <rounding.h>
// F-weighted
//static const float or3Fc[3] = {1.623f, -0.982f, 0.109f};
static const float or9Fc[9] = {2.412f, -3.370f, 3.937f, -4.174f, 3.353f, -2.205f,
1.281f, -0.569f, 0.0847f};
// modified-E weighted
//static const float or2MEc[2] = {1.537f, -0.8367f};
//static const float or3MEc[3] = {1.652f, -1.049f, 0.1382f};
//static const float or9MEc[9] = {1.662f, -1.263f, 0.4827f, -0.2913f, 0.1268f, -
0.1124f, 0.03252f, -0.01265f, -0.03524f};
// improved-E weighted
//static const float or5IEc[5] = {2.033f, -2.165f, 1.959f, -1.590f, 0.6149f};
//static const float or9IEc[9] = {2.847f, -4.685f, 6.214f, -7.184f, 6.639f, -5.032f,
3.263f, -1.632f, 0.4191f};
// Much faster than C's % operator (anyway, x will never be > 2*n in this case so
this is very simple and fast)
// Tell the compiler to try to inline as much as possible, since it makes it run so
much faster since these functions get called at least once per sample
__inline my_mod(int x, int n)
{
if(x > n) x-=n;
return x;
}
{
int r;
__asm {
fld f
fistp r
}
return r;
}
a = a * 140359821 + 1;
return a * (1.0f / 0xFFFFFFFF);
}
class NS9dither16
{
public:
NS9dither16();
~NS9dither16();
short processS(float samp);
int processI(float samp);
void reset();
private:
int order;
int HistPos;
float* c; // Coeffs
float* EH; // Error History
};
__inline NS9dither16::NS9dither16()
{
order = 9;
// Start at top of error history - you can make it start anywhere from 0-8 if
you really want to
HistPos=8;
__inline NS9dither16::~NS9dither16()
{
//if(c) delete [] c; // if you don't have _aligned_free
//if(EH) delete [] EH; // if you don't have _aligned_free
if(c) _aligned_free(c); // don't really need "if," since it is OK to free
null pointer, but still...
if(EH) _aligned_free(EH);
}
// Force inline because VC++.NET doesn't inline for some reason (VC++ 6 does)
__forceinline short NS9dither16::processS(float samp)
{
int output;
return ltos(output);
}
{
//samp -= c[x] * EH[(HistPos+x) % order];
samp -= c[x] * EH[HistPos+x];
}*/
// Unrolled loop for faster execution
/*samp -= c[0]*EH[HistPos] + c[1]*EH[HistPos+1] + c[2]*EH[HistPos+2] +
c[3]*EH[HistPos+3] + c[4]*EH[HistPos+4] + c[5]*EH[HistPos+5]
+
c[6]*EH[HistPos+6] + c[7]*EH[HistPos+7] +
c[8]*EH[HistPos+8];*/
// This arrangement seems to execute 3 clock cycles faster on a P-III
samp -= c[8]*EH[HistPos+8] + c[7]*EH[HistPos+7] + c[6]*EH[HistPos+6] +
c[5]*EH[HistPos+5] + c[4]*EH[HistPos+4] + c[3]*EH[HistPos+3]
+
c[2]*EH[HistPos+1] + c[1]*EH[HistPos+1] + c[0]*EH[HistPos];
return output;
}
Notch filter
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);
/* 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;
/* Check clipping */
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
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
Where
out[n] = in[n]*a0 + in[n-1]*a1 + out[n-1]*b1;
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
Notes :
This formula can be found from a patent related to parallel combs structure. The formula places
the first echoes coming out of parallel combs to uniformly distributed sequence. If T_ ,...,T_n
are the delay lines in increasing order, the formula can be derived by setting T_(k-1)/T_k =
Constant and T_n/(2*T_1) = Constant, where 2*T_1 is the echo coming just after the echo T_n.
I figured this out myself as it is not told in the patent. The formula is not the best which one can
come up. I use a search method to find echo sequences which are uniform enough for long
enough time. The formula is uniform for a short time only.
The formula doesn't work good for series allpass and FDN structures, for which a similar
formula can be derived with the same idea. The search method works for these structures as
well.
Notes :
This code shows what the difference is betwee FM and PM.
The code is NOT optimised, nor should it be used like this.
It is an <b>EXAMPLE</b>
Linked files
////////////////////////////////////////////////////
//
// this code was NEVER MEANT TO BE USED.
//
// use as EXPLANATION ONLY for the difference between
// Phase Modulation and Frequency Modulation.
// there are MANY ways to speed this code up.
//
// [email protected] | [email protected]
//
// ps:
// we use the 'previous' value of the phase in all the algo's to make sure that
// the first call to the getSampleXX() function returns the wave at phase 'zero'
//
////////////////////////////////////////////////////
#include "math.h";
#define Pi 3.141592f
class SimpleOscillator
{
SimpleOscillator(const float sampleRate = 44100.f, const long tableSize =
4096)
{
this->tableSize = tableSize;
this->sampleRate = sampleRate;
phase = 0.f;
makeTable();
}
~SimpleOscillator()
{
delete [] table;
}
return lookup(lookupPhase);
}
// frequency modulation
// the fm input should be in HZ.
//
// example:
// osc1.getSampleFM(440.f, osc2.getSample(0.5f) * 5.f )
// would give a signal where the frequency of the signal is
// modulated between 435hz and 445hz at a 0.5hz rate
//
float generateSampleFM(const float frequency, const float fm)
{
float lookupPhase = phase;
return lookup(lookupPhase);
}
// phase modulation
//
// a phase mod value of 1.f will increase the "phase" of the wave by a full
cycle
// i.e. calling getSamplePM(440.f,1.f) will result in the "same" wave as
getSamplePM(440.f,0.f)
//
float generateSamplePM(const float frequency, const float pm)
{
float lookupPhase = phase + (float)tableSize * pm;
wrap(lookupPhase)
return lookup(lookupPhase);
}
//
float lookup(const float phase)
{
return table[(long)phase];
}
// wrap around
//
void wrap(float &in)
{
while(in < 0.f)
in += (float)tableSize;
return in;
}
private:
float sampleRate;
float phase;
float *table;
long tableSize;
void makeTable()
{
table = new float[tableSize];
for(long i=0;i<tableSize;i++)
{
float x = Pi * 2.f * (float)i / (float)tableSize;
table[i] = (float)sin(x);
}
}
Phaser code
Notes :
(see linked file)
Linked files
/*
Date: Mon, 24 Aug 1998 07:02:40 -0700
Reply-To: music-dsp
Originator: [email protected]
Sender: music-dsp
Precedence: bulk
From: "Ross Bencina" <[email protected]>
To: Multiple recipients of list <music-dsp>
Subject: Re: Phaser revisited [code included]
X-Comment: Music Hackers Unite! http://shoko.calarts.edu/~glmrboy/musicdsp/music-
dsp.html
Status: RO
Hi again,
Ross B.
*/
/*
class: Phaser
implemented by: Ross Bencina <[email protected]>
date: 24/8/98
Thanks to:
*/
class Phaser{
public:
Phaser() //initialise to some usefull defaults...
: _fb( .7f )
, _lfoPhase( 0.f )
, _depth( 1.f )
, _zm1( 0.f )
{
Range( 440.f, 1600.f );
Rate( .5f );
}
//calculate output
float y = _alps[0].Update(
_alps[1].Update(
_alps[2].Update(
_alps[3].Update(
_alps[4].Update(
_alps[5].Update( inSamp + _zm1 * _fb ))))));
_zm1 = y;
return y;
}
private:
float _a1, _zm1;
};
AllpassDelay _alps[6];
float _zm1;
};
Notes :
(see linked file)
Linked files
---
[email protected]
http://www.abel.co.uk/~maxim/
Polyphase Filters
Linked files
x0=0.0;
x1=0.0;
x2=0.0;
y0=0.0;
y1=0.0;
y2=0.0;
};
CAllPassFilter::~CAllPassFilter()
{
};
//shuffle outputs
y2=y1;
y1=y0;
//allpass filter 1
const double output=x2+((input-y2)*a);
y0=output;
return output;
};
allpassfilter[i]=new CAllPassFilter(coefficient[i]);
}
numfilters=N;
};
CAllPassFilterCascade::~CAllPassFilterCascade()
{
delete[] allpassfilter;
};
int i=0;
do
{
output=allpassfilter[i]->process(output);
i++;
} while (i<numfilters);
return output;
};
double b_coefficients[6]=
{0.13654762463195771
,0.42313861743656667
,0.6775400499741616
,0.839889624849638
,0.9315419599631839
,0.9878163707328971
};
filter_a=new CAllPassFilterCascade(a_coefficients,6);
filter_b=new CAllPassFilterCascade(b_coefficients,6);
}
else if (order==10) //rejection=86dB, transition band=0.01
{
double a_coefficients[5]=
{0.051457617441190984
,0.35978656070567017
,0.6725475931034693
,0.8590884928249939
,0.9540209867860787
};
double b_coefficients[5]=
{0.18621906251989334
,0.529951372847964
,0.7810257527489514
,0.9141815687605308
,0.985475023014907
};
filter_a=new CAllPassFilterCascade(a_coefficients,5);
filter_b=new CAllPassFilterCascade(b_coefficients,5);
}
else if (order==8) //rejection=69dB, transition band=0.01
{
double a_coefficients[4]=
{0.07711507983241622
,0.4820706250610472
,0.7968204713315797
,0.9412514277740471
};
double b_coefficients[4]=
{0.2659685265210946
,0.6651041532634957
,0.8841015085506159
,0.9820054141886075
};
filter_a=new CAllPassFilterCascade(a_coefficients,4);
filter_b=new CAllPassFilterCascade(b_coefficients,4);
}
else if (order==6) //rejection=51dB, transition band=0.01
{
double a_coefficients[3]=
{0.1271414136264853
,0.6528245886369117
,0.9176942834328115
};
double b_coefficients[3]=
{0.40056789819445626
,0.8204163891923343
,0.9763114515836773
};
filter_a=new CAllPassFilterCascade(a_coefficients,3);
filter_b=new CAllPassFilterCascade(b_coefficients,3);
}
else if (order==4) //rejection=53dB,transition band=0.05
{
double a_coefficients[2]=
{0.12073211751675449
,0.6632020224193995
};
double b_coefficients[2]=
{0.3903621872345006
,0.890786832653497
};
filter_a=new CAllPassFilterCascade(a_coefficients,2);
filter_b=new CAllPassFilterCascade(b_coefficients,2);
}
filter_a=new CAllPassFilterCascade(&a_coefficients,1);
filter_b=new CAllPassFilterCascade(&b_coefficients,1);
}
}
else //softer slopes, more attenuation and less stopband ripple
{
if (order==12) //rejection=150dB, transition band=0.05
{
double a_coefficients[6]=
{0.01677466677723562
,0.13902148819717805
,0.3325011117394731
,0.53766105314488
,0.7214184024215805
,0.8821858402078155
};
double b_coefficients[6]=
{0.06501319274445962
,0.23094129990840923
,0.4364942348420355
,0.06329609551399348
,0.80378086794111226
,0.9599687404800694
};
filter_a=new CAllPassFilterCascade(a_coefficients,6);
filter_b=new CAllPassFilterCascade(b_coefficients,6);
}
else if (order==10) //rejection=133dB, transition band=0.05
{
double a_coefficients[5]=
{0.02366831419883467
,0.18989476227180174
,0.43157318062118555
,0.6632020224193995
,0.860015542499582
};
double b_coefficients[5]=
{0.09056555904993387
,0.3078575723749043
,0.5516782402507934
,0.7652146863779808
,0.95247728378667541
};
filter_a=new CAllPassFilterCascade(a_coefficients,5);
filter_b=new CAllPassFilterCascade(b_coefficients,5);
}
else if (order==8) //rejection=106dB, transition band=0.05
{
double a_coefficients[4]=
{0.03583278843106211
,0.2720401433964576
,0.5720571972357003
,0.827124761997324
};
double b_coefficients[4]=
{0.1340901419430669
,0.4243248712718685
,0.7062921421386394
,0.9415030941737551
};
filter_a=new CAllPassFilterCascade(a_coefficients,4);
filter_b=new CAllPassFilterCascade(b_coefficients,4);
}
else if (order==6) //rejection=80dB, transition band=0.05
{
double a_coefficients[3]=
{0.06029739095712437
,0.4125907203610563
,0.7727156537429234
};
double b_coefficients[3]=
{0.21597144456092948
,0.6043586264658363
,0.9238861386532906
};
filter_a=new CAllPassFilterCascade(a_coefficients,3);
filter_b=new CAllPassFilterCascade(b_coefficients,3);
}
double b_coefficients[2]=
{0.28382934487410993
,0.8344118914807379
};
filter_a=new CAllPassFilterCascade(a_coefficients,2);
filter_b=new CAllPassFilterCascade(b_coefficients,2);
}
filter_a=new CAllPassFilterCascade(&a_coefficients,1);
filter_b=new CAllPassFilterCascade(&b_coefficients,1);
}
}
oldout=0.0;
};
CHalfBandFilter::~CHalfBandFilter()
{
delete filter_a;
delete filter_b;
};
return output;
}
class CAllPassFilter
{
public:
private:
double a;
double x0;
double x1;
double x2;
double y0;
double y1;
double y2;
};
class CAllPassFilterCascade
{
public:
CAllPassFilterCascade(const double* coefficients, int N);
~CAllPassFilterCascade();
private:
CAllPassFilter** allpassfilter;
int numfilters;
};
class CHalfBandFilter
{
public:
CHalfBandFilter(const int order, const bool steep);
~CHalfBandFilter();
private:
CAllPassFilterCascade* filter_a;
CAllPassFilterCascade* filter_b;
double oldout;
};
pow(x,4) approximation
Notes :
Very hacked, but it gives a rough estimate of x**4 by modifying exponent and mantissa.
Code :
float p4fast(float in)
{
long *lp,l;
lp=(long *)(&in);
l=*lp;
l-=0x3F800000l; /* un-bias */
l<<=2; /* **4 */
l+=0x3F800000l; /* bias */
*lp=l;
/* compiler will read this from memory since & operator had been used */
return in;
}
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):
= j * (2/T) * tan(w*T/2)
or
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.
Pseudo-Random generator
Notes :
This can be used to generate random numeric sequences or to synthesise a white noise audio
signal.
If you only use some of the bits, use the most significant bits by shifting right.
Do not just mask off the low bits.
Code :
/* Calculate pseudo-random 32 bit number based on linear congruential
method. */
unsigned long GenerateRandomNumber( void )
{
/* Change this for different random sequences. */
static unsigned long randSeed = 22222;
randSeed = (randSeed * 196314165) + 907633515;
return randSeed;
}
Pulsewidth modulation
Notes :
Take an upramping sawtooth and its inverse, a downramping sawtooth. Adding these two
waves with a well defined delay between 0 and period (1/f)
results in a square wave with a duty cycle ranging from 0 to 100%.
RBJ-Audio-EQ-Cookbook
Notes :
(see linkfile)
A superb collection of filters used in a lot of (commercial) plugins and effects.
There's also a very interesting paper linked about biquad EQ filters.
Notes :
(see linkfile)
Code to read the .WA! (compressed .WAV) parts of GigaSampler .GIG files.
For related info on reading .GIG files see http://www.linuxdj.com/evo
Resonant filter
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);
Notes :
Hard to calculate coefficients, easy to process algorithm
Code :
resofreq = pole frequency
amp = magnitude at pole frequency (approx)
double pi = 3.141592654;
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++) {
/* Check clipping */
temp = vibrapos;
if (temp > 32767) {
temp = 32767;
} else if (temp < -32768) temp = -32768;
Linked files
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/**************************************************************************
*************************************************************************/
typedef struct {
unsigned int length; /* size of filter */
float *history; /* pointer to history in filter */
float *coef; /* pointer to coefficients of filter */
} FILTER;
typedef struct {
double a0, a1, a2; /* numerator coefficients */
double b0, b1, b2; /* denominator coefficients */
} BIQUAD;
void szxform(
double *a0, double *a1, double *a2, /* numerator coefficients */
double *b0, double *b1, double *b2, /* denominator coefficients */
double fc, /* Filter cutoff frequency */
double fs, /* sampling rate */
double *k, /* overall gain factor */
/*
* --------------------------------------------------------------------
*
* iir_filter - Perform IIR filtering sample by sample on floats
*
* Implements cascaded direct form II second order sections.
* Requires FILTER structure for history and coefficients.
* The length in the filter structure specifies the number of sections.
* The size of the history array is 2*iir->length.
* The size of the coefficient array is 4*iir->length + 1 because
* the first coefficient is the overall scale factor for the filter.
* Returns one output sample for each input sample. Allocates history
* array if not previously allocated.
*
* float iir_filter(float input,FILTER *iir)
*
* float input new float input sample
* FILTER *iir pointer to FILTER structure
*
* Returns float value giving the current output.
*
* Allocation errors cause an error message and a call to exit.
* --------------------------------------------------------------------
*/
float iir_filter(input,iir)
float input; /* new input sample */
FILTER *iir; /* pointer to FILTER structure */
{
unsigned int i;
float *hist1_ptr,*hist2_ptr,*coef_ptr;
float output,new_hist,history1,history2;
if(!iir->history) {
iir->history = (float *) calloc(2*iir->length,sizeof(float));
if(!iir->history) {
printf("\nUnable to allocate history array in iir_filter\n");
exit(1);
}
}
*hist2_ptr++ = *hist1_ptr;
*hist1_ptr++ = new_hist;
hist1_ptr++;
hist2_ptr++;
}
return(output);
}
/*
* --------------------------------------------------------------------
*
* main()
*
* Example main function to show how to update filter coefficients.
* We create a 4th order filter (24 db/oct roloff), consisting
* of two second order sections.
* --------------------------------------------------------------------
*/
int main()
{
FILTER iir;
float *coef;
double fs, fc; /* Sampling frequency, cutoff frequency */
double Q; /* Resonance > 1.0 < 1000 */
unsigned nInd;
double a0, a1, a2, b0, b1, b2;
double k; /* overall gain factor */
/*
* Setup filter s-domain coefficients
*/
/* Section 1 */
ProtoCoef[0].a0 = 1.0;
ProtoCoef[0].a1 = 0;
ProtoCoef[0].a2 = 0;
ProtoCoef[0].b0 = 1.0;
ProtoCoef[0].b1 = 0.765367;
ProtoCoef[0].b2 = 1.0;
/* Section 2 */
ProtoCoef[1].a0 = 1.0;
ProtoCoef[1].a1 = 0;
ProtoCoef[1].a2 = 0;
ProtoCoef[1].b0 = 1.0;
ProtoCoef[1].b1 = 1.847759;
ProtoCoef[1].b2 = 1.0;
/*
* Allocate array of z-domain coefficients for each filter section
* plus filter gain variable
*/
iir.coef = (float *) calloc(4 * iir.length + 1, sizeof(float));
if (!iir.coef)
{
printf("Unable to allocate coef array, exiting\n");
exit(1);
}
Q = 1; /* Resonance */
fc = 5000; /* Filter cutoff (Hz) */
fs = 44100; /* Sampling frequency (Hz) */
/*
* Compute z-domain coefficients for each biquad section
* for new Cutoff Frequency and Resonance
*/
for (nInd = 0; nInd < iir.length; nInd++)
{
a0 = ProtoCoef[nInd].a0;
a1 = ProtoCoef[nInd].a1;
a2 = ProtoCoef[nInd].a2;
b0 = ProtoCoef[nInd].b0;
b1 = ProtoCoef[nInd].b1 / Q; /* Divide by resonance or Q
*/
b2 = ProtoCoef[nInd].b2;
szxform(&a0, &a1, &a2, &b0, &b1, &b2, fc, fs, &k, coef);
coef += 4; /* Point to next filter
section */
}
iir.coef[0] = k;
Reposting bilinear.c just in case the other one was not the latest version.
#include <math.h>
void prewarp(double *a0, double *a1, double *a2, double fc, double fs);
void bilinear(
double a0, double a1, double a2, /* numerator coefficients */
double b0, double b1, double b2, /* denominator coefficients */
double *k, /* overall gain factor */
double fs, /* sampling rate */
float *coef); /* pointer to 4 iir coefficients */
/*
* ----------------------------------------------------------
* Pre-warp the coefficients of a numerator or denominator.
* Note that a0 is assumed to be 1, so there is no wrapping
* of it.
* ----------------------------------------------------------
*/
void prewarp(
double *a0, double *a1, double *a2,
double fc, double fs)
{
double wp, pi;
pi = 4.0 * atan(1.0);
wp = 2.0 * fs * tan(pi * fc / fs);
/*
* ----------------------------------------------------------
* bilinear()
*
* Transform the numerator and denominator coefficients
* of s-domain biquad section into corresponding
* z-domain coefficients.
*
* Store the 4 IIR coefficients in array pointed by coef
* in following order:
* beta1, beta2 (denominator)
* alpha1, alpha2 (numerator)
*
* Arguments:
* a0-a2 - s-domain numerator coefficients
* b0-b2 - s-domain denominator coefficients
* k - filter gain factor. initially set to 1
* and modified by each biquad section in such
* a way, as to make it the coefficient by
* which to multiply the overall filter gain
* in order to achieve a desired overall filter
gain,
* specified in initial value of k.
* fs - sampling rate (Hz)
* coef - array of z-domain coefficients to be filled in.
*
* Return:
* On return, set coef z-domain coefficients
* ----------------------------------------------------------
*/
void bilinear(
double a0, double a1, double a2, /* numerator coefficients */
double b0, double b1, double b2, /* denominator coefficients */
double *k, /* overall gain factor */
double fs, /* sampling rate */
float *coef /* pointer to 4 iir coefficients */
)
{
double ad, bd;
/* Denominator */
*coef++ = (2. * b0 - 8. * b2 * fs * fs)
/ bd; /* beta1 */
*coef++ = (4. * b2 * fs * fs - 2. * b1 * fs + b0)
/ bd; /* beta2 */
/* Nominator */
*coef++ = (2. * a0 - 8. * a2 * fs * fs)
/ ad; /* alpha1 */
*coef = (4. * a2 * fs * fs - 2. * a1 * fs + a0)
/ ad; /* alpha2 */
}
/*
* ----------------------------------------------------------
* Transform from s to z domain using bilinear transform
* with prewarp.
*
* Arguments:
* For argument description look at bilinear()
*
* coef - pointer to array of floating point coefficients,
* corresponding to output of bilinear transofrm
* (z domain).
*
* Note: frequencies are in Hz.
* ----------------------------------------------------------
*/
void szxform(
double *a0, double *a1, double *a2, /* numerator coefficients */
The filter consists of two 2nd order secions since highest s power is 2.
Now we can take the coefficients, or the numbers by which s is multiplied
and plug them into a standard formula to be used by bilinear transform.
a2 * s^2 + a1 * s + a0
H(s) = ----------------------
b2 * s^2 + b1 * s + b0
0 + 0 + 1 0 + 0 + 1
-- ----------------------- * --------------------------
1 + ((1/Q) * 0.765367) + 1 1 + ((1/Q) * 1.847759) + 1
Section 1:
a2 = 0; a1 = 0; a0 = 1;
b2 = 1; b1 = 0.5176387; b0 = 1;
Section 2:
a2 = 0; a1 = 0; a0 = 1;
b2 = 1; b1 = 1.847759; b0 = 1;
Upon return from each call, the k argument will be set to a value,
by which to multiply our actual signal in order for the gain
to be one. On second call to szxform() we provide k that was
changed by the previous section. During actual audio filtering
function iir_filter() will use this k
Summary:
*
* ----------------------------------------------------------
*References:
*Van Valkenburg, "Analog Filter Design"
*Oxford University Press 1982
*ISBN 0-19-510734-9
*
*C Language Algorithms for Digital Signal Processing
*Paul Embree, Bruce Kimble
*Prentice Hall, 1991
*ISBN 0-13-133406-9
*
*Digital Filter Designer's Handbook
*With C++ Algorithms
*Britton Rorabaugh
*McGraw Hill, 1997
*ISBN 0-07-053806-9
* ----------------------------------------------------------
Reverberation techniques
Notes :
* Parallel comb filters, followed by series allpass filters. This was the original design by
Schroeder, and was extended by Moorer. Has a VERY metallic sound for sharp transients.
* Several allpass filters in serie (also proposed by Schroeder). Also suffers from metallic sound.
* 2nd-order comb and allpass filters (described by Moorer). Not supposed to give much of an
advantage over first order sections.
* Nested allpass filters, where an allpass filter will replace the delay line in another allpass
filter. Pioneered by Gardner. Haven't heard the results.
* Strange allpass amp delay line based structure in Jon Dattorro article (JAES). Four allpass
filters are used as an input to a cool "figure-8" feedback loop, where four allpass reverberators
are used in series with
a few delay lines. Outputs derived from various taps in structure. Supposedly based on a
Lexicon reverb design. Modulating delay lines are used in some of the allpass structures to
"spread out" the eigentones.
* Waveguide-based reverbs, where the reverb structure is based upon the junction of many
waveguides. Julius Smith developed these. Recently, these have been shown to be essentially
equivalent to the feedback delay network reverbs. Also sound very nice. Modulating delay lines
and scattering values used to spread out eigentones.
* Convolution-based reverbs, where the sound to be reverbed is convolved with the impulse
response of a room, or with exponentially-decaying white noise. Supposedly the best sound, but
very computationally expensive, and not very flexible.
* FIR-based reverbs. Essentially the same as convolution. Probably not used, but shorter FIR
filters are probably used in combination with many of the above techniques, to provide early
reflections.
Saturation
Type : Waveshaper
References : Posted by Bram
Notes :
when the input is below a certain threshold (t) these functions return the input, if it goes over
that threshold, they return a soft shaped saturation.
Neigther claims to be fast ;-)
Code :
float saturate(float x, float t)
{
if(fabs(x)<t)
return x
else
{
if(x > 0.f);
return t + (1.f-t)*tanh((x-t)/(1-t));
else
return -(t + (1.f-t)*tanh((-x-t)/(1-t)));
}
}
or
float sigmoid(x)
{
if(fabs(x)<1)
return x*(1.5f - 0.5f*x*x);
else
return x > 0.f ? 1.f : -1.f;
}
SawSin
Notes :
This simple peak follower will give track the peaks of a signal. It will rise rapidly when the
input is rising, and then decay exponentially when the input drops. It can be used to drive VU
meters, or used in an automatic gain control circuit.
Code :
// halfLife = time in seconds for output to decay to half value after an
impulse
Notes :
Code for approximation of cos, sin, tan and inv sin, etc.
Surprisingly accurate and very usable.
Code :
Linked files
//----------------------------------------------------------------------
Real Math::FastSin0 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = 7.61e-03f;
fResult *= fASqr;
fResult -= 1.6605e-01f;
fResult *= fASqr;
fResult += 1.0f;
fResult *= fAngle;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastSin1 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = -2.39e-08f;
fResult *= fASqr;
fResult += 2.7526e-06f;
fResult *= fASqr;
fResult -= 1.98409e-04f;
fResult *= fASqr;
fResult += 8.3333315e-03f;
fResult *= fASqr;
fResult -= 1.666666664e-01f;
fResult *= fASqr;
fResult += 1.0f;
fResult *= fAngle;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastCos0 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = 3.705e-02f;
fResult *= fASqr;
fResult -= 4.967e-01f;
fResult *= fASqr;
fResult += 1.0f;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastCos1 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = -2.605e-07f;
fResult *= fASqr;
fResult += 2.47609e-05f;
fResult *= fASqr;
fResult -= 1.3888397e-03f;
fResult *= fASqr;
fResult += 4.16666418e-02f;
fResult *= fASqr;
fResult -= 4.999999963e-01f;
fResult *= fASqr;
fResult += 1.0f;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastTan0 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = 2.033e-01f;
fResult *= fASqr;
fResult += 3.1755e-01f;
fResult *= fASqr;
fResult += 1.0f;
fResult *= fAngle;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastTan1 (Real fAngle)
{
Real fASqr = fAngle*fAngle;
Real fResult = 9.5168091e-03f;
fResult *= fASqr;
fResult += 2.900525e-03f;
fResult *= fASqr;
fResult += 2.45650893e-02f;
fResult *= fASqr;
fResult += 5.33740603e-02f;
fResult *= fASqr;
fResult += 1.333923995e-01f;
fResult *= fASqr;
fResult += 3.333314036e-01f;
fResult *= fASqr;
fResult += 1.0f;
fResult *= fAngle;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastInvSin (Real fValue)
{
Real fRoot = Math::Sqrt(1.0f-fValue);
Real fResult = -0.0187293f;
fResult *= fValue;
fResult += 0.0742610f;
fResult *= fValue;
fResult -= 0.2121144f;
fResult *= fValue;
fResult += 1.5707288f;
fResult = HALF_PI - fRoot*fResult;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastInvCos (Real fValue)
{
Real fRoot = Math::Sqrt(1.0f-fValue);
Real fResult = -0.0187293f;
fResult *= fValue;
fResult += 0.0742610f;
fResult *= fValue;
fResult -= 0.2121144f;
fResult *= fValue;
fResult += 1.5707288f;
fResult *= fRoot;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastInvTan0 (Real fValue)
{
Real fVSqr = fValue*fValue;
Real fResult = 0.0208351f;
fResult *= fVSqr;
fResult -= 0.085133f;
fResult *= fVSqr;
fResult += 0.180141f;
fResult *= fVSqr;
fResult -= 0.3302995f;
fResult *= fVSqr;
fResult += 0.999866f;
fResult *= fValue;
return fResult;
}
//----------------------------------------------------------------------
Real Math::FastInvTan1 (Real fValue)
{
Real fVSqr = fValue*fValue;
Real fResult = 0.0028662257f;
fResult *= fVSqr;
fResult -= 0.0161657367f;
fResult *= fVSqr;
fResult += 0.0429096138f;
fResult *= fVSqr;
fResult -= 0.0752896400f;
fResult *= fVSqr;
fResult += 0.1065626393f;
fResult *= fVSqr;
fResult -= 0.1420889944f;
fResult *= fVSqr;
fResult += 0.1999355085f;
fResult *= fVSqr;
fResult -= 0.3333314528f;
fResult *= fVSqr;
fResult += 1.0f;
fResult *= fValue;
return fResult;
}
//----------------------------------------------------------------------
Sine calculation
Notes :
Code from JSyn for a sine wave generator based on a Taylor Expansion. It is not as efficient as
the filter methods, but it has linear frequency control and is, therefore, suitable for FM or other
time varying applications where accurate frequency is needed. The sine generated is accurate to
at least 16 bits.
Code :
for(i=0; i < nSamples ; i++)
{
//Generate sawtooth phasor to provide phase for sine generation
IncrementWrapPhase(phase, freqPtr[i]);
//Wrap phase back into region where results are more accurate
x = yp * PI;
x2 = x*x;
Type : Pitch Scaling (often incorrectly referred to as "Pitch Shifting") using the Fourier
transform
References : Posted by sms[AT]dspdimension.com
Linked file : http://www.dspdimension.com
Code :
See above web site
Soft saturation
Type : waveshaper
References : Posted by Bram de Jong
Notes :
This only works for positive values of x. a should be in the range 0..1
Code :
x < a:
f(x) = x
x > a:
f(x) = a + (x-a)/(1+((x-a)/(1-a))^2)
x > 1:
f(x) = (a+1)/2
Square Waves
Notes :
One way to do a square wave:
You need two buzz generators (see Dodge & Jerse, or the Csound source code, for
implementation details). One of the buzz generators runs at the desired square wave frequency,
while the second buzz generator is exactly one octave above this pitch. Subtract the higher
octave buzz generator's output from the lower buzz generator's output - the result should be a
signal with all odd harmonics, all at equal amplitude. Filter the resultant signal (maybe integrate
it). Voila, a bandlimited square wave! Well, I think it should work...
The one question I have with the above technique is whether it produces a waveform that truly
resembles a square wave in the time domain. Even if the number of harmonics, and the relative
ratio of the harmonics, is identical to an "ideal" bandwidth-limited square wave, it may have an
entirely different waveshape. No big deal, unless the signal is processed by a nonlinearity, in
which case the results of the nonlinear processing will be far different than the processing of a
waveform that has a similar shape to a square wave.
State variable
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
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 = MIN(0.25, 2.0*sin(PI*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;
}
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.
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.
This we knew already: We can split a convolution into four convolutions of halved size.
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:
= (A1*B2 + A2*B1)*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).
= 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 :)
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;
if (size == 1)
{
out[0] = in1[0] * in2[0];
return;
}
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];
}
}
}
}
}
Type : Goertzel
References : Posted by [email protected]
Linked file : http://www.ii.uib.no/~espenr/tonedetect.zip
Notes :
Goertzel is basically DFT of parts of a spectrum not the total spectrum as you normally do with
FFT. So if you just want to check out the power for some frequencies this could be better. Is
good for DTFM detection I've heard.
The WNk isn't calculated 100% correctly, but it seems to work so ;) Yeah and the code is C++
so you might have to do some small adjustment to compile it as C.
Code :
/** Tone detect by Goertzel algorithm
*
* This program basically searches for tones (sines) in a sample and reports
the different dB it finds for
* different frequencies. Can easily be extended with some thresholding to
report true/false on detection.
* I'm far from certain goertzel it implemented 100% correct, but it works
:)
*
* Hint, the SAMPLERATE, BUFFERSIZE, FREQUENCY, NOISE and SIGNALVOLUME all
affects the outcome of the reported dB. Tweak
* em to find the settings best for your application. Also, seems to be
pretty sensitive to noise (whitenoise anyway) which
* is a bit sad. Also I don't know if the goertzel really likes float values
for the frequency ... And using 44100 as
* samplerate for detecting 6000 Hz tone is kinda silly I know :)
*
* Written by: Espen Riskedal, [email protected], july-2002
*/
#include <iostream>
#include <cmath>
#include <cstdlib>
using std::rand;
// math stuff
using std::cos;
using std::abs;
using std::exp;
using std::log10;
// iostream stuff
using std::cout;
using std::endl;
#define PI 3.14159265358979323844
// change the defines if you want to
#define SAMPLERATE 44100
#define BUFFERSIZE 8820
#define FREQUENCY 6000
#define NOISE 0.05
#define SIGNALVOLUME 0.8
/** The Goertzel algorithm computes the k-th DFT coefficient of the input
signal using a second-order filter.
* http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html.
* Basiclly it just does a DFT of the frequency we want to check, and none
of the others (FFT calculates for all frequencies).
*/
float goertzel(float *x, int N, float frequency, int samplerate) {
float Skn, Skn1, Skn2;
Skn = Skn1 = Skn2 = 0;
return tone;
}
return 0;
}
Notes :
k >= 1 is the "clipping hardness". 1 gives a smooth clipping, and a high value gives
hardclipping.
Don't set k too high, because the formula use the pow() function, which use exp() and would
overflow easily. 100 seems to be a reasonable value for "hardclipping"
Code :
f (x) = sign (x) * pow (atan (pow (abs (x), k)), (1 / k));
Linked files
Notes :
(see linkfile)
Filters included are:
presence
shelvelowpass
2polebp
peaknotch
peaknotch2
Linked files
/*
* Presence and Shelve filters as given in
* James A. Moorer
* The manifold joys of conformal mapping:
* applications to digital filtering in the studio
* JAES, Vol. 31, No. 11, 1983 November
*/
double bw2angle(a,bw)
double a,bw;
{
double T,d,sn,cs,mag,delta,theta,tmp,a2,a4,asnd;
T = tan(2.0*PI*bw);
a2 = a*a;
a4 = a2*a2;
d = 2.0*a2*T;
sn = (1.0 + a4)*T;
cs = (1.0 - a4);
mag = sqrt(sn*sn + cs*cs);
d /= mag;
delta = atan2(sn,cs);
asnd = asin(d);
theta = 0.5*(PI - asnd - delta);
tmp = 0.5*(asnd-delta);
if ((tmp > 0.0) && (tmp < theta)) theta = tmp;
return(theta/(2.0*PI));
}
void presence(cf,boost,bw,a0,a1,a2,b1,b2)
double cf,boost,bw,*a0,*a1,*a2,*b1,*b2;
{
double a,A,F,xfmbw,C,tmp,alphan,alphad,b0,recipb0,asq,F2,a2plus1,ma2plus1;
a = tan(PI*(cf-0.25));
asq = a*a;
A = pow(10.0,boost/20.0);
if ((boost < 6.0) && (boost > -6.0)) F = sqrt(A);
else if (A > 1.0) F = A/sqrt(2.0);
else F = A*sqrt(2.0);
xfmbw = bw2angle(a,bw);
C = 1.0/tan(2.0*PI*xfmbw);
F2 = F*F;
tmp = A*A - F2;
if (fabs(tmp) <= SPN) alphad = C;
else alphad = sqrt(C*C*(F2-1.0)/tmp);
alphan = A*alphad;
b0 = a2plus1 + alphad*ma2plus1;
*b2 = a2plus1 - alphad*ma2plus1;
recipb0 = 1.0/b0;
*a0 *= recipb0;
*a1 *= recipb0;
*a2 *= recipb0;
*b1 = *a1;
*b2 *= recipb0;
}
void shelve(cf,boost,a0,a1,a2,b1,b2)
double cf,boost,*a0,*a1,*a2,*b1,*b2;
{
double a,A,F,tmp,b0,recipb0,asq,F2,gamma2,siggam2,gam2p1;
double gamman,gammad,ta0,ta1,ta2,tb0,tb1,tb2,aa1,ab1;
a = tan(PI*(cf-0.25));
asq = a*a;
A = pow(10.0,boost/20.0);
if ((boost < 6.0) && (boost > -6.0)) F = sqrt(A);
else if (A > 1.0) F = A/sqrt(2.0);
else F = A*sqrt(2.0);
F2 = F*F;
gamma2 = gamman*gamman;
gam2p1 = 1.0 + gamma2;
siggam2 = 2.0*sqrt(2.0)/2.0*gamman;
ta0 = gam2p1 + siggam2;
ta1 = -2.0*(1.0 - gamma2);
ta2 = gam2p1 - siggam2;
gamma2 = gammad*gammad;
gam2p1 = 1.0 + gamma2;
siggam2 = 2.0*sqrt(2.0)/2.0*gammad;
tb0 = gam2p1 + siggam2;
tb1 = -2.0*(1.0 - gamma2);
tb2 = gam2p1 - siggam2;
aa1 = a*ta1;
*a0 = ta0 + aa1 + asq*ta2;
*a1 = 2.0*a*(ta0+ta2)+(1.0+asq)*ta1;
*a2 = asq*ta0 + aa1 + ta2;
ab1 = a*tb1;
b0 = tb0 + ab1 + asq*tb2;
*b1 = 2.0*a*(tb0+tb2)+(1.0+asq)*tb1;
*b2 = asq*tb0 + ab1 + tb2;
recipb0 = 1.0/b0;
*a0 *= recipb0;
*a1 *= recipb0;
*a2 *= recipb0;
*b1 *= recipb0;
*b2 *= recipb0;
}
void initfilter(f)
filter *f;
{
f->x1 = 0.0;
f->x2 = 0.0;
f->y1 = 0.0;
f->y2 = 0.0;
f->y = 0.0;
}
void setfilter_presence(f,freq,boost,bw)
filter *f;
double freq,boost,bw;
{
presence(freq/(double)SR,boost,bw/(double)SR,
&f->cx,&f->cx1,&f->cx2,&f->cy1,&f->cy2);
f->cy1 = -f->cy1;
f->cy2 = -f->cy2;
}
void setfilter_shelve(f,freq,boost)
filter *f;
double freq,boost;
{
shelve(freq/(double)SR,boost,
&f->cx,&f->cx1,&f->cx2,&f->cy1,&f->cy2);
f->cy1 = -f->cy1;
f->cy2 = -f->cy2;
}
void setfilter_shelvelowpass(f,freq,boost)
filter *f;
double freq,boost;
{
double gain;
gain = pow(10.0,boost/20.0);
shelve(freq/(double)SR,boost,
&f->cx,&f->cx1,&f->cx2,&f->cy1,&f->cy2);
f->cx /= gain;
f->cx1 /= gain;
f->cx2 /= gain;
f->cy1 = -f->cy1;
f->cy2 = -f->cy2;
}
/*
* As in ''An introduction to digital filter theory'' by Julius O. Smith
* and in Moore's book; I use the normalized version in Moore's book.
*/
void setfilter_2polebp(f,freq,R)
filter *f;
double freq,R;
{
double theta;
theta = 2.0*PI*freq/(double)SR;
f->cx = 1.0-R;
f->cx1 = 0.0;
f->cx2 = -(1.0-R)*R;
f->cy1 = 2.0*R*cos(theta);
f->cy2 = -R*R;
}
/*
* As in
* Stanley A. White
* Design of a digital biquadratic peaking or notch filter
* for digital audio equalization
* JAES, Vol. 34, No. 6, 1986 June
*/
void setfilter_peaknotch(f,freq,M,bw)
filter *f;
double freq,M,bw;
{
double w0,p,om,ta,d;
w0 = 2.0*PI*freq;
if ((1.0/sqrt(2.0) < M) && (M < sqrt(2.0))) {
fprintf(stderr,"peaknotch filter: 1/sqrt(2) < M < sqrt(2)\n");
exit(-1);
}
if (M <= 1.0/sqrt(2.0)) p = sqrt(1.0-2.0*M*M);
if (sqrt(2.0) <= M) p = sqrt(M*M-2.0);
om = 2.0*PI*bw;
ta = tan(om/((double)SR*2.0));
d = p+ta;
f->cx = (p+M*ta)/d;
f->cx1 = -2.0*p*cos(w0/(double)SR)/d;
f->cx2 = (p-M*ta)/d;
f->cy1 = 2.0*p*cos(w0/(double)SR)/d;
f->cy2 = -(p-ta)/d;
}
/*
* Some JAES's article on ladder filter.
* freq (Hz), gdb (dB), bw (Hz)
*/
void setfilter_peaknotch2(f,freq,gdb,bw)
filter *f;
double freq,gdb,bw;
{
double k,w,bwr,abw,gain;
k = pow(10.0,gdb/20.0);
w = 2.0*PI*freq/(double)SR;
bwr = 2.0*PI*bw/(double)SR;
abw = (1.0-tan(bwr/2.0))/(1.0+tan(bwr/2.0));
gain = 0.5*(1.0+k+abw-k*abw);
f->cx = 1.0*gain;
f->cx1 = gain*(-2.0*cos(w)*(1.0+abw))/(1.0+k+abw-k*abw);
f->cx2 = gain*(abw+k*abw+1.0-k)/(abw-k*abw+1.0+k);
f->cy1 = 2.0*cos(w)/(1.0+tan(bwr/2.0));
f->cy2 = -abw;
}
double applyfilter(f,x)
filter *f;
double x;
{
f->x = x;
f->y = f->cx * f->x + f->cx1 * f->x1 + f->cx2 * f->x2
+ f->cy1 * f->y1 + f->cy2 * f->y2;
f->x2 = f->x1;
f->x1 = f->x;
f->y2 = f->y1;
f->y1 = f->y;
return(f->y);
}
Notes :
C code and project file for MSVC6 for a bandwidth-limited saw/square (with PWM) generator
using MinBLEPS.
This code is based on Eli's MATLAB MinBLEP code and uses his original minblep.mat file.
Instead of keeping a list of all active MinBLEPS, the output of each MinBLEP is stored in a
buffer, in which all consequent MinBLEPS and the waveform output are added together. This
optimization makes it fast enough to be used realtime.
Produces slight aliasing when sweeping high frequencies. I don't know wether Eli's original
code does the same, because I don't have MATLAB. Any help would be appreciated.
The project name is 'hardsync', because it's easy to generate hardsync using MinBLEPS.
Code :
WaveShaper
Type : waveshaper
References : Posted by Bram de Jong
Notes :
where x (in [-1..1] will be distorted and a is a distortion parameter that goes from 1 to infinity
The equation is valid for positive and negativ values.
If a is 1, it results in a slight distortion and with bigger a's the signal get's more funky.
Code :
f(x,a) = x*(abs(x) + a)/(x^2 + (a-1)*abs(x) + 1)
Waveshaper
Type : waveshaper
References : Posted by Jon Watte
Notes :
A favourite of mine is using a sin() function instead.
This will have the "unfortunate" side effect of removing
odd harmonics if you take it to the extreme: a triangle
wave gets mapped to a pure sine wave.
This will work with a going from .1 or so to a= 5 and bigger!
The mathematical limits for a = 0 actually turns it into a linear
function at that point, but unfortunately FPUs aren't that good
with calculus :-) Once a goes above 1, you start getting clipping
in addition to the "soft" wave shaping. It starts getting into
more of an effect and less of a mastering tool, though :-)
Code :
(input: a == "overdrive amount")
z = M_PI * a;
s = 1/sin(z)
b = 1/a
if (x > b)
f(x) = 1
else
f(x) = sin(z*x)*s
Waveshaper
Notes :
amount should be in [-1..1[ Plot it and stand back in astonishment! ;)
Code :
x = input in [-1..1]
y = output
k = 2*amount/(1-amount);
f(x) = (1+k)*x/(1+k*abs(x))
Notes :
> The other question; what's a 'waveshaper' algorithm. Is it simply another
> word for distortion?
Code :
float waveshape_distort( float in ) {
return 1.5f * in - 0.5f * in *in * in;
}
Waveshaper :: Gloubi-boulga
Notes :
Multiply input by gain before processing
Code :
const double x = input * 0.686306;
const double a = 1 + exp (sqrt (fabs (x)) * -0.75);
output = (exp (x) - exp (-x * a)) / (exp (x) + exp (-x));
Wavetable Synthesis
Notes :
Wavetable sythesis AES paper by RBJ.
Weird synthesis
Notes :
(quoted from Andy's mail...)
What I've done in a soft-synth I've been working on is used what I've termed Fooglers, no
reason, just liked the name :) Anyway all I've done is use a *VERY* short delay line of 256
samples and then use 2 controllable taps into the delay with High Frequency Damping, and a
feedback parameter.
Using a tiny fixed delay size of approx. 4.8ms (really 256 samples/1k memory with floats)
means this costs, in terms of cpu consumption practically nothing, and the filter is a real simple
1 pole low-pass filter. Maybe not DSP'litically correct but all I wanted was to avoid the high
frequencies trashing the delay line when high feedbacks (99%->99.9%) are used (when the fun
starts ;).
I've been getting some really sexy sounds out of this idea, and of course you can have the delay
line tuneable if you choose to use fractional taps, but I'm happy with it as it is.. 1 nice simple,
yet powerful addition to the base oscillators.
In reality you don't need 2 taps, but I found that using 2 added that extra element of funkiness...
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