# Musicdsp.org¶

Welcome to musicdsp.org.

Musicdsp.org is a collection of algorithms, thoughts and snippets, gathered for the music dsp community. Most of this data was gathered by and for the people of the splendid Music-DSP mailing list at http://sites.music.columbia.edu/cmc/music-dsp/

Important

## Synthesis¶

• Author or source: Ross Bencina, Olli Niemitalo, …
• Type: waveform generation
• Created: 2002-01-17 01:01:39
notes
```Ross Bencina: original source code poster
Olli Niemitalo: UpdateWithCubicInterpolation
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84``` ```//this code is meant as an EXAMPLE //uncomment if you need an FM oscillator //define FM_OSCILLATOR /* members are: float phase; int TableSize; float sampleRate; float *table, dtable0, dtable1, dtable2, dtable3; ->these should be filled as folows... (remember to wrap around!!!) table[i] = the wave-shape dtable0[i] = table[i+1] - table[i]; dtable1[i] = (3.f*(table[i]-table[i+1])-table[i-1]+table[i+2])/2.f dtable2[i] = 2.f*table[i+1]+table[i-1]-(5.f*table[i]+table[i+2])/2.f dtable3[i] = (table[i+1]-table[i-1])/2.f */ float Oscillator::UpdateWithoutInterpolation(float frequency) { int i = (int) phase; phase += (sampleRate/(float TableSize)/frequency; if(phase >= (float)TableSize) phase -= (float)TableSize; #ifdef FM_OSCILLATOR if(phase < 0.f) phase += (float)TableSize; #endif return table[i] ; } float Oscillator::UpdateWithLinearInterpolation(float frequency) { int i = (int) phase; float alpha = phase - (float) i; phase += (sampleRate/(float)TableSize)/frequency; if(phase >= (float)TableSize) phase -= (float)TableSize; #ifdef FM_OSCILLATOR if(phase < 0.f) phase += (float)TableSize; #endif /* dtable0[i] = table[i+1] - table[i]; //remember to wrap around!!! */ return table[i] + dtable0[i]*alpha; } float Oscillator::UpdateWithCubicInterpolation( float frequency ) { int i = (int) phase; float alpha = phase - (float) i; phase += (sampleRate/(float)TableSize)/frequency; if(phase >= (float)TableSize) phase -= (float)TableSize; #ifdef FM_OSCILLATOR if(phase < 0.f) phase += (float)TableSize; #endif /* //remember to wrap around!!! dtable1[i] = (3.f*(table[i]-table[i+1])-table[i-1]+table[i+2])/2.f dtable2[i] = 2.f*table[i+1]+table[i-1]-(5.f*table[i]+table[i+2])/2.f dtable3[i] = (table[i+1]-table[i-1])/2.f */ return ((dtable1[i]*alpha + dtable2[i])*alpha + dtable3[i])*alpha+table[i]; } ```

### AM Formantic Synthesis¶

• Author or source: Paul Sernine
• Created: 2006-07-05 20:14:14
notes
```Here is another tutorial from Doc Rochebois.
It performs formantic synthesis without filters and without grains. Instead, it uses
"double carrier amplitude modulation" to pitch shift formantic waveforms. Just beware the
phase relationships to avoid interferences. Some patches of the DX7 used the same trick
but phase interferences were a problem. Here, Thierry Rochebois avoids them by using
cosine-phased waveforms.

Various formantic waveforms are precalculated and put in tables, they correspond to
different formant widths.
The runtime uses many intances (here 4) of these and pitch shifts them with double
carriers (to preserve the harmonicity of the signal).

This is a tutorial code, it can be optimized in many ways.
Have Fun

Paul
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137``` ```// FormantsAM.cpp // Thierry Rochebois' "Formantic Synthesis by Double Amplitude Modulation" // Based on a tutorial by Thierry Rochebois. // Comments by Paul Sernine. // The spectral content of the signal is obtained by adding amplitude modulated formantic // waveforms. The amplitude modulations spectraly shift the formantic waveforms. // Continuous spectral shift, without losing the harmonic structure, is obtained // by using crossfaded double carriers (multiple of the base frequency). // To avoid unwanted interference artifacts, phase relationships must be of the // "cosine type". // The output is a 44100Hz 16bit stereo PCM file. #include #include #include //Approximates cos(pi*x) for x in [-1,1]. inline float fast_cos(const float x) { float x2=x*x; return 1+x2*(-4+2*x2); } //Length of the table #define L_TABLE (256+1) //The last entry of the table equals the first (to avoid a modulo) //Maximal formant width #define I_MAX 64 //Table of formants float TF[L_TABLE*I_MAX]; //Formantic function of width I (used to fill the table of formants) float fonc_formant(float p,const float I) { float a=0.5f; int hmax=int(10*I)>L_TABLE/2?L_TABLE/2:int(10*I); float phi=0.0f; for(int h=1;hI_MAX-2?I_MAX-2:i; // width limitation float P=(L_TABLE-1)*(p+1)*0.5f; // phase normalisation int P0=(int)P; float fP=P-P0; // Integer and fractional int I0=(int)i; float fI=i-I0; // parts of the phase (p) and width (i). int i00=P0+L_TABLE*I0; int i10=i00+L_TABLE; //bilinear interpolation. return (1-fI)*(TF[i00] + fP*(TF[i00+1]-TF[i00])) + fI*(TF[i10] + fP*(TF[i10+1]-TF[i10])); } // Double carrier. // h : position (float harmonic number) // p : phase float porteuse(const float h,const float p) { float h0=floor(h); //integer and float hf=h-h0; //decimal part of harmonic number. // modulos pour ramener p*h0 et p*(h0+1) dans [-1,1] float phi0=fmodf(p* h0 +1+1000,2.0f)-1.0f; float phi1=fmodf(p*(h0+1)+1+1000,2.0f)-1.0f; // two carriers. float Porteuse0=fast_cos(phi0); float Porteuse1=fast_cos(phi1); // crossfade between the two carriers. return Porteuse0+hf*(Porteuse1-Porteuse0); } int main() { //Formant table for various french vowels (you can add your own) float F1[]={ 730, 200, 400, 250, 190, 350, 550, 550, 450}; float A1[]={ 1.0f, 0.5f, 1.0f, 1.0f, 0.7f, 1.0f, 1.0f, 0.3f, 1.0f}; float F2[]={ 1090, 2100, 900, 1700, 800, 1900, 1600, 850, 1100}; float A2[]={ 2.0f, 0.5f, 0.7f, 0.7f,0.35f, 0.3f, 0.5f, 1.0f, 0.7f}; float F3[]={ 2440, 3100, 2300, 2100, 2000, 2500, 2250, 1900, 1500}; float A3[]={ 0.3f,0.15f, 0.2f, 0.4f, 0.1f, 0.3f, 0.7f, 0.2f, 0.2f}; float F4[]={ 3400, 4700, 3000, 3300, 3400, 3700, 3200, 3000, 3000}; float A4[]={ 0.2f, 0.1f, 0.2f, 0.3f, 0.1f, 0.1f, 0.3f, 0.2f, 0.3f}; float f0,dp0,p0=0.0f; int F=7; //number of the current formant preset float f1,f2,f3,f4,a1,a2,a3,a4; f1=f2=f3=f4=100.0f;a1=a2=a3=a4=0.0f; init_formant(); FILE *f=fopen("sortie.pcm","wb"); for(int ns=0;ns<10*44100;ns++) { if(0==(ns%11025)){F++;F%=8;} //formant change f0=12*powf(2.0f,4-4*ns/(10*44100.0f)); //sweep f0*=(1.0f+0.01f*sinf(ns*0.0015f)); //vibrato dp0=f0*(1/22050.0f); float un_f0=1.0f/f0; p0+=dp0; //phase increment p0-=2*(p0>1); { //smoothing of the commands. float r=0.001f; f1+=r*(F1[F]-f1);f2+=r*(F2[F]-f2);f3+=r*(F3[F]-f3);f4+=r*(F4[F]-f4); a1+=r*(A1[F]-a1);a2+=r*(A2[F]-a2);a3+=r*(A3[F]-a3);a4+=r*(A4[F]-a4); } //The f0/fn coefficients stand for a -3dB/oct spectral enveloppe float out= a1*(f0/f1)*formant(p0,100*un_f0)*porteuse(f1*un_f0,p0) +0.7f*a2*(f0/f2)*formant(p0,120*un_f0)*porteuse(f2*un_f0,p0) + a3*(f0/f3)*formant(p0,150*un_f0)*porteuse(f3*un_f0,p0) + a4*(f0/f4)*formant(p0,300*un_f0)*porteuse(f4*un_f0,p0); short s=short(15000.0f*out); fwrite(&s,2,1,f);fwrite(&s,2,1,f); //fichier raw pcm stereo } fclose(f); return 0; } ```

• Date: 2007-04-24 12:04:12
• By: Baltazar
```Quite interesting and efficient for an algo that does not use any filter ;-)
```
• Date: 2007-08-14 11:30:14
• By: phoenix-69
```Very funny sound !
```
• Date: 2008-08-19 20:51:30
• By: Wait.
```What header files are you including?
```

### Alias-free waveform generation with analog filtering¶

notes
```(see linkfile)
```

### Another LFO class¶

notes
```This LFO uses an unsigned 32-bit phase and increment whose 8 Most Significant Bits adress
a Look-up table while the 24 Least Significant Bits are used as the fractionnal part.
Note: As the phase overflow automatically, the index is always in the range 0-255.

It performs linear interpolation, but it is easy to add other types of interpolation.

Don't know how good it could be as an oscillator, but I found it good enough for a LFO.
BTW there is also different kind of waveforms.

Modifications:
We could use phase on 64-bit or change the proportion of bits used by the index and the
fractionnal part.
```

```This type of oscillator is know as a numerically controlled oscillator(nco) or phase accumulator sythesiser. Integrated circuits that implement it in hardware are available such as the AD7008 from Analog Devices.
The frequency resolution is very high and is = (SampleRate)/32^2. So if clocked at 44.1Khz the frequency resolution would be 0.00001026Hz!
As you said the output waveform can be whatever shape you choose to put in the lookup table. The phase register is already in saw tooth form.

Regards,
Tony
```
```It works great!
Here's a Delphi version I just knocked up. Both VCL and KOL supported.

code:
unit PALFO;
//
// purpose: LUT based LFO
// Remarks: Translated from c++ sources by Remy Mueller, www.musicdsp.org

interface
uses
{\$IFDEF KOL}
Windows, Kol,KolMath;
{\$ELSE}
Windows, math;
{\$ENDIF}

const
k1Div24lowerBits = 1/(1 shl 24);

WFStrings:array[0..4] of string =
('triangle','sinus', 'sawtooth', 'square', 'exponent');

type
Twaveform = (triangle, sinus, sawtooth, square, exponent);

{\$IFDEF KOL}
PPaLfo = ^TPALfo;
TPaLfo = object(TObj)
{\$ELSE}
TPaLfo = class
{\$ENDIF}
private
FTable:array[0..256] of Single;// 1 more for linear interpolation
FPhase,
FInc:Single;
FRate: Single;
FSampleRate: Single;
FWaveForm: TWaveForm;
procedure SetRate(const Value: Single);
procedure SetSampleRate(const Value: Single);
procedure SetWaveForm(const Value: TWaveForm);
public
{\$IFNDEF KOL}
constructor create(SampleRate:Single);virtual;
{\$ENDIF}
// increments the phase and outputs the new LFO value.
// return the new LFO value between [-1;+1]
function WaveformName:String;
function Tick:Single;
// The rate in Hz
property Rate:Single read FRate write SetRate;
// The Samplerate
property SampleRate:Single read FSampleRate write SetSampleRate;
property WaveForm:TWaveForm read FWaveForm write SetWaveForm;
end;

{\$IFDEF KOL}
function NewPaLfo(aSamplerate:Single):PPaLfo;
{\$ENDIF}

implementation

{ TPaLfo }
{\$IFDEF KOL}
function NewPaLfo(aSamplerate:Single):PPaLfo;
begin
New(Result,Create);
with Result^ do
begin
FPhase:=0;
Finc:=0;
FSamplerate:=aSamplerate;
SetWaveform(triangle);
FRate:=1;
end;
end;
{\$ELSE}
constructor TPaLfo.create(SampleRate: Single);
begin
inherited create;
FPhase:=0;
Finc:=0;
FSamplerate:=aSamplerate;
SetWaveform(triangle);
FRate:=1;
end;
{\$ENDIF}

procedure TPaLfo.SetRate(const Value: Single);
begin
FRate := Value;
// the rate in Hz is converted to a phase increment with the following formula
// f[ inc = (256*rate/samplerate) * 2^24]
Finc :=  (256 * Frate / Fsamplerate) * (1 shl 24);
end;

procedure TPaLfo.SetSampleRate(const Value: Single);
begin
FSampleRate := Value;
end;

procedure TPaLfo.SetWaveForm(const Value: TWaveForm);
var
i:integer;
begin
FWaveForm := Value;
Case Fwaveform of
sinus:
for i:=0 to 256 do
FTable[i] := sin(2*pi*(i/256));
triangle:
begin
for i:=0 to 63 do
begin
FTable[i] := i / 64;
FTable[i+64] :=(64-i) / 64;
FTable[i+128] := - i / 64;
FTable[i+192] := - (64-i) / 64;
end;
FTable := 0;
end;
sawtooth:
begin
for i:=0 to 255 do
FTable[i] := 2*(i/255) - 1;
FTable := -1;
end;
square:
begin
for i:=0 to 127 do
begin
FTable[i]     :=  1;
FTable[i+128] := -1;
end;
FTable := 1;
end;
exponent:
begin
// symetric exponent similar to triangle
for i:=0 to 127 do
begin
FTable[i] := 2 * ((exp(i/128) - 1) / (exp(1) - 1)) - 1  ;
FTable[i+128] := 2 * ((exp((128-i)/128) - 1) / (exp(1) - 1)) - 1  ;
end;
FTable := -1;
end;
end;
end;

function TPaLfo.WaveformName:String;
begin
result:=WFStrings[Ord(Fwaveform)];
end;

function TPaLfo.Tick: Single;
var
i:integer;
frac:Single;
begin
// the 8 MSB are the index in the table in the range 0-255
i := Pinteger(Fphase)^ shr 24;
// and the 24 LSB are the fractionnal part
frac := (PInteger(Fphase)^ and \$00FFFFFF) * k1Div24lowerBits;
// increment the phase for the next tick
Fphase :=FPhase + Finc; // the phase overflows itself
Result:= Ftable[i]*(1-frac) + Ftable[i+1]* frac; // linear interpolation
end;

end.
```
```Oops,

This one is correct:

code:
unit PALFO;
//
// purpose: LUT based LFO
// Remarks: Translated from c++ sources by Remy Mueller, www.musicdsp.org

interface
uses
{\$IFDEF KOL}
Windows, Kol,KolMath;
{\$ELSE}
Windows, math;
{\$ENDIF}

const
k1Div24lowerBits = 1/(1 shl 24);

WFStrings:array[0..4] of string =
('triangle','sinus', 'sawtooth', 'square', 'exponent');

type
Twaveform = (triangle, sinus, sawtooth, square, exponent);

{\$IFDEF KOL}
PPaLfo = ^TPALfo;
TPaLfo = object(TObj)
{\$ELSE}
TPaLfo = class
{\$ENDIF}
private
FTable:array[0..256] of Single;// 1 more for linear interpolation
FPhase,
FInc:dword;
FRate: Single;
FSampleRate: Single;
FWaveForm: TWaveForm;
procedure SetRate(const Value: Single);
procedure SetSampleRate(const Value: Single);
procedure SetWaveForm(const Value: TWaveForm);
public
{\$IFNDEF KOL}
constructor create(SampleRate:Single);virtual;
{\$ENDIF}
// increments the phase and outputs the new LFO value.
// return the new LFO value between [-1;+1]
function WaveformName:String;
function Tick:Single;
// The rate in Hz
property Rate:Single read FRate write SetRate;
// The Samplerate
property SampleRate:Single read FSampleRate write SetSampleRate;
property WaveForm:TWaveForm read FWaveForm write SetWaveForm;
end;

{\$IFDEF KOL}
function NewPaLfo(aSamplerate:Single):PPaLfo;
{\$ENDIF}

implementation

{ TPaLfo }
{\$IFDEF KOL}
function NewPaLfo(aSamplerate:Single):PPaLfo;
begin
New(Result,Create);
with Result^ do
begin
FPhase:=0;
FSamplerate:=aSamplerate;
SetWaveform(sinus);
Rate:=1;
end;
end;
{\$ELSE}
constructor TPaLfo.create(SampleRate: Single);
begin
inherited create;
FPhase:=0;
FSamplerate:=aSamplerate;
SetWaveform(sinus);
FRate:=1;
end;
{\$ENDIF}

procedure TPaLfo.SetRate(const Value: Single);
begin
FRate := Value;
// the rate in Hz is converted to a phase increment with the following formula
// f[ inc = (256*rate/samplerate) * 2^24]
Finc :=  round((256 * Frate / Fsamplerate) * (1 shl 24));
end;

procedure TPaLfo.SetSampleRate(const Value: Single);
begin
FSampleRate := Value;
end;

procedure TPaLfo.SetWaveForm(const Value: TWaveForm);
var
i:integer;
begin
FWaveForm := Value;
Case Fwaveform of
sinus:
for i:=0 to 256 do
FTable[i] := sin(2*pi*(i/256));
triangle:
begin
for i:=0 to 63 do
begin
FTable[i] := i / 64;
FTable[i+64] :=(64-i) / 64;
FTable[i+128] := - i / 64;
FTable[i+192] := - (64-i) / 64;
end;
FTable := 0;
end;
sawtooth:
begin
for i:=0 to 255 do
FTable[i] := 2*(i/255) - 1;
FTable := -1;
end;
square:
begin
for i:=0 to 127 do
begin
FTable[i]     :=  1;
FTable[i+128] := -1;
end;
FTable := 1;
end;
exponent:
begin
// symetric exponent similar to triangle
for i:=0 to 127 do
begin
FTable[i] := 2 * ((exp(i/128) - 1) / (exp(1) - 1)) - 1  ;
FTable[i+128] := 2 * ((exp((128-i)/128) - 1) / (exp(1) - 1)) - 1  ;
end;
FTable := -1;
end;
end;
end;

function TPaLfo.WaveformName:String;
begin
result:=WFStrings[Ord(Fwaveform)];
end;

function TPaLfo.Tick: Single;
var
i:integer;
frac:Single;
begin
// the 8 MSB are the index in the table in the range 0-255
i := Fphase shr 24;
// and the 24 LSB are the fractionnal part
frac := (Fphase and \$00FFFFFF) * k1Div24lowerBits;
// increment the phase for the next tick
Fphase :=FPhase + Finc; // the phase overflows itself
Result:= Ftable[i]*(1-frac) + Ftable[i+1]* frac; // linear interpolation
end;

end.
```

### Another cheap sinusoidal LFO¶

notes
```Some pseudo code for a easy to calculate LFO.

You can even make a rough triangle wave out of this by substracting the output of 2 of
these with different phases.

PJ
```
code
 ```1 2 3 4 5 6 7``` ```r = the rate 0..1 -------------- p += r if(p > 1) p -= 2; out = p*(1-abs(p)); -------------- ```

```Slick! I like it!

Sincerely,
Frederick Umminger
```
```Great! just what I wanted a fast trinagle lfo. :D

float rate = 0..1;
float phase1 = 0;
float phase2 = 0.1f;

---------------
phase1 += rate;
if (phase1>1) phase1 -= 2;

phase2 += rate;
if (phase2>1) phase2 -= 2;

out = (phase1*(1-abs(phase1)) - phase2*(1-abs(phase2))) * 10;
---------------
```
```Nice! If you want the output range to be between -1..1 then use:

--------------
p += r
if(p > 2) p -= 4;
out = p*(2-abs(p));
--------------
```
```A better way of making a triangle LFO (out range is -1..1):

rate = 0..1;
p = -1;
{
p += rate;
if (p>1) p -= 4.0f;
out = abs(-(abs(p)-2))-1;
}
```
```/* this goes from -1 to +1 */
#include <iostream>
#include <math.h>

using namespace std;
int main(int argc, char *argv[]) {

//r = the rate 0..1
float r = 0.5f;
float p = 0.f;
float result=0.f;
//--------------
for(int i=1;i<=2048;i++){
p += r;
if(p > 1) p -= 2;
result = 4*p*(1-fabs(p));
cout << result;
cout <<"\r";
}
}
```

### Antialiased square generator¶

• Author or source: Paul Sernine
• Type: 1st April edition
• Created: 2006-04-01 12:46:23
notes
```It is based on a code by Thierry Rochebois, obfuscated by me.
It generates a 16bit MONO raw pcm file. Have Fun.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21``` ```//sqrfish.cpp #include #include //obfuscation P.Sernine int main() {float ccc,cccc=0,CC=0,cc=0,CCCC, CCC,C,c; FILE *CCCCCCC=fopen("sqrfish.pcm", "wb" ); int ccccc= 0; float CCCCC=6.89e-6f; for(int CCCCCC=0;CCCCCC<1764000;CCCCCC++ ){ if(!(CCCCCC%7350)){if(++ccccc>=30){ ccccc =0; CCCCC*=2;}CCC=1;}ccc=CCCCC*expf(0.057762265f* "aiakahiafahadfaiakahiahafahadf"[ccccc]);CCCC =0.75f-1.5f*ccc;cccc+=ccc;CCC*=0.9999f;cccc-= 2*(cccc>1);C=cccc+CCCC*CC; c=cccc+CCCC*cc; C -=2*(C>1);c-=2*(c>1);C+=2*(C<-1); c+=1+2 *(c<-1);c-=2*(c>1);C=C*C*(2 *C*C-4); c=c*c*(2*c*c-4); short cccccc=short(15000.0f* CCC*(C-c )*CCC);CC=0.5f*(1+C+CC);cc=0.5f*(1+ c+cc); fwrite(&cccccc,2,1,CCCCCCC);} //algo by Thierry Rochebois fclose(CCCCCCC); return 0000000;} ```

### Arbitary shaped band-limited waveform generation (using oversampling and low-pass filtering)¶

