-
-
Save vojd/42bac90414650f831c57b00b95964d31 to your computer and use it in GitHub Desktop.
Python implementation of "Cookbook formulae for audio EQ biquad filter coefficients"
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -*- coding: utf-8 -*- | |
| """ | |
| Created on Thu Mar 28 18:51:00 2013 | |
| UNFINISHED AND BUGGY | |
| Python/SciPy implementation of the filters described in | |
| "Cookbook formulae for audio EQ biquad filter coefficients" | |
| by Robert Bristow-Johnson | |
| http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt | |
| These functions will output analog or digital transfer functions, deriving | |
| the latter using the bilinear transform, as is done in the reference. | |
| Overall gain parameters are not included. | |
| "BLT frequency warping has been taken into account for | |
| both significant frequency relocation (this is the normal "prewarping" that | |
| is necessary when using the BLT) and for bandwidth readjustment (since the | |
| bandwidth is compressed when mapped from analog to digital using the BLT)." | |
| TODO: combine lowpass and highpass? and bandpass? | |
| TODO: generate analog poles/zeros prototypes and convert them or output them directly? | |
| TODO: Use ordinary frequency instead of rad/s for analog filters? angular | |
| matches scipy, but these are usually used in audio. Compare with CSound | |
| functions, etc. | |
| TODO: sane defaults for Q for all filters | |
| TODO: Try to think of better names than "outer", "constantq", "skirt", etc | |
| TODO: Bandwidth is wrong for high-frequency peaking digital filters, | |
| despite using the equations in the cookbook | |
| TODO: functions should accept Q, BW, or S, since these are not trivially | |
| derived otherwise? | |
| Q (the EE kind of definition, except for peakingEQ in which A*Q is | |
| the classic EE Q. That adjustment in definition was made so that | |
| a boost of N dB followed by a cut of N dB for identical Q and | |
| f0/Fs results in a precisely flat unity gain filter or "wire".) | |
| _or_ BW, the bandwidth in octaves (between -3 dB frequencies for BPF | |
| and notch or between midpoint (dBgain/2) gain frequencies for | |
| peaking EQ) | |
| _or_ S, a "shelf slope" parameter (for shelving EQ only). When S = 1, | |
| the shelf slope is as steep as it can be and remain monotonically | |
| increasing or decreasing gain with frequency. The shelf slope, in | |
| dB/octave, remains proportional to S for all other values for a | |
| fixed f0/Fs and dBgain. | |
| Then compute a few intermediate variables: | |
| alpha = sin(w0)/(2*Q) (case: Q) | |
| = sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0) ) (case: BW) | |
| = sin(w0)/2 * sqrt( (A + 1/A)*(1/S - 1) + 2 ) (case: S) | |
| FYI: The relationship between bandwidth and Q is | |
| 1/Q = 2*sinh(ln(2)/2*BW*w0/sin(w0)) (digital filter w BLT) | |
| or 1/Q = 2*sinh(ln(2)/2*BW) (analog filter prototype) | |
| """ | |
| from __future__ import division | |
| from math import pi, tan, sinh | |
| from math import log as ln | |
| from cmath import sqrt | |
| import numpy as np | |
| from scipy.signal import tf2zpk, tf2ss, lp2lp, bilinear | |
| def _transform(b, a, Wn, analog, output): | |
| """ | |
| Shift prototype filter to desired frequency, convert to digital with | |
| pre-warping, and return in various formats. | |
| """ | |
| Wn = np.asarray(Wn) | |
| if not analog: | |
| if np.any(Wn < 0) or np.any(Wn > 1): | |
| raise ValueError("Digital filter critical frequencies " | |
| "must be 0 <= Wn <= 1") | |
| fs = 2.0 | |
| warped = 2 * fs * tan(pi * Wn / fs) | |
| else: | |
| warped = Wn | |
| # Shift frequency | |
| b, a = lp2lp(b, a, wo=warped) | |
| # Find discrete equivalent if necessary | |
| if not analog: | |
| b, a = bilinear(b, a, fs=fs) | |
| # Transform to proper out type (pole-zero, state-space, numer-denom) | |
| if output in ('zpk', 'zp'): | |
| return tf2zpk(b, a) | |
| elif output in ('ba', 'tf'): | |
| return b, a | |
| elif output in ('ss', 'abcd'): | |
| return tf2ss(b, a) | |
| elif output in ('sos'): | |
| raise NotImplementedError('second-order sections not yet implemented') | |
| else: | |
| raise ValueError('Unknown output type {0}'.format(output)) | |
| def lowpass(Wn, Q=1/sqrt(2), analog=False, output='ba'): | |
| """ | |
| Generic biquad lowpass filter design | |
| Design a 2nd-order analog or digital lowpass filter with variable Q and | |
| return the filter coefficients. | |
| Analog prototype: H(s) = 1 / (s**2 + s/Q + 1) | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Corner frequency of the filter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * 1/sqrt(2) (default) is a Butterworth filter, with maximally-flat | |
| passband | |
| * 1/sqrt(3) is a Bessel filter, with maximally-flat group delay. | |
| * 1/2 is a Linkwitz-Riley filter, used to make lowpass and highpass | |
| sections that sum flat to unity gain. | |
| analog : bool, optional | |
| When True, return an analog filter, otherwise a digital filter is | |
| returned. | |
| output : {'ba', 'zpk', 'ss'}, optional | |
| Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or | |
| state-space ('ss'). | |
| Default is 'ba'. | |
| Returns | |
| ------- | |
| b, a : ndarray, ndarray | |
| Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. | |
| Only returned if ``output='ba'``. | |
| z, p, k : ndarray, ndarray, float | |
| Zeros, poles, and system gain of the IIR filter transfer | |
| function. Only returned if ``output='zpk'``. | |
| """ | |
| # H(s) = 1 / (s**2 + s/Q + 1) | |
| b = np.array([1]) | |
| a = np.array([1, 1/Q, 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def highpass(Wn, Q=1/sqrt(2), analog=False, output='ba'): | |
| """ | |
| Generic biquad highpass filter design | |
| Design a 2nd-order analog or digital highpass filter with variable Q and | |
| return the filter coefficients. | |
| Analog prototype: H(s) = s**2 / (s**2 + s/Q + 1) | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Corner frequency of the filter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * 1/sqrt(2) (default) is a Butterworth filter, with maximally-flat | |
| passband | |
| * 1/sqrt(3) is a Bessel filter, with maximally-flat group delay. | |
| * 1/2 is a Linkwitz-Riley filter, used to make lowpass and highpass | |
| sections that sum flat to unity gain. | |
| analog : bool, optional | |
| When True, return an analog filter, otherwise a digital filter is | |
| returned. | |
| output : {'ba', 'zpk', 'ss'}, optional | |
| Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or | |
| state-space ('ss'). | |
| Default is 'ba'. | |
| Returns | |
| ------- | |
| b, a : ndarray, ndarray | |
| Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. | |
| Only returned if ``output='ba'``. | |
| z, p, k : ndarray, ndarray, float | |
| Zeros, poles, and system gain of the IIR filter transfer | |
| function. Only returned if ``output='zpk'``. | |
| """ | |
| # H(s) = s**2 / (s**2 + s/Q + 1) | |
| b = np.array([1, 0, 0]) | |
| a = np.array([1, 1/Q, 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def bandpass(Wn, Q=1, type='skirt', analog=False, output='ba'): | |
| """ | |
| Biquad bandpass filter design | |
| Design a 2nd-order analog or digital bandpass filter with variable Q and | |
| return the filter coefficients. | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Center frequency of the filter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * sqrt(2) is 1 octave wide | |
| type : {'skirt', 'peak'}, optional | |
| The type of filter. | |
| ``skirt`` | |
| Type 1 (default), has a constant skirt gain, with peak gain = Q | |
| Transfer function: H(s) = s / (s**2 + s/Q + 1) | |
| ``peak`` | |
| Type 2, has a constant peak gain of 0 dB, and the skirt changes | |
| with the Q. | |
| Transfer function: H(s) = (s/Q) / (s**2 + s/Q + 1) | |
| analog : bool, optional | |
| When True, return an analog filter, otherwise a digital filter is | |
| returned. | |
| output : {'ba', 'zpk', 'ss'}, optional | |
| Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or | |
| state-space ('ss'). | |
| Default is 'ba'. | |
| Returns | |
| ------- | |
| b, a : ndarray, ndarray | |
| Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. | |
| Only returned if ``output='ba'``. | |
| z, p, k : ndarray, ndarray, float | |
| Zeros, poles, and system gain of the IIR filter transfer | |
| function. Only returned if ``output='zpk'``. | |
| """ | |
| if type in (1, 'skirt'): | |
| # H(s) = s / (s**2 + s/Q + 1) | |
| b = np.array([0, 1, 0]) | |
| elif type in (2, 'peak'): | |
| # H(s) = (s/Q) / (s**2 + s/Q + 1) | |
| b = np.array([0, 1/Q, 0]) | |
| else: | |
| raise ValueError('"%s" is not a known bandpass type' % type) | |
| a = np.array([1, 1/Q, 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def notch(Wn, Q=10, analog=False, output='ba'): | |
| """ | |
| Biquad notch filter design | |
| Design a 2nd-order analog or digital notch filter with variable Q and | |
| return the filter coefficients. | |
| The notch differs from a peaking cut filter in that the gain at the | |
| notch center frequency is 0, or -Inf dB. | |
| Transfer function: H(s) = (s**2 + 1) / (s**2 + s/Q + 1) | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Center frequency of the filter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * sqrt(2) is 1 octave wide | |
| analog : bool, optional | |
| When True, return an analog filter, otherwise a digital filter is | |
| returned. | |
| output : {'ba', 'zpk', 'ss'}, optional | |
| Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or | |
| state-space ('ss'). | |
| Default is 'ba'. | |
| Returns | |
| ------- | |
| b, a : ndarray, ndarray | |
| Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. | |
| Only returned if ``output='ba'``. | |
| z, p, k : ndarray, ndarray, float | |
| Zeros, poles, and system gain of the IIR filter transfer | |
| function. Only returned if ``output='zpk'``. | |
| """ | |
| # H(s) = (s**2 + 1) / (s**2 + s/Q + 1) | |
| b = np.array([1, 0, 1]) | |
| a = np.array([1, 1/Q, 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def allpass(Wn, Q=1, analog=False, output='ba'): | |
| """ | |
| Biquad allpass filter design | |
| Design a 2nd-order analog or digital allpass filter with variable Q and | |
| return the filter coefficients. | |
| Transfer function: H(s) = (s**2 - s/Q + 1) / (s**2 + s/Q + 1) | |
| Wn is center frequency | |
| """ | |
| # H(s) = (s**2 - s/Q + 1) / (s**2 + s/Q + 1) | |
| b = np.array([1, -1/Q, 1]) | |
| a = np.array([1, 1/Q, 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def peaking(Wn, dBgain, Q=None, BW=None, type='half', analog=False, output='ba'): | |
| """ | |
| Biquad peaking filter design | |
| Design a 2nd-order analog or digital peaking filter with variable Q and | |
| return the filter coefficients. Used in graphic or parametric EQs. | |
| Transfer function: H(s) = (s**2 + s*(Az/Q) + 1) / (s**2 + s/(Ap*Q) + 1) | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Center frequency of the filter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| dBgain : float | |
| The gain at the center frequency, in dB. Positive for boost, | |
| negative for cut. | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * Q = sqrt(2) (default) produces a bandwidth of 1 octave | |
| ftype : {'half', 'constant'}, optional | |
| Where on the curve to measure the bandwidth of the filter. | |
| ``half`` | |
| Bandwidth is defined using the points on the curve at which the | |
| gain in dB is half of the peak gain. This is the method used in | |
| "Cookbook formulae for audio EQ biquad filter coefficients" | |
| ``constant`` | |
| Bandwidth is defined using the points -3 dB down from the peak | |
| gain (or +3 dB up from the cut gain), maintaining constant Q | |
| regardless of center frequency or boost gain. This is | |
| symmetrical in dB, so that a boost and cut with identical | |
| parameters sum to unity gain. | |
| This is the method used in "Constant-Q" hardware equalizers. | |
| [ref: http://www.rane.com/note101.html] | |
| Klark Teknik calls this "symmetrical Q" http://www.klarkteknik.com/faq-06.php | |
| constant Q asymmetrical | |
| constant Q for both boost and cut, which makes them asymmetrical (not implemented) | |
| Half-gain Hybrid | |
| Defined symmetrical at half gain point except for 3 dB or less (not implemented) | |
| analog : bool, optional | |
| When True, return an analog filter, otherwise a digital filter is | |
| returned. | |
| output : {'ba', 'zpk', 'ss'}, optional | |
| Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or | |
| state-space ('ss'). | |
| Default is 'ba'. | |
| Notes | |
| ----- | |
| Due to bilinear transform, this is always 0 dB at fs/2, but it would be | |
| better if the curve fell off symmetrically. | |
| Orfanidis describes a digital filter that more accurately matches the | |
| analog filter, but it is far more complicated. | |
| Orfanidis, Sophocles J., "Digital Parametric Equalizer Design with | |
| Prescribed Nyquist-Frequency Gain" | |
| """ | |
| if Q is None and BW is None: | |
| BW = 1 # octave | |
| if Q is None: | |
| # w0 = Wn | |
| # Q = 1/(2*sinh(ln(2)/2*BW*w0/sin(w0))) # digital filter w BLT | |
| Q = 1/(2*sinh(ln(2)/2*BW)) # analog filter prototype | |
| # TODO: In testing, neither of these is even close to correct near | |
| # fs/2, and the difference between them is very small | |
| if type in ('half'): | |
| A = 10.0**(dBgain/40.0) # for peaking and shelving EQ filters only | |
| Az = A | |
| Ap = A | |
| elif type in ('constantq'): | |
| A = 10.0**(dBgain/20.0) | |
| if dBgain > 0: # boost | |
| Az = A | |
| Ap = 1 | |
| else: # cut | |
| Az = 1 | |
| Ap = A | |
| else: | |
| raise ValueError('"%s" is not a known peaking type' % type) | |
| # H(s) = (s**2 + s*(Az/Q) + 1) / (s**2 + s/(Ap*Q) + 1) | |
| b = np.array([1, Az/Q, 1]) | |
| a = np.array([1, 1/(Ap*Q), 1]) | |
| return _transform(b, a, Wn, analog, output) | |
| def shelf(Wn, dBgain, S=1, btype='low', ftype='half', analog=False, output='ba'): | |
| """ | |
| Biquad shelving filter design | |
| Design a 2nd-order analog or digital shelving filter with variable slope | |
| and return the filter coefficients. | |
| Parameters | |
| ---------- | |
| Wn : float | |
| Turnover frequency of the filter, defined by the `ftype` parameter. | |
| For digital filters, `Wn` is normalized from 0 to 1, where 1 is the | |
| Nyquist frequency, pi radians/sample. (`Wn` is thus in | |
| half-cycles / sample.) | |
| For analog filters, `Wn` is an angular frequency (e.g. rad/s). | |
| dBgain : float | |
| The gain at the center frequency, in dB. Positive for boost, | |
| negative for cut. | |
| Q : float | |
| Quality factor of the filter. Examples: | |
| * Q fdsafda | |
| ftype : {'half', 'outer', 'inner'}, optional | |
| fpoint? | |
| fdef? | |
| Definition of the filter's turnover frequency | |
| ``half`` | |
| Wn is defined as the point on the curve at which the | |
| gain in dB is half of the shelf gain, or midway between the | |
| filter's pole and zero. This method is used in | |
| "Cookbook formulae for audio EQ biquad filter coefficients" | |
| ``outer`` | |
| Wn is defined as the point 3 dB up or down from the shelf's | |
| plateau. | |
| This is symmetrical in dB, so that a boost and cut with identical | |
| parameters sum to unity gain. | |
| This is defined using the location of the outer pole or zero of | |
| the filter (the lower of the two for a low shelf, higher of the | |
| two for a high shelf), so will not be exactly 3 dB at lower shelf | |
| gains. This method is used in ____ hardware audio equalizers. | |
| ``inner`` | |
| Wn is defined as the point 3 dB up or down from unity gain. | |
| This is symmetrical in dB, so that a boost and cut with identical | |
| parameters sum to unity gain. | |
| btype : {'low', 'high'}, optional | |
| Band type of the filter, low shelf or high shelf. | |
| ftype is the meaning of f, either midpoint of slope, fstop or fturnover | |
| turnover frequency at large boost/cuts, this is 3 dB away from unity gain | |
| stop frequency at large boost/cuts, this is 3 dB away from plateau | |
| tonmeister defines outer as fstop and inner as fturnover | |
| as does http://www.soundonsound.com/sos/dec05/articles/qa1205_3.htm | |
| Understanding Audio defines turnover as outer | |
| as does ems.music.utexas.edu/dwnld/mus329j10/Filter%20Basics.ppt | |
| also calls it knee | |
| R is transition ratio fstop/fturnover. at R=1, fstop = fturnover | |
| If the transition ratio is less than 1, then the filter is a low shelving filter. If the transition ratio is greater than 1, then the filter is a high shelving filter. | |
| highShelf: H(s) = A * (A*s**2 + (sqrt(A)/Q)*s + 1)/( s**2 + (sqrt(A)/Q)*s + A) | |
| lowShelf: H(s) = A * ( s**2 + (sqrt(A)/Q)*s + A)/(A*s**2 + (sqrt(A)/Q)*s + 1) | |
| 2*sqrt(A)*alpha = sin(w0) * sqrt( (A**2 + 1)*(1/S - 1) + 2*A ) | |
| is a handy intermediate variable for shelving EQ filters. | |
| The relationship between shelf slope and Q is | |
| 1/Q = sqrt((A + 1/A)*(1/S - 1) + 2) | |
| f0 shelf midpoint frequency | |
| _or_ S, a "shelf slope" parameter (for shelving EQ only). When S = 1, | |
| the shelf slope is as steep as it can be and remain monotonically | |
| increasing or decreasing gain with frequency. The shelf slope, in | |
| dB/octave, remains proportional to S for all other values for a | |
| fixed f0/Fs and dBgain. | |
| """ | |
| Q = None # TODO: Maybe this should be a function parameter? | |
| if ftype in ('mid', 'half'): | |
| A = 10.0**(dBgain/40.0) # for peaking and shelving EQ filters only | |
| if Q is None: | |
| Q = 1/sqrt((A + 1/A)*(1/S - 1) + 2) | |
| Az = A | |
| Ap = A | |
| elif ftype in ('outer'): | |
| A = 10.0**(dBgain/20.0) | |
| if Q is None: | |
| Q = 1/sqrt((A + 1/A)*(1/S - 1) + 2) | |
| if dBgain > 0: # boost | |
| Az = A | |
| Ap = 1 | |
| else: # cut | |
| Az = 1 | |
| Ap = A | |
| elif ftype in ('inner'): | |
| A = 10.0**(dBgain/20.0) | |
| if Q is None: | |
| Q = 1/sqrt((A + 1/A)*(1/S - 1) + 2) | |
| if dBgain > 0: # boost | |
| Az = 1 | |
| Ap = A | |
| else: # cut | |
| Az = A | |
| Ap = 1 | |
| else: | |
| raise ValueError('"%s" is not a known shelf type' % ftype) | |
| if btype == 'low': | |
| # H(s) = A * ( s**2 + (sqrt(A)/Q)*s + A)/(A*s**2 + (sqrt(A)/Q)*s + 1) | |
| b = Ap * np.array([1, sqrt(Az)/Q, Az]) | |
| a = np.array([Ap, sqrt(Ap)/Q, 1]) | |
| elif btype == 'high': | |
| # H(s) = A * (A*s**2 + (sqrt(A)/Q)*s + 1)/( s**2 + (sqrt(A)/Q)*s + A) | |
| b = Ap * np.array([Az, sqrt(Az)/Q, 1]) | |
| a = np.array([1, sqrt(Ap)/Q, Ap]) | |
| else: | |
| raise ValueError('"%s" is not a known shelf type' % btype) | |
| return _transform(b, a, Wn, analog, output) | |
| if __name__ == "__main__": | |
| from scipy.signal import freqs, freqz | |
| from numpy import log10 | |
| from matplotlib.pyplot import (title, grid, show, plot, xlabel, | |
| ylabel, xscale, figure, yticks, xlim, | |
| margins) | |
| for ftype in ('half', 'constantq'): | |
| figure() | |
| for boost in range(-24, 0, 3): | |
| b, a = peaking(10, dBgain=boost, Q=sqrt(2), type=ftype, analog=True) | |
| w, h = freqs(b, a, 10000) | |
| plot(w, 20 * log10(abs(h)), 'r', alpha=0.5) | |
| for boost in range(0, 25, 3): | |
| b, a = peaking(10, dBgain=boost, Q=sqrt(2), type=ftype, analog=True) | |
| w, h = freqs(b, a, 10000) | |
| plot(w, 20 * log10(abs(h)), 'b', alpha=0.5) | |
| xscale('log') | |
| title(ftype + ' frequency response') | |
| xlim(0.1, 1000) | |
| xlabel('Frequency [radians / second]') | |
| ylabel('Amplitude [dB]') | |
| yticks(range(-24, 25, 3)) | |
| margins(0, 0.1) | |
| grid(True, color = '0.7', linestyle='-', which='major', axis='both') | |
| grid(True, color = '0.9', linestyle='-', which='minor', axis='both') | |
| show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment