added UnitCircle with sine and cosine functions

This commit is contained in:
Ahmet Inan 2019-01-18 16:56:58 +01:00
commit 15843df336
3 changed files with 118 additions and 0 deletions

View file

@ -72,3 +72,9 @@ Some everyday helpers:
When working with [Analog-to-digital](https://en.wikipedia.org/wiki/Analog-to-digital_converter) and [Digital-to-analog](https://en.wikipedia.org/wiki/Digital-to-analog_converter) converters, we often face the ugly truth, that we can't always have a precise [Sampling](https://en.wikipedia.org/wiki/Sampling_(signal_processing)) rate.
But if we can estimate the Sampling frequency offset, we can correct it by [Resampling](https://en.wikipedia.org/wiki/Sample-rate_conversion) the sampled data.
### [unit_circle.hh](unit_circle.hh)
Sometimes we only need [trigonometric functions](https://en.wikipedia.org/wiki/Trigonometric_functions) that stay on the [unit circle](https://en.wikipedia.org/wiki/Unit_circle):
* [sine function](https://en.wikipedia.org/wiki/Sine)
* [cosine function](https://en.wikipedia.org/wiki/Trigonometric_functions#cosine)

35
tests/unit_circle_test.cc Normal file
View file

@ -0,0 +1,35 @@
/*
Test for the UnitCircle functions
Copyright 2019 Ahmet Inan <inan@aicodix.de>
*/
#include <cassert>
#include <iostream>
#include <limits>
#include <cmath>
#include "unit_circle.hh"
#include "const.hh"
template <typename TYPE, typename PRECISE>
void test(int MAX)
{
for (int N = 1; N <= MAX; ++N) {
for (int n = 0; n < N; ++n) {
PRECISE x = DSP::Const<PRECISE>::TwoPi() * PRECISE(n) / PRECISE(N);
PRECISE cos_err = std::abs(DSP::UnitCircle<TYPE>::cos(n, N) - std::cos(x));
assert(cos_err < std::numeric_limits<TYPE>::epsilon());
PRECISE sin_err = std::abs(DSP::UnitCircle<TYPE>::sin(n, N) - std::sin(x));
assert(sin_err < std::numeric_limits<TYPE>::epsilon());
}
}
}
int main()
{
test<float, double>(5000);
test<double, long double>(5000);
std::cerr << "UnitCircle functions test passed!" << std::endl;
return 0;
}

77
unit_circle.hh Normal file
View file

@ -0,0 +1,77 @@
/*
Trigonometric functions on the unit circle
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>
class UnitCircle
{
static constexpr TYPE
a1 = 1.0,
a2 = -0.5,
a3 = -1.66666666666666307295E-1,
a4 = 4.16666666666665929218E-2,
a5 = 8.33333333332211858878E-3,
a6 = -1.38888888888730564116E-3,
a7 = -1.98412698295895385996E-4,
a8 = 2.48015872888517045348E-5,
a9 = 2.75573136213857245213E-6,
a10 = -2.75573141792967388112E-7,
a11 = -2.50507477628578072866E-8,
a12 = 2.08757008419747316778E-9,
a13 = 1.58962301576546568060E-10,
a14 = -1.13585365213876817300E-11;
static constexpr TYPE cosine_kernel(TYPE xx)
{
return a1 + xx * (xx * (xx * (xx * (xx * (xx * (xx * a14 + a12) + a10) + a8) + a6) + a4) + a2);
}
static constexpr TYPE sine_kernel(TYPE x, TYPE xx)
{
return x * (xx * (xx * (xx * (xx * (xx * (xx * a13 + a11) + a9) + a7) + a5) + a3) + a1);
}
static constexpr TYPE cosine_approx(TYPE x)
{
return cosine_kernel(x * x);
}
static constexpr TYPE sine_approx(TYPE x)
{
return sine_kernel(x, x * x);
}
public:
static constexpr TYPE rad(int n, int N)
{
return Const<TYPE>::TwoPi() * TYPE(n) / TYPE(N);
}
static constexpr TYPE cos(int n, int N)
{
return
8*n < 1*N ? cosine_approx(rad(n, N)) :
8*n < 3*N ? -sine_approx(rad(4*n-N, 4*N)) :
8*n < 5*N ? -cosine_approx(rad(2*n-N, 2*N)) :
8*n < 7*N ? sine_approx(rad(4*n-3*N, 4*N)) :
cosine_approx(rad(n-N, N));
}
static constexpr TYPE sin(int n, int N)
{
return
8*n < 1*N ? sine_approx(rad(n, N)) :
8*n < 3*N ? cosine_approx(rad(4*n-N, 4*N)) :
8*n < 5*N ? -sine_approx(rad(2*n-N, 2*N)) :
8*n < 7*N ? -cosine_approx(rad(4*n-3*N, 4*N)) :
sine_approx(rad(n-N, N));
}
};
}