-// ==================================================================================
-// Copyright (c) 2016 HiFi-LoFi
-//
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is furnished
-// to do so, subject to the following conditions:
-//
-// The above copyright notice and this permission notice shall be included in
-// all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
-// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
-// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
-// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
-// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-// ==================================================================================
-
-#include "AudioFFT.h"
-
-#include <cassert>
-#include <cmath>
-#include <cstring>
-
-
-#if defined(AUDIOFFT_APPLE_ACCELERATE)
- #define AUDIOFFT_APPLE_ACCELERATE_USED
- #include <Accelerate/Accelerate.h>
- #include <vector>
-#elif defined (AUDIOFFT_FFTW3)
- #define AUDIOFFT_FFTW3_USED
- #include <fftw3.h>
-#else
- #if !defined(AUDIOFFT_OOURA)
- #define AUDIOFFT_OOURA
- #endif
- #define AUDIOFFT_OOURA_USED
- #include <vector>
-#endif
-
-
-namespace audiofft
-{
-
- namespace details
- {
-
- static bool IsPowerOf2(size_t val)
- {
- return (val == 1 || (val & (val-1)) == 0);
- }
-
-
- template<typename TypeDest, typename TypeSrc>
- void ConvertBuffer(TypeDest* dest, const TypeSrc* src, size_t len)
- {
- for (size_t i=0; i<len; ++i)
- {
- dest[i] = static_cast<TypeDest>(src[i]);
- }
- }
-
-
- template<typename TypeDest, typename TypeSrc, typename TypeFactor>
- void ScaleBuffer(TypeDest* dest, const TypeSrc* src, const TypeFactor factor, size_t len)
- {
- for (size_t i=0; i<len; ++i)
- {
- dest[i] = static_cast<TypeDest>(static_cast<TypeFactor>(src[i]) * factor);
- }
- }
-
-
- // ================================================================
-
-
-#ifdef AUDIOFFT_OOURA_USED
-
- /**
- * @internal
- * @class OouraFFT
- * @brief FFT implementation based on the great radix-4 routines by Takuya Ooura
- */
- class OouraFFT : public AudioFFTImpl
- {
- public:
- OouraFFT() :
- AudioFFTImpl(),
- _size(0),
- _ip(),
- _w(),
- _buffer()
- {
- }
-
- virtual void init(size_t size) override
- {
- if (_size != size)
- {
- _ip.resize(2 + static_cast<int>(std::sqrt(static_cast<double>(size))));
- _w.resize(size / 2);
- _buffer.resize(size);
- _size = size;
-
- const int size4 = static_cast<int>(_size) / 4;
- makewt(size4, _ip.data(), _w.data());
- makect(size4, _ip.data(), _w.data() + size4);
- }
- }
-
- virtual void fft(const float* data, float* re, float* im) override
- {
- // Convert into the format as required by the Ooura FFT
- ConvertBuffer(&_buffer[0], data, _size);
-
- rdft(static_cast<int>(_size), +1, _buffer.data(), _ip.data(), _w.data());
-
- // Convert back to split-complex
- {
- double* b = &_buffer[0];
- double* bEnd = b + _size;
- float *r = re;
- float *i = im;
- while (b != bEnd)
- {
- *(r++) = static_cast<float>(*(b++));
- *(i++) = static_cast<float>(-(*(b++)));
- }
- }
- const size_t size2 = _size / 2;
- re[size2] = -im[0];
- im[0] = 0.0;
- im[size2] = 0.0;
- }
-
- virtual void ifft(float* data, const float* re, const float* im) override
- {
- // Convert into the format as required by the Ooura FFT
- {
- double* b = &_buffer[0];
- double* bEnd = b + _size;
- const float *r = re;
- const float *i = im;
- while (b != bEnd)
- {
- *(b++) = static_cast<double>(*(r++));
- *(b++) = -static_cast<double>(*(i++));
- }
- _buffer[1] = re[_size / 2];
- }
-
- rdft(static_cast<int>(_size), -1, _buffer.data(), _ip.data(), _w.data());
-
- // Convert back to split-complex
- ScaleBuffer(data, &_buffer[0], 2.0 / static_cast<double>(_size), _size);
- }
-
- private:
- size_t _size;
- std::vector<int> _ip;
- std::vector<double> _w;
- std::vector<double> _buffer;
-
- void rdft(int n, int isgn, double *a, int *ip, double *w)
- {
- int nw = ip[0];
- int nc = ip[1];
-
- if (isgn >= 0)
- {
- if (n > 4)
- {
- bitrv2(n, ip + 2, a);
- cftfsub(n, a, w);
- rftfsub(n, a, nc, w + nw);
- }
- else if (n == 4)
- {
- cftfsub(n, a, w);
- }
- double xi = a[0] - a[1];
- a[0] += a[1];
- a[1] = xi;
- }
- else
- {
- a[1] = 0.5 * (a[0] - a[1]);
- a[0] -= a[1];
- if (n > 4)
- {
- rftbsub(n, a, nc, w + nw);
- bitrv2(n, ip + 2, a);
- cftbsub(n, a, w);
- }
- else if (n == 4)
- {
- cftfsub(n, a, w);
- }
- }
- }
-
-
- /* -------- initializing routines -------- */
-
- void makewt(int nw, int *ip, double *w)
- {
- int j, nwh;
- double delta, x, y;
-
- ip[0] = nw;
- ip[1] = 1;
- if (nw > 2) {
- nwh = nw >> 1;
- delta = atan(1.0) / nwh;
- w[0] = 1;
- w[1] = 0;
- w[nwh] = cos(delta * nwh);
- w[nwh + 1] = w[nwh];
- if (nwh > 2) {
- for (j = 2; j < nwh; j += 2) {
- x = cos(delta * j);
- y = sin(delta * j);
- w[j] = x;
- w[j + 1] = y;
- w[nw - j] = y;
- w[nw - j + 1] = x;
- }
- bitrv2(nw, ip + 2, w);
- }
- }
- }
-
-
- void makect(int nc, int *ip, double *c)
- {
- int j, nch;
- double delta;
-
- ip[1] = nc;
- if (nc > 1) {
- nch = nc >> 1;
- delta = atan(1.0) / nch;
- c[0] = cos(delta * nch);
- c[nch] = 0.5 * c[0];
- for (j = 1; j < nch; j++) {
- c[j] = 0.5 * cos(delta * j);
- c[nc - j] = 0.5 * sin(delta * j);
- }
- }
- }
-
-
- /* -------- child routines -------- */
-
-
- void bitrv2(int n, int *ip, double *a)
- {
- int j, j1, k, k1, l, m, m2;
- double xr, xi, yr, yi;
-
- ip[0] = 0;
- l = n;
- m = 1;
- while ((m << 3) < l) {
- l >>= 1;
- for (j = 0; j < m; j++) {
- ip[m + j] = ip[j] + l;
- }
- m <<= 1;
- }
- m2 = 2 * m;
- if ((m << 3) == l) {
- for (k = 0; k < m; k++) {
- for (j = 0; j < k; j++) {
- j1 = 2 * j + ip[k];
- k1 = 2 * k + ip[j];
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- j1 += m2;
- k1 += 2 * m2;
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- j1 += m2;
- k1 -= m2;
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- j1 += m2;
- k1 += 2 * m2;
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- }
- j1 = 2 * k + m2 + ip[k];
- k1 = j1 + m2;
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- }
- } else {
- for (k = 1; k < m; k++) {
- for (j = 0; j < k; j++) {
- j1 = 2 * j + ip[k];
- k1 = 2 * k + ip[j];
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- j1 += m2;
- k1 += m2;
- xr = a[j1];
- xi = a[j1 + 1];
- yr = a[k1];
- yi = a[k1 + 1];
- a[j1] = yr;
- a[j1 + 1] = yi;
- a[k1] = xr;
- a[k1 + 1] = xi;
- }
- }
- }
- }
-
-
- void cftfsub(int n, double *a, double *w)
- {
- int j, j1, j2, j3, l;
- double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
- l = 2;
- if (n > 8) {
- cft1st(n, a, w);
- l = 8;
- while ((l << 2) < n) {
- cftmdl(n, l, a, w);
- l <<= 2;
- }
- }
- if ((l << 2) == n) {
- for (j = 0; j < l; j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = a[j + 1] + a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = a[j + 1] - a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- a[j2] = x0r - x2r;
- a[j2 + 1] = x0i - x2i;
- a[j1] = x1r - x3i;
- a[j1 + 1] = x1i + x3r;
- a[j3] = x1r + x3i;
- a[j3 + 1] = x1i - x3r;
- }
- } else {
- for (j = 0; j < l; j += 2) {
- j1 = j + l;
- x0r = a[j] - a[j1];
- x0i = a[j + 1] - a[j1 + 1];
- a[j] += a[j1];
- a[j + 1] += a[j1 + 1];
- a[j1] = x0r;
- a[j1 + 1] = x0i;
- }
- }
- }
-
-
- void cftbsub(int n, double *a, double *w)
- {
- int j, j1, j2, j3, l;
- double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
- l = 2;
- if (n > 8) {
- cft1st(n, a, w);
- l = 8;
- while ((l << 2) < n) {
- cftmdl(n, l, a, w);
- l <<= 2;
- }
- }
- if ((l << 2) == n) {
- for (j = 0; j < l; j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = -a[j + 1] - a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = -a[j + 1] + a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i - x2i;
- a[j2] = x0r - x2r;
- a[j2 + 1] = x0i + x2i;
- a[j1] = x1r - x3i;
- a[j1 + 1] = x1i - x3r;
- a[j3] = x1r + x3i;
- a[j3 + 1] = x1i + x3r;
- }
- } else {
- for (j = 0; j < l; j += 2) {
- j1 = j + l;
- x0r = a[j] - a[j1];
- x0i = -a[j + 1] + a[j1 + 1];
- a[j] += a[j1];
- a[j + 1] = -a[j + 1] - a[j1 + 1];
- a[j1] = x0r;
- a[j1 + 1] = x0i;
- }
- }
- }
-
-
- void cft1st(int n, double *a, double *w)
- {
- int j, k1, k2;
- double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
- double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
- x0r = a[0] + a[2];
- x0i = a[1] + a[3];
- x1r = a[0] - a[2];
- x1i = a[1] - a[3];
- x2r = a[4] + a[6];
- x2i = a[5] + a[7];
- x3r = a[4] - a[6];
- x3i = a[5] - a[7];
- a[0] = x0r + x2r;
- a[1] = x0i + x2i;
- a[4] = x0r - x2r;
- a[5] = x0i - x2i;
- a[2] = x1r - x3i;
- a[3] = x1i + x3r;
- a[6] = x1r + x3i;
- a[7] = x1i - x3r;
- wk1r = w[2];
- x0r = a[8] + a[10];
- x0i = a[9] + a[11];
- x1r = a[8] - a[10];
- x1i = a[9] - a[11];
- x2r = a[12] + a[14];
- x2i = a[13] + a[15];
- x3r = a[12] - a[14];
- x3i = a[13] - a[15];
- a[8] = x0r + x2r;
- a[9] = x0i + x2i;
- a[12] = x2i - x0i;
- a[13] = x0r - x2r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[10] = wk1r * (x0r - x0i);
- a[11] = wk1r * (x0r + x0i);
- x0r = x3i + x1r;
- x0i = x3r - x1i;
- a[14] = wk1r * (x0i - x0r);
- a[15] = wk1r * (x0i + x0r);
- k1 = 0;
- for (j = 16; j < n; j += 16) {
- k1 += 2;
- k2 = 2 * k1;
- wk2r = w[k1];
- wk2i = w[k1 + 1];
- wk1r = w[k2];
- wk1i = w[k2 + 1];
- wk3r = wk1r - 2 * wk2i * wk1i;
- wk3i = 2 * wk2i * wk1r - wk1i;
- x0r = a[j] + a[j + 2];
- x0i = a[j + 1] + a[j + 3];
- x1r = a[j] - a[j + 2];
- x1i = a[j + 1] - a[j + 3];
- x2r = a[j + 4] + a[j + 6];
- x2i = a[j + 5] + a[j + 7];
- x3r = a[j + 4] - a[j + 6];
- x3i = a[j + 5] - a[j + 7];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- x0r -= x2r;
- x0i -= x2i;
- a[j + 4] = wk2r * x0r - wk2i * x0i;
- a[j + 5] = wk2r * x0i + wk2i * x0r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[j + 2] = wk1r * x0r - wk1i * x0i;
- a[j + 3] = wk1r * x0i + wk1i * x0r;
- x0r = x1r + x3i;
- x0i = x1i - x3r;
- a[j + 6] = wk3r * x0r - wk3i * x0i;
- a[j + 7] = wk3r * x0i + wk3i * x0r;
- wk1r = w[k2 + 2];
- wk1i = w[k2 + 3];
- wk3r = wk1r - 2 * wk2r * wk1i;
- wk3i = 2 * wk2r * wk1r - wk1i;
- x0r = a[j + 8] + a[j + 10];
- x0i = a[j + 9] + a[j + 11];
- x1r = a[j + 8] - a[j + 10];
- x1i = a[j + 9] - a[j + 11];
- x2r = a[j + 12] + a[j + 14];
- x2i = a[j + 13] + a[j + 15];
- x3r = a[j + 12] - a[j + 14];
- x3i = a[j + 13] - a[j + 15];
- a[j + 8] = x0r + x2r;
- a[j + 9] = x0i + x2i;
- x0r -= x2r;
- x0i -= x2i;
- a[j + 12] = -wk2i * x0r - wk2r * x0i;
- a[j + 13] = -wk2i * x0i + wk2r * x0r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[j + 10] = wk1r * x0r - wk1i * x0i;
- a[j + 11] = wk1r * x0i + wk1i * x0r;
- x0r = x1r + x3i;
- x0i = x1i - x3r;
- a[j + 14] = wk3r * x0r - wk3i * x0i;
- a[j + 15] = wk3r * x0i + wk3i * x0r;
- }
- }
-
-
- void cftmdl(int n, int l, double *a, double *w)
- {
- int j, j1, j2, j3, k, k1, k2, m, m2;
- double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
- double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
- m = l << 2;
- for (j = 0; j < l; j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = a[j + 1] + a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = a[j + 1] - a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- a[j2] = x0r - x2r;
- a[j2 + 1] = x0i - x2i;
- a[j1] = x1r - x3i;
- a[j1 + 1] = x1i + x3r;
- a[j3] = x1r + x3i;
- a[j3 + 1] = x1i - x3r;
- }
- wk1r = w[2];
- for (j = m; j < l + m; j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = a[j + 1] + a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = a[j + 1] - a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- a[j2] = x2i - x0i;
- a[j2 + 1] = x0r - x2r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[j1] = wk1r * (x0r - x0i);
- a[j1 + 1] = wk1r * (x0r + x0i);
- x0r = x3i + x1r;
- x0i = x3r - x1i;
- a[j3] = wk1r * (x0i - x0r);
- a[j3 + 1] = wk1r * (x0i + x0r);
- }
- k1 = 0;
- m2 = 2 * m;
- for (k = m2; k < n; k += m2) {
- k1 += 2;
- k2 = 2 * k1;
- wk2r = w[k1];
- wk2i = w[k1 + 1];
- wk1r = w[k2];
- wk1i = w[k2 + 1];
- wk3r = wk1r - 2 * wk2i * wk1i;
- wk3i = 2 * wk2i * wk1r - wk1i;
- for (j = k; j < l + k; j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = a[j + 1] + a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = a[j + 1] - a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- x0r -= x2r;
- x0i -= x2i;
- a[j2] = wk2r * x0r - wk2i * x0i;
- a[j2 + 1] = wk2r * x0i + wk2i * x0r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[j1] = wk1r * x0r - wk1i * x0i;
- a[j1 + 1] = wk1r * x0i + wk1i * x0r;
- x0r = x1r + x3i;
- x0i = x1i - x3r;
- a[j3] = wk3r * x0r - wk3i * x0i;
- a[j3 + 1] = wk3r * x0i + wk3i * x0r;
- }
- wk1r = w[k2 + 2];
- wk1i = w[k2 + 3];
- wk3r = wk1r - 2 * wk2r * wk1i;
- wk3i = 2 * wk2r * wk1r - wk1i;
- for (j = k + m; j < l + (k + m); j += 2) {
- j1 = j + l;
- j2 = j1 + l;
- j3 = j2 + l;
- x0r = a[j] + a[j1];
- x0i = a[j + 1] + a[j1 + 1];
- x1r = a[j] - a[j1];
- x1i = a[j + 1] - a[j1 + 1];
- x2r = a[j2] + a[j3];
- x2i = a[j2 + 1] + a[j3 + 1];
- x3r = a[j2] - a[j3];
- x3i = a[j2 + 1] - a[j3 + 1];
- a[j] = x0r + x2r;
- a[j + 1] = x0i + x2i;
- x0r -= x2r;
- x0i -= x2i;
- a[j2] = -wk2i * x0r - wk2r * x0i;
- a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
- x0r = x1r - x3i;
- x0i = x1i + x3r;
- a[j1] = wk1r * x0r - wk1i * x0i;
- a[j1 + 1] = wk1r * x0i + wk1i * x0r;
- x0r = x1r + x3i;
- x0i = x1i - x3r;
- a[j3] = wk3r * x0r - wk3i * x0i;
- a[j3 + 1] = wk3r * x0i + wk3i * x0r;
- }
- }
- }
-
-
- void rftfsub(int n, double *a, int nc, double *c)
- {
- int j, k, kk, ks, m;
- double wkr, wki, xr, xi, yr, yi;
-
- m = n >> 1;
- ks = 2 * nc / m;
- kk = 0;
- for (j = 2; j < m; j += 2) {
- k = n - j;
- kk += ks;
- wkr = 0.5 - c[nc - kk];
- wki = c[kk];
- xr = a[j] - a[k];
- xi = a[j + 1] + a[k + 1];
- yr = wkr * xr - wki * xi;
- yi = wkr * xi + wki * xr;
- a[j] -= yr;
- a[j + 1] -= yi;
- a[k] += yr;
- a[k + 1] -= yi;
- }
- }
-
-
- void rftbsub(int n, double *a, int nc, double *c)
- {
- int j, k, kk, ks, m;
- double wkr, wki, xr, xi, yr, yi;
-
- a[1] = -a[1];
- m = n >> 1;
- ks = 2 * nc / m;
- kk = 0;
- for (j = 2; j < m; j += 2) {
- k = n - j;
- kk += ks;
- wkr = 0.5 - c[nc - kk];
- wki = c[kk];
- xr = a[j] - a[k];
- xi = a[j + 1] + a[k + 1];
- yr = wkr * xr + wki * xi;
- yi = wkr * xi - wki * xr;
- a[j] -= yr;
- a[j + 1] = yi - a[j + 1];
- a[k] += yr;
- a[k + 1] = yi - a[k + 1];
- }
- a[m + 1] = -a[m + 1];
- }
-
- OouraFFT(const OouraFFT&) = delete;
- OouraFFT& operator=(const OouraFFT&) = delete;
- };
-
- std::unique_ptr<AudioFFTImpl> MakeAudioFFTImpl()
- {
- return std::unique_ptr<OouraFFT>(new OouraFFT());
- }
-
-
-#endif // AUDIOFFT_OOURA_USED
-
-
- // ================================================================
-
-
-#ifdef AUDIOFFT_APPLE_ACCELERATE_USED
-
-
- /**
- * @internal
- * @class AppleAccelerateFFT
- * @brief FFT implementation using the Apple Accelerate framework internally
- */
- class AppleAccelerateFFT : public AudioFFTImpl
- {
- public:
- AppleAccelerateFFT() :
- AudioFFTImpl(),
- _size(0),
- _powerOf2(0),
- _fftSetup(0),
- _re(),
- _im()
- {
- }
-
- virtual ~AppleAccelerateFFT()
- {
- init(0);
- }
-
- virtual void init(size_t size) override
- {
- if (_fftSetup)
- {
- vDSP_destroy_fftsetup(_fftSetup);
- _size = 0;
- _powerOf2 = 0;
- _fftSetup = 0;
- _re.clear();
- _im.clear();
- }
-
- if (size > 0)
- {
- _size = size;
- _powerOf2 = 0;
- while ((1 << _powerOf2) < _size)
- {
- ++_powerOf2;
- }
- _fftSetup = vDSP_create_fftsetup(_powerOf2, FFT_RADIX2);
- _re.resize(_size / 2);
- _im.resize(_size / 2);
- }
- }
-
- virtual void fft(const float* data, float* re, float* im) override
- {
- const size_t size2 = _size / 2;
- DSPSplitComplex splitComplex;
- splitComplex.realp = re;
- splitComplex.imagp = im;
- vDSP_ctoz(reinterpret_cast<const COMPLEX*>(data), 2, &splitComplex, 1, size2);
- vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_FORWARD);
- const float factor = 0.5f;
- vDSP_vsmul(re, 1, &factor, re, 1, size2);
- vDSP_vsmul(im, 1, &factor, im, 1, size2);
- re[size2] = im[0];
- im[0] = 0.0f;
- im[size2] = 0.0f;
- }
-
- virtual void ifft(float* data, const float* re, const float* im) override
- {
- const size_t size2 = _size / 2;
- ::memcpy(_re.data(), re, size2 * sizeof(float));
- ::memcpy(_im.data(), im, size2 * sizeof(float));
- _im[0] = re[size2];
- DSPSplitComplex splitComplex;
- splitComplex.realp = _re.data();
- splitComplex.imagp = _im.data();
- vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_INVERSE);
- vDSP_ztoc(&splitComplex, 1, reinterpret_cast<COMPLEX*>(data), 2, size2);
- const float factor = 1.0f / static_cast<float>(_size);
- vDSP_vsmul(data, 1, &factor, data, 1, _size);
- }
-
- private:
- size_t _size;
- size_t _powerOf2;
- FFTSetup _fftSetup;
- std::vector<float> _re;
- std::vector<float> _im;
-
- AppleAccelerateFFT(const AppleAccelerateFFT&) = delete;
- AppleAccelerateFFT& operator=(const AppleAccelerateFFT&) = delete;
- };
-
-
- std::unique_ptr<AudioFFTImpl> MakeAudioFFTImpl()
- {
- return std::unique_ptr<AppleAccelerateFFT>(new AppleAccelerateFFT());
- }
-
-
-#endif // AUDIOFFT_APPLE_ACCELERATE_USED
-
-
- // ================================================================
-
-
-#ifdef AUDIOFFT_FFTW3_USED
-
-
- /**
- * @internal
- * @class FFTW3FFT
- * @brief FFT implementation using FFTW3 internally (see fftw.org)
- */
- class FFTW3FFT : public AudioFFTImpl
- {
- public:
- FFTW3FFT() :
- AudioFFTImpl(),
- _size(0),
- _complexSize(0),
- _planForward(0),
- _planBackward(0),
- _data(0),
- _re(0),
- _im(0)
- {
- }
-
- virtual ~FFTW3FFT()
- {
- init(0);
- }
-
- virtual void init(size_t size) override
- {
- if (_size != size)
- {
- if (_size > 0)
- {
- fftwf_destroy_plan(_planForward);
- fftwf_destroy_plan(_planBackward);
- _planForward = 0;
- _planBackward = 0;
- _size = 0;
- _complexSize = 0;
-
- if (_data)
- {
- fftwf_free(_data);
- _data = 0;
- }
-
- if (_re)
- {
- fftwf_free(_re);
- _re = 0;
- }
-
- if (_im)
- {
- fftwf_free(_im);
- _im = 0;
- }
- }
-
- if (size > 0)
- {
- _size = size;
- _complexSize = AudioFFT::ComplexSize(_size);
- const size_t complexSize = AudioFFT::ComplexSize(_size);
- _data = reinterpret_cast<float*>(fftwf_malloc(_size * sizeof(float)));
- _re = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));
- _im = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));
-
- fftw_iodim dim;
- dim.n = static_cast<int>(size);
- dim.is = 1;
- dim.os = 1;
- _planForward = fftwf_plan_guru_split_dft_r2c(1, &dim, 0, 0, _data, _re, _im, FFTW_MEASURE);
- _planBackward = fftwf_plan_guru_split_dft_c2r(1, &dim, 0, 0, _re, _im, _data, FFTW_MEASURE);
- }
- }
- }
-
- virtual void fft(const float* data, float* re, float* im) override
- {
- ::memcpy(_data, data, _size * sizeof(float));
- fftwf_execute_split_dft_r2c(_planForward, _data, _re, _im);
- ::memcpy(re, _re, _complexSize * sizeof(float));
- ::memcpy(im, _im, _complexSize * sizeof(float));
- }
-
- void ifft(float* data, const float* re, const float* im)
- {
- ::memcpy(_re, re, _complexSize * sizeof(float));
- ::memcpy(_im, im, _complexSize * sizeof(float));
- fftwf_execute_split_dft_c2r(_planBackward, _re, _im, _data);
- ScaleBuffer(data, _data, 1.0f / static_cast<float>(_size), _size);
- }
-
- private:
- size_t _size;
- size_t _complexSize;
- fftwf_plan _planForward;
- fftwf_plan _planBackward;
- float* _data;
- float* _re;
- float* _im;
-
- FFTW3FFT(const FFTW3FFT&) = delete;
- FFTW3FFT& operator=(const FFTW3FFT&) = delete;
- };
-
-
- std::unique_ptr<AudioFFTImpl> MakeAudioFFTImpl()
- {
- return std::unique_ptr<FFTW3FFT>(new FFTW3FFT());
- }
-
-
-#endif // AUDIOFFT_FFTW3_USED
-
- } // End of namespace details
-
-
- // =============================================================
-
-
- AudioFFT::AudioFFT() :
- _impl(details::MakeAudioFFTImpl())
- {
- }
-
-
- void AudioFFT::init(size_t size)
- {
- assert(details::IsPowerOf2(size));
- _impl->init(size);
- }
-
-
- void AudioFFT::fft(const float* data, float* re, float* im)
- {
- _impl->fft(data, re, im);
- }
-
-
- void AudioFFT::ifft(float* data, const float* re, const float* im)
- {
- _impl->ifft(data, re, im);
- }
-
-
- size_t AudioFFT::ComplexSize(size_t size)
- {
- return (size / 2) + 1;
- }
-
-} // End of namespace
+// ==================================================================================\r
+// Copyright (c) 2017 HiFi-LoFi\r
+//\r
+// Permission is hereby granted, free of charge, to any person obtaining a copy\r
+// of this software and associated documentation files (the "Software"), to deal\r
+// in the Software without restriction, including without limitation the rights\r
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\r
+// copies of the Software, and to permit persons to whom the Software is furnished\r
+// to do so, subject to the following conditions:\r
+//\r
+// The above copyright notice and this permission notice shall be included in\r
+// all copies or substantial portions of the Software.\r
+//\r
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\r
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS\r
+// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR\r
+// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER\r
+// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION\r
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
+// ==================================================================================\r
+\r
+#include "AudioFFT.h"\r
+\r
+#include <cassert>\r
+#include <cmath>\r
+#include <cstring>\r
+\r
+\r
+#if defined(AUDIOFFT_APPLE_ACCELERATE)\r
+ #define AUDIOFFT_APPLE_ACCELERATE_USED\r
+ #include <Accelerate/Accelerate.h>\r
+ #include <vector>\r
+#elif defined (AUDIOFFT_FFTW3)\r
+ #define AUDIOFFT_FFTW3_USED\r
+ #include <fftw3.h>\r
+#else\r
+ #if !defined(AUDIOFFT_OOURA)\r
+ #define AUDIOFFT_OOURA\r
+ #endif\r
+ #define AUDIOFFT_OOURA_USED\r
+ #include <vector>\r
+#endif\r
+\r
+\r
+namespace audiofft\r
+{\r
+\r
+ namespace detail\r
+ {\r
+\r
+ class AudioFFTImpl\r
+ {\r
+ public:\r
+ AudioFFTImpl() = default;\r
+ AudioFFTImpl(const AudioFFTImpl&) = delete;\r
+ AudioFFTImpl& operator=(const AudioFFTImpl&) = delete;\r
+ virtual ~AudioFFTImpl() = default;\r
+ virtual void init(size_t size) = 0;\r
+ virtual void fft(const float* data, float* re, float* im) = 0;\r
+ virtual void ifft(float* data, const float* re, const float* im) = 0;\r
+ };\r
+\r
+\r
+ constexpr bool IsPowerOf2(size_t val)\r
+ {\r
+ return (val == 1 || (val & (val-1)) == 0);\r
+ }\r
+\r
+\r
+ template<typename TypeDest, typename TypeSrc>\r
+ void ConvertBuffer(TypeDest* dest, const TypeSrc* src, size_t len)\r
+ {\r
+ for (size_t i=0; i<len; ++i)\r
+ {\r
+ dest[i] = static_cast<TypeDest>(src[i]);\r
+ }\r
+ }\r
+\r
+\r
+ template<typename TypeDest, typename TypeSrc, typename TypeFactor>\r
+ void ScaleBuffer(TypeDest* dest, const TypeSrc* src, const TypeFactor factor, size_t len)\r
+ {\r
+ for (size_t i=0; i<len; ++i)\r
+ {\r
+ dest[i] = static_cast<TypeDest>(static_cast<TypeFactor>(src[i]) * factor);\r
+ }\r
+ }\r
+\r
+ } // End of namespace detail\r
+\r
+\r
+ // ================================================================\r
+\r
+\r
+#ifdef AUDIOFFT_OOURA_USED\r
+\r
+ /**\r
+ * @internal\r
+ * @class OouraFFT\r
+ * @brief FFT implementation based on the great radix-4 routines by Takuya Ooura\r
+ */\r
+ class OouraFFT : public detail::AudioFFTImpl\r
+ {\r
+ public:\r
+ OouraFFT() :\r
+ detail::AudioFFTImpl(),\r
+ _size(0),\r
+ _ip(),\r
+ _w(),\r
+ _buffer()\r
+ {\r
+ }\r
+\r
+ OouraFFT(const OouraFFT&) = delete;\r
+ OouraFFT& operator=(const OouraFFT&) = delete;\r
+\r
+ virtual void init(size_t size) override\r
+ {\r
+ if (_size != size)\r
+ {\r
+ _ip.resize(2 + static_cast<int>(std::sqrt(static_cast<double>(size))));\r
+ _w.resize(size / 2);\r
+ _buffer.resize(size);\r
+ _size = size;\r
+\r
+ const int size4 = static_cast<int>(_size) / 4;\r
+ makewt(size4, _ip.data(), _w.data());\r
+ makect(size4, _ip.data(), _w.data() + size4);\r
+ }\r
+ }\r
+\r
+ virtual void fft(const float* data, float* re, float* im) override\r
+ {\r
+ // Convert into the format as required by the Ooura FFT\r
+ detail::ConvertBuffer(_buffer.data(), data, _size);\r
+\r
+ rdft(static_cast<int>(_size), +1, _buffer.data(), _ip.data(), _w.data());\r
+\r
+ // Convert back to split-complex\r
+ {\r
+ double* b = _buffer.data();\r
+ double* bEnd = b + _size;\r
+ float *r = re;\r
+ float *i = im;\r
+ while (b != bEnd)\r
+ {\r
+ *(r++) = static_cast<float>(*(b++));\r
+ *(i++) = static_cast<float>(-(*(b++)));\r
+ }\r
+ }\r
+ const size_t size2 = _size / 2;\r
+ re[size2] = -im[0];\r
+ im[0] = 0.0;\r
+ im[size2] = 0.0;\r
+ }\r
+\r
+ virtual void ifft(float* data, const float* re, const float* im) override\r
+ {\r
+ // Convert into the format as required by the Ooura FFT\r
+ {\r
+ double* b = _buffer.data();\r
+ double* bEnd = b + _size;\r
+ const float *r = re;\r
+ const float *i = im;\r
+ while (b != bEnd)\r
+ {\r
+ *(b++) = static_cast<double>(*(r++));\r
+ *(b++) = -static_cast<double>(*(i++));\r
+ }\r
+ _buffer[1] = re[_size / 2];\r
+ }\r
+\r
+ rdft(static_cast<int>(_size), -1, _buffer.data(), _ip.data(), _w.data());\r
+\r
+ // Convert back to split-complex\r
+ detail::ScaleBuffer(data, _buffer.data(), 2.0 / static_cast<double>(_size), _size);\r
+ }\r
+\r
+ private:\r
+ size_t _size;\r
+ std::vector<int> _ip;\r
+ std::vector<double> _w;\r
+ std::vector<double> _buffer;\r
+\r
+ void rdft(int n, int isgn, double *a, int *ip, double *w)\r
+ {\r
+ int nw = ip[0];\r
+ int nc = ip[1];\r
+\r
+ if (isgn >= 0)\r
+ {\r
+ if (n > 4)\r
+ {\r
+ bitrv2(n, ip + 2, a);\r
+ cftfsub(n, a, w);\r
+ rftfsub(n, a, nc, w + nw);\r
+ }\r
+ else if (n == 4)\r
+ {\r
+ cftfsub(n, a, w);\r
+ }\r
+ double xi = a[0] - a[1];\r
+ a[0] += a[1];\r
+ a[1] = xi;\r
+ }\r
+ else\r
+ {\r
+ a[1] = 0.5 * (a[0] - a[1]);\r
+ a[0] -= a[1];\r
+ if (n > 4)\r
+ {\r
+ rftbsub(n, a, nc, w + nw);\r
+ bitrv2(n, ip + 2, a);\r
+ cftbsub(n, a, w);\r
+ }\r
+ else if (n == 4)\r
+ {\r
+ cftfsub(n, a, w);\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ /* -------- initializing routines -------- */\r
+\r
+ void makewt(int nw, int *ip, double *w)\r
+ {\r
+ int j, nwh;\r
+ double delta, x, y;\r
+\r
+ ip[0] = nw;\r
+ ip[1] = 1;\r
+ if (nw > 2) {\r
+ nwh = nw >> 1;\r
+ delta = atan(1.0) / nwh;\r
+ w[0] = 1;\r
+ w[1] = 0;\r
+ w[nwh] = cos(delta * nwh);\r
+ w[nwh + 1] = w[nwh];\r
+ if (nwh > 2) {\r
+ for (j = 2; j < nwh; j += 2) {\r
+ x = cos(delta * j);\r
+ y = sin(delta * j);\r
+ w[j] = x;\r
+ w[j + 1] = y;\r
+ w[nw - j] = y;\r
+ w[nw - j + 1] = x;\r
+ }\r
+ bitrv2(nw, ip + 2, w);\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ void makect(int nc, int *ip, double *c)\r
+ {\r
+ int j, nch;\r
+ double delta;\r
+\r
+ ip[1] = nc;\r
+ if (nc > 1) {\r
+ nch = nc >> 1;\r
+ delta = atan(1.0) / nch;\r
+ c[0] = cos(delta * nch);\r
+ c[nch] = 0.5 * c[0];\r
+ for (j = 1; j < nch; j++) {\r
+ c[j] = 0.5 * cos(delta * j);\r
+ c[nc - j] = 0.5 * sin(delta * j);\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ /* -------- child routines -------- */\r
+\r
+\r
+ void bitrv2(int n, int *ip, double *a)\r
+ {\r
+ int j, j1, k, k1, l, m, m2;\r
+ double xr, xi, yr, yi;\r
+\r
+ ip[0] = 0;\r
+ l = n;\r
+ m = 1;\r
+ while ((m << 3) < l) {\r
+ l >>= 1;\r
+ for (j = 0; j < m; j++) {\r
+ ip[m + j] = ip[j] + l;\r
+ }\r
+ m <<= 1;\r
+ }\r
+ m2 = 2 * m;\r
+ if ((m << 3) == l) {\r
+ for (k = 0; k < m; k++) {\r
+ for (j = 0; j < k; j++) {\r
+ j1 = 2 * j + ip[k];\r
+ k1 = 2 * k + ip[j];\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ j1 += m2;\r
+ k1 += 2 * m2;\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ j1 += m2;\r
+ k1 -= m2;\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ j1 += m2;\r
+ k1 += 2 * m2;\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ }\r
+ j1 = 2 * k + m2 + ip[k];\r
+ k1 = j1 + m2;\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ }\r
+ } else {\r
+ for (k = 1; k < m; k++) {\r
+ for (j = 0; j < k; j++) {\r
+ j1 = 2 * j + ip[k];\r
+ k1 = 2 * k + ip[j];\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ j1 += m2;\r
+ k1 += m2;\r
+ xr = a[j1];\r
+ xi = a[j1 + 1];\r
+ yr = a[k1];\r
+ yi = a[k1 + 1];\r
+ a[j1] = yr;\r
+ a[j1 + 1] = yi;\r
+ a[k1] = xr;\r
+ a[k1 + 1] = xi;\r
+ }\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ void cftfsub(int n, double *a, double *w)\r
+ {\r
+ int j, j1, j2, j3, l;\r
+ double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;\r
+\r
+ l = 2;\r
+ if (n > 8) {\r
+ cft1st(n, a, w);\r
+ l = 8;\r
+ while ((l << 2) < n) {\r
+ cftmdl(n, l, a, w);\r
+ l <<= 2;\r
+ }\r
+ }\r
+ if ((l << 2) == n) {\r
+ for (j = 0; j < l; j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = a[j + 1] + a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = a[j + 1] - a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ a[j2] = x0r - x2r;\r
+ a[j2 + 1] = x0i - x2i;\r
+ a[j1] = x1r - x3i;\r
+ a[j1 + 1] = x1i + x3r;\r
+ a[j3] = x1r + x3i;\r
+ a[j3 + 1] = x1i - x3r;\r
+ }\r
+ } else {\r
+ for (j = 0; j < l; j += 2) {\r
+ j1 = j + l;\r
+ x0r = a[j] - a[j1];\r
+ x0i = a[j + 1] - a[j1 + 1];\r
+ a[j] += a[j1];\r
+ a[j + 1] += a[j1 + 1];\r
+ a[j1] = x0r;\r
+ a[j1 + 1] = x0i;\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ void cftbsub(int n, double *a, double *w)\r
+ {\r
+ int j, j1, j2, j3, l;\r
+ double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;\r
+\r
+ l = 2;\r
+ if (n > 8) {\r
+ cft1st(n, a, w);\r
+ l = 8;\r
+ while ((l << 2) < n) {\r
+ cftmdl(n, l, a, w);\r
+ l <<= 2;\r
+ }\r
+ }\r
+ if ((l << 2) == n) {\r
+ for (j = 0; j < l; j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = -a[j + 1] - a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = -a[j + 1] + a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i - x2i;\r
+ a[j2] = x0r - x2r;\r
+ a[j2 + 1] = x0i + x2i;\r
+ a[j1] = x1r - x3i;\r
+ a[j1 + 1] = x1i - x3r;\r
+ a[j3] = x1r + x3i;\r
+ a[j3 + 1] = x1i + x3r;\r
+ }\r
+ } else {\r
+ for (j = 0; j < l; j += 2) {\r
+ j1 = j + l;\r
+ x0r = a[j] - a[j1];\r
+ x0i = -a[j + 1] + a[j1 + 1];\r
+ a[j] += a[j1];\r
+ a[j + 1] = -a[j + 1] - a[j1 + 1];\r
+ a[j1] = x0r;\r
+ a[j1 + 1] = x0i;\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ void cft1st(int n, double *a, double *w)\r
+ {\r
+ int j, k1, k2;\r
+ double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;\r
+ double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;\r
+\r
+ x0r = a[0] + a[2];\r
+ x0i = a[1] + a[3];\r
+ x1r = a[0] - a[2];\r
+ x1i = a[1] - a[3];\r
+ x2r = a[4] + a[6];\r
+ x2i = a[5] + a[7];\r
+ x3r = a[4] - a[6];\r
+ x3i = a[5] - a[7];\r
+ a[0] = x0r + x2r;\r
+ a[1] = x0i + x2i;\r
+ a[4] = x0r - x2r;\r
+ a[5] = x0i - x2i;\r
+ a[2] = x1r - x3i;\r
+ a[3] = x1i + x3r;\r
+ a[6] = x1r + x3i;\r
+ a[7] = x1i - x3r;\r
+ wk1r = w[2];\r
+ x0r = a[8] + a[10];\r
+ x0i = a[9] + a[11];\r
+ x1r = a[8] - a[10];\r
+ x1i = a[9] - a[11];\r
+ x2r = a[12] + a[14];\r
+ x2i = a[13] + a[15];\r
+ x3r = a[12] - a[14];\r
+ x3i = a[13] - a[15];\r
+ a[8] = x0r + x2r;\r
+ a[9] = x0i + x2i;\r
+ a[12] = x2i - x0i;\r
+ a[13] = x0r - x2r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[10] = wk1r * (x0r - x0i);\r
+ a[11] = wk1r * (x0r + x0i);\r
+ x0r = x3i + x1r;\r
+ x0i = x3r - x1i;\r
+ a[14] = wk1r * (x0i - x0r);\r
+ a[15] = wk1r * (x0i + x0r);\r
+ k1 = 0;\r
+ for (j = 16; j < n; j += 16) {\r
+ k1 += 2;\r
+ k2 = 2 * k1;\r
+ wk2r = w[k1];\r
+ wk2i = w[k1 + 1];\r
+ wk1r = w[k2];\r
+ wk1i = w[k2 + 1];\r
+ wk3r = wk1r - 2 * wk2i * wk1i;\r
+ wk3i = 2 * wk2i * wk1r - wk1i;\r
+ x0r = a[j] + a[j + 2];\r
+ x0i = a[j + 1] + a[j + 3];\r
+ x1r = a[j] - a[j + 2];\r
+ x1i = a[j + 1] - a[j + 3];\r
+ x2r = a[j + 4] + a[j + 6];\r
+ x2i = a[j + 5] + a[j + 7];\r
+ x3r = a[j + 4] - a[j + 6];\r
+ x3i = a[j + 5] - a[j + 7];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ x0r -= x2r;\r
+ x0i -= x2i;\r
+ a[j + 4] = wk2r * x0r - wk2i * x0i;\r
+ a[j + 5] = wk2r * x0i + wk2i * x0r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[j + 2] = wk1r * x0r - wk1i * x0i;\r
+ a[j + 3] = wk1r * x0i + wk1i * x0r;\r
+ x0r = x1r + x3i;\r
+ x0i = x1i - x3r;\r
+ a[j + 6] = wk3r * x0r - wk3i * x0i;\r
+ a[j + 7] = wk3r * x0i + wk3i * x0r;\r
+ wk1r = w[k2 + 2];\r
+ wk1i = w[k2 + 3];\r
+ wk3r = wk1r - 2 * wk2r * wk1i;\r
+ wk3i = 2 * wk2r * wk1r - wk1i;\r
+ x0r = a[j + 8] + a[j + 10];\r
+ x0i = a[j + 9] + a[j + 11];\r
+ x1r = a[j + 8] - a[j + 10];\r
+ x1i = a[j + 9] - a[j + 11];\r
+ x2r = a[j + 12] + a[j + 14];\r
+ x2i = a[j + 13] + a[j + 15];\r
+ x3r = a[j + 12] - a[j + 14];\r
+ x3i = a[j + 13] - a[j + 15];\r
+ a[j + 8] = x0r + x2r;\r
+ a[j + 9] = x0i + x2i;\r
+ x0r -= x2r;\r
+ x0i -= x2i;\r
+ a[j + 12] = -wk2i * x0r - wk2r * x0i;\r
+ a[j + 13] = -wk2i * x0i + wk2r * x0r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[j + 10] = wk1r * x0r - wk1i * x0i;\r
+ a[j + 11] = wk1r * x0i + wk1i * x0r;\r
+ x0r = x1r + x3i;\r
+ x0i = x1i - x3r;\r
+ a[j + 14] = wk3r * x0r - wk3i * x0i;\r
+ a[j + 15] = wk3r * x0i + wk3i * x0r;\r
+ }\r
+ }\r
+\r
+\r
+ void cftmdl(int n, int l, double *a, double *w)\r
+ {\r
+ int j, j1, j2, j3, k, k1, k2, m, m2;\r
+ double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;\r
+ double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;\r
+\r
+ m = l << 2;\r
+ for (j = 0; j < l; j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = a[j + 1] + a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = a[j + 1] - a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ a[j2] = x0r - x2r;\r
+ a[j2 + 1] = x0i - x2i;\r
+ a[j1] = x1r - x3i;\r
+ a[j1 + 1] = x1i + x3r;\r
+ a[j3] = x1r + x3i;\r
+ a[j3 + 1] = x1i - x3r;\r
+ }\r
+ wk1r = w[2];\r
+ for (j = m; j < l + m; j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = a[j + 1] + a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = a[j + 1] - a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ a[j2] = x2i - x0i;\r
+ a[j2 + 1] = x0r - x2r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[j1] = wk1r * (x0r - x0i);\r
+ a[j1 + 1] = wk1r * (x0r + x0i);\r
+ x0r = x3i + x1r;\r
+ x0i = x3r - x1i;\r
+ a[j3] = wk1r * (x0i - x0r);\r
+ a[j3 + 1] = wk1r * (x0i + x0r);\r
+ }\r
+ k1 = 0;\r
+ m2 = 2 * m;\r
+ for (k = m2; k < n; k += m2) {\r
+ k1 += 2;\r
+ k2 = 2 * k1;\r
+ wk2r = w[k1];\r
+ wk2i = w[k1 + 1];\r
+ wk1r = w[k2];\r
+ wk1i = w[k2 + 1];\r
+ wk3r = wk1r - 2 * wk2i * wk1i;\r
+ wk3i = 2 * wk2i * wk1r - wk1i;\r
+ for (j = k; j < l + k; j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = a[j + 1] + a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = a[j + 1] - a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ x0r -= x2r;\r
+ x0i -= x2i;\r
+ a[j2] = wk2r * x0r - wk2i * x0i;\r
+ a[j2 + 1] = wk2r * x0i + wk2i * x0r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[j1] = wk1r * x0r - wk1i * x0i;\r
+ a[j1 + 1] = wk1r * x0i + wk1i * x0r;\r
+ x0r = x1r + x3i;\r
+ x0i = x1i - x3r;\r
+ a[j3] = wk3r * x0r - wk3i * x0i;\r
+ a[j3 + 1] = wk3r * x0i + wk3i * x0r;\r
+ }\r
+ wk1r = w[k2 + 2];\r
+ wk1i = w[k2 + 3];\r
+ wk3r = wk1r - 2 * wk2r * wk1i;\r
+ wk3i = 2 * wk2r * wk1r - wk1i;\r
+ for (j = k + m; j < l + (k + m); j += 2) {\r
+ j1 = j + l;\r
+ j2 = j1 + l;\r
+ j3 = j2 + l;\r
+ x0r = a[j] + a[j1];\r
+ x0i = a[j + 1] + a[j1 + 1];\r
+ x1r = a[j] - a[j1];\r
+ x1i = a[j + 1] - a[j1 + 1];\r
+ x2r = a[j2] + a[j3];\r
+ x2i = a[j2 + 1] + a[j3 + 1];\r
+ x3r = a[j2] - a[j3];\r
+ x3i = a[j2 + 1] - a[j3 + 1];\r
+ a[j] = x0r + x2r;\r
+ a[j + 1] = x0i + x2i;\r
+ x0r -= x2r;\r
+ x0i -= x2i;\r
+ a[j2] = -wk2i * x0r - wk2r * x0i;\r
+ a[j2 + 1] = -wk2i * x0i + wk2r * x0r;\r
+ x0r = x1r - x3i;\r
+ x0i = x1i + x3r;\r
+ a[j1] = wk1r * x0r - wk1i * x0i;\r
+ a[j1 + 1] = wk1r * x0i + wk1i * x0r;\r
+ x0r = x1r + x3i;\r
+ x0i = x1i - x3r;\r
+ a[j3] = wk3r * x0r - wk3i * x0i;\r
+ a[j3 + 1] = wk3r * x0i + wk3i * x0r;\r
+ }\r
+ }\r
+ }\r
+\r
+\r
+ void rftfsub(int n, double *a, int nc, double *c)\r
+ {\r
+ int j, k, kk, ks, m;\r
+ double wkr, wki, xr, xi, yr, yi;\r
+\r
+ m = n >> 1;\r
+ ks = 2 * nc / m;\r
+ kk = 0;\r
+ for (j = 2; j < m; j += 2) {\r
+ k = n - j;\r
+ kk += ks;\r
+ wkr = 0.5 - c[nc - kk];\r
+ wki = c[kk];\r
+ xr = a[j] - a[k];\r
+ xi = a[j + 1] + a[k + 1];\r
+ yr = wkr * xr - wki * xi;\r
+ yi = wkr * xi + wki * xr;\r
+ a[j] -= yr;\r
+ a[j + 1] -= yi;\r
+ a[k] += yr;\r
+ a[k + 1] -= yi;\r
+ }\r
+ }\r
+\r
+\r
+ void rftbsub(int n, double *a, int nc, double *c)\r
+ {\r
+ int j, k, kk, ks, m;\r
+ double wkr, wki, xr, xi, yr, yi;\r
+\r
+ a[1] = -a[1];\r
+ m = n >> 1;\r
+ ks = 2 * nc / m;\r
+ kk = 0;\r
+ for (j = 2; j < m; j += 2) {\r
+ k = n - j;\r
+ kk += ks;\r
+ wkr = 0.5 - c[nc - kk];\r
+ wki = c[kk];\r
+ xr = a[j] - a[k];\r
+ xi = a[j + 1] + a[k + 1];\r
+ yr = wkr * xr + wki * xi;\r
+ yi = wkr * xi - wki * xr;\r
+ a[j] -= yr;\r
+ a[j + 1] = yi - a[j + 1];\r
+ a[k] += yr;\r
+ a[k + 1] = yi - a[k + 1];\r
+ }\r
+ a[m + 1] = -a[m + 1];\r
+ }\r
+ };\r
+\r
+\r
+ /**\r
+ * @internal\r
+ * @brief Concrete FFT implementation\r
+ */\r
+ typedef OouraFFT AudioFFTImplementation;\r
+\r
+\r
+#endif // AUDIOFFT_OOURA_USED\r
+\r
+\r
+ // ================================================================\r
+\r
+\r
+#ifdef AUDIOFFT_APPLE_ACCELERATE_USED\r
+\r
+\r
+ /**\r
+ * @internal\r
+ * @class AppleAccelerateFFT\r
+ * @brief FFT implementation using the Apple Accelerate framework internally\r
+ */\r
+ class AppleAccelerateFFT : public detail::AudioFFTImpl\r
+ {\r
+ public:\r
+ AppleAccelerateFFT() :\r
+ detail::AudioFFTImpl(),\r
+ _size(0),\r
+ _powerOf2(0),\r
+ _fftSetup(0),\r
+ _re(),\r
+ _im()\r
+ {\r
+ }\r
+\r
+ AppleAccelerateFFT(const AppleAccelerateFFT&) = delete;\r
+ AppleAccelerateFFT& operator=(const AppleAccelerateFFT&) = delete;\r
+\r
+ virtual ~AppleAccelerateFFT()\r
+ {\r
+ init(0);\r
+ }\r
+\r
+ virtual void init(size_t size) override\r
+ {\r
+ if (_fftSetup)\r
+ {\r
+ vDSP_destroy_fftsetup(_fftSetup);\r
+ _size = 0;\r
+ _powerOf2 = 0;\r
+ _fftSetup = 0;\r
+ _re.clear();\r
+ _im.clear();\r
+ }\r
+\r
+ if (size > 0)\r
+ {\r
+ _size = size;\r
+ _powerOf2 = 0;\r
+ while ((1 << _powerOf2) < _size)\r
+ {\r
+ ++_powerOf2;\r
+ }\r
+ _fftSetup = vDSP_create_fftsetup(_powerOf2, FFT_RADIX2);\r
+ _re.resize(_size / 2);\r
+ _im.resize(_size / 2);\r
+ }\r
+ }\r
+\r
+ virtual void fft(const float* data, float* re, float* im) override\r
+ {\r
+ const size_t size2 = _size / 2;\r
+ DSPSplitComplex splitComplex;\r
+ splitComplex.realp = re;\r
+ splitComplex.imagp = im;\r
+ vDSP_ctoz(reinterpret_cast<const COMPLEX*>(data), 2, &splitComplex, 1, size2);\r
+ vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_FORWARD);\r
+ const float factor = 0.5f;\r
+ vDSP_vsmul(re, 1, &factor, re, 1, size2);\r
+ vDSP_vsmul(im, 1, &factor, im, 1, size2);\r
+ re[size2] = im[0];\r
+ im[0] = 0.0f;\r
+ im[size2] = 0.0f;\r
+ }\r
+\r
+ virtual void ifft(float* data, const float* re, const float* im) override\r
+ {\r
+ const size_t size2 = _size / 2;\r
+ ::memcpy(_re.data(), re, size2 * sizeof(float));\r
+ ::memcpy(_im.data(), im, size2 * sizeof(float));\r
+ _im[0] = re[size2];\r
+ DSPSplitComplex splitComplex;\r
+ splitComplex.realp = _re.data();\r
+ splitComplex.imagp = _im.data();\r
+ vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_INVERSE);\r
+ vDSP_ztoc(&splitComplex, 1, reinterpret_cast<COMPLEX*>(data), 2, size2);\r
+ const float factor = 1.0f / static_cast<float>(_size);\r
+ vDSP_vsmul(data, 1, &factor, data, 1, _size);\r
+ }\r
+\r
+ private:\r
+ size_t _size;\r
+ size_t _powerOf2;\r
+ FFTSetup _fftSetup;\r
+ std::vector<float> _re;\r
+ std::vector<float> _im;\r
+ };\r
+\r
+\r
+ /**\r
+ * @internal\r
+ * @brief Concrete FFT implementation\r
+ */\r
+ typedef AppleAccelerateFFT AudioFFTImplementation;\r
+\r
+\r
+#endif // AUDIOFFT_APPLE_ACCELERATE_USED\r
+\r
+\r
+ // ================================================================\r
+\r
+\r
+#ifdef AUDIOFFT_FFTW3_USED\r
+\r
+\r
+ /**\r
+ * @internal\r
+ * @class FFTW3FFT\r
+ * @brief FFT implementation using FFTW3 internally (see fftw.org)\r
+ */\r
+ class FFTW3FFT : public detail::AudioFFTImpl\r
+ {\r
+ public:\r
+ FFTW3FFT() :\r
+ detail::AudioFFTImpl(),\r
+ _size(0),\r
+ _complexSize(0),\r
+ _planForward(0),\r
+ _planBackward(0),\r
+ _data(0),\r
+ _re(0),\r
+ _im(0)\r
+ {\r
+ }\r
+\r
+ FFTW3FFT(const FFTW3FFT&) = delete;\r
+ FFTW3FFT& operator=(const FFTW3FFT&) = delete;\r
+\r
+ virtual ~FFTW3FFT()\r
+ {\r
+ init(0);\r
+ }\r
+\r
+ virtual void init(size_t size) override\r
+ {\r
+ if (_size != size)\r
+ {\r
+ if (_size > 0)\r
+ {\r
+ fftwf_destroy_plan(_planForward);\r
+ fftwf_destroy_plan(_planBackward);\r
+ _planForward = 0;\r
+ _planBackward = 0;\r
+ _size = 0;\r
+ _complexSize = 0;\r
+\r
+ if (_data)\r
+ {\r
+ fftwf_free(_data);\r
+ _data = 0;\r
+ }\r
+\r
+ if (_re)\r
+ {\r
+ fftwf_free(_re);\r
+ _re = 0;\r
+ }\r
+\r
+ if (_im)\r
+ {\r
+ fftwf_free(_im);\r
+ _im = 0;\r
+ }\r
+ }\r
+\r
+ if (size > 0)\r
+ {\r
+ \r
+ _size = size;\r
+ _complexSize = AudioFFT::ComplexSize(_size);\r
+ const size_t complexSize = AudioFFT::ComplexSize(_size);\r
+ _data = reinterpret_cast<float*>(fftwf_malloc(_size * sizeof(float)));\r
+ _re = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));\r
+ _im = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));\r
+ fftwf_set_timelimit(0.01);\r
+ fftw_iodim dim;\r
+ dim.n = static_cast<int>(size);\r
+ dim.is = 1;\r
+ dim.os = 1;\r
+ _planForward = fftwf_plan_guru_split_dft_r2c(1, &dim, 0, 0, _data, _re, _im, FFTW_MEASURE);\r
+ _planBackward = fftwf_plan_guru_split_dft_c2r(1, &dim, 0, 0, _re, _im, _data, FFTW_MEASURE);\r
+ }\r
+ }\r
+ }\r
+\r
+ virtual void fft(const float* data, float* re, float* im) override\r
+ {\r
+ ::memcpy(_data, data, _size * sizeof(float));\r
+ fftwf_execute_split_dft_r2c(_planForward, _data, _re, _im);\r
+ ::memcpy(re, _re, _complexSize * sizeof(float));\r
+ ::memcpy(im, _im, _complexSize * sizeof(float));\r
+ }\r
+\r
+ virtual void ifft(float* data, const float* re, const float* im) override\r
+ {\r
+ ::memcpy(_re, re, _complexSize * sizeof(float));\r
+ ::memcpy(_im, im, _complexSize * sizeof(float));\r
+ fftwf_execute_split_dft_c2r(_planBackward, _re, _im, _data);\r
+ detail::ScaleBuffer(data, _data, 1.0f / static_cast<float>(_size), _size);\r
+ }\r
+\r
+ private:\r
+ size_t _size;\r
+ size_t _complexSize;\r
+ fftwf_plan _planForward;\r
+ fftwf_plan _planBackward;\r
+ float* _data;\r
+ float* _re;\r
+ float* _im;\r
+ };\r
+\r
+\r
+ /**\r
+ * @internal\r
+ * @brief Concrete FFT implementation\r
+ */\r
+ typedef FFTW3FFT AudioFFTImplementation;\r
+\r
+\r
+#endif // AUDIOFFT_FFTW3_USED\r
+\r
+\r
+ // =============================================================\r
+\r
+\r
+ AudioFFT::AudioFFT() :\r
+ _impl(new AudioFFTImplementation())\r
+ {\r
+ }\r
+\r
+\r
+ AudioFFT::~AudioFFT()\r
+ {\r
+ }\r
+\r
+\r
+ void AudioFFT::init(size_t size)\r
+ {\r
+ assert(detail::IsPowerOf2(size));\r
+ _impl->init(size);\r
+ }\r
+\r
+\r
+ void AudioFFT::fft(const float* data, float* re, float* im)\r
+ {\r
+ _impl->fft(data, re, im);\r
+ }\r
+\r
+\r
+ void AudioFFT::ifft(float* data, const float* re, const float* im)\r
+ {\r
+ _impl->ifft(data, re, im);\r
+ }\r
+\r
+\r
+ size_t AudioFFT::ComplexSize(size_t size)\r
+ {\r
+ return (size / 2) + 1;\r
+ }\r
+\r
+} // End of namespace\r
-// ==================================================================================
-// Copyright (c) 2016 HiFi-LoFi
-//
-// Permission is hereby granted, free of charge, to any person obtaining a copy
-// of this software and associated documentation files (the "Software"), to deal
-// in the Software without restriction, including without limitation the rights
-// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-// copies of the Software, and to permit persons to whom the Software is furnished
-// to do so, subject to the following conditions:
-//
-// The above copyright notice and this permission notice shall be included in
-// all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
-// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
-// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
-// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
-// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
-// ==================================================================================
-
-#ifndef _AUDIOFFT_H
-#define _AUDIOFFT_H
-
-
-/**
-* AudioFFT provides real-to-complex/complex-to-real FFT routines.
-*
-* Features:
-*
-* - Real-complex FFT and complex-real inverse FFT for power-of-2-sized real data.
-*
-* - Uniform interface to different FFT implementations (currently Ooura, FFTW3 and Apple Accelerate).
-*
-* - Complex data is handled in "split-complex" format, i.e. there are separate
-* arrays for the real and imaginary parts which can be useful for SIMD optimizations
-* (split-complex arrays have to be of length (size/2+1) representing bins from DC
-* to Nyquist frequency).
-*
-* - Output is "ready to use" (all scaling etc. is already handled internally).
-*
-* - No allocations/deallocations after the initialization which makes it usable
-* for real-time audio applications (that's what I wrote it for and using it).
-*
-*
-* How to use it in your project:
-*
-* - Add the .h and .cpp file to your project - that's all.
-*
-* - To get extra speed, you can link FFTW3 to your project and define
-* AUDIOFFT_FFTW3 (however, please check whether your project suits the
-* according license).
-*
-* - To get the best speed on Apple platforms, you can link the Apple
-* Accelerate framework to your project and define
-* AUDIOFFT_APPLE_ACCELERATE (however, please check whether your
-* project suits the according license).
-*
-*
-* Remarks:
-*
-* - AudioFFT is not intended to be the fastest FFT, but to be a fast-enough
-* FFT suitable for most audio applications.
-*
-* - AudioFFT uses the quite liberal MIT license.
-*
-*
-* Example usage:
-* @code
-* #include "AudioFFT.h"
-*
-* void Example()
-* {
-* const size_t fftSize = 1024; // Needs to be power of 2!
-*
-* std::vector<float> input(fftSize, 0.0f);
-* std::vector<float> re(audiofft::AudioFFT::ComplexSize(fftSize));
-* std::vector<float> im(audiofft::AudioFFT::ComplexSize(fftSize));
-* std::vector<float> output(fftSize);
-*
-* audiofft::AudioFFT fft;
-* fft.init(1024);
-* fft.fft(input.data(), re.data(), im.data());
-* fft.ifft(output.data(), re.data(), im.data());
-* }
-* @endcode
-*/
-
-
-#include <cstddef>
-#include <memory>
-
-
-namespace audiofft
-{
-
- namespace details
- {
-
- class AudioFFTImpl
- {
- public:
- AudioFFTImpl() = default;
- virtual ~AudioFFTImpl() = default;
- virtual void init(size_t size) = 0;
- virtual void fft(const float* data, float* re, float* im) = 0;
- virtual void ifft(float* data, const float* re, const float* im) = 0;
-
- private:
- AudioFFTImpl(const AudioFFTImpl&) = delete;
- AudioFFTImpl& operator=(const AudioFFTImpl&) = delete;
- };
- }
-
-
- // ======================================================
-
-
- /**
- * @class AudioFFT
- * @brief Performs 1D FFTs
- */
- class AudioFFT
- {
- public:
- /**
- * @brief Constructor
- */
- AudioFFT();
-
- /**
- * @brief Initializes the FFT object
- * @param size Size of the real input (must be power 2)
- */
- void init(size_t size);
-
- /**
- * @brief Performs the forward FFT
- * @param data The real input data (has to be of the length as specified in init())
- * @param re The real part of the complex output (has to be of length as returned by ComplexSize())
- * @param im The imaginary part of the complex output (has to be of length as returned by ComplexSize())
- */
- void fft(const float* data, float* re, float* im);
-
- /**
- * @brief Performs the inverse FFT
- * @param data The real output data (has to be of the length as specified in init())
- * @param re The real part of the complex input (has to be of length as returned by ComplexSize())
- * @param im The imaginary part of the complex input (has to be of length as returned by ComplexSize())
- */
- void ifft(float* data, const float* re, const float* im);
-
- /**
- * @brief Calculates the necessary size of the real/imaginary complex arrays
- * @param size The size of the real data
- * @return The size of the real/imaginary complex arrays
- */
- static size_t ComplexSize(size_t size);
-
- private:
- std::unique_ptr<details::AudioFFTImpl> _impl;
-
- AudioFFT(const AudioFFT&) = delete;
- AudioFFT& operator=(const AudioFFT&) = delete;
- };
-
-
- /**
- * @deprecated
- * @brief Let's keep an AudioFFTBase type around for now because it has been here already in the 1st version in order to avoid breaking existing code.
- */
- typedef AudioFFT AudioFFTBase;
-
-} // End of namespace
-
-#endif // Header guard
+// ==================================================================================\r
+// Copyright (c) 2017 HiFi-LoFi\r
+//\r
+// Permission is hereby granted, free of charge, to any person obtaining a copy\r
+// of this software and associated documentation files (the "Software"), to deal\r
+// in the Software without restriction, including without limitation the rights\r
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\r
+// copies of the Software, and to permit persons to whom the Software is furnished\r
+// to do so, subject to the following conditions:\r
+//\r
+// The above copyright notice and this permission notice shall be included in\r
+// all copies or substantial portions of the Software.\r
+//\r
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\r
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS\r
+// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR\r
+// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER\r
+// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION\r
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.\r
+// ==================================================================================\r
+\r
+#ifndef _AUDIOFFT_H\r
+#define _AUDIOFFT_H\r
+\r
+\r
+/**\r
+* AudioFFT provides real-to-complex/complex-to-real FFT routines.\r
+*\r
+* Features:\r
+*\r
+* - Real-complex FFT and complex-real inverse FFT for power-of-2-sized real data.\r
+*\r
+* - Uniform interface to different FFT implementations (currently Ooura, FFTW3 and Apple Accelerate).\r
+*\r
+* - Complex data is handled in "split-complex" format, i.e. there are separate\r
+* arrays for the real and imaginary parts which can be useful for SIMD optimizations\r
+* (split-complex arrays have to be of length (size/2+1) representing bins from DC\r
+* to Nyquist frequency).\r
+*\r
+* - Output is "ready to use" (all scaling etc. is already handled internally).\r
+*\r
+* - No allocations/deallocations after the initialization which makes it usable\r
+* for real-time audio applications (that's what I wrote it for and using it).\r
+*\r
+*\r
+* How to use it in your project:\r
+*\r
+* - Add the .h and .cpp file to your project - that's all.\r
+*\r
+* - To get extra speed, you can link FFTW3 to your project and define\r
+* AUDIOFFT_FFTW3 (however, please check whether your project suits the\r
+* according license).\r
+*\r
+* - To get the best speed on Apple platforms, you can link the Apple\r
+* Accelerate framework to your project and define\r
+* AUDIOFFT_APPLE_ACCELERATE (however, please check whether your\r
+* project suits the according license).\r
+*\r
+*\r
+* Remarks:\r
+*\r
+* - AudioFFT is not intended to be the fastest FFT, but to be a fast-enough\r
+* FFT suitable for most audio applications.\r
+*\r
+* - AudioFFT uses the quite liberal MIT license.\r
+*\r
+*\r
+* Example usage:\r
+* @code\r
+* #include "AudioFFT.h"\r
+*\r
+* void Example()\r
+* {\r
+* const size_t fftSize = 1024; // Needs to be power of 2!\r
+*\r
+* std::vector<float> input(fftSize, 0.0f);\r
+* std::vector<float> re(audiofft::AudioFFT::ComplexSize(fftSize));\r
+* std::vector<float> im(audiofft::AudioFFT::ComplexSize(fftSize));\r
+* std::vector<float> output(fftSize);\r
+*\r
+* audiofft::AudioFFT fft;\r
+* fft.init(1024);\r
+* fft.fft(input.data(), re.data(), im.data());\r
+* fft.ifft(output.data(), re.data(), im.data());\r
+* }\r
+* @endcode\r
+*/\r
+\r
+\r
+#include <cstddef>\r
+#include <memory>\r
+\r
+\r
+namespace audiofft\r
+{\r
+\r
+ namespace detail\r
+ {\r
+ class AudioFFTImpl;\r
+ }\r
+\r
+\r
+ // =============================================================\r
+\r
+\r
+ /**\r
+ * @class AudioFFT\r
+ * @brief Performs 1D FFTs\r
+ */\r
+ class AudioFFT\r
+ {\r
+ public:\r
+ /**\r
+ * @brief Constructor\r
+ */\r
+ AudioFFT();\r
+\r
+ AudioFFT(const AudioFFT&) = delete;\r
+ AudioFFT& operator=(const AudioFFT&) = delete;\r
+\r
+ /**\r
+ * @brief Destructor\r
+ */\r
+ ~AudioFFT();\r
+\r
+ /**\r
+ * @brief Initializes the FFT object\r
+ * @param size Size of the real input (must be power 2)\r
+ */\r
+ void init(size_t size);\r
+\r
+ /**\r
+ * @brief Performs the forward FFT\r
+ * @param data The real input data (has to be of the length as specified in init())\r
+ * @param re The real part of the complex output (has to be of length as returned by ComplexSize())\r
+ * @param im The imaginary part of the complex output (has to be of length as returned by ComplexSize())\r
+ */\r
+ void fft(const float* data, float* re, float* im);\r
+\r
+ /**\r
+ * @brief Performs the inverse FFT\r
+ * @param data The real output data (has to be of the length as specified in init())\r
+ * @param re The real part of the complex input (has to be of length as returned by ComplexSize())\r
+ * @param im The imaginary part of the complex input (has to be of length as returned by ComplexSize())\r
+ */\r
+ void ifft(float* data, const float* re, const float* im);\r
+\r
+ /**\r
+ * @brief Calculates the necessary size of the real/imaginary complex arrays\r
+ * @param size The size of the real data\r
+ * @return The size of the real/imaginary complex arrays\r
+ */\r
+ static size_t ComplexSize(size_t size);\r
+\r
+ private:\r
+ std::unique_ptr<detail::AudioFFTImpl> _impl;\r
+ };\r
+\r
+\r
+ /**\r
+ * @deprecated\r
+ * @brief Let's keep an AudioFFTBase type around for now because it has been here already in the 1st version in order to avoid breaking existing code.\r
+ */\r
+ typedef AudioFFT AudioFFTBase;\r
+\r
+} // End of namespace\r
+\r
+#endif // Header guard\r
--- /dev/null
+#include "ConvolverThreadPool.h"
+#include "FFTConvolver.h"
+#include "config.h"
+
+ConvolverThreadPool::ConvolverThreadPool()
+ : _convolvers(), _threads(), _taskQueue(), _queueMutex(), _condition(), _completionCV(),
+ _stop(false), _activeTasks(0) {}
+
+ConvolverThreadPool::~ConvolverThreadPool() {
+ shutdown();
+
+ // Delete all convolver instances
+ for (size_t i = 0; i < _convolvers.size(); ++i) {
+ delete _convolvers[i];
+ }
+ _convolvers.clear();
+}
+
+bool ConvolverThreadPool::init(size_t numThreads, size_t numConvolvers) {
+ // Clean up any existing threads first
+ shutdown();
+
+ // Validate parameters
+ if (numThreads == 0 || numConvolvers == 0) {
+ return false;
+ }
+
+ // Delete any existing convolvers
+ for (size_t i = 0; i < _convolvers.size(); ++i) {
+ delete _convolvers[i];
+ }
+ _convolvers.clear();
+
+ // Reset state
+ _stop = false;
+ _activeTasks = 0;
+
+ // Create convolver instances (but don't initialize them yet)
+ _convolvers.resize(numConvolvers);
+ for (size_t i = 0; i < numConvolvers; ++i) {
+ _convolvers[i] = new FFTConvolver();
+ }
+
+ // Create worker threads
+ _threads.reserve(numThreads);
+ for (size_t i = 0; i < numThreads; i++) {
+#ifndef COMPILE_FOR_OSX
+ pthread_setname_np(pthread_self(), "convolver");
+#endif
+ _threads.emplace_back([this]() { workerThread(); });
+ }
+
+ return true;
+}
+
+// Mutex to protect FFTW plan creation (required)
+std::mutex fftwMutex;
+
+bool ConvolverThreadPool::initConvolver(size_t convolverId, size_t blockSize, const Sample *ir,
+ size_t irLen) {
+ if (convolverId >= _convolvers.size()) {
+ return false;
+ }
+
+ // Make sure no tasks are using this convolver
+ waitForAll();
+ std::lock_guard<std::mutex> lock(fftwMutex);
+ return _convolvers[convolverId]->init(blockSize, ir, irLen);
+}
+
+bool ConvolverThreadPool::initAllConvolvers(size_t blockSize, const Sample *ir, size_t irLen) {
+ // Make sure no tasks are running
+ waitForAll();
+
+ for (size_t i = 0; i < _convolvers.size(); ++i) {
+ if (!_convolvers[i]->init(blockSize, ir, irLen)) {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+void ConvolverThreadPool::processAsync(size_t convolverId, const Sample *input, Sample *output,
+ size_t len) {
+ assert(convolverId < _convolvers.size());
+
+ {
+ std::lock_guard<std::mutex> lock(_queueMutex);
+
+ // Create the task
+ _taskQueue.push([this, convolverId, input, output, len]() {
+ _convolvers[convolverId]->process(input, output, len);
+ });
+
+ ++_activeTasks;
+ }
+
+ // Wake up one worker thread
+ _condition.notify_one();
+}
+
+void ConvolverThreadPool::waitForAll() {
+ std::unique_lock<std::mutex> lock(_queueMutex);
+ _completionCV.wait(lock, [this]() { return _taskQueue.empty() && _activeTasks == 0; });
+}
+
+void ConvolverThreadPool::clearState(size_t convolverId) {
+ assert(convolverId < _convolvers.size());
+
+ // Make sure no tasks are running before clearing state
+ waitForAll();
+
+ // _convolvers[convolverId]->clearState();
+}
+
+void ConvolverThreadPool::clearAllStates() {
+ // Make sure no tasks are running before clearing state
+ waitForAll();
+
+ for (size_t i = 0; i < _convolvers.size(); ++i) {
+ _convolvers[i]->clearState();
+ }
+}
+
+void ConvolverThreadPool::workerThread() {
+ while (true) {
+ std::function<void()> task;
+
+ // Wait for a task or stop signal
+ {
+ std::unique_lock<std::mutex> lock(_queueMutex);
+ _condition.wait(lock, [this]() { return _stop || !_taskQueue.empty(); });
+
+ // Exit if stopping and no more tasks
+ if (_stop && _taskQueue.empty()) {
+ return;
+ }
+
+ // Get the next task
+ task = std::move(_taskQueue.front());
+ _taskQueue.pop();
+ }
+
+ // Execute the task (outside the lock for better parallelism)
+ task();
+
+ // Mark task as complete
+ {
+ std::lock_guard<std::mutex> lock(_queueMutex);
+ --_activeTasks;
+ _completionCV.notify_all();
+ }
+ }
+}
+
+void ConvolverThreadPool::shutdown() {
+ // Signal threads to stop
+ {
+ std::lock_guard<std::mutex> lock(_queueMutex);
+ _stop = true;
+ }
+ _condition.notify_all();
+
+ // Wait for all threads to finish
+ for (auto &thread : _threads) {
+ if (thread.joinable()) {
+ thread.join();
+ }
+ }
+
+ // Clear thread vector
+ _threads.clear();
+
+ // Clear remaining tasks - properly this time
+ {
+ std::lock_guard<std::mutex> lock(_queueMutex);
+ while (!_taskQueue.empty()) {
+ _taskQueue.pop();
+ }
+ }
+
+ // Reset state
+ _activeTasks = 0;
+ _stop = false;
+}
\ No newline at end of file
--- /dev/null
+#ifndef CONVOLVER_THREAD_POOL_H
+#define CONVOLVER_THREAD_POOL_H
+
+#include "FFTConvolver.h"
+#include <cassert>
+#include <condition_variable>
+#include <functional>
+#include <mutex>
+#include <queue>
+#include <thread>
+#include <vector>
+
+// Use the fftconvolver namespace
+using fftconvolver::FFTConvolver;
+using fftconvolver::Sample;
+
+class ConvolverThreadPool {
+private:
+ std::vector<FFTConvolver *> _convolvers;
+ std::vector<std::thread> _threads;
+ std::queue<std::function<void()>> _taskQueue;
+ std::mutex _queueMutex;
+ std::condition_variable _condition;
+ std::condition_variable _completionCV;
+ bool _stop;
+ size_t _activeTasks;
+
+public:
+ ConvolverThreadPool();
+ ~ConvolverThreadPool();
+
+ // Initialize the thread pool (creates convolvers but doesn't initialize them)
+ // numThreads: number of worker threads (level of parallelism)
+ // numConvolvers: total number of convolver instances
+ bool init(size_t numThreads, size_t numConvolvers);
+
+ // Initialize a specific convolver with its IR
+ // convolverId: which convolver to initialize
+ // blockSize: block size for convolution
+ // ir: impulse response data
+ // irLen: length of impulse response
+ bool initConvolver(size_t convolverId, size_t blockSize, const Sample *ir, size_t irLen);
+
+ // Initialize all convolvers with the same IR
+ bool initAllConvolvers(size_t blockSize, const Sample *ir, size_t irLen);
+
+ // Queue a convolution task (non-blocking)
+ void processAsync(size_t convolverId, const Sample *input, Sample *output, size_t len);
+
+ // Wait for all queued tasks to complete (blocking)
+ void waitForAll();
+
+ // Clear the state of a specific convolver
+ void clearState(size_t convolverId);
+
+ // Clear the state of all convolvers
+ void clearAllStates();
+
+ // Get the number of convolvers
+ size_t getNumConvolvers() const { return _convolvers.size(); }
+
+ // Get the number of threads
+ size_t getNumThreads() const { return _threads.size(); }
+
+ void shutdown();
+
+private:
+ void workerThread();
+ // void shutdown();
+};
+
+#endif // CONVOLVER_THREAD_POOL_H
\ No newline at end of file
_inputBufferFill = 0;
}
+void FFTConvolver::clearState()
+{
+ if (_segCount == 0)
+ {
+ return; // Not initialized
+ }
+
+ _inputBuffer.setZero();
+ _inputBufferFill = 0;
+ _overlap.setZero();
+
+ for (size_t i = 0; i < _segCount; ++i)
+ {
+ _segments[i]->setZero();
+ }
+
+ _preMultiplied.setZero();
+ _conv.setZero();
+ _current = 0;
+}
bool FFTConvolver::init(size_t blockSize, const Sample* ir, size_t irLen)
{
// ==================================================================================
-// Copyright (c) 2012 HiFi-LoFi
+// Copyright (c) 2017 HiFi-LoFi
//
-// This 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 3 of the License, or
-// (at your option) any later version.
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is furnished
+// to do so, subject to the following conditions:
//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
//
-// You should have received a copy of the GNU General Public License
-// along with this program. If not, see <http://www.gnu.org/licenses/>.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// ==================================================================================
#ifndef _FFTCONVOLVER_FFTCONVOLVER_H
*/
void reset();
+ /**
+ * @brief Clears audio history
+ */
+
+ void clearState();
+
private:
size_t _blockSize;
size_t _segSize;
// ==================================================================================
-// Copyright (c) 2012 HiFi-LoFi
+// Copyright (c) 2017 HiFi-LoFi
//
-// This 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 3 of the License, or
-// (at your option) any later version.
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is furnished
+// to do so, subject to the following conditions:
//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
//
-// You should have received a copy of the GNU General Public License
-// along with this program. If not, see <http://www.gnu.org/licenses/>.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// ==================================================================================
#include "Utilities.h"
// ==================================================================================
-// Copyright (c) 2012 HiFi-LoFi
+// Copyright (c) 2017 HiFi-LoFi
//
-// This 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 3 of the License, or
-// (at your option) any later version.
+// Permission is hereby granted, free of charge, to any person obtaining a copy
+// of this software and associated documentation files (the "Software"), to deal
+// in the Software without restriction, including without limitation the rights
+// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the Software is furnished
+// to do so, subject to the following conditions:
//
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
+// The above copyright notice and this permission notice shall be included in
+// all copies or substantial portions of the Software.
//
-// You should have received a copy of the GNU General Public License
-// along with this program. If not, see <http://www.gnu.org/licenses/>.
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
// ==================================================================================
#ifndef _FFTCONVOLVER_UTILITIES_H
_size = size;
}
}
- if (_data)
- setZero();
+ setZero();
}
size_t size() const
#include "convolver.h"
+#include "ConvolverThreadPool.h"
#include "FFTConvolver.h"
#include "Utilities.h"
#include <pthread.h>
#include <sndfile.h>
+#include <unistd.h>
extern "C" void _warn(const char *filename, const int linenumber, const char *format, ...);
extern "C" void _debug(const char *filename, const int linenumber, int level, const char *format,
#define warn(...) _warn(__FILE__, __LINE__, __VA_ARGS__)
#define debug(...) _debug(__FILE__, __LINE__, __VA_ARGS__)
-fftconvolver::FFTConvolver convolver_l;
-fftconvolver::FFTConvolver convolver_r;
+// Create and initialize the thread pool
+ConvolverThreadPool pool;
-// always lock use this when accessing the playing conn value
-pthread_mutex_t convolver_lock = PTHREAD_MUTEX_INITIALIZER;
+void convolver_pool_init(size_t numThreads, size_t numConvolvers) {
+ if (!pool.init(numThreads, numConvolvers)) {
+ debug(1, "failed to initialize thread pool!");
+ } else {
+ debug(1, "thread pool initialized with %u threads and %u convolvers.", numThreads,
+ numConvolvers);
+ }
+}
-int convolver_init(const char *filename, int max_length) {
+void convolver_pool_closedown() {
+ pool.shutdown(); // Just shutdown, don't delete
+ debug(3, "thread pool shut down");
+}
+
+int convolver_init(const char *filename, unsigned char channel_count,
+ double max_length_in_seconds, size_t block_size) {
+ debug(3, "convolver_init");
int success = 0;
- SF_INFO info;
+ SF_INFO info = {}; // Zero everything, including format
if (filename) {
SNDFILE *file = sf_open(filename, SFM_READ, &info);
if (file) {
-
- if (info.samplerate == 44100) {
- if ((info.channels == 1) || (info.channels == 2)) {
- const size_t size = info.frames > max_length ? max_length : info.frames;
- float buffer[size * info.channels];
-
+ size_t max_length = (size_t)(max_length_in_seconds * info.samplerate);
+ const size_t size =
+ (unsigned int)info.frames > max_length ? max_length : (unsigned int)info.frames;
+ float *buffer = (float*)malloc(sizeof(float) * size * info.channels);
+ if (buffer != NULL) {
+ // float buffer[size * info.channels];
+ float *abuffer = (float*)malloc(sizeof(float) * size);
+ if (abuffer != NULL) {
size_t l = sf_readf_float(file, buffer, size);
if (l != 0) {
- pthread_mutex_lock(&convolver_lock);
- convolver_l.reset(); // it is possible that init could be called more than once
- convolver_r.reset(); // so it could be necessary to remove all previous settings
-
+ unsigned int cc;
if (info.channels == 1) {
- convolver_l.init(352, buffer, size);
- convolver_r.init(352, buffer, size);
- } else {
- // deinterleave
- float buffer_l[size];
- float buffer_r[size];
-
- unsigned int i;
- for (i = 0; i < size; ++i) {
- buffer_l[i] = buffer[2 * i + 0];
- buffer_r[i] = buffer[2 * i + 1];
+ for (cc = 0; cc < channel_count; cc++) {
+ if (!pool.initConvolver(cc, block_size, buffer, size)) {
+ debug(1, "new convolver failed to initialize convolver %u ", cc);
+ }
+ }
+ } else if (info.channels == channel_count) {
+ // we have to deinterleave the ir file channels for each convolver
+ // float abuffer[size];
+ for (cc = 0; cc < channel_count; cc++) {
+ unsigned int i;
+ for (i = 0; i < size; ++i) {
+ abuffer[i] = buffer[channel_count * i + cc];
+ }
+ if (!pool.initConvolver(cc, block_size, abuffer, size)) {
+ debug(1, "new convolver failed to initialize convolver %u ", cc);
+ }
}
-
- convolver_l.init(352, buffer_l, size);
- convolver_r.init(352, buffer_r, size);
}
- pthread_mutex_unlock(&convolver_lock);
success = 1;
}
- debug(1,
+ debug(2,
"convolution impulse response filter initialized from \"%s\" with %d channel%s and "
"%d samples",
filename, info.channels, info.channels == 1 ? "" : "s", size);
+ sf_close(file);
+ free((void*)abuffer);
} else {
- warn("Convolution impulse response filter file \"%s\" contains %d channels. Only 1 or 2 "
- "is supported.",
- filename, info.channels);
+ debug(1, "failed to init convolvers because insufficient memory was available");
}
+ free((void*)buffer);
} else {
- warn("Convolution impulse response filter file \"%s\" sample rate is %d Hz. Only 44100 Hz "
- "is supported.",
- filename, info.samplerate);
+ warn("failed to init convolvers because insufficient memory was available");
}
+ } else {
+ warn("Convolution impulse response filter file \"%s\" can not be opened. Please check that "
+ "it exists, is a valid sound file and has appropriate access permissions.",
+ filename);
+ }
+ }
+ return success;
+}
+
+void convolver_process(unsigned int channel, float *data, int length) {
+ pool.processAsync(channel, data, data, length);
+}
+
+void convolver_wait_for_all() { pool.waitForAll(); }
+
+void convolver_clear_state() {
+ pool.clearAllStates();
+}
+
+const unsigned int max_channels = 8;
+fftconvolver::FFTConvolver convolvers[max_channels];
+
+// fftconvolver::FFTConvolver convolver_l;
+// fftconvolver::FFTConvolver convolver_r;
+
+// always lock use this when accessing the playing conn value
+/*
+
+pthread_mutex_t convolver_lock = PTHREAD_MUTEX_INITIALIZER;
+
+int convolver_init(const char *filename, unsigned char channel_count, double max_length_in_seconds,
+ size_t block_size) {
+ debug(1, "convolver_init");
+ int success = 0;
+ SF_INFO info;
+ if (filename) {
+ SNDFILE *file = sf_open(filename, SFM_READ, &info);
+ if (file) {
+ size_t max_length = (size_t)(max_length_in_seconds * info.samplerate);
+ const size_t size =
+ (unsigned int)info.frames > max_length ? max_length : (unsigned int)info.frames;
+ float buffer[size * info.channels];
+
+ size_t l = sf_readf_float(file, buffer, size);
+ if (l != 0) {
+ pthread_mutex_lock(&convolver_lock);
+
+ unsigned int cc;
+ for (cc = 0; cc < channel_count; cc++) {
+ convolvers[cc].reset();
+ }
+
+ if (info.channels == 1) {
+ for (cc = 0; cc < channel_count; cc++) {
+ convolvers[cc].init(block_size, buffer, size);
+ }
+ } else if (info.channels == channel_count) {
+ // we have to deinterleave the ir file channels for each convolver
+ for (cc = 0; cc < channel_count; cc++) {
+ float abuffer[size];
+ unsigned int i;
+ for (i = 0; i < size; ++i) {
+ abuffer[i] = buffer[channel_count * i + cc];
+ }
+ convolvers[cc].init(block_size, abuffer, size);
+ }
+ }
+ pthread_mutex_unlock(&convolver_lock);
+ success = 1;
+ }
+ debug(2,
+ "convolution impulse response filter initialized from \"%s\" with %d channel%s and "
+ "%d samples",
+ filename, info.channels, info.channels == 1 ? "" : "s", size);
sf_close(file);
} else {
warn("Convolution impulse response filter file \"%s\" can not be opened. Please check that "
}
return success;
}
+void convolver_reset() {
+ debug(1, "convolver_reset");
+ pthread_mutex_lock(&convolver_lock);
+ unsigned int cc;
+ for (cc = 0; cc < max_channels; cc++) {
+ convolvers[cc].reset();
+ }
+ // convolver_l.reset(); // it is possible that init could be called more than once
+ // convolver_r.reset(); // so it could be necessary to remove all previous settings
+ pthread_mutex_unlock(&convolver_lock);
+}
+
+void convolver_clear_state() {
+ debug(1, "convolver_clear_state");
+ pthread_mutex_lock(&convolver_lock);
+ unsigned int cc;
+ for (cc = 0; cc < max_channels; cc++) {
+ convolvers[cc].clearState();
+ }
+ // convolver_l.reset(); // it is possible that init could be called more than once
+ // convolver_r.reset(); // so it could be necessary to remove all previous settings
+ pthread_mutex_unlock(&convolver_lock);
+}
+
+void convolver_process(unsigned int channel, float *data, int length) {
+ pthread_mutex_lock(&convolver_lock);
+ convolvers[channel].process(data, data, length);
+ pthread_mutex_unlock(&convolver_lock);
+ usleep(100);
+}
void convolver_process_l(float *data, int length) {
pthread_mutex_lock(&convolver_lock);
convolver_r.process(data, data, length);
pthread_mutex_unlock(&convolver_lock);
}
+*/
\ No newline at end of file
#ifdef __cplusplus
extern "C" {
#endif
+
+ #include <stddef.h>
-int convolver_init(const char* file, int max_length);
-void convolver_process_l(float* data, int length);
-void convolver_process_r(float* data, int length);
+// int convolver_init(const char* file, unsigned char channel_count, double max_length_in_seconds, size_t block_size);
+void convolver_reset();
+//void convolver_clear_state();
+// void convolver_process(unsigned int channel, float *data, int length);
+// void convolver_process_l(float* data, int length);
+// void convolver_process_r(float* data, int length);
+
+void convolver_pool_init(size_t numThreads, size_t numConvolvers);
+void convolver_pool_closedown();
+int convolver_init(const char* file, unsigned char channel_count, double max_length_in_seconds, size_t block_size);
+void convolver_process(unsigned int channel, float *data, int length);
+void convolver_clear_state();
+void convolver_wait_for_all();
#ifdef __cplusplus
}
endif
if USE_CONVOLUTION
-shairport_sync_SOURCES += FFTConvolver/AudioFFT.cpp FFTConvolver/FFTConvolver.cpp FFTConvolver/Utilities.cpp FFTConvolver/convolver.cpp
+shairport_sync_SOURCES += FFTConvolver/AudioFFT.cpp FFTConvolver/FFTConvolver.cpp FFTConvolver/Utilities.cpp FFTConvolver/convolver.cpp FFTConvolver/ConvolverThreadPool.cpp
AM_CXXFLAGS += -std=c++11
endif
#include <netinet/in.h>
#endif
+#ifdef CONFIG_CONVOLUTION
+#include <ctype.h>
+#include <sndfile.h>
+#endif
+
#ifdef CONFIG_OPENSSL
#include <openssl/aes.h> // needed for older AES stuff
#include <openssl/bio.h> // needed for BIO_new_mem_buf
}
#endif
-
-int config_lookup_non_empty_string(const config_t * the_config, const char * path, const char ** value) {
- int response = config_lookup_string(the_config, path, value);
+int config_lookup_non_empty_string(const config_t *cfg, const char *path, const char **value) {
+ int response = config_lookup_string(cfg, path, value);
if (response == CONFIG_TRUE) {
if ((value != NULL) && ((*value == NULL) || (*value[0] == 0))) {
warn("The \"%s\" parameter is an empty string and has been ignored.", path);
#endif
return ret;
}
+
+#ifdef CONFIG_CONVOLUTION
+/* Parse comma-separated filenames with optional quotes
+ * Returns array of ir_file_info_t structs (caller must free both array and filenames)
+ * count is set to number of filenames found
+ * Returns NULL on error
+ */
+ir_file_info_t *parse_ir_filenames(const char *input, unsigned int *file_count) {
+ if (!input || !file_count)
+ return NULL;
+
+ *file_count = 0;
+ unsigned int capacity = 10;
+ ir_file_info_t *files = malloc(capacity * sizeof(ir_file_info_t));
+ if (!files)
+ return NULL;
+
+ const char *p = input;
+
+ while (*p) {
+ /* Skip whitespace before filename */
+ while (isspace((unsigned char)*p))
+ p++;
+ if (!*p)
+ break;
+
+ /* Check if we need to resize array */
+ if (*file_count >= capacity) {
+ capacity *= 2;
+ ir_file_info_t *temp = realloc(files, capacity * sizeof(ir_file_info_t));
+ if (!temp) {
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+ files = temp;
+ }
+
+ /* Parse one filename */
+ char quote_char = 0;
+ char *buffer = NULL;
+ size_t buf_len = 0;
+ size_t buf_cap = 64;
+
+ if (*p == '"' || *p == '\'') {
+ /* Quoted filename */
+ quote_char = *p;
+ p++;
+
+ buffer = malloc(buf_cap);
+ if (!buffer) {
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+
+ /* Parse quoted string with escape handling */
+ while (*p && *p != quote_char) {
+ if (*p == '\\' && *(p + 1)) {
+ /* Escape sequence */
+ p++;
+ if (buf_len >= buf_cap - 1) {
+ buf_cap *= 2;
+ char *temp = realloc(buffer, buf_cap);
+ if (!temp) {
+ free(buffer);
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+ buffer = temp;
+ }
+ buffer[buf_len++] = *p++;
+ } else {
+ if (buf_len >= buf_cap - 1) {
+ buf_cap *= 2;
+ char *temp = realloc(buffer, buf_cap);
+ if (!temp) {
+ free(buffer);
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+ buffer = temp;
+ }
+ buffer[buf_len++] = *p++;
+ }
+ }
+ buffer[buf_len] = '\0';
+ if (*p == quote_char)
+ p++; /* Skip closing quote */
+
+ files[*file_count].samplerate = 0;
+ // files[*file_count].evaluation = ev_unchecked;
+ files[*file_count].filename = buffer;
+ (*file_count)++;
+ } else {
+ /* Unquoted filename - read until comma or end, handle escapes */
+ buffer = malloc(buf_cap);
+ if (!buffer) {
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+
+ while (*p && *p != ',') {
+ if (*p == '\\' && *(p + 1)) {
+ /* Escape sequence */
+ p++;
+ if (buf_len >= buf_cap - 1) {
+ buf_cap *= 2;
+ char *temp = realloc(buffer, buf_cap);
+ if (!temp) {
+ free(buffer);
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+ buffer = temp;
+ }
+ buffer[buf_len++] = *p++;
+ } else {
+ if (buf_len >= buf_cap - 1) {
+ buf_cap *= 2;
+ char *temp = realloc(buffer, buf_cap);
+ if (!temp) {
+ free(buffer);
+ for (unsigned int i = 0; i < *file_count; i++)
+ free(files[i].filename);
+ free(files);
+ return NULL;
+ }
+ buffer = temp;
+ }
+ buffer[buf_len++] = *p++;
+ }
+ }
+
+ /* Trim trailing whitespace */
+ while (buf_len > 0 && isspace((unsigned char)buffer[buf_len - 1])) {
+ buf_len--;
+ }
+ buffer[buf_len] = '\0';
+
+ files[*file_count].samplerate = 0;
+ files[*file_count].channels = 0;
+ // files[*file_count].evaluation = ev_unchecked;
+ files[*file_count].filename = buffer;
+ (*file_count)++;
+ }
+
+ /* Skip comma and whitespace */
+ while (isspace((unsigned char)*p))
+ p++;
+ if (*p == ',') {
+ p++;
+ while (isspace((unsigned char)*p))
+ p++;
+ }
+ }
+
+ return files;
+}
+
+/* Do a quick sanity check on the files -- see if they can be opened as sound files */
+void sanity_check_ir_files(const int option_print_level, ir_file_info_t *files,
+ unsigned int count) {
+ if (files != NULL) {
+ debug(option_print_level, "convolution impulse response files: %d found.", count);
+ for (unsigned int i = 0; i < count; i++) {
+
+ SF_INFO sfinfo = {};
+ // sfinfo.format = 0;
+
+ SNDFILE *file = sf_open(files[i].filename, SFM_READ, &sfinfo);
+ if (file) {
+ // files[i].evaluation = ev_okay;
+ files[i].samplerate = sfinfo.samplerate;
+ files[i].channels = sfinfo.channels;
+ debug(
+ option_print_level,
+ "convolution impulse response file \"%s\": %" PRId64 " frames (%.1f seconds), %d channel%s at %d frames per second.",
+ files[i].filename,
+ sfinfo.frames,
+ (float) sfinfo.frames / sfinfo.samplerate,
+ sfinfo.channels,
+ sfinfo.channels == 1 ? "" : "s",
+ sfinfo.samplerate);
+ sf_close(file);
+ } else {
+ // files[i].evaluation = ev_invalid;
+ debug(option_print_level, "convolution impulse response file \"%s\" %s", files[i].filename,
+ sf_strerror(NULL));
+ warn("Error accessing the convolution impulse response file \"%s\". %s", files[i].filename,
+ sf_strerror(NULL));
+ }
+ }
+ } else {
+ debug(option_print_level, "no convolution impulse response files found.");
+ }
+}
+
+/* Free the array returned by parse_filenames */
+void free_ir_filenames(ir_file_info_t *files, unsigned int file_count) {
+ if (!files)
+ return;
+ for (unsigned int i = 0; i < file_count; i++) {
+ free(files[i].filename);
+ }
+ free(files);
+}
+#endif
#define SAFAMILY sa_family
#endif
+#if defined(CONFIG_CONVOLUTION)
+typedef enum { ev_unchecked, ev_okay, ev_invalid } ir_file_evaluation;
+
+typedef struct {
+ unsigned int samplerate; // initialized to 0, will be filter frame rate
+ unsigned int channels;
+ char *filename; // the parsed filename
+} ir_file_info_t;
+#endif
+
#if defined(CONFIG_DBUS_INTERFACE) || defined(CONFIG_MPRIS_INTERFACE)
#include <glib.h>
typedef enum {
DBT_default = 0,
- DBT_system, // use the system bus
- DBT_session, // use the session bus
+ DBT_system, // use the system bus
+ DBT_session, // use the session bus
} dbus_message_bus_t;
#endif
uint32_t channel_set;
#ifdef CONFIG_CONVOLUTION
- int convolution;
- int convolver_valid;
- char *convolution_ir_file;
+ int convolution_enabled;
+ unsigned int convolution_rate; // 0 means the convolver has never been initialised, so ignore
+ // convolver_valid.
+ // but if this is the same as the current rate and convolver_valid is false, it means that an
+ // attempt to initialise the convolver has failed.
+ size_t convolution_block_size;
+ unsigned int convolution_ir_file_count;
+ ir_file_info_t *convolution_ir_files; // NULL or an array of information about all the impulse
+ // response files loaded
+ int convolution_ir_files_updated; // set to true if the convolution_ir_files are changed. Cleared
+ // when the convolver has been initialised
+ int convolver_valid; // set to true if the convolver can be initialised
+ unsigned int convolution_threads; // number of threads in the convolver thread pool
float convolution_gain;
- int convolution_max_length;
+ double convolution_max_length_in_seconds;
#endif
- int loudness;
+ int loudness_enabled;
float loudness_reference_volume_db;
int alsa_use_hardware_mute;
double alsa_maximum_stall_time;
extern uint64_t minimum_dac_queue_size;
-int config_lookup_non_empty_string(const config_t * the_config, const char * path, const char ** value);
+int config_lookup_non_empty_string(const config_t *cfg, const char *path, const char **value);
int config_set_lookup_bool(config_t *cfg, char *where, int *dst);
int check_string_or_list_setting(config_setting_t *setting, const char *item);
int check_int_or_list_setting(config_setting_t *setting, const int item);
char *bnprintf(char *buffer, ssize_t max_bytes, const char *format, ...);
+#ifdef CONFIG_CONVOLUTION
+
+/* Parse comma-separated filenames with optional quotes from the input string
+ * Returns array of ir_file_info_t structs (caller must free both array and filenames)
+ * count is set to number of filenames found
+ * Returns NULL on error
+ */
+ir_file_info_t *parse_ir_filenames(const char *input, unsigned int *file_count);
+// Access: files[i].filename, files[i].rate, files[i].evaluation
+
+/* Do a quick sanity check on the files -- see if they can be opened as sound files */
+void sanity_check_ir_files(const int option_print_level, ir_file_info_t *files, unsigned int count);
+
+/* Free the array returned by parse_filenames */
+void free_ir_filenames(ir_file_info_t *files, unsigned int file_count);
+
+#endif
+
#ifdef CONFIG_USE_GIT_VERSION_STRING
extern char git_version_string[];
#endif
else
AC_CHECK_LIB([sndfile], [sf_open], , AC_MSG_ERROR(Convolution support requires the sndfile library -- libsndfile1-dev suggested!))
fi
+
fi
AM_CONDITIONAL([USE_CONVOLUTION], [test "x$with_convolution" = "xyes"])
}
#ifdef CONFIG_CONVOLUTION
-gboolean notify_convolution_callback(ShairportSync *skeleton,
- __attribute__((unused)) gpointer user_data) {
+gboolean notify_convolution_enabled_callback(ShairportSync *skeleton,
+ __attribute__((unused)) gpointer user_data) {
// debug(1, "\"notify_convolution_callback\" called.");
- if (shairport_sync_get_convolution(skeleton)) {
- debug(1, ">> activate convolution filter");
- config.convolution = 1;
- config.convolver_valid =
- convolver_init(config.convolution_ir_file, config.convolution_max_length);
+ if (shairport_sync_get_convolution_enabled(skeleton)) {
+ debug(1, ">> activate convolution impulse response filter");
+ config.convolution_enabled = 1;
} else {
debug(1, ">> deactivate convolution impulse response filter");
- config.convolution = 0;
+ config.convolution_enabled = 0;
+ convolver_clear_state();
}
return TRUE;
}
#else
-gboolean notify_convolution_callback(__attribute__((unused)) ShairportSync *skeleton,
- __attribute__((unused)) gpointer user_data) {
+gboolean notify_convolution_enabled_callback(__attribute__((unused)) ShairportSync *skeleton,
+ __attribute__((unused)) gpointer user_data) {
+ warn(">> Convolution support is not built in to this build of Shairport Sync.");
+ return TRUE;
+}
+#endif
+
+#ifdef CONFIG_CONVOLUTION
+gboolean notify_convolution_maximum_length_in_seconds_callback(ShairportSync *skeleton,
+ __attribute__((unused))
+ gpointer user_data) {
+
+ gdouble th = shairport_sync_get_convolution_maximum_length_in_seconds(skeleton);
+ if ((th >= 0.0) && (th <= 15.0)) {
+ debug(1, ">> set convolution maximum length in seconds to %f.", th);
+ config.convolution_max_length_in_seconds = th;
+ } else {
+ debug(1, ">> invalid convolution gain: %f. Ignored.", th);
+ shairport_sync_set_convolution_maximum_length_in_seconds(
+ skeleton, config.convolution_max_length_in_seconds);
+ }
+ return TRUE;
+}
+#else
+gboolean notify_convolution_maximum_length_in_seconds_callback(__attribute__((unused))
+ ShairportSync *skeleton,
+ __attribute__((unused))
+ gpointer user_data) {
warn(">> Convolution support is not built in to this build of Shairport Sync.");
return TRUE;
}
__attribute__((unused)) gpointer user_data) {
gdouble th = shairport_sync_get_convolution_gain(skeleton);
- if ((th <= 0.0) && (th >= -100.0)) {
+ if ((th <= 18.0) && (th >= -60.0)) {
debug(1, ">> set convolution gain to %f.", th);
config.convolution_gain = th;
} else {
}
#endif
#ifdef CONFIG_CONVOLUTION
-gboolean notify_convolution_impulse_response_file_callback(ShairportSync *skeleton,
- __attribute__((unused))
- gpointer user_data) {
- char *th = (char *)shairport_sync_get_convolution_impulse_response_file(skeleton);
- if ((config.convolution_ir_file == 0) || (config.convolver_valid == 0) ||
- (strcmp(config.convolution_ir_file, th) != 0)) {
- if (config.convolution_ir_file != NULL) {
- debug(1, ">> freeing current configuration impulse response filter file \"%s\".",
- config.convolution_ir_file);
- free(config.convolution_ir_file);
- }
- config.convolution_ir_file = strdup(th);
- debug(1, ">> set configuration impulse response filter file to \"%s\".",
- config.convolution_ir_file);
- config.convolver_valid =
- convolver_init(config.convolution_ir_file, config.convolution_max_length);
+gboolean notify_convolution_impulse_response_files_callback(ShairportSync *skeleton,
+ __attribute__((unused))
+ gpointer user_data) {
+ char *th = (char *)shairport_sync_get_convolution_impulse_response_files(skeleton);
+ if (config.convolution_ir_files != NULL) {
+ debug(1, ">> freeing current configuration impulse response filter files.");
+ free_ir_filenames(config.convolution_ir_files, config.convolution_ir_file_count);
+ config.convolution_ir_files = NULL;
+ config.convolution_ir_file_count = 0;
}
+ config.convolution_ir_files = parse_ir_filenames(th, &config.convolution_ir_file_count);
+ sanity_check_ir_files(1, config.convolution_ir_files, config.convolution_ir_file_count);
+ debug(1, ">> setting %d configuration impulse response filter%s",
+ config.convolution_ir_file_count, config.convolution_ir_file_count == 1 ? "" : "s");
+ config.convolution_ir_files_updated = 1;
return TRUE;
}
#else
-gboolean notify_convolution_impulse_response_file_callback(__attribute__((unused))
- ShairportSync *skeleton,
- __attribute__((unused))
- gpointer user_data) {
+gboolean notify_convolution_impulse_response_files_callback(__attribute__((unused))
+ ShairportSync *skeleton,
+ __attribute__((unused))
+ gpointer user_data) {
__attribute__((unused)) char *th =
(char *)shairport_sync_get_convolution_impulse_response_file(skeleton);
return TRUE;
}
#endif
-gboolean notify_loudness_callback(ShairportSync *skeleton,
- __attribute__((unused)) gpointer user_data) {
+gboolean notify_loudness_enabled_callback(ShairportSync *skeleton,
+ __attribute__((unused)) gpointer user_data) {
// debug(1, "\"notify_loudness_callback\" called.");
- if (shairport_sync_get_loudness(skeleton)) {
+ if (shairport_sync_get_loudness_enabled(skeleton)) {
debug(1, ">> activate loudness filter");
- config.loudness = 1;
+ config.loudness_enabled = 1;
} else {
debug(1, ">> deactivate loudness filter");
- config.loudness = 0;
+ config.loudness_enabled = 0;
}
return TRUE;
}
static void on_dbus_name_acquired(GDBusConnection *connection, const gchar *name,
__attribute__((unused)) gpointer user_data) {
-
+ const char *str = NULL;
debug(2, "Shairport Sync native D-Bus interface \"%s\" acquired on the %s bus.", name,
(dbus_bus_type == G_BUS_TYPE_SESSION) ? "session" : "system");
G_CALLBACK(notify_volume_control_profile_callback), NULL);
g_signal_connect(shairportSyncSkeleton, "notify::disable-standby",
G_CALLBACK(notify_disable_standby_callback), NULL);
- g_signal_connect(shairportSyncSkeleton, "notify::convolution",
- G_CALLBACK(notify_convolution_callback), NULL);
+ g_signal_connect(shairportSyncSkeleton, "notify::convolution-enabled",
+ G_CALLBACK(notify_convolution_enabled_callback), NULL);
g_signal_connect(shairportSyncSkeleton, "notify::convolution-gain",
G_CALLBACK(notify_convolution_gain_callback), NULL);
- g_signal_connect(shairportSyncSkeleton, "notify::convolution-impulse-response-file",
- G_CALLBACK(notify_convolution_impulse_response_file_callback), NULL);
- g_signal_connect(shairportSyncSkeleton, "notify::loudness", G_CALLBACK(notify_loudness_callback),
- NULL);
+ g_signal_connect(shairportSyncSkeleton, "notify::convolution-maximum-length-in-seconds",
+ G_CALLBACK(notify_convolution_maximum_length_in_seconds_callback), NULL);
+ g_signal_connect(shairportSyncSkeleton, "notify::convolution-impulse-response-files",
+ G_CALLBACK(notify_convolution_impulse_response_files_callback), NULL);
+ g_signal_connect(shairportSyncSkeleton, "notify::loudness-enabled",
+ G_CALLBACK(notify_loudness_enabled_callback), NULL);
g_signal_connect(shairportSyncSkeleton, "notify::loudness-threshold",
G_CALLBACK(notify_loudness_threshold_callback), NULL);
g_signal_connect(shairportSyncSkeleton, "notify::drift-tolerance",
shairport_sync_set_disable_standby(SHAIRPORT_SYNC(shairportSyncSkeleton), TRUE);
}
- if (config.loudness == 0) {
- debug(1, ">> loudness set to \"off\"");
- shairport_sync_set_loudness(SHAIRPORT_SYNC(shairportSyncSkeleton), FALSE);
+ if (config.loudness_enabled == 0) {
+ debug(1, ">> loudness_enabled is false");
+ shairport_sync_set_loudness_enabled(SHAIRPORT_SYNC(shairportSyncSkeleton), FALSE);
} else {
- debug(1, ">> loudness set to \"on\"");
- shairport_sync_set_loudness(SHAIRPORT_SYNC(shairportSyncSkeleton), TRUE);
+ debug(1, ">> loudness_enabled is true");
+ shairport_sync_set_loudness_enabled(SHAIRPORT_SYNC(shairportSyncSkeleton), TRUE);
}
#ifdef CONFIG_CONVOLUTION
- if (config.convolution == 0) {
- debug(1, ">> convolution set to \"off\"");
- shairport_sync_set_convolution(SHAIRPORT_SYNC(shairportSyncSkeleton), FALSE);
+ if (config.convolution_enabled == 0) {
+ debug(1, ">> convolution_enabled is false");
+ shairport_sync_set_convolution_enabled(SHAIRPORT_SYNC(shairportSyncSkeleton), FALSE);
+ } else {
+ debug(1, ">> convolution_enabled is true");
+ shairport_sync_set_convolution_enabled(SHAIRPORT_SYNC(shairportSyncSkeleton), TRUE);
+ }
+
+ if ((config.cfg != NULL) &&
+ (config_lookup_non_empty_string(config.cfg, "dsp.convolution_ir_files", &str))) {
+ shairport_sync_set_convolution_impulse_response_files(SHAIRPORT_SYNC(shairportSyncSkeleton),
+ str);
} else {
- debug(1, ">> convolution set to \"on\"");
- shairport_sync_set_convolution(SHAIRPORT_SYNC(shairportSyncSkeleton), TRUE);
+ shairport_sync_set_convolution_impulse_response_files(SHAIRPORT_SYNC(shairportSyncSkeleton),
+ NULL);
}
- if (config.convolution_ir_file)
- shairport_sync_set_convolution_impulse_response_file(SHAIRPORT_SYNC(shairportSyncSkeleton),
- config.convolution_ir_file);
-// else
-// shairport_sync_set_convolution_impulse_response_file(SHAIRPORT_SYNC(shairportSyncSkeleton),
-// NULL);
+ shairport_sync_set_convolution_maximum_length_in_seconds(
+ SHAIRPORT_SYNC(shairportSyncSkeleton), config.convolution_max_length_in_seconds);
#endif
shairport_sync_set_service_name(SHAIRPORT_SYNC(shairportSyncSkeleton), config.service_name);
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:DriftTolerance variant:double:0.001
# Is Loudness Enabled:
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:LoudnessThreshold
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:LoudnessEnabled
# Enable Loudness Filter
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:Loudness variant:boolean:true
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:LoudnessEnabled variant:boolean:true
# Get Loudness Threshold
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:LoudnessThreshold
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:LoudnessThreshold variant:double:-15.0
# Is Convolution enabled:
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:Convolution
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:ConvolutionEnabled
# Enable Convolution
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:Convolution variant:boolean:true
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:ConvolutionEnabled variant:boolean:true
# Get Convolution Gain:
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:ConvolutionGain
-# Set Convolution Gain -- the gain applied before convolution is applied -- to -10.0 dB
+# Set Convolution Gain -- the gain applied after convolution is applied -- to -10.0 dB
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:ConvolutionGain variant:double:-10
-# Get Convolution Impulse Response File:
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:ConvolutionImpulseResponseFile
+# Get Convolution Impulse Response Files:
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:ConvolutionImpulseResponseFiles
+
+# Set Convolution Impulse Response Files:
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:ConvolutionImpulseResponseFiles variant:string:"'/home/pi/filters/Sennheiser HD 205 minimum phase 48000Hz.wav', '/home/pi/filters/Sennheiser HD 205 minimum phase 44100Hz.wav'"
+
+# Get Convolution Impulse Response File Maximum Length:
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:ConvolutionMaximumLengthInSeconds
-# Set Convolution Impulse Response File:
-dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:ConvolutionImpulseResponseFile variant:string:"/etc/shairport-sync/boom.wav"
+# Set Convolution Impulse Response File Maximum Length:
+dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Set string:org.gnome.ShairportSync string:ConvolutionMaximumLengthInSeconds variant:double:1
# Get the Protocol Shairport Sync was built for -- AirPlay or AirPlay 2:
dbus-send --print-reply --system --dest=org.gnome.ShairportSync /org/gnome/ShairportSync org.freedesktop.DBus.Properties.Get string:org.gnome.ShairportSync string:Protocol
#include "common.h"
#include <math.h>
-loudness_processor loudness_r;
-loudness_processor loudness_l;
+#define MAXCHANNELS 8
-void _loudness_set_volume(loudness_processor *p, float volume) {
+// loudness_processor_dynamic loudness_r;
+// loudness_processor_dynamic loudness_l;
+
+loudness_processor_dynamic loudness_filters[MAXCHANNELS];
+
+// if the rate or the loudness_reference_volume_db change, recalculate the
+// loudness volume parameters
+
+static int loudness_fix_volume_parameter = 0;
+static unsigned int loudness_rate_parameter = 0;
+static float loudness_volume_reference_parameter = 0.0;
+static loudness_processor_static lps = {0.0, 0.0, 0.0, 0.0, 0.0};
+
+void _loudness_set_volume(loudness_processor_static *p, float volume, unsigned int sample_rate) {
float gain = -(volume - config.loudness_reference_volume_db) * 0.5;
- if (gain < 0)
+ if (gain < 0) {
gain = 0;
+ }
+ debug(2, "Volume: %.1f dB - Loudness gain @10Hz: %.1f dB", volume, gain);
float Fc = 10.0;
float Q = 0.5;
// Formula from http://www.earlevel.com/main/2011/01/02/biquad-formulas/
- float Fs = 44100.0;
+ float Fs = sample_rate * 1.0;
float K = tan(M_PI * Fc / Fs);
float V = pow(10.0, gain / 20.0);
p->b2 = (1 - 1 / Q * K + K * K) * norm;
}
-float loudness_process(loudness_processor *p, float i0) {
- float o0 = p->a0 * i0 + p->a1 * p->i1 + p->a2 * p->i2 - p->b1 * p->o1 - p->b2 * p->o2;
+float loudness_process(loudness_processor_dynamic *p, float i0) {
+ float o0 = lps.a0 * i0 + lps.a1 * p->i1 + lps.a2 * p->i2 - lps.b1 * p->o1 - lps.b2 * p->o2;
p->o2 = p->o1;
p->o1 = o0;
return o0;
}
-void loudness_set_volume(float volume) {
- float gain = -(volume - config.loudness_reference_volume_db) * 0.5;
- if (gain < 0)
- gain = 0;
+void loudness_update(rtsp_conn_info *conn) {
+ // first, see if loudness can be enabled
+ int do_loudness = config.loudness_enabled;
+ if ((config.output->parameters != NULL) && (config.output->parameters()->volume_range != NULL)) {
+ do_loudness = 0; // if we are using external (usually hardware) volume controls.
+ }
+
+ if (do_loudness) {
+ // check the volume parameters
+ if ((conn->fix_volume != loudness_fix_volume_parameter) ||
+ (conn->input_rate != loudness_rate_parameter) ||
+ (config.loudness_reference_volume_db != loudness_volume_reference_parameter)) {
+ debug(1, "update loudness parameters");
+ float new_volume = 20 * log10((double)conn->fix_volume / 65536);
+ _loudness_set_volume(&lps, new_volume, conn->input_rate);
+ // _loudness_set_volume(&loudness_r, new_volume, conn->input_rate);
+ loudness_fix_volume_parameter = conn->fix_volume;
+ loudness_rate_parameter = conn->input_rate;
+ loudness_volume_reference_parameter = config.loudness_reference_volume_db;
+ }
+ }
+ conn->do_loudness = do_loudness;
+}
- debug(3, "Volume: %.1f dB - Loudness gain @10Hz: %.1f dB", volume, gain);
- _loudness_set_volume(&loudness_l, volume);
- _loudness_set_volume(&loudness_r, volume);
+void loudness_process_blocks(float *fbufs, unsigned int channel_length,
+ unsigned int number_of_channels, float gain) {
+ unsigned int channel_number, sample_index;
+ float *sample_pointer = fbufs;
+ for (channel_number = 0; channel_number < number_of_channels; channel_number++) {
+ for (sample_index = 0; sample_index < channel_length; sample_index++) {
+ *sample_pointer = loudness_process(&loudness_filters[channel_number], *sample_pointer * gain);
+ sample_pointer++;
+ }
+ }
}
+
+void loudness_reset() {
+ unsigned int i;
+ for (i = 0; i < MAXCHANNELS; i++) {
+ loudness_filters[i].i1 = 0.0;
+ loudness_filters[i].i2 = 0.0;
+ loudness_filters[i].o1 = 0.0;
+ loudness_filters[i].o2 = 0.0;
+ }
+}
\ No newline at end of file
#pragma once
+#include "player.h"
#include <stdio.h>
typedef struct {
float a0, a1, a2, b1, b2;
+} loudness_processor_static;
+
+typedef struct {
float i1, i2, o1, o2;
-} loudness_processor;
+} loudness_processor_dynamic;
-extern loudness_processor loudness_r;
-extern loudness_processor loudness_l;
+// extern loudness_processor_dynamic loudness_r;
+// extern loudness_processor_dynamic loudness_l;
void loudness_set_volume(float volume);
-float loudness_process(loudness_processor *p, float sample);
+float loudness_process(loudness_processor_dynamic *p, float sample);
+void loudness_update(rtsp_conn_info *conn);
+
+void loudness_reset();
+void loudness_process_blocks(float *fbufs, unsigned int channel_length,
+ unsigned int number_of_channels, float gain);
\ No newline at end of file
<property name='Active' type='b' access='read'/>
<property name="DisableStandby" type="b" access="readwrite" />
<property name="DisableStandbyMode" type="s" access="readwrite" />
- <property name="Loudness" type="b" access="readwrite" />
+ <property name="LoudnessEnabled" type="b" access="readwrite" />
<property name="LoudnessThreshold" type="d" access="readwrite" />
- <property name="Convolution" type="b" access="readwrite" />
+ <property name="ConvolutionEnabled" type="b" access="readwrite" />
<property name="ConvolutionGain" type="d" access="readwrite" />
- <property name="ConvolutionImpulseResponseFile" type="s" access="readwrite" />
+ <property name="ConvolutionImpulseResponseFiles" type="s" access="readwrite" />
+ <property name="ConvolutionMaximumLengthInSeconds" type="d" access="readwrite" />
<property name="DriftTolerance" type="d" access="readwrite" />
<property name='Volume' type='d' access='readwrite'/>
<method name="DropSession"/>
#if LIBAVUTIL_VERSION_MAJOR >= 57
- AVChannelLayout output_channel_layout;
+ AVChannelLayout output_channel_layout = {0};
av_opt_get_chlayout(swr, "out_chlayout", 0, &output_channel_layout);
conn->resampler_output_channels = output_channel_layout.nb_channels;
for (c = 0; c < 64; c++) {
int64_t hyper_sample = sample;
int result = 0;
- if ((config.loudness != 0) && (conn->input_rate == 44100)) {
+ if (conn->do_loudness != 0) {
hyper_sample <<=
32; // Do not apply volume as it has already been done with the Loudness DSP filter
} else {
// we need to set up the output device to correspond to
// the input format w.r.t. rate, depth and channels
// because we'll be sending silence before the first real frame.
+ debug(3, "reset loudness filters.");
+ loudness_reset();
#ifdef CONFIG_FFMPEG
// Set up the output chain, including the software resampler.
debug(2, "set up the output chain to %s for FFmpeg.",
uint64_t time_of_last_metadata_progress_update =
0; // the assignment is to stop a compiler warning...
#endif
+ double highest_convolver_output_db = 0.0;
uint64_t previous_frames_played = 0; // initialised to avoid a "possibly uninitialised" warning
uint64_t previous_raw_measurement_time =
0; // initialised to avoid a "possibly uninitialised" warning
malloc(sps_format_sample_size(
FORMAT_FROM_ENCODED_FORMAT(config.current_output_configuration)) *
CHANNELS_FROM_ENCODED_FORMAT(config.current_output_configuration) *
- ((inframe->length) * conn->output_sample_ratio + +INTERPOLATION_LIMIT));
+ ((inframe->length) * conn->output_sample_ratio + INTERPOLATION_LIMIT));
if (conn->outbuf == NULL)
die("Failed to allocate memory for an output buffer.");
// now, before outputting anything to the output device, check the stats
- uint32_t stats_logging_interval_in_frames = 8 * RATE_FROM_ENCODED_FORMAT(config.current_output_configuration);
- if ((stats_logging_interval_in_frames != 0) && (frames_since_last_stats_logged > stats_logging_interval_in_frames)) {
+ uint32_t stats_logging_interval_in_frames =
+ 8 * RATE_FROM_ENCODED_FORMAT(config.current_output_configuration);
+ if ((stats_logging_interval_in_frames != 0) &&
+ (frames_since_last_stats_logged > stats_logging_interval_in_frames)) {
// here, calculate the input and output frame rates, where possible, even if
// statistics have not been requested this is to calculate them in case they are
gap_to_fix = -inframe->timestamp_gap; // this is frames at the input rate
int64_t gap_to_fix_ns = (gap_to_fix * 1000000000) / conn->input_rate;
gap_to_fix = (gap_to_fix_ns *
- RATE_FROM_ENCODED_FORMAT(config.current_output_configuration)) /
- 1000000000; // this is frames at the output rate
- // debug(3, "due to timstamp gap of %d frames, skip %" PRId64 " output frames.", inframe->timestamp_gap, gap_to_fix);
+ RATE_FROM_ENCODED_FORMAT(config.current_output_configuration)) /
+ 1000000000; // this is frames at the output rate
+ // debug(3, "due to timstamp gap of %d frames, skip %" PRId64 " output
+ // frames.", inframe->timestamp_gap, gap_to_fix);
}
}
-
-
-
if (gap_to_fix > 0) {
// debug(1, "drop %u frames, timestamp: %u, skipping_frames_at_start_of_play is
sync_error_ns = 0; // don't invoke any sync checking
sync_error = 0;
}
- }
+ }
}
// debug(1, "frames_to_skip: %u.", frames_to_skip);
// don't do any sync error calculations if you're skipping frames
// Apply DSP here
- // check the state of loudness and convolution flags here and don't change them
- // for the frame
-
- int do_loudness = ((config.loudness != 0) && (conn->input_rate == 44100));
-
-#ifdef CONFIG_CONVOLUTION
- int do_convolution = 0;
- if ((config.convolution) && (config.convolver_valid) && (conn->input_rate == 44100))
- do_convolution = 1;
-
- // we will apply the convolution gain if convolution is enabled, even if there
- // is no valid convolution happening
-
- int convolution_is_enabled = 0;
- if (config.convolution)
- convolution_is_enabled = 1;
-#endif
+ loudness_update(conn);
- if (do_loudness
+ if (conn->do_loudness
#ifdef CONFIG_CONVOLUTION
- || convolution_is_enabled
+ || config.convolution_enabled
#endif
) {
- int32_t *tbuf32 = (int32_t *)conn->tbuf;
- float fbuf_l[inbuflength];
- float fbuf_r[inbuflength];
+
+ float(*fbufs)[1024] = malloc(conn->input_num_channels * sizeof(*fbufs));
+ // debug(1, "size of array allocated is %d bytes.", conn->input_num_channels *
+ // sizeof(*fbufs));
+ int32_t *tbuf32 = conn->tbuf;
// Deinterleave, and convert to float
- int i;
- for (i = 0; i < inbuflength; ++i) {
- fbuf_l[i] = tbuf32[2 * i];
- fbuf_r[i] = tbuf32[2 * i + 1];
+ unsigned int i, j;
+ for (i = 0; i < inframe->length; i++) {
+ for (j = 0; j < conn->input_num_channels; j++) {
+ fbufs[j][i] = tbuf32[conn->input_num_channels * i + j];
+ }
}
#ifdef CONFIG_CONVOLUTION
// Apply convolution
- if (do_convolution) {
- convolver_process_l(fbuf_l, inbuflength);
- convolver_process_r(fbuf_r, inbuflength);
- }
- if (convolution_is_enabled) {
+ // First, have we got the right convolution setup?
+
+ static int convolver_is_valid = 0;
+ static size_t current_convolver_block_size = 0;
+ static unsigned int current_convolver_rate = 0;
+ static unsigned int current_convolver_channels = 0;
+ static double current_convolver_maximum_length_in_seconds = 0;
+
+ if (config.convolution_enabled) {
+ if (
+ // if any of these are true, we need to create a new convolver
+ // (conn->convolver_is_valid == 0) ||
+ (current_convolver_block_size != inframe->length) ||
+ (current_convolver_rate != conn->input_rate) ||
+ !((current_convolver_channels == 1) ||
+ (current_convolver_channels == conn->input_num_channels)) ||
+ (current_convolver_maximum_length_in_seconds !=
+ config.convolution_max_length_in_seconds) ||
+ (config.convolution_ir_files_updated == 1)) {
+
+ // look for a convolution ir file with a matching rate and channel count
+
+ convolver_is_valid = 0; // declare any current convolver as invalid
+ current_convolver_block_size = inframe->length;
+ current_convolver_rate = conn->input_rate;
+ current_convolver_channels = conn->input_num_channels;
+ current_convolver_maximum_length_in_seconds =
+ config.convolution_max_length_in_seconds;
+ config.convolution_ir_files_updated = 0;
+ debug(2, "try to initialise a %u/%u convolver.", current_convolver_rate,
+ current_convolver_channels);
+ char *convolver_file_found = NULL;
+ unsigned int ir = 0;
+ while ((ir < config.convolution_ir_file_count) &&
+ (convolver_file_found == NULL)) {
+ if ((config.convolution_ir_files[ir].samplerate ==
+ current_convolver_rate) &&
+ (config.convolution_ir_files[ir].channels ==
+ current_convolver_channels)) {
+ convolver_file_found = config.convolution_ir_files[ir].filename;
+ } else {
+ ir++;
+ }
+ }
+
+ // if no luck, try for a single-channel IR file
+ if (convolver_file_found == NULL) {
+ current_convolver_channels = 1;
+ ir = 0;
+ while ((ir < config.convolution_ir_file_count) &&
+ (convolver_file_found == NULL)) {
+ if ((config.convolution_ir_files[ir].samplerate ==
+ current_convolver_rate) &&
+ (config.convolution_ir_files[ir].channels ==
+ current_convolver_channels)) {
+ convolver_file_found = config.convolution_ir_files[ir].filename;
+ } else {
+ ir++;
+ }
+ }
+ }
+ if (convolver_file_found != NULL) {
+ // we have an apparently suitable convolution ir file, so lets initialise
+ // a convolver
+ convolver_is_valid = convolver_init(
+ convolver_file_found, conn->input_num_channels,
+ config.convolution_max_length_in_seconds, inframe->length);
+ convolver_wait_for_all();
+ // if (convolver_is_valid)
+ // debug(1, "convolver_init for %u channels was successful.", conn->input_num_channels);
+ // convolver_is_valid = convolver_init(
+ // convolver_file_found, conn->input_num_channels,
+ // config.convolution_max_length_in_seconds, inframe->length);
+ }
+
+ if (convolver_is_valid == 0)
+ debug(1, "can not initialise a %u/%u convolver.", current_convolver_rate,
+ conn->input_num_channels);
+ else
+ debug(1, "convolver: \"%s\".", convolver_file_found);
+ }
+ if (convolver_is_valid != 0) {
+ for (j = 0; j < conn->input_num_channels; j++) {
+ // convolver_process(j, fbufs[j], inframe->length);
+ convolver_process(j, fbufs[j], inframe->length);
+ }
+ convolver_wait_for_all();
+ }
+
+ // apply convolution gain even if no convolution is done...
float gain = pow(10.0, config.convolution_gain / 20.0);
- for (i = 0; i < inbuflength; ++i) {
- fbuf_l[i] *= gain;
- fbuf_r[i] *= gain;
+ for (i = 0; i < inframe->length; ++i) {
+ for (j = 0; j < conn->input_num_channels; j++) {
+ float output_level_db = 0.0;
+ if (fbufs[j][i] < 0.0)
+ output_level_db = 20 * log10(fbufs[j][i] / (float)INT32_MIN * 1.0);
+ else
+ output_level_db = 20 * log10(fbufs[j][i] / (float)INT32_MAX);
+ if (output_level_db > highest_convolver_output_db) {
+ highest_convolver_output_db = output_level_db;
+ if ((highest_convolver_output_db + config.convolution_gain) > 0.0)
+ warn("clipping %.1f dB with convolution gain set to %.1f dB!", highest_convolver_output_db + config.convolution_gain, config.convolution_gain);
+ }
+ fbufs[j][i] *= gain;
+ }
}
}
+
#endif
+ if (conn->do_loudness) {
+ loudness_process_blocks((float *)fbufs, inframe->length,
+ conn->input_num_channels,
+ (float)conn->fix_volume / 65536);
+ }
- if (do_loudness) {
- // Apply volume and loudness
- // Volume must be applied here because the loudness filter will increase the
- // signal level and it would saturate the int32_t otherwise
- float gain = conn->fix_volume / 65536.0f;
- // float gain_db = 20 * log10(gain);
- // debug(1, "Applying soft volume dB: %f k: %f", gain_db, gain);
-
- for (i = 0; i < inbuflength; ++i) {
- fbuf_l[i] = loudness_process(&loudness_l, fbuf_l[i] * gain);
- fbuf_r[i] = loudness_process(&loudness_r, fbuf_r[i] * gain);
+ // Interleave and convert back to int32_t
+ for (i = 0; i < inframe->length; i++) {
+ for (j = 0; j < conn->input_num_channels; j++) {
+ tbuf32[conn->input_num_channels * i + j] = fbufs[j][i];
}
}
- // Interleave and convert back to int32_t
- for (i = 0; i < inbuflength; ++i) {
- tbuf32[2 * i] = fbuf_l[i];
- tbuf32[2 * i + 1] = fbuf_r[i];
+ if (fbufs != NULL) {
+ free(fbufs);
+ fbufs = NULL;
}
}
+ // }
#ifdef CONFIG_SOXR
double t = config.audio_backend_buffer_interpolation_threshold_in_seconds *
RATE_FROM_ENCODED_FORMAT(config.current_output_configuration);
software_attenuation, temp_fix_volume);
conn->fix_volume = temp_fix_volume;
-
- // if (config.loudness)
- loudness_set_volume(software_attenuation / 100);
}
if (conn != NULL)
debug(3, "Connection %d: AirPlay Volume set to %.3f, Output Level set to: %.2f dB.",
void player_flush(uint32_t timestamp, rtsp_conn_info *conn) {
debug(3, "player_flush");
do_flush(timestamp, conn);
+#ifdef CONFIG_CONVOLUTION
+ convolver_clear_state();
+#endif
#ifdef CONFIG_METADATA
// only send a flush metadata message if the first packet has been seen -- it's a bogus message
// otherwise
}
free(pt);
// reset_anchor_info(conn); // say the clock is no longer valid
+#ifdef CONFIG_CONVOLUTION
+ convolver_clear_state();
+#endif
response = 0; // deleted
} else {
debug(2, "Connection %d: no player thread.", conn->connection_number);
uint32_t flush_rtp_timestamp;
uint64_t time_of_last_audio_packet;
seq_t ab_read, ab_write;
+
+ int do_loudness; // if loudness is requested and there is no external mixer
#ifdef CONFIG_MBEDTLS
mbedtls_aes_context dctx;
#include <sys/types.h>
#include <time.h>
#include <unistd.h>
+
#ifdef CONFIG_AIRPLAY_2
// #include "plist_xml_strings.h"
#include "ptp-utilities.h"
#include <sodium.h>
#endif
+#ifdef CONFIG_CONVOLUTION
+#include "FFTConvolver/convolver.h"
+#endif
+
+
struct Nvll {
char *name;
double value;
// 0; // This may be set to 1 by a flush, so don't zero it during start.
packets_played_in_this_sequence = 0;
new_buffer_needed = 0;
+#ifdef CONFIG_CONVOLUTION
+ convolver_clear_state();
+#endif
}
if ((play_enabled == 0) && (conn->ap2_play_enabled != 0)) {
#endif
#endif
+#ifdef CONFIG_CONVOLUTION
+#include "FFTConvolver/convolver.h"
+#endif
+
+
#ifdef CONFIG_DBUS_INTERFACE
#include "dbus-service.h"
#endif
debug(2, "Connection %d: SETRATEANCHORI Pause playing.", conn->connection_number);
conn->ap2_play_enabled = 0;
activity_monitor_signify_activity(0);
+#ifdef CONFIG_CONVOLUTION
+ // convolver_clear_state();
+#endif
+
// reset_anchor_info(conn);
#ifdef CONFIG_METADATA
send_ssnc_metadata('paus', NULL, 0, 1); // pause -- contains cancellation points
// --with-convolution
dsp =
{
-
-// NOTE: convolution and loudness filters currently work only for 44100 input and are disabled for other frame rates.
-
//////////////////////////////////////////
-// This convolution filter can be used to apply almost any correction to the audio signal, like frequency and phase correction.
+// This convolution filter can be used with inpulse responses files
+// to apply almost any correction to the audio signal, like frequency and phase correction.
// For example you could measure (with a good microphone and a sweep-sine) the frequency response of your speakers + room,
// and apply a correction to get a flat response curve.
+// Long impulse response files require lots of floating-point processing!
//////////////////////////////////////////
//
-// convolution = "no"; // Set this to "yes" to activate the convolution filter.
-// convolution_ir_file = "impulse.wav"; // Impulse Response file to be convolved to the audio stream
-// convolution_gain = -4.0; // Static gain applied to prevent clipping during the convolution process
-// convolution_max_length = 44100; // Truncate the input file to this length in order to save CPU.
-
-
+// convolution_enabled = "no"; // Set this to "yes" to activate the convolution filter.
+// convolution_thread_pool_size = 1; // Number of CPU threads that can work on convolution at the same time.
+// convolution_ir_files = "<comma-separated list of full path names to impulse response files, for sample rates of 44100 and 48000 -- stereo or mono only>";
+// convolution_gain = -4.0; // Static gain applied after convolution, between -60 and +18 dB. Useful for preventing clipping or amplifying low convolution output levels
+// convolution_max_length_in_seconds = 1.0; // Truncate the input file to this length in order to save CPU.
+//
//////////////////////////////////////////
// This loudness filter is used to compensate for human ear non linearity.
-// When the volume decreases, our ears loose more sentisitivity in the low range frequencies than in the mid range ones.
-// This filter aims at compensating for this loss, applying a variable gain to low frequencies depending on the volume.
+// When the volume decreases, our ears lose more sensitivity to low frequencies than to mid range ones.
+// This filter aims at compensating for this loss, applying a small extra gain to low frequencies, depending on the volume.
// More info can be found here: https://en.wikipedia.org/wiki/Equal-loudness_contour
-// For this filter to work properly, you should disable (or set to a fix value) all other volume control and only let shairport-sync control your volume.
-// The setting "loudness_reference_volume_db" should be set at the volume reported by shairport-sync when listening to music at a normal listening volume.
+//
+// NOTE: To use this filter, Shairport Sync must _not_ be controlling a hardware mixer.
+//
+// The setting "loudness_reference_volume_db" is the highest level at which the loudness filter will take effect. Above this level the filter is off.
//////////////////////////////////////////
//
-// loudness = "no"; // Set this to "yes" to activate the loudness filter
-// loudness_reference_volume_db = -20.0; // Above this level the filter will have no effect anymore. Below this level it will gradually boost the low frequencies.
-
+// loudness_enabled = "no"; // Set this to "yes" to activate the loudness filter (only works if you are _not_ using hardware volume control)
+// loudness_reference_volume_db = -16.0; // Above this level the filter will have no effect. Below this level it will gradually boost the low frequencies.
};
// How to deal with metadata, including artwork
#ifdef CONFIG_AIRPLAY_2
int has_fltp_capable_aac_decoder(void) {
-
// return 1 if the AAC decoder advertises fltp decoding capability, which
// is needed for decoding Buffered Audio streams
+ debug(3, "checking availability of an fltp-capable aac decoder");
int has_capability = 0;
const AVCodec *codec = avcodec_find_decoder(AV_CODEC_ID_AAC);
if (codec != NULL) {
- const enum AVSampleFormat *p = codec->sample_fmts;
- if (p != NULL) {
- while ((has_capability == 0) && (*p != AV_SAMPLE_FMT_NONE)) {
- if (*p == AV_SAMPLE_FMT_FLTP)
+ const enum AVSampleFormat *formats = NULL;
+#if LIBAVCODEC_VERSION_INT >= AV_VERSION_INT(61, 13, 100)
+ // New API (FFmpeg 7.1+) for getting formats
+ debug(3, "getting sample formats the new way");
+ int format_count = 0;
+ if ((avcodec_get_supported_config(NULL, codec, AV_CODEC_CONFIG_SAMPLE_FORMAT, 0,
+ (const void **)&formats, &format_count) < 0) ||
+ (format_count == 0))
+ formats = NULL; // not clear if the returned pointer is nulled on error or on zero items
+#else
+ debug(3, "getting sample formats the old way");
+ // older API
+ formats = codec->sample_fmts;
+#endif
+ if (formats != NULL) {
+ while ((has_capability == 0) && (*formats != AV_SAMPLE_FMT_NONE)) {
+ if (*formats == AV_SAMPLE_FMT_FLTP) {
has_capability = 1;
- p++;
+ }
+ formats++;
}
}
+ } else {
+ debug(3, "no AAC codec found.");
}
return has_capability;
}
#endif
#ifdef CONFIG_CONVOLUTION
- config.convolution_max_length = 8192;
+ config.convolution_max_length_in_seconds = 1.0;
+ config.convolution_gain = -4.0;
+ config.convolution_threads = 1; // This is to merely to minimise potential power supply noise some
+ // CPUs make switching cores on and off. E.g. Pi 3.
#endif
- config.loudness_reference_volume_db = -20;
+ config.loudness_reference_volume_db = -16;
#ifdef CONFIG_METADATA_HUB
config.cover_art_cache_dir = "/tmp/shairport-sync/.cache/coverart";
#endif
// config_setting_t *setting;
- const char *str = 0;
+ const char *str = NULL;
int value = 0;
double dvalue = 0.0;
if (config_file_real_path == NULL) {
debug(2, "can't resolve the configuration file \"%s\".", config.configfile);
} else {
- debug(2, "looking for configuration file at full path \"%s\"", config_file_real_path);
+ debug(1, "looking for configuration file at full path \"%s\"", config_file_real_path);
/* Read the file. If there is an error, report it and exit. */
if (config_read_file(&config_file_stuff, config_file_real_path)) {
config_set_auto_convert(&config_file_stuff,
}
#ifdef CONFIG_CONVOLUTION
+
if (config_lookup_string(config.cfg, "dsp.convolution", &str)) {
if (strcasecmp(str, "no") == 0)
- config.convolution = 0;
+ config.convolution_enabled = 0;
else if (strcasecmp(str, "yes") == 0) {
- config.convolution = 1;
-#ifdef CONFIG_AIRPLAY_2
- inform("Note that convolution currently only works on audio received at 44,100 frames "
- "per second.");
-#endif
+ config.convolution_enabled = 1;
+ }
+ warn("the \"dsp\" \"convolution\" setting is deprecated and will be removed due to its "
+ "potential ambiguity. Please use \"convolution_enabled\" instead.");
+ }
+
+ if (config_lookup_string(config.cfg, "dsp.convolution_enabled", &str)) {
+ if (strcasecmp(str, "no") == 0)
+ config.convolution_enabled = 0;
+ else if (strcasecmp(str, "yes") == 0) {
+ config.convolution_enabled = 1;
} else
- die("Invalid dsp.convolution setting \"%s\". It should be \"yes\" or \"no\"", str);
+ die("Invalid dsp.convolution_enabled setting \"%s\". It should be \"yes\" or \"no\"",
+ str);
+ }
+
+ if (config_lookup_int(config.cfg, "dsp.convolution_thread_pool_size", &value)) {
+ if ((value >= 1) && (value <= 64)) {
+ config.convolution_threads = value;
+ } else {
+ warn("Invalid value \"%u\" for \"convolution_thread_pool_size\". It must be between 1 and 64."
+ "The default of %u will be used instead.",
+ value, config.convolution_threads);
+ }
}
if (config_lookup_float(config.cfg, "dsp.convolution_gain", &dvalue)) {
config.convolution_gain = dvalue;
- if (dvalue > 10 || dvalue < -50)
- die("Invalid value \"%f\" for dsp.convolution_gain. It should be between -50 and +10 dB",
+ if (dvalue > 18 || dvalue < -60)
+ die("Invalid value \"%f\" for dsp.convolution_gain. It should be between -60 and +18 dB",
dvalue);
}
if (config_lookup_int(config.cfg, "dsp.convolution_max_length", &value)) {
- config.convolution_max_length = value;
-
+ config.convolution_max_length_in_seconds = (double)value / 44100;
+ warn("the \"dsp\" \"convolution_max_length\" setting is deprecated, as it assumes a fixed "
+ "sample rate of 44,100. It will be removed. "
+ "Please use convolution_max_length_in_seconds instead.");
if (value < 1 || value > 200000)
die("dsp.convolution_max_length must be within 1 and 200000");
}
+ if (config_lookup_float(config.cfg, "dsp.convolution_max_length_in_seconds", &dvalue)) {
+
+ if (dvalue > 20 || dvalue < 0) {
+ warn("Invalid value \"%f\" for dsp.convolution_max_length_in_seconds -- ignored. It "
+ "should be between 0 and 20. It is set to %f.1.",
+ dvalue, config.convolution_max_length_in_seconds);
+ } else {
+ config.convolution_max_length_in_seconds = dvalue;
+ }
+ }
+
if (config_lookup_non_empty_string(config.cfg, "dsp.convolution_ir_file", &str)) {
- config.convolution_ir_file = strdup(str);
- config.convolver_valid =
- convolver_init(config.convolution_ir_file, config.convolution_max_length);
+ warn(
+ "the \"dsp\" \"convolution_ir_file\" setting is deprecated and will be removed. Please "
+ "use \"convolution_ir_files\" instead, which allows multiple comma-separated files.");
+ config.convolution_ir_files = parse_ir_filenames(str, &config.convolution_ir_file_count);
}
- if (config.convolution && config.convolution_ir_file == NULL) {
- warn("Convolution enabled but no convolution_ir_file provided");
+ if (config_lookup_non_empty_string(config.cfg, "dsp.convolution_ir_files", &str)) {
+ config.convolution_ir_files = parse_ir_filenames(str, &config.convolution_ir_file_count);
}
#endif
+
if (config_lookup_string(config.cfg, "dsp.loudness", &str)) {
if (strcasecmp(str, "no") == 0)
- config.loudness = 0;
+ config.loudness_enabled = 0;
else if (strcasecmp(str, "yes") == 0) {
- config.loudness = 1;
-#ifdef CONFIG_AIRPLAY_2
- inform("Note that the loudness filter currently only works on audio received at 44,100 "
- "frames per second.");
-#endif
+ config.loudness_enabled = 1;
+ }
+ warn("the \"dsp\" \"loudness\" setting is deprecated and will be removed due to its "
+ "potential ambiguity. Please use \"loudness_enabled\" instead.");
+ }
+
+ if (config_lookup_string(config.cfg, "dsp.loudness_enabled", &str)) {
+ if (strcasecmp(str, "no") == 0)
+ config.loudness_enabled = 0;
+ else if (strcasecmp(str, "yes") == 0) {
+ config.loudness_enabled = 1;
} else
- die("Invalid dsp.loudness \"%s\". It should be \"yes\" or \"no\"", str);
+ die("Invalid dsp.loudness_enabled \"%s\". It should be \"yes\" or \"no\"", str);
}
if (config_lookup_float(config.cfg, "dsp.loudness_reference_volume_db", &dvalue)) {
dvalue);
}
- if (config.loudness == 1 &&
+ if (config.loudness_enabled == 1 &&
config_lookup_non_empty_string(config.cfg, "alsa.mixer_control_name", &str))
- die("Loudness activated but hardware volume is active. You must remove "
- "\"alsa.mixer_control_name\" to use the loudness filter.");
+ die("The loudness filter is activated but cannot be used because the volume is being "
+ "controlled by a hardware mixer. "
+ "You must not use a hardware mixer when using the loudness filter.");
} else {
if (config_error_type(&config_file_stuff) == CONFIG_ERR_FILE_IO)
#endif
#ifdef CONFIG_CONVOLUTION
- if (config.convolution_ir_file)
- free(config.convolution_ir_file);
+ if (config.convolution_ir_files) {
+ free_ir_filenames(config.convolution_ir_files, config.convolution_ir_file_count);
+ config.convolution_ir_files = NULL;
+ config.convolution_ir_file_count = 0;
+ }
+ convolver_pool_closedown();
#endif
if (config.regtype)
#endif
#ifdef CONFIG_CONVOLUTION
- debug(option_print_level, "convolution is %d.", config.convolution);
- debug(option_print_level, "convolution IR file is \"%s\"", config.convolution_ir_file);
- debug(option_print_level, "convolution max length %d", config.convolution_max_length);
+ debug(option_print_level, "convolution_enabled is %s.",
+ config.convolution_enabled != 0 ? "true" : "false");
+ debug(option_print_level, "convolution maximum length is %f seconds.",
+ config.convolution_max_length_in_seconds);
debug(option_print_level, "convolution gain is %f", config.convolution_gain);
+ sanity_check_ir_files(option_print_level, config.convolution_ir_files,
+ config.convolution_ir_file_count);
#endif
- debug(option_print_level, "loudness is %d.", config.loudness);
+ debug(option_print_level, "loudness_enabled is %s.",
+ config.loudness_enabled != 0 ? "true" : "false");
debug(option_print_level, "loudness reference level is %f", config.loudness_reference_volume_db);
#ifdef CONFIG_SOXR
send_ssnc_metadata('svna', config.service_name, strlen(config.service_name), 1);
#endif
+#ifdef CONFIG_CONVOLUTION
+ convolver_pool_init(config.convolution_threads, 8); // 8 channels
+#endif
activity_monitor_start();
debug(3, "create an RTSP listener");
named_pthread_create(&rtsp_listener_thread, NULL, &rtsp_listen_loop, NULL, "bonjour");
projects / thirdparty / shairport-sync.git / commitdiff