code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16``` ```Arbitary shaped band-limited waveform generation (using oversampling and low-pass filtering) There are many articles about band-limited waveform synthesis techniques, that provide correct and fast methods for generating classic analogue waveforms, such as saw, pulse, and triangle wave. However, generating arbitary shaped band-limited waveforms, such as the "sawsin" shape (found in this source-code archive), seems to be quite hard using these techniques. My analogue waveforms are generated in a _very_ high sampling rate (actually it's 1.4112 GHz for 44.1 kHz waveforms, using 32x oversampling). Using this sample-rate, the amplitude of the aliasing harmonics are negligible (the base analogue waveforms has exponentially decreasing harmonics amplitudes). Using a 511-tap windowed sync FIR filter (with Blackman-Harris window, and 12 kHz cutoff frequency) the harmonics above 20 kHz are killed, the higher harmonics (that cause the sharp overshoot at step response) are dampened. The filtered signal downsampled to 44.1 kHz contains the audible (non-aliased) harmonics only. This waveform synthesis is performed for wavetables of 4096, 2048, 1024, ... 8, 4, 2 samples. The real-time signal is interpolated from these waveform-tables, using Hermite-(cubic-)interpolation for the waveforms, and linear interpolation between the two wavetables near the required note. This procedure is quite time-consuming, but the whole waveform (or, in my implementation, the whole waveform-set) can be precalculated (or saved at first launch of the synth) and reloaded at synth initialization. I don't know if this is a theoretically correct solution, but the waveforms sound good (no audible aliasing). Please let me know if I'm wrong... ```

```Why can't you use fft/ifft
to synthesis directly wavetables of 2048,1024,..?
It'd be not so
"time consuming" comparing to FIR filtering.
Further cubic interpolation still might give you audible
distortion in some cases.
--Alex.
```
```What should I use instead of cubic interpolation? (I had already some aliasing problems with cubic interpolation, but that can be solved by oversampling 4x the realtime signal generation)
Is this theory of generating waves from wavetables of 4096, 2084, ... 8, 4, 2 samples wrong?
```
```I think tablesize should not vary
depending on tone (4096,2048...)
and you'd better stay with the same table size for all notes (for example 4096, 4096...).

To avoid interpolation noise
(it's NOT caused by aliasing)
try to increase wavetable size
and be sure that waveform spectrum has
steep roll off
(don't forget Gibbs phenomena as well).
```
```you say that the higher harmonics (that cause the sharp overshoot at step response) are dampened.
How ? Or is it a result of the filtering ?
```
```Yes. The FIR-filter cutoff is set to 12 kHz, so it dampens the audible frequencies too. This way the frequncies above 20 kHz are about -90 dB (don't remember exactly, but killing all harmonics above 20 kHz was the main reason to set the cutoff to 12 kHz).

Anyway, as Alex suggested, FFT/IFFT seems to be a better solution to this problem.
```

### Audiable alias free waveform gen using width sine¶

notes
```Warning, my english abilities is terribly limited.

How ever, the other day when finally understanding what bandlimited wave creation is (i am
a noobie, been doing DSP stuf on and off for a half/year) it hit me i can implement one
little part in my synths. It's all about the freq (that i knew), very simple you can
reduce alias (the alias that you can hear that is) extremely by keeping track of your
frequence, the way i solved it is using a factor, afact = 1 - sin(f*2PI). This means you
can do audiable alias free synthesis without very complex algorithms or very huge tables,
even though the sound becomes kind of low-filtered.
Propably something like this is mentioned b4, but incase it hasn't this is worth looking
up

The psuedo code describes it more.

// Druttis
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24``` ```f := freq factor, 0 - 0.5 (0 to half samplingrate) afact(f) = 1 - sin(f*2PI) t := time (0 to ...) ph := phase shift (0 to 1) fm := freq mod (0 to 1) sine(t,f,ph,fm) = sin((t*f+ph)*2PI + 0.5PI*fm*afact(f)) fb := feedback (0 to 1) (1 max saw) saw(t,f,ph,fm,fb) = sine(t,f,ph,fb*sine(t-1,f,ph,fm)) pm := pulse mod (0 to 1) (1 max pulse) pw := pulse width (0 to 1) (1 square) pulse(t,f,ph,fm,fb,pm,pw) = saw(t,f,ph,fm,fb) - (t,f,ph+0.5*pw,fm,fb) * pm I am not completely sure about fm for saw & pulse since i cant test that atm. but it should work :) otherwise just make sure fm are 0 for saw & pulse. As you can see the saw & pulse wave are very variable. // Druttis ```

```Um, reading it I can see a big flaw...

afact(f) = 1 - sin(f*2PI) is not correct!

should be

afact(f) = 1 - sqrt(f * 2 / sr)

where sr := samplingrate
f should be exceed half sr
```
```f has already be divided by sr, right ? So it should become :

afact (f) = 1 - sqrt (f  * 2)

And i see a typo (saw forgotten in the second expression) :

pulse(t,f,ph,fm,fb,pm,pw) = saw(t,f,ph,fm,fb) - saw(t,f,ph+0.5*pw,fm,fb) * pm

However I haven't checked the formula.
```
```Hi Lauent,
I'm new to that DSP stuff and can't get the key to
what'S the meaning of afact? - Can you explain please!? - Thanks in advice!
```
```I've been playing around with this for some time. Expect a major update in a while, as soon as I know how to describe it :)
```
```afact is used as an amplitude factor for fm or fb depending on the carrier frequency. The higher frequency the lower afact. It's not completely resolving the problem with aliasing but it is a cheap way that dramatically reduces it.
```

### Band Limited waveforms my way¶

• Author or source: Anton Savov (gb.liam@ottna)
• Type: classic Sawtooth example
• Created: 2009-06-22 18:09:08
notes
```This is my <ugly> C++ code for generating a single cycle of a Sawtooth in a table
normaly i create my "fundamental" table big enough to hold on around 20-40Hz in the
current Sampling rate
also, i create the table twice as big, i do "mip-maps" then
so the size should be a power of two, say 1024 for 44100Hz = 44100/1024 = ~43.066Hz
then the mip-maps are with decreasing sizes (twice) 512, 256, 128, 64, 32, 16, 8, 4, and 2

if the "gibbs" effect is what i think it is - then i have a simple solution
here is my crappy code:
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35``` ```int sz = 1024; // the size of the table int i = 0; float *table; // pointer to the table double scale = 1.0; double pd; // phase double omega = 1.0 / (double)(sz); while (i < sz) { double amp = scale; double x = 0.0; // the sample double h = 1; // harmonic number (starts from 1) double dd; // fix high frequency "ring" pd = (double)(i) / (double)(sz); // calc phase double hpd = pd; // phase of the harmonic while (true) // start looping for this sample { if ((omega * h) < 0.5) // harmonic frequency is in range? { dd = cos(omega * h * 2 * pi); x = x + (amp * dd * sin(hpd * 2 * pi)); h = h + 1; hpd = pd * h; amp = 1.0 / h; } else { break; } } table[i] = x; ++i; } the peaks are around +/- 0.8 a square can be generated by just changing h = h+2; the peaks would be +/- 0.4 any bugs/improvements? ```

```excuse me, there is a typo

amp = scale / h;
```
• Date: 2009-09-20 10:34:18
• By: antto mail bg
```for even smoother edges:
dd = cos(sin(omega * h * pi) * 0.5 * pi);
no visual ringing, smooth waveform
```
• Date: 2009-09-25 04:15:31
• By: antto mail bg
```to get +/- 1.0 amplitude: scale = 1.25
```

### Bandlimited sawtooth synthesis¶

notes
```This is working code for synthesizing a bandlimited sawtooth waveform. The algorithm is
DSF BLIT + leaky integrator. Includes driver code.

There are two parameters you may tweak:

1) Desired attenuation at nyquist. A low value yields a duller sawtooth but gets rid of
those annoying CLICKS when sweeping the frequency up real high. Must be strictly less than
1.0!

2) Integrator leakiness/cut off. Affects the shape of the waveform to some extent, esp. at
the low end. Ideally you would want to set this low, but too low a setting will give you
problems with DC.

Have fun!
/Emanuel Landeholm

```

```there is no need to use a butterworth design for a simple leaky integrator, in this case actually the
variable curcps can be used directly in a simple: leak += curcps * (blit - leak);

this produces a nearly perfect saw shape in almost all cases
```
```The square wave type will be able to be generated from this source.
Please teach if it is possible.
```

### Bandlimited waveform generation¶

notes
```(see linkfile)
```

```The code to reduce the gibbs effect causes aliasing if a transition is made from wavetable A with x partials to wavetable B with y partials.

The aliasing can clearly be seen in a spectral view.

The problem is, that the volume modification for partial N is different depending on the number of partials the wavetable row contains
```

### Bandlimited waveforms synopsis.¶

notes
```(see linkfile)
```

• Date: 2005-11-15 20:07:11
```The abs(sin) method from the Lane CMJ paper is not bandlimited! It's basically just a crappy method for BLIT.

You forgot to mention Eli Brandt's minBLEP method. It's the best! You just have to know how to properly generate a nice minblep table... (slightly dilated, see Stilson and Smith BLIT paper, at the end regarding table implementation issues)
```

### Bandlimited waveforms…¶

• Author or source: Paul Kellet
• Created: 2002-01-17 00:56:04
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.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34``` ```float p=0.0f; //current position float dp=1.0f; //change in postion per sample float pmax; //maximum position float x; //position in sinc function float leak=0.995f; //leaky integrator float dc; //dc offset float saw; //output //set frequency... pmax = 0.5f * getSampleRate() / freqHz; dc = -0.498f/pmax; //for each sample... p += dp; if(p < 0.0f) { p = -p; dp = -dp; } else if(p > pmax) { p = pmax + pmax - p; dp = -dp; } x= pi * p; if(x < 0.00001f) x=0.00001f; //don't divide by 0 saw = leak*saw + dc + (float)sin(x)/(x); ```

```Hi,
Has anyone managed to implement this in a VST?
If anyone could mail me and talk me through it I'd be very grateful.  Yes, I'm a total newbie and yes, I'm after a quick-fix solution...we all have to start somewhere, eh?

As it stands, where I should be getting a sawtooth I'm getting a full-on and inaudible signal...!

Even a small clue would be nice.
Cheers,
A
```
```this sounds quite nice, maybe going to use it an LV 2 plugin
```
• Date: 2016-01-17 13:24:11
• By: pvdmeer [atorsomething] gmail [point] com
```this is really amazing, and easily hacked into a lut-based algo. i'll try windowing it too, but it already looks like aliasing is well within acceptable levels.
```

### Butterworth¶

code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37``` ```First calculate the prewarped digital frequency: K = tan(Pi * Frequency / Samplerate); Now calc some intermediate variables: (see 'Factors of Polynoms' at http://en.wikipedia.org/wiki/Butterworth_filter, especially if you want a higher order like 48dB/Oct) a = 0.76536686473 * Q * K; b = 1.84775906502 * Q * K; K = K*K; (to optimize it a little bit) Calculate the first biquad: A0 = (K+a+1); A1 = 2*(1-K); A2 =(a-K-1); B0 = K; B1 = 2*B0; B2 = B0; Calculate the second biquad: A3 = (K+b+1); A4 = 2*(1-K); A5 = (b-K-1); B3 = K; B4 = 2*B3; B5 = B3; Then calculate the output as follows: Stage1 = B0*Input + State0; State0 = B1*Input + A1/A0*Stage1 + State1; State1 = B2*Input + A2/A0*Stage1; Output = B3*Stage1 + State2; State2 = B4*Stage1 + A4/A3*Output + State2; State3 = B5*Stage1 + A5/A3*Output; ```

```If you have a Q factor different than 1, then filter won't be a Butterworth filter (in terms of maximally flat passpand). So, your filter is a kind of a tweaked Butterworth filter with added resonance.

Highpass version should be:

A0 = (K+a+1);
A1 = 2*(1-K);
A2 =(a-K-1);
B0 = 1;
B1 = -2;
B2 = 1;

A3 = (K+b+1);
A4 = 2*(1-K);
A5 = (b-K-1);
B3 = 1;
B4 = -2;
B5 = 1;

The rest is the same. You might want to leave out B0, B2, B3 and B5 completely, because they all equal to 1, and optimize the highpass loop as:

Stage1 = Input + State0;
State0 = B1*Input + A1/A0*Stage1 + State1;
State1 = Input + A2/A0*Stage1;

Output = Stage1 + State2;
State2 = B4*Stage1 + A4/A3*Output + State2;
State3 = Stage1 + A5/A3*Output;

Anyone confirms this code works? (Being too lazy to throw this into a compiler...)

Cheers,
Peter
```
```And of course it is not a good idea to do divisions in the process loop, because they are very heavy, so the best is to precalculate A1/A0, A2/A0, A4/A3 and A5/A3 after the calculation of coefficients:

inv_A0 = 1.0/A0;
A1A0 = A1 * inv_A0;
A2A0 = A2 * inv_A0;

inv_A3 = 1.0/A3;
A4A3 = A4 * inv_A3;
A5A3 = A5 * inv_A3;

(The above should be faster than writing

A1A0 = A1/A0;
A2A0 = A2/A0;
A4A3 = A4/A3;
A5A3 = A5/A3;

but I think some compilers do this optimization automatically.)

Then the lowpass process loop becomes

Stage1 = B0*Input + State0;
State0 = B1*Input + A1A0*Stage1 + State1;
State1 = B2*Input + A2A0*Stage1;

Output = B3*Stage1 + State2;
State2 = B4*Stage1 + A4A3*Output + State2;
State3 = B5*Stage1 + A5A3*Output;

Much faster, isn't it?
```
```Once you figured it out, it's even possible to do higher order butterworth shelving filters. Here's an example of an 8th order lowshelf.

First we start as usual prewarping the cutoff frequency:

K = tan(fW0*0.5);

Then we settle up the Coefficient V:

V = Power(GainFactor,-1/4)-1;

Finally here's the loop to calculate the filter coefficients:

for i = 0 to 3
{
cm = cos(PI*(i*2+1) / (2*8) );

B[3*i+0] = 1/ ( 1 + 2*K*cm + K*K + 2*V*K*K + 2*V*K*cm + V*V*K*K);
B[3*i+1] = 2 * ( 1 - K*K - 2*V*K*K - V*V*K*K);
B[3*i+2] = (-1 + 2*K*cm - K*K - 2*V*K*K + 2*V*K*cm - V*V*K*K);
A[3*i+0] = ( 1-2*K*cm+K*K);
A[3*i+1] = 2*(-1 +K*K);
A[3*i+2] = ( 1+2*K*cm+K*K);
}
```
```Hmm... interesting. I guess the phase response/group delay gets quite funky, which is generally unwanted for an equalizer.

I think the 1/ is not necessary for the first B coefficient! (of course you divide all the other coeffs with the inverse of that coeff at the end...)

I guess the next will be Chebyshev shelving filters ;)

BTW did you check whether my 4 pole highpass Butterworth code is correct?

Peter
```
```The 1/ is of course an error here. It's left of my own implementation, where I divide directly. Also I think A and B is exchanged.

I've already nearly done all the different filter types (except elliptic filters), but I won't post too much here.
The highpass maybe a highpass, but not the exact complementary. At least my (working) implementation looks different.

The lowpass<->highpass transform is to replace s with 1/s and by doing this, more than one sign is changing.
```
```Different authors tend to mix up A and B coeffs.

If I take this lowpass derived by bilinear transform and change B0 and B2 to 1, and B1 to -2 then I get a perfect highpass. At least that's what I see in FilterExplorer. Probably you could get the same by replacing s with 1/s and deriving it by bilinear transform.

Well, there are many ways to get a filter working, for example if I replace tan(PI*w) with 1.0/tan(PI*w), inverse the sign of B1, and replace A1 = 2*(1-K) with A1 = 2*(K-1), I also get the same highpass.

Well, the reason for the sign inversion is that the coeffs you named B1 and B2 here are responsible for the locations of the zeroes. B1 is responsible for the angle, and B2 is for the radius, so if you invert B1 then the zeroes get on the opposite side of the unit circle, so you get a highpass filter. You then need to adjust the gain coefficient (B0) so that the passband gain is 1. Well, this is not a very precise explanation, but this is the reason why this works.
```
```Ok, you're right. I've been doing too much stuff these days that I missed that simple thing. Your version is even numerical better, because there is less potential of annihilation. Thanks for that.

Btw. the group delay really get a little bit 'funky', i've also noticed that, but for not too high orders it doesn't hurt that much.
```
```Well, it isn't a big problem unless you start modulating the filter very fast... then you get this strange pitch-shifting effect ;)

Well, I read sometimes that when you do EQing, phase is a very important factor. I guess that's why ppl sell a lot of linear phase EQ plugins. Or just the marketing? Don't know, haven't compared linear and non-linear phase stuff very much..
```
```State2 = B4*Stage1 + A4/A3*Output + State2;
State2 = B4*Stage1 + A4/A3*Output + State3;
```

### C# Oscilator class¶

• Author or source: neotec
• Type: Sine, Saw, Variable Pulse, Triangle, C64 Noise
• Created: 2007-01-08 10:49:36
notes
```Parameters:

Pitch: The Osc's pitch in Cents [0 - 14399] startig at A -> 6.875Hz
Pulsewidth: [0 - 65535] -> 0% to 99.99%
Value: The last Output value, a set to this property 'syncs' the Oscilator
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194``` ```public class SynthOscilator { public enum OscWaveformType { SAW, PULSE, TRI, NOISE, SINE } public int Pitch { get { return this._Pitch; } set { this._Pitch = this.MinMax(0, value, 14399); this.OscStep = WaveSteps[this._Pitch]; } } public int PulseWidth { get { return this._PulseWidth; } set { this._PulseWidth = this.MinMax(0, value, 65535); } } public OscWaveformType Waveform { get { return this._WaveForm; } set { this._WaveForm = value; } } public int Value { get { return this._Value; } set { this._Value = 0; this.OscNow = 0; } } private int _Pitch; private int _PulseWidth; private int _Value; private OscWaveformType _WaveForm; private int OscNow; private int OscStep; private int ShiftRegister; public const double BaseFrequence = 6.875; public const int SampleRate = 44100; public static int[] WaveSteps = new int; public static int[] SineTable = new int; public SynthOscilator() { if (WaveSteps.Length == 0) this.CalcSteps(); if (SineTable.Length == 0) this.CalcSine(); this._Pitch = 7200; this._PulseWidth = 32768; this._WaveForm = OscWaveformType.SAW; this.ShiftRegister = 0x7ffff8; this.OscNow = 0; this.OscStep = WaveSteps[this._Pitch]; this._Value = 0; } private void CalcSteps() { WaveSteps = new int; for (int i = 0; i < 14400; i++) { double t0, t1, t2; t0 = Math.Pow(2.0, (double)i / 1200.0); t1 = BaseFrequence * t0; t2 = (t1 * 65536.0) / (double)this.SampleRate; WaveSteps[i] = (int)Math.Round(t2 * 4096.0); } } private void CalcSine() { SineTable = new int; double s = Math.PI / 32768.0; for (int i = 0; i < 65536; i++) { double v = Math.Sin((double)i * s) * 32768.0; int t = (int)Math.Round(v) + 32768; if (t < 0) t = 0; else if (t > 65535) t = 65535; SineTable[i] = t; } } public override int Run() { int ret = 32768; int osc = this.OscNow >> 12; switch (this._WaveForm) { case OscWaveformType.SAW: ret = osc; break; case OscWaveformType.PULSE: if (osc < this.PulseWidth) ret = 65535; else ret = 0; break; case OscWaveformType.TRI: if (osc < 32768) ret = osc << 1; else ret = 131071 - (osc << 1); break; case OscWaveformType.NOISE: ret = ((this.ShiftRegister & 0x400000) >> 11) | ((this.ShiftRegister & 0x100000) >> 10) | ((this.ShiftRegister & 0x010000) >> 7) | ((this.ShiftRegister & 0x002000) >> 5) | ((this.ShiftRegister & 0x000800) >> 4) | ((this.ShiftRegister & 0x000080) >> 1) | ((this.ShiftRegister & 0x000010) << 1) | ((this.ShiftRegister & 0x000004) << 2); ret <<= 4; break; case OscWaveformType.SINE: ret = SynthTools.SineTable[osc]; break; default: break; } this.OscNow += this.OscStep; if (this.OscNow > 0xfffffff) { int bit0 = ((this.ShiftRegister >> 22) ^ (this.ShiftRegister >> 17)) & 0x1; this.ShiftRegister <<= 1; this.ShiftRegister &= 0x7fffff; this.ShiftRegister |= bit0; } this.OscNow &= 0xfffffff; this._Value = ret - 32768; return this._Value; } public int MinMax(int a, int b, int c) { if (b < a) return a; else if (b > c) return c; else return b; } } ```

### C++ gaussian noise generation¶

notes
```References :
Tobybears delphi noise generator was the basis. Simply converted it to C++.
http://www.musicdsp.org/archive.php?classid=0#129
The output is in noise.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24``` ```/* Include requisits */ #include #include /* Generate a new random seed from system time - do this once in your constructor */ srand(time(0)); /* Setup constants */ const static int q = 15; const static float c1 = (1 << q) - 1; const static float c2 = ((int)(c1 / 3)) + 1; const static float c3 = 1.f / c1; /* random number in range 0 - 1 not including 1 */ float random = 0.f; /* the white noise */ float noise = 0.f; for (int i = 0; i < numSamples; i++) { random = ((float)rand() / (float)(RAND_MAX + 1)); noise = (2.f * ((random * c2) + (random * c2) + (random * c2)) - 3.f * (c2 - 1.f)) * c3; } ```

```What's the difference between the much simpler noise generator:

randSeed = (randSeed * 196314165) + 907633515;      out=((int)randSeed)*0.0000000004656612873077392578125f;

and this one? they both sound the same to my ears...
```
```How can you change the variance (sigma)?
```
```This is NOT a good code to generate Gaussian Noice. Look into:
(random * c2) + (random * c2) + (random * c2)
It is all nonsense! The reason of adding three numbers it the Central Limit Theorem to aproximate Gaussian distribution. But the random numbers inside must differ, which is not the case. The code on original link http://www.musicdsp.org/archive.php?classid=0#129 is correct.
```

### Chebyshev waveshaper (using their recursive definition)¶

• Author or source: mdsp
• Type: chebyshev
• Created: 2005-01-10 18:03:29
notes
```someone asked for it on kvr-audio.

I use it in an unreleased additive synth.
There's no oversampling needed in my case since I feed it with a pure sinusoid and I
control the order to not have frequencies above Fs/2. Otherwise you should oversample by
the order you'll use in the function or bandlimit the signal before the waveshaper. unless
you really want that aliasing effect... :)

I hope the code is self-explaining, otherwise there's plenty of sites explaining chebyshev
polynoms and their applications.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20``` ```float chebyshev(float x, float A[], int order) { // To = 1 // T1 = x // Tn = 2.x.Tn-1 - Tn-2 // out = sum(Ai*Ti(x)) , i C {1,..,order} float Tn_2 = 1.0f; float Tn_1 = x; float Tn; float out = A*Tn_1; for(int n=2;n<=order;n++) { Tn = 2.0f*x*Tn_1 - Tn_2; out += A[n-1]*Tn; Tn_2 = Tn_1; Tn_1 = Tn; } return out; } ```

• Date: 2005-01-10 18:10:12
• By: mdsp
```BTW you can achieve an interresting effect by feeding back the ouput in the input. it adds a kind of interresting pitched noise to the signal.

I think VirSyn is using something similar in microTERA.
```
```Hi, it was me that asked about this on KvR. It seems that it is possible to use such a waveshaper on a non-sinusoidal input without oversampling; split the input signal into bands, and use the highest frequency in each band to determine which order polynomials to send each band to. The idea about feeding back the output to the input occured to me as well, good to know that such an effect might be interesting... If I come across any other points of interest while coding this plugin, I'll be glad to mention them on here.

Dan
```
• Date: 2005-01-12 11:48:00
• By: mdsp
```of course you can use it on non sinusoidal input, but you won't achieve the same result.

if you express your input as a sum of sinusoids of frequencies [f0 f1 f2 ...] and use the chebyshev polynom of order 2 you won't have 2*[f0 f1 f2...] as the resulting frequencies.
As it's a nonlinear function you can't use the superposition theorem anymore.

beware that chebyshev polynoms are sensitive to
the range of your input. Your sinusoid has to have a gain exaclty equal to 1 in order to work as expected.

that's a nice trick but it has it's limits.
```

### Cubic polynomial envelopes¶

• Author or source: Andy Mucho
• Type: envellope generation
• Created: 2002-01-17 00:59:16
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17``` ```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

```I have try this and it works fine, but what hell is bigs?????

bye bye

float time = (float)pRect.Width();                 //time in sampleframes
float startlevel = (float)pRect.Height();  //max h vedi ma 1.0
float midlevel = 500.f;
float endlevel = 0.f;

float k = startlevel + endlevel - (midlevel * 2);
float r = startlevel;
float s = (endlevel - startlevel - (2 * k)) / time;
float t= (2 * k) / (time * time);

float bigr = r;
float bigs = s + t;
float bigt = 2 * t;

for(int i=0;i<time;i++)
{
bigr = bigr + bigs;
bigs = bigs + bigt;                                                     //bigs? co'è
pDC->SetPixel(i,(int)bigr,RGB (0, 0, 0));
}
```
```the method uses a technique called forward differencing, which is based on the fact that a successive values of an polynomial function can be calculated using only additions instead of evaluating the whole polynomial which would require huge amount of multiplications.

Actually the method presented here uses only a quadratic curve, not cubic. The number of the variables in the adder is N+1, where N is the order of the polynomial to be generated. In this example we have only three, thus second order function. For linear we would have two variables: the current value and the constant adder.

The trickiest part is to set up the adder variables...

```

### DSF (super-set of BLIT)¶

• Author or source: David Lowenfels
• Type: matlab code
• Created: 2003-04-02 23:59:24
notes
```Discrete Summation Formula ala Moorer

computes equivalent to sum{k=0:N-1}(a^k * sin(beta + k*theta))
modified from Emanuel Landeholm's C code
output should never clip past [-1,1]

If using for BLIT synthesis for virtual analog:
N = maxN;
a = attn_at_Nyquist ^ (1/maxN); %hide top harmonic popping in and out when sweeping
frequency
beta = pi/2;
num = 1 - a^N * cos(N*theta) - a*( cos(theta) - a^N * cos(N*theta - theta) ); %don't waste
time on beta

You can also get growing harmonics if a > 1, but the min statement in the code must be
removed, and the scaling will be weird.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26``` ```function output = dsf( freq, a, H, samples, beta) %a = rolloff coeffecient %H = number of harmonic overtones (fundamental not included) %beta = harmonic phase shift samplerate = 44.1e3; freq = freq/samplerate; %normalize frequency bandlimit = samplerate / 2; %Nyquist maxN = 1 + floor( bandlimit / freq ); %prevent aliasing N = min(H+2,maxN); theta = 2*pi * phasor(freq, samples); epsilon = 1e-6; a = min(a, 1-epsilon); %prevent divide by zero num = sin(beta) - a*sin(beta-theta) - a^N*sin(beta + N*theta) + a^(N+1)*sin(beta+(N-1)*theta); den = (1 + a * ( a - 2*cos(theta) )); output = 2*(num ./ den - 1) * freq; %subtract by one to remove DC, scale by freq to normalize output = output * maxN/N; %OPTIONAL: rescale to give louder output as rolloff increases function out = phasor(normfreq, samples); out = mod( (0:samples-1)*normfreq , 1); out = out * 2 - 1; %make bipolar ```

• Date: 2003-04-03 15:05:42
• By: David Lowenfels
```oops, there's an error in this version. frequency should not be normalized until after the maxN calculation is done.
```

### Direct pink noise synthesis with auto-correlated generator¶

notes
```Canonical C++ class with minimum system dependencies, BUT you must provide your own
uniform random number generator. Accurate range is a little over 9 octaves, degrading
gracefully beyond this. Estimated deviations +-0.25 dB from ideal 1/f curve in range.
Scaled to fit signed 16-bit range.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152``` ```// Pink noise class using the autocorrelated generator method. // Method proposed and described by Larry Trammell "the RidgeRat" -- // see http://home.earthlink.net/~ltrammell/tech/newpink.htm // There are no restrictions. // // ------------------------------------------------------------------ // // This is a canonical, 16-bit fixed-point implementation of the // generator in 32-bit arithmetic. There are only a few system // dependencies. // // -- access to an allocator 'malloc' for operator new // -- access to definition of 'size_t' // -- assumes 32-bit two's complement arithmetic // -- assumes long int is 32 bits, short int is 16 bits // -- assumes that signed right shift propagates the sign bit // // It needs a separate URand class to provide uniform 16-bit random // numbers on interval [1,65535]. The assumed class must provide // methods to query and set the current seed value, establish a // scrambled initial seed value, and evaluate uniform random values. // // // ----------- header ----------------------------------------------- // pinkgen.h #ifndef _pinkgen_h_ #define _pinkgen_h_ 1 #include #include // You must provide the uniform random generator class. #ifndef _URand_h_ #include "URand.h" #endif class PinkNoise { private: // Coefficients (fixed) static long int const pA; static short int const pPSUM; // Internal pink generator state long int contrib; // stage contributions long int accum; // combined generators void internal_clear( ); // Include a UNoise component URand ugen; public: PinkNoise( ); PinkNoise( PinkNoise & ); ~PinkNoise( ); void * operator new( size_t ); void pinkclear( ); short int pinkrand( ); } ; #endif // ----------- implementation --------------------------------------- // pinkgen.cpp #include "pinkgen.h" // Static class data long int const PinkNoise::pA = { 14055, 12759, 10733, 12273, 15716 }; short int const PinkNoise::pPSUM = { 22347, 27917, 29523, 29942, 30007 }; // Clear generator to a zero state. void PinkNoise::pinkclear( ) { int i; for (i=0; i<5; ++i) { contrib[i]=0L; } accum = 0L; } // PRIVATE, clear generator and also scramble the internal // uniform generator seed. void PinkNoise::internal_clear( ) { pinkclear(); ugen.seed(0); // Randomizes the seed! } // Constructor. Guarantee that initial state is cleared // and uniform generator scrambled. PinkNoise::PinkNoise( ) { internal_clear(); } // Copy constructor. Preserve generator state from the source // object, including the uniform generator seed. PinkNoise::PinkNoise( PinkNoise & Source ) { int i; for (i=0; i<5; ++i) contrib[i]=Source.contrib[i]; accum = Source.accum; ugen.seed( Source.ugen.seed( ) ); } // Operator new. Just fetch required object storage. void * PinkNoise::operator new( size_t size ) { return malloc(size); } // Destructor. No special action required. PinkNoise::~PinkNoise( ) { /* NIL */ } // Coding artifact for convenience #define UPDATE_CONTRIB(n) \ { \ accum -= contrib[n]; \ contrib[n] = (long)randv * pA[n]; \ accum += contrib[n]; \ break; \ } // Evaluate next randomized 'pink' number with uniform CPU loading. short int PinkNoise::pinkrand( ) { short int randu = ugen.urand() & 0x7fff; // U[0,32767] short int randv = (short int) ugen.urand(); // U[-32768,32767] // Structured block, at most one update is performed while (1) { if (randu < pPSUM) UPDATE_CONTRIB(0); if (randu < pPSUM) UPDATE_CONTRIB(1); if (randu < pPSUM) UPDATE_CONTRIB(2); if (randu < pPSUM) UPDATE_CONTRIB(3); if (randu < pPSUM) UPDATE_CONTRIB(4); break; } return (short int) (accum >> 16); } // ----------- application ------------------------------- short int pink_signal; void example(void) { PinkNoise pinkgen; int i; for (i=0; i<1024; ++i) pink_signal[i] = pinkgen.pinkrand(); } ```

### Discrete Summation Formula (DSF)¶

• Author or source: Stylson, Smith and others… (Alexander Kritov)
• Created: 2002-02-10 12:43:30
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14``` ```double DSF (double x, // input double a, // a<1.0 double N, // N

• Date: 2002-11-08 11:21:19
• By: dfl[AT]ccrma.stanford.edu
```According to Stilson + Smith, this should be

double s1 = pow(a,N+1.0)*sin((N-1.0)*x+fi);
^
!

Could be a typo though?
```
• Date: 2003-03-14 17:01:46
• By: Alex
```yepp..
```
• Date: 2003-03-20 04:20:51
• By: TT
```So what is wrong about "double" up there?
For DSF, do we have to update the phase (fi input) at every sample?
Another question is what's the input x supposed to represent? Thanks!
```
• Date: 2003-04-01 01:45:47
• By: David Lowenfels
```input x should be the phase, and fi is the initial phase I guess? Seems redundant to me.
There is nothing wrong with the double, there is a sign typo in the original posting.
```
```I'm not so sure that there is a sign typo. (I know--I'm five years late to this party.)

The author of this code just seems to have an off-by-one definition of N. If you expand it all out, it looks like Stilson & Smith's paper, except you have N here where S&S had N+1, and you have N-1 where S&S had N.

I think the code is equivalent. You just have to understand how to choose N to avoid aliasing.

I don't have it working yet, but that's how it looks to me as I prepare a DSF oscillator. More later.
```
• Date: 2008-11-02 11:47:07
• By: mysterious T
```Got it working nicely, but it took a few minutes to pluck it apart. Had to correct it for my pitch scheme, too. But it's quite amazing! Funny concept, though, it's like a generator with a built in filter. It holds up into very high pitches, too, in terms of aliasing, as far as I can tell... ehm...and without any further oversampling (so far).

Really, really nice! I was looking for a way to give my sinus an edge! ;)
```

### Drift generator¶

notes
```I use this drift to modulate any sound parameter of my synth.
It is very effective if it slightly modulates amplitude or frequency of an FM modulator.
It is based on an incremental random variable, sine-warped.
I like it because it is "continuous" (as opposite to "sample and hold"), and I can set
variation rate and max variation.
It can go to upper or lower constraint (+/- max drift) but it gradually decreases rate of
variation when approaching to the limit.
I use it exactly as an LFO (-1.f .. +1.f)
I use a table for sin instead of sin() function because this way I can change random
distribution, by selecting a different curve (different table) from sine...

I hope that it is clear ... (sigh... :-)
Bye!!!
P.S. Thank you for help in previous submission ;-)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15``` ```const int kSamples //Number of samples in fSinTable below float fSinTable[kSamples] // Tabulated sin() [0 - 2pi[ amplitude [-1.f .. 1.f] float fWhere// Index float fRate // Max rate of variation float fLimit //max or min value float fDrift // Output //I assume that random() is a number from 0.f to 1.f, otherwise scale it fWhere += fRate * random() //I update this drift in a long-term cycle, so I don't care of branches if (fWhere >= 1.f) fWhere -= 1.f else if (fWhere < 0.f) sWhere += 1.f fDrift = fLimit * fSinTable[(long) (fWhere * kSamples)] ```

```...sorry...
random() must be in [-1..+1] !!!
```

### Easy noise generation¶

notes
```Easy noise generation,
in .hpp,
b_noise = 19.1919191919191919191919191919191919191919;

alternatively, the number 19 below can be replaced with a number of your choice, to get
that particular flavour of noise.

Regards,
Ove Karlsen.
```
code
 ```1 2 3 4 5 6 7 8``` ``` b_noise = b_noise * b_noise; int i_noise = b_noise; b_noise = b_noise - i_noise; double b_noiseout; b_noiseout = b_noise - 0.5; b_noise = b_noise + 19; ```

```This is quite a good PRNG! The numbers it generates exhibit a slight a pattern (obviously, since it's not very sophisticated) but they seem quite usable! The real FFT spectrum is very flat and "white" with just one or two aberrant spikes while the imaginary spectrum is almost perfect (as is the case with most PRNGs). Very nice! Either that or I need more practice with MuPad...
```
```Alternatively you can do:

double b_noiselast = b_noise;
b_noise = b_noise + 19;
b_noise = b_noise * b_noise;
b_noise = b_noise + ((-b_noise + b_noiselast) * 0.5);
int i_noise = b_noise;
b_noise = b_noise - i_noise;

This will remove the patterning.
```
```>>b_noise = b_noise + ((-b_noise + b_noiselast) * 0.5);

That seems to reduce to just:

b_noise=(b_noise+b_noiselast) * 0.5;
```
```Hi, is this integer? Please do not disturb the forum, rather send me an email.

B.i.T
```
```The line is written like that, so you can change 0.5, to for instance 0.19.
```
```>>The line is written like that, so you can change 0.5, to for instance 0.19.

OK. Why would I do that? What's that number control?
```
```It controls the patterning. I usually write my algorithms tweakable.

You could try even lower aswell, maybe 1e-19.
```

### Fast Exponential Envelope Generator¶

• Author or source: Christian Schoenebeck
• Created: 2005-03-03 14:44:11
notes
```The naive way to implement this would be to use a exp() call for each point
of the envelope. Unfortunately exp() is quite a heavy function for most
CPUs, so here is a numerical, much faster way to compute an exponential
envelope (performance gain measured in benchmark: about factor 100 with a
Intel P4, gcc -O3 --fast-math -march=i686 -mcpu=i686).

Note: you can't use a value of 0.0 for levelEnd. Instead you have to use an
appropriate, very small value (e.g. 0.001 should be sufficiently small
enough).
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17``` ```const float sampleRate = 44100; float coeff; float currentLevel; void init(float levelBegin, float levelEnd, float releaseTime) { currentLevel = levelBegin; coeff = (log(levelEnd) - log(levelBegin)) / (releaseTime * sampleRate); } inline void calculateEnvelope(int samplePoints) { for (int i = 0; i < samplePoints; i++) { currentLevel += coeff * currentLevel; // do something with 'currentLevel' here ... } } ```

```is there a typo in the runtime equation? or am i missing something in the implementation?
```
• Date: 2005-03-04 15:13:56
• By: schoenebeck (@) software ( minus ) engineering.org
```Why should there be a typo?

Here is my benchmark code btw:
http://stud.fh-heilbronn.de/~cschoene/studienarbeit/benchmarks/exp.cpp
```
```ok, i think i get it. this can only work on blocks of samples, right? not per-sample calc?

i was confused because i could not find the input sample(s) in the runtime code. but now i see that the equation does not take an input; it merely generates a defined envelope accross the number of samples. my bad.
```
• Date: 2005-03-04 19:16:29
• By: schoenebeck (@) software ( minus ) engineering.org
```Well, the code above is only meant to show the principle. Of course you
would adjust it for your application. The question if you are calculating
on a per-sample basis or applying the envelope to a block of samples
within a tight loop doesn't really matter; it would just mean an
adjustment of the interface of the execution code, which is trivial.
```
```This is not working for long envelopes because of numerical accury problems. Try calculating is over 10 seconds @ 192KHz to see what I mean: it drifts.
I have an equivalent system that permits to have linear to log and to exp curves with a simple parameter. I may submit it one of these days...

Sebastien Metrot
--
http://www.usbsounds.com
```
• Date: 2005-03-05 13:48:12
• By: schoenebeck (@) software ( minus ) engineering.org
```No, here is a test app which shows the introduced drift:
http://stud.fh-heilbronn.de/~cschoene/studienarbeit/benchmarks/expaccuracy.cpp

Even with an envelope duration of 30s, which is really quite long, a sample
rate of 192kHz and single-precision floating point calculation I get this
result:

Calculated sample points: 5764846
Demanded duration: 30.000000 s
Actual duration: 30.025240 s

So the envelope just drifts about 25ms for that long envelope!
```
```I believe you are seeing unrealistic results with this test because on x86 the fpu's internal format is 80bits and your compiler probably optimises this cases quite easily. Try doing the same test, calculating the same envelope, but by breaking the calculation in blocks of 256 or 512 samples at a time and then storing in memory the temp values for the next block. In this case you may see diferent results and a much bigger drift (that's my experience with the same algo).
Anyway my algo is a bit diferent as it permits to change the curent type with a parameter, this makes the formula looks like
value = value * coef + contant;
May be this leads to more calculation errors :).
```
• Date: 2005-03-09 13:33:43
• By: schoenebeck (@) software ( minus ) engineering.org
```And again... no! :)

Replace the C equation by:

asm volatile (
"mulss %%xmm1,%%xmm0  # coeff *= currentLevel\n\t"
"addss %%xmm0,%%xmm1  # currentLevel += coeff * currentLevel\n\t"
"movss %%xmm1,%0      # store currentLevel\n\t"
: "=m" (currentLevel) /* %0 */
: "m" (coeff),        /* %1 */
"m" (currentLevel)  /* %2 */
);

This is a SSE1 assembly implementation. The SSE registers are only 32 bit
large by guarantee. And this is the result I get:

Calculated sample points: 5764845
Demanded duration: 30.000000 s
Actual duration: 30.025234 s

So this result differs just 1 sample point from the x86 FPU solution! So
believe me, this numerical solution is safe!

(Of course the assembly code above is NOT meant as optimization, it's just
to demonstrate the accuracy even for 32 bit / single precision FP
calculation)
```
• Date: 2005-03-23 22:42:06
• By: m (at) mindplay (dot) dk
```in my tests, the following code produced the exact same results, and saves one operation (the addition) per sample - so it should be faster:

const float sampleRate = 44100;
float coeff;
float currentLevel;

void init(float levelBegin, float levelEnd, float releaseTime) {
currentLevel = levelBegin;
coeff = exp(log(levelEnd)) /
(releaseTime * sampleRate);
}

inline void calculateEnvelope(int samplePoints) {
for (int i = 0; i < samplePoints; i++) {
currentLevel *= coeff;
// do something with 'currentLevel' here
...
}
}

...

Also, assuming that your startLevel is 1.0, to calculate an appropriate endLevel, you can use something like:

endLevel = 10 ^ dB/20;

where dB is your endLevel in decibels (and must be a negative value of course) - for amplitude envelopes, -90 dB should be a suitable level for "near inaudible"...
```
• Date: 2005-03-31 14:45:51
• By: schoenebeck (@) software ( minus ) engineering.org
```              Sorry, you are right of course; that simplification of the execution
equation works here because we are calculating all points with linear
discretization. But you will agree that your init() function is not good,
because exp(log(x)) == x and it's not generalized at all. Usually you might
have more than one exp segment in your EG and maybe even have an exp attack
segment. So we arrive at the following solution:

const float sampleRate = 44100;
float coeff;
float currentLevel;

void init(float levelBegin, float levelEnd, float releaseTime) {
currentLevel = levelBegin;
coeff = 1.0f + (log(levelEnd) - log(levelBegin)) /
(releaseTime * sampleRate);
}

inline void calculateEnvelope(int samplePoints) {
for (int i = 0; i < samplePoints; i++) {
currentLevel *= coeff;
// do something with 'currentLevel' here
...
}
}

You can use a dB conversion for both startLevel and endLevel of course.
```
• Date: 2006-03-10 01:53:44
• By: na
```i would say that calculation of coeff is still wrong. It should be :
coeff = pow( levelEnd / levelBegin, 1 / N );
```
• Date: 2006-03-10 02:23:29
• By: na[ eldar # starman # ee]
```or coeff = exp(log(levelEnd/levelBegin) /
(releaseTime * sampleRate) );
not sure but it looks computationally more expensive
```
• Date: 2006-11-26 15:44:04
• By: hc.xmg@.i.l.e
```what's about?
coeff = 1.0f + (log(levelEnd) - log(levelBegin)) /
(releaseTime * sampleRate - 1);
```
• Date: 2006-11-26 15:55:12
• By: hc.xmg@.i.l.e
```sorry for the double post. and i'm now almost sure, that it should be:
coeff = 1.0f + (log(levelEnd) - log(levelBegin)) /
(releaseTime * sampleRate + 1);
```

### Fast LFO in Delphi…¶

notes
```[from Didier's mail...]
[see attached zip file too!]

I was working on a flanger, & needed an LFO for it. I first used a Sin(), but it was too
slow, then tried a big wavetable, but it wasn't accurate enough.
I then checked the alternate sine generators from your web site, & while they're good,
they all can drift, so you're also wasting too much CPU in branching for the drift checks.
So I made a quick & easy linear LFO, then a sine-like version of it. Can be useful for
LFO's, not to output as sound.
If has no branching & is rather simple. 2 Abs() but apparently they're fast. In all cases
faster than a Sin()
It's in delphi, but if you understand it you can translate it if you want.
It uses a 32bit integer counter that overflows, & a power for the sine output.
If you don't know delphi, \$ is for hex (h at the end in c++?), Single is 32bit float,
integer is 32bit integer (signed, normally).
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75``` ```unit Unit1; interface uses Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs, StdCtrls, ExtCtrls, ComCtrls; type TForm1 = class(TForm) PaintBox1: TPaintBox; Bevel1: TBevel; procedure PaintBox1Paint(Sender: TObject); private { Private declarations } public { Public declarations } end; var Form1: TForm1; implementation {\$R *.DFM} procedure TForm1.PaintBox1Paint(Sender: TObject); var n,Pos,Speed:Integer; Output,Scale,HalfScale,PosMul:Single; OurSpeed,OurScale:Single; begin OurSpeed:=100; // 100 samples per cycle OurScale:=100; // output in -100..100 Pos:=0; // position in our linear LFO Speed:=Round(\$100000000/OurSpeed); // --- triangle LFO --- Scale:=OurScale*2; PosMul:=Scale/\$80000000; // loop for n:=0 to 299 do Begin // inc our 32bit integer LFO pos & let it overflow. It will be seen as signed when read by the math unit Pos:=Pos+Speed; Output:=Abs(Pos*PosMul)-OurScale; // visual Paintbox1.Canvas.Pixels[n,Round(100+Output)]:=clRed; End; // --- sine-like LFO --- Scale:=Sqrt(OurScale*4); PosMul:=Scale/\$80000000; HalfScale:=Scale/2; // loop for n:=0 to 299 do Begin // inc our 32bit integer LFO pos & let it overflow. It will be seen as signed when read by the math unit Pos:=Pos+Speed; Output:=Abs(Pos*PosMul)-HalfScale; Output:=Output*(Scale-Abs(Output)); // visual Paintbox1.Canvas.Pixels[n,Round(100+Output)]:=clBlue; End; end; end. ```

```LFO Class...

unit Unit1;

interface

uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls, ExtCtrls, ComCtrls;

type
TLFOType = (lfoTriangle,lfoSine);
TLFO = class(TObject)
private
iSpeed     : Integer;
fSpeed     : Single;
fMax, fMin : Single;
fValue     : Single;
fPos       : Integer;
fType      : TLFOType;
fScale     : Single;
fPosMul    : Single;
fHalfScale : Single;
function GetValue:Single;
procedure SetType(tt: TLFOType);
procedure SetMin(v:Single);
procedure SetMax(v:Single);
public
{ Public declarations }
constructor Create;
published
property Speed:Single read FSpeed Write FSpeed;
property Min:Single read FMin write SetMin;
property Max:Single read FMax Write SetMax;
property LFO:TLFOType read fType Write SetType;
end;

TForm1 = class(TForm)
Bevel1: TBevel;
PaintBox1: TPaintBox;
procedure PaintBox1Paint(Sender: TObject);
private
{ Private declarations }
public
{ Public declarations }
end;

var
Form1: TForm1;

implementation

{\$R *.DFM}

constructor TLFO.Create;
begin
fSpeed:=100;
fMax:=1;
fMin:=0;
fValue:=0;
fPos:=0;
iSpeed:=Round(\$100000000/fSpeed);
fType:=lfoTriangle;
fScale:=fMax-((fMin+fMax)/2);
end;

procedure TLFO.SetType(tt: TLFOType);
begin
fType:=tt;
if fType = lfoSine then
begin
fPosMul:=(Sqrt(fScale*2))/\$80000000;
fHalfScale:=(Sqrt(fScale*2))/2;
end
else
begin
fPosMul:=fScale/\$80000000;
end;
end;

procedure TLFO.SetMin(v: Single);
begin
fMin:=v;
fScale:=fMax-((fMin+fMax)/2);
end;

procedure TLFO.SetMax(v: Single);
begin
fMax:=v;
fScale:=fMax-((fMin+fMax)/2);
end;

function TLFO.GetValue:Single;
begin
if fType = lfoSine then
begin
Result:=Abs(fPos*fPosMul)-fHalfScale;
Result:=Result*(fHalfScale*2-Abs(Result))*2;
Result:=Result+((fMin+fMax)/2);
end
else
begin
Result:=Abs(fPos*(2*fPosMul))+fMin;
end;
fPos:=fPos+iSpeed;
end;

procedure TForm1.PaintBox1Paint(Sender: TObject);
var n    : Integer;
LFO1 : TLFO;
begin

LFO1:=TLFO.Create;
LFO1.Min:=-100;
LFO1.Max:=100;
LFO1.Speed:=100;
LFO1.LFO:=lfoTriangle;

// --- triangle LFO ---
for n:=0 to 299 do Paintbox1.Canvas.Pixels[n,Round(100+LFO1.Value)]:=clRed;

LFO1.LFO:=lfoSine;
// --- sine-like LFO ---
for n:=0 to 299 do Paintbox1.Canvas.Pixels[n,Round(100+LFO1.Value)]:=clBlue;
end;

end.
```
```Ups, i forgot something:

TLFO = class(TObject)
private
...
procedure SetSpeed(v:Single);
public
...
published
...
property Speed:Single read FSpeed Write SetSpeed;
...
end;

...

constructor TLFO.Create;
begin
...
Speed:=100;
...
end;

procedure TLFO.SetSpeed(v:Single);
begin
fSpeed:=v;
iSpeed:=Round(\$100000000/fSpeed);
end;

...
```

### Fast SIN approximation for usage in e.g. additive synthesizers¶

• Author or source: neotec
• Created: 2008-12-09 11:56:33
notes
```This code presents 2 'fastsin' functions. fastsin2 is less accurate than fastsin. In fact
it's a simple taylor series, but optimized for integer phase.

phase is in 0 -> (2^32)-1 range and maps to 0 -> ~2PI

I get about 55000000 fastsin's per second on my P4,3.2GHz which would give a nice Kawai K5
emulation using 64 harmonics and 8->16 voices.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32``` ```float fastsin(UINT32 phase) { const float frf3 = -1.0f / 6.0f; const float frf5 = 1.0f / 120.0f; const float frf7 = -1.0f / 5040.0f; const float frf9 = 1.0f / 362880.0f; const float f0pi5 = 1.570796327f; float x, x2, asin; UINT32 tmp = 0x3f800000 | (phase >> 7); if (phase & 0x40000000) tmp ^= 0x007fffff; x = (*((float*)&tmp) - 1.0f) * f0pi5; x2 = x * x; asin = ((((frf9 * x2 + frf7) * x2 + frf5) * x2 + frf3) * x2 + 1.0f) * x; return (phase & 0x80000000) ? -asin : asin; } float fastsin2(UINT32 phase) { const float frf3 = -1.0f / 6.0f; const float frf5 = 1.0f / 120.0f; const float frf7 = -1.0f / 5040.0f; const float f0pi5 = 1.570796327f; float x, x2, asin; UINT32 tmp = 0x3f800000 | (phase >> 7); if (phase & 0x40000000) tmp ^= 0x007fffff; x = (*((float*)&tmp) - 1.0f) * f0pi5; x2 = x * x; asin = (((frf7 * x2 + frf5) * x2 + frf3) * x2 + 1.0f) * x; return (phase & 0x80000000) ? -asin : asin; } ```

• Date: 2008-12-09 12:11:11
• By: neotec
```PS: To use this as an OSC you'll need the following vars/equ's:

UINT32 phase = 0;
UINT32 step = frequency * powf(2.0f, 32.0f) / samplerate;

Then it's just:
...
out = fastsin(phase);
phase += step;
...
```
```Woah! Seven multiplies, on top of those adds and memory lookup. Is this really all that fast?
```

### Fast Whitenoise Generator¶

notes
```This is Whitenoise... :o)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18``` ```float g_fScale = 2.0f / 0xffffffff; int g_x1 = 0x67452301; int g_x2 = 0xefcdab89; void whitenoise( float* _fpDstBuffer, // Pointer to buffer unsigned int _uiBufferSize, // Size of buffer float _fLevel ) // Noiselevel (0.0 ... 1.0) { _fLevel *= g_fScale; while( _uiBufferSize-- ) { g_x1 ^= g_x2; *_fpDstBuffer++ = g_x2 * _fLevel; g_x2 += g_x1; } } ```

```Works well! Kinda fast! The spectrum looks completely flat in an FFT analyzer.
```
```As I said! :-)
Take care
```
```I'm now waiting for pink and brown. :-)
```
```To get pink noise, you can apply a 3dB/Oct filter, for example the pink noise filter in the Filters section.

To get brown noise, apply an one pole LP filter to get a 6dB/oct slope.

Peter
```
```Yeah, I know how to do it with a filter. I was just looking to see if this guy had anything else clever up his sleeve.

I'm currently using this great stuff:

vellocet.com/dsp/noise/VRand.html
```
```I compiled it, but I get some grainyness that a unisgned long LC algorithm does not give me...  am I the only one?

pa
```
```Did you do everything right? It works here.
```
```I've noticed that my code is similar to a so called "feedback shift register" as used in the Commodore C64 Soundchip 6581 called SID for noise generation.

en.wikipedia.org/wiki/Linear_feedback_shift_register
en.wikipedia.org/wiki/MOS_Technology_SID
www.cc65.org/mailarchive/2003-06/3156.html
```
```SID noise! cool.
```

### Fast sine and cosine calculation¶

• Author or source: Lot’s or references… Check Julius O. SMith mainly
• Type: waveform generation
• Created: 2002-01-17 00:54:32
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13``` ```init: float a = 2.f*(float)sin(Pi*frequency/samplerate); float s; s = 0.5f; s = 0.f; loop: s = s - a*s; s = s + a*s; output_sine = s; output_cosine = s ```

• Date: 2003-04-05 10:52:49
• By: DFL
```Yeah, this is a cool trick! :)

FYI you can set s to whatever amplitude of sinewave you desire. With 0.5, you will get +/- 0.5
```
```After a while it may drift, so you should resync it as follows:
const float tmp=1.5f-0.5f*(s*s+s*s);
s*=tmp; s*=tmp;

This assumes you set s to 1.0 initially.

'Tick
```
• Date: 2003-04-08 09:19:40
• By: DFL
```Just to expalin the above "resync" equation
(3-x)/2 is an approximation of 1/sqrt(x)
So the above is actually renormalizing the complex magnitude.
[ sin^2 (x) + cos^2(x) = 1 ]
```
• Date: 2003-05-15 08:26:22
• By: nigel
```This is the Chamberlin state variable filter specialized for infinite Q oscillation. A few things to note:

Like the state variable filter, the upper frequency limit for stability is around one-sixth the sample rate.

The waveform symmetry is very pure at low frequencies, but gets skewed as you get near the upper limit.

For low frequencies, sin(n) is very close to n, so the calculation for "a" can be reduced to a = 2*Pi*frequency/samplerate.

You shouldn't need to resync the oscillator--for fixed point and IEEE floating point, errors cancel exactly, so the osciallator runs forever without drifting in amplitude or frequency.
```
```I made a nice little console 'game' using your cordic sinewave approximation. Download it at http://users.pandora.be/antipro/Other/Ascillator.zip (includes source code). Just for oldschool fun :).
```
• Date: 2004-12-22 16:52:20
• By: hplus
```Note that the peaks of the waveforms will actually be between samples, and the functions will be phase offset by one half sample's worth. If you need exact phase, you can compensate by interpolating using cubic hermite interpolation.
```
• Date: 2007-07-24 20:33:12
• By: more on that topic…
```... can be found in Jon Datorro, Effect Design, Part 3, a paper that can be easily found in the web.

Funny, this is just a complex multiply that is optimized for small angles (low frequencies)

When the CPU rounding mode is set to nearest, it should be stable, at least for small frequencies.
```
```More on that can be found in Jon Datorro, Effect Design, Part 3, a paper that can be easily found in the web.

Funny, this is just a complex multiply that is optimized for small angles (low frequencies)

When the CPU rounding mode is set to nearest, it should be stable, at least for small frequencies.
```
```How do I set a particular phase for this? I've tried setting s = cos(phase) and s = sin(phase), but that didn't seem to be accurate enough.

Thanks
```

### Fast sine wave calculation¶

• Author or source: James McCartney in Computer Music Journal, also the Julius O. Smith paper
• Type: waveform generation
• Created: 2002-01-17 00:52:33
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22``` ```Variables: ip = phase of the first output sample in radians w = freq*pi / samplerate b1 = 2.0 * cos(w) Init: y1=sin(ip-w) y2=sin(ip-2*w) Loop: y0 = b1*y1 - y2 y2 = y1 y1 = y0 output is in y0 (y0 = sin(ip + n*freq*pi / samplerate), n= 0, 1, 2, ... I *think*) Later note by James McCartney: if you unroll such a loop by 3 you can even eliminate the assigns!! y0 = b1*y1 - y2 y2 = b1*y0 - y1 y1 = b1*y2 - y0 ```

```try using this to make sine waves with frequency less that 1. I did and it gives very rough half triangle-like waves. Is there any way to fix this? I want to use a sine generated for LFO so I need one that works for low frequencies.
```
```looks like the formula has gotten munged.
w = freq * twopi / samplerate
```

### Fast square wave generator¶

• Author or source: Wolfgang (ed.tfoxen@redienhcsw)
• Type: NON-bandlimited osc…
• Created: 2002-02-10 12:46:22
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15``` ```Idea: use integer overflow to avoid conditional jumps. // init: typedef unsigned long ui32; float sampleRate = 44100.0f; // whatever float freq = 440.0f; // 440 Hz float one = 1.0f; ui32 intOver = 0L; ui32 intIncr = (ui32)(4294967296.0 / hostSampleRate / freq)); // loop: (*((ui32 *)&one)) &= 0x7FFFFFFF; // mask out sign bit (*((ui32 *)&one)) |= (intOver & 0x80000000); intOver += intIncr; ```

```So, how would I get the output into a float variable like square_out, for instance?
```
```In response to lancej, yo can declare a union with a float and a int and operate the floatas as here using the int part of the union.

If I remeber correctly the value for -1.f = 0xBF800000 and the value for 1.f = 0x3F800000, note the 0x80000000 difference between them that is the sign.
```

### Gaussian White Noise¶

notes
```SOURCE:

Steven W. Smith:
The Scientist and Engineer's Guide to Digital Signal Processing
http://www.dspguide.com
```
code
 ```1 2 3 4 5 6``` ```#define PI 3.1415926536f float R1 = (float) rand() / (float) RAND_MAX; float R2 = (float) rand() / (float) RAND_MAX; float X = (float) sqrt( -2.0f * log( R1 )) * cos( 2.0f * PI * R2 ); ```

```The previous one seems better for me, since it requires only a rand, half log and half sqrt per sample.
Actually, I used that one, but I can't remember where I found it, too. Maybe on Knuth's book.
```

### Gaussian White noise¶

• Author or source: Alexey Menshikov
• Created: 2002-08-01 01:13:47
notes
```Code I use sometimes, but don't remember where I ripped it from.

- Alexey Menshikov
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26``` ```#define ranf() ((float) rand() / (float) RAND_MAX) float ranfGauss (int m, float s) { static int pass = 0; static float y2; float x1, x2, w, y1; if (pass) { y1 = y2; } else { do { x1 = 2.0f * ranf () - 1.0f; x2 = 2.0f * ranf () - 1.0f; w = x1 * x1 + x2 * x2; } while (w >= 1.0f); w = (float)sqrt (-2.0 * log (w) / w); y1 = x1 * w; y2 = x2 * w; } pass = !pass; return ( (y1 * s + (float) m)); } ```

• Date: 2004-01-29 15:41:35
• By: davidchristenATgmxDOTnet
```White Noise does !not! consist of uniformly distributed values. Because in white noise, the power of the frequencies are uniformly distributed. The values must be normal (or gaussian) distributed. This is achieved by the Box-Muller Transformation. This function is the polar form of the Box-Muller Transformation. It is faster and numeriacally more stable than the basic form. The basic form is coded in the other (second) post.
Detailed information on this topic:
http://www.taygeta.com/random/gaussian.html
http://www.eece.unm.edu/faculty/bsanthan/EECE-541/white2.pdf

Cheers David
```
```I'm trying to implement this in C#, but y2 isn't initialized. Is this a typo?
```
```@nick: Way to late, but y2 will always be initialized as in the first run "pass" is 0 (i.e. false). The C# compiler just can't prove it.
```
```David is wrong. The distribution of the sample values is irrelevant. 'white' simply describes the spectrum. Any series of sequentially independent random values -- whatever their distribution -- will have a white spectrum.
```

### Generator¶

• Author or source: Paul Sernine
• Type: antialiased sawtooth
• Created: 2006-02-23 22:38:56
notes
```This code generates a swept antialiasing sawtooth in a raw 16bit pcm file.
It is based on the quad differentiation of a 5th order polynomial. The polynomial
harmonics (and aliased harmonics) decay at 5*6 dB per oct. The differenciators correct the
spectrum and waveform, while aliased harmonics are still attenuated.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47``` ```/* clair.c Examen Partiel 2b T.Rochebois 02/03/98 */ #include #include main() { double phase=0,dphase,freq,compensation; double aw0=0,aw1=0,ax0=0,ax1=0,ay0=0,ay1=0,az0=0,az1=0,sortie; short aout; int sr=44100; //sample rate (Hz) double f_debut=55.0;//start freq (Hz) double f_fin=sr/6.0;//end freq (Hz) double octaves=log(f_fin/f_debut)/log(2.0); double duree=50.0; //duration (s) int i; FILE* f; f=fopen("saw.pcm","wb"); for(i=0;i1.0) phase-=2.0; //phase wrapping (-1,+1) //polynomial calculation (extensive continuity at -1 +1) // 7 1 3 1 5 //P(x) = --- x - -- x + --- x // 360 36 120 aw0=phase*(7.0/360.0 + phase*phase*(-1/36.0 + phase*phase*(1/120.0))); // quad differentiation (first order high pass filters) ax0=aw1-aw0; ay0=ax1-ax0; az0=ay1-ay0; sortie=az1-az0; //compensation of the attenuation of the quad differentiator //this can be calculated at "control rate" and linearly //interpolated at sample rate. compensation=1.0/(dphase*dphase*dphase*dphase); // compensation and output aout=(short)(15000.0*compensation*sortie); fwrite(&aout,1,2,f); //old memories of differentiators aw1=aw0; ax1=ax0; ay1=ay0; az1=az0; } fclose(f); } ```

```More infos and discussions in the KVR thread:
http://www.kvraudio.com/forum/viewtopic.php?t=123498
```
```nice but i prefer the fishy algo, it generates less alias.
bonaveture rosignol
```

### Inverted parabolic envelope¶

• Author or source: James McCartney
• Type: envellope generation
• Created: 2002-01-17 00:57:43
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19``` ```dur = duration in samples midlevel = amplitude at midpoint beglevel = beginning and ending level (typically zero) amp = midlevel - beglevel; rdur = 1.0 / dur; rdur2 = rdur * rdur; level = beglevel; slope = 4.0 * amp * (rdur - rdur2); curve = -8.0 * amp * rdur2; ... for (i=0; i

```This parabola approximation seems more like a linear than a parab/expo envelope... or i'm mistaking something but i tryed everything and is only linear.
```
```slope is linear, but 'slope' is a function of 'curve'. If you imagine you threw a ball upwards, think of 'curve' as the gravity, 'slope' as the vertical velocity, and 'level' as the vertical displacement.
```
• Date: 2005-01-17 07:39:28
• By: asynth(at)io(dot)com
```This is not an approximation of a parabola, it IS a parabola.
This entry has become corrupted since it was first posted. Should be:

for (i=0; i<dur; ++i) {
out = level;
level += slope;
slope += curve;
}
```

### Matlab/octave code for minblep table generation¶

notes
```When I tested this code, it was running with each function in a separate file... so it
might need some tweaking (endfunction statements?) if you try and run it all as one file.

Enjoy!

PS There's a C++ version by Daniel Werner here.
http://www.experimentalscene.com/?type=2&id=1
Not sure if it the output is any different than my version.
(eg no thresholding in minphase calculation)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137``` ```% Octave/Matlab code to generate a minblep table for bandlimited synthesis %% original minblep technique described by Eli Brandt: %% http://www.cs.cmu.edu/~eli/L/icmc01/hardsync.html % (c) David Lowenfels 2004 % you may use this code freely to generate your tables, % but please send me a free copy of the software that you % make with it, or at least send me an email to say hello % and put my name in the software credits :) % (IIRC: mps and clipdb functions are from Julius Smith) % usage: % fc = dilation factor % Nzc = number of zero crossings % omega = oversampling factor % thresh = dB threshold for minimum phase calc mbtable = minblep( fc, Nzc, omega, thresh ); mblen = length( mbtable ); save -binary mbtable.mat mbtable ktable nzc mblen; ********************************************* function [out] = minblep( fc, Nzc, omega, thresh ) out = filter( 1, [1 -1], minblip( fc, Nzc, omega, thresh ) ); len = length( out ); normal = mean( out( floor(len*0.7):len ) ) out = out / normal; %% normalize %% now truncate so it ends at proper phase cycle for minimum discontinuity thresh = 1e-6; for i = len:-1:len-1000 % pause a = out(i) - thresh - 1; b = out(i-1) - thresh - 1; % i if( (abs(a) < thresh) & (a > b) ) break; endif endfor %out = out'; out = out(1:i); ********************************************* function [out] = minblip( fc, Nzc, omega, thresh ) if (nargin < 4 ) thresh = -100; end if (nargin < 3 ) omega = 64; end if (nargin < 2 ) Nzc = 16; end if (nargin < 1 ) fc = 0.9; end blip = sinctable( omega, Nzc, fc ); %% length(blip) must be nextpow2! (if fc < 1 ); mag = fft( blip ); out = real( ifft( mps( mag, thresh ) ) ); ********************************************* function [sm] = mps(s, thresh) % [sm] = mps(s) % create minimum-phase spectrum sm from complex spectrum s if (nargin < 2 ) thresh = -100; endif s = clipdb(s, thresh); sm = exp( fft( fold( ifft( log( s ))))); ********************************************* function [clipped] = clipdb(s,cutoff) % [clipped] = clipdb(s,cutoff) % Clip magnitude of s at its maximum + cutoff in dB. % Example: clip(s,-100) makes sure the minimum magnitude % of s is not more than 100dB below its maximum magnitude. % If s is zero, nothing is done. as = abs(s); mas = max(as(:)); if mas==0, return; end if cutoff >= 0, return; end thresh = mas*10^(cutoff/20); % db to linear toosmall = find(as < thresh); clipped = s; clipped(toosmall) = thresh; ********************************************* function [out, phase] = sinctable( omega, Nzc, fc ) if (nargin < 3 ) fc = 1.0 %% cutoff frequency end %if if (nargin < 2 ) Nzc = 16 %% number of zero crossings end %if if (nargin < 1 ) omega = 64 %% oversampling factor end %if Nzc = Nzc / fc %% This ensures more flatness at the ends. phase = linspace( -Nzc, Nzc, Nzc*omega*2 ); %sinc = sin( pi * fc * phase) ./ (pi * fc * phase); num = sin( pi*fc*phase ); den = pi*fc*phase; len = length( phase ); sinc = zeros( len, 1 ); %sinc = num ./ den; for i=1:len if ( den(i) ~= 0 ) sinc(i) = num(i) / den(i); else sinc(i) = 1; end end %for out = sinc; window = blackman( len ); out = out .* window; ```

notes
```Please see the full description of the algorithm with public domain c++ code here:
```
code
 ```1 2 3 4``` ```It's here: http://zynaddsubfx.sourceforge.net/doc/PADsynth/PADsynth.htm You may copy it (everything is public domain). Paul ```

```Impressed at first hearing! Well documented.
```
```1. Isn't this plain additive synthesis
2. Isn't this the algorithm used by the waldorf microwave synths?
```
```1. Nope. This is not a plain additive synthesis. It's a special kind :-P Read the doc again :)
2. No way.. this is NOT even close to waldord microwave synths. Google for this :)
```

### PRNG for non-uniformly distributed values from trigonometric identity¶

• Author or source: moc.liamg@321tiloen
• Type: pseudo-random number generator
• Created: 2009-09-02 07:16:14
notes
```a method, which generates random numbers in the [-1,+1] range, while having a probability
density function with less concentration of values near zero for sin().
you can use an approximation of sin() and/or experiment with such an equation for
different distributions. using tan() will accordingly invert the pdf graph i.e. more
concentration near zero, but the output range will be also affected.

http://www.stat.wisc.edu/~larget/math496/random2.html

regards
lubomir
```
code
 ```1 2 3 4 5``` ```//init x=y=1; //sampleloop y=sin((x+=1)*y); ```

### Parabolic shaper¶

notes
```This function can be used for oscillators or shaper.
it can be driven by a phase accumulator or an audio input.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```Function Parashape(inp:single):single; var fgh,tgh:single; begin fgh := inp ; fgh := 0.25-f_abs(fgh) ; tgh := fgh ; tgh := 1-2*f_abs(tgh); fgh := fgh*8; result := fgh*tgh ; end; // f_abs is the function of ddsputils unit. ```

### Phase modulation Vs. Frequency modulation¶

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>

```

### Phase modulation Vs. Frequency modulation II¶

• Author or source: James McCartney
• Created: 2003-12-01 08:24:26
notes
```The difference between FM & PM in a digital oscillator is that FM is added to the
frequency before the phase integration, while PM is added to the phase after the phase
integration. Phase integration is when the old phase for the oscillator is added to the
current frequency (in radians per sample) to get the new phase for the oscillator. The
equivalent PM modulator to obtain the same waveform as FM is the integral of the FM
modulator. Since the integral of sine waves are inverted cosine waves this is no problem.
In modulators with multiple partials, the equivalent PM modulator will have different
relative partial amplitudes. For example, the integral of a square wave is a triangle
wave; they have the same harmonic content, but the relative partial amplitudes are
different. These differences make no difference since we are not trying to exactly
recreate FM, but real (or nonreal) instruments.

The reason PM is better is because in PM and FM there can be non-zero energy produced at 0
Hz, which in FM will produce a shift in pitch if the FM wave is used again as a modulator,
however in PM the DC component will only produce a phase shift. Another reason PM is
better is that the modulation index (which determines the number of sidebands produced and
which in normal FM is calculated as the modulator amplitude divided by frequency of
modulator) is not dependant on the frequency of the modulator, it is always equal to the
amplitude of the modulator in radians. The benefit of solving the DC frequency shift
problem, is that cascaded carrier-modulator pairs and feedback modulation are possible.
The simpler calculation of modulation index makes it easier to have voices keep the same
harmonic structure throughout all pitches.

The basic mathematics of phase modulation are available in any text on electronic
communication theory.

Below is some C code for a digital oscillator that implements FM,PM,and AM. It illustrates
the difference in implementation of FM & PM. It is only meant as an example, and not as an
efficient implementation.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38``` ```/* Example implementation of digital oscillator with FM, PM, & AM */ #define PI 3.14159265358979 #define RADIANS_TO_INDEX (512.0 / (2.0 * PI)) typedef struct{ /* oscillator data */ double freq; /* oscillator frequency in radians per sample */ double phase; /* accumulated oscillator phase in radians */ double wavetable; /* waveform lookup table */ } OscilRec; /* oscil - compute 1 sample of oscillator output whose freq. phase and * wavetable are in the OscilRec structure pointed to by orec. */ double oscil(orec, fm, pm, am) OscilRec *orec; /* pointer to the oscil's data */ double fm; /* frequency modulation input in radians per sample */ double pm; /* phase modulation input in radians */ double am; /* amplitude modulation input in any units you want */ { long tableindex; /* index into wavetable */ double instantaneous_freq; /* oscillator freq + freq modulation */ double instantaneous_phase; /* oscillator phase + phase modulation */ double output; /* oscillator output */ instantaneous_freq = orec->freq + fm; /* get instantaneous freq */ orec->phase += instantaneous_freq; /* accumulate phase */ instantaneous_phase = orec->phase + pm; /* get instantaneous phase */ /* convert to lookup table index */ tableindex = RADIANS_TO_INDEX * instantaneous_phase; tableindex &= 511; /* make it mod 512 === eliminate multiples of 2*k*PI */ output = orec->wavetable[tableindex] * am; /* lookup and mult by am input */ return (output); /* return oscillator output */ } ```

```As the PM/FM is "iterative", won't this code produce different results at different sampling rates? How can this be prevented?

```

### Pseudo-Random generator¶

• Author or source: Hal Chamberlain, “Musical Applications of Microprocessors” (Phil Burk)
• Type: Linear Congruential, 32bit
• Created: 2002-01-17 03:11:50
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
 ```1 2 3 4 5 6 7 8``` ```/* 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; } ```

notes
```This is written in Actionscript 3.0 (Flash9). You can listen to the example at
http://lab.andre-michelle.com/playing-with-pulse-harmonics
It allows to morph between a sinus like quadratic function and an ordinary pulse width
with adjustable pulse width. Note that the slope while morphing is always zero at the edge
points of the waveform. It is not just distorsion.
```
code
 `1` ```http://lab.andre-michelle.com/swf/f9/pulsequad/PulseQuad.as ```

### Pulsewidth modulation¶

• Author or source: Steffan Diedrichsen
• Type: waveform generation
• Created: 2002-01-15 21:29:52
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%.
```

### Quick & Dirty Sine¶

• Author or source: MisterToast
• Type: Sine Wave Synthesis
• Created: 2007-01-08 10:50:10
notes
```This is proof of concept only (but code works--I have it in my synth now).

Note that x must come in as 0<x<=4096. If you want to scale it to something else (like
0<x<=2*M_PI), do it in the call. Or do the math to scale the constants properly.

There's not much noise in here. A few little peaks here and there. When the signal is at
-20dB, the worst noise is at around -90dB.

For speed, you can go all floats without much difference. You can get rid of that unitary
negate pretty easily, as well. A couple other tricks can speed it up further--I went for
clarity in the code.

The result comes out a bit shy of the range -1<x<1. That is, the peak is something like
0.999.

Where did this come from? I'm experimenting with getting rid of my waveform tables, which
require huge amounts of memory. Once I had the Hamming anti-ringing code in, it looked
like all my waveforms were smooth enough to approximate with curves. So I started with
sine. Pulled my table data into Excel and then threw the data into a curve-fitting
application.

This would be fine for a synth. The noise is low enough that you could easily get away
with it. Ideal for a low-memory situation. My final code will be a bit harder to
understand, as I'll break the curve up and curve-fit smaller sections.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24``` ```float xSin(double x) { //x is scaled 0<=x<4096 const double A=-0.015959964859; const double B=217.68468676; const double C=0.000028716332164; const double D=-0.0030591066066; const double E=-7.3316892871734489e-005; double y; bool negate=false; if (x>2048) { negate=true; x-=2048; } if (x>1024) x=2048-x; if (negate) y=-((A+x)/(B+C*x*x)+D*x-E); else y=(A+x)/(B+C*x*x)+D*x-E; return (float)y; } ```

```Improved version:

float xSin(double x)
{
//x is scaled 0<=x<4096
const double A=-0.40319426317E-08;
const double B=0.21683205691E+03;
const double C=0.28463350538E-04;
const double D=-0.30774648337E-02;
double y;

bool negate=false;
if (x>2048)
{
negate=true;
x-=2048;
}
if (x>1024)
x=2048-x;
y=(A+x)/(B+C*x*x)+D*x;
if (negate)
return (float)(-y);
else
return (float)y;
}
```
```%This is Matlab code. you can convert it to C
%All it take to make a high quality sine
%wave is 1 multiply and one subtract.
%You first have to initialize the 2 unit delays
% and the coefficient

Fs = 48000;      %Sample rate
oscfreq = 1000.0; %Oscillator frequency in Hz
c1 = 2 * cos(2 * pi * oscfreq / Fs);
%Initialize the unit delays
d1 = sin(2 * pi * oscfreq / Fs);
d2 = 0;
%Initialization done here is the oscillator loop
% which generates a sinewave
for j=1:100
output = d1;        %This is the sine value
fprintf(1, '%f\n', output);
%one multiply and one subtract is all it takes
d0 = d1 * c1 - d2;
d2 = d1;   %Shift the unit delays
d1 = d0;
end
```
```Hi,

Can I use this code in a GPL2 or GPL3 licensed program (a soft synth project called Snarl)? In other words, will you grant permission for me to re-license your code? And what name should I write down as copyright holder in the headers?

Thanks,
Juuso Alasuutari
```
```Juuso,

Absolutely!

Toast
```

### Quick and dirty sine generator¶

notes
```this is part of my library, although I've seen a lot of sine generators, I've never seen
the simplest one, so I try to do it,
tell me something, I've try it and work so tell me something about it
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26``` ```PSPsample PSPsin1::doOsc(int numCh) { double x=0; double t=0; if(m_time[numCh]>m_sampleRate) //re-init cycle m_time[numCh]=0; if(m_time[numCh]>0) { t =(double)(((double)m_time[numCh])/(double)m_sampleRate); x=(m_2PI *(double)(t)*m_freq); } else x=0; PSPsample r=(PSPsample) sin(x+m_phase)*m_amp; m_time[numCh]++; return r; } ```

```isn't the sin() function a little bit heavyweight?  Since this is based upon slices of time, would it not be much more processor efficient to use a state variable filter that is self oscillating?

The operation:
t =(double)(((double)m_time[numCh])/(double)m_sampleRate);

also seems a little bit much, since t could be calculated by adding an interval value, which would eliminate the divide (needs more clocks).  The divide would then only need to be done once.
An FDIV may take 39 clock cycles minimum(depending on the operands), whilst an FADD is far faster (3 clocks).  An FMUL is comparable to an add, which would be a predominant instruction if using the SVF method.

FSIN may take between 16-126 clock cylces.

(clock cycle info nabbed from: http://www.singlix.com/trdos/pentium.txt)
```
```See also the fun with sinusoids page:
http://www.audiomulch.com/~rossb/code/sinusoids/
```
```For audio generation, sines are expensive i think, they are so perfect and take up more processing. it's rare to find a synth that sounds nicer with a sine compared to a parabol wave. My favourite parabolic wave is simply triangle wave with x*x with one of the half periods flipped. x*x is a very fast!!!
```
```hmm... x*x second half flipped...
very cool ! i'll give it a try!!
```

### RBJ Wavetable 101¶

notes
```see linked file
```

### Randja compressor¶

• Author or source: randja
• Type: compressor
• Created: 2009-09-29 07:37:05
notes
```I had found this code on the internet then made some improvements (speed) and now I post
it here for others to see.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111``` ```#include #define max(a,b) (a>b?a:b) class compressor { private: float threshold; float attack, release, envelope_decay; float output; float transfer_A, transfer_B; float env, gain; public: compressor() { threshold = 1.f; attack = release = envelope_decay = 0.f; output = 1.f; transfer_A = 0.f; transfer_B = 1.f; env = 0.f; gain = 1.f; } void set_threshold(float value) { threshold = value; transfer_B = output * pow(threshold,-transfer_A); } void set_ratio(float value) { transfer_A = value-1.f; transfer_B = output * pow(threshold,-transfer_A); } void set_attack(float value) { attack = exp(-1.f/value); } void et_release(float value) { release = exp(-1.f/value); envelope_decay = exp(-4.f/value); /* = exp(-1/(0.25*value)) */ } void set_output(float value) { output = value; transfer_B = output * pow(threshold,-transfer_A); } void reset() { env = 0.f; gain = 1.f; } __forceinline void process(float *input_left, float *input_right,float *output_left, float *output_right, int frames) { float det, transfer_gain; for(int i=0; i= env ? det : det+envelope_decay*(env-det); transfer_gain = env > threshold ? pow(env,transfer_A)*transfer_B:output; gain = transfer_gain < gain ? transfer_gain+attack *(gain-transfer_gain): transfer_gain+release*(gain-transfer_gain); output_left[i] = input_left[i] * gain; output_right[i] = input_right[i] * gain; } } __forceinline void process(double *input_left, double *input_right, double *output_left, double *output_right,int frames) { double det, transfer_gain; for(int i=0; i= env ? det : det+envelope_decay*(env-det); transfer_gain = env > threshold ? pow(env,transfer_A)*transfer_B:output; gain = transfer_gain < gain ? transfer_gain+attack *(gain-transfer_gain): transfer_gain+release*(gain-transfer_gain); output_left[i] = input_left[i] * gain; output_right[i] = input_right[i] * gain; } } }; ```

```env = det >= env ? det : det+envelope_decay*(env-det);

I have found it is good to add either a timed peak hold-before-release or something like a "peak magnetism" with an attack/release time that is within a cycle of the lowest expected frequency in the side chain...often ~80Hz is a good reference point... for example here is a snippet to illustrate the idea...please keep in mind this will not work as show because of some obvious things needing to be added...
class compressor {
float patk;  //"seek" slope
float ipatk;
float peakdet;
float target;
float envelope;  //the final value used for gain control
int hold;  //time to hold peak
int holdcnt;

patk = SAMPLE_RATE/80.0f;
ipatk = 1.0f - patk;
.
.
.
}
///now the envelope following function
float compressor::envelope_tracker(float input)
{
float rect = fabs(input);
holdcnt++;
if(rect>peakdet) {
peakdet = rect;
if(holdcnt>hold) //the hold is really optional
{
target = rect;
holdcnt = 0;
}
}
peakdet*=ipatk;  //sort of like capacitor+diode discharge.

if(envelope<target) envelope+=patk;
else envelope -= patk;  //sort of like a heat seeking missile ;)
//for the seek function this is really nothing more than a feedback control system, so a more elaborate higher order system could be used...
//the linear interpolation is good to use because
///it is simple while the main attack & release
//functions filter off the corners.
//In either case it is an improvement to a leaky
//peak detector.

}

Another thing I saw in the form of an analog circuit was a "round robin" peak detector which had 4 peak detectors being reset by a JFET to ground and was clocked by simple binary counter IC.  The full cycle from 0 to 15 on the binary counter was set at about 80Hz, then the JFET gates were reset every /4 of that period.  That way the envelope tracker pretty much grabbed every peak during the 1/80 period, then the "release" time was 1/80th second.  This could be implemented digitally very easily ;)
```
```Hi randja,

I am trying to use your compressor.
Could you explain in what unity shall be given the values of threshold, ratio, release, attack and output ?
Have you a good set of values for these parameters to make it work ?

```
```@vmeserette:

To quote the original code that randja "found" on the Internet and then "made some improvements" (he only added __forceinline):
"
How to use it:
--------------
free_comp.cpp and free_comp.hpp implement a simple C++ class called free_comp.
(People are expected to know what to do with it! If not seek help on a beginner
programming forum!)
The free_comp class implements the following methods:
void    set_threshold(float value);
Sets the threshold level; 'value' must be a _positive_ value
representing the amplitude as a floating point figure which should be
1.0 at 0dBFS

void    set_ratio(float value);
Sets the ratio; 'value' must be in range [0.0; 1.0] with 0.0
representing a oo:1 ratio, 0.5 a 2:1 ratio; 1.0 a 1:1 ratio and so on

void    set_attack(float value);
Sets the attack time; 'value' gives the attack time in _samples_

void    set_release(float value);
Sets the release time; 'value' gives the release time in _samples_

void    set_output(float value);
Sets the output gain; 'value' represents the gain, where 0dBFS is 1.0
(see set_threshold())

void    reset();
Resets the internal state of the compressor (gain reduction)

void    process(float *input_left, float *input_right,
float *output_left, float *output_right,
int frames, int skip);
void    process(double *input_left, double *input_right,
double *output_left, double *output_right,
int frames, int skip);
Processes a stereo stream of length 'frames' from either two arrays of
floats or arrays of doubles 'input_left' and 'input_right' then puts
the processed data in 'output_left' and 'output_right'.
'input_{left,right}' and 'output_{left,right}' may be the same location
in which case the algorithm will work in place. '{input,output}_left'
and '{input,output}_right' can also point to the same data, in which
case the algorithm works in mono (although if you process a lot of mono
data it will yield more performance if you modify the source to make the
algorithm mono in the first place).
The 'skip' parameter allows for processing of interleaved as well as two
separate contiguous streams. For two separate streams this value should
be 1, for interleaved stereo it should be 2 (but it can also have other
values than that to process specific channels in an interleaved audio
stream, though if you do that it is highly recommended to study the
source first to check whether it yields the expected behaviour or not).
```

### Rossler and Lorenz Oscillators¶

notes
```The Rossler and Lorenz functions are iterated chaotic systems - they trace smooth curves
that never repeat the same way twice. Lorenz is "unpitched", having no distinct peaks in
its spectrum -- similar to pink noise. Rossler exhibits definite spectral peaks against a

Time-domain and frequency spectrum of these two functions, as well as other info, can be
found at:

http://www.physics.emory.edu/~weeks/research/tseries1.html

These functions might be useful in simulating "analog drift."
```
code
 ```1 2``` ```Available on the web at: http://www.tinygod.com/code/BLorenzOsc.zip ```

```A Delphi/pascal version for VCL, KOL, Kylix and Freepascal on my website:
http://members.chello.nl/t.koning8/loro_sc.pas

Nice work!
```

### SawSin¶

• Author or source: Alexander Kritov
• Type: Oscillator shape
• Created: 2002-02-10 12:40:59
code
 ```1 2 3 4 5 6 7 8``` ```double sawsin(double x) { double t = fmod(x/(2*M_PI),(double)1.0); if (t>0.5) return -sin(x); if (t<=0.5) return (double)2.0*t-1.0; } ```

### Simple Time Stretching-Granular Synthesizer¶

• Author or source: Harry-Chris
• Created: 2008-12-18 07:08:20
notes
```Matlab function that implements crude time stretching - granulizing function, by overlap
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35``` ```function y = gran_func(x, w, H,H2, Fs, tr_amount) % x -> input signal % w -> Envelope - Window Vector % H1 -> Original Hop Size % H2 -> Synthesis Hop Size % Fs -> Sample Rate % str_amount -> time stretching factor M = length(w); pin = 1; pend = length(x) - M; y = zeros(1, floor( str_amount * length(x)) +M); count = 1; idx = 1; while pin < pend input = x(pin : pin+M-1) .* w'; y(idx : idx + M - 1) = y(idx : idx + M - 1) + input; pin = pin + H; count = count + 1; idx = idx + H2; end ```

### Sine calculation¶

• Author or source: Phil Burk
• Type: waveform generation, Taylor approximation of sin()
• Created: 2002-01-17 00:57:01
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28``` ```for(i=0; i < nSamples ; i++) { //Generate sawtooth phasor to provide phase for sine generation IncrementWrapPhase(phase, freqPtr[i]); //Wrap phase back into region where results are more accurate if(phase > 0.5) yp = 1.0 - phase; else { if(phase < -0.5) yp = -1.0 - phase; else yp = phase; } x = yp * PI; x2 = x*x; //Taylor expansion out to x**9/9! factored into multiply-adds fastsin = x*(x2*(x2*(x2*(x2*(1.0/362880.0) - (1.0/5040.0)) + (1.0/120.0)) - (1.0/6.0)) + 1.0); outPtr[i] = fastsin * amplPtr[i]; } ```

### Smooth random LFO Generator¶

• Author or source: Rob Belcham
• Created: 2009-06-30 08:31:24
notes
```I've been after a random LFO that's suitable for modulating a delay line for ages ( e.g
for chorus / reverb modulation) , so after i rolled my own, i thought i'd better make it
my first contribution to the music-dsp community.

My aim was to achive a sinusoidal based random but smooth waveform with a frequency
control that has no discontinuities and stays within a -1:1 range. If you listen to it, it
sounds quite like brown noise, or wind through a microphone (at rate = 100Hz for example)

It's written as a matlab m function, so shouldn't be too hard to port to C.

The oscillator generates a random level stepped waveform with random time spent at each
step (within bounds). These levels are linearly interpolated between and used to drive the
frequency of a sinewave. To achive amplitude variation, at each zero crossing a new random
amplitude scale factor is generated. The amplitude coefficient is ramped to this value
with a simple exponential.

An example call would be,
t = 4; Fs = 44100;
y = random_lfo(100, t*Fs, Fs);
axis([0, t*Fs, -1, 1]);
plot(y)

Enjoy !
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50``` ```% Random LFO Generator % creates a random sinusoidal waveform with no discontinuities % rate = average rate in Hz % N = run length in samples % Fs = sample frequency in Hz function y = random_lfo(rate, N, Fs) step_freq_scale = Fs / (1*rate); min_Cn = 0.1 * step_freq_scale; An = 0; lastA = 0; Astep = 0; y = zeros(1,N); % output x = 0; % sine phase lastSign = 0; amp_scale = 0.6; new_amp_scale = 0.6; amp_scale_ramp = exp(1000/Fs)-1; for (n=1:N) if (An == 0) || (An>=Cn) % generate a new random freq scale factor Cn = floor(step_freq_scale * rand()); % limit to prevent rapid transitions Cn = max(Cn, min_Cn); % generate new value & step coefficient newA = 0.1 + 0.9*rand(); Astep = (newA - lastA) / Cn; A = lastA; lastA = newA; % reset counter An = 0; end An = An + 1; % generate output y(n) = sin(x) * amp_scale; % ramp amplitude amp_scale = amp_scale + ( new_amp_scale - amp_scale ) * amp_scale_ramp; sin_inc = 2*pi*rate*A/Fs; A = A + Astep; % increment phase x = x + sin_inc; if (x >= 2*pi) x = x - 2*pi; end % scale at each zero crossing if (sign(y(n)) ~= 0) && (sign(y(n)) ~= lastSign) lastSign = sign(y(n)); new_amp_scale = 0.25 + 0.75*rand(); end; end; ```

### Square Waves¶

• Author or source: Sean Costello
• Type: waveform generation
• Created: 2002-01-15 21:26:38
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.
```

```Actually, I don't think this would work...
The proper way to do it is subtract a phase shifted buzz (aka BLIT) at the same frequency. This is equivalent to comb filtering, which will notch out the even harmonics.
```
```The above comment is correct, and my concept is inaccurate. My technique may have produced a signal with the proper harmonic structure, but it has been nearly 10 years since I wrote the post, so I can't remember what I was working with.

DFL's technique can be implemented with two buzz generators, or with a single buzz generator in conjunction with a fractional delay, where the delay controls the amount of phase shift.
```

### Trammell Pink Noise (C++ class)¶

code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48``` ```#ifndef _PinkNoise_H #define _PinkNoise_H // Technique by Larry "RidgeRat" Trammell 3/2006 // http://home.earthlink.net/~ltrammell/tech/pinkalg.htm // implementation and optimization by David Lowenfels #include #include #define PINK_NOISE_NUM_STAGES 3 class PinkNoise { public: PinkNoise() { srand ( time(NULL) ); // initialize random generator clear(); } void clear() { for( size_t i=0; i< PINK_NOISE_NUM_STAGES; i++ ) state[ i ] = 0.0; } float tick() { static const float RMI2 = 2.0 / float(RAND_MAX); // + 1.0; // change for range [0,1) static const float offset = A + A + A; // unrolled loop float temp = float( rand() ); state = P * (state - temp) + temp; temp = float( rand() ); state = P * (state - temp) + temp; temp = float( rand() ); state = P * (state - temp) + temp; return ( A*state + A*state + A*state )*RMI2 - offset; } protected: float state[ PINK_NOISE_NUM_STAGES ]; static const float A[ PINK_NOISE_NUM_STAGES ]; static const float P[ PINK_NOISE_NUM_STAGES ]; }; const float PinkNoise::A[] = { 0.02109238, 0.07113478, 0.68873558 }; // rescaled by (1+P)/(1-P) const float PinkNoise::P[] = { 0.3190, 0.7756, 0.9613 }; #endif ```

```Many thanks to David Lowenfels for posting this implementation of the early experimental version. I recommend switching to the new algorithm form described in 'newpink.htm' -- better range to 9+ octaves, better accuracy to +-0.25 dB, and leveled computational loading.  So where is MY submission to the archive? Um...  well, it's coming... if he doesn't beat me to the punch again and post his code first! -- Larry Trammell (the RidgeRat)
```

### Waveform generator using MinBLEPS¶

notes
```C code and project file for MSVC6 for a bandwidth-limited saw/square (with PWM) generator
using MinBLEPS.

This code is based on Eli's MATLAB MinBLEP code and uses his original minblep.mat file.
Instead of keeping a list of all active MinBLEPS, the output of each MinBLEP is stored in
a buffer, in which all consequent MinBLEPS and the waveform output are added together.
This optimization makes it fast enough to be used realtime.

Produces slight aliasing when sweeping high frequencies. I don't know wether Eli's
original code does the same, because I don't have MATLAB. Any help would be appreciated.

The project name is 'hardsync', because it's easy to generate hardsync using MinBLEPS.
```
code
 `1` ``` ```

```http://www.slack.net/~ant/bl-synth/windowed-impulse/

This page also describes a similar algorithm for generating waves. Could the aliasing be due to the fact that the blep only occurs after the discontinuity? On this page the blep also occurs in the opposite direction as well, leading up to the discontinuity.
```
• Date: 2008-02-11 18:42:15
• By: kernel[@}audiospillage.com
```The sawtooth is a nice oscillator but I can't seem to get the square wave to work properly.  Anyone else had any luck with this?  Also, it's worth noting that the code assumes it is running on a little endian architecture.
```
```I have written GPLv3 C++ source code for a MinBLEP oscillator and also public domain C++ source code for generating the MinBLEP without MatLab.

http://www.experimentalscene.com/articles/minbleps.php - Article and Code

http://www.experimentalscene.com/source.php - Look in DarkWave / latest version / CoreMachines / VCO.cpp
```

### Wavetable Synthesis¶

notes
```Wavetable sythesis AES paper by RBJ.
```

### Weird synthesis¶

• Author or source: Andy M00cho
• Created: 2002-01-17 00:55:09
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...
```

```Andy:
4.8 m.sec.fixed. What are the variables in the two controllable taps and is the 6dB filter variable frequency wise?
Phil
```
```What you have there is the core of a physical modelling algorithm.  I have done virtually the same thing to model plucked string instruments in Reaktor.  It's amazingly realistic.  See http://joeorgren.com
```

## Analysis¶

### Beat Detector Class¶

notes
```This class was designed for a VST plugin. Basically, it's just a 2nd order LP filter,
followed by an enveloppe detector (thanks Bram), feeding a Schmitt trigger. The rising
edge detector provides a 1-sample pulse each time a beat is detected. Code is self
documented...
Note : The class uses a fixed comparison level, you may need to change it.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89``` ```// ***** BEATDETECTOR.H ***** #ifndef BeatDetectorH #define BeatDetectorH class TBeatDetector { private: float KBeatFilter; // Filter coefficient float Filter1Out, Filter2Out; float BeatRelease; // Release time coefficient float PeakEnv; // Peak enveloppe follower bool BeatTrigger; // Schmitt trigger output bool PrevBeatPulse; // Rising edge memory public: bool BeatPulse; // Beat detector output TBeatDetector(); ~TBeatDetector(); virtual void setSampleRate(float SampleRate); virtual void AudioProcess (float input); }; #endif // ***** BEATDETECTOR.CPP ***** #include "BeatDetector.h" #include "math.h" #define FREQ_LP_BEAT 150.0f // Low Pass filter frequency #define T_FILTER 1.0f/(2.0f*M_PI*FREQ_LP_BEAT) // Low Pass filter time constant #define BEAT_RTIME 0.02f // Release time of enveloppe detector in second TBeatDetector::TBeatDetector() // Beat detector constructor { Filter1Out=0.0; Filter2Out=0.0; PeakEnv=0.0; BeatTrigger=false; PrevBeatPulse=false; setSampleRate(44100); } TBeatDetector::~TBeatDetector() { // Nothing specific to do... } void TBeatDetector::setSampleRate (float sampleRate) // Compute all sample frequency related coeffs { KBeatFilter=1.0/(sampleRate*T_FILTER); BeatRelease=(float)exp(-1.0f/(sampleRate*BEAT_RTIME)); } void TBeatDetector::AudioProcess (float input) // Process incoming signal { float EnvIn; // Step 1 : 2nd order low pass filter (made of two 1st order RC filter) Filter1Out=Filter1Out+(KBeatFilter*(input-Filter1Out)); Filter2Out=Filter2Out+(KBeatFilter*(Filter1Out-Filter2Out)); // Step 2 : peak detector EnvIn=fabs(Filter2Out); if (EnvIn>PeakEnv) PeakEnv=EnvIn; // Attack time = 0 else { PeakEnv*=BeatRelease; PeakEnv+=(1.0f-BeatRelease)*EnvIn; } // Step 3 : Schmitt trigger if (!BeatTrigger) { if (PeakEnv>0.3) BeatTrigger=true; } else { if (PeakEnv<0.15) BeatTrigger=false; } // Step 4 : rising edge detector BeatPulse=false; if ((BeatTrigger)&&(!PrevBeatPulse)) BeatPulse=true; PrevBeatPulse=BeatTrigger; } ```

```// Nice work!
//Here's a Delphi and freepascal version:
unit beattrigger;

interface

type
TBeatDetector = class
private
KBeatFilter,               // Filter coefficient
Filter1Out,
Filter2Out,
BeatRelease,               // Release time coefficient
PeakEnv:single;            // Peak enveloppe follower
BeatTrigger,               // Schmitt trigger output
PrevBeatPulse:Boolean;     // Rising edge memory
public
BeatPulse:Boolean;            // Beat detector output
constructor Create;
procedure setSampleRate(SampleRate:single);
procedure AudioProcess (input:single);
end;

function fabs(value:single):Single;

implementation

const
FREQ_LP_BEAT = 150.0;                    // Low Pass filter frequency
T_FILTER = 1.0/(2.0 * PI*FREQ_LP_BEAT);  // Low Pass filter time constant
BEAT_RTIME = 0.02;   // Release time of enveloppe detector in second

constructor TBeatDetector.create;
// Beat detector constructor
begin
inherited;
Filter1Out:=0.0;
Filter2Out:=0.0;
PeakEnv:=0.0;
BeatTrigger:=false;
PrevBeatPulse:=false;
setSampleRate(44100);
end;

procedure TBeatDetector.setSampleRate (sampleRate:single);
// Compute all sample frequency related coeffs
begin
KBeatFilter:=1.0/(sampleRate*T_FILTER);
BeatRelease:= exp(-1.0/(sampleRate*BEAT_RTIME));
end;

function fabs(value:single):Single;
asm
fld value
fabs
fwait
end;

procedure  TBeatDetector.AudioProcess (input:single);
var
EnvIn:Single;
// Process incoming signal
begin
// Step 1 : 2nd order low pass filter (made of two 1st order RC filter)
Filter1Out:=Filter1Out+(KBeatFilter*(input-Filter1Out));
Filter2Out:=Filter2Out+(KBeatFilter*(Filter1Out-Filter2Out));
// Step 2 : peak detector
EnvIn:=fabs(Filter2Out);
if EnvIn>PeakEnv then PeakEnv:=EnvIn  // Attack time = 0
else
begin
PeakEnv:=PeakEnv*BeatRelease;
PeakEnv:=PeakEnv+(1.0-BeatRelease)*EnvIn;
end;
// Step 3 : Schmitt trigger
if not BeatTrigger then
begin
if PeakEnv>0.3 then BeatTrigger:=true;
end
else
begin
if PeakEnv<0.15 then BeatTrigger:=false;
end;

// Step 4 : rising edge detector
BeatPulse:=false;
if (BeatTrigger = true ) and( not PrevBeatPulse) then
BeatPulse:=true;
PrevBeatPulse:=BeatTrigger;
end;

end.
```

### Coefficients for Daubechies wavelets 1-38¶

• Author or source: Computed by Kazuo Hatano, Compiled and verified by Olli Niemitalo
• Type: wavelet transform
• Created: 2002-01-17 02:00:43
• Linked files: `daub.h`.

### DFT¶

• Author or source: Andy Mucho
• Type: fourier transform
• Created: 2002-01-17 01:59:38
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26``` ```AnalyseWaveform(float *waveform, int framesize) { float aa[MaxPartials]; float bb[MaxPartials]; for(int i=0;i

### Envelope detector¶

• Author or source: Bram
• Created: 2002-04-12 21:37:18
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23``` ```//attack and release in seconds float ga = (float) exp(-1/(SampleRate*attack)); float gr = (float) exp(-1/(SampleRate*release)); float envelope=0; for(...) { //get your data into 'input' EnvIn = std::abs(input); if(envelope < EnvIn) { envelope *= ga; envelope += (1-ga)*EnvIn; } else { envelope *= gr; envelope += (1-gr)*EnvIn; } //envelope now contains.........the envelope ;) } ```

```// Slightly faster version of the envelope follower using one multiply form.

// attTime and relTime is in seconds

float ga = exp(-1.0f/(sampleRate*attTime));
float gr = exp(-1.0f/(sampleRate*relTime));

float envOut = 0.0f;

for( ... )
{
// get your data into 'input'

envIn = fabs(input);

if( envOut < envIn )
envOut = envIn + ga * (envOut - envIn);
else
envOut = envIn + gr * (envOut - envIn);

// envOut now contains the envelope
}
```

### Envelope follower with different attack and release¶

• Author or source: Bram
• Created: 2003-01-15 00:21:39
notes
```xxxx_in_ms is xxxx in milliseconds ;-)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13``` ```init:: attack_coef = exp(log(0.01)/( attack_in_ms * samplerate * 0.001)); release_coef = exp(log(0.01)/( release_in_ms * samplerate * 0.001)); envelope = 0.0; loop:: tmp = fabs(in); if(tmp > envelope) envelope = attack_coef * (envelope - tmp) + tmp; else envelope = release_coef * (envelope - tmp) + tmp; ```

```// the expressions of the form:

xxxx_coef = exp(log(0.01)/( xxxx_in_ms * samplerate * 0.001));

// can be simplified a little bit to:

xxxx_coef = pow(0.01, 1.0/( xxxx_in_ms * samplerate * 0.001));
```
```Here the definition of the attack/release time is the time for the envelope to fall from 100% to 1%.
In the other version, the definition is for the envelope to fall from 100% to 36.7%. So in this one
the envelope is about 4.6 times faster.
```

### 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)."
```

```Thank you very much, this was useful, and it worked right out of the box,
so to speak. It's very efficient, and the algorithm is readable. It also includes some
very useful functions.
```

### FFT classes in C++ and Object Pascal¶

• Author or source: Laurent de Soras (Object Pascal translation by Frederic Vanmol)
• Type: Real-to-Complex FFT and Complex-to-Real IFFT
• Created: 2002-02-14 02:09:26
• Linked files: `FFTReal.zip`.
notes
```(see linkfile)
```

• Author or source: Timo H Tossavainen
• Type: wavelet transform
• Created: 2002-01-17 01:54:52
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37``` ```void inline wht_bfly (long& a, long& b) { long tmp = a; a += b; b = tmp - b; } // just a integer log2 int inline l2 (long x) { int l2; for (l2 = 0; x > 0; x >>=1) { ++ l2; } return (l2); } //////////////////////////////////////////// // Fast in-place Walsh-Hadamard Transform // //////////////////////////////////////////// void FWHT (std::vector& data) { const int log2 = l2 (data.size()) - 1; for (int i = 0; i < log2; ++i) { for (int j = 0; j < (1 << log2); j += 1 << (i+1)) { for (int k = 0; k < (1<

### Frequency response from biquad coefficients¶

notes
```Here is a formula for plotting the frequency response of a biquad filter. Depending on the
coefficients that you have, you might have to use negative values for the b- coefficients.
```
code
 ```1 2 3 4 5``` ```//w = frequency (0 < w < PI) //square(x) = x*x y = 20*log((sqrt(square(a0*square(cos(w))-a0*square(sin(w))+a1*cos(w)+a2)+square(2*a0*cos(w)*sin(w)+a1*(sin(w))))/ sqrt(square(square(cos(w))- square(sin(w))+b1*cos(w)+b2)+square(2* cos(w)*sin(w)+b1*(sin(w)))))); ```

```this formula can have roundoff errors with frequencies close to zero... (especially a problem
with high samplerate filters)

here is a better formula:

20*log10[|H(e^jw)|] =
10*log10[ (b0+b1+b2)^2 - 4*(b0*b1 + 4*b0*b2 + b1*b2)*phi + 16*b0*b2*phi^2 ]
-10*log10[ (a0+a1+a2)^2 - 4*(a0*a1 + 4*a0*a2 + a1*a2)*phi + 16*a0*a2*phi^2 ]

where phi = sin^2(w/2)
```

### Java FFT¶

• Author or source: Loreno Heer
• Type: FFT Analysis
• Created: 2003-11-25 17:38:15
notes
```May not work correctly ;-)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142``` ```// WTest.java /* Copyright (C) 2003 Loreno Heer, (helohe at bluewin dot ch) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ public class WTest{ private static double[] sin(double step, int size){ double f = 0; double[] ret = new double[size]; for(int i = 0; i < size; i++){ ret[i] = Math.sin(f); f += step; } return ret; } private static double[] add(double[] a, double[] b){ double[] c = new double[a.length]; for(int i = 0; i < a.length; i++){ c[i] = a[i] + b[i]; } return c; } private static double[] sub(double[] a, double[] b){ double[] c = new double[a.length]; for(int i = 0; i < a.length; i++){ c[i] = a[i] - b[i]; } return c; } private static double[] add(double[] a, double b){ double[] c = new double[a.length]; for(int i = 0; i < a.length; i++){ c[i] = a[i] + b; } return c; } private static double[] cp(double[] a, int size){ double[] c = new double[size]; for(int i = 0; i < size; i++){ c[i] = a[i]; } return c; } private static double[] mul(double[] a, double b){ double[] c = new double[a.length]; for(int i = 0; i < a.length; i++){ c[i] = a[i] * b; } return c; } private static void print(double[] value){ for(int i = 0; i < value.length; i++){ System.out.print(i + "," + value[i] + "\n"); } System.out.println(); } private static double abs(double[] a){ double c = 0; for(int i = 0; i < a.length; i++){ c = ((c * i) + Math.abs(a[i])) / (i + 1); } return c; } private static double[] fft(double[] a, int min, int max, int step){ double[] ret = new double[(max - min) / step]; int i = 0; for(int d = min; d < max; d = d + step){ double[] f = sin(fc(d), a.length); double[] dif = sub(a, f); ret[i] = 1 - abs(dif); i++; } return ret; } private static double[] fft_log(double[] a){ double[] ret = new double; int i = 0; for(double d = 0; d < 15.5; d = d + 0.01){ double[] f = sin(fc(Math.pow(2,d)), a.length); double[] dif = sub(a, f); ret[i] = Math.abs(1 - abs(dif)); i++; } return ret; } private static double fc(double d){ return d * Math.PI / res; } private static void print_log(double[] value){ for(int i = 0; i < value.length; i++){ System.out.print(Math.pow(2,((double)i/100d)) + "," + value[i] + "\n"); } System.out.println(); } public static void main(String[] args){ double[] f_0 = sin(fc(440), sample_length); // res / pi =>14005 //double[] f_1 = sin(.02, sample_length); double[] f_2 = sin(fc(520), sample_length); //double[] f_3 = sin(.25, sample_length); //double[] f = add( add( add(f_0, f_1), f_2), f_3); double[] f = add(f_0, f_2); //print(f); double[] d = cp(f,1000); print_log(fft_log(d)); } static double length = .2; // sec static int res = 44000; // resoultion (pro sec) static int sample_length = res; // resoultion } ```

### LPC analysis (autocorrelation + Levinson-Durbin recursion)¶

notes
```The autocorrelation function implements a warped autocorrelation, so that frequency
resolution can be specified by the variable 'lambda'.  Levinson-Durbin recursion
calculates autoregression coefficients a and reflection coefficients (for lattice filter
implementation) K.  Comments for Levinson-Durbin function implement matlab version of the
same function.

No optimizations.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80``` ```//find the order-P autocorrelation array, R, for the sequence x of length L and warping of lambda //wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0); wAutocorrelate(float * x, unsigned int L, float * R, unsigned int P, float lambda) { double * dl = new double [L]; double * Rt = new double [L]; double r1,r2,r1t; R=0; Rt=0; r1=0; r2=0; r1t=0; for(unsigned int k=0; k

```// Blind Object Pascal Translation:
// --------------------------------

unit Levinson;

interface

type
TDoubleArray = array of Double;
TSingleArray = array of Single;

implementation

//find the P-order autocorrelation array, R, for the sequence x of length L and warping of lambda
procedure Autocorrelate(x,R : TSingleArray; P : Integer; lambda : Single; l: Integer = -1);
var dl,Rt     : TDoubleArray;
r1,r2,r1t : Double;
k,i       : Integer;
begin
// Initialization
if l=-1 then l:=Length(x);
SetLength(dl,l);
SetLength(Rt,l);
R:=0;
Rt:=0;
r1:=0;
r2:=0;
r1t:=0;

for k:=0 to l-1 do
begin
Rt:=Rt+x[k]*x[k];
dl[k]:=r1-lambda*(x[k]-r2);
r1:= x[k];
r2:= dl[k];
end;

for i:=1 to P do
begin
Rt[i]:=0;
r1:=0;
r2:=0;
for k:=0 to L-1 do
begin
Rt[i]:=Rt[i]+dl[k]*x[k];
r1t:= dl[k];
dl[k]:=r1-lambda*(r1t-r2);
r1:=r1t;
r2:=dl[k];
end;
end;

for i:=1 to P do R[i]:=Rt[i];
setlength(Rt,0);
setlength(dl,0);
end;

// Calculate the Levinson-Durbin recursion for the autocorrelation sequence
// R of length P+1 and return the autocorrelation coefficients a and reflection coefficients K
procedure LevinsonRecursion(P : Integer; R,A,K : TSingleArray);
var Am1       : TDoubleArray;
i,j,s,m   : Integer;
km,Em1,Em : Double;
err       : Double;
begin
SetLength(Am1,62);
if (R=0.0) then
begin
for i:=1 to P do
begin
K[i]:=0.0;
A[i]:=0.0;
end;
end
else
begin
for j:=0 to P do
begin
A:=0;
Am1:=0;
end;
A:=1;
Am1:=1;
km:=0;
Em1:=R;
for m:=1 to P do
begin
err:=0.0;
for j:=1 to m-1 do err:=err+Am1[j]*R[m-j];
km:=(R[m]-err)/Em1;
K[m-1]:=-km;
A[m]:=km;
for j:=1 to m-1 do A[j]:=Am1[j]-km*Am1[m-j];
Em:=(1-km*km)*Em1;
for s:=0 to P do Am1[s]:=A[s];
Em1:=Em;
end;
end;
end;

end.
```

• Author or source: Wilfried Welti
• Created: 2002-01-17 03:08:11
notes
```use add_value with all values which enter the look-ahead area,
and remove_value with all value which leave this area. to get
the maximum value in the look-ahead area, use get_max_value.
in the very beginning initialize the table with zeroes.

If you always want to know the maximum amplitude in
problem. very primitive approach using a look-up table
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63``` ```void lookup_add(unsigned section, unsigned size, unsigned value) { if (section==value) lookup[section]++; else { size >>= 1; if (value>section) { lookup[section]++; lookup_add(section+size,size,value); } else lookup_add(section-size,size,value); } } void lookup_remove(unsigned section, unsigned size, unsigned value) { if (section==value) lookup[section]--; else { size >>= 1; if (value>section) { lookup[section]--; lookup_remove(section+size,size,value); } else lookup_remove(section-size,size,value); } } unsigned lookup_getmax(unsigned section, unsigned size) { unsigned max = lookup[section] ? section : 0; size >>= 1; if (size) if (max) { max = lookup_getmax((section+size),size); if (!max) max=section; } else max = lookup_getmax((section-size),size); return max; } void add_value(unsigned value) { lookup_add(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value); } void remove_value(unsigned value) { lookup_remove(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1, value); } unsigned get_max_value() { return lookup_getmax(LOOKUP_VALUES>>1, LOOKUP_VALUES>>1); } ```

### Magnitude and phase plot of arbitrary IIR function, up to 5th order¶

• Author or source: George Yohng
• Type: magnitude and phase at any frequency
• Created: 2002-08-01 00:43:57
notes
```Amplitude and phase calculation of IIR equation
run at sample rate "sampleRate" at frequency "F".

AMPLITUDE
-----------
cf_mag(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------

PHASE
-----------
cf_phi(F,sampleRate,
a0,a1,a2,a3,a4,a5,
b0,b1,b2,b3,b4,b5)
-----------

If you need a frequency diagram, draw a plot for
F=0...sampleRate/2

If you need amplitude in dB, use cf_lin2db(cf_mag(.......))

Set b0=-1 if you have such function:

y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
+ b1*y[n-1] + b2*y[n-2] + b3*y[n-3] + b4*y[n-4] + b5*y[n-5];

Set b0=1  if you have such function:

y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] + a3*x[n-3] + a4*x[n-4] + a5*x[n-5] +
- b1*y[n-1] - b2*y[n-2] - b3*y[n-3] - b4*y[n-4] - b5*y[n-5];

Do not try to reverse engineer these formulae - they don't give any sense
other than they are derived from transfer function, and they work. :)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73``` ```/* C file can be downloaded from http://www.yohng.com/dsp/cfsmp.c */ #define C_PI 3.14159265358979323846264 double cf_mag(double f,double rate, double a0,double a1,double a2,double a3,double a4,double a5, double b0,double b1,double b2,double b3,double b4,double b5) { return sqrt((a0*a0 + a1*a1 + a2*a2 + a3*a3 + a4*a4 + a5*a5 + 2*(a0*a1 + a1*a2 + a2*a3 + a3*a4 + a4*a5)*cos((2*f*C_PI)/rate) + 2*(a0*a2 + a1*a3 + a2*a4 + a3*a5)*cos((4*f*C_PI)/rate) + 2*a0*a3*cos((6*f*C_PI)/rate) + 2*a1*a4*cos((6*f*C_PI)/rate) + 2*a2*a5*cos((6*f*C_PI)/rate) + 2*a0*a4*cos((8*f*C_PI)/rate) + 2*a1*a5*cos((8*f*C_PI)/rate) + 2*a0*a5*cos((10*f*C_PI)/rate))/ (b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 + 2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) + 2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) + 2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) + 2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) + 2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate))); } double cf_phi(double f,double rate, double a0,double a1,double a2,double a3,double a4,double a5, double b0,double b1,double b2,double b3,double b4,double b5) { atan2((a0*b0 + a1*b1 + a2*b2 + a3*b3 + a4*b4 + a5*b5 + (a0*b1 + a1*(b0 + b2) + a2*(b1 + b3) + a5*b4 + a3*(b2 + b4) + a4*(b3 + b5))*cos((2*f*C_PI)/rate) + ((a0 + a4)*b2 + (a1 + a5)*b3 + a2*(b0 + b4) + a3*(b1 + b5))*cos((4*f*C_PI)/rate) + a3*b0*cos((6*f*C_PI)/rate) + a4*b1*cos((6*f*C_PI)/rate) + a5*b2*cos((6*f*C_PI)/rate) + a0*b3*cos((6*f*C_PI)/rate) + a1*b4*cos((6*f*C_PI)/rate) + a2*b5*cos((6*f*C_PI)/rate) + a4*b0*cos((8*f*C_PI)/rate) + a5*b1*cos((8*f*C_PI)/rate) + a0*b4*cos((8*f*C_PI)/rate) + a1*b5*cos((8*f*C_PI)/rate) + (a5*b0 + a0*b5)*cos((10*f*C_PI)/rate))/ (b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 + 2*((b0*b1 + b1*b2 + b3*(b2 + b4) + b4*b5)*cos((2*f*C_PI)/rate) + (b2*(b0 + b4) + b3*(b1 + b5))*cos((4*f*C_PI)/rate) + (b0*b3 + b1*b4 + b2*b5)*cos((6*f*C_PI)/rate) + (b0*b4 + b1*b5)*cos((8*f*C_PI)/rate) + b0*b5*cos((10*f*C_PI)/rate))), ((a1*b0 + a3*b0 + a5*b0 - a0*b1 + a2*b1 + a4*b1 - a1*b2 + a3*b2 + a5*b2 - a0*b3 - a2*b3 + a4*b3 - a1*b4 - a3*b4 + a5*b4 - a0*b5 - a2*b5 - a4*b5 + 2*(a3*b1 + a5*b1 - a0*b2 + a4*(b0 + b2) - a1*b3 + a5*b3 + a2*(b0 - b4) - a0*b4 - a1*b5 - a3*b5)*cos((2*f*C_PI)/rate) + 2*(a3*b0 + a4*b1 + a5*(b0 + b2) - a0*b3 - a1*b4 - a0*b5 - a2*b5)* cos((4*f*C_PI)/rate) + 2*a4*b0*cos((6*f*C_PI)/rate) + 2*a5*b1*cos((6*f*C_PI)/rate) - 2*a0*b4*cos((6*f*C_PI)/rate) - 2*a1*b5*cos((6*f*C_PI)/rate) + 2*a5*b0*cos((8*f*C_PI)/rate) - 2*a0*b5*cos((8*f*C_PI)/rate))*sin((2*f*C_PI)/rate))/ (b0*b0 + b1*b1 + b2*b2 + b3*b3 + b4*b4 + b5*b5 + 2*(b0*b1 + b1*b2 + b2*b3 + b3*b4 + b4*b5)*cos((2*f*C_PI)/rate) + 2*(b0*b2 + b1*b3 + b2*b4 + b3*b5)*cos((4*f*C_PI)/rate) + 2*b0*b3*cos((6*f*C_PI)/rate) + 2*b1*b4*cos((6*f*C_PI)/rate) + 2*b2*b5*cos((6*f*C_PI)/rate) + 2*b0*b4*cos((8*f*C_PI)/rate) + 2*b1*b5*cos((8*f*C_PI)/rate) + 2*b0*b5*cos((10*f*C_PI)/rate))); } double cf_lin2db(double lin) { if (lin<9e-51) return -1000; /* prevent invalid operation */ return 20*log10(lin); } ```

• Date: 2004-01-02 08:46:35
• By: Rob
```% Actually it is simpler to simply take the zero-padded b and a coefficients and do real->complex
% FFT like this (matlab code):

H_complex=fft(b,N)./fft(a,N);
phase=angle(H_complex);
Magn=abs(H_complex);

% This will give you N/2 points from 0 to pi angle freq (or 0 to nyquist freq).
% /Rob
```
```// Here are the formulas if you only have a biquad. But i am not sure, if maybe the
// phase is shifted with pi/2...

20*Log10(
sqrt(
(a0*a0+a1*a1+a2*a2+
2*(a0*a1+a1*a2)*cos(w)+
2*(a0*a2)* cos(2*w)
)
/
(
1 + b1*b1 + b2*b2 +
2*(b1 + b1*b2)*cos(w)+
2*b2*cos(2*w)
)
)
)

ArcTan2(
(
a0+a1*b1+a2*b2+
(a0*b1+a1*(1+b2)+a2*b1)*cos(w)+
(a0*b2+a2)*cos(2*w)
)
/
(
1+b1*b1+b2*b2+
2*
(
(b1+b1*b2)*cos(w)+ b2*cos(2*w)
)
)
,
(
(
a1-a0*b1+a2*b1-a1*b2+
2*(-a0*b2+a2)*cos(w)
)*sin(w)
/
(
1+b1*b1+b2*b2+
2*(b1 + b1*b2)*cos(w)+
2*b2*cos(2*w)
)
)
)
```
```// Recursive Delphi Code with arbitrary order:

unit Plot;

interface

type TArrayOfDouble = Array of Double;

function MagnitudeCalc(f,rate : Double; a,b : TArrayOfDouble): Double;

implementation

uses Math;

function MulVectCalc(const v: TArrayOfDouble; const Z, N : Integer) : Double;
begin
if N=0
then result:=0
else result:=(v[N-1]*v[N-1+Z])+MulVectCalc(v,Z,N-1);
end;

function MagCascadeCalc(const v: TArrayOfDouble; const w : double; N, Order : Integer ): Double;
begin
if N=1
then result:=(MulVectCalc(v,0,Order))
else result:=((MulVectCalc(v,N-1,1+Order-N)*(2*cos((N-1)*w))+MagCascadeCalc(v, w, N-1, Order )));
end;

function MagnitudeCalc(f,rate : Double; a,b : TArrayOfDouble): Double;
var w : Double;
begin
w:=(2*f*pi)/rate;
end;

end.
```
```function CalcMagPart(w: Double; C : TDoubleArray):Double;
var i,j,l : Integer;
temp  : Double;
begin
l:=Length(C);
temp:=0;
for j:=0 to l-1
do temp:=temp+C[j]*C[j];
result:=temp;
for i:=1 to l-1 do
begin
temp:=0;
for j:=0 to l-i-1
do temp:=temp+C[j]*C[j+i];
result:=Result+2*temp*cos(i*w);
end;
end;

function CalcMagnitude_dB(const f,rate: Double; const A,B: TDoubleArray): Double;
var w : Double;
begin
w:=(2*f*pi)/rate;
result:=10*log10(CalcMagPart(w,A)/CalcMagPart(w,B));
end;

// Here's a really fast function for an arbitrary IIR with high order without stack overflows
// or recursion. And specially for John without sqrt.
```

### Measuring interpollation noise¶

• Author or source: Jon Watte
• Created: 2002-01-17 02:00:09
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.
```

### QFT and DQFT (double precision) classes¶

notes
```Since it's a Visual C++ project (though it has relatively portable C++) I
guess the main audience are PC users.  As such I'm including a zip file.
Some PC users wouldn't know what to do with a tgz file.

The QFT and DQFT (double precision) classes supply the following functions:

1. Real valued FFT and inverse FFT functions.  Note that separate arraysare used
for real and imaginary component of the resulting spectrum.

2. Decomposition of a spectrum into a separate spectrum of the evensamples
and a spectrum of the odd samples.  This can be useful for buildingfilter banks.

3. Reconstituting a spectrum from separate spectrums of the even samples
and odd samples.  This can be useful for building filter banks.

4. A discrete Sin transform (a QFT decomposes an FFT into a DST and DCT).

5. A discrete Cos transfrom.

6. Since a QFT does it's last stage calculating from the outside in thelast part
can be left unpacked and only calculated as needed in the case wherethe entire
spectrum isn't needed (I used this for calculating correlations andconvolutions
where I only needed half of the results).
ReverseNoUnpack()
UnpackStep()
and NegUnpackStep()
implement this functionality

NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseNoUnpack() does not.

7. Also if you only want the first half of the results you can call ReverseHalf()

NOTE Reverse() normalizes its results (divides by one half the blocklength), but
ReverseHalf() does not.

8. QFT is less numerically stable than regular FFTs.  With singleprecision calculations,
a block length of 2^15 brings the accuracy down to being barelyaccurate enough.
At that size, single precision calculations tested sound files wouldoccasionally have
a sample off by 2, and a couple off by 1 per block. Full volume whitenoise would
generate
a few samples off by as much as 6 per block at the end, beginning and middle.

No matter what the inputs the errors are always at the same positions in the block.
There some sort of cancelation that gets more delicate as the block size gets bigger.

For the sake of doing convolutions and the like where the forward transform is
done only once for one of the inputs, I created a AccurateForward() function.
It uses a regular FFT algorithm for blocks larger than 2^12, and decomposes into even
and
odd FFTs recursively.

In any case you can always use the double precision routines to get more accuracy.
DQFT even has routines that take floats as inputs and return double precision
spectrum outputs.

As for portability:

#define _USE_ASM

If you comment those define out, then what's left is C++ with no assembly language.

2. There is unnecessary windows specific code in "criticalSection.h"
I used a critical section because objects are not reentrant (each object has
permanent scratch pad memory), but obviously critical sections are operating
system specific.  In any case that code can easily be taken out.

If you look at my code and see that there's an a test built in the examples
that makes sure that the results are in the ballpark of being right. It
wasn't that I expected the answers to be far off, it was that I uncommenting
the "no assembly language" versions of some routines and I wanted to make
sure that they weren't broken.
```

### Simple peak follower¶

• Author or source: Phil Burk
• Type: amplitude analysis
• Created: 2002-01-17 01:57:19
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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25``` ```// halfLife = time in seconds for output to decay to half value after an impulse static float output = 0.0; float scalar = pow( 0.5, 1.0/(halfLife * sampleRate))); if( input < 0.0 ) input = -input; /* Absolute value. */ if ( input >= output ) { /* When we hit a peak, ride the peak to the top. */ output = input; } else { /* Exponential decay of output when signal is low. */ output = output * scalar; /* ** When current gets close to 0.0, set current to 0.0 to prevent FP underflow ** which can cause a severe performance degradation due to a flood ** of interrupts. */ if( output < VERY_SMALL_FLOAT ) output = 0.0; } ```

```#ifndef VERY_SMALL_FLOAT
#define VERY_SMALL_FLOAT 1.0e-30F
#endif
```

### Tone detection with Goertzel¶

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
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107``` ```/** Tone detect by Goertzel algorithm * * This program basically searches for tones (sines) in a sample and reports the different dB it finds for * different frequencies. Can easily be extended with some thresholding to report true/false on detection. * I'm far from certain goertzel it implemented 100% correct, but it works :) * * Hint, the SAMPLERATE, BUFFERSIZE, FREQUENCY, NOISE and SIGNALVOLUME all affects the outcome of the reported dB. Tweak * em to find the settings best for your application. Also, seems to be pretty sensitive to noise (whitenoise anyway) which * is a bit sad. Also I don't know if the goertzel really likes float values for the frequency ... And using 44100 as * samplerate for detecting 6000 Hz tone is kinda silly I know :) * * Written by: Espen Riskedal, espenr@ii.uib.no, july-2002 */ #include #include #include using std::rand; // math stuff using std::cos; using std::abs; using std::exp; using std::log10; // iostream stuff using std::cout; using std::endl; #define PI 3.14159265358979323844 // change the defines if you want to #define SAMPLERATE 44100 #define BUFFERSIZE 8820 #define FREQUENCY 6000 #define NOISE 0.05 #define SIGNALVOLUME 0.8 /** The Goertzel algorithm computes the k-th DFT coefficient of the input signal using a second-order filter. * http://ptolemy.eecs.berkeley.edu/papers/96/dtmf_ict/www/node3.html. * Basiclly it just does a DFT of the frequency we want to check, and none of the others (FFT calculates for all frequencies). */ float goertzel(float *x, int N, float frequency, int samplerate) { float Skn, Skn1, Skn2; Skn = Skn1 = Skn2 = 0; for (int i=0; i 1) tone[i]= 2*cos(A)*tone[i-1] - tone[i-2]; else if (i > 0) tone[i] = 2*cos(A)*tone[i-1] - (cos(A)); else tone[i] = 2*cos(A)*cos(A) - cos(2*A); } for (int i=0; i dB: " << power(goertzel(tone, BUFFERSIZE, freq, SAMPLERATE)) << endl; } delete tone; return 0; } ```

```// yet untested Delphi translation of the algorithm:

function Goertzel(Buffer:array of Single; frequency, samplerate: single):single;
var Skn, Skn1, Skn2 : Single;
i               : Integer;
temp1, temp2    : Single;
begin
skn:=0;
skn1:=0;
skn2:=0;
temp1:=2*PI*frequency/samplerate;
temp2:=Cos(temp1);
for i:=0 to Length(Buffer) do
begin
Skn2 = Skn1;
Skn1 = Skn;
Skn = 2*temp2*Skn1 - Skn2 + Buffer[i];
end;
Result:=(Skn - exp(-temp1)*Skn1);
end;

// Maybe someone can use it...
//
// Christian
```

### Tone detection with Goertzel (x86 ASM)¶

notes
```This is an "assemblified" version of the Goertzel Tone Detector. It is about 2 times
faster than the original code.

The code has been tested and it works fine.

Hope you can use it. I'm gonna try to build a Tuner (as VST-Plugin). I hope, that this
will work :-\ If anyone is intrested, please let me know.

Christian
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54``` ```function Goertzel_x87(Buffer :Psingle; BLength:Integer; frequency: Single; samplerate: Single):Single; asm mov ecx,BLength mov eax,Buffer fld x2 fldpi fmulp fmul frequency fdiv samplerate fld st(0) fcos fld x2 fmulp fxch st(1) fldz fsub st(0),st(1) fstp st(1) fldl2e fmul fld st(0) frndint fsub st(1),st(0) fxch st(1) f2xm1 fld1 fadd fscale fstp st(1) fldz fldz fldz @loopStart: fxch st(1) fxch st(2) fstp st(0) fld st(3) fmul st(0),st(1) fsub st(0),st(2) fld [eax].Single faddp add eax,4 loop @loopStart @loopEnd: fxch st(3) fmulp st(2), st(0) fsub st(0),st(1) fstp result ffree st(2) ffree st(1) ffree st(0) end; ```

```// Here's a variant on the theme that compensates for harmonics:

Function Goertzel(.Buffer: array of double; frequency, samplerate: double):.double;
var
Qkn, Qkn1, Qkn2, Wkn, Mk: double;
i: integer;
begin
Qkn:=0; Qkn1:=0;
Wkn:=2*.PI*.frequency/samplerate;
Mk:=2*.Cos(.Wkn);
for i:=0 to High(.Buffer) do begin
Qkn2: = Qkn1; Qkn1: = Qkn;
Qkn  : = Buffer[.i ] + Mk*.Qkn1 - Qkn2;
end;
Result: = sqrt(.Qkn*.Qkn + Qkn1*.Qkn1 - Mk*.Qkn*.Qkn1);
end;

// Posted on www.delphimaster.ru by Jeer
```

### Vintage VU meters tutorial¶

notes
```Here is a short tutorial about vintage-styled VU meters:

http://neolit123.blogspot.com/2009/03/designing-analog-vu-meter-in-dsp.html
```

## Filters¶

### 1 pole LPF for smooth parameter changes¶

• Author or source: moc.liamg@odiugoiz
• Type: 1-pole LPF class
• Created: 2008-09-22 20:27:06
notes
```This is a very simple class that I'm using in my plugins for smoothing parameter changes
that directly affect audio stream.
It's a 1-pole LPF, very easy on CPU.

Change the value of variable "a" (0~1) for slower or a faster response.

Of course you can also use it as a lowpass filter for audio signals.
```
code
 ```1 2 3 4 5 6 7 8 9``` ```class CParamSmooth { public: CParamSmooth() { a = 0.99f; b = 1.f - a; z = 0; }; ~CParamSmooth(); inline float Process(float in) { z = (in * b) + (z * a); return z; } private: float a, b, z; }; ```

```I've used this a lot, but here's an important thing: it won't work the same for multiple sample rates, so if you set a to 0.9995 for example, this will be lower if the sample rate is higher than you intended. I fix it by doing this:

/*must compensate this factor for sample rate change*/

float srCompensate;
srCompensate = sr/44100.0f;
float compensated_a;
compensated_a = powf(a, (1.0f/srCompensate));
b = 1.0f-compensated_a;

Then if you start with a built in value (or range) designed for 44100hz, it will scale up with the sample rate so you will get the same amount of smoothing. I'm not sure if this is mathematically correct, but I came up with it very quickly and it works a charm for me.

Good work, very useful and easy to use filter.
```
```*edit, a won't be LOWER if the sample rate changes, but it won't have the same effect.
```
```New version, now you can specify the speed response of the parameter in ms. and sampling rate:

class CParamSmooth
{
public:
CParamSmooth(float smoothingTimeInMs, float samplingRate)
{
const float c_twoPi = 6.283185307179586476925286766559f;

a = exp(-c_twoPi / (smoothingTimeInMs * 0.001f * samplingRate));
b = 1.0f - a;
z = 0.0f;
}

~CParamSmooth()
{

}

inline float process(float in)
{
z = (in * b) + (z * a);
return z;
}

private:
float a;
float b;
float z;
};
```

### 1-RC and C filter¶

notes
```This filter is called 1-RC and C since it uses these two parameters. C and R correspond to
raw cutoff and inverted resonance, and have a range from 0 to 1.
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12``` ```//Parameter calculation //cutoff and resonance are from 0 to 127 c = pow(0.5, (128-cutoff) / 16.0); r = pow(0.5, (resonance+24) / 16.0); //Loop: v0 = (1-r*c)*v0 - (c)*v1 + (c)*input; v1 = (1-r*c)*v1 + (c)*v0; output = v1; ```

• Date: 2005-01-13 18:25:57
• By: yes
```input is not in 0 - 1 range.

for cutoff i guess 128.

for reso the same ?
```
```Nice. This is very similar to a state variable filter in many ways. Relationship between c and frequency:

c = 2*sin(pi*freq/samplerate)

You can approximate this (tuning error towards nyquist):

c = 2*pi*freq/samplerate

Relationship between r and q factor:

r = 1/q

This filter has stability issues for high r values. State variable filter stability limits seem to work fine here. It can also be oversampled for better stability and wider frequency range (use 0.5*original frequency):

//Loop:

v0 = (1-r*c)*v0 - c*v1 + c*input;
v1 = (1-r*c)*v1 + c*v0;
tmp = v1;

v0 = (1-r*c)*v0 - c*v1 + c*input;
v1 = (1-r*c)*v1 + c*v0;
output = (tmp+v1)*0.5;

-- peter schoffhauzer
```

### 18dB/oct resonant 3 pole LPF with tanh() dist¶

notes
```Implementation in CSound and Sync Modular...
```

### 1st and 2nd order pink noise filters¶

notes
```Here are some new lower-order pink noise filter coefficients.

These have approximately equiripple error in decibels from 20hz to 20khz at a 44.1khz
sampling rate.

1st order, ~ +/- 3 dB error (not recommended!)
num = [0.05338071119116  -0.03752455712906]
den = [1.00000000000000  -0.97712493947102]

2nd order, ~ +/- 0.9 dB error
num = [ 0.04957526213389  -0.06305581334498   0.01483220320740 ]
den = [ 1.00000000000000  -1.80116083982126   0.80257737639225 ]
```

### 3 Band Equaliser¶

• Author or source: Neil C
• Created: 2006-08-29 20:34:25
notes
```Simple 3 band equaliser with adjustable low and high frequencies ...

Fairly fast algo, good quality output (seems to be accoustically transparent with all
gains set to 1.0)

How to use ...

1. First you need to declare a state for your eq

EQSTATE eq;

2. Now initialise the state (we'll assume your output frequency is 48Khz)

set_3band_state(eq,880,5000,480000);

Your EQ bands are now as follows (approximatley!)

low  band = 0Hz to 880Hz
mid  band = 880Hz to 5000Hz
high band = 5000Hz to 24000Hz

3. Set the gains to some values ...

eq.lg = 1.5;  // Boost bass by 50%
eq.mg = 0.75; // Cut mid by 25%
eq.hg = 1.0;  // Leave high band alone

4. You can now EQ some samples

out_sample = do_3band(eq,in_sample)

Have fun and mail me if any problems ...  etanza at lycos dot co dot uk

Neil C / Etanza Systems, 2006 :)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203``` ```First the header file .... //--------------------------------------------------------------------------- // // 3 Band EQ :) // // EQ.H - Header file for 3 band EQ // // (c) Neil C / Etanza Systems / 2K6 // // Shouts / Loves / Moans = etanza at lycos dot co dot uk // // This work is hereby placed in the public domain for all purposes, including // use in commercial applications. // // The author assumes NO RESPONSIBILITY for any problems caused by the use of // this software. // //---------------------------------------------------------------------------- #ifndef __EQ3BAND__ #define __EQ3BAND__ // ------------ //| Structures | // ------------ typedef struct { // Filter #1 (Low band) double lf; // Frequency double f1p0; // Poles ... double f1p1; double f1p2; double f1p3; // Filter #2 (High band) double hf; // Frequency double f2p0; // Poles ... double f2p1; double f2p2; double f2p3; // Sample history buffer double sdm1; // Sample data minus 1 double sdm2; // 2 double sdm3; // 3 // Gain Controls double lg; // low gain double mg; // mid gain double hg; // high gain } EQSTATE; // --------- //| Exports | // --------- extern void init_3band_state(EQSTATE* es, int lowfreq, int highfreq, int mixfreq); extern double do_3band(EQSTATE* es, double sample); #endif // #ifndef __EQ3BAND__ //--------------------------------------------------------------------------- Now the source ... //---------------------------------------------------------------------------- // // 3 Band EQ :) // // EQ.C - Main Source file for 3 band EQ // // (c) Neil C / Etanza Systems / 2K6 // // Shouts / Loves / Moans = etanza at lycos dot co dot uk // // This work is hereby placed in the public domain for all purposes, including // use in commercial applications. // // The author assumes NO RESPONSIBILITY for any problems caused by the use of // this software. // //---------------------------------------------------------------------------- // NOTES : // // - Original filter code by Paul Kellet (musicdsp.pdf) // // - Uses 4 first order filters in series, should give 24dB per octave // // - Now with P4 Denormal fix :) //---------------------------------------------------------------------------- // ---------- //| Includes | // ---------- #include #include "eq.h" // ----------- //| Constants | // ----------- static double vsa = (1.0 / 4294967295.0); // Very small amount (Denormal Fix) // --------------- //| Initialise EQ | // --------------- // Recommended frequencies are ... // // lowfreq = 880 Hz // highfreq = 5000 Hz // // Set mixfreq to whatever rate your system is using (eg 48Khz) void init_3band_state(EQSTATE* es, int lowfreq, int highfreq, int mixfreq) { // Clear state memset(es,0,sizeof(EQSTATE)); // Set Low/Mid/High gains to unity es->lg = 1.0; es->mg = 1.0; es->hg = 1.0; // Calculate filter cutoff frequencies es->lf = 2 * sin(M_PI * ((double)lowfreq / (double)mixfreq)); es->hf = 2 * sin(M_PI * ((double)highfreq / (double)mixfreq)); } // --------------- //| EQ one sample | // --------------- // - sample can be any range you like :) // // Note that the output will depend on the gain settings for each band // (especially the bass) so may require clipping before output, but you // knew that anyway :) double do_3band(EQSTATE* es, double sample) { // Locals double l,m,h; // Low / Mid / High - Sample Values // Filter #1 (lowpass) es->f1p0 += (es->lf * (sample - es->f1p0)) + vsa; es->f1p1 += (es->lf * (es->f1p0 - es->f1p1)); es->f1p2 += (es->lf * (es->f1p1 - es->f1p2)); es->f1p3 += (es->lf * (es->f1p2 - es->f1p3)); l = es->f1p3; // Filter #2 (highpass) es->f2p0 += (es->hf * (sample - es->f2p0)) + vsa; es->f2p1 += (es->hf * (es->f2p0 - es->f2p1)); es->f2p2 += (es->hf * (es->f2p1 - es->f2p2)); es->f2p3 += (es->hf * (es->f2p2 - es->f2p3)); h = es->sdm3 - es->f2p3; // Calculate midrange (signal - (low + high)) m = es->sdm3 - (h + l); // Scale, Combine and store l *= es->lg; m *= es->mg; h *= es->hg; // Shuffle history buffer es->sdm3 = es->sdm2; es->sdm2 = es->sdm1; es->sdm1 = sample; // Return result return(l + m + h); } //---------------------------------------------------------------------------- ```

```Great Thanks!
I have one problem the below:
double  f2p0;     // Poles ...
double  f2p1;
double  f2p2;
double  f2p3;
that I want to know the starting value
```
```yuri:

The invocation of memset() during the initialization method sets all the the members of the struct to zero.
```
```This is great -- I want to develop a compressor/limiter/expander and have been looking long and hard for bandpass / eq filtering code.  Here it is!

I am sure we could easily expand this into an x band eq.

Thanks!
```
• Date: 2007-07-05 06:49:45
• By: tom tom
```Hi !

I've just transposed your code under Delphi.

It works well if the gain is under 1, but if i put gain > 1 i get clipping (annoying sound clips), even at 1.1;

Is it normal ?

I convert my smallint (44100 16 bits) to double before process, and convert the obtained value back to smallint with clipping (if < -32768 i set it to -32768, and if > 32768 i set it to 32768).

What did i do wrong ?

Regards

Tom
```
```Hi.

Maybe the answer is quite easy. The upper limit is 32767 not 32768.

Regards

Herbert
```
```Hi, Can U send me a full source code for this 3 band state eq from start to end ??
I really need it for my study in school.
I hope you can send me, to my email.

thanks you.
regard

angga
```
```How can I expand this 3 Band EQ into X Band EQ..?!

Anybody answer me, or email me..
```
```For more bands, you could take the low-pass and repeat the process on that.
```
```This is a great little filter, I am using it in an application but when I first started playing with it I noticed some problems. The mid range didn't seem to be calculated properly, a friend of mine who knows more about dsp than I do took a quick look at it and suggested the following change:

m          = es->sdm3 - (h + l);

Should be:

m          = sample - (h + l);

I've tested it with this small fix and everything works perfectly now. Just thought I'd bring this to your attention... Thanks for a great code snippet!
```
```What problems were you getting? Doesn't removing the delay cause phase problems?
```
```Hi Great Stuff,

how to create 6 band equalizer, is any algorithm for 6 band same like 3 band, please help me if any one

```
```How to extend this to 6 band equalizer?
```
```hello! thanks for your code!!
i tried to use the code in my project of guitar distortions in real time (in C) and i could'nt, i'm in linux using jack audio server, and it starts to have x-runs everytime i turn on the equalizer. do you any idea of how solving this? (from 5 to 50 milliseconds o x-runs)
i was thinking of coding it in assembler but i don't know if that would be the solution.
excuse me for my english, i'm from argentina and it's been a while since i last wrote in this language!
mariano
```
```Dev c++ can not run it? any suggesstions, how to run it?
```
```This example is exactly what i've been looking for. This little piece of code executes faster then the one I have been using before.

FYI,

I will use the code in a 3 band compressor / limiter / clipper for FM broadcasting.

Thanks for sharing.

Best regards,

Vincent Bruinink.
```
```I've got this filtering audio on iOS by running my sample through do_3band in the render callback. However, I'm getting a fair amount of distortion. Here's an example with my EQ3Band gains all set to 1.0:

Here's the code for my implementation:

https://github.com/tassock/mixerhost/commit/4b8b87028bfffe352ed67609f747858059a3e89b

I assume others using this aren't having this same distortion issue? If so, what sort of audio sample formats are you using (big/little endian, float/integer samples, etc). Thanks!
```
```hi,

Lucie
```
```      I've made an implementation of 3 band Equalizer to read a wave file,apply filtering and then save wave  outputfile.
With your code I have a lot of distortion,I think the problem maybe coefficient calculation:

es->f1p0  += (es->lf * (sample   - es->f1p0))+vsa ;
...

Can anyone resolve distortion?
```
```Works fine for me on iOS. Maybe you feed a interleaved stereo signal with the same EQSTATE instance (you'll need one EQSTATE for each channel)?
```
```It's all about WHERE you init you EQ.  Try a little :) If you don't find out yourself, I'll help you.
```
```I added this code to my xcode project. But where i can pass the values and method. May be its funny question for you guyz but please help me. Its loking good to me but in xcode , i am using avaudio player for play sound and making sound app.
Thanks
```
```The distortion will occur of you are trying to adapt this to stereo and do so by simply adding an outer loop per channel. Doing so will cause the filter values to compound and cause distortion. Instead, you will need to duplicate all the filter values and keep them separate from the other channel.
```
```Any good x-band equalizer equivalents for C# that I can use?
```
```s'cool thanks!
have portet to c#
works fine!
```
```Converted Neil's C Code 3 band equalizer to Delphi class, for those who are interested.

1. Create instance of class
Public
eq:TEQ;

2. On form create

eq := TEq.Create;

//Initialize
init_3band_state(880,5000,44100,50,-25,0);

//init_3band_state(lowfreq,highfreq,mixfreq:integer;BassGain,MidGain,HighGain:Double);

3. process: pass Raw 16Bit PCM to eq

eq.Equalize(const Data: Pointer; DataSize: DWORD);

Works like a charm form me

***********************************************
TEq=Class
private
lf,f1p0,f1p1,f1p2,f1p3:double;
hf,f2p0,f2p1,f2p2,f2p3:double;
sdm1,sdm2,sdm3:double;
vsa:double;
lg,mg,hg:double;
Function do_3band(sample:Smallint):Smallint;
public
constructor create;
destructor destroy;override;
procedure init_3band_state(lowfreq,highfreq,mixfreq:integer;BassGain,MidGain,HighGain:Double);
procedure Equalize(const Data: Pointer; DataSize: DWORD);
end;

constructor TEQ.create;
begin
inherited create;
vsa := (1.0 / 4294967295.0);
lg := 1.0;
mg := 1.0;
hg := 1.0;

end;

destructor TEQ.destroy;
begin
inherited destroy;
end;

procedure TEQ.init_3band_state(lowfreq,highfreq,mixfreq:integer;BassGain,MidGain,HighGain:Double);
begin

{(880,5000,44100,1.5,0.75,1.0)
eq.lg = 1.5; // Boost bass by 50%
eq.mg = 0.75; // Cut mid by 25%
eq.hg = 1.0; // Leave high band alone }

lg := 1+(BassGain/100);
mg := 1+(MidGain/100);
hg := 1+(HighGain/100);

// Calculate filter cutoff frequencies

lf := 2 * sin(PI * (lowfreq / mixfreq));
hf := 2 * sin(PI * (highfreq / mixfreq));
end;

Function TEQ.do_3band(sample:Smallint):Smallint;
var l,m,h:double;
res:integer;
begin

// Filter #1 (lowpass)

f1p0 := f1p0 + (lf * (sample - f1p0)) + vsa;
f1p1 := f1p1 + (lf * (f1p0 - f1p1));
f1p2 := f1p2 + (lf * (f1p1 - f1p2));
f1p3 := f1p3 + (lf * (f1p2 - f1p3));

l := f1p3;

// Filter #2 (highpass)

f2p0 := f2p0 + (hf * (sample - f2p0)) + vsa;
f2p1 := f2p1 + (hf * (f2p0 - f2p1));
f2p2 := f2p2 + (hf * (f2p1 - f2p2));
f2p3 := f2p3 + (hf * (f2p2 - f2p3));

h := sdm3 - f2p3;

// Calculate midrange (signal - (low + high))

m := sdm3 - (h + l);

// Scale, Combine and store

l := l * lg;
m := m * mg;
h := h * hg;

// Shuffle history buffer

sdm3 := sdm2;
sdm2 := sdm1;
sdm1 := sample;

// Return result
res := trunc(l+m+h);
if res > 32767 then res := 32767 else if res < -32768 then res := -32768;

result := res;
end;

procedure TEQ.Equalize(const Data: Pointer; DataSize: DWORD);
var pSample: PSmallInt;
begin
pSample := Data;
while DataSize > 0 do
begin
pSample^ := do_3band(pSample^);
Inc(pSample);
Dec(DataSize, 2);
end;
end;
```

### 303 type filter with saturation¶

notes
```I posted a filter to the Csound mailing list a couple of weeks ago that has a 303 flavor
to it. It basically does wacky distortions to the sound. I used Runge-Kutta for the diff
eq. simulation though which makes it somewhat sluggish.

This is a CSound score!!
```

notes
```Original from T. Lossius - ttblue project

Optimized version in pseudo-code.

[! The filter is unstable for fast automation changes in the lower frequency range.
Parameter interpolation and/or oversampling should fix this. !]

The sum of the Linkwitz-Riley (Butterworth squared) HP and LP outputs, will result an all-
pass filter at Fc and flat magnitude response - close to ideal for crossovers.

Lubomir I. Ivanov
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63``` ```//----------------------------------------- // [code] //----------------------------------------- //fc -> cutoff frequency //pi -> 3.14285714285714 //srate -> sample rate //================================================ // shared for both lp, hp; optimizations here //================================================ wc=2*pi*fc; wc2=wc*wc; wc3=wc2*wc; wc4=wc2*wc2; k=wc/tan(pi*fc/srate); k2=k*k; k3=k2*k; k4=k2*k2; sqrt2=sqrt(2); sq_tmp1=sqrt2*wc3*k; sq_tmp2=sqrt2*wc*k3; a_tmp=4*wc2*k2+2*sq_tmp1+k4+2*sq_tmp2+wc4; b1=(4*(wc4+sq_tmp1-k4-sq_tmp2))/a_tmp; b2=(6*wc4-8*wc2*k2+6*k4)/a_tmp; b3=(4*(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp; b4=(k4-2*sq_tmp1+wc4-2*sq_tmp2+4*wc2*k2)/a_tmp; //================================================ // low-pass //================================================ a0=wc4/a_tmp; a1=4*wc4/a_tmp; a2=6*wc4/a_tmp; a3=a1; a4=a0; //===================================================== // high-pass //===================================================== a0=k4/a_tmp; a1=-4*k4/a_tmp; a2=6*k4/a_tmp; a3=a1; a4=a0; //===================================================== // sample loop - same for lp, hp //===================================================== tempx=input; tempy=a0*tempx+a1*xm1+a2*xm2+a3*xm3+a4*xm4-b1*ym1-b2*ym2-b3*ym3-b4*ym4; xm4=xm3; xm3=xm2; xm2=xm1; xm1=tempx; ym4=ym3; ym3=ym2; ym2=ym1; ym1=tempy; output=tempy; ```

```LR2 with DFII:

//------------------------------
// LR2
// fc -> cutoff frequency
// pi -> 3.14285714285714
// srate -> sample rate
//------------------------------
fpi = pi*fc;
wc = 2*fpi;
wc2 = wc*wc;
wc22 = 2*wc2;
k = wc/tan(fpi/srate);
k2 = k*k;
k22 = 2*k2;
wck2 = 2*wc*k;
tmpk = (k2+wc2+wck2);
//b shared
b1 = (-k22+wc22)/tmpk;
b2 = (-wck2+k2+wc2)/tmpk;
//---------------
// low-pass
//---------------
a0_lp = (wc2)/tmpk;
a1_lp = (wc22)/tmpk;
a2_lp = (wc2)/tmpk;
//----------------
// high-pass
//----------------
a0_hp = (k2)/tmpk;
a1_hp = (-k22)/tmpk;
a2_hp = (k2)/tmpk;

//=========================
// sample loop, in -> input
//=========================
//---lp
lp_out = a0_lp*in + lp_xm0;
lp_xm0 = a1_lp*in - b1*lp_out + lp_xm1;
lp_xm1 = a2_lp*in - b2*lp_out;
//---hp
hp_out = a0_hp*in + hp_xm0;
hp_xm0 = a1_hp*in - b1*hp_out + hp_xm1;
hp_xm1 = a2_hp*in - b2*hp_out;

// the two are with 180 degrees phase shift,
// so you need to invert the phase of one.
out = lp_out + hp_out*(-1);

//result is allpass at Fc
```
```I've converted this Linkwits Riley 4 into intrinsics. It's set up with the cross over point and the sample rate. The function 'ProcessSplit' returns the low and high parts. It uses _mm_malloc to align the variables to 16 bytes, as putting them into the class as __m128 vars doesn't guarantee alignment.
Enjoy!  :)

//-----------------------------------------------
//-----------------------------------------------
//-----------------------------------------------

#pragma once

#include <xmmintrin.h>

{
__m128 *ab;
__m128 *al;
__m128 *ah;
__m128 *xm;
__m128 *yml;
__m128 *ymh;

float a0l;
float a0h;

public:

void ResetSplit();

__inline void FIL_Linkwitz_Riley4::ProcessSplit(const float in, float &low, float &high)
{
__m128 m1;

m1 = _mm_sub_ps(_mm_mul_ps(*al, *xm), _mm_mul_ps(*ab, *yml));
low  = a0l * in + m1.m128_f32 + m1.m128_f32 + m1.m128_f32 + m1.m128_f32;

m1 = _mm_sub_ps(_mm_mul_ps(*ah, *xm), _mm_mul_ps(*ab, *ymh));
high = a0h * in + m1.m128_f32 + m1.m128_f32 + m1.m128_f32 + m1.m128_f32;

*xm  = _mm_shuffle_ps(*xm, *xm, _MM_SHUFFLE(2,1,0,0));
(*xm).m128_f32 = in;
*yml = _mm_shuffle_ps(*yml, *yml, _MM_SHUFFLE(2,1,0,0));
(*yml).m128_f32 = low;
*ymh = _mm_shuffle_ps(*ymh, *ymh, _MM_SHUFFLE(2,1,0,0));
(*ymh).m128_f32 = high;
}

};

//-----------------------------------------------
//-----------------------------------------------
//-----------------------------------------------

#include <math.h>

{
ab  = (__m128*)_mm_malloc(16, 16);
al  = (__m128*)_mm_malloc(16, 16);
ah  = (__m128*)_mm_malloc(16, 16);
xm  = (__m128*)_mm_malloc(16, 16);
yml = (__m128*)_mm_malloc(16, 16);
ymh = (__m128*)_mm_malloc(16, 16);

float wc  = 2.0f * PI * fc;
float wc2 = wc*wc;
float wc3 = wc2*wc;
float wc4 = wc2*wc2;
float k  = wc / tanf(PI * fc / srate);
float k2 = k*k;
float k3 = k2*k;
float k4 = k2*k2;
float sqrt2   = sqrtf(2.0f);
float sq_tmp1 = sqrt2 *wc3 * k;
float sq_tmp2 = sqrt2 *wc * k3;
float a_tmp       = 4.0f * wc2 * k2 + 2.0f * sq_tmp1 + k4 + 2.0f * sq_tmp2 + wc4;

(*ab).m128_f32 = (4.0f *(wc4+sq_tmp1-k4-sq_tmp2))/a_tmp;
(*ab).m128_f32 = (6.0f *wc4-8*wc2*k2+6*k4)/a_tmp;
(*ab).m128_f32 = (4.0f *(wc4-sq_tmp1+sq_tmp2-k4))/a_tmp;
(*ab).m128_f32 = (k4 -2.0f * sq_tmp1 + wc4 - 2.0f * sq_tmp2 + 4.0f * wc2 * k2) / a_tmp;

//================================================
// low-pass
//================================================
a0l       = wc4/a_tmp;
(*al).m128_f32 = 4.0f * wc4 / a_tmp;
(*al).m128_f32 = 6.0f * wc4 / a_tmp;
(*al).m128_f32 = (*al).m128_f32;
(*al).m128_f32 = a0l;

//=====================================================
// high-pass
//=====================================================
a0h       = k4 / a_tmp;
(*ah).m128_f32 = -4.0f * k4 / a_tmp;
(*ah).m128_f32 =  6.0f * k4 / a_tmp;
(*ah).m128_f32 = (*ah).m128_f32;
(*ah).m128_f32 = a0h;

ResetSplit();
}

{
_mm_free((void*)ab);
_mm_free((void*)al);
_mm_free((void*)ah);
_mm_free((void*)xm);
_mm_free((void*)yml);
_mm_free((void*)ymh);
}

{
// Reset history...
*xm  = _mm_set1_ps(0.0f);
*yml = _mm_set1_ps(0.0f);
*ymh = _mm_set1_ps(0.0f);
}
//-----------------------------------------------
//-----------------------------------------------
//-----------------------------------------------
```
```I've no idea why I'm accessing those pointers like that! But never mind. :)
```
```Your pi value is wrong:
pi -> 3.14285714285714

It should be 3.1415692 ect.
```
```Or even 3.1415926535!
LOL.
```
```I don't think this is unstable for changes in frequency.  It's unstable for low frequencies.

Here's my implementation.

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <assert.h>

class LRCrossoverFilter { // LR4 crossover filter
private:
struct filterCoefficents {
float a0, a1, a2, a3, a4;
} lpco, hpco;

float b1co, b2co, b3co, b4co;

struct {
float xm1 = 0.0f;
float xm2 = 0.0f;
float xm3 = 0.0f;
float xm4 = 0.0f;
float ym1 = 0.0f, ym2 = 0.0f, ym3 = 0.0f, ym4 = 0.0f;
} hptemp, lptemp;

float coFreqRunningAv = 100.0f;
public:
void setup(float crossoverFrequency, float sr);
void processBlock(float * in, float * outHP, float * outLP, int numSamples);
void dumpCoefficents(struct filterCoefficents x) {
std::cout << "a0: " << x.a0 << "\n";
std::cout << "a1: " << x.a1 << "\n";
std::cout << "a2: " << x.a2 << "\n";
std::cout << "a3: " << x.a3 << "\n";
std::cout << "a4: " << x.a4 << "\n";
}
void dumpInformation() {
std::cout << "-----\nfrequency: "<< coFreqRunningAv << "\n";
std::cout << "lpco:\n";
dumpCoefficents(lpco);
std::cout << "hpco:\n";
dumpCoefficents(hpco);
std::cout << "bco:\nb1: ";
std::cout << b1co << "\nb2: " << b2co << "\nb3: " <<  b3co << "\nb4: " << b4co << "\n";
}

};

void LRCrossoverFilter::setup(float crossoverFrequency, float sr) {

const float pi = 3.141f;

coFreqRunningAv = crossoverFrequency;

float cowc=2*pi*coFreqRunningAv;
float cowc2=cowc*cowc;
float cowc3=cowc2*cowc;
float cowc4=cowc2*cowc2;

float cok=cowc/tan(pi*coFreqRunningAv/sr);
float cok2=cok*cok;
float cok3=cok2*cok;
float cok4=cok2*cok2;
float sqrt2=sqrt(2);
float sq_tmp1 = sqrt2 * cowc3 * cok;
float sq_tmp2 = sqrt2 * cowc * cok3;
float a_tmp = 4*cowc2*cok2 + 2*sq_tmp1 + cok4 + 2*sq_tmp2+cowc4;

b1co=(4*(cowc4+sq_tmp1-cok4-sq_tmp2))/a_tmp;

b2co=(6*cowc4-8*cowc2*cok2+6*cok4)/a_tmp;

b3co=(4*(cowc4-sq_tmp1+sq_tmp2-cok4))/a_tmp;

b4co=(cok4-2*sq_tmp1+cowc4-2*sq_tmp2+4*cowc2*cok2)/a_tmp;

//================================================
// low-pass
//================================================
lpco.a0=cowc4/a_tmp;
lpco.a1=4*cowc4/a_tmp;
lpco.a2=6*cowc4/a_tmp;
lpco.a3=lpco.a1;
lpco.a4=lpco.a0;

//=====================================================
// high-pass
//=====================================================
hpco.a0=cok4/a_tmp;
hpco.a1=-4*cok4/a_tmp;
hpco.a2=6*cok4/a_tmp;
hpco.a3=hpco.a1;
hpco.a4=hpco.a0;

}

void LRCrossoverFilter::processBlock(float * in, float * outHP, float * outLP, int numSamples) {

float tempx, tempy;
for (int i = 0; i<numSamples; i++) {
tempx=in[i];

// High pass

tempy = hpco.a0*tempx +
hpco.a1*hptemp.xm1 +
hpco.a2*hptemp.xm2 +
hpco.a3*hptemp.xm3 +
hpco.a4*hptemp.xm4 -
b1co*hptemp.ym1 -
b2co*hptemp.ym2 -
b3co*hptemp.ym3 -
b4co*hptemp.ym4;

hptemp.xm4=hptemp.xm3;
hptemp.xm3=hptemp.xm2;
hptemp.xm2=hptemp.xm1;
hptemp.xm1=tempx;
hptemp.ym4=hptemp.ym3;
hptemp.ym3=hptemp.ym2;
hptemp.ym2=hptemp.ym1;
hptemp.ym1=tempy;
outHP[i]=tempy;

assert(tempy<10000000);

// Low pass

tempy = lpco.a0*tempx +
lpco.a1*lptemp.xm1 +
lpco.a2*lptemp.xm2 +
lpco.a3*lptemp.xm3 +
lpco.a4*lptemp.xm4 -
b1co*lptemp.ym1 -
b2co*lptemp.ym2 -
b3co*lptemp.ym3 -
b4co*lptemp.ym4;

lptemp.xm4=lptemp.xm3; // these are the same as hptemp and could be optimised away
lptemp.xm3=lptemp.xm2;
lptemp.xm2=lptemp.xm1;
lptemp.xm1=tempx;
lptemp.ym4=lptemp.ym3;
lptemp.ym3=lptemp.ym2;
lptemp.ym2=lptemp.ym1;
lptemp.ym1=tempy;
outLP[i] = tempy;

assert(!isnan(outLP[i]));
}
}

int main () {
LRCrossoverFilter filter;
float data;
float lp, hp;

filter.setup(50.0, 44100.0f);
filter.dumpInformation();

for (int i = 0; i<2000; i++) {
data[i] = sinf(i/100.f);
}
filter.processBlock(data, hp, lp, 2000);

}

I'll try and fix it, but this kind of work is new to me, so all suggestions appreciated (Including "You Fool, you've copied the code wrong"). cheers!
```
```I tried this code for a crossover - firstly the SSE intrinsics version then the full original version. Both have problems with the HPF output.
With a crossover frequency of 200Hz and a pure sine tone input (any pitch) I get loud (-16dBFS)low frequency noise in the HPF output. This noise level reduces as the crossover frequency increases but it is unusable in its current state.
Can anyone post a solution for this problem?
Thanks.....Chris
```
```I also tried the LR2 code, this works better but there is still low frequency noise (-56dBFS & Xover 200Hz) in the HPF output.
Seems there is a fundamental problem with the HPF coefficients in this code :(
The LF Noise for both LR2 and LR4 appears to be a modulating DC offset - maybe that can guide the Filter Gurus to identify and solve the problem.
Cheers.....Chris
```

### Another 4-pole lowpass…¶

• Author or source: ten.xmg@zlipzzuf
• Type: 4-pole LP/HP
• Created: 2004-09-06 08:40:52
notes
```Vaguely based on the Stilson/Smith Moog paper, but going in a rather different direction
from others I've seen here.

The parameters are peak frequency and peak magnitude (g below); both are reasonably
accurate for magnitudes above 1. DC gain is 1.

The filter has some undesirable properties - e.g. it's unstable for low peak freqs if
implemented in single precision (haven't been able to cleanly separate it into biquads or
onepoles to see if that helps), and it responds so strongly to parameter changes that it's
not advisable to update the coefficients much more rarely than, say, every eight samples
during sweeps, which makes it somewhat expensive.

I like the sound, however, and the accuracy is nice to have, since many filters are not
very strong in that respect.

I haven't looked at the HP again for a while, but IIRC it had approximately the same good
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64``` ```double coef; double d; double omega; //peak freq double g; //peak mag // calculating coefficients: double k,p,q,a; double a0,a1,a2,a3,a4; k=(4.0*g-3.0)/(g+1.0); p=1.0-0.25*k;p*=p; // LP: a=1.0/(tan(0.5*omega)*(1.0+p)); p=1.0+a; q=1.0-a; a0=1.0/(k+p*p*p*p); a1=4.0*(k+p*p*p*q); a2=6.0*(k+p*p*q*q); a3=4.0*(k+p*q*q*q); a4= (k+q*q*q*q); p=a0*(k+1.0); coef=p; coef=4.0*p; coef=6.0*p; coef=4.0*p; coef=p; coef=-a1*a0; coef=-a2*a0; coef=-a3*a0; coef=-a4*a0; // or HP: a=tan(0.5*omega)/(1.0+p); p=a+1.0; q=a-1.0; a0=1.0/(p*p*p*p+k); a1=4.0*(p*p*p*q-k); a2=6.0*(p*p*q*q+k); a3=4.0*(p*q*q*q-k); a4= (q*q*q*q+k); p=a0*(k+1.0); coef=p; coef=-4.0*p; coef=6.0*p; coef=-4.0*p; coef=p; coef=-a1*a0; coef=-a2*a0; coef=-a3*a0; coef=-a4*a0; // per sample: out=coef*in+d; d=coef*in+coef*out+d; d=coef*in+coef*out+d; d=coef*in+coef*out+d; d=coef*in+coef*out; ```

```Yet untested object pascal translation:

unit T4PoleUnit;

interface

type TFilterType=(ftLowPass, ftHighPass);
T4Pole=class(TObject)
private
fGain     : Double;
fFreq     : Double;
fSR       : Single;
protected
fCoeffs     : array[0..8] of Double;
d           : array[0..3] of Double;
fFilterType : TFilterType;
procedure SetGain(s:Double);
procedure SetFrequency(s:Double);
procedure SetFilterType(v:TFilterType);
procedure Calc;
public
constructor Create;
function Process(s:single):single;
published
property Gain: Double read fGain write SetGain;
property Frequency: Double read fFreq write SetFrequency;
property SampleRate: Single read fSR write fSR;
property FilterType: TFilterType read fFilterType write SetFilterType;
end;

implementation

uses math;

const kDenorm = 1.0e-25;

constructor T4Pole.Create;
begin
inherited create;
fFreq:=1000;
fSR:=44100;
Calc;
end;

procedure T4Pole.SetFrequency(s:Double);
begin
fFreq:=s;
Calc;
end;

procedure T4Pole.SetGain(s:Double);
begin
fGain:=s;
Calc;
end;

procedure T4Pole.SetFilterType(v:TFilterType);
begin
fFilterType:=v;
Calc;
end;

procedure T4Pole.Calc;
var k,p,q,b,s : Double;
a         : array[0..4] of Double;
begin
fGain:=1;
if fFilterType=ftLowPass
then s:=1
else s:=-1;
// calculating coefficients:
k:=(4.0*fGain-3.0)/(fGain+1.0);
p:=1.0-0.25*k;
p:=p*p;

if fFilterType=ftLowPass
then b:=1.0/(tan(pi*fFreq/fSR)*(1.0+p))
else b:=tan(pi*fFreq/fSR)/(1.0+p);
p:=1.0+b;
q:=s*(1.0-b);

a := 1.0/(  k+p*p*p*p);
a := 4.0*(s*k+p*p*p*q);
a := 6.0*(  k+p*p*q*q);
a := 4.0*(s*k+p*q*q*q);
a :=     (  k+q*q*q*q);
p    := a*(k+1.0);

fCoeffs:=p;
fCoeffs:=4.0*p*s;
fCoeffs:=6.0*p;
fCoeffs:=4.0*p*s;
fCoeffs:=p;
fCoeffs:=-a*a;
fCoeffs:=-a*a;
fCoeffs:=-a*a;
fCoeffs:=-a*a;
end;

function T4Pole.Process(s:single):single;
begin
Result:=fCoeffs*s+d;
d:=fCoeffs*s+fCoeffs*Result+d;
d:=fCoeffs*s+fCoeffs*Result+d;
d:=fCoeffs*s+fCoeffs*Result+d;
d:=fCoeffs*s+fCoeffs*Result;
end;

end.
```
```so bad that this filter is so unstable. i tested it and is has a really nice sound. but frequencies below 200 hz are not possible. :-(
```

### Bass Booster¶

• Author or source: Johny Dupej
• Type: LP and SUM
• Created: 2006-08-11 12:47:34
notes
```This function adds a low-passed signal to the original signal. The low-pass has a quite
wide response.

Params:
selectivity - frequency response of the LP (higher value gives a steeper one) [70.0 to
140.0 sounds good]
ratio - how much of the filtered signal is mixed to the original
gain2 - adjusts the final volume to handle cut-offs (might be good to set dynamically)
```
code
 ``` 1 2 3 4 5 6 7 8 9 10 11 12``` ```#define saturate(x) __min(__max(-1.0,x),1.0) float BassBoosta(float sample) { static float selectivity, gain1, gain2, ratio, cap; gain1 = 1.0/(selectivity + 1.0); cap= (sample + cap*selectivity )*gain1; sample = saturate((sample + cap*ratio)*gain2); return sample; } ```

```Implementation of the RBJ cookbook, in C.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821 1822 1823 1824 1825 1826 1827 1828 1829 1830 1831 1832 1833 1834 1835 1836 1837 1838 1839 1840 1841 1842 1843 1844 1845 1846 1847 1848 1849 1850 1851 1852 1853 1854 1855 1856 1857 1858 1859 1860 1861 1862 1863 1864 1865 1866 1867 1868 1869 1870 1871 1872 1873 1874 1875 1876 1877 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 2020 2021 2022 2023 2024 2025 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045 2046 2047 2048 2049 2050 2051 2052 2053 2054 2055 2056 2057 2058 2059 2060 2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 2112 2113 2114 2115 2116 2117 2118 2119 2120 2121 2122 2123 2124 2125 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142 2143 2144 2145 2146 2147 2148 2149 2150 2151 2152 2153 2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187 2188 2189 2190 2191 2192 2193 2194 2195 2196 2197 2198 2199 2200 2201 2202 2203 2204 2205 2206 2207 2208 2209 2210 2211 2212 2213 2214 2215 2216 2217 2218 2219 2220 2221 2222 2223 2224 2225 2226 2227 2228 2229 2230 2231 2232 2233 2234 2235 2236 2237 2238 2239 2240 2241 2242 2243 2244 2245 2246 2247 2248 2249 2250 2251 2252 2253 2254 2255 2256 2257 2258 2259 2260 2261 2262 2263 2264 2265 2266 2267 2268 2269 2270 2271 2272 2273 2274 2275 2276 2277 2278 2279 2280 2281 2282 2283 2284 2285 2286 2287 2288 2289 2290 2291 2292 2293 2294 2295 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 2309 2310 2311 2312 2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 2325 2326 2327 2328 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 2339 2340 2341 2342 2343 2344 2345 2346 2347 2348 2349 2350 2351 2352 2353 2354 2355 2356 2357 2358 2359 2360 2361 2362 2363 2364 2365 2366 2367 2368 2369 2370 2371 2372 2373 2374 2375 2376 2377 2378 2379 2380 2381 2382 2383 2384 2385 2386 2387 2388 2389 2390 2391 2392 2393 2394 2395 2396 2397 2398 2399 2400 2401 2402 2403 2404 2405 2406 2407 2408 2409 2410 2411 2412 2413 2414 2415 2416 2417 2418 2419 2420 2421 2422 2423 2424 2425 2426 2427 2428 2429 2430 2431 2432 2433 2434 2435 2436 2437 2438 2439 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453 2454 2455 2456 2457 2458 2459 2460 2461 2462 2463 2464 2465 2466 2467 2468 2469 2470 2471 2472 2473 2474 2475 2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 2486 2487 2488 2489 2490 2491 2492 2493 2494 2495 2496 2497 2498 2499 2500 2501 2502 2503 2504 2505 2506``` ```/* "A Collection of Useful C++ Classes for Digital Signal Processing" By Vincent Falco Please direct all comments to either the music-dsp mailing list or the DSP and Plug-in Development forum: http://music.columbia.edu/cmc/music-dsp/ http://www.kvraudio.com/forum/viewforum.php?f=33 http://www.kvraudio.com/forum/ Support is provided for performing N-order Dsp floating point filter operations on M-channel data with a caller specified floating point type. The implementation breaks a high order IIR filter down into a series of cascaded second order stages. Tests conclude that numerical stability is maintained even at higher orders. For example the Butterworth low pass filter is stable at up to 53 poles. Processing functions are provided to use either Direct Form I or Direct Form II of the filter transfer function. Direct Form II is slightly faster but can cause discontinuities in the output if filter parameters are changed during processing. Direct Form I is slightly slower, but maintains fidelity even when parameters are changed during processing. To support fast parameter changes, filters provide two functions for adjusting parameters. A high accuracy Setup() function, and a faster form called SetupFast() that uses approximations for trigonometric functions. The approximations work quite well and should be suitable for most applications. Channels are stored in an interleaved format with M samples per frame arranged contiguously. A single class instance can process all M channels simultaneously in an efficient manner. A 'skip' parameter causes the processing function to advance by skip additional samples in the destination buffer in between every frame. Through manipulation of the skip paramter it is possible to exclude channels from processing (for example, only processing the left half of stereo interleaved data). For multichannel data which is not interleaved, it will be necessary to instantiate multiple instance of the filter and set skip=0. There are a few other utility classes and functions included that may prove useful. Classes: Complex CascadeStages Biquad BiquadLowPass BiquadHighPass BiquadBandPass1 BiquadBandPass2 BiquadBandStop BiquadAllPass BiquadPeakEq BiquadLowShelf BiquadHighShelf PoleFilter Butterworth ButterLowPass ButterHighPass ButterBandPass ButterBandStop Chebyshev1 Cheby1LowPass Cheby1HighPass Cheby1BandPass Cheby1BandStop Chebyshev2 Cheby2LowPass Cheby2HighPass Cheby2BandPass Cheby2BandStop EnvelopeFollower AutoLimiter Functions: zero() copy() mix() scale() interleave() deinterleave() Order for PoleFilter derived classes is specified in the number of poles, except for band pass and band stop filters, for which the number of pole pairs is specified. For some filters there are two versions of Setup(), the one called SetupFast() uses approximations to trigonometric functions for speed. This is an option if you are doing frequent parameter changes to the filter. There is an example function at the bottom that shows how to use the classes. Filter ideas are based on a java applet (http://www.falstad.com/dfilter/) developed by Paul Falstad. All of this code was written by the author Vincent Falco except where marked. -------------------------------------------------------------------------------- License: MIT License (http://www.opensource.org/licenses/mit-license.php) Copyright (c) 2009 by Vincent Falco Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * */ /* To Do: - Shelving, peak, all-pass for Butterworth, Chebyshev, and Elliptic. The Biquads have these versions and I would like the others to have them as well. It would also be super awesome if higher order filters could have a "Q" parameter for resonance but I'm not expecting miracles, it would require redistributing the poles and zeroes. But if there's a research paper or code out there...I could incorporate it. - Denormal testing and fixing I'd like to know if denormals are a problem. And if so, it would be nice to have a small function that can reproduce the denormal problem. This way I can test the fix under all conditions. I will include the function as a "unit test" object in the header file so anyone can verify its correctness. But I'm a little lost. - Optimized template specializations for stages=1, channels={1,2} There are some pretty obvious optimizations I am saving for "the end". I don't want to do them until the code is finalized. - Optimized template specializations for SSE, other special instructions - Optimized trigonometric functions for fast parameter changes - Elliptic curve based filter coefficients - More utility functions for manipulating sample buffers - Need fast version of pow( 10, x ) */ #ifndef __DSP_FILTER__ #define __DSP_FILTER__ #include #include #include #include #include //#define DSP_USE_STD_COMPLEX #ifdef DSP_USE_STD_COMPLEX #include #endif #define DSP_SSE3_OPTIMIZED #ifdef DSP_SSE3_OPTIMIZED //#include //#include #include #endif namespace Dsp { //-------------------------------------------------------------------------- // WARNING: Here there be templates //-------------------------------------------------------------------------- //-------------------------------------------------------------------------- // // Configuration // //-------------------------------------------------------------------------- // Regardless of the type of sample that the filter operates on (e.g. // float or double), all calculations are performed using double (or // better) for stability and accuracy. This controls the underlying // type used for calculations: typedef double CalcT; typedef int Int32; // Must be 32 bits typedef __int64 Int64; // Must be 64 bits // This is used to prevent denormalization. const CalcT vsa=1.0 / 4294967295.0; // for CalcT as float // These constants are so important, I made my own copy. If you improve // the resolution of CalcT be sure to add more significant digits to these. const CalcT kPi =3.1415926535897932384626433832795; const CalcT kPi_2 =1.57079632679489661923; const CalcT kLn2 =0.693147180559945309417; const CalcT kLn10 =2.30258509299404568402; //-------------------------------------------------------------------------- // Some functions missing from template inline T acosh( T x ) { return log( x+::sqrt(x*x-1) ); } //-------------------------------------------------------------------------- // // Fast Trigonometric Functions // //-------------------------------------------------------------------------- // Three approximations for both sine and cosine at a given angle. // The faster the routine, the larger the error. // From http://lab.polygonal.de/2007/07/18/fast-and-accurate-sinecosine-approximation/ // Tuned for maximum pole stability. r must be in the range 0..pi // This one breaks down considerably at the higher angles. It is // included only for educational purposes. inline void fastestsincos( CalcT r, CalcT *sn, CalcT *cs ) { const CalcT c=0.70710678118654752440; // sqrt(2)/2 CalcT v=(2-4*c)*r*r+c; if(r kPi) x -= 2*kPi; //compute sine if (x < 0) *sn = 1.27323954 * x + 0.405284735 * x * x; else *sn = 1.27323954 * x - 0.405284735 * x * x; //compute cosine: sin(x + PI/2) = cos(x) x += kPi_2; if (x > kPi ) x -= 2*kPi; if (x < 0) *cs = 1.27323954 * x + 0.405284735 * x * x; else *cs = 1.27323954 * x - 0.405284735 * x * x; } // Slower than ::fastersincos() but still faster than // sin(), and with the best accuracy of these routines. inline void fastsincos( CalcT x, CalcT *sn, CalcT *cs ) { CalcT s, c; //always wrap input angle to -PI..PI if (x < -kPi) x += 2*kPi; else if (x > kPi) x -= 2*kPi; //compute sine if (x < 0) { s = 1.27323954 * x + .405284735 * x * x; if (s < 0) s = .225 * (s * -s - s) + s; else s = .225 * (s * s - s) + s; } else { s = 1.27323954 * x - 0.405284735 * x * x; if (s < 0) s = .225 * (s * -s - s) + s; else s = .225 * (s * s - s) + s; } *sn=s; //compute cosine: sin(x + PI/2) = cos(x) x += kPi_2; if (x > kPi ) x -= 2*kPi; if (x < 0) { c = 1.27323954 * x + 0.405284735 * x * x; if (c < 0) c = .225 * (c * -c - c) + c; else c = .225 * (c * c - c) + c; } else { c = 1.27323954 * x - 0.405284735 * x * x; if (c < 0) c = .225 * (c * -c - c) + c; else c = .225 * (c * c - c) + c; } *cs=c; } // Faster approximations to sqrt() // From http://ilab.usc.edu/wiki/index.php/Fast_Square_Root // The faster the routine, the more error in the approximation. // Log Base 2 Approximation // 5 times faster than sqrt() inline float fastsqrt1( float x ) { union { Int32 i; float x; } u; u.x = x; u.i = (Int32(1)<<29) + (u.i >> 1) - (Int32(1)<<22); return u.x; } inline double fastsqrt1( double x ) { union { Int64 i; double x; } u; u.x = x; u.i = (Int64(1)<<61) + (u.i >> 1) - (Int64(1)<<51); return u.x; } // Log Base 2 Approximation with one extra Babylonian Step // 2 times faster than sqrt() inline float fastsqrt2( float x ) { float v=fastsqrt1( x ); v = 0.5f * (v + x/v); // One Babylonian step return v; } inline double fastsqrt2(const double x) { double v=fastsqrt1( x ); v = 0.5f * (v + x/v); // One Babylonian step return v; } // Log Base 2 Approximation with two extra Babylonian Steps // 50% faster than sqrt() inline float fastsqrt3( float x ) { float v=fastsqrt1( x ); v = v + x/v; v = 0.25f* v + x/v; // Two Babylonian steps return v; } inline double fastsqrt3(const double x) { double v=fastsqrt1( x ); v = v + x/v; v = 0.25 * v + x/v; // Two Babylonian steps return v; } //-------------------------------------------------------------------------- // // Complex // //-------------------------------------------------------------------------- #ifdef DSP_USE_STD_COMPLEX template inline std::complex polar( const T &m, const T &a ) { return std::polar( m, a ); } template inline T norm( const std::complex &c ) { return std::norm( c ); } template inline T abs( const std::complex &c ) { return std::abs(c); } template inline std::complex addmul( const std::complex &c, T v, const std::complex &c1 ) { return std::complex( c.real()+v*c1.real(), c.imag()+v*c1.imag() ); } template inline T arg( const std::complex &c ) { return std::arg( c ); } template inline std::complex recip( const std::complex &c ) { T n=1.0/Dsp::norm(c); return std::complex( n*c.real(), n*c.imag() ); } template inline std::complex sqrt( const std::complex &c ) { return std::sqrt( c ); } typedef std::complex Complex; #else //-------------------------------------------------------------------------- // "Its always good to have a few extra wheels in case one goes flat." template struct ComplexT { ComplexT(); ComplexT( T r_, T i_=0 ); template ComplexT( const ComplexT &c ); T imag ( void ) const; T real ( void ) const; ComplexT & neg ( void ); ComplexT & conj ( void ); template ComplexT & add ( const ComplexT &c ); template ComplexT & sub ( const ComplexT &c ); template ComplexT & mul ( const ComplexT &c ); template ComplexT & div ( const ComplexT &c ); template ComplexT & addmul ( T v, const ComplexT &c ); ComplexT operator- ( void ) const; ComplexT operator+ ( T v ) const; ComplexT operator- ( T v ) const; ComplexT operator* ( T v ) const; ComplexT operator/ ( T v ) const; ComplexT & operator+= ( T v ); ComplexT & operator-= ( T v ); ComplexT & operator*= ( T v ); ComplexT & operator/= ( T v ); template ComplexT operator+ ( const ComplexT &c ) const; template ComplexT operator- ( const ComplexT &c ) const; template ComplexT operator* ( const ComplexT &c ) const; template ComplexT operator/ ( const ComplexT &c ) const; template ComplexT & operator+= ( const ComplexT &c ); template ComplexT & operator-= ( const ComplexT &c ); template ComplexT & operator*= ( const ComplexT &c ); template ComplexT & operator/= ( const ComplexT &c ); private: ComplexT & add ( T v ); ComplexT & sub ( T v ); ComplexT & mul ( T c, T d ); ComplexT & mul ( T v ); ComplexT & div ( T v ); T r; T i; }; //-------------------------------------------------------------------------- template inline ComplexT::ComplexT() { } template inline ComplexT::ComplexT( T r_, T i_ ) { r=r_; i=i_; } template template inline ComplexT::ComplexT( const ComplexT &c ) { r=c.r; i=c.i; } template inline T ComplexT::imag( void ) const { return i; } template inline T ComplexT::real( void ) const { return r; } template inline ComplexT &ComplexT::neg( void ) { r=-r; i=-i; return *this; } template inline ComplexT &ComplexT::conj( void ) { i=-i; return *this; } template inline ComplexT &ComplexT::add( T v ) { r+=v; return *this; } template inline ComplexT &ComplexT::sub( T v ) { r-=v; return *this; } template inline ComplexT &ComplexT::mul( T c, T d ) { T ac=r*c; T bd=i*d; // must do i first i=(r+i)*(c+d)-(ac+bd); r=ac-bd; return *this; } template inline ComplexT &ComplexT::mul( T v ) { r*=v; i*=v; return *this; } template inline ComplexT &ComplexT::div( T v ) { r/=v; i/=v; return *this; } template template inline ComplexT &ComplexT::add( const ComplexT &c ) { r+=c.r; i+=c.i; return *this; } template template inline ComplexT &ComplexT::sub( const ComplexT &c ) { r-=c.r; i-=c.i; return *this; } template template inline ComplexT &ComplexT::mul( const ComplexT &c ) { return mul( c.r, c.i ); } template template inline ComplexT &ComplexT::div( const ComplexT &c ) { T s=1.0/norm(c); return mul( c.r*s, -c.i*s ); } template inline ComplexT ComplexT::operator-( void ) const { return ComplexT(*this).neg(); } template inline ComplexT ComplexT::operator+( T v ) const { return ComplexT(*this).add( v ); } template inline ComplexT ComplexT::operator-( T v ) const { return ComplexT(*this).sub( v ); } template inline ComplexT ComplexT::operator*( T v ) const { return ComplexT(*this).mul( v ); } template inline ComplexT ComplexT::operator/( T v ) const { return ComplexT(*this).div( v ); } template inline ComplexT &ComplexT::operator+=( T v ) { return add( v ); } template inline ComplexT &ComplexT::operator-=( T v ) { return sub( v ); } template inline ComplexT &ComplexT::operator*=( T v ) { return mul( v ); } template inline ComplexT &ComplexT::operator/=( T v ) { return div( v ); } template template inline ComplexT ComplexT::operator+( const ComplexT &c ) const { return ComplexT(*this).add(c); } template template inline ComplexT ComplexT::operator-( const ComplexT &c ) const { return ComplexT(*this).sub(c); } template template inline ComplexT ComplexT::operator*( const ComplexT &c ) const { return ComplexT(*this).mul(c); } template template inline ComplexT ComplexT::operator/( const ComplexT &c ) const { return ComplexT(*this).div(c); } template template inline ComplexT &ComplexT::operator+=( const ComplexT &c ) { return add( c ); } template template inline ComplexT &ComplexT::operator-=( const ComplexT &c ) { return sub( c ); } template template inline ComplexT &ComplexT::operator*=( const ComplexT &c ) { return mul( c ); } template template inline ComplexT &ComplexT::operator/=( const ComplexT &c ) { return div( c ); } //-------------------------------------------------------------------------- template inline ComplexT polar( const T &m, const T &a ) { return ComplexT( m*cos(a), m*sin(a) ); } template inline T norm( const ComplexT &c ) { return c.real()*c.real()+c.imag()*c.imag(); } template inline T abs( const ComplexT &c ) { return ::sqrt( c.real()*c.real()+c.imag()*c.imag() ); } template inline ComplexT addmul( const ComplexT &c, T v, const ComplexT &c1 ) { return ComplexT( c.real()+v*c1.real(), c.imag()+v*c1.imag() ); } template inline T arg( const ComplexT &c ) { return atan2( c.imag(), c.real() ); } template inline ComplexT recip( const ComplexT &c ) { T n=1.0/norm(c); return ComplexT( n*c.real(), -n*c.imag() ); } template inline ComplexT sqrt( const ComplexT &c ) { return polar( ::sqrt(abs(c)), arg(c)*0.5 ); } //-------------------------------------------------------------------------- typedef ComplexT Complex; #endif //-------------------------------------------------------------------------- // // Numerical Analysis // //-------------------------------------------------------------------------- // Implementation of Brent's Method provided by // John D. Cook (http://www.johndcook.com/) // The return value of Minimize is the minimum of the function f. // The location where f takes its minimum is returned in the variable minLoc. // Notation and implementation based on Chapter 5 of Richard Brent's book // "Algorithms for Minimization Without Derivatives". template CalcT BrentMinimize ( TFunction& f, // [in] objective function to minimize CalcT leftEnd, // [in] smaller value of bracketing interval CalcT rightEnd, // [in] larger value of bracketing interval CalcT epsilon, // [in] stopping tolerance CalcT& minLoc // [out] location of minimum ) { CalcT d, e, m, p, q, r, tol, t2, u, v, w, fu, fv, fw, fx; static const CalcT c = 0.5*(3.0 - ::sqrt(5.0)); static const CalcT SQRT_DBL_EPSILON = ::sqrt(DBL_EPSILON); CalcT& a = leftEnd; CalcT& b = rightEnd; CalcT& x = minLoc; v = w = x = a + c*(b - a); d = e = 0.0; fv = fw = fx = f(x); int counter = 0; loop: counter++; m = 0.5*(a + b); tol = SQRT_DBL_EPSILON*::fabs(x) + epsilon; t2 = 2.0*tol; // Check stopping criteria if (::fabs(x - m) > t2 - 0.5*(b - a)) { p = q = r = 0.0; if (::fabs(e) > tol) { // fit parabola r = (x - w)*(fx - fv); q = (x - v)*(fx - fw); p = (x - v)*q - (x - w)*r; q = 2.0*(q - r); (q > 0.0) ? p = -p : q = -q; r = e; e = d; } if (::fabs(p) < ::fabs(0.5*q*r) && p < q*(a - x) && p < q*(b - x)) { // A parabolic interpolation step d = p/q; u = x + d; // f must not be evaluated too close to a or b if (u - a < t2 || b - u < t2) d = (x < m) ? tol : -tol; } else { // A golden section step e = (x < m) ? b : a; e -= x; d = c*e; } // f must not be evaluated too close to x if (::fabs(d) >= tol) u = x + d; else if (d > 0.0) u = x + tol; else u = x - tol; fu = f(u); // Update a, b, v, w, and x if (fu <= fx) { (u < x) ? b = x : a = x; v = w; fv = fw; w = x; fw = fx; x = u; fx = fu; } else { (u < x) ? a = u : b = u; if (fu <= fw || w == x) { v = w; fv = fw; w = u; fw = fu; } else if (fu <= fv || v == x || v == w) { v = u; fv = fu; } } goto loop; // Yes, the dreaded goto statement. But the code // here is faithful to Brent's orginal pseudocode. } return fx; } //-------------------------------------------------------------------------- // // Infinite Impulse Response Filters // //-------------------------------------------------------------------------- // IIR filter implementation using multiple second-order stages. class CascadeFilter { public: // Process data in place using Direct Form I // skip is added after each frame. // Direct Form I is more suitable when the filter parameters // are changed often. However, it is slightly slower. template void ProcessI( size_t frames, T *dest, int skip=0 ); // Process data in place using Direct Form II // skip is added after each frame. // Direct Form II is slightly faster than Direct Form I, // but changing filter parameters on stream can result // in discontinuities in the output. It is best suited // for a filter whose parameters are set only once. template void ProcessII( size_t frames, T *dest, int skip=0 ); // Convenience function that just calls ProcessI. // Feel free to change the implementation. template void Process( size_t frames, T *dest, int skip=0 ); // Determine response at angular frequency (0<=w<=kPi) Complex Response( CalcT w ); // Clear the history buffer. void Clear( void ); protected: struct Hist; struct Stage; // for m_nchan==2 #ifdef DSP_SSE3_OPTIMIZED template void ProcessISSEStageStereo( size_t frames, T *dest, Stage *stage, Hist *h, int skip ); template void ProcessISSEStereo( size_t frames, T *dest, int skip ); #endif protected: void Reset ( void ); void Normalize ( CalcT scale ); void SetAStage ( CalcT x1, CalcT x2 ); void SetBStage ( CalcT x0, CalcT x1, CalcT x2 ); void SetStage ( CalcT a1, CalcT a2, CalcT b0, CalcT b1, CalcT b2 ); protected: struct Hist { CalcT v; }; struct Stage { CalcT a; // a unused CalcT b; void Reset( void ); }; struct ResponseFunctor { CascadeFilter *m_filter; CalcT operator()( CalcT w ); ResponseFunctor( CascadeFilter *filter ); }; int m_nchan; int m_nstage; Stage * m_stagep; Hist * m_histp; }; //-------------------------------------------------------------------------- template void CascadeFilter::ProcessI( size_t frames, T *dest, int skip ) { #ifdef DSP_SSE3_OPTIMIZED if( m_nchan==2 ) ProcessISSEStereo( frames, dest, skip ); else #endif while( frames-- ) { Hist *h=m_histp; for( int j=m_nchan;j;j-- ) { CalcT in=CalcT(*dest); Stage *s=m_stagep; for( int i=m_nstage;i;i--,h++,s++ ) { CalcT out; out=s->b*in + s->b*h->v + s->b*h->v + s->a*h->v + s->a*h->v; h->v=h->v; h->v=in; h->v=h->v; h->v=out; in=out; } *dest++=T(in); } dest+=skip; } } // A good compiler already produces code that is optimized even for // the general case. The only way to make it go faster is to // to implement it in assembler or special instructions. Like this: #ifdef DSP_SSE3_OPTIMIZED // ALL SSE OPTIMIZATIONS ASSUME CalcT as double template inline void CascadeFilter::ProcessISSEStageStereo( size_t frames, T *dest, Stage *s, Hist *h, int skip ) { assert( m_nchan==2 ); #if 1 CalcT b0=s->b; __m128d m0=_mm_loadu_pd( &s->a ); // a1 , a2 __m128d m1=_mm_loadu_pd( &s->b ); // b1 , b2 __m128d m2=_mm_loadu_pd( &h.v ); // h->v , h->v __m128d m3=_mm_loadu_pd( &h.v ); // h->v , h->v __m128d m4=_mm_loadu_pd( &h.v ); // h->v , h->v __m128d m5=_mm_loadu_pd( &h.v ); // h->v , h->v while( frames-- ) { CalcT in, b0in, out; __m128d m6; __m128d m7; in=CalcT(*dest); b0in=b0*in; m6=_mm_mul_pd ( m1, m2 ); // b1*h->v , b2*h->v m7=_mm_mul_pd ( m0, m3 ); // a1*h->v , a2*h->v m6=_mm_add_pd ( m6, m7 ); // b1*h->v + a1*h->v, b2*h->v + a2*h->v m7=_mm_load_sd( &b0in ); // b0*in , 0 m6=_mm_add_sd ( m6, m7 ); // b1*h->v + a1*h->v + in*b0 , b2*h->v + a2*h->v + 0 m6=_mm_hadd_pd( m6, m7 ); // b1*h->v + a1*h->v + in*b0 + b2*h->v + a2*h->v, in*b0 _mm_store_sd( &out, m6 ); m6=_mm_loadh_pd( m6, &in ); // out , in m2=_mm_shuffle_pd( m6, m2, _MM_SHUFFLE2( 0, 1 ) ); // h->v=in , h->v=h->v m3=_mm_shuffle_pd( m6, m3, _MM_SHUFFLE2( 0, 0 ) ); // h->v=out, h->v=h->v *dest++=T(out); in=CalcT(*dest); b0in=b0*in; m6=_mm_mul_pd ( m1, m4 ); // b1*h->v , b2*h->v m7=_mm_mul_pd ( m0, m5 ); // a1*h->v , a2*h->v m6=_mm_add_pd ( m6, m7 ); // b1*h->v + a1*h->v, b2*h->v + a2*h->v m7=_mm_load_sd( &b0in ); // b0*in , 0 m6=_mm_add_sd ( m6, m7 ); // b1*h->v + a1*h->v + in*b0 , b2*h->v + a2*h->v + 0 m6=_mm_hadd_pd( m6, m7 ); // b1*h->v + a1*h->v + in*b0 + b2*h->v + a2*h->v, in*b0 _mm_store_sd( &out, m6 ); m6=_mm_loadh_pd( m6, &in ); // out , in m4=_mm_shuffle_pd( m6, m4, _MM_SHUFFLE2( 0, 1 ) ); // h->v=in , h->v=h->v m5=_mm_shuffle_pd( m6, m5, _MM_SHUFFLE2( 0, 0 ) ); // h->v=out, h->v=h->v *dest++=T(out); dest+=skip; } // move history from registers back to state _mm_storeu_pd( &h.v, m2 ); _mm_storeu_pd( &h.v, m3 ); _mm_storeu_pd( &h.v, m4 ); _mm_storeu_pd( &h.v, m5 ); #else // Template-specialized version from which the assembly was modeled CalcT a1=s->a; CalcT a2=s->a; CalcT b0=s->b; CalcT b1=s->b; CalcT b2=s->b; while( frames-- ) { CalcT in, out; in=CalcT(*dest); out=b0*in+b1*h.v+b2*h.v +a1*h.v+a2*h.v; h.v=h.v; h.v=in; h.v=h.v; h.v=out; in=out; *dest++=T(in); in=CalcT(*dest); out=b0*in+b1*h.v+b2*h.v +a1*h.v+a2*h.v; h.v=h.v; h.v=in; h.v=h.v; h.v=out; in=out; *dest++=T(in); dest+=skip; } #endif } // Note there could be a loss of accuracy here. Unlike the original version // of Process...() we are applying each stage to all of the input data. // Since the underlying type T could be float, the results from this function // may be different than the unoptimized version. However, it is much faster. template void CascadeFilter::ProcessISSEStereo( size_t frames, T *dest, int skip ) { assert( m_nchan==2 ); Stage *s=m_stagep; Hist *h=m_histp; for( int i=m_nstage;i;i--,h+=2,s++ ) { ProcessISSEStageStereo( frames, dest, s, h, skip ); } } #endif template void CascadeFilter::ProcessII( size_t frames, T *dest, int skip ) { while( frames-- ) { Hist *h=m_histp; for( int j=m_nchan;j;j-- ) { CalcT in=CalcT(*dest); Stage *s=m_stagep; for( int i=m_nstage;i;i--,h++,s++ ) { CalcT d2=h->v=h->v; CalcT d1=h->v=h->v; CalcT d0=h->v= in+s->a*d1 + s->a*d2; in=s->b*d0 + s->b*d1 + s->b*d2; } *dest++=T(in); } dest+=skip; } } template inline void CascadeFilter::Process( size_t frames, T *dest, int skip ) { ProcessI( frames, dest, skip ); } inline Complex CascadeFilter::Response( CalcT w ) { Complex ch( 1 ); Complex cbot( 1 ); Complex czn1=polar( 1., -w ); Complex czn2=polar( 1., -2*w ); Stage *s=m_stagep; for( int i=m_nstage;i;i-- ) { Complex ct( s->b ); Complex cb( 1 ); ct=addmul( ct, s->b, czn1 ); cb=addmul( cb, -s->a, czn1 ); ct=addmul( ct, s->b, czn2 ); cb=addmul( cb, -s->a, czn2 ); ch*=ct; cbot*=cb; s++; } return ch/cbot; } inline void CascadeFilter::Clear( void ) { ::memset( m_histp, 0, m_nstage*m_nchan*sizeof(m_histp) ); } inline void CascadeFilter::Stage::Reset( void ) { a=0; a=0; b=1; b=0; b=0; } inline void CascadeFilter::Reset( void ) { Stage *s=m_stagep; for( int i=m_nstage;i;i--,s++ ) s->Reset(); } // Apply scale factor to stage coefficients. inline void CascadeFilter::Normalize( CalcT scale ) { // We are throwing the normalization into the first // stage. In theory it might be nice to spread it around // to preserve numerical accuracy. Stage *s=m_stagep; s->b*=scale; s->b*=scale; s->b*=scale; } inline void CascadeFilter::SetAStage( CalcT x1, CalcT x2 ) { Stage *s=m_stagep; for( int i=m_nstage;i;i-- ) { if( s->a==0 && s->a==0 ) { s->a=x1; s->a=x2; s=0; break; } if( s->a==0 && x2==0 ) { s->a=-s->a*x1; s->a+=x1; s=0; break; } s++; } assert( s==0 ); } inline void CascadeFilter::SetBStage( CalcT x0, CalcT x1, CalcT x2 ) { Stage *s=m_stagep; for( int i=m_nstage;i;i-- ) { if( s->b==0 && s->b==0 ) { s->b=x0; s->b=x1; s->b=x2; s=0; break; } if( s->b==0 && x2==0 ) { // (b0 + z b1)(x0 + z x1) = (b0 x0 + (b1 x0+b0 x1) z + b1 x1 z^2) s->b=s->b*x1; s->b=s->b*x0+s->b*x1; s->b*=x0; s=0; break; } s++; } assert( s==0 ); } // Optimized version for Biquads inline void CascadeFilter::SetStage( CalcT a1, CalcT a2, CalcT b0, CalcT b1, CalcT b2 ) { assert( m_nstage==1 ); Stage *s=&m_stagep; s->a=a1; s->a=a2; s->b=b0; s->b=b1; s->b=b2; } inline CalcT CascadeFilter::ResponseFunctor::operator()( CalcT w ) { return -Dsp::abs(m_filter->Response( w )); } inline CascadeFilter::ResponseFunctor::ResponseFunctor( CascadeFilter *filter ) { m_filter=filter; } //-------------------------------------------------------------------------- template class CascadeStages : public CascadeFilter { public: CascadeStages(); private: Hist m_hist [stages*channels]; Stage m_stage [stages]; }; //-------------------------------------------------------------------------- template CascadeStages::CascadeStages( void ) { m_nchan=channels; m_nstage=stages; m_stagep=m_stage; m_histp=m_hist; Clear(); } //-------------------------------------------------------------------------- // // Biquad Second Order IIR Filters // //-------------------------------------------------------------------------- // Formulas from http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt template class Biquad : public CascadeStages<1, channels> { protected: void Setup( const CalcT a, const CalcT b ); }; //-------------------------------------------------------------------------- template inline void Biquad::Setup( const CalcT a, const CalcT b ) { Reset(); // transform Biquad coefficients CalcT ra0=1/a; SetAStage( -a*ra0, -a*ra0 ); SetBStage( b*ra0, b*ra0, b*ra0 ); } //-------------------------------------------------------------------------- template class BiquadLowPass : public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template inline void BiquadLowPass::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT a0 = 1 / ( 1 + alph ); CalcT b1 = 1 - cs; CalcT b0 = a0 * b1 * 0.5; CalcT a1 = 2 * cs; CalcT a2 = alph - 1; SetStage( a1*a0, a2*a0, b0, b1*a0, b0 ); } template void BiquadLowPass::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadLowPass::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- template class BiquadHighPass : public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template inline void BiquadHighPass::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT a0 = -1 / ( 1 + alph ); CalcT b1 = -( 1 + cs ); CalcT b0 = a0 * b1 * -0.5; CalcT a1 = -2 * cs; CalcT a2 = 1 - alph; SetStage( a1*a0, a2*a0, b0, b1*a0, b0 ); } template void BiquadHighPass::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadHighPass::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- // Constant skirt gain, peak gain=Q template class BiquadBandPass1 : public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template inline void BiquadBandPass1::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT a0 = -1 / ( 1 + alph ); CalcT b0 = a0 * ( sn * -0.5 ); CalcT a1 = -2 * cs; CalcT a2 = 1 - alph; SetStage( a1*a0, a2*a0, b0, 0, -b0 ); } template void BiquadBandPass1::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadBandPass1::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- // Constant 0dB peak gain template class BiquadBandPass2 : public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template inline void BiquadBandPass2::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT b0 = -alph; CalcT b2 = alph; CalcT a0 = -1 / ( 1 + alph ); CalcT a1 = -2 * cs; CalcT a2 = 1 - alph; SetStage( a1*a0, a2*a0, b0*a0, 0, b2*a0 ); } template void BiquadBandPass2::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadBandPass2::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- template class BiquadBandStop : public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template inline void BiquadBandStop::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT a0 = 1 / ( 1 + alph ); CalcT b1 = a0 * ( -2 * cs ); CalcT a2 = alph - 1; SetStage( -b1, a2*a0, a0, b1, a0 ); } template void BiquadBandStop::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadBandStop::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- template class BiquadAllPass: public Biquad { public: void Setup ( CalcT normFreq, CalcT q ); void SetupFast ( CalcT normFreq, CalcT q ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT q ); }; //-------------------------------------------------------------------------- template void BiquadAllPass::SetupCommon( CalcT sn, CalcT cs, CalcT q ) { CalcT alph = sn / ( 2 * q ); CalcT b2 = 1 + alph; CalcT a0 = 1 / b2; CalcT b0 =( 1 - alph ) * a0; CalcT b1 = -2 * cs * a0; SetStage( -b1, -b0, b0, b1, b2*a0 ); } template void BiquadAllPass::Setup( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); SetupCommon( sn, cs, q ); } template void BiquadAllPass::SetupFast( CalcT normFreq, CalcT q ) { CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); SetupCommon( sn, cs, q ); } //-------------------------------------------------------------------------- template class BiquadPeakEq: public Biquad { public: void Setup ( CalcT normFreq, CalcT dB, CalcT bandWidth ); void SetupFast ( CalcT normFreq, CalcT dB, CalcT bandWidth ); protected: void SetupCommon ( CalcT sn, CalcT cs, CalcT alph, CalcT A ); }; //-------------------------------------------------------------------------- template inline void BiquadPeakEq::SetupCommon( CalcT sn, CalcT cs, CalcT alph, CalcT A ) { CalcT t=alph*A; CalcT b0 = 1 - t; CalcT b2 = 1 + t; t=alph/A; CalcT a0 = 1 / ( 1 + t ); CalcT a2 = t - 1; CalcT b1 = a0 * ( -2 * cs ); CalcT a1 = -b1; SetStage( a1, a2*a0, b0*a0, b1, b2*a0 ); } template void BiquadPeakEq::Setup( CalcT normFreq, CalcT dB, CalcT bandWidth ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); CalcT alph = sn * sinh( kLn2/2 * bandWidth * w0/sn ); SetupCommon( sn, cs, alph, A ); } template void BiquadPeakEq::SetupFast( CalcT normFreq, CalcT dB, CalcT bandWidth ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); CalcT alph = sn * sinh( kLn2/2 * bandWidth * w0/sn ); SetupCommon( sn, cs, alph, A ); } //-------------------------------------------------------------------------- template class BiquadLowShelf : public Biquad { public: void Setup ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 ); void SetupFast ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 ); protected: void SetupCommon ( CalcT cs, CalcT A, CalcT sa ); }; //-------------------------------------------------------------------------- template inline void BiquadLowShelf::SetupCommon( CalcT cs, CalcT A, CalcT sa ) { CalcT An = A-1; CalcT Ap = A+1; CalcT Ancs = An*cs; CalcT Apcs = Ap*cs; CalcT b0 = A * (Ap - Ancs + sa ); CalcT b2 = A * (Ap - Ancs - sa ); CalcT b1 = 2 * A * (An - Apcs); CalcT a2 = sa - (Ap + Ancs); CalcT a0 = 1 / (Ap + Ancs + sa ); CalcT a1 = 2 * (An + Apcs); SetStage( a1*a0, a2*a0, b0*a0, b1*a0, b2*a0 ); } template void BiquadLowShelf::Setup( CalcT normFreq, CalcT dB, CalcT shelfSlope ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); CalcT al = sn / 2 * ::sqrt( (A + 1/A) * (1/shelfSlope - 1) + 2 ); CalcT sa = 2 * ::sqrt( A ) * al; SetupCommon( cs, A, sa ); } // This could be optimized further template void BiquadLowShelf::SetupFast( CalcT normFreq, CalcT dB, CalcT shelfSlope ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); CalcT al = sn / 2 * fastsqrt1( (A + 1/A) * (1/shelfSlope - 1) + 2 ); CalcT sa = 2 * fastsqrt1( A ) * al; SetupCommon( cs, A, sa ); } //-------------------------------------------------------------------------- template class BiquadHighShelf : public Biquad { public: void Setup ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 ); void SetupFast ( CalcT normFreq, CalcT dB, CalcT shelfSlope=1.0 ); protected: void SetupCommon ( CalcT cs, CalcT A, CalcT sa ); }; //-------------------------------------------------------------------------- template void BiquadHighShelf::SetupCommon( CalcT cs, CalcT A, CalcT sa ) { CalcT An = A-1; CalcT Ap = A+1; CalcT Ancs = An*cs; CalcT Apcs = Ap*cs; CalcT b0 = A * (Ap + Ancs + sa ); CalcT b1 = -2 * A * (An + Apcs); CalcT b2 = A * (Ap + Ancs - sa ); CalcT a0 = (Ap - Ancs + sa ); CalcT a2 = Ancs + sa - Ap; CalcT a1 = -2 * (An - Apcs); SetStage( a1/a0, a2/a0, b0/a0, b1/a0, b2/a0 ); } template void BiquadHighShelf::Setup( CalcT normFreq, CalcT dB, CalcT shelfSlope ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT cs = cos(w0); CalcT sn = sin(w0); CalcT alph = sn / 2 * ::sqrt( (A + 1/A) * (1/shelfSlope - 1) + 2 ); CalcT sa = 2 * ::sqrt( A ) * alph; SetupCommon( cs, A, sa ); } template void BiquadHighShelf::SetupFast( CalcT normFreq, CalcT dB, CalcT shelfSlope ) { CalcT A = pow( 10, dB/40 ); CalcT w0 = 2 * kPi * normFreq; CalcT sn, cs; fastsincos( w0, &sn, &cs ); CalcT alph = sn / 2 * fastsqrt1( (A + 1/A) * (1/shelfSlope - 1) + 2 ); CalcT sa = 2 * fastsqrt1( A ) * alph; SetupCommon( cs, A, sa ); } //-------------------------------------------------------------------------- // // General N-Pole IIR Filter // //-------------------------------------------------------------------------- template class PoleFilter : public CascadeStages { public: PoleFilter(); virtual int CountPoles ( void )=0; virtual int CountZeroes ( void )=0; virtual Complex GetPole ( int i )=0; virtual Complex GetZero ( int i )=0; protected: virtual Complex GetSPole ( int i, CalcT wc )=0; protected: // Determines the method of obtaining // unity gain coefficients in the passband. enum Hint { // No normalizating hintNone, // Use Brent's method to find the maximum hintBrent, // Use the response at a given frequency hintPassband }; Complex BilinearTransform ( const Complex &c ); Complex BandStopTransform ( int i, const Complex &c ); Complex BandPassTransform ( int i, const Complex &c ); Complex GetBandStopPole ( int i ); Complex GetBandStopZero ( int i ); Complex GetBandPassPole ( int i ); Complex GetBandPassZero ( int i ); void Normalize ( void ); void Prepare ( void ); virtual void BrentHint ( CalcT *w0, CalcT *w1 ); virtual CalcT PassbandHint( void ); protected: Hint m_hint; int m_n; CalcT m_wc; CalcT m_wc2; }; //-------------------------------------------------------------------------- template inline PoleFilter::PoleFilter( void ) { m_hint=hintNone; } template inline Complex PoleFilter::BilinearTransform( const Complex &c ) { return (c+1.)/(-c+1.); } template inline Complex PoleFilter::BandStopTransform( int i, const Complex &c ) { CalcT a=cos((m_wc+m_wc2)*.5) / cos((m_wc-m_wc2)*.5); CalcT b=tan((m_wc-m_wc2)*.5); Complex c2(0); c2=addmul( c2, 4*(b*b+a*a-1), c ); c2+=8*(b*b-a*a+1); c2*=```