An Introduction to Biquad filters

This is a quick introduction to biquad filters. I will explain what they do and how we could implement one in c++. We wont go too deep into the mathematical part of it, because that would exceed the scope.

What is a biquad filter?

The term biquad is short for bi-quadratic, and is a common name for a two-pole, two-zero digital filter. It is a IIR filter, that can replace FIR filter, that use a greater amount of hardware resources. Biquad filters come in several forms.

Direct form 1
Direct form 2
Transposed direct form 1
Transposed direct form 2

The Direct from 1 is the most straightforward one and also the one we are going to create. It works best in a fixed point processor because it only has a single summation point.

Biquad filter direct form 1 — Direct form 1

All of the forms each come with their advantages and disadvantages. For example the Direct form 2 would be best for floating points.

Anyways, one of the reasons why biquad filters can be important is because they are fast and can create different filters depending on their coefficients.

peak
high/low/all pass
high/low shelve
band pass
band reject
smoothing

The behaviour of the filter is determined by its coefficients, which are

B1, B2
A0, A1, A2

It also needs two states Z1 and Z2, which are the last two processed samples.

How to create a biquad in code

We are going to use C++ for this one.

Biquad.h

First we will want to have an enum for convinience, to switch between different filter types.

enum BiquadType
{
    None,
    Lowpass,
    Highpass,
    Peak,
    Bandpass,
    Notch,
    Lowshelf,
    Highshelf,
};

Next we will create the Class Biquad. There is nothing scary about that part. We just need the sampleRate of our audio input we are using, the type of filter and create variables for the coefficients and states.


class Biquad
{
public:
    ~Biquad() = default;
    Biquad(const float sampleRate, const BiquadType bTypeParam, const float frequency
        ,const float qualityFactor, const float peakGain)
        :m_sampleRate(sampleRate)
        , type(bTypeParam)
        , in_z1(0.0f)
        , in_z2(0.0f)
        , coeff_a0(0.0f)
        , coeff_a1(0.0f)
        , coeff_a2(0.0f)
        , coeff_b1(0.0f)
        , coeff_b2(0.0f)
    {
        ComputeCoeff(frequency, qualityFactor, peakGain);
    }

    void ComputeCoeff(const float frequency, const float qualityFactor, const float peakGain);
    float ProcessSample(float const sample);

private:
    float m_sampleRate;
    BiquadType type;

    float in_z1;
    float in_z2;

    float coeff_a0;
    float coeff_a1;
    float coeff_a2;
    float coeff_b1;
    float coeff_b2;
};

The biggest chunck of work is calculating all the coefficients in ComputeCoeff().

ProcessSample() will be called on every sample that we pipe through that filter and will apply the right volume to each sample.

Biquad.cpp

As said the hardest part is the coefficient calculation, so we are going to tackle that one first.

To get a better grasp of that topic Will C. Pirkle has a very good section in his book Designing Audio Effect Plugins in C++ about Biquad Filters. To cut it short a biquad filter can become unstable. To prevent that the two poles need to stay inside of the unit circle. The book also already provides us with functions to calculate the right coefficient for each filter type, so we are just going to take those because otherwise that would be too much for the context of this post . In code that represents as follows.

