Dima353/openal-soft
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Use a FIR filter for the UHJ all-pass

Browse Source
This commit is contained in:
Chris Robinson
2020-05-08 15:38:14 -07:00
parent 12bb5a47cd
commit 1c320c3532
3 changed files with 141 additions and 82 deletions
Download Patch File Download Diff File Expand all files Collapse all files
alc/alc.cpp
+3
View File
@@ -2084,6 +2084,9 @@ static ALCenum UpdateDeviceParams(ALCdevice *device, const int *attrList)
device->SourcesMax, device->NumMonoSources, device->NumStereoSources,
device->AuxiliaryEffectSlotMax, device->NumAuxSends);
if(Uhj2Encoder *uhj{device->Uhj_Encoder.get()})
device->FixedLatency += nanoseconds{seconds{uhj->sFilterSize}} / device->Frequency;
/* Enable the stablizer only for formats that have front-left, front-right,
* and front-center outputs.
*/
alc/uhjfilter.cpp
+116 -55
View File
@@ -3,50 +3,99 @@
#include "uhjfilter.h"
#ifdef HAVE_SSE_INTRINSICS
#include <xmmintrin.h>
#endif
#include <algorithm>
#include <iterator>
#include "AL/al.h"
#include "alcomplex.h"
#include "alnumeric.h"
#include "opthelpers.h"
namespace {
/* This is the maximum number of samples processed for each inner loop
* iteration. */
#define MAX_UPDATE_SAMPLES 128
using complex_d = std::complex<double>;
constexpr std::array<float,4> Filter1CoeffSqr{{
0.479400865589f, 0.876218493539f, 0.976597589508f, 0.997499255936f
}};
constexpr std::array<float,4> Filter2CoeffSqr{{
0.161758498368f, 0.733028932341f, 0.945349700329f, 0.990599156685f
}};
void allpass_process(al::span<AllPassState,4> state, float *dst, const float *src,
const std::array<float,4> &coeffs, const size_t todo)
std::array<float,Uhj2Encoder::sFilterSize> GenerateFilter()
{
const std::array<float,4> aa{coeffs};
std::array<std::array<float,2>,4> z{{state[0].z, state[1].z, state[2].z, state[3].z}};
auto proc_sample = [aa,&z](float sample) noexcept -> float
/* Some notes on this filter construction.
*
* An impulse in the frequency domain is represented by a continuous series
* of +1,-1 values, with a 0 imaginary term. Consequently, that impulse
* with a +90 degree phase offset would be represented by 0s with imaginary
* terms that alternate between +1,-1. Converting that to the time domain
* results in a FIR filter that can be convolved with the incoming signal
* to apply a wide-band 90-degree phase shift.
*
* A particularly notable aspect of the time-domain filter response is that
* every other coefficient is 0. This allows doubling the effective size of
* the filter, by only storing the non-0 coefficients and double-stepping
* over the input to apply it.
*
* Additionally, the resulting filter is independent of the sample rate.
* The same filter can be applied regardless of the device's sample rate
* and achieve the same effect, although a lower rate allows the filter to
* cover more time and improve the results.
*/
constexpr complex_d c0{0.0, 1.0};
constexpr complex_d c1{0.0, -1.0};
constexpr size_t half_size{32768};
/* Generate a frequency domain impulse with a +90 degree phase offset. Keep
* the latter half clear for converting to the time domain.
*/
auto fftBuffer = std::vector<complex_d>(half_size*2, complex_d{});
for(size_t i{0};i < half_size;i += 2)
{
for(size_t i{0};i < 4;++i)
fftBuffer[i ] = c0;
fftBuffer[i+1] = c1;
}
complex_fft(fftBuffer, 1.0);
/* Reverse and truncate the filter to a usable size, and store only the
* non-0 terms. Should this be windowed?
*/
std::array<float,Uhj2Encoder::sFilterSize> ret;
auto fftiter = fftBuffer.data() + half_size + (Uhj2Encoder::sFilterSize-1);
for(float &coeff : ret)
{
coeff = static_cast<float>(fftiter->real() / half_size);
fftiter -= 2;
}
return ret;
}
const auto PShiftCoeffs = GenerateFilter();
void allpass_process(al::span<float> dst, const float *RESTRICT src)
{
for(float &output : dst)
{
#ifdef HAVE_SSE_INTRINSICS
constexpr size_t todo{PShiftCoeffs.size()>>2};
__m128 r4{_mm_setzero_ps()};
for(size_t i{0};i < todo;i+=4)
{
const float output{sample*aa[i] + z[i][0]};
z[i][0] = z[i][1];
z[i][1] = output*aa[i] - sample;
sample = output;
const __m128 coeffs{_mm_load_ps(&PShiftCoeffs[i])};
const __m128 s{_mm_setr_ps(src[i*2], src[i*2 + 2], src[i*2 + 4], src[i*2 + 6])};
r4 = _mm_add_ps(r4, _mm_mul_ps(s, coeffs));
}
return sample;
};
std::transform(src, src+todo, dst, proc_sample);
state[0].z = z[0];
state[1].z = z[1];
state[2].z = z[2];
state[3].z = z[3];
r4 = _mm_add_ps(r4, _mm_shuffle_ps(r4, r4, _MM_SHUFFLE(0, 1, 2, 3)));
r4 = _mm_add_ps(r4, _mm_movehl_ps(r4, r4));
float ret{_mm_cvtss_f32(r4)};
#else
float ret{0.0f};
for(size_t i{0};i < PShiftCoeffs.size();++i)
ret += src[i*2] * PShiftCoeffs[i];
#endif
output += ret;
++src;
}
}
} // namespace
@@ -73,40 +122,52 @@ void allpass_process(al::span<AllPassState,4> state, float *dst, const float *sr
*/
void Uhj2Encoder::encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
FloatBufferLine *InSamples, const size_t SamplesToDo)
const FloatBufferLine *InSamples, const size_t SamplesToDo)
{
ASSUME(SamplesToDo > 0);
const auto winput = al::assume_aligned<16>(InSamples[0].cbegin());
const auto xinput = al::assume_aligned<16>(InSamples[1].cbegin());
const auto yinput = al::assume_aligned<16>(InSamples[2].cbegin());
/* D = 0.6554516*Y */
std::transform(yinput, yinput+SamplesToDo, mTemp.begin(),
[](const float y) noexcept -> float { return 0.6554516f*y; });
/* NOTE: Filter1 requires a 1 sample delay for the final output, so take
* the last processed sample from the previous run as the first output
* sample.
*/
mSide[0] = mLastY;
allpass_process(mFilter1_Y, mSide.data()+1, mTemp.data(), Filter1CoeffSqr, SamplesToDo);
mLastY = mSide[SamplesToDo];
/* D += j(-0.3420201*W + 0.5098604*X) */
std::transform(winput, winput+SamplesToDo, xinput, mTemp.begin(),
[](const float w, const float x) noexcept -> float
{ return -0.3420201f*w + 0.5098604f*x; });
allpass_process(mFilter2_WX, mTemp.data(), mTemp.data(), Filter2CoeffSqr, SamplesToDo);
for(size_t i{0};i < SamplesToDo;++i)
mSide[i] += mTemp[i];
const float *RESTRICT winput{al::assume_aligned<16>(InSamples[0].data())};
const float *RESTRICT xinput{al::assume_aligned<16>(InSamples[1].data())};
const float *RESTRICT yinput{al::assume_aligned<16>(InSamples[2].data())};
/* S = 0.9396926*W + 0.1855740*X */
std::transform(winput, winput+SamplesToDo, xinput, mTemp.begin(),
std::transform(winput, winput+SamplesToDo, xinput, mMid.begin(),
[](const float w, const float x) noexcept -> float
{ return 0.9396926f*w + 0.1855740f*x; });
mMid[0] = mLastWX;
allpass_process(mFilter1_WX, mMid.data()+1, mTemp.data(), Filter1CoeffSqr, SamplesToDo);
mLastWX = mMid[SamplesToDo];
/* D = 0.6554516*Y */
std::transform(yinput, yinput+SamplesToDo, mSide.begin(),
[](const float y) noexcept -> float { return 0.6554516f*y; });
/* Apply a delay to the non-filtered signal to align with the filter delay. */
if LIKELY(SamplesToDo >= sFilterSize)
{
auto buffer_end = mMid.begin() + SamplesToDo;
auto delay_end = std::rotate(mMid.begin(), buffer_end - sFilterSize, buffer_end);
std::swap_ranges(mMid.begin(), delay_end, mMidDelay.begin());
buffer_end = mSide.begin() + SamplesToDo;
delay_end = std::rotate(mSide.begin(), buffer_end - sFilterSize, buffer_end);
std::swap_ranges(mSide.begin(), delay_end, mSideDelay.begin());
}
else
{
auto buffer_end = mMid.begin() + SamplesToDo;
auto delay_start = std::swap_ranges(mMid.begin(), buffer_end, mMidDelay.begin());
std::rotate(mMidDelay.begin(), delay_start, mMidDelay.end());
buffer_end = mSide.begin() + SamplesToDo;
delay_start = std::swap_ranges(mSide.begin(), buffer_end, mSideDelay.begin());
std::rotate(mSideDelay.begin(), delay_start, mSideDelay.end());
}
/* D += j(-0.3420201*W + 0.5098604*X) */
auto tmpiter = std::copy(mSideHistory.cbegin(), mSideHistory.cend(), mTemp.begin());
std::transform(winput, winput+SamplesToDo, xinput, tmpiter,
[](const float w, const float x) noexcept -> float
{ return -0.3420201f*w + 0.5098604f*x; });
std::copy_n(mTemp.cbegin()+SamplesToDo, mSideHistory.size(), mSideHistory.begin());
allpass_process({mSide.data(), SamplesToDo}, mTemp.data());
/* Left = (S + D)/2.0 */
float *RESTRICT left{al::assume_aligned<16>(LeftOut.data())};
alc/uhjfilter.h
+22 -27
View File
@@ -7,10 +7,6 @@
#include "almalloc.h"
struct AllPassState {
std::array<float,2> z{{0.0f, 0.0f}};
};
/* Encoding 2-channel UHJ from B-Format is done as:
*
* S = 0.9396926*W + 0.1855740*X
@@ -21,36 +17,35 @@ struct AllPassState {
*
* where j is a wide-band +90 degree phase shift.
*
* The phase shift is done using a Hilbert transform, described here:
* https://web.archive.org/web/20060708031958/http://www.biochem.oulu.fi/~oniemita/dsp/hilbert/
* It works using 2 sets of 4 chained filters. The first filter chain produces
* a phase shift of varying magnitude over a wide range of frequencies, while
* the second filter chain produces a phase shift 90 degrees ahead of the
* first over the same range.
*
* Combining these two stages requires the use of three filter chains. S-
* channel output uses a Filter1 chain on the W and X channel mix, while the D-
* channel output uses a Filter1 chain on the Y channel plus a Filter2 chain on
* the W and X channel mix. This results in the W and X input mix on the D-
* channel output having the required +90 degree phase shift relative to the
* other inputs.
* The phase shift is done using a FIR filter derived from an FFT'd impulse
* with the desired shift.
*/
struct Uhj2Encoder {
alignas(16) std::array<float,BUFFERSIZE> mTemp;
alignas(16) std::array<float,BUFFERSIZE+1> mMid;
alignas(16) std::array<float,BUFFERSIZE+1> mSide;
/* A particular property of the filter allows it to cover nearly twice its
* length, so the filter size is also the effective delay (despite being
* center-aligned).
*/
constexpr static size_t sFilterSize{128};
AllPassState mFilter1_Y[4];
AllPassState mFilter2_WX[4];
AllPassState mFilter1_WX[4];
float mLastY{0.0f}, mLastWX{0.0f};
/* Delays for the unfiltered signal. */
alignas(16) std::array<float,sFilterSize> mMidDelay;
alignas(16) std::array<float,sFilterSize> mSideDelay;
/* Encodes a 2-channel UHJ (stereo-compatible) signal from a B-Format input
/* History for the FIR filter. */
alignas(16) std::array<float,sFilterSize*2 - 1> mSideHistory;
alignas(16) std::array<float,BUFFERSIZE + sFilterSize*2> mTemp;
alignas(16) std::array<float,BUFFERSIZE> mMid;
alignas(16) std::array<float,BUFFERSIZE> mSide;
/**
* Encodes a 2-channel UHJ (stereo-compatible) signal from a B-Format input
* signal. The input must use FuMa channel ordering and scaling.
*/
void encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut, FloatBufferLine *InSamples,
const size_t SamplesToDo);
void encode(FloatBufferLine &LeftOut, FloatBufferLine &RightOut,
const FloatBufferLine *InSamples, const size_t SamplesToDo);
DEF_NEWDEL(Uhj2Encoder)
};