diff --git a/window.hh b/window.hh index ef65291..006c57d 100644 --- a/window.hh +++ b/window.hh @@ -12,101 +12,95 @@ Copyright 2018 Ahmet Inan namespace DSP { +template +struct Window +{ + virtual TYPE operator () (int, int) = 0; +}; + template -class Rect +class Taps { TYPE w[TAPS]; public: - Rect() + Taps(Window *func) { for (int n = 0; n < TAPS; ++n) - w[n] = TYPE(1); + w[n] = (*func)(n, TAPS); } inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } inline operator const TYPE * () const { return w; } }; -template -class Hann +template +struct Rect : public Window { - TYPE w[TAPS]; -public: - Hann() - { - for (int n = 0; n < TAPS; ++n) - w[n] = TYPE(0.5) * (TYPE(1) - std::cos(Const::TwoPi() * TYPE(n) / TYPE(TAPS - 1))); - } - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } + TYPE operator () (int n, int N) { return n >= 0 && n < N ? 1 : 0; } }; -template -class Hamming +template +struct Hann : public Window { - TYPE w[TAPS]; -public: - Hamming() + TYPE operator () (int n, int N) { - for (int n = 0; n < TAPS; ++n) - w[n] = TYPE(0.54) - TYPE(0.46) * std::cos(Const::TwoPi() * TYPE(n) / TYPE(TAPS - 1)); + return TYPE(0.5) * (TYPE(1) - std::cos(Const::TwoPi() * TYPE(n) / TYPE(N - 1))); } - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } }; -template -class Lanczos +template +struct Hamming : public Window { - TYPE w[TAPS]; - TYPE sinc(TYPE x) + TYPE operator () (int n, int N) + { + return TYPE(0.54) - TYPE(0.46) * std::cos(Const::TwoPi() * TYPE(n) / TYPE(N - 1)); + } +}; + +template +class Lanczos : public Window +{ + static TYPE sinc(TYPE x) { return TYPE(0) == x ? TYPE(1) : std::sin(Const::Pi() * x) / (Const::Pi() * x); } public: - Lanczos() + TYPE operator () (int n, int N) { - for (int n = 0; n < TAPS; ++n) - w[n] = sinc(TYPE(2 * n) / TYPE(TAPS - 1) - TYPE(1)); + return sinc(TYPE(2 * n) / TYPE(N - 1) - TYPE(1)); } - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } }; -template -class Blackman +template +class Blackman : public Window { - TYPE w[TAPS]; + TYPE a0, a1, a2; public: - Blackman(TYPE a0, TYPE a1, TYPE a2) - { - for (int n = 0; n < TAPS; ++n) - w[n] = a0 - a1 * std::cos(Const::TwoPi() * TYPE(n) / TYPE(TAPS - 1)) + a2 * std::cos(Const::FourPi() * TYPE(n) / TYPE(TAPS - 1)); - } + Blackman(TYPE a0, TYPE a1, TYPE a2) : a0(a0), a1(a1), a2(a2) {} Blackman(TYPE a) : Blackman((TYPE(1) - a) / TYPE(2), TYPE(0.5), a / TYPE(2)) {} // "exact Blackman" Blackman() : Blackman(TYPE(7938) / TYPE(18608), TYPE(9240) / TYPE(18608), TYPE(1430) / TYPE(18608)) {} - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } -}; - -template -class Gauss -{ - TYPE w[TAPS]; -public: - Gauss(TYPE o) + TYPE operator () (int n, int N) { - for (int n = 0; n < TAPS; ++n) - w[n] = std::exp(- TYPE(0.5) * std::pow((TYPE(n) - TYPE(TAPS - 1) / TYPE(2)) / (o * TYPE(TAPS - 1) / TYPE(2)), TYPE(2))); + return a0 - a1 * std::cos(Const::TwoPi() * TYPE(n) / TYPE(N - 1)) + a2 * std::cos(Const::FourPi() * TYPE(n) / TYPE(N - 1)); } - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } }; -template -class Kaiser +template +class Gauss : public Window { - TYPE w[TAPS]; + TYPE o; +public: + Gauss(TYPE o) : o(o) {} + TYPE operator () (int n, int N) + { + return std::exp(- TYPE(0.5) * std::pow((TYPE(n) - TYPE(N - 1) / TYPE(2)) / (o * TYPE(N - 1) / TYPE(2)), TYPE(2))); + } +}; + +template +class Kaiser : public Window +{ + TYPE a; /* i0() implements the zero-th order modified Bessel function of the first kind: https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions:_I%CE%B1,_K%CE%B1 @@ -115,7 +109,7 @@ class Kaiser We obviously can't use the factorial here, so let's get rid of it: $= 1 + \left(\frac{x}{2 \cdot 1}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\right)^2 + \left(\frac{x}{2 \cdot 1}\cdot \frac{x}{2 \cdot 2}\cdot \frac{x}{2 \cdot 3}\right)^2 + .. = 1 + \sum_{m=1}^\infty \left(\prod_{n=1}^m \frac{x}{2n}\right)^2$ */ - TYPE i0(TYPE x) + static TYPE i0(TYPE x) { Kahan sum(1.0); TYPE val = 1.0; @@ -130,13 +124,11 @@ class Kaiser return sum(); } public: - Kaiser(TYPE a) + Kaiser(TYPE a) : a(a) {} + TYPE operator () (int n, int N) { - for (int n = 0; n < TAPS; ++n) - w[n] = i0(Const::Pi() * a * std::sqrt(TYPE(1) - std::pow(TYPE(2 * n) / TYPE(TAPS - 1) - TYPE(1), TYPE(2)))) / i0(Const::Pi() * a); + return i0(Const::Pi() * a * std::sqrt(TYPE(1) - std::pow(TYPE(2 * n) / TYPE(N - 1) - TYPE(1), TYPE(2)))) / i0(Const::Pi() * a); } - inline TYPE operator () (int n) { return n >= 0 && n < TAPS ? w[n] : 0; } - inline operator const TYPE * () const { return w; } }; }