void Biquad::ComputeCoeff(const float frequency, const float qualityFactor, const float peakGain)
{
	double const factorV = pow(10.0, (fabs(peakGain) / 20.0));
	double const factorK = tan(M_PI * (frequency / m_sampleRate));
	double const factorKSpuare = factorK * factorK;

	switch (type)
	{
		case Lowpass:
		{
			double const norm = 1.0 / (1.0 + factorK / qualityFactor + factorKSpuare);
			coeff_a0 = static_cast<float>(factorKSpuare * norm);
			coeff_a1 = static_cast<float>(2.0 * coeff_a0);
			coeff_a2 = coeff_a0;
			coeff_b1 = static_cast<float>(2.0 * (factorKSpuare - 1.0) * norm);
			coeff_b2 = static_cast<float>((1.0 - factorK / qualityFactor + factorKSpuare) 	
				    * norm);
			break;
		}
		case Highpass:
		{
			double const norm = 1.0 / (1.0 + factorK / qualityFactor + factorKSpuare);
			coeff_a0 = static_cast<float>(1 * norm);
			coeff_a1 = static_cast<float>(-2.0f * coeff_a0);
			coeff_a2 = coeff_a0;
			coeff_b1 = static_cast<float>(2 * (factorKSpuare - 1) * norm);
			coeff_b2 = static_cast<float>(norm * (1 - factorK / qualityFactor + 	
				    factorKSpuare));
			break;
		}
		case Peak:
		{
			double const x1 = factorK / qualityFactor;
			double const x2 = factorK * factorV / qualityFactor;

			if (peakGain >= 0)
			{
				double const norm = 1 / (1 + x1 + factorKSpuare);
				coeff_a0 = static_cast<float>((1 + x2 + factorKSpuare) * norm);
				coeff_a1 = static_cast<float>(2 * (factorKSpuare - 1) * norm);
				coeff_a2 = static_cast<float>((1 - x2 + factorKSpuare) * norm);
				coeff_b1 = coeff_a1;
				coeff_b2 = static_cast<float>((1 - x1 + factorKSpuare) * norm);
			}
			else
			{
				double const norm = 1 / (1 + x2 + factorKSpuare);
				coeff_a0 = static_cast<float>((1 + x1 + factorKSpuare) * norm);
				coeff_a1 = static_cast<float>((2 * (factorKSpuare - 1)) * norm);
				coeff_a2 = static_cast<float>((1 - x1 + factorKSpuare) * norm);
				coeff_b1 = coeff_a1;
				coeff_b2 = static_cast<float>((1 - x2 + factorKSpuare) * norm);
		}
		break;
	}
	case Bandpass:
	{
		double const norm = 1 / (1 + factorK / qualityFactor + factorKSpuare);
		coeff_a0 = static_cast<float>(factorK / qualityFactor * norm);
		coeff_a1 = 0;
		coeff_a2 = -coeff_a0;
		coeff_b1 = static_cast<float>(2 * (factorKSpuare - 1) * norm);
		coeff_b2 = static_cast<float>((1 - factorK / qualityFactor + factorKSpuare) *norm);
		break;
	}
	case Notch:
	{
		double const norm = 1 / (1 + factorK / qualityFactor + factorKSpuare);
		coeff_a0 = static_cast<float>((1 + factorKSpuare) * norm);
		coeff_a1 = static_cast<float>(2 * factorKSpuare - 1) * norm;
		coeff_a2 = coeff_a0;
		coeff_b1 = coeff_a1;
		coeff_b2 = static_cast<float>((1 - factorK / qualityFactor + factorKSpuare) *norm);
		break;
	}
	case Highshelf:
	{
		double const x1 = sqrt(2 * factorV) * factorK;
		double const x2 = sqrt(2) * factorK;
		if (peakGain >= 0)
		{
			double const norm = 1 / (1 + x2 + factorKSpuare);
			coeff_a0 = static_cast<float>((factorV + x1 + factorKSpuare) * norm);
			coeff_a1 = static_cast<float>(2.0 * (factorKSpuare - factorV) * norm);
			coeff_a2 = static_cast<float>((factorV - x1 + factorKSpuare) * norm);
			coeff_b1 = static_cast<float>(2.0 * (factorKSpuare - 1) * norm);
			coeff_b2 = static_cast<float>((1.0 - x2 + factorKSpuare) * norm);
		}
		else
		{
			double const norm = 1.0 / (factorV + x1 + factorKSpuare);
			coeff_a0 = static_cast<float>((1.0 + x1 + factorKSpuare) * norm);
			coeff_a1 = static_cast<float>(2.0 * (factorKSpuare - 1.0) * norm);
			coeff_a2 = static_cast<float>((1.0 - x1 + factorKSpuare) * norm);
			coeff_b1 = static_cast<float>(2.0 * (factorKSpuare - factorV) * norm);
			coeff_b2 = static_cast<float>((factorV - x1 + factorKSpuare) * norm);
		}
		break;
	}
	case Lowshelf:
	{
		double const x3 = sqrt(2) * factorK;
		double const x1 = sqrt(2.0 * factorV) * factorK;
		double const x2 = factorV * factorKSpuare;
		if (peakGain >= 0)
		{
			double const norm = 1.0 / (1.0 + x3 + factorKSpuare);
			coeff_a0 = static_cast<float>((1.0 + x1 + x2) * norm);
			coeff_a1 = static_cast<float>(2.0 * (x2 - 1.0) * norm);
			coeff_a2 = static_cast<float>((1.0 - x1 + x2) * norm);
			coeff_b1 = static_cast<float>(2.0 * (factorKSpuare - 1.0) * norm);
			coeff_b2 = static_cast<float>((1.0 - x3 + factorKSpuare) * norm);
		}
		else
		{
			double const norm = 1.0 / (1.0 + x1 + x2);
			coeff_a0 = static_cast<float>((1.0 + x3 + factorKSpuare) * norm);
			coeff_a1 = static_cast<float>(2.0 * (factorKSpuare - 1.0) * norm);
			coeff_a2 = static_cast<float>((1.0 - x3 + factorKSpuare) * norm);
			coeff_b1 = static_cast<float>(2.0 * (x2 - 1.0) * norm);
			coeff_b2 = static_cast<float>((1.0 - x1 + x2) * norm);
		}
		break;
	}
	}
}

Now we are just left with the processing part, which is fairly small and way easier to digest.

float Biquad::ProcessSample(float const sample)
{
	if (type != BiquadType::None)
	{
		float const fCurOut = sample * coeff_a0 + in_z1;
		in_z1 = sample * coeff_a1 + in_z2 - coeff_b1 * fCurOut;
		in_z2 = sample * coeff_a2 - coeff_b2 * fCurOut;
		return fCurOut;
	}
	else
	{
		return sample;
	}
}

Conclusion

Biquad filters are a great option to FIR filters to create runtime effects, as once the filters are set, it is a relative cheap calculation. With that we can create custom filterbanks that do frequency shaping of the signal. It even allows us to create a spatializer that simulates a HRTF with it, but thats a topic for a later post.

Here is a website where you can play around with values and see how a biquad filter behaves.

An introduction to Biquad filters

What is a biquad filter?

How to create a biquad in code

Biquad.h

Biquad.cpp

Conclusion

Recent Posts

GCRadix