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