<!DOCTYPE html>
<html>
  <head>
    <title>
      Test Tail Time for IIRFilter
    </title>
    <script src="../../imported/w3c/web-platform-tests/resources/testharness.js"></script>
    <script src="../../resources/testharnessreport.js"></script>
    <script src="../resources/audit-util.js"></script>
    <script src="../resources/audit.js"></script>
  </head>
  <body>
    <script id="layout-test-code">
      let audit = Audit.createTaskRunner();
      let renderQuantumFrames = 128;

      // Must be a power of two to eliminate round-off differences between thsi
      // JS code and the WebAudio implementation.  Otherwise, the sample rate is
      // arbitrary.
      let sampleRate = 16384;

      // Fairly arbitrary, but should be long enough so that the node propagates
      // silence before the end of the offline context.
      let renderDuration = 1;


      audit.define('1-pole tail', (task, should) => {
        let pole = 0.99;
        let IIROptions = {feedforward: [1], feedback: [1, -pole]};
        // For the given filter, we can actually compute where the tail
        // begins.  The impulse response for the 1-pole filter is h(n) =
        // a^n, where a = 0.9.  The tail here starts when a^n < eps =
        // 1/32768.  So n > log(eps)/log(a), or 98.7.  Round that up to the
        // nearest render quantum frames.
        let tail = Math.ceil(Math.log(1 / 32768) / Math.log(pole));

        runTest(should, IIROptions, tail, '1-pole').then(() => task.done());
      });

      audit.define('2 real pole test', (task, should) => {
        // Simple example of a 2-pole IIR filter where both poles are real.
        // We arbitrarily select a pole at 9.99 and one at -0.5. The IIRFilter
        // is then
        //   1 / ((z-0.99) * (z + 0.5))
        //     = 1/(z^2-0.49z-0.495)
        //     = z^-2/(1-0.49/z-0.495/z^2)
        let IIROptions = {feedforward: [0, 0, 1], feedback: [1, -0.49, -0.495]};

        // For this particular filter, we can analytically compute the impulse
        // response using partical fractios:
        //
        //   1 / ((z-0.99) * (z + 0.5))
        //   = 1/(-0.5-0.99)/(z + 0.5) - 1/(-0.5-0.99)/(z - 0.99)
        //   = 1/1.49*(1/(z-0.99) - 1/(z+0.5))
        //   = 1/1.49*[1/z*sum(.99^n/z^n,n,0,inf)
        //       - 1/z*sum((-0.5)^n/z^n,n,0,inf)]
        //   = 1/1.49/z*sum((0.99^n-(-0.5)^n)/z^n)
        //
        // So the tail begins when 1/1.49*(0.99^n-(-0.5)^n) < 1/32768.  This can
        // be solved numerically to give n = 995.
        let tail = 995;
        tail = renderQuantumFrames * Math.ceil(tail / renderQuantumFrames);

        runTest(should, IIROptions, tail, '2 real poles')
            .then(() => task.done());
      });

      audit.define('2 complex poles', (task, should) => {
        // Simple example of a 2-pole IIR filter where both poles are complex
        // conjugates.  In this case, the poles will be r*exp(+/-i*theta) where
        // r = 0.99 and theta = 0.01.  The filter is then
        //
        //  1/(z^2-2*r*cos(theta) + r^2)
        //    = z^(-2)/(1-2*r*cos(theta)/z + r^2/z^2)
        let r = 0.99;
        let theta = 0.01;
        let IIROptions = {
          feedforward: [0, 0, 1],
          feedback: [1, -2 * r * Math.cos(theta), r * r]
        };

        // Again, we can use partial fractions as for 2 real pole case to get an
        // analytically solution for the impulse response.  For simplicity, let
        // p1 = r*exp(i*theta), p2 = r*exp(-i*theta).  Then:
        //
        //  1/(z^2-2*r*cos(theta) + r^2)
        //    = 1/(z-p1)/(z-p2)
        //    = 1/(p2-p1)*[1/(z-p2) - 1/(z-p1)]
        //    = 1/(p2-p1)*[1/z*sum(p2^n/z^n) - 1/z*sum(p1^n/z^n)]
        //    = 1/(p2-p1)/z*sum((p2^n-p1^n)/z^n)
        //
        // So the tail begins when
        //   1/32768 > |1/(p2-p1)*(p2^n-p1^n)|
        //           = 1/(r*sin(theta))*|r^n*(exp(-i*theta*n)-exp(i*theta*n))|
        //           = 1/(2*r*sin(theta))*(2*r^n*|sin(theta*n)|);
        //           = r^(n-1)*|sin(theta*n)|/sin(theta)
        //
        // This can be solved numerically to for n;
        let tail = 1474.256;
        tail = renderQuantumFrames * Math.ceil(tail / renderQuantumFrames);

        runTest(should, IIROptions, tail, '2 complex poles')
            .then(() => task.done());
      });

      audit.define('repeated poles', (task, should) => {
        // Two repeated roots.  Let p be the repeated pole.  Then the filter is
        //
        //    1/(z-p)^2
        //      = z^(-2)/(1-p/z)^2
        //      = z^(-2)/(1-2*p/z+p*p/z^2)

        let pole = 0.99;
        let IIROptions = {
          feedforward: [0, 0, 1],
          feedback: [1, -2 * pole, pole * pole]
        };

        // We can analytically compute the impulse response of this filter to be
        //
        //   1/z^2*sum(p^n*(n+1)/z^n, n, 0, inf)
        //     = sum(p^n*(n+1)/z^(n+2), n, 0, inf)
        //     = 1/p^2*sum((p^k*(k-1))/z^k,k,2,inf))
        //
        // Therefore the tail starts when p^(k-2)*(k-1) < 1/32768.  We can solve
        // this numerically to be 1781.213;

        let tail = 1781.213;
        runTest(should, IIROptions, tail, '2 repeated poles')
            .then(() => task.done());

      });

      audit.define('4-th order', (task, should) => {
        // Test consistency of tail times between a 4-th order direct IIR filter
        // and the equivalent cascade of second-order sections.  The first
        // channel of the output is the cascaded biquad, and the second channel
        // is the 4-th order equivalent.
        let context =
            new OfflineAudioContext(2, renderDuration * sampleRate, sampleRate);

        let src = new AudioBufferSourceNode(
            context, {buffer: createImpulseBuffer(context, 1)});

        // This is a 4-th order lowpass elliptic filter designed using
        // http://rtoy.github.io/webaudio-hacks/more/filter-design/filter-design.html.
        // The sample rate is 16384 Hz with a passband at 3600 Hz with a 0.25 dB
        // attenuation, and a stopband at 4800 Hz, with a stopband attenuation
        // of 30 dB. (Nothing really special except that this gives a 4-th order
        // filter).

        let f0 = context.createIIRFilter(
            [0.6410686464424084, 0.2607836369670137, 0.6410686464424084],
            [1, -0.2287413068432929, 0.7716622366951231]);
        let f1 = context.createIIRFilter(
            [0.21283904239866536, 0.3184888523034876, 0.21283904239866536],
            [1, -0.4686913542990081, 0.21285829139982618]);

        // The poles for f0 are 0.1143706534216465 +/- 0.8709658950447078*i or
        // 0.8784430753868592*%e^(+/-1.440228658066206*%i),
        //
        // The poles for f1 are 0.2343456771495041 +/- 0.3974171548903829*i or
        // 0.4613656807780854*%e^(+/-1.038005727602151*%i.
        //
        // Thus, the tail time for f0 is approximately 80, but this is an
        // approximation since we didn't include the affect of the numerator.
        // Round this up to the next render to get an actual tail time of 128.
        //
        // Similarly, for f0, the tail time is 14.3. Thus, the actual tail time
        // is alos 128 for this filter.
        //
        // Since these biquads are cascaded, the total tail time for both is the
        // sum or 256 frames.  However, the tail actually ends two render quanta
        // after this for a total of 512 frames.

        let biquadTailEnd = 512;

        // The equivalent 4-th order filter created multiplying the f0 and f1
        // coefficients together appropriately.
        let f = context.createIIRFilter(
            [
              0.136444436820611, 0.259678157018493, 0.355945554878375,
              0.259678157018493, 0.136444436820611
            ],
            [
              1.000000000000000, -0.697432661142301, 1.091729600983457,
              -0.410360902525266, 0.164254705240692
            ]);

        let merger = context.createChannelMerger(2);
        merger.connect(context.destination);

        src.connect(f0).connect(f1).connect(merger, 0, 0);
        src.connect(f).connect(merger, 0, 1);

        src.start();

        context.startRendering()
            .then(renderedBuffer => {
              // c0 = cascaded biquads
              // c1 = 4-th order filter
              let c0 = renderedBuffer.getChannelData(0);
              let c1 = renderedBuffer.getChannelData(1);

              // Sanity check: The two filters should have the same output
              // within some rounding error.
              should(
                  c0.slice(0, biquadTailEnd),
                  'Filter outputs[0:' + (biquadTailEnd - 1) + ']')
                  .beCloseToArray(
                      c1.slice(0, biquadTailEnd),
                      {absoluteThreshold: 1.4902e-8});
              should(
                  c0.slice(biquadTailEnd),
                  'Filter outputs[' + biquadTailEnd + ':]')
                  .beEqualToArray(c1.slice(biquadTailEnd));

              // Verify that after the tail time, the outputs are zero, and not
              // before for both the biquads and 4-th order filters.
              should(
                  c0.slice(0, biquadTailEnd),
                  'cascaded biquad output[0:' + (biquadTailEnd - 1) + ']')
                  .notBeConstantValueOf(0);
              should(
                  c0.slice(biquadTailEnd),
                  'cascaded biquad output[' + biquadTailEnd + ':]')
                  .beConstantValueOf(0);

              should(
                  c1.slice(0, biquadTailEnd),
                  '4-th order output[0:' + (biquadTailEnd - 1) + ']')
                  .notBeConstantValueOf(0);
              should(
                  c1.slice(biquadTailEnd),
                  '4-th order output[' + biquadTailEnd + ':]')
                  .beConstantValueOf(0);
            })
            .then(() => task.done());
      });

      function runTest(should, IIROptions, tailFrames, prefix) {
        let context =
            new OfflineAudioContext(1, renderDuration * sampleRate, sampleRate);

        let src = new AudioBufferSourceNode(
            context, {buffer: createImpulseBuffer(context, 1)});

        let iir = new IIRFilterNode(context, IIROptions);

        src.connect(iir).connect(context.destination);

        src.start();

        return context.startRendering().then(renderedBuffer => {
          let audio = renderedBuffer.getChannelData(0);

          // Round up the tailFrames to the nearest render quantum.
          let tailQuantum =
              renderQuantumFrames * Math.ceil(tailFrames / renderQuantumFrames);
          let tailEndFrame = tailQuantum + 2 * renderQuantumFrames;

          should(tailEndFrame, prefix + ': tail end frame')
              .beLessThanOrEqualTo(context.length);

          // Clamp to the render duration so we don't go off the end.
          tailEndFrame = Math.min(tailEndFrame, context.length);

          for (let k = 0; k < tailEndFrame; k += renderQuantumFrames) {
            should(
                audio.slice(k, k + renderQuantumFrames),
                prefix + ': output[' + k + ':' + (k + renderQuantumFrames - 1) +
                    ']')
                .notBeConstantValueOf(0);
          }

          if (tailEndFrame < context.length) {
            // All frames after should be zero because we're propagating
            // silence.
            should(
                audio.slice(tailEndFrame),
                'output[' + tailEndFrame + ':' + (context.length - 1) + ']')
                .beConstantValueOf(0);
          }
        });
      }

      audit.run();
    </script>
  </body>
</html>