diff options
Diffstat (limited to 'alc/uhjfilter.cpp')
| -rw-r--r-- | alc/uhjfilter.cpp | 161 |
1 files changed, 111 insertions, 50 deletions
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())}; |