Draw frequency response of filter from transfer function


#1

Hi,
I’m using the gamma synthesis library (https://github.com/LancePutnam/Gamma) for my audio processing. It’s a very nice and easy to use library. The only problem I’m having: Drawing the frequency response of the filter. Right now I’m doing it very inefficently by filtering an impulse, applying a fourier transform and interpolating over each frequency band amplitude. To get smooth results I need a pretty big fft buffer. It’s really not efficent at all.

In the documention of the filter I found a link with the transfer functions of all the different filter types and a lot of extra information. I have absolutely no clue of how transfer functions work. But as far as I can see, it should be possible to calculate the frequency responses directly from the information.

The said document can be found here: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Could somebody help me transform this into a plotting function?
Help would be very much appreciated!

Thanks!


#2

The transfer function is exactly what you want. Given:

H(z) = (...) / (...)

You can get the complex frequency response at some angle w by doing H(w).

This is the code for evaluating the transfer function of a standard biquad (RBJ’s filters are biquads). It is basically expanding and reformulating the transfer function equation and optimizing it for real output. Given a biquad consisting of a0, a1, a2, b1, b2, the transfer function can be evaluated at normalized frequency w0 like so:

		auto const piw0 = w0 * cpl::simd::consts<T>::pi;
		auto const cosw = std::cos(piw0);
		auto const sinw = std::sin(piw0);

		auto square = [](auto z) { return z * z; };

		auto const numerator = sqrt(square(a0*square(cosw) - a0*square(sinw) + a1*cosw + a2) + square(2 * a0*cosw*sinw + a1*(sinw)));
		auto const denominator = sqrt(square(square(cosw) - square(sinw) + b1*cosw + b2) + square(2 * cosw*sinw + b1*(sinw)));

		/*auto const w = std::sin(cpl::simd::consts<T>::pi * w / 2);
		auto const phi = w * w;
		auto const phi2 = phi * phi;

		// black magic. supposedly gives better precision than direct evaluation of H(z)
		auto const numerator = 10 * std::log10((b0 + b1 + b2) * (b0 + b1 + b2) - 4 * (b0 * b1 + 4 * b0 * b2 + b1 * b2) * phi + 16 * b0 * b2 * phi2);
		auto const denominator = -10 * std::log10((a0 + a1 + a2) * (a0 + a1 + a2) - 4 * (a0 * a1 + 4 * a0 * a2 + a1 * a2) * phi + 16 * a0 * a2 * phi2);
		return numerator + denominator; */

		return 20 * std::log10(numerator / denominator);

The last part calculates it to decibels, but can obviously be removed.


#3

Ok, so I tried working on this, but there’s some things I don’t really understand. Please bear with me.

You’re saying I can evaluate the function at “normalized frequency w0”. I looked it up and the definition varied. I use an interval from 0-1. Is this correct?

Also you’re saying “given a biquad consisting of a0, a1, a2, b1, b2”. Where has b0 gone in this example? Did you use it to normalize the other coefficients?

Right now my function looks like this. The coefficients are from the lowpass filter. filter.freq() is the cutoff frequency and the filter.res() the resonance.

float calcMagnitude(gam::Biquad<> filter, float freq) {
auto w = 2*M_PI*filter.freq()/44100.f;
auto alpha = sin(w)/(2*filter.res());

auto b0 = (1-cos(w))/2.f;
auto b1 = 1-cos(w);
auto b2 = (1-cos(w))/2.f;
auto a0 = 1+alpha;
auto a1 = -2*cos(w);
auto a2 = 1-alpha;

auto const piw0 = freq * M_PI;
auto const cosw = std::cos(piw0);
auto const sinw = std::sin(piw0);

auto square = [](double z) { return z * z; };

auto const numerator = sqrt(square(a0*square(cosw) - a0*square(sinw) + a1*cosw + a2) + square(2 * a0*cosw*sinw + a1*(sinw)));
auto const denominator = sqrt(square(square(cosw) - square(sinw) + b1*cosw + b2) + square(2 * cosw*sinw + b1*(sinw)));

return numerator / denominator;
}

I call it with freq values from 0 - 1 in 0.001 steps, since that for me is “normalized frequency”. (I know, I shouldn’t recalculate the coefficents every time. This is just for testing.)

Right now it’s not really working. The results are totally of. Help is very much appreciated!


#4

Nevermind, it was a drawing error on my side. Works like a charm! Thank you very much!