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 FMD{1..5}, atan and atan2 approximations

Browse source
This commit is contained in:
Ahmet Inan 2019-02-06 14:29:07 +01:00
commit 99baec64a5
3 changed files with 217 additions and 0 deletions
Download patch file Download diff file Expand all files Collapse all files

8
README.md
View file

@ -44,6 +44,14 @@ The following [infinite impulse response](https://en.wikipedia.org/wiki/Infinite
[Numerically controlled oscillator](https://en.wikipedia.org/wiki/Numerically_controlled_oscillator) implemented using a [phasor](https://en.wikipedia.org/wiki/Phasor) and [complex multiplication](https://en.wikipedia.org/wiki/Complex_number#Multiplication) instead of a [lookup table](https://en.wikipedia.org/wiki/Lookup_table).
### [fmd.hh](fmd.hh)
[Frequency modulation](https://en.wikipedia.org/wiki/Frequency_modulation) [demodulation](https://en.wikipedia.org/wiki/Demodulation) with and without [atan2](https://en.wikipedia.org/wiki/Atan2).
### [atan2.hh](atan2.hh)
[atan](https://en.wikipedia.org/wiki/Inverse_trigonometric_functions) and [atan2](https://en.wikipedia.org/wiki/Atan2).
### [const.hh](const.hh)
Some constants we need

74
atan2.hh Normal file
View file

@ -0,0 +1,74 @@
/*
atan and atan2 approximations
Constants below lifted from the Cephes Mathematical Library:
https://www.netlib.org/cephes/cmath.tgz
Copyright 2019 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
#include "const.hh"
namespace DSP {
template <typename TYPE>
TYPE atan_kernel(TYPE x)
{
static const TYPE
P0 = -8.750608600031904122785E-1,
P1 = -1.615753718733365076637E1,
P2 = -7.500855792314704667340E1,
P3 = -1.228866684490136173410E2,
P4 = -6.485021904942025371773E1,
Q0 = 2.485846490142306297962E1,
Q1 = 1.650270098316988542046E2,
Q2 = 4.328810604912902668951E2,
Q3 = 4.853903996359136964868E2,
Q4 = 1.945506571482613964425E2;
TYPE xx = x * x;
return x + x * xx *
(xx*(xx*(xx*(xx*P0+P1)+P2)+P3)+P4) /
(xx*(xx*(xx*(xx*(xx+Q0)+Q1)+Q2)+Q3)+Q4);
}
template <typename TYPE>
TYPE atan_approx(TYPE x)
{
static const TYPE
MOREBITS = 6.123233995736765886130E-17;
if (x <= TYPE(0.66))
return atan_kernel(x);
if (x > TYPE(2.41421356237309504880))
return Const<TYPE>::HalfPi() + MOREBITS -
atan_kernel(TYPE(1) / x);
return Const<TYPE>::FourthPi() + TYPE(0.5) * MOREBITS +
atan_kernel((x - TYPE(1)) / (x + TYPE(1)));
}
template <typename TYPE>
TYPE atan(TYPE x)
{
return x >= TYPE(0) ? atan_approx(x) : -atan_approx(-x);
}
template <typename TYPE>
TYPE atan2(TYPE y, TYPE x)
{
if (y == TYPE(0))
return x >= TYPE(0) ?
copysign(TYPE(0), y) :
copysign(Const<TYPE>::Pi(), y);
if (x == TYPE(0))
return copysign(Const<TYPE>::HalfPi(), y);
TYPE phase = atan(y / x);
if (x < TYPE(0))
return phase > TYPE(0) ?
phase - Const<TYPE>::Pi() :
phase + Const<TYPE>::Pi();
return phase;
}
}

135
fmd.hh Normal file
View file

@ -0,0 +1,135 @@
/*
FM Demodulation
Copyright 2019 Ahmet Inan <inan@aicodix.de>
*/
#pragma once
namespace DSP {
template <typename TYPE>
class FMD1
{
typedef TYPE complex_type;
typedef typename complex_type::value_type value_type;
value_type prev;
value_type scale;
public:
constexpr FMD1() : prev(0), scale(0)
{
}
void bandwidth(value_type bw)
{
scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
}
value_type operator()(complex_type input)
{
value_type phase = arg(input);
value_type delta = phase - prev;
prev = phase;
return scale * delta;
}
};
template <typename TYPE>
class FMD2
{
typedef TYPE complex_type;
typedef typename complex_type::value_type value_type;
complex_type prev;
value_type scale;
public:
constexpr FMD2() : prev(1, 0), scale(0)
{
}
void bandwidth(value_type bw)
{
scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
}
value_type operator()(complex_type input)
{
value_type phase = arg(input / prev);
prev = input;
return scale * phase;
}
};
template <typename TYPE>
class FMD3
{
typedef TYPE complex_type;
typedef typename complex_type::value_type value_type;
complex_type prev;
value_type scale;
public:
constexpr FMD3() : prev(1, 0), scale(0)
{
}
void bandwidth(value_type bw)
{
scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
}
value_type operator()(complex_type input)
{
value_type phase = arg(input * conj(prev));
prev = input;
return scale * phase;
}
};
template <typename TYPE>
class FMD4
{
typedef TYPE complex_type;
typedef typename complex_type::value_type value_type;
complex_type prev;
value_type scale;
public:
constexpr FMD4() : prev(1, 0), scale(0)
{
}
void bandwidth(value_type bw)
{
scale = value_type(1) / (bw * DSP::Const<value_type>::Pi());
}
value_type operator()(complex_type input)
{
value_type phase =
prev.real() * input.imag() -
prev.imag() * input.real();
phase /= norm(input);
prev = input;
return scale * phase;
}
};
template <typename TYPE>
class FMD5
{
typedef TYPE complex_type;
typedef typename complex_type::value_type value_type;
complex_type prev1, prev2;
value_type scale;
public:
constexpr FMD5() : prev1(1, 0), prev2(1, 0), scale(0)
{
}
void bandwidth(value_type bw)
{
scale = value_type(1) / (bw * DSP::Const<value_type>::TwoPi());
}
value_type operator()(complex_type input)
{
complex_type diff = input - prev2;
value_type phase =
prev1.real() * diff.imag() -
prev1.imag() * diff.real();
phase /= norm(prev1);
prev2 = prev1; prev1 = input;
return scale * phase;
}
};
}