SOS blowing up with K-weighting filter coefficients

Questions
1

I’m trying to implement the “K” frequency weighting filter as described in ITU-R BS.1770-4 (p. 3-5) using SOS.ar. This filter is a simple cascade of two second-order filter sections (biquads), namely a high shelving filter and a high-pass filter. The biquad coefficients in the paper are given for a fixed sample rate (fs=48000).
I’ve made sure

  • that my server is running at that exact sample rate,
  • that the coefficients are given in the correct order (the feedback/feedforward coefficient naming is reversed in the paper, so SOS.ar’s a0, a1, a2 are called b0, b1, b2 in the paper), and
  • that these exact coefficients are working as expected in Max’s biquad~ object,

but SOS.ar still blows up. The same thing happens with LTI.ar from sc3-plugins. I’ve tested this with both the high shelf and lowpass coefficients given in the paper. Here’s my code, make sure to mute your server outs before executing the second block:

s = s ?? { Server.default }; 
s.sampleRate; // -> 48000.0, no issues there

// !!!
// WARNING: MUTE BEFORE EXECUTING THIS BLOCK
// !!!
(
// valid filter coefficients for a sample rate of 48000 Hz
// given in ITU-R BS.1770-4
var a0, a1, a2, b1, b2;

// values from paper
// shelf
b1 = 1.69065929318241.neg;
b2 = 0.73248077421585;
a0 = 1.53512485958697;
a1 = 2.69169618940638.neg;
a2 = 1.19839281085285;

// hp
// b1 = 1.99004745483398.neg;
// b2 = 0.99007225036621;
// a0 = 1.0;
// a1 = 2.0.neg;
// a2 = 1.0;

{
	SOS.ar(WhiteNoise.ar, a0, a1, a2, b1, b2) * -30.dbamp ! 2;
}.play
)

Does anyone know why this happens? I know FluidLoudness.kr exists, but I want to have a stab at rolling my own “perceptually” weighted filters for other purporses (e.g. to use in a compressor before the envelope calculation stage).

2

In Figure 3, there are minus signs for a1 and a2. Therefore their filter has the following difference equation (using their a’s and b’s):

y[t] = b0 * x[t] + b1 * x[t-1] + b2 * x[t-2] - a1 * y[t-1] - a2 * y[t-2]

The SOS help file, in addition to swapping the a’s and b’s, inverts the sign of the y[t] coefficients on the right-hand side. Therefore, you need to flip the sign of your b1’s and b2’s.

ITU’s/Max’s convention is that way because when using the Z-transform, the transfer function becomes

H(z) = (b0 + b1 z^-1 + b2 z^-2) / (1 + a1 z^-1 + a2 z^-2)

and all the signs end up positive.

1 Like
3

I did not know this, but it makes sense! Thanks for the explanation. I noticed the minus sign pretty much right after posting the OP, that somehow always happens… (figuring out the solution to a problem right after asking on a forum)

4

I’ve ported this piece of code by Samuel Gaehwiler (Klangfreund) to convert a set of biquad coeffs given for a fixed sample rate into any other sample rate (e.g. from fs=48000 to fs=44100), if anyone needs it.

