core/uhjfilter.cpp

   1
   2 #include "config.h"
   3
   4 #include "uhjfilter.h"
   5
   6 #include <algorithm>
   7 #include <iterator>
   8
   9 #include "alcomplex.h"
  10 #include "alnumeric.h"
  11 #include "opthelpers.h"
  12 #include "phase_shifter.h"
  13
  14
  15 UhjQualityType UhjDecodeQuality{UhjQualityType::Default};
  16 UhjQualityType UhjEncodeQuality{UhjQualityType::Default};
  17
  18
  19 namespace {
  20
  21 const PhaseShifterT<UhjLength256> PShiftLq{};
  22 const PhaseShifterT<UhjLength512> PShiftHq{};
  23
  24 template<size_t N>
  25 struct GetPhaseShifter;
  26 template<>
  27 struct GetPhaseShifter<UhjLength256> { static auto& Get() noexcept { return PShiftLq; } };
  28 template<>
  29 struct GetPhaseShifter<UhjLength512> { static auto& Get() noexcept { return PShiftHq; } };
  30
  31
  32 constexpr float square(float x) noexcept
  33 { return x*x; }
  34
  35 /* Filter coefficients for the 'base' all-pass IIR, which applies a frequency-
  36  * dependent phase-shift of N degrees. The output of the filter requires a 1-
  37  * sample delay.
  38  */
  39 constexpr std::array<float,4> Filter1Coeff{{
  40     square(0.6923878f), square(0.9360654322959f), square(0.9882295226860f),
  41     square(0.9987488452737f)
  42 }};
  43 /* Filter coefficients for the offset all-pass IIR, which applies a frequency-
  44  * dependent phase-shift of N+90 degrees.
  45  */
  46 constexpr std::array<float,4> Filter2Coeff{{
  47     square(0.4021921162426f), square(0.8561710882420f), square(0.9722909545651f),
  48     square(0.9952884791278f)
  49 }};
  50
  51 } // namespace
  52
  53 void UhjAllPassFilter::process(const al::span<const float,4> coeffs,
  54     const al::span<const float> src, const bool updateState, float *RESTRICT dst)
  55 {
  56     auto state = mState;
  57
  58     auto proc_sample = [&state,coeffs](float x) noexcept -> float
  59     {
  60         for(size_t i{0};i < 4;++i)
  61         {
  62             const float y{x*coeffs[i] + state[i].z[0]};
  63             state[i].z[0] = state[i].z[1];
  64             state[i].z[1] = y*coeffs[i] - x;
  65             x = y;
  66         }
  67         return x;
  68     };
  69     std::transform(src.begin(), src.end(), dst, proc_sample);
  70     if(updateState) LIKELY mState = state;
  71 }
  72
  73
  74 /* Encoding UHJ from B-Format is done as:
  75  *
  76  * S = 0.9396926*W + 0.1855740*X
  77  * D = j(-0.3420201*W + 0.5098604*X) + 0.6554516*Y
  78  *
  79  * Left = (S + D)/2.0
  80  * Right = (S - D)/2.0
  81  * T = j(-0.1432*W + 0.6512*X) - 0.7071068*Y
  82  * Q = 0.9772*Z
  83  *
  84  * where j is a wide-band +90 degree phase shift. 3-channel UHJ excludes Q,
  85  * while 2-channel excludes Q and T.
  86  *
  87  * The phase shift is done using a linear FIR filter derived from an FFT'd
  88  * impulse with the desired shift.
  89  */
  90
  91 template<size_t N>
  92 void UhjEncoder<N>::encode(float *LeftOut, float *RightOut,
  93     const al::span<const float*const,3> InSamples, const size_t SamplesToDo)
  94 {
  95     const auto &PShift = GetPhaseShifter<N>::Get();
  96
  97     ASSUME(SamplesToDo > 0);
  98
  99     const float *RESTRICT winput{al::assume_aligned<16>(InSamples[0])};
 100     const float *RESTRICT xinput{al::assume_aligned<16>(InSamples[1])};
 101     const float *RESTRICT yinput{al::assume_aligned<16>(InSamples[2])};
 102
 103     std::copy_n(winput, SamplesToDo, mW.begin()+sFilterDelay);
 104     std::copy_n(xinput, SamplesToDo, mX.begin()+sFilterDelay);
 105     std::copy_n(yinput, SamplesToDo, mY.begin()+sFilterDelay);
 106
 107     /* S = 0.9396926*W + 0.1855740*X */
 108     for(size_t i{0};i < SamplesToDo;++i)
 109         mS[i] = 0.9396926f*mW[i] + 0.1855740f*mX[i];
 110
 111     /* Precompute j(-0.3420201*W + 0.5098604*X) and store in mD. */
 112     std::transform(winput, winput+SamplesToDo, xinput, mWX.begin() + sWXInOffset,
 113         [](const float w, const float x) noexcept -> float
 114         { return -0.3420201f*w + 0.5098604f*x; });
 115     PShift.process({mD.data(), SamplesToDo}, mWX.data());
 116
 117     /* D = j(-0.3420201*W + 0.5098604*X) + 0.6554516*Y */
 118     for(size_t i{0};i < SamplesToDo;++i)
 119         mD[i] = mD[i] + 0.6554516f*mY[i];
 120
 121     /* Copy the future samples to the front for next time. */
 122     std::copy(mW.cbegin()+SamplesToDo, mW.cbegin()+SamplesToDo+sFilterDelay, mW.begin());
 123     std::copy(mX.cbegin()+SamplesToDo, mX.cbegin()+SamplesToDo+sFilterDelay, mX.begin());
 124     std::copy(mY.cbegin()+SamplesToDo, mY.cbegin()+SamplesToDo+sFilterDelay, mY.begin());
 125     std::copy(mWX.cbegin()+SamplesToDo, mWX.cbegin()+SamplesToDo+sWXInOffset, mWX.begin());
 126
 127     /* Apply a delay to the existing output to align with the input delay. */
 128     auto *delayBuffer = mDirectDelay.data();
 129     for(float *buffer : {LeftOut, RightOut})
 130     {
 131         float *distbuf{al::assume_aligned<16>(delayBuffer->data())};
 132         ++delayBuffer;
 133
 134         float *inout{al::assume_aligned<16>(buffer)};
 135         auto inout_end = inout + SamplesToDo;
 136         if(SamplesToDo >= sFilterDelay) LIKELY
 137         {
 138             auto delay_end = std::rotate(inout, inout_end - sFilterDelay, inout_end);
 139             std::swap_ranges(inout, delay_end, distbuf);
 140         }
 141         else
 142         {
 143             auto delay_start = std::swap_ranges(inout, inout_end, distbuf);
 144             std::rotate(distbuf, delay_start, distbuf + sFilterDelay);
 145         }
 146     }
 147
 148     /* Combine the direct signal with the produced output. */
 149
 150     /* Left = (S + D)/2.0 */
 151     float *RESTRICT left{al::assume_aligned<16>(LeftOut)};
 152     for(size_t i{0};i < SamplesToDo;i++)
 153         left[i] += (mS[i] + mD[i]) * 0.5f;
 154     /* Right = (S - D)/2.0 */
 155     float *RESTRICT right{al::assume_aligned<16>(RightOut)};
 156     for(size_t i{0};i < SamplesToDo;i++)
 157         right[i] += (mS[i] - mD[i]) * 0.5f;
 158 }
 159
 160 /* This encoding implementation uses two sets of four chained IIR filters to
 161  * produce the desired relative phase shift. The first filter chain produces a
 162  * phase shift of varying degrees over a wide range of frequencies, while the
 163  * second filter chain produces a phase shift 90 degrees ahead of the first
 164  * over the same range. Further details are described here:
 165  *
 166  * https://web.archive.org/web/20060708031958/http://www.biochem.oulu.fi/~oniemita/dsp/hilbert/
 167  *
 168  * 2-channel UHJ output requires the use of three filter chains. The S channel
 169  * output uses a Filter1 chain on the W and X channel mix, while the D channel
 170  * output uses a Filter1 chain on the Y channel plus a Filter2 chain on the W
 171  * and X channel mix. This results in the W and X input mix on the D channel
 172  * output having the required +90 degree phase shift relative to the other
 173  * inputs.
 174  */
 175 void UhjEncoderIIR::encode(float *LeftOut, float *RightOut,
 176     const al::span<const float *const, 3> InSamples, const size_t SamplesToDo)
 177 {
 178     ASSUME(SamplesToDo > 0);
 179
 180     const float *RESTRICT winput{al::assume_aligned<16>(InSamples[0])};
 181     const float *RESTRICT xinput{al::assume_aligned<16>(InSamples[1])};
 182     const float *RESTRICT yinput{al::assume_aligned<16>(InSamples[2])};
 183
 184     /* S = 0.9396926*W + 0.1855740*X */
 185     std::transform(winput, winput+SamplesToDo, xinput, mTemp.begin(),
 186         [](const float w, const float x) noexcept { return 0.9396926f*w + 0.1855740f*x; });
 187     mFilter1WX.process(Filter1Coeff, {mTemp.data(), SamplesToDo}, true, mS.data()+1);
 188     mS[0] = mDelayWX; mDelayWX = mS[SamplesToDo];
 189
 190     /* Precompute j(-0.3420201*W + 0.5098604*X) and store in mWX. */
 191     std::transform(winput, winput+SamplesToDo, xinput, mTemp.begin(),
 192         [](const float w, const float x) noexcept { return -0.3420201f*w + 0.5098604f*x; });
 193     mFilter2WX.process(Filter2Coeff, {mTemp.data(), SamplesToDo}, true, mWX.data());
 194
 195     /* Apply filter1 to Y and store in mD. */
 196     mFilter1Y.process(Filter1Coeff, {yinput, SamplesToDo}, SamplesToDo, mD.data()+1);
 197     mD[0] = mDelayY; mDelayY = mD[SamplesToDo];
 198
 199     /* D = j(-0.3420201*W + 0.5098604*X) + 0.6554516*Y */
 200     for(size_t i{0};i < SamplesToDo;++i)
 201         mD[i] = mWX[i] + 0.6554516f*mD[i];
 202
 203     /* Apply the base filter to the existing output to align with the processed
 204      * signal.
 205      */
 206     mFilter1Direct[0].process(Filter1Coeff, {LeftOut, SamplesToDo}, true, mTemp.data()+1);
 207     mTemp[0] = mDirectDelay[0]; mDirectDelay[0] = mTemp[SamplesToDo];
 208
 209     /* Left = (S + D)/2.0 */
 210     float *RESTRICT left{al::assume_aligned<16>(LeftOut)};
 211     for(size_t i{0};i < SamplesToDo;i++)
 212         left[i] = (mS[i] + mD[i])*0.5f + mTemp[i];
 213
 214     mFilter1Direct[1].process(Filter1Coeff, {RightOut, SamplesToDo}, true, mTemp.data()+1);
 215     mTemp[0] = mDirectDelay[1]; mDirectDelay[1] = mTemp[SamplesToDo];
 216
 217     /* Right = (S - D)/2.0 */
 218     float *RESTRICT right{al::assume_aligned<16>(RightOut)};
 219     for(size_t i{0};i < SamplesToDo;i++)
 220         right[i] = (mS[i] - mD[i])*0.5f + mTemp[i];
 221 }
 222
 223
 224 /* Decoding UHJ is done as:
 225  *
 226  * S = Left + Right
 227  * D = Left - Right
 228  *
 229  * W = 0.981532*S + 0.197484*j(0.828331*D + 0.767820*T)
 230  * X = 0.418496*S - j(0.828331*D + 0.767820*T)
 231  * Y = 0.795968*D - 0.676392*T + j(0.186633*S)
 232  * Z = 1.023332*Q
 233  *
 234  * where j is a +90 degree phase shift. 3-channel UHJ excludes Q, while 2-
 235  * channel excludes Q and T.
 236  */
 237 template<size_t N>
 238 void UhjDecoder<N>::decode(const al::span<float*> samples, const size_t samplesToDo,
 239     const bool updateState)
 240 {
 241     static_assert(sInputPadding <= sMaxPadding, "Filter padding is too large");
 242
 243     const auto &PShift = GetPhaseShifter<N>::Get();
 244
 245     ASSUME(samplesToDo > 0);
 246
 247     {
 248         const float *RESTRICT left{al::assume_aligned<16>(samples[0])};
 249         const float *RESTRICT right{al::assume_aligned<16>(samples[1])};
 250         const float *RESTRICT t{al::assume_aligned<16>(samples[2])};
 251
 252         /* S = Left + Right */
 253         for(size_t i{0};i < samplesToDo+sInputPadding;++i)
 254             mS[i] = left[i] + right[i];
 255
 256         /* D = Left - Right */
 257         for(size_t i{0};i < samplesToDo+sInputPadding;++i)
 258             mD[i] = left[i] - right[i];
 259
 260         /* T */
 261         for(size_t i{0};i < samplesToDo+sInputPadding;++i)
 262             mT[i] = t[i];
 263     }
 264
 265     float *RESTRICT woutput{al::assume_aligned<16>(samples[0])};
 266     float *RESTRICT xoutput{al::assume_aligned<16>(samples[1])};
 267     float *RESTRICT youtput{al::assume_aligned<16>(samples[2])};
 268
 269     /* Precompute j(0.828331*D + 0.767820*T) and store in xoutput. */
 270     auto tmpiter = std::copy(mDTHistory.cbegin(), mDTHistory.cend(), mTemp.begin());
 271     std::transform(mD.cbegin(), mD.cbegin()+samplesToDo+sInputPadding, mT.cbegin(), tmpiter,
 272         [](const float d, const float t) noexcept { return 0.828331f*d + 0.767820f*t; });
 273     if(updateState) LIKELY
 274         std::copy_n(mTemp.cbegin()+samplesToDo, mDTHistory.size(), mDTHistory.begin());
 275     PShift.process({xoutput, samplesToDo}, mTemp.data());
 276
 277     /* W = 0.981532*S + 0.197484*j(0.828331*D + 0.767820*T) */
 278     for(size_t i{0};i < samplesToDo;++i)
 279         woutput[i] = 0.981532f*mS[i] + 0.197484f*xoutput[i];
 280     /* X = 0.418496*S - j(0.828331*D + 0.767820*T) */
 281     for(size_t i{0};i < samplesToDo;++i)
 282         xoutput[i] = 0.418496f*mS[i] - xoutput[i];
 283
 284     /* Precompute j*S and store in youtput. */
 285     tmpiter = std::copy(mSHistory.cbegin(), mSHistory.cend(), mTemp.begin());
 286     std::copy_n(mS.cbegin(), samplesToDo+sInputPadding, tmpiter);
 287     if(updateState) LIKELY
 288         std::copy_n(mTemp.cbegin()+samplesToDo, mSHistory.size(), mSHistory.begin());
 289     PShift.process({youtput, samplesToDo}, mTemp.data());
 290
 291     /* Y = 0.795968*D - 0.676392*T + j(0.186633*S) */
 292     for(size_t i{0};i < samplesToDo;++i)
 293         youtput[i] = 0.795968f*mD[i] - 0.676392f*mT[i] + 0.186633f*youtput[i];
 294
 295     if(samples.size() > 3)
 296     {
 297         float *RESTRICT zoutput{al::assume_aligned<16>(samples[3])};
 298         /* Z = 1.023332*Q */
 299         for(size_t i{0};i < samplesToDo;++i)
 300             zoutput[i] = 1.023332f*zoutput[i];
 301     }
 302 }
 303
 304 void UhjDecoderIIR::decode(const al::span<float*> samples, const size_t samplesToDo,
 305     const bool updateState)
 306 {
 307     static_assert(sInputPadding <= sMaxPadding, "Filter padding is too large");
 308
 309     ASSUME(samplesToDo > 0);
 310
 311     {
 312         const float *RESTRICT left{al::assume_aligned<16>(samples[0])};
 313         const float *RESTRICT right{al::assume_aligned<16>(samples[1])};
 314
 315         /* S = Left + Right */
 316         for(size_t i{0};i < samplesToDo;++i)
 317             mS[i] = left[i] + right[i];
 318
 319         /* D = Left - Right */
 320         for(size_t i{0};i < samplesToDo;++i)
 321             mD[i] = left[i] - right[i];
 322     }
 323
 324     float *RESTRICT woutput{al::assume_aligned<16>(samples[0])};
 325     float *RESTRICT xoutput{al::assume_aligned<16>(samples[1])};
 326     float *RESTRICT youtput{al::assume_aligned<16>(samples[2])};
 327
 328     /* Precompute j(0.828331*D + 0.767820*T) and store in xoutput. */
 329     std::transform(mD.cbegin(), mD.cbegin()+samplesToDo, youtput, mTemp.begin(),
 330         [](const float d, const float t) noexcept { return 0.828331f*d + 0.767820f*t; });
 331     mFilter2DT.process(Filter2Coeff, {mTemp.data(), samplesToDo}, updateState, xoutput);
 332
 333     /* Apply filter1 to S and store in mTemp. */
 334     mTemp[0] = mDelayS;
 335     mFilter1S.process(Filter1Coeff, {mS.data(), samplesToDo}, updateState, mTemp.data()+1);
 336     if(updateState) LIKELY mDelayS = mTemp[samplesToDo];
 337
 338     /* W = 0.981532*S + 0.197484*j(0.828331*D + 0.767820*T) */
 339     for(size_t i{0};i < samplesToDo;++i)
 340         woutput[i] = 0.981532f*mTemp[i] + 0.197484f*xoutput[i];
 341     /* X = 0.418496*S - j(0.828331*D + 0.767820*T) */
 342     for(size_t i{0};i < samplesToDo;++i)
 343         xoutput[i] = 0.418496f*mTemp[i] - xoutput[i];
 344
 345
 346     /* Apply filter1 to (0.795968*D - 0.676392*T) and store in mTemp. */
 347     std::transform(mD.cbegin(), mD.cbegin()+samplesToDo, youtput, youtput,
 348         [](const float d, const float t) noexcept { return 0.795968f*d - 0.676392f*t; });
 349     mTemp[0] = mDelayDT;
 350     mFilter1DT.process(Filter1Coeff, {youtput, samplesToDo}, updateState, mTemp.data()+1);
 351     if(updateState) LIKELY mDelayDT = mTemp[samplesToDo];
 352
 353     /* Precompute j*S and store in youtput. */
 354     mFilter2S.process(Filter2Coeff, {mS.data(), samplesToDo}, updateState, youtput);
 355
 356     /* Y = 0.795968*D - 0.676392*T + j(0.186633*S) */
 357     for(size_t i{0};i < samplesToDo;++i)
 358         youtput[i] = mTemp[i] + 0.186633f*youtput[i];
 359
 360
 361     if(samples.size() > 3)
 362     {
 363         float *RESTRICT zoutput{al::assume_aligned<16>(samples[3])};
 364
 365         /* Apply filter1 to Q and store in mTemp. */
 366         mTemp[0] = mDelayQ;
 367         mFilter1Q.process(Filter1Coeff, {zoutput, samplesToDo}, updateState, mTemp.data()+1);
 368         if(updateState) LIKELY mDelayQ = mTemp[samplesToDo];
 369
 370         /* Z = 1.023332*Q */
 371         for(size_t i{0};i < samplesToDo;++i)
 372             zoutput[i] = 1.023332f*mTemp[i];
 373     }
 374 }
 375
 376
 377 /* Super Stereo processing is done as:
 378  *
 379  * S = Left + Right
 380  * D = Left - Right
 381  *
 382  * W = 0.6098637*S - 0.6896511*j*w*D
 383  * X = 0.8624776*S + 0.7626955*j*w*D
 384  * Y = 1.6822415*w*D - 0.2156194*j*S
 385  *
 386  * where j is a +90 degree phase shift. w is a variable control for the
 387  * resulting stereo width, with the range 0 <= w <= 0.7.
 388  */
 389 template<size_t N>
 390 void UhjStereoDecoder<N>::decode(const al::span<float*> samples, const size_t samplesToDo,
 391     const bool updateState)
 392 {
 393     static_assert(sInputPadding <= sMaxPadding, "Filter padding is too large");
 394
 395     const auto &PShift = GetPhaseShifter<N>::Get();
 396
 397     ASSUME(samplesToDo > 0);
 398
 399     {
 400         const float *RESTRICT left{al::assume_aligned<16>(samples[0])};
 401         const float *RESTRICT right{al::assume_aligned<16>(samples[1])};
 402
 403         for(size_t i{0};i < samplesToDo+sInputPadding;++i)
 404             mS[i] = left[i] + right[i];
 405
 406         /* Pre-apply the width factor to the difference signal D. Smoothly
 407          * interpolate when it changes.
 408          */
 409         const float wtarget{mWidthControl};
 410         const float wcurrent{(mCurrentWidth < 0.0f) ? wtarget : mCurrentWidth};
 411         if(wtarget == wcurrent || !updateState)
 412         {
 413             for(size_t i{0};i < samplesToDo+sInputPadding;++i)
 414                 mD[i] = (left[i] - right[i]) * wcurrent;
 415             mCurrentWidth = wcurrent;
 416         }
 417         else
 418         {
 419             const float wstep{(wtarget - wcurrent) / static_cast<float>(samplesToDo)};
 420             float fi{0.0f};
 421             for(size_t i{0};i < samplesToDo;++i)
 422             {
 423                 mD[i] = (left[i] - right[i]) * (wcurrent + wstep*fi);
 424                 fi += 1.0f;
 425             }
 426             for(size_t i{samplesToDo};i < samplesToDo+sInputPadding;++i)
 427                 mD[i] = (left[i] - right[i]) * wtarget;
 428             mCurrentWidth = wtarget;
 429         }
 430     }
 431
 432     float *RESTRICT woutput{al::assume_aligned<16>(samples[0])};
 433     float *RESTRICT xoutput{al::assume_aligned<16>(samples[1])};
 434     float *RESTRICT youtput{al::assume_aligned<16>(samples[2])};
 435
 436     /* Precompute j*D and store in xoutput. */
 437     auto tmpiter = std::copy(mDTHistory.cbegin(), mDTHistory.cend(), mTemp.begin());
 438     std::copy_n(mD.cbegin(), samplesToDo+sInputPadding, tmpiter);
 439     if(updateState) LIKELY
 440         std::copy_n(mTemp.cbegin()+samplesToDo, mDTHistory.size(), mDTHistory.begin());
 441     PShift.process({xoutput, samplesToDo}, mTemp.data());
 442
 443     /* W = 0.6098637*S - 0.6896511*j*w*D */
 444     for(size_t i{0};i < samplesToDo;++i)
 445         woutput[i] = 0.6098637f*mS[i] - 0.6896511f*xoutput[i];
 446     /* X = 0.8624776*S + 0.7626955*j*w*D */
 447     for(size_t i{0};i < samplesToDo;++i)
 448         xoutput[i] = 0.8624776f*mS[i] + 0.7626955f*xoutput[i];
 449
 450     /* Precompute j*S and store in youtput. */
 451     tmpiter = std::copy(mSHistory.cbegin(), mSHistory.cend(), mTemp.begin());
 452     std::copy_n(mS.cbegin(), samplesToDo+sInputPadding, tmpiter);
 453     if(updateState) LIKELY
 454         std::copy_n(mTemp.cbegin()+samplesToDo, mSHistory.size(), mSHistory.begin());
 455     PShift.process({youtput, samplesToDo}, mTemp.data());
 456
 457     /* Y = 1.6822415*w*D - 0.2156194*j*S */
 458     for(size_t i{0};i < samplesToDo;++i)
 459         youtput[i] = 1.6822415f*mD[i] - 0.2156194f*youtput[i];
 460 }
 461
 462 void UhjStereoDecoderIIR::decode(const al::span<float*> samples, const size_t samplesToDo,
 463     const bool updateState)
 464 {
 465     static_assert(sInputPadding <= sMaxPadding, "Filter padding is too large");
 466
 467     ASSUME(samplesToDo > 0);
 468
 469     {
 470         const float *RESTRICT left{al::assume_aligned<16>(samples[0])};
 471         const float *RESTRICT right{al::assume_aligned<16>(samples[1])};
 472
 473         for(size_t i{0};i < samplesToDo;++i)
 474             mS[i] = left[i] + right[i];
 475
 476         /* Pre-apply the width factor to the difference signal D. Smoothly
 477          * interpolate when it changes.
 478          */
 479         const float wtarget{mWidthControl};
 480         const float wcurrent{(mCurrentWidth < 0.0f) ? wtarget : mCurrentWidth};
 481         if(wtarget == wcurrent || !updateState)
 482         {
 483             for(size_t i{0};i < samplesToDo;++i)
 484                 mD[i] = (left[i] - right[i]) * wcurrent;
 485             mCurrentWidth = wcurrent;
 486         }
 487         else
 488         {
 489             const float wstep{(wtarget - wcurrent) / static_cast<float>(samplesToDo)};
 490             float fi{0.0f};
 491             for(size_t i{0};i < samplesToDo;++i)
 492             {
 493                 mD[i] = (left[i] - right[i]) * (wcurrent + wstep*fi);
 494                 fi += 1.0f;
 495             }
 496             mCurrentWidth = wtarget;
 497         }
 498     }
 499
 500     float *RESTRICT woutput{al::assume_aligned<16>(samples[0])};
 501     float *RESTRICT xoutput{al::assume_aligned<16>(samples[1])};
 502     float *RESTRICT youtput{al::assume_aligned<16>(samples[2])};
 503
 504     /* Apply filter1 to S and store in mTemp. */
 505     mTemp[0] = mDelayS;
 506     mFilter1S.process(Filter1Coeff, {mS.data(), samplesToDo}, updateState, mTemp.data()+1);
 507     if(updateState) LIKELY mDelayS = mTemp[samplesToDo];
 508
 509     /* Precompute j*D and store in xoutput. */
 510     mFilter2D.process(Filter2Coeff, {mD.data(), samplesToDo}, updateState, xoutput);
 511
 512     /* W = 0.6098637*S - 0.6896511*j*w*D */
 513     for(size_t i{0};i < samplesToDo;++i)
 514         woutput[i] = 0.6098637f*mTemp[i] - 0.6896511f*xoutput[i];
 515     /* X = 0.8624776*S + 0.7626955*j*w*D */
 516     for(size_t i{0};i < samplesToDo;++i)
 517         xoutput[i] = 0.8624776f*mTemp[i] + 0.7626955f*xoutput[i];
 518
 519     /* Precompute j*S and store in youtput. */
 520     mFilter2S.process(Filter2Coeff, {mS.data(), samplesToDo}, updateState, youtput);
 521
 522     /* Apply filter1 to D and store in mTemp. */
 523     mTemp[0] = mDelayD;
 524     mFilter1D.process(Filter1Coeff, {mD.data(), samplesToDo}, updateState, mTemp.data()+1);
 525     if(updateState) LIKELY mDelayD = mTemp[samplesToDo];
 526
 527     /* Y = 1.6822415*w*D - 0.2156194*j*S */
 528     for(size_t i{0};i < samplesToDo;++i)
 529         youtput[i] = 1.6822415f*mTemp[i] - 0.2156194f*youtput[i];
 530 }
 531
 532
 533 template struct UhjEncoder<UhjLength256>;
 534 template struct UhjDecoder<UhjLength256>;
 535 template struct UhjStereoDecoder<UhjLength256>;
 536
 537 template struct UhjEncoder<UhjLength512>;
 538 template struct UhjDecoder<UhjLength512>;
 539 template struct UhjStereoDecoder<UhjLength512>;