Compute the tail time for an IIRFilter from its coefficients

third_party/WebKit/Source/platform/audio/IIRFilter.cpp

 // Copyright 2016 The Chromium Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 #include "platform/audio/IIRFilter.h"
 
 #include <complex>
+#include "platform/audio/AudioUtilities.h"
+#include "platform/audio/VectorMath.h"
 #include "platform/wtf/MathExtras.h"
 
 namespace blink {
 
 // The length of the memory buffers for the IIR filter.  This MUST be a power of
 // two and must be greater than the possible length of the filter coefficients.
 const int kBufferLength = 32;
 static_assert(kBufferLength >= IIRFilter::kMaxOrder + 1,
               "Internal IIR buffer length must be greater than maximum IIR "
               "Filter order.");
 
 IIRFilter::IIRFilter(const AudioDoubleArray* feedforward,
                      const AudioDoubleArray* feedback)
     : buffer_index_(0), feedback_(feedback), feedforward_(feedforward) {
   // These are guaranteed to be zero-initialized.
   x_buffer_.Allocate(kBufferLength);
   y_buffer_.Allocate(kBufferLength);
 }
 
 IIRFilter::~IIRFilter() {}
 
 void IIRFilter::Reset() {
   x_buffer_.Zero();
   y_buffer_.Zero();
+  buffer_index_ = 0;
 }
 
 static std::complex<double> EvaluatePolynomial(const double* coef,
                                                std::complex<double> z,
                                                int order) {
   // Use Horner's method to evaluate the polynomial P(z) = sum(coef[k]*z^k, k,
   // 0, order);
   std::complex<double> result = 0;
 
   for (int k = order; k >= 0; --k)
@@ skipping 90 matching lines @@
         feedforward_->Data(), z_recip, feedforward_->size() - 1);
     std::complex<double> denominator =
         EvaluatePolynomial(feedback_->Data(), z_recip, feedback_->size() - 1);
     std::complex<double> response = numerator / denominator;
     mag_response[k] = static_cast<float>(abs(response));
     phase_response[k] =
         static_cast<float>(atan2(imag(response), real(response)));
   }
 }
 
+double IIRFilter::TailTime(double sample_rate) {
+  // The maximum tail time.  This is somewhat arbitrary, but we're assuming that
+  // no one is going to expect the IIRFilter to produce an output after this
+  // much time after the inputs have stopped.
+  const double kMaxTailTime = 60;

[hongchan, 2017/05/12 16:24:05]

Let's pull these constants out of this function - probably move them to the top
of this file.

[Raymond Toy, 2017/05/12 19:42:51]

But they're not relevant to any other part of the code.

+
+  // If the maximum amplitude of the impulse response is less than this, we
+  // assume that we've reached the tail of the response.  Currently, this means
+  // that the impulse is less than 1 bit of a 16-bit PCM value.
+  const float kMaxTailAmplitude = 1 / 32768.0;

[hongchan, 2017/05/12 16:24:05]

Ditto.

+
+  // How to compute the tail time?  We're going to filter an impulse
+  // for |kMaxTailTime| seconds, in blocks of kRenderQuantumFrames at
+  // a time.  The maximum magnitude of this block is saved.  After all
+  // of the samples have been computed, find the last block with a
+  // maximum magnitude greater than |kMaxTaileAmplitude|.  That block
+  // index + 1 will be the tail time.  We don't need to be
+  // super-accurate in computing the tail time since we process on
+  // blocks, block accuracy is good enough, and the value just needs
+  // to be larger than the "real" tail time, so we don't prematurely
+  // zero out the output of the node.
+
+  // Number of render quanta needed to reach the max tail time.
+  int number_of_blocks = std::ceil(sample_rate * kMaxTailTime /
+                                   AudioUtilities::kRenderQuantumFrames);
+
+  // Input and output buffers for filtering.
+  AudioFloatArray input(AudioUtilities::kRenderQuantumFrames);
+  AudioFloatArray output(AudioUtilities::kRenderQuantumFrames);
+
+  // Array to hold the max magnitudes
+  AudioFloatArray magnitudes(number_of_blocks);
+
+  // Create the impulse input signal.
+  input[0] = 1;
+
+  // Process the first block and get the max magnitude of the output.
+  Process(input.Data(), output.Data(), AudioUtilities::kRenderQuantumFrames);
+  VectorMath::Vmaxmgv(output.Data(), 1, &magnitudes[0],
+                      AudioUtilities::kRenderQuantumFrames);
+
+  // Process the rest of the signal, getting the max magnitude of the
+  // output for each block.
+  input[0] = 0;
+
+  for (int k = 1; k < number_of_blocks; ++k) {
+    Process(input.Data(), output.Data(), AudioUtilities::kRenderQuantumFrames);
+    VectorMath::Vmaxmgv(output.Data(), 1, &magnitudes[k],
+                        AudioUtilities::kRenderQuantumFrames);
+  }
+
+  // Done computing the impulse response; reset the state so the actual node
+  // starts in the expected initial state.
+  Reset();
+
+  // Find the last block with amplitude greater than the threshold.
+  int index = number_of_blocks - 1;
+  for (int k = number_of_blocks - 1; k >= 0; --k) {

[hongchan, 2017/05/12 16:24:04]

|number_of_blocks -1| can be |index| here.

[Raymond Toy, 2017/05/12 19:42:51]

Done.

+    if (magnitudes[k] > kMaxTailAmplitude) {
+      index = k;
+      break;
+    }
+  }
+
+  // The magnitude first become lower than the threshold at the next block.
+  // Compute the corresponding time value value; that's the tail time.
+  double time =
+      (index + 1) * AudioUtilities::kRenderQuantumFrames / sample_rate;
+
+  return time;

[hongchan, 2017/05/12 16:24:05]

We don't need |time| here. Right?

[Raymond Toy, 2017/05/12 19:42:51]

Done.

+}
+
 }  // namespace blink