Mirror/aicodix___dsp
1
0
Fork 0
mirror of https://github.com/aicodix/dsp.git synced 2026-04-27 14:30:36 +00:00
Code Issues Projects Releases Packages Wiki Activity

added EMA and Biquad

Browse source
This commit is contained in:
Ahmet Inan 2019-02-04 10:46:53 +01:00
commit 639204fed8
3 changed files with 171 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

12
README.md
View file

@ -28,6 +28,18 @@ Implemented are the following [finite impulse response](https://en.wikipedia.org
* [high-pass filter](https://en.wikipedia.org/wiki/High-pass_filter)
* [band-pass filter](https://en.wikipedia.org/wiki/Band-pass_filter)
### [ema.hh](ema.hh)
The [exponential moving average](https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average) is an [infinite impulse response](https://en.wikipedia.org/wiki/Infinite_impulse_response) [low-pass filter](https://en.wikipedia.org/wiki/Low-pass_filter).
There is also support for cascading, to improve [roll-off](https://en.wikipedia.org/wiki/Roll-off) while a correction factor helps to keep the same [cutoff frequency](https://en.wikipedia.org/wiki/Cutoff_frequency).
### [biquad.hh](biquad.hh)
The following [infinite impulse response](https://en.wikipedia.org/wiki/Infinite_impulse_response) [digital biquad filter](https://en.wikipedia.org/wiki/Digital_biquad_filter) implementations are available:
* [second-order Butterworth low pass filter](https://en.wikipedia.org/wiki/Butterworth_filter)
* [2n-order Butterworth cascade of second-order low pass filters](https://en.wikipedia.org/wiki/Butterworth_filter)
### [const.hh](const.hh)
Some constants we need

71
biquad.hh Normal file
View file

@ -0,0 +1,71 @@
/*
Digital biquad filter
Copyright 2019 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include "const.hh"
#include "unit_circle.hh"
namespace DSP {
template <typename TYPE, typename VALUE>
class Biquad
{
TYPE x1, x2, y1, y2;
VALUE b0a0, b1a0, b2a0, a1a0, a2a0;
public:
constexpr Biquad() : x1(0), x2(0), y1(0), y2(0),
b0a0(1), b1a0(0), b2a0(0), a1a0(0), a2a0(0)
{
}
void lowpass(int n, int N, VALUE Q = Const<VALUE>::InvSqrtTwo())
{
VALUE alpha = UnitCircle<VALUE>::sin(n, N) / (VALUE(2) * Q),
cn = UnitCircle<VALUE>::cos(n, N),
b0 = (VALUE(1) - cn) / VALUE(2),
b1 = VALUE(1) - cn,
b2 = (VALUE(1) - cn) / VALUE(2),
a0 = VALUE(1) + alpha,
a1 = -VALUE(2) * cn,
a2 = VALUE(1) - alpha;
b0a0 = b0 / a0;
b1a0 = b1 / a0;
b2a0 = b2 / a0;
a1a0 = a1 / a0;
a2a0 = a2 / a0;
}
TYPE operator()(TYPE x0)
{
TYPE y0 = b0a0*x0 + b1a0*x1 + b2a0*x2
- a1a0*y1 - a2a0*y2;
x2 = x1; x1 = x0;
y2 = y1; y1 = y0;
return y0;
}
};
template <typename TYPE, typename VALUE, int NUM = 1>
class BiquadCascade
{
static const int ORDER = 2 * NUM;
Biquad<TYPE, VALUE> cascade[NUM];
public:
void lowpass(int n, int N)
{
for (int i = 0; i < NUM; ++i)
cascade[i].lowpass(n, N, VALUE(1) / (VALUE(2) *
UnitCircle<VALUE>::cos(2*i+1, 4*ORDER)));
}
TYPE operator()(TYPE input)
{
for (int i = 0; i < NUM; ++i)
input = cascade[i](input);
return input;
}
};
}

88
ema.hh Normal file
View file

@ -0,0 +1,88 @@
/*
Exponential moving average
Copyright 2019 Ahmet Inan <inan@aicodix.de>
*/
#include "utils.hh"
#include "unit_circle.hh"
#pragma once
namespace DSP {
template <typename TYPE, typename VALUE>
class EMA
{
TYPE prev;
VALUE alpha;
static VALUE cutoff_alpha(int n, int N)
{
VALUE x = UnitCircle<VALUE>::cos(n, N);
return x-VALUE(1)+sqrt(x*(x-VALUE(4))+VALUE(3));
}
public:
EMA() : prev(0), alpha(1)
{
}
void cutoff(int n, int N)
{
alpha = cutoff_alpha(n, N);
}
TYPE operator()(TYPE input)
{
return prev = lerp(alpha, prev, input);
}
};
template <typename TYPE, typename VALUE, int ORDER = 1>
class EMACascade
{
TYPE prev[ORDER];
VALUE alpha;
static constexpr VALUE cascade_factor(int n)
{
// sqrt(e(1/n*l(2))-1)
return
n == 2 ? .6435942529055826247354434374182098089242027424440076511561520093520L :
n == 3 ? .5098245285339585980866134620415975171004662971729213347866667412748L :
n == 4 ? .4349794420460822959023617403105260887838801269492020945122292703957L :
n == 5 ? .3856142567346739149847213965799930737314794827929705449459805476567L :
n == 6 ? .3499457790992384314823879950215014587750945565870355960474164484675L :
n == 7 ? .3226290651410878970708971817580301448745259798477682170908870012270L :
n == 8 ? .3008450309798346315492310951214492771360484456002599192820777432131L :
n == 9 ? .2829482972069387977799123584440613415063420296903386876512703264517L :
n == 10 ? .2679056970956257290860485974661891235831324413057223464282157983007L :
n == 11 ? .2550315459702243937103098780342370981449477556708628034249027758407L :
n == 12 ? .2438505574307659124524083486415529630644763241448445300317390218311L :
n == 13 ? .2340215299532217691188557628487431016851279721568792383213019527102L :
n == 14 ? .2252923404228812762935182605055552181633447787610798690912140601327L :
n == 15 ? .2174721196397982365782144895603276387766971846333067133969338153904L :
n == 16 ? .2104133608576552353576644943161359426249259240013118976757629219108L :
1;
}
static VALUE cutoff_alpha(int n, int N)
{
VALUE x = UnitCircle<VALUE>::cos(n, N);
return x-VALUE(1)+sqrt(x*(x-VALUE(4))+VALUE(3));
}
public:
EMACascade() : alpha(1)
{
for (int i = 0; i < ORDER; ++i)
prev[i] = TYPE(0);
}
void cutoff(int n, int N)
{
alpha = cutoff_alpha(nearbyint(n / cascade_factor(ORDER)), N);
}
TYPE operator()(TYPE input)
{
for (int i = 0; i < ORDER; ++i)
prev[i] = input = lerp(alpha, prev[i], input);
return input;
}
};
}