From f68d2e06bcfb81efda107d3b4c3aa7dbc2d73bc2 Mon Sep 17 00:00:00 2001 From: Robin Gareus Date: Thu, 6 Oct 2016 00:16:44 +0200 Subject: update qm-dsp library --- libs/qm-dsp/base/KaiserWindow.cpp | 65 ++++++++++++++++++++++ libs/qm-dsp/base/KaiserWindow.h | 105 ++++++++++++++++++++++++++++++++++++ libs/qm-dsp/base/Pitch.cpp | 2 +- libs/qm-dsp/base/Pitch.h | 4 ++ libs/qm-dsp/base/SincWindow.cpp | 45 ++++++++++++++++ libs/qm-dsp/base/SincWindow.h | 61 +++++++++++++++++++++ libs/qm-dsp/base/Window.h | 110 ++++++++++++++++++++++++-------------- 7 files changed, 351 insertions(+), 41 deletions(-) create mode 100644 libs/qm-dsp/base/KaiserWindow.cpp create mode 100644 libs/qm-dsp/base/KaiserWindow.h create mode 100644 libs/qm-dsp/base/SincWindow.cpp create mode 100644 libs/qm-dsp/base/SincWindow.h (limited to 'libs/qm-dsp/base') diff --git a/libs/qm-dsp/base/KaiserWindow.cpp b/libs/qm-dsp/base/KaiserWindow.cpp new file mode 100644 index 0000000000..4fe838e090 --- /dev/null +++ b/libs/qm-dsp/base/KaiserWindow.cpp @@ -0,0 +1,65 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP library + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#include "KaiserWindow.h" + +#include "maths/MathUtilities.h" + +KaiserWindow::Parameters +KaiserWindow::parametersForTransitionWidth(double attenuation, + double transition) +{ + Parameters p; + p.length = 1 + (attenuation > 21.0 ? + ceil((attenuation - 7.95) / (2.285 * transition)) : + ceil(5.79 / transition)); + p.beta = (attenuation > 50.0 ? + 0.1102 * (attenuation - 8.7) : + attenuation > 21.0 ? + 0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) : + 0); + return p; +} + +static double besselTerm(double x, int i) +{ + if (i == 0) { + return 1; + } else { + double f = MathUtilities::factorial(i); + return pow(x/2, i*2) / (f*f); + } +} + +static double bessel0(double x) +{ + double b = 0.0; + for (int i = 0; i < 20; ++i) { + b += besselTerm(x, i); + } + return b; +} + +void +KaiserWindow::init() +{ + double denominator = bessel0(m_beta); + bool even = (m_length % 2 == 0); + for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) { + double k = double(2*i) / double(m_length-1) - 1.0; + m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator); + } + for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) { + m_window.push_back(m_window[int(m_length/2) - i - 1]); + } +} diff --git a/libs/qm-dsp/base/KaiserWindow.h b/libs/qm-dsp/base/KaiserWindow.h new file mode 100644 index 0000000000..f16a4b6c16 --- /dev/null +++ b/libs/qm-dsp/base/KaiserWindow.h @@ -0,0 +1,105 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP library + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#ifndef KAISER_WINDOW_H +#define KAISER_WINDOW_H + +#include +#include + +/** + * Kaiser window: A windower whose bandwidth and sidelobe height + * (signal-noise ratio) can be specified. These parameters are traded + * off against the window length. + */ +class KaiserWindow +{ +public: + struct Parameters { + int length; + double beta; + }; + + /** + * Construct a Kaiser windower with the given length and beta + * parameter. + */ + KaiserWindow(Parameters p) : m_length(p.length), m_beta(p.beta) { init(); } + + /** + * Construct a Kaiser windower with the given attenuation in dB + * and transition width in samples. + */ + static KaiserWindow byTransitionWidth(double attenuation, + double transition) { + return KaiserWindow + (parametersForTransitionWidth(attenuation, transition)); + } + + /** + * Construct a Kaiser windower with the given attenuation in dB + * and transition bandwidth in Hz for the given samplerate. + */ + static KaiserWindow byBandwidth(double attenuation, + double bandwidth, + double samplerate) { + return KaiserWindow + (parametersForBandwidth(attenuation, bandwidth, samplerate)); + } + + /** + * Obtain the parameters necessary for a Kaiser window of the + * given attenuation in dB and transition width in samples. + */ + static Parameters parametersForTransitionWidth(double attenuation, + double transition); + + /** + * Obtain the parameters necessary for a Kaiser window of the + * given attenuation in dB and transition bandwidth in Hz for the + * given samplerate. + */ + static Parameters parametersForBandwidth(double attenuation, + double bandwidth, + double samplerate) { + return parametersForTransitionWidth + (attenuation, (bandwidth * 2 * M_PI) / samplerate); + } + + int getLength() const { + return m_length; + } + + const double *getWindow() const { + return m_window.data(); + } + + void cut(double *src) const { + cut(src, src); + } + + void cut(const double *src, double *dst) const { + for (int i = 0; i < m_length; ++i) { + dst[i] = src[i] * m_window[i]; + } + } + +private: + int m_length; + double m_beta; + std::vector m_window; + + void init(); +}; + +#endif diff --git a/libs/qm-dsp/base/Pitch.cpp b/libs/qm-dsp/base/Pitch.cpp index e86d6792b7..9ba8b4b3f7 100644 --- a/libs/qm-dsp/base/Pitch.cpp +++ b/libs/qm-dsp/base/Pitch.cpp @@ -39,7 +39,7 @@ Pitch::getPitchForFrequency(float frequency, midiPitch = midiPitch + 1; centsOffset = -(100.0 - centsOffset); } - + if (centsOffsetReturn) *centsOffsetReturn = centsOffset; return midiPitch; } diff --git a/libs/qm-dsp/base/Pitch.h b/libs/qm-dsp/base/Pitch.h index 7b1d6c96aa..5b55ecc4a3 100644 --- a/libs/qm-dsp/base/Pitch.h +++ b/libs/qm-dsp/base/Pitch.h @@ -15,6 +15,10 @@ #ifndef _PITCH_H_ #define _PITCH_H_ +/** + * Convert between musical pitch (i.e. MIDI pitch number) and + * fundamental frequency. + */ class Pitch { public: diff --git a/libs/qm-dsp/base/SincWindow.cpp b/libs/qm-dsp/base/SincWindow.cpp new file mode 100644 index 0000000000..5d31a631e3 --- /dev/null +++ b/libs/qm-dsp/base/SincWindow.cpp @@ -0,0 +1,45 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP library + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#include "SincWindow.h" + +#include + +void +SincWindow::init() +{ + if (m_length < 1) { + return; + } else if (m_length < 2) { + m_window.push_back(1); + return; + } else { + + int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2); + int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2); + double m = 2 * M_PI / m_p; + + for (int i = 0; i < n0; ++i) { + double x = ((m_length / 2) - i) * m; + m_window.push_back(sin(x) / x); + } + + m_window.push_back(1.0); + + for (int i = 1; i < n1; ++i) { + double x = i * m; + m_window.push_back(sin(x) / x); + } + } +} + diff --git a/libs/qm-dsp/base/SincWindow.h b/libs/qm-dsp/base/SincWindow.h new file mode 100644 index 0000000000..bb35d90c20 --- /dev/null +++ b/libs/qm-dsp/base/SincWindow.h @@ -0,0 +1,61 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP library + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#ifndef SINC_WINDOW_H +#define SINC_WINDOW_H + +#include + +/** + * A window containing values of the sinc function, i.e. sin(x)/x with + * sinc(0) == 1, with x == 0 at the centre. + */ +class SincWindow +{ +public: + /** + * Construct a windower of the given length, containing the values + * of sinc(x) with x=0 in the middle, i.e. at sample (length-1)/2 + * for odd or (length/2)+1 for even length, such that the distance + * from -pi to pi (the nearest zero crossings either side of the + * peak) is p samples. + */ + SincWindow(int length, double p) : m_length(length), m_p(p) { init(); } + + int getLength() const { + return m_length; + } + + const double *getWindow() const { + return m_window.data(); + } + + void cut(double *src) const { + cut(src, src); + } + + void cut(const double *src, double *dst) const { + for (int i = 0; i < m_length; ++i) { + dst[i] = src[i] * m_window[i]; + } + } + +private: + int m_length; + double m_p; + std::vector m_window; + + void init(); +}; + +#endif diff --git a/libs/qm-dsp/base/Window.h b/libs/qm-dsp/base/Window.h index 02ca426d66..93909115c0 100644 --- a/libs/qm-dsp/base/Window.h +++ b/libs/qm-dsp/base/Window.h @@ -18,6 +18,7 @@ #include #include #include +#include enum WindowType { RectangularWindow, @@ -25,18 +26,28 @@ enum WindowType { HammingWindow, HanningWindow, BlackmanWindow, - GaussianWindow, - ParzenWindow + BlackmanHarrisWindow, + + FirstWindow = RectangularWindow, + LastWindow = BlackmanHarrisWindow }; +/** + * Various shaped windows for sample frame conditioning, including + * cosine windows (Hann etc) and triangular and rectangular windows. + */ template class Window { public: /** - * Construct a windower of the given type. + * Construct a windower of the given type and size. + * + * Note that the cosine windows are periodic by design, rather + * than symmetrical. (A window of size N is equivalent to a + * symmetrical window of size N+1 with the final element missing.) */ - Window(WindowType type, size_t size) : m_type(type), m_size(size) { encache(); } + Window(WindowType type, int size) : m_type(type), m_size(size) { encache(); } Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); } Window &operator=(const Window &w) { if (&w == this) return *this; @@ -46,79 +57,98 @@ public: return *this; } virtual ~Window() { delete[] m_cache; } - + void cut(T *src) const { cut(src, src); } void cut(const T *src, T *dst) const { - for (size_t i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i]; + for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i]; } WindowType getType() const { return m_type; } - size_t getSize() const { return m_size; } + int getSize() const { return m_size; } + + std::vector getWindowData() const { + std::vector d; + for (int i = 0; i < m_size; ++i) { + d.push_back(m_cache[i]); + } + return d; + } protected: WindowType m_type; - size_t m_size; + int m_size; T *m_cache; - + void encache(); }; template void Window::encache() { - size_t n = m_size; + int n = m_size; T *mult = new T[n]; - size_t i; + int i; for (i = 0; i < n; ++i) mult[i] = 1.0; switch (m_type) { - + case RectangularWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * 0.5; + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * 0.5; } break; - + case BartlettWindow: - for (i = 0; i < n/2; ++i) { - mult[i] = mult[i] * (i / T(n/2)); - mult[i + n/2] = mult[i + n/2] * (1.0 - (i / T(n/2))); + if (n == 2) { + mult[0] = mult[1] = 0; // "matlab compatible" + } else if (n == 3) { + mult[0] = 0; + mult[1] = mult[2] = 2./3.; + } else if (n > 3) { + for (i = 0; i < n/2; ++i) { + mult[i] = mult[i] * (i / T(n/2)); + mult[i + n - n/2] = mult[i + n - n/2] * (1.0 - (i / T(n/2))); + } } break; - + case HammingWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n)); + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n)); + } } break; - + case HanningWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n)); + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n)); + } } break; - + case BlackmanWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n) - + 0.08 * cos(4 * M_PI * i / n)); + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n) + + 0.08 * cos(4 * M_PI * i / n)); + } } break; - - case GaussianWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * exp((-1.0 / (n*n)) * ((T(2*i) - n) * - (T(2*i) - n))); - } - break; - - case ParzenWindow: - for (i = 0; i < n; ++i) { - mult[i] = mult[i] * (1.0 - fabs((T(2*i) - n) / T(n + 1))); + + case BlackmanHarrisWindow: + if (n > 1) { + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.35875 + - 0.48829 * cos(2 * M_PI * i / n) + + 0.14128 * cos(4 * M_PI * i / n) + - 0.01168 * cos(6 * M_PI * i / n)); + } } break; } - + m_cache = mult; } -- cgit v1.2.3