Mike Brady
f30dd7aee6
Enable the AirPlay 2 build to operate an optional classic-AirPlay-only service or to gracefully degrade to classic AirPlay, as follows:
1. Add a new command-line option '-—service-type=<type>' and an equivalent configuration entry 'service_type = "<type>"' in the 'general' section of the configuration file,
where <type> can be "auto", "classic" or "airplay2":
1. "auto" (default) means that the service will be AirPlay 2 if NQPTP is running. If NQPTP is not running, classic AirPlay service will be provided instead, and in that case,
"(Classic)" will be appended to the default AirPlay service name visible to AirPlay clients like Apple Music, for example "RaspberryPi3B (Classic)".
2. "classic" means the service will be classic AirPlay (aka AirPlay 1).
3. "airplay2" means the service will be the modern AirPlay 2. In this case, as distinct from "auto", if NQPTP is not running, Shairport Sync will log an error and terminate.
2. In the systemd service file, NQPTP is now a "Want" rather than a "Require". If it's present, then it will be launched before Shairport Sync. If it's absent, Shairport Sync will launch anyway.
3. Improve the delivery of input format changes and emit 'sdsc' metadata when changes occur.
4. Emit new format information in the log if statistics is enabled.
Docker Changes
1. Support for linux/arm/v6 has been dropped, as Docker is no longer supported.
2. NQPTP is not started in the AirPlay 2 Docker image if '--service-type=classic' or '--service-type=airplay1' is in the command line options at the end of the docker run command.
The purpose is to ensure that ports 319 and 320 are left alone when the AirPlay 2 image is set to provide Classic service only.
Note that setting the configuration file 'service_type' to 'classic' will not prevent NQPTP from starting up -- you must use the command line option.
3. A new "dev" target has been added. It is a large image containing the custom-built FFmpeg library, NQPTP, Avahi and D-Bus along with the Shairport Sync source and
all necessary development tools. When started, Avahi, D-Bus and NQPTP are all installed and running. The bash shell has also been added and is entered.
Stability Improvements
Reorganise session preemption to fully terminate the existing session before starting a new one.
Don't delay closing the event port to wait for it to be closed at the client end.
Add a safe_socket_close() function to ensure sockets are fully closed. Use -1 to designate closed rather than 0, to prevent attempts to reclose sockets, causing mayhem.
Re-order FFmpeg decommissioning during teardown.
910 lines
25 KiB
C++
910 lines
25 KiB
C++
// ==================================================================================
|
|
// Copyright (c) 2017 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 detail {
|
|
|
|
class AudioFFTImpl {
|
|
public:
|
|
AudioFFTImpl() = default;
|
|
AudioFFTImpl(const AudioFFTImpl &) = delete;
|
|
AudioFFTImpl &operator=(const AudioFFTImpl &) = delete;
|
|
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;
|
|
};
|
|
|
|
constexpr 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);
|
|
}
|
|
}
|
|
|
|
} // End of namespace detail
|
|
|
|
// ================================================================
|
|
|
|
#ifdef AUDIOFFT_OOURA_USED
|
|
|
|
/**
|
|
* @internal
|
|
* @class OouraFFT
|
|
* @brief FFT implementation based on the great radix-4 routines by Takuya Ooura
|
|
*/
|
|
class OouraFFT : public detail::AudioFFTImpl {
|
|
public:
|
|
OouraFFT() : detail::AudioFFTImpl(), _size(0), _ip(), _w(), _buffer() {}
|
|
|
|
OouraFFT(const OouraFFT &) = delete;
|
|
OouraFFT &operator=(const OouraFFT &) = delete;
|
|
|
|
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
|
|
detail::ConvertBuffer(_buffer.data(), data, _size);
|
|
|
|
rdft(static_cast<int>(_size), +1, _buffer.data(), _ip.data(), _w.data());
|
|
|
|
// Convert back to split-complex
|
|
{
|
|
double *b = _buffer.data();
|
|
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.data();
|
|
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
|
|
detail::ScaleBuffer(data, _buffer.data(), 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];
|
|
}
|
|
};
|
|
|
|
/**
|
|
* @internal
|
|
* @brief Concrete FFT implementation
|
|
*/
|
|
typedef OouraFFT AudioFFTImplementation;
|
|
|
|
#endif // AUDIOFFT_OOURA_USED
|
|
|
|
// ================================================================
|
|
|
|
#ifdef AUDIOFFT_APPLE_ACCELERATE_USED
|
|
|
|
/**
|
|
* @internal
|
|
* @class AppleAccelerateFFT
|
|
* @brief FFT implementation using the Apple Accelerate framework internally
|
|
*/
|
|
class AppleAccelerateFFT : public detail::AudioFFTImpl {
|
|
public:
|
|
AppleAccelerateFFT()
|
|
: detail::AudioFFTImpl(), _size(0), _powerOf2(0), _fftSetup(0), _re(), _im() {}
|
|
|
|
AppleAccelerateFFT(const AppleAccelerateFFT &) = delete;
|
|
AppleAccelerateFFT &operator=(const AppleAccelerateFFT &) = delete;
|
|
|
|
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;
|
|
};
|
|
|
|
/**
|
|
* @internal
|
|
* @brief Concrete FFT implementation
|
|
*/
|
|
typedef AppleAccelerateFFT AudioFFTImplementation;
|
|
|
|
#endif // AUDIOFFT_APPLE_ACCELERATE_USED
|
|
|
|
// ================================================================
|
|
|
|
#ifdef AUDIOFFT_FFTW3_USED
|
|
|
|
/**
|
|
* @internal
|
|
* @class FFTW3FFT
|
|
* @brief FFT implementation using FFTW3 internally (see fftw.org)
|
|
*/
|
|
class FFTW3FFT : public detail::AudioFFTImpl {
|
|
public:
|
|
FFTW3FFT()
|
|
: detail::AudioFFTImpl(), _size(0), _complexSize(0), _planForward(0), _planBackward(0),
|
|
_data(0), _re(0), _im(0) {}
|
|
|
|
FFTW3FFT(const FFTW3FFT &) = delete;
|
|
FFTW3FFT &operator=(const FFTW3FFT &) = delete;
|
|
|
|
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)));
|
|
fftwf_set_timelimit(0.01);
|
|
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));
|
|
}
|
|
|
|
virtual void ifft(float *data, const float *re, const float *im) override {
|
|
::memcpy(_re, re, _complexSize * sizeof(float));
|
|
::memcpy(_im, im, _complexSize * sizeof(float));
|
|
fftwf_execute_split_dft_c2r(_planBackward, _re, _im, _data);
|
|
detail::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;
|
|
};
|
|
|
|
/**
|
|
* @internal
|
|
* @brief Concrete FFT implementation
|
|
*/
|
|
typedef FFTW3FFT AudioFFTImplementation;
|
|
|
|
#endif // AUDIOFFT_FFTW3_USED
|
|
|
|
// =============================================================
|
|
|
|
AudioFFT::AudioFFT() : _impl(new AudioFFTImplementation()) {}
|
|
|
|
AudioFFT::~AudioFFT() {}
|
|
|
|
void AudioFFT::init(size_t size) {
|
|
assert(detail::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; }
|
|
|
|
} // namespace audiofft
|