Line data Source code
1 : // Copyright 2016 The Chromium Authors. All rights reserved.
2 : // Use of this source code is governed by a BSD-style license that can be
3 : // found in the LICENSE file.
4 :
5 : #include "IIRFilter.h"
6 :
7 : #include "DenormalDisabler.h"
8 :
9 : #include <mozilla/Assertions.h>
10 :
11 : #include <complex>
12 :
13 : namespace blink {
14 :
15 : // The length of the memory buffers for the IIR filter. This MUST be a power of two and must be
16 : // greater than the possible length of the filter coefficients.
17 : const int kBufferLength = 32;
18 : static_assert(kBufferLength >= IIRFilter::kMaxOrder + 1,
19 : "Internal IIR buffer length must be greater than maximum IIR Filter order.");
20 :
21 0 : IIRFilter::IIRFilter(const AudioDoubleArray* feedforwardCoef, const AudioDoubleArray* feedbackCoef)
22 : : m_bufferIndex(0)
23 : , m_feedback(feedbackCoef)
24 0 : , m_feedforward(feedforwardCoef)
25 : {
26 0 : m_xBuffer.SetLength(kBufferLength);
27 0 : m_yBuffer.SetLength(kBufferLength);
28 0 : reset();
29 0 : }
30 :
31 0 : IIRFilter::~IIRFilter()
32 : {
33 0 : }
34 :
35 0 : void IIRFilter::reset()
36 : {
37 0 : memset(m_xBuffer.Elements(), 0, m_xBuffer.Length() * sizeof(double));
38 0 : memset(m_yBuffer.Elements(), 0, m_yBuffer.Length() * sizeof(double));
39 0 : }
40 :
41 0 : static std::complex<double> evaluatePolynomial(const double* coef, std::complex<double> z, int order)
42 : {
43 : // Use Horner's method to evaluate the polynomial P(z) = sum(coef[k]*z^k, k, 0, order);
44 0 : std::complex<double> result = 0;
45 :
46 0 : for (int k = order; k >= 0; --k)
47 0 : result = result * z + std::complex<double>(coef[k]);
48 :
49 0 : return result;
50 : }
51 :
52 0 : void IIRFilter::process(const float* sourceP, float* destP, size_t framesToProcess)
53 : {
54 : // Compute
55 : //
56 : // y[n] = sum(b[k] * x[n - k], k = 0, M) - sum(a[k] * y[n - k], k = 1, N)
57 : //
58 : // where b[k] are the feedforward coefficients and a[k] are the feedback coefficients of the
59 : // filter.
60 :
61 : // This is a Direct Form I implementation of an IIR Filter. Should we consider doing a
62 : // different implementation such as Transposed Direct Form II?
63 0 : const double* feedback = m_feedback->Elements();
64 0 : const double* feedforward = m_feedforward->Elements();
65 :
66 0 : MOZ_ASSERT(feedback);
67 0 : MOZ_ASSERT(feedforward);
68 :
69 : // Sanity check to see if the feedback coefficients have been scaled appropriately. It must
70 : // be EXACTLY 1!
71 0 : MOZ_ASSERT(feedback[0] == 1);
72 :
73 0 : int feedbackLength = m_feedback->Length();
74 0 : int feedforwardLength = m_feedforward->Length();
75 0 : int minLength = std::min(feedbackLength, feedforwardLength);
76 :
77 0 : double* xBuffer = m_xBuffer.Elements();
78 0 : double* yBuffer = m_yBuffer.Elements();
79 :
80 0 : for (size_t n = 0; n < framesToProcess; ++n) {
81 : // To help minimize roundoff, we compute using double's, even though the filter coefficients
82 : // only have single precision values.
83 0 : double yn = feedforward[0] * sourceP[n];
84 :
85 : // Run both the feedforward and feedback terms together, when possible.
86 0 : for (int k = 1; k < minLength; ++k) {
87 0 : int n = (m_bufferIndex - k) & (kBufferLength - 1);
88 0 : yn += feedforward[k] * xBuffer[n];
89 0 : yn -= feedback[k] * yBuffer[n];
90 : }
91 :
92 : // Handle any remaining feedforward or feedback terms.
93 0 : for (int k = minLength; k < feedforwardLength; ++k)
94 0 : yn += feedforward[k] * xBuffer[(m_bufferIndex - k) & (kBufferLength - 1)];
95 :
96 0 : for (int k = minLength; k < feedbackLength; ++k)
97 0 : yn -= feedback[k] * yBuffer[(m_bufferIndex - k) & (kBufferLength - 1)];
98 :
99 : // Save the current input and output values in the memory buffers for the next output.
100 0 : m_xBuffer[m_bufferIndex] = sourceP[n];
101 0 : m_yBuffer[m_bufferIndex] = yn;
102 :
103 0 : m_bufferIndex = (m_bufferIndex + 1) & (kBufferLength - 1);
104 :
105 : // Avoid introducing a stream of subnormals
106 0 : destP[n] = WebCore::DenormalDisabler::flushDenormalFloatToZero(yn);
107 0 : MOZ_ASSERT(destP[n] == 0.0 || fabs(destP[n]) > FLT_MIN, "output should not be subnormal");
108 : }
109 0 : }
110 :
111 0 : void IIRFilter::getFrequencyResponse(int nFrequencies, const float* frequency, float* magResponse, float* phaseResponse)
112 : {
113 : // Evaluate the z-transform of the filter at the given normalized frequencies from 0 to 1. (One
114 : // corresponds to the Nyquist frequency.)
115 : //
116 : // The z-tranform of the filter is
117 : //
118 : // H(z) = sum(b[k]*z^(-k), k, 0, M) / sum(a[k]*z^(-k), k, 0, N);
119 : //
120 : // The desired frequency response is H(exp(j*omega)), where omega is in [0, 1).
121 : //
122 : // Let P(x) = sum(c[k]*x^k, k, 0, P) be a polynomial of order P. Then each of the sums in H(z)
123 : // is equivalent to evaluating a polynomial at the point 1/z.
124 :
125 0 : for (int k = 0; k < nFrequencies; ++k) {
126 : // zRecip = 1/z = exp(-j*frequency)
127 0 : double omega = -M_PI * frequency[k];
128 0 : std::complex<double> zRecip = std::complex<double>(cos(omega), sin(omega));
129 :
130 0 : std::complex<double> numerator = evaluatePolynomial(m_feedforward->Elements(), zRecip, m_feedforward->Length() - 1);
131 0 : std::complex<double> denominator = evaluatePolynomial(m_feedback->Elements(), zRecip, m_feedback->Length() - 1);
132 : // Strangely enough, using complex division:
133 : // e.g. Complex response = numerator / denominator;
134 : // fails on our test machines, yielding infinities and NaNs, so we do
135 : // things the long way here.
136 0 : double n = norm(denominator);
137 0 : double r = (real(numerator)*real(denominator) + imag(numerator)*imag(denominator)) / n;
138 0 : double i = (imag(numerator)*real(denominator) - real(numerator)*imag(denominator)) / n;
139 0 : std::complex<double> response = std::complex<double>(r, i);
140 :
141 0 : magResponse[k] = static_cast<float>(abs(response));
142 0 : phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
143 : }
144 0 : }
145 :
146 0 : bool IIRFilter::buffersAreZero()
147 : {
148 0 : double* xBuffer = m_xBuffer.Elements();
149 0 : double* yBuffer = m_yBuffer.Elements();
150 :
151 0 : for (size_t k = 0; k < m_feedforward->Length(); ++k) {
152 0 : if (xBuffer[(m_bufferIndex - k) & (kBufferLength - 1)] != 0.0) {
153 0 : return false;
154 : }
155 : }
156 :
157 0 : for (size_t k = 0; k < m_feedback->Length(); ++k) {
158 0 : if (fabs(yBuffer[(m_bufferIndex - k) & (kBufferLength - 1)]) >= FLT_MIN) {
159 0 : return false;
160 : }
161 : }
162 :
163 0 : return true;
164 : }
165 :
166 : } // namespace blink
|
