A question about the teensy Audio Library Biquad filter header file source file

Status
Not open for further replies.

yihanhu

Member
Code:
void setCoefficients(uint32_t stage, const int *coefficients);
	void setCoefficients(uint32_t stage, const double *coefficients) {
		int coef[5];
		// 1073741824 bytes = 1 GB
		// Google! Prob has something to do with memory allocation 
		// Ask MUS171 Professor
		coef[0] = coefficients[0] * 1073741824.0;
		coef[1] = coefficients[1] * 1073741824.0;
		coef[2] = coefficients[2] * 1073741824.0;
		coef[3] = coefficients[3] * 1073741824.0;
		coef[4] = coefficients[4] * 1073741824.0;
		setCoefficients(stage, coef);
	}

	// Compute common filter functions
	// http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
	void setLowpass(uint32_t stage, float frequency, float q = 0.7071) {
		int coef[5];
		double w0 = frequency * (2 * 3.141592654 / AUDIO_SAMPLE_RATE_EXACT);
		double sinW0 = sin(w0);
		double alpha = sinW0 / ((double)q * 2.0);
		double cosW0 = cos(w0);
		double scale = 1073741824.0 / (1.0 + alpha);
		/* b0 */ coef[0] = ((1.0 - cosW0) / 2.0) * scale;
		/* b1 */ coef[1] = (1.0 - cosW0) * scale;
		/* b2 */ coef[2] = coef[0];
		/* a1 */ coef[3] = (-2.0 * cosW0) * scale;
		/* a2 */ coef[4] = (1.0 - alpha) * scale;
		setCoefficients(stage, coef);
	}

I understand that we are using bilinear transformation to obtain the coefficients. But I am not quite sure about the scaling factor here which is 1073741824.0. If someone could explain the reason that we include such number in the code, I would very much appreciate it. Thank you!!
 
Code:
void setCoefficients(uint32_t stage, const int *coefficients);
	void setCoefficients(uint32_t stage, const double *coefficients) {
		int coef[5];
		// 1073741824 bytes = 1 GB
		// Google! Prob has something to do with memory allocation 
		// Ask MUS171 Professor
		coef[0] = coefficients[0] * 1073741824.0;
		coef[1] = coefficients[1] * 1073741824.0;
		coef[2] = coefficients[2] * 1073741824.0;
		coef[3] = coefficients[3] * 1073741824.0;
		coef[4] = coefficients[4] * 1073741824.0;
		setCoefficients(stage, coef);
	}

	// Compute common filter functions
	// http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
	void setLowpass(uint32_t stage, float frequency, float q = 0.7071) {
		int coef[5];
		double w0 = frequency * (2 * 3.141592654 / AUDIO_SAMPLE_RATE_EXACT);
		double sinW0 = sin(w0);
		double alpha = sinW0 / ((double)q * 2.0);
		double cosW0 = cos(w0);
		double scale = 1073741824.0 / (1.0 + alpha);
		/* b0 */ coef[0] = ((1.0 - cosW0) / 2.0) * scale;
		/* b1 */ coef[1] = (1.0 - cosW0) * scale;
		/* b2 */ coef[2] = coef[0];
		/* a1 */ coef[3] = (-2.0 * cosW0) * scale;
		/* a2 */ coef[4] = (1.0 - alpha) * scale;
		setCoefficients(stage, coef);
	}

I understand that we are using bilinear transformation to obtain the coefficients. But I am not quite sure about the scaling factor here which is 1073741824.0. If someone could explain the reason that we include such number in the code, I would very much appreciate it. Thank you!!

FWIW, is 1073741824 is 0x40000000 (or 0x1UL << 30). I don't know why it is used.
 
To get good performance on older Teensy's, most of the audio code uses fixed-point instead of floating-point math.

In this case, biquad coefs are converted to fixed-point with sign bit, one integer bit, and 30 fraction bits. This effectively represents [-2.0, +2.0] with enough fraction bits to handle low-cutoff high-Q filters where coefs are very sensitive to quantization. DSP people usually call this Q30.

You can see how it works by following the scaling of the biquad loop:
Code:
            in2 = *data;
            sum = signed_multiply_accumulate_32x16b(sum, b0, in2);
            sum = signed_multiply_accumulate_32x16t(sum, b1, bprev);
            sum = signed_multiply_accumulate_32x16b(sum, b2, bprev);
            sum = signed_multiply_accumulate_32x16t(sum, a1, aprev);
            sum = signed_multiply_accumulate_32x16b(sum, a2, aprev);
            out2 = signed_saturate_rshift(sum, 16, 14);
            sum &= 0x3FFF;

You can consider the input samples to be Q0 (no fraction bits).
Each multiply-add does sum += (a * b) >> 16, so the scaling is (Q30 * Q0) >> 16 = Q14 and sum keeps 14 fraction bits.
The final shift discards the 14 fraction bits and saturates the output to 16 integer bits.

A problem with just throwing away the fraction bits, is output is reused by the feedback path (aprev in the code above). A low-freq high-Q filter recirculates for a long time, and the accumulated quantization noise ruins the SNR. One solution is to keep fraction bits in the feedback path, but that requires a 32x32 multiply instead of 32x16 SIMD multiply which requires fewer loads.

You get almost as good result by using first-order quantization-error feedback. The above code implements this via "fraction saving." It is important to note how a signed fixed-point shift (sum >> 14) implements "truncation towards -infinity". That is, 1.xxx truncates to 1 and -1.xxx truncates to -2 so quantization error is guaranteed to be a positive fraction. The (sum &= 0x3FFF) extracts this fraction and saves it for the next iteration. Then by replacing the first multiply with a multiply-add, it feeds back into the next sum essentially for free.

Hope this helps.
 
Status
Not open for further replies.
Back
Top