Forum Rule: Always post complete source code & details to reproduce any issue!
Results 1 to 3 of 3

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

  1. #1
    Junior Member
    Join Date
    Feb 2018
    Posts
    9

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

    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!!

  2. #2
    Senior Member+ MichaelMeissner's Avatar
    Join Date
    Nov 2012
    Location
    Ayer Massachussetts
    Posts
    2,952
    Quote Originally Posted by yihanhu View Post
    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.

  3. #3
    Junior Member
    Join Date
    Apr 2014
    Location
    Seattle, WA
    Posts
    12
    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.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •