mirror of
https://github.com/LMMS/lmms.git
synced 2026-04-10 00:58:36 -04:00
226 lines
7.2 KiB
C++
226 lines
7.2 KiB
C++
/*
|
|
* HilbertTransform.h
|
|
*
|
|
* Copyright (c) 2025 Lost Robot <r94231/at/gmail/dot/com>
|
|
*
|
|
* This file is part of LMMS - https://lmms.io
|
|
*
|
|
* This program is free software; you can redistribute it and/or
|
|
* modify it under the terms of the GNU General Public
|
|
* License as published by the Free Software Foundation; either
|
|
* version 2 of the License, or (at your option) any later version.
|
|
*
|
|
* 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.
|
|
*
|
|
* You should have received a copy of the GNU General Public
|
|
* License along with this program (see COPYING); if not, write to the
|
|
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
|
|
* Boston, MA 02110-1301 USA.
|
|
*
|
|
* -------------------------------------------------------------------------
|
|
* Full credit for the Hilbert IIR coefficients and general usage goes to:
|
|
*
|
|
* Copyright (c) 2024 Geraint Luff / Signalsmith Audio Ltd.
|
|
* Released under the 0BSD (Zero-Clause BSD) License.
|
|
*
|
|
* https://github.com/Signalsmith-Audio/hilbert-iir/blob/main/hilbert.h
|
|
* -------------------------------------------------------------------------
|
|
*/
|
|
|
|
#ifndef LMMS_HILBERT_TRANSFORM_H
|
|
#define LMMS_HILBERT_TRANSFORM_H
|
|
#include <cmath>
|
|
|
|
#ifdef __SSE2__
|
|
#include <emmintrin.h>
|
|
#endif
|
|
|
|
namespace lmms
|
|
{
|
|
|
|
template<int Channels>
|
|
struct HilbertIIRFloat
|
|
{
|
|
static constexpr int order = 12;
|
|
|
|
alignas(16) static constexpr float baseCoeffsR[order] = {
|
|
-0.000224352093802f, 0.0107500557815f, -0.0456795873917f,
|
|
0.11282500582f, -0.208067578452f, 0.28717837501f,
|
|
-0.254675294431f, 0.0481081835026f, 0.227861357867f,
|
|
-0.365411839137f, 0.280729061131f, -0.0935061787728f
|
|
};
|
|
alignas(16) static constexpr float baseCoeffsI[order] = {
|
|
0.00543499018201f, -0.0173890685681f, 0.0229166931429f,
|
|
0.00278413661237f, -0.104628958675f, 0.33619239719f,
|
|
-0.683033899655f, 0.954061589374f, -0.891273574569f,
|
|
0.525088317271f, -0.155131206606f, 0.00512245855404f
|
|
};
|
|
alignas(16) static constexpr float basePolesR[order] = {
|
|
-0.00495335976478f, -0.017859491302f, -0.0413714373155f,
|
|
-0.0882148408885f, -0.17922965812f, -0.338261800753f,
|
|
-0.557688699732f, -0.735157736148f, -0.719057381172f,
|
|
-0.517871025209f, -0.280197469471f, -0.0852751354531f
|
|
};
|
|
alignas(16) static constexpr float basePolesI[order] = {
|
|
0.0092579876872f, 0.0273493725543f, 0.0744756910287f,
|
|
0.178349677457f, 0.39601340223f, 0.829229533354f,
|
|
1.61298538328f, 2.79987398682f, 4.16396166128f,
|
|
5.29724826804f, 5.99598602388f, 6.3048492377f
|
|
};
|
|
static constexpr float baseDirect = 0.000262057212648f;
|
|
|
|
alignas(16) float coeffsR[order];
|
|
alignas(16) float coeffsI[order];
|
|
alignas(16) float polesR[order];
|
|
alignas(16) float polesI[order];
|
|
alignas(16) float stateR[Channels][order];
|
|
alignas(16) float stateI[Channels][order];
|
|
float direct;
|
|
|
|
HilbertIIRFloat(float sampleRate = 48000.0f, float passbandGain = 2.0f)
|
|
{
|
|
const float freqFactor = std::fmin(0.46f, 20000.0f / sampleRate);
|
|
const float coeffScale = freqFactor * passbandGain;
|
|
direct = baseDirect * 2.0f * passbandGain * freqFactor;
|
|
for (int i = 0; i < order; ++i)
|
|
{
|
|
coeffsR[i] = baseCoeffsR[i] * coeffScale;
|
|
coeffsI[i] = baseCoeffsI[i] * coeffScale;
|
|
const float a = basePolesR[i] * freqFactor;
|
|
const float b = basePolesI[i] * freqFactor;
|
|
const float ea = std::exp(a);
|
|
polesR[i] = ea * std::cos(b);
|
|
polesI[i] = ea * std::sin(b);
|
|
}
|
|
reset();
|
|
}
|
|
|
|
inline void reset()
|
|
{
|
|
for (int ch = 0; ch < Channels; ++ch)
|
|
{
|
|
for (int i = 0; i < order; ++i)
|
|
{
|
|
stateR[ch][i] = stateI[ch][i] = 0.0f;
|
|
}
|
|
}
|
|
}
|
|
|
|
inline void processReal(float x, int channel, float *out)
|
|
{
|
|
float *sR = stateR[channel], *sI = stateI[channel];
|
|
#ifdef __SSE2__
|
|
const __m128 vx = _mm_set1_ps(x);
|
|
__m128 sumR = _mm_setzero_ps(), sumI = _mm_setzero_ps();
|
|
for (int i = 0; i < order; i += 4)
|
|
{
|
|
__m128 vr = _mm_load_ps(&sR[i]);
|
|
__m128 vi = _mm_load_ps(&sI[i]);
|
|
__m128 vpr = _mm_load_ps(&polesR[i]);
|
|
__m128 vpi = _mm_load_ps(&polesI[i]);
|
|
__m128 vcr = _mm_load_ps(&coeffsR[i]);
|
|
__m128 vci = _mm_load_ps(&coeffsI[i]);
|
|
|
|
__m128 rpr = _mm_mul_ps(vr, vpr);
|
|
__m128 impi = _mm_mul_ps(vi, vpi);
|
|
__m128 xcr = _mm_mul_ps(vx, vcr);
|
|
__m128 nr = _mm_add_ps(_mm_sub_ps(rpr, impi), xcr);
|
|
|
|
__m128 rpi = _mm_mul_ps(vr, vpi);
|
|
__m128 impr = _mm_mul_ps(vi, vpr);
|
|
__m128 xci = _mm_mul_ps(vx, vci);
|
|
__m128 ni = _mm_add_ps(_mm_add_ps(rpi, impr), xci);
|
|
|
|
_mm_store_ps(&sR[i], nr);
|
|
_mm_store_ps(&sI[i], ni);
|
|
|
|
sumR = _mm_add_ps(sumR, nr);
|
|
sumI = _mm_add_ps(sumI, ni);
|
|
}
|
|
float tmpR[4], tmpI[4];
|
|
_mm_storeu_ps(tmpR, sumR);
|
|
_mm_storeu_ps(tmpI, sumI);
|
|
out[0] = x * direct + (tmpR[0] + tmpR[1] + tmpR[2] + tmpR[3]);
|
|
out[1] = (tmpI[0] + tmpI[1] + tmpI[2] + tmpI[3]);
|
|
#else
|
|
float sumR = 0.0f, sumI = 0.0f;
|
|
for (int i = 0; i < order; ++i)
|
|
{
|
|
const float r = sR[i], im = sI[i], pr = polesR[i], pi = polesI[i];
|
|
const float nr = r * pr - im * pi + x * coeffsR[i];
|
|
const float ni = r * pi + im * pr + x * coeffsI[i];
|
|
sR[i] = nr; sI[i] = ni;
|
|
sumR += nr; sumI += ni;
|
|
}
|
|
out[0] = x * direct + sumR;
|
|
out[1] = sumI;
|
|
#endif
|
|
}
|
|
|
|
inline void processComplex(const float *x, int channel, float *out)
|
|
{
|
|
const float xr = x[0], xi = x[1];
|
|
float *sR = stateR[channel], *sI = stateI[channel];
|
|
#ifdef __SSE2__
|
|
const __m128 vxr = _mm_set1_ps(xr), vxi = _mm_set1_ps(xi);
|
|
__m128 sumR = _mm_setzero_ps(), sumI = _mm_setzero_ps();
|
|
for (int i = 0; i < order; i += 4)
|
|
{
|
|
__m128 vr = _mm_load_ps(&sR[i]);
|
|
__m128 vi = _mm_load_ps(&sI[i]);
|
|
__m128 vpr = _mm_load_ps(&polesR[i]);
|
|
__m128 vpi = _mm_load_ps(&polesI[i]);
|
|
__m128 vcr = _mm_load_ps(&coeffsR[i]);
|
|
__m128 vci = _mm_load_ps(&coeffsI[i]);
|
|
|
|
__m128 xrcr = _mm_mul_ps(vxr, vcr);
|
|
__m128 xici = _mm_mul_ps(vxi, vci);
|
|
__m128 xrci = _mm_mul_ps(vxr, vci);
|
|
__m128 xicr = _mm_mul_ps(vxi, vcr);
|
|
|
|
__m128 rpr = _mm_mul_ps(vr, vpr);
|
|
__m128 impi = _mm_mul_ps(vi, vpi);
|
|
__m128 rpi = _mm_mul_ps(vr, vpi);
|
|
__m128 impr = _mm_mul_ps(vi, vpr);
|
|
|
|
__m128 nr = _mm_add_ps(_mm_sub_ps(rpr, impi), _mm_sub_ps(xrcr, xici));
|
|
__m128 ni = _mm_add_ps(_mm_add_ps(rpi, impr), _mm_add_ps(xrci, xicr));
|
|
|
|
_mm_store_ps(&sR[i], nr);
|
|
_mm_store_ps(&sI[i], ni);
|
|
|
|
sumR = _mm_add_ps(sumR, nr);
|
|
sumI = _mm_add_ps(sumI, ni);
|
|
}
|
|
float tmpR[4], tmpI[4];
|
|
_mm_storeu_ps(tmpR, sumR);
|
|
_mm_storeu_ps(tmpI, sumI);
|
|
const float sr = tmpR[0] + tmpR[1] + tmpR[2] + tmpR[3];
|
|
const float si = tmpI[0] + tmpI[1] + tmpI[2] + tmpI[3];
|
|
out[0] = xr * direct + sr;
|
|
out[1] = xi * direct + si;
|
|
#else
|
|
float sumR = 0.0f, sumI = 0.0f;
|
|
for (int i = 0; i < order; ++i)
|
|
{
|
|
const float r = sR[i], im = sI[i];
|
|
const float pr = polesR[i], pi = polesI[i], cr = coeffsR[i], ci = coeffsI[i];
|
|
const float xrcr = xr * cr, xici = xi * ci, xrci = xr * ci, xicr = xi * cr;
|
|
const float nr = r * pr - im * pi + (xrcr - xici);
|
|
const float ni = r * pi + im * pr + (xrci + xicr);
|
|
sR[i] = nr; sI[i] = ni;
|
|
sumR += nr; sumI += ni;
|
|
}
|
|
out[0] = xr * direct + sumR;
|
|
out[1] = xi * direct + sumI;
|
|
#endif
|
|
}
|
|
};
|
|
|
|
} // namespace lmms
|
|
|
|
#endif // LMMS_HILBERT_TRANSFORM_H
|