aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--alc/alc.cpp3
-rw-r--r--alc/uhjfilter.cpp161
-rw-r--r--alc/uhjfilter.h49
3 files changed, 136 insertions, 77 deletions
diff --git a/alc/alc.cpp b/alc/alc.cpp
index 040bffa6..6b4bde3f 100644
--- a/alc/alc.cpp
+++ b/alc/alc.cpp
@@ -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.
*/
diff --git a/alc/uhjfilter.cpp b/alc/uhjfilter.cpp
index 13d44130..99737cc9 100644
--- a/alc/uhjfilter.cpp
+++ b/alc/uhjfilter.cpp
@@ -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>;
+
+std::array<float,Uhj2Encoder::sFilterSize> GenerateFilter()
+{
+ /* 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)
+ {
+ 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();
-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)
+void allpass_process(al::span<float> dst, const float *RESTRICT src)
{
- 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
+ for(float &output : dst)
{
- for(size_t i{0};i < 4;++i)
+#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());
+ 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, mMid.begin(),
+ [](const float w, const float x) noexcept -> float
+ { return 0.9396926f*w + 0.1855740f*x; });
/* D = 0.6554516*Y */
- std::transform(yinput, yinput+SamplesToDo, mTemp.begin(),
+ std::transform(yinput, yinput+SamplesToDo, mSide.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];
+
+ /* 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) */
- std::transform(winput, winput+SamplesToDo, xinput, mTemp.begin(),
+ 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; });
- allpass_process(mFilter2_WX, mTemp.data(), mTemp.data(), Filter2CoeffSqr, SamplesToDo);
- for(size_t i{0};i < SamplesToDo;++i)
- mSide[i] += mTemp[i];
-
- /* S = 0.9396926*W + 0.1855740*X */
- std::transform(winput, winput+SamplesToDo, xinput, mTemp.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];
+ 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())};
diff --git a/alc/uhjfilter.h b/alc/uhjfilter.h
index 0593cdb9..db8e55ec 100644
--- a/alc/uhjfilter.h
+++ b/alc/uhjfilter.h
@@ -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};
+
+ /* Delays for the unfiltered signal. */
+ alignas(16) std::array<float,sFilterSize> mMidDelay;
+ alignas(16) std::array<float,sFilterSize> mSideDelay;
+
+ /* History for the FIR filter. */
+ alignas(16) std::array<float,sFilterSize*2 - 1> mSideHistory;
+
+ alignas(16) std::array<float,BUFFERSIZE + sFilterSize*2> mTemp;
- AllPassState mFilter1_Y[4];
- AllPassState mFilter2_WX[4];
- AllPassState mFilter1_WX[4];
- float mLastY{0.0f}, mLastWX{0.0f};
+ 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
+ /**
+ * 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)
};
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-30 06:26:15 +0000