diff --git a/README.md b/README.md index 7847b75..24e2fc5 100644 --- a/README.md +++ b/README.md @@ -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) + diff --git a/tests/unit_circle_test.cc b/tests/unit_circle_test.cc new file mode 100644 index 0000000..fc2d9a1 --- /dev/null +++ b/tests/unit_circle_test.cc @@ -0,0 +1,35 @@ +/* +Test for the UnitCircle functions + +Copyright 2019 Ahmet Inan +*/ + +#include +#include +#include +#include +#include "unit_circle.hh" +#include "const.hh" + +template +void test(int MAX) +{ + for (int N = 1; N <= MAX; ++N) { + for (int n = 0; n < N; ++n) { + PRECISE x = DSP::Const::TwoPi() * PRECISE(n) / PRECISE(N); + PRECISE cos_err = std::abs(DSP::UnitCircle::cos(n, N) - std::cos(x)); + assert(cos_err < std::numeric_limits::epsilon()); + PRECISE sin_err = std::abs(DSP::UnitCircle::sin(n, N) - std::sin(x)); + assert(sin_err < std::numeric_limits::epsilon()); + } + } +} + +int main() +{ + test(5000); + test(5000); + std::cerr << "UnitCircle functions test passed!" << std::endl; + return 0; +} + diff --git a/unit_circle.hh b/unit_circle.hh new file mode 100644 index 0000000..0df6b85 --- /dev/null +++ b/unit_circle.hh @@ -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 +*/ + +#pragma once + +#include "const.hh" + +namespace DSP { + +template +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::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)); + } +}; + +} +