GCC Code Coverage Report

Directory: src/
File: src/LowPassFilter.cpp
Date: 2026-09-16 13:46:10
Exec Total Coverage
Lines: 106 106 100.0%
Branches: 70 106 66.0%
Line Branch Exec Source
1 /*****************************************************************************
2 *
3 * Copyright (C) 2021 - 2022 Jonathan Grahl <jonathan.grahl@igh.de>
4 * 2021 - 2026 Florian Pose <florian.pose@igh.de>
5 *
6 * This file is part of the reta library (realtime-automation).
7 *
8 * The reta library is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Lesser General Public License as published
10 * by the Free Software Foundation, either version 3 of the License, or (at
11 * your option) any later version.
12 *
13 * The reta library is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
16 * General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public License
19 * along with the reta library. If not, see <http://www.gnu.org/licenses/>.
20 *
21 ****************************************************************************/
22
23 #include "reta/LowPassFilter.h"
24
25 #include "Base.h"
26
27 #include <cmath> // cos(), sin()
28 #include <sstream>
29 #include <stdexcept>
30
31 using std::function;
32 using std::invalid_argument;
33 using std::make_unique;
34 using std::shared_ptr;
35 using std::string;
36 using std::stringstream;
37 using std::vector;
38
39 using namespace Reta;
40
41 /****************************************************************************/
42
43 39 struct RETA_LOCAL LowPassFilter::Impl : public Base
44 {
45 // Direct Form II Transposed biquad section (RBJ audio cookbook). For
46 // order == 1, a1/b0 are set up as a plain single-pole filter and b1/b2/a2
47 // stay 0, so the same process() reduces to a first-order filter without
48 // any special-casing in update().
49 62 struct Biquad
50 {
51 double b0 {1.0}, b1 {0.0}, b2 {0.0}, a1 {0.0}, a2 {0.0};
52 double z1 {0.0}, z2 {0.0};
53
54 322549 double process(double in)
55 {
56 322549 double out = b0 * in + z1;
57 322549 z1 = b1 * in - a1 * out + z2;
58 322549 z2 = b2 * in - a2 * out;
59 322549 return out;
60 }
61 };
62
63 Impl(shared_ptr<Task>,
64 const string &,
65 double,
66 unsigned int,
67 unsigned int);
68
69 void
70 update(const vector<double> &,
71 const vector<uint8_t> & = vector<uint8_t>());
72 void updateCoefficients(unsigned int, double);
73
74 shared_ptr<Task> task;
75
76 // order == 1: single real pole, one section
77 // Even order >= 2: (order / 2) RBJ Butterworth biquad sections
78 unsigned int order;
79
80 vector<double> tau;
81 vector<double> output;
82
83 // Set (order >= 2 only) whenever tau had to be clamped to stay clear
84 // of the Nyquist frequency, see updateCoefficients().
85 vector<uint8_t> tauLimited;
86
87 // tau the cascade coefficients were last computed for, per channel, so
88 // that the (comparatively expensive) trigonometry in
89 // updateCoefficients() only runs when tau actually changes.
90 vector<double> stageTau;
91 vector<vector<Biquad>> stages;
92 };
93
94 /****************************************************************************/
95
96 41 LowPassFilter::Impl::Impl(
97 shared_ptr<Task> task,
98 const string &prefix,
99 double initTau,
100 unsigned int width,
101 41 unsigned int filterOrder) :
102 Base {prefix},
103 task {task},
104 order {filterOrder},
105 tau(width, initTau),
106 output(width, 0.0),
107 tauLimited(width, 0),
108 stageTau(width, 0.0),
109
5/10
✓ Branch 4 taken 41 times.
✗ Branch 5 not taken.
✓ Branch 9 taken 41 times.
✗ Branch 10 not taken.
✓ Branch 14 taken 41 times.
✗ Branch 15 not taken.
✓ Branch 19 taken 41 times.
✗ Branch 20 not taken.
✓ Branch 24 taken 41 times.
✗ Branch 25 not taken.
 43 stages(width)
110 {
111
2/2
✓ Branch 1 taken 40 times.
✓ Branch 2 taken 1 times.
 41 checkZeroWidth(width);
112
113
5/6
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 35 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 4 times.
 40 if (order == 0 || (order > 1 && order % 2 != 0)) {
114
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 stringstream str;
115 str << __PRETTY_FUNCTION__ << " : " << prefix
116
5/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
 1 << " Invalid filter order " << order
117
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 << " (must be 1, or an even number >= 2)!";
118
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 throw invalid_argument(str.str());
119 }
120
121
2/2
✓ Branch 3 taken 62 times.
✓ Branch 4 taken 39 times.
 101 for (auto &s : stages) {
122
3/4
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 58 times.
✓ Branch 3 taken 62 times.
✗ Branch 4 not taken.
 62 s.resize(order == 1 ? 1 : order / 2);
123 }
124
125
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 39 pdserv *pdserv {task->getPdServ()};
126
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 39 pdtask *pdtask {task->getPdTask()};
127
128
1/2
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
 78 pdserv_parameter(
129
1/2
✓ Branch 2 taken 39 times.
✗ Branch 3 not taken.
 78 pdserv, (prefix + "/Tau").c_str(), 0666, pd_double_T, tau.data(),
130 tau.size(), NULL, NULL, NULL);
131
1/2
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
 78 pdserv_signal(
132
1/2
✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
 78 pdtask, 1, (prefix + "/Output").c_str(), pd_double_T,
133 39 output.data(), output.size(), NULL);
134
1/2
✓ Branch 3 taken 39 times.
✗ Branch 4 not taken.
 78 pdserv_signal(
135
1/2
✓ Branch 1 taken 39 times.
✗ Branch 2 not taken.
 78 pdtask, 1, (prefix + "/TauLimited").c_str(), pd_boolean_T,
136 39 tauLimited.data(), tauLimited.size(), NULL);
137 39 }
138
139 /****************************************************************************/
140
141 13 void LowPassFilter::Impl::updateCoefficients(unsigned int i, double t)
142 {
143 13 double period {task->getPeriod()};
144
145
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 6 times.
 13 if (order == 1) {
146 // Discretized single-pole filter, equivalent to Euler-integrating
147 // dy/dt = (x - y) / tau:
148 // y[n] = y[n-1] + (period / tau) * (x[n] - y[n-1]).
149 7 Biquad &s {stages.at(i).at(0)};
150 7 s.b0 = period / t;
151 7 s.b1 = 0.0;
152 7 s.b2 = 0.0;
153 7 s.a1 = period / t - 1.0;
154 7 s.a2 = 0.0;
155
156 7 stageTau.at(i) = t;
157 7 return;
158 }
159
160 6 double sampleRate {1.0 / period};
161 6 double nyquist {sampleRate / 2.0};
162 6 double cutoffHz {1.0 / (2.0 * M_PI * t)};
163 6 tauLimited.at(i) = cutoffHz > nyquist * 0.99;
164
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 5 times.
 6 if (tauLimited.at(i)) {
165 // stay clear of the Nyquist limit, where cos(w0) approaches -1
166 1 cutoffHz = nyquist * 0.99;
167 }
168
169 6 double w0 {2.0 * M_PI * cutoffHz / sampleRate};
170 6 double cosw0 {cos(w0)};
171 6 double sinw0 {sin(w0)};
172
173 6 unsigned int nSections {order / 2};
174
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 6 times.
 12 for (unsigned int k = 1; k <= nSections; k++) {
175 // Per-section Q from the Butterworth pole angles, so that the
176 // cascade of these biquads reproduces the N-th order Butterworth
177 // response (maximally flat passband).
178 6 double theta {M_PI * (2.0 * k - 1.0) / (2.0 * order)};
179 6 double Q {1.0 / (2.0 * cos(theta))};
180 6 double alpha {sinw0 / (2.0 * Q)};
181
182 // Standard bilinear-transform lowpass biquad ("RBJ cookbook").
183 6 double a0 {1.0 + alpha};
184 6 Biquad &s {stages.at(i).at(k - 1)};
185 6 s.b0 = ((1.0 - cosw0) / 2.0) / a0;
186 6 s.b1 = (1.0 - cosw0) / a0;
187 6 s.b2 = s.b0;
188 6 s.a1 = (-2.0 * cosw0) / a0;
189 6 s.a2 = (1.0 - alpha) / a0;
190 }
191
192 6 stageTau.at(i) = t;
193 }
194
195 /****************************************************************************/
196
197 321566 void LowPassFilter::Impl::update(
198 const vector<double> &input,
199 const vector<uint8_t> &reset)
200 {
201
2/2
✓ Branch 1 taken 323569 times.
✓ Branch 2 taken 321566 times.
 645135 for (unsigned int i = 0; i < output.size(); i++) {
202
6/6
✓ Branch 1 taken 323549 times.
✓ Branch 2 taken 20 times.
✓ Branch 4 taken 322549 times.
✓ Branch 5 taken 1000 times.
✓ Branch 6 taken 322549 times.
✓ Branch 7 taken 1020 times.
 323569 if (tau.at(i) > 0.0 && !reset.at(i)) {
203
2/2
✓ Branch 2 taken 13 times.
✓ Branch 3 taken 322536 times.
 322549 if (tau.at(i) != stageTau.at(i)) {
204 13 updateCoefficients(i, tau.at(i));
205 }
206
207 322549 double out {input.at(i)};
208
3/4
✓ Branch 1 taken 322549 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 322549 times.
✓ Branch 7 taken 322549 times.
 645098 for (auto &s : stages.at(i)) {
209 322549 out = s.process(out);
210 }
211 322549 output.at(i) = out;
212 }
213 else {
214 1020 output.at(i) = input.at(i);
215
216 // Drop the filter memory so that it starts fresh from the
217 // current input once tau becomes positive again.
218
3/4
✓ Branch 1 taken 1020 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 1020 times.
✓ Branch 7 taken 1020 times.
 2040 for (auto &s : stages.at(i)) {
219 1020 s.z1 = s.z2 = 0.0;
220 }
221 }
222 }
223 321566 }
224
225 /****************************************************************************/
226
227 29 LowPassFilter::LowPassFilter(
228 shared_ptr<Task> task,
229 const std::string &prefix,
230 double initTau,
231 29 unsigned int width) :
232
2/2
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 1 times.
 30 LowPassFilter {task, prefix, initTau, width, 1}