(
// this block of code converts a set of biquad coefficients given for a fixed sample rate
// to a new sample rate. this is used to convert the "k" frequency weighting filter coefficients
// given in ITU-R BS.1770-4, which are given for a fixed fs=48k, into any other sample rate.

// set of valid filter coefficients for a sample rate of 48000 Hz
var a1_48k, a2_48k, b0_48k, b1_48k, b2_48k;

// calculated filter coefficients for the used sample rate
var a1, a2, b0, b1, b2;

// filter parameters which appear in the transfer function
// these are used to calculate the new coefficients
var k, q, vh, vb, vl;

// old/new sample rates
var oldSampleRate = 48000, newSampleRate = 44100;

// function which carries out computation of filter parameters
var calcFilterParams;

// print coefficients, plot magnitude response
var verbose = true, plot = true;

// play unfiltered/filtered audio to compare
var play = false;

// values from paper
// high shelving filter
// gain = 4dB, f_c ~= 1681.9 Hz
a1_48k =  1.69065929318241.neg;
a2_48k =  0.73248077421585;
b0_48k =  1.53512485958697;
b1_48k =  2.69169618940638.neg;
b2_48k =  1.19839281085285;

// high-pass filter
// f_c ~= 37.5 Hz, Q = 0.5
// a1_48k =  1.99004745483398.neg;
// a2_48k =  0.99007225036621;
// b0_48k =  1.0;
// b1_48k =  2.0.neg;
// b2_48k =  1.0;


calcFilterParams = { |a1, a2, b0, b1, b2|
	var kOverQ = (2 - (2 * a2)) / (a2 - a1 + 1);
	k = sqrt((a1 + a2 + 1) / (a2 - a1 + 1));
	
	q = k / kOverQ;
	vb = (b0 - b2) / (1 - a2);
	vh = (b0 - b1 + b2) / (a2 - a1 + 1);
	vl = (b0 + b1 + b2) / (a1 + a2 + 1);
};

// run filter param calculation
calcFilterParams.(a1_48k, a2_48k, b0_48k, b1_48k, b2_48k);

// calculate filter coefficients for new sample rate
if (newSampleRate == oldSampleRate) {
	#a1, a2, b0, b1, b2 = [a1_48k, a2_48k, b0_48k, b1_48k, b2_48k];
} {
	var kSq, commonFactor;
	
	k = tan(atan(k) * oldSampleRate / newSampleRate);
	kSq = k.squared;
	commonFactor = 1 / (1 + (k/q) + (k.squared));

	a1 = 2 * (kSq - 1) * commonFactor;
	a2 = (1 - (k/q) + kSq) * commonFactor;
	b0 = (vh + (vb*k/q) + (vl*kSq)) * commonFactor;
	b1 = 2 * (vl*kSq - vh) * commonFactor;
	b2 = (vh - (vb*k/q) + (vl*kSq)) * commonFactor;
	
	if (verbose) {
		// compare coefficients
		"old coefficients (fs = 48k)".postln;
		[a1_48k, a2_48k, b0_48k, b1_48k, b2_48k].debug;
		"".postln;
		"new coefficients (fs = %)\n".postf(newSampleRate);
		[a1, a2, b0, b1, b2].debug;
		
		if (plot) {
			// plot magnitude response
			// requires MagnitudeResponse quark (https://scsynth.org/t/magnituderesponseview-filterkernels/5907)
			
			// old coeffs
			// SOS.magnitudeResponse([b0_48k, b1_48k, b2_48k, a1_48k.neg, a2_48k.neg], 2048, 20, 2e4, oldSampleRate).plot
			
			// new coeffs
			SOS.magnitudeResponse([b0, b1, b2, a1.neg, a2.neg], 2048, 20, 2e4, newSampleRate).plot
		};
	};
};


// compare unfiltered/filtered audio (using white noise)
if (play) {
	{
		var in = WhiteNoise.ar * -20.dbamp;
		var filtered = SOS.ar(in, b0_48k, b1_48k, b2_48k, a1_48k.neg, a2_48k.neg); // a1/a2 need to be negated for SOS.ar
		Select.ar(MouseX.kr > 0.5, [in, filtered]) ! 2
	}.play;
	// }.scopeResponse; // magnitude reponse
};
)
1 Like
5

Hey there! I have wondered about the same thing, I’m working in max atm. But I don’t get the right results with biquad~ or filtergraph~ in max as specified in the ITU pdf. I tried switching the a and b coefficients, which at least resulted in a similar freq response for the low cut filter as described in the specification, well it’s a low cut at least :slight_smile: But I can’t get the high shelve to work.

Also tried inverting the signs of my b coefficients, but it’s still blowing up the filter.
Any advice would be appreciated. Especially from @bovil43810 since you made it work in max :slight_smile:

Thanks!