A low-pass filter is a filter which removes high frequencies from a signal. Low-pass filters are very useful in digital audio applications. For example:
- Removing artifacts when resampling.
- Emulating the response of a speaker cabinet.
- Removing harsh sounding frequencies introduced by wave shaping effects such as overdrive.
In this post, we will be looking at how to design a digital low-pass filter.
This is a problem that was famously tackled by Stephen Butterworth. As he put it:
An ideal electrical filter should not only completely reject the unwanted frequencies but should also have uniform sensitivity for the wanted frequencies.
Butterworth discovered a formula for designing low-pass filters of arbitrary order. As we increase the order, we get progressively closer to Butterworth’s idea of an ideal filter. So a low order filter will suffer from slow frequency roll off. A higher order filter will have a faster roll off (at the expense of slightly higher computational expense).
Creating a Butterworth filter in Matlab
One way to design the filter is to use Matlab (or Octave). The command is as follows:
[b, a] = butter(N, Wc)
where N is the order of the IIR filter and Wc is the relative cutoff frequency. The value of Wc should be between 0 and 1, where 1 corresponds to the Nyquist rate, i.e. half the sample rate.
The resulting arrays a and b give the IIR coefficients required to implement the filter. We’ll give more details on these in the following sections.
The cutoff frequency is the frequency at which the magnitude response of the filter is .
For example, if the sample rate is 44.1kHz and we want the cutoff frequency to be 8kHz, we could generate a second-order filter as follows:
fc = 8000;
fs = 44100;
[b,a] = butter(2, fc/(fs/2))
This would give the following output:
b =
0.1772 0.3545 0.1772
a =
1.0000 -0.5087 0.2177
Note that (following convention) is always unity.
Implementing the IIR filter in code
Once we know the coefficients and , we can write a formula to process the samples of the input signal in order to calculate the output signal .
The Butterworth low-pass filter is a type of IIR filter. Butterworth filters have the same feedforward and feedback order (say ). In this case, the difference formula can be written as:
Note that we have used the convention that is unity.
Designing the filter in C++
We have seen how to design the filter in Matlab and we have established a formula that we can code in C++ to implement the filter.
That’s fine if we know in advance what the cutoff frequency is. But if the frequency can vary at runtime (e.g. we allow the user to select the frequency via the user interface) then we have a harder problem. Unlike Matlab, C++ doesn’t have a standard library for designing filters.
At this point, the following post on dsprelated.com proved to be very useful:
In this post, the author helpfully recreates Matlab’s built-in butter() function in Matlab code. So we just need to translate the sample Matlab code into C++ and we are done.
My C++ implementation is presented below. Whilst Matlab natively supports complex polynomials, in C++ we have a little more work to do. For complex numbers, we can use the standard library’s std::complex class. However, for polynomials, we need to recreate some of Matlab’s built-in functions, namely the poly() function and the multiplication operator.
ButterworthFilter.h
#pragma once
#include <vector>
class IIR_Coeffs
{
public:
std::vector<double> a;
std::vector<double> b;
};
IIR_Coeffs butter_synth(int N, double fc, double fs);
ButterworthFilter.cpp
#include "ButterworthFilter.h"
#include <complex>
#include <cmath>
#include <numeric>
#include <assert.h>
using namespace std::literals::complex_literals;
using Polynomial = std::vector<std::complex<double>>;
Polynomial operator*(const Polynomial& p, const Polynomial& q)
{
size_t n = p.size() + q.size() - 1;
Polynomial result(n);
for (size_t i = 0; i < p.size(); i++)
for (size_t j = 0; j < q.size(); j++)
result[i + j] += p[i] * q[j];
return result;
}
// Emulate Matlab 'poly' function. Calculate the coefficients of the polynomial
// with the specified roots.
static std::vector<std::complex<double>> poly(std::vector<std::complex<double>> roots)
{
Polynomial result{1.0};
for (auto root : roots)
{
Polynomial factor({-root, 1.0});
result = result * factor;
}
// Matlab returns the highest order coefficients first.
std::reverse(result.begin(), result.end());
return result;
}
// Emulate Matlab 'sum' function.
static std::complex<double> sum(const std::vector<std::complex<double>>& v)
{
return std::accumulate(v.begin(), v.end(), 0.0i);
}
// Taken from here : https://www.dsprelated.com/showarticle/1119.php
// The blog post gives Matlab source code which I have converted to C++ below:
IIR_Coeffs butter_synth(int N, double fc, double fs)
{
constexpr double pi = 3.14159265358979323846;
IIR_Coeffs coeffs;
std::vector<std::complex<double>> pa(N);
std::vector<std::complex<double>> p(N);
std::vector<std::complex<double>> q(N, -1.0);
assert(fc < fs / 2); // Cutoff frequency must be less that fs/2
// I. Find poles of analog filter
for (int i = 0; i < N; i++)
{
int k = i + 1;
double theta = (2 * k - 1) * pi / (2 * N);
pa[i] = -sin(theta) + 1.0i * cos(theta);
}
// II. Scale poles in frequency
double Fc = fs / pi * tan(pi * fc / fs);
for (size_t i=0; i<pa.size(); i++)
pa[i] *= 2 * pi * Fc;
// III. Find coeffs of digital filter poles and zeros in the z plane
for (size_t i = 0; i < N; i++)
p[i] = (1.0 + pa[i] / (2 * fs)) / (1.0 - pa[i] / (2 * fs));
auto a = poly(p);
for (size_t i = 0; i < a.size(); i++)
a[i] = a[i].real();
auto b = poly(q);
auto K = sum(a) / sum(b);
for (size_t i = 0; i < b.size(); i++)
b[i] *= K;
for (auto coeff : a)
coeffs.a.push_back(coeff.real());
for (auto coeff : b)
coeffs.b.push_back(coeff.real());
return coeffs;
}
Conclusion
Whilst the C++ implementation is not as succinct as the original Matlab code, it could be very useful for developing DSP applications such as audio plugins. This is an area I hope to return to in future posts.