233 28 {}
234
235 /****************************************************************************/
236
237 41 LowPassFilter::LowPassFilter(
238 shared_ptr<Task> task,
239 const std::string &prefix,
240 double initTau,
241 unsigned int width,
242 41 unsigned int order) :
243 41 impl {make_unique<Impl>(task, prefix, initTau, width, order)}
244 39 {}
245
246 /****************************************************************************/
247
248 39 LowPassFilter::~LowPassFilter()
249 39 {}
250
251 /****************************************************************************/
252
253 2 void LowPassFilter::setTau(double newTau, unsigned int i)
254 {
255 2 impl->tau.at(i) = newTau;
256 2 }
257
258 /****************************************************************************/
259
260 320554 void LowPassFilter::update(double input, bool reset)
261 {
262
3/6
✓ Branch 3 taken 320554 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 320554 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 320554 times.
✗ Branch 11 not taken.
 320554 impl->update(vector<double> {input}, vector<uint8_t> {reset});
263 320554 }
264
265 /****************************************************************************/
266
267 12 void LowPassFilter::update(
268 const vector<double> &input,
269 const vector<uint8_t> &reset)
270 {
271 12 impl->update(input, reset);
272 12 }
273
274 /****************************************************************************/
275
276 1000 void LowPassFilter::update(
277 function<double(unsigned int)> input,
278 function<bool(unsigned int)> reset)
279 {
280
1/2
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
 2000 vector<double> inputVector(impl->output.size());
281
1/2
✓ Branch 4 taken 1000 times.
✗ Branch 5 not taken.
 2000 vector<uint8_t> resetVector(impl->output.size());
282
283
2/2
✓ Branch 2 taken 3000 times.
✓ Branch 3 taken 1000 times.
 4000 for (unsigned int i = 0; i < impl->output.size(); i++) {
284
2/4
✓ Branch 1 taken 3000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3000 times.
✗ Branch 5 not taken.
 3000 inputVector.at(i) = input(i);
285
2/4
✓ Branch 1 taken 3000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3000 times.
✗ Branch 5 not taken.
 3000 resetVector.at(i) = reset(i);
286 }
287
288
1/2
✓ Branch 2 taken 1000 times.
✗ Branch 3 not taken.
 1000 impl->update(inputVector, resetVector);
289 1000 }
290
291 /****************************************************************************/
292
293 316070 double LowPassFilter::getOutput(unsigned int i) const
294 {
295 316070 return impl->output.at(i);
296 }
297
298 /****************************************************************************/
299
300 4 bool LowPassFilter::isTauLimited(unsigned int i) const
301 {
302 4 return impl->tauLimited.at(i);
303 }
304
305 /****************************************************************************/
306