Digital Filters

I present some basic material on digital filters, along with interactive plots to visualize their properties and help in filter design. The compilation is essentially my reference book, but maybe you’ll find it useful too.

1. Basics
1.1. Basic Relations
1.2. Biquads
2. First-order Low Pass Filter
2.1. Continuous Time
2.2. Backward Difference, LPF1  (with interactive plot)
2.3. Bilinear, LPF2  (with interactive plot)
3. CIC Filter, Average  (with interactive plot)
4. Second-order Low Pass Filter
4.1. Continuous Time  (with interactive plot)
4.2. Backward Difference, LPF3
4.3. Bilinear, LPF4  (with interactive plot)
4.4. Cascaded First-order Low Pass Filters, LPF5
4.5. Butterworth, LPF6
4.6. Biquad
5. Second-order Inverse Chebyshev Low Pass Filter
5.1. Continuous Time
5.2. Biquad, LPF7  (with interactive plot)
6. Second-order Peaking and Notch Filter
6.1. Continuous Time
6.2. Biquad
7. Design Tools
7.1. CIC + Biquad Low Pass Filter
7.2. CIC + Inverse Chebyshev Biquad Low Pass Filter
Appendix


1. Basics [1|2|3|4|5|6|7|A]

1.1. Basic Relations
For a continuous-time transfer function H_s(s) the frequency response H^{(s)}_\omega(\omega) is calculated by substituting s = j \omega, where \omega is the real frequency. The frequency response H^{(z)}_\omega(\omega) of a discrete-time transfer function H_z(z) is calculated by substituting

z = \exp(j \omega \Delta t) ,

where \omega is the real frequency, and \Delta t the sampling time step. The sampling rate is

f_{sampling} = \dfrac{1}{\Delta t} .

Whenever the domain is obvious from the context, the super- and/or subscripts to the transfer function are dropped, and H(\omega), H(s) and H(z) written for simplicity.

The response of the discrete-time transfer function at zero frequency is given by H_\omega(0) = H_z(1), and that at the Nyquist frequency f_{Ny} = \frac{1}{2} f_{sampling} by H_\omega(\omega_{Ny}) = H_z(-1).

In order to derive H(z) from a continuous-time transfer function H(s), usually the variable s is substituted by different approximations. Typical choices are:

backward difference: s = \dfrac{1}{\Delta t} \left( 1 - z^{-1} \right) , Eq. (1.2)    
bilinear transformation, expansion of ln(z): s = \dfrac{2}{\Delta t} \dfrac{ 1 - z^{-1} }{  1 + z^{-1} } , Eq. (1.3)

These approximations represent non-linear mappings from continuous-time to discrete-time. For the bilinear transform, the effect is completely captured by warping the frequency axis according to

\omega^{(z)} = 2 \arctan\left( \dfrac{1}{2} \omega^{(s)} \Delta t \right) f_{sampling} ,

where \omega^{(s)} and \omega^{(z)} are the continuous-time and discrete-time frequencies, respectively.

The phase delay t_{PD}(\omega) is defined as

t_{PD}(\omega) =  -\dfrac{\phi(\omega)}{\omega} .

Its meaning is clarified by the relation \cos( \omega t + \phi ) = \cos[ \omega( t - t_{PD})]. For a linear-phase filter holds \phi(\omega) = -\tau \omega, corresponding to a constant phase delay of t_{PD}(\omega) = \tau.

The group delay t_{GD}(\omega) is defined as

t_{GD}(\omega) =  -\dfrac{d \phi(\omega)}{d \omega} .

1.2. Biquads
A common implementation of second-order filters uses a biquad. The transfer function is

H(z) = \dfrac{ a_0 + a_1 z^{-1} + a_2 z^{-2} }{ 1 + b_1 z^{-1} + b_2 z^{-2} } ,

and the difference equation reads

y_n = -b_1 y_{n-1} -b_2 y_{n-2} + a_0 x_n + a_1 x_{n-1} + a_2 x_{n-2} .

The biquad is quite general and can produce many different filters. For a low-pass filter, with unit zero-frequency gain or H_\omega(0) = H_z(1) = 1, the coefficients are restricted to a_0 + a_1 + a_2 = 1 + b_1 + b_2.

Often the bilinear transform is used to generate the biquad coefficients from a continuous-time transfer function H(s). The continuous-time frequency is typically normalized to some characteristic time T, or characteristic frequency \Omega = 1 / T, i.e., H(s) is written in terms of the variable Ts. One finds:

H(s) = \dfrac{ A_0 + A_1 (Ts)^1 + A_2 (Ts)^2 }{ B_0 + B_1 (Ts)^1 + B_2 (Ts)^2 } ,

H(z) = \dfrac{ \left( A_0 \Delta t^2 + 2 A_1 T \Delta t + 4 A_2 T^2 \right) + 2 \left( A_0 \Delta t^2 - 4 A_2 T^2 \right) z^{-1} + \left( A_0 \Delta t^2 - 2 A_1 T \Delta t + 4 A_2 T^2 \right) z^{-2} }{ \left( B_0 \Delta t^2 + 2 B_1 T \Delta t + 4 B_2 T^2 \right) + 2 \left( B_0 \Delta t^2 - 4 B_2 T^2 \right) z^{-1} + \left( B_0 \Delta t^2 - 2 B_1 T \Delta t + 4 B_2 T^2 \right) z^{-2} } .

The coefficients in H(z) should be normalized by C = B_0 \Delta t^2 + 2 B_1 T \Delta t + 4 B_2 T^2 to yield the biquad form used in this article. The value of T may be corrected with the warping formula, to match the response of the discrete-time filter with the response of the continuous-time filter at e.g. the characteristic frequency: T' = \frac{1}{2} \Delta t /  \tan\left( \frac{1}{2} \Omega \Delta t \right).

The code for implementing biquads is quite straight forward.

Biquad Direct Form I (preferred for fixed point):

  datatype y, num;
  num = a0 * x + a1 * regX1 + a2 * regX2;
  y = num - b1 * regY1 - b2 * regY2;
  regX2 = regX1;
  regX1 = x;
  regY2 = regY1;
  regY1 = y;
  return y; //the return value usually will have fewer bits than datatype

Biquad Direct Form II (preferred for floating point):

  datatype y;
  reg0 = x - b1 * reg1 - b2 * reg2;
  y = a0 * reg0 + a1 * reg1 + a2 * reg2;
  reg2 = reg1;
  reg1 = reg0;
  return y;

Biquad Direct Form II Transposed (supposedly preferred most for floating point):

  datatype y;
  y = a0 * x + reg1;
  reg1 = a1 * x - b1 * y + reg2;
  reg2 = a2 * x - b2 * y;
  return y;

In some cases the coefficients can be represented by integers resulting in more efficient fixed point code. This is filter dependent, and thus discussed in the corresponding chapters.

References I like:


2. First-order Low Pass Filter [1|2|3|4|5|6|7|A]

2.1. First-order Continuous Time
The continuous-time transfer function of a first-order low-pass filter (LPF) reads:

H(s) = \dfrac{1}{1 + T s} .

The amplitude is

|H(\omega)| = \dfrac{1}{\sqrt{1 + \left(T \omega \right)^2}} ,

and the cut-off frequency is

f_g = \dfrac{1}{2 \pi T} .

By using either the backward difference or the bilinear transformation for substituting s, two different discrete-time LPFs are derived.

2.2. First-order Backward Difference, LPF1
With a suitable parameterization, the discrete-time transfer function is obtained as

H(z) = \dfrac{\alpha}{1 - \left(1-\alpha\right) z^{-1}} ,

\alpha = \dfrac{\Delta t}{T + \Delta t} .

The difference equation becomes

y_n = y_{n-1} + \alpha \left( x_n - y_{n-1} \right) ,

which is probably the most frequently used version of a first-order LPF.

The amplitude frequency dependence is given by

|H(\omega)| = \dfrac{ \alpha }{ \sqrt{ 1 + \left(1-\alpha\right)^2 - 2 \left(1-\alpha\right) \cos(\omega \Delta t) } } .

It exhibits neither a zero nor a pole in the entire frequency range. The cutoff frequency f^{(s)}_g of the corresponding continuous-time filter is

f^{(s)}_g = \dfrac{1}{2\pi} \dfrac{ \alpha }{ 1 - \alpha } f_{sampling} .

However, the frequency axis is warped, and the -3 dB cutoff frequency f_g of the discrete-time filter is shifted to lower values, according to

f_g = \dfrac{1}{\pi} \arcsin\left( \dfrac{\alpha}{2\sqrt{1-\alpha}} \right) f_{sampling} \approx \dfrac{1}{2\pi} \dfrac{\alpha}{\sqrt{1-\alpha}} f_{sampling} .

The filter constant maybe chosen such that it is represented exactly in integer arithmetic, that is

\alpha = \dfrac{1}{N}

where N is an integer with N>1. The filter is then represented by the equations

f_n = f_{n-1} - f_{n-1}/N + x_n ,

y_n = f_n / N .

The following settings can be realized:

N T f^{(s)}_g f_g
\Delta t f_{s} f_{s}
f_s Hz ms Hz Hz

   log   phase
First-order Backward Difference LPF: Discrete-time first-order LPF derived from backward difference, for integer inverse filter constant.

2.3. First-order Bilinear, LPF2
With a parameterization ensuring H(1)=1 and H(-1) = 0, the discrete-time transfer function is obtained as

H(z) = \dfrac{ \dfrac{1}{2} \alpha \left( 1 + z^{-1} \right) }{ 1 - \left(1-\alpha \right) z^{-1} } ,

\alpha = \dfrac{2 \Delta t}{2 T + \Delta t} .

The difference equation becomes

y_n = y_{n-1} - \alpha y_{n-1} + \dfrac{1}{2} \alpha \left( x_n + x_{n-1} \right) .

The amplitude frequency dependence is given by

|H(\omega)| = \sqrt{\dfrac{ \dfrac{1}{2} \alpha^2 \left(1 + \cos(\omega \Delta t) \right) }{  1 + \left(1-\alpha\right)^2 - 2 \left(1-\alpha\right) \cos(\omega \Delta t) }}.

It exhibits a zero at \omega \Delta t = \pi or the Nyquist frequency f_{Ny} = \frac{1}{2} f_{sampling}. It has no poles in the entire frequency range. The cutoff frequency f^{(s)}_g of the corresponding continuous-time filter is

f^{(s)}_g = \dfrac{1}{2 \pi} \dfrac{ 2 \alpha }{ 2 - \alpha } f_{sampling} .

However, the frequency axis is warped, and the -3 dB cutoff frequency f_g of the discrete-time filter is shifted to lower values, according to

f_g = \dfrac{1}{\pi} \arctan\left( \dfrac{ \alpha }{ 2-\alpha } \right) f_{sampling} .

Note: This can be derived directly from the above equations, or from using the bilinear transform warping formula given in Chapter 1.1.

It is often convenient to chose the filter constant such that it is represented exactly in integer arithmetic. Two options emerge depending on how the 1/2 factor is handled. For instance:

\dfrac{\alpha}{2} = \dfrac{1}{N} ,

where N is an integer with N>2. The filter is then represented by the equations

f_n = f_{n-1} - 2 f_{n-1}/N + x_n + x_{n-1} ,

y_n = f_n / N .

The following settings can be realized:

N T f^{(s)}_g f_g
\Delta t f_{s} f_{s}
f_s Hz ms Hz Hz

   log   phase
First-order Bilinear LPF: Discrete-time first-order LPF derived from bilinear transformation, for integer inverse filter constant.


3. CIC Filter, Average [1|2|3|4|5|6|7|A]

Up to a normalization constant, the (one-stage) CIC and moving average filters are equivalent. They are FIR filters with a low-pass characteristic, and have no continuous-time counterpart.

The discrete-time transfer function of a moving average reads:

H(z) = \dfrac{1}{M} \sum\limits_{m=0}^{M-1} z^{-m} = \dfrac{1}{M} \dfrac{1 - z^{-M}}{1 - z^{-1}} ,

where M is the number of samples averaged. The difference equation reads

y_n = \dfrac{1}{M} \sum\limits_{m=0}^{M-1} x_{n-m} = y_{n-1} + \dfrac{1}{M} \left( x_n - x_{n-M} \right) ,

where the last expression points to the recursive calculation of the moving average.

The amplitude frequency dependence is

|H(\omega)| = \dfrac{1}{M} \sqrt{\dfrac{ 1 - \cos(M \omega \Delta t) }{ 1 - \cos(\omega \Delta t)  }} = \dfrac{1}{M} \left| \dfrac{ \sin\left(\dfrac{1}{2} M \omega \Delta t\right) }{ \sin\left( \dfrac{1}{2}\omega \Delta t\right)  } \right| .

It exhibits zeros at M \omega \Delta t = 2 \pi k or f = \dfrac{k}{M} f_{sampling}, with k integer and k \neq 0,M,2M,.... It has no poles. The frequency of the first zero is

f_0 = \dfrac{1}{M} f_{sampling} ,

and represents the passband edge. The moving average filter has quite significant attenuation in the passband, and the -3 dB frequency f_g lies below half of the frequency of the first zero:

f_g = 0.5 f_0     for M = 2 ,
f_g = 0.44295 f_0   for M \to \infty .

The impulse response of the moving average is symmetric, and the filter is therefore a linear-phase filter, with a delay of (M-1)/2 samples or t_{D} = \frac{1}{2} (M-1) \Delta t. The phase at the -3 dB frequency f_g is thus roughly -\pi/2, and at the first zero f_0 equal to -(1-1/M)\pi.

References I like:
M f_g f_0
f_s f_s
f_s Hz Hz Hz

   log   phase
CIC Filter: Discrete-time moving average filter, for averaging length M.


4. Second-order Low Pass Filter [1|2|3|4|5|6|7|A]

4.1. Second-order Continuous Time
The continuous-time transfer function of a second-order low-pass filter can be written as

H(s) = \dfrac{1}{1 + 2 d T s + T^2 s^2} ,

but many alternative parameterizations are in use. The characteristic frequency is

f_c = \dfrac{1}{2 \pi T} ,

and d is the damping parameter, where d>1 corresponds to an over-damped and d<1 to an under-damped system. Accordingly, for d \ge 1 the transfer function exhibits two real poles, and for d<1 a complex conjugated pair of poles. The quality factor is given by Q = 1 / (2d).

The amplitude and phase are

|H(\omega)| = \dfrac{1}{\sqrt{ 1 + 2 \left( 2 d^2 - 1 \right) (T \omega)^2 + (T \omega)^4 }} ,    \varphi(\omega) = \arctan\left( \dfrac{2 d T\omega}{ 1 - (T \omega)^2} \right) ,

and the cut-off frequency is

f_g = \sqrt{ \sqrt{ 1 + \left(1 - 2d^2\right)^2 } +  \left(1 - 2d^2\right) }\, f_c .

For d < 1/\sqrt{2}, the amplitude exceeds unity in the pass band, with a maximum at

f_{max} = \sqrt{ 1 - 2d^2 }\, f_c    |H|_{max} = \dfrac{1}{ 2d\sqrt{1 - d^2} } ,

and passes again through unity gain at the frequency

f_1 = \sqrt{ 2 \left(1 - 2d^2\right) }\, f_c .

d f_g f_{max} |H|_{max} f_1
1 0.6436 f_c
1/\sqrt{2} f_c
1/2 1.2720 f_c 0.7071 f_c 1.249 \mbox{ dB} f_c
f_c f_c dB f_c
f_c Hz Hz Hz Hz

   log   phase
Continuous Time Second-order LPF: Continuous-time second-order LPF, for specified damping d, in relation to the characteristic frequency f_c.

4.2. Second-order Backward Difference, LPF3
The discrete-time transfer function assumes the following structure

H(z) = \dfrac{\alpha}{1 - \left(1 - \alpha + \beta \right) z^{-1} + \beta z^{-2} } ,

with two parameters \alpha and \beta, which are related to d and T as

\alpha = \dfrac{ \Delta t^2 }{ T^2 +2dT\Delta t + \Delta t^2 } ,    \beta = \dfrac{ T^2 }{ T^2 +2dT\Delta t + \Delta t^2 } ,

The difference equation becomes

y_n = y_{n-1} + \beta \left( y_{n-1} - y_{n-2} \right) + \alpha \left( x_n - y_{n-1} \right) .

4.3. Second-order Bilinear, LPF4
With a paramterization ensuring H(1)=1 and H(-1) = 0, the discrete-time transfer function assumes the following structure

H(z) = \dfrac{ \dfrac{1}{4} \alpha \left( 1 + 2 z^{-1} + z^{-2} \right) }{ 1 - \left(2 - \alpha - \beta \right) z^{-1} + \left(1 - \beta \right) z^{-2} } ,

with two parameters \alpha and \beta, which are related to d and T as

\alpha = \dfrac{ 4 \Delta t^2 }{ 4 T^2 + 4 d T \Delta t + \Delta t^2 } ,     \beta = \dfrac{ 8 d T \Delta t  }{ 4 T^2 + 4 d T \Delta t + \Delta t^2 } ,     \dfrac{\alpha}{\beta} = \dfrac{1}{2 d} \dfrac{ \Delta t }{ T } ,

or then inverted

d = \dfrac{1}{2} \dfrac{\beta}{\alpha} \dfrac{\Delta t}{T} ,    \dfrac{\Delta t}{T} = \sqrt{\dfrac{ 4 \alpha  }{ 4 - 2 \beta - \alpha }} .

The difference equation becomes

y_n = y_{n-1} - \alpha y_{n-1} + \left(1 - \beta\right) \left( y_{n-1} - y_{n-2} \right) + \dfrac{1}{4} \alpha \left( x_n + 2 x_{n-1} + x_{n-2} \right) .

The amplitude frequency dependence exhibits a zero at \omega \Delta t = \pi or the Nyquist frequency f_{Ny} = \frac{1}{2} f_{sampling}. The frequency axis is warped, and the -3 dB cutoff frequency f_g of the discrete-time filter is shifted to lower values, according to

f_g = \dfrac{1}{\pi} \arctan( \pi f^{(s)}_g \Delta t) f_{sampling} ,

where f^{(s)}_g is the cutoff frequency of the continuous time filter given in Chapter 4.1. Analogous relations hold for the the frequencies of the maximum and the unity gain crossing, f_{max} and f_1. The height of the maximum is not affected by the warping, and is identical to |H|_{max}.

d T f_c f^{(s)}_g f_g
\Delta t f_{s} f_{s} f_{s}
f_s Hz ms Hz Hz Hz
f_{max} = Hz    |H|_{max} = dB    f_1 = Hz

   log   phase
Second-order Bilinear LPF: Discrete-time second-order LPF derived from bilinear transformation.

4.4. Cascaded First-order Low Pass Filters, LPF5
For d = 1, the second-order LPF reduces to a cascade of two identical first-order LPFs, each with time constant T.

Only the bilinear filter is considered. The constants \alpha and \beta reduce to \alpha = \alpha_1^2 and \beta = (1- \alpha_1)^2, where \alpha_1 is the \alpha-constant of the first-order LPF (LPF2), and the discrete-time transfer function becomes

H(z) = \left( \dfrac{ \dfrac{1}{2} \alpha_1 \left( 1 + z^{-1} \right) }{ 1 - \left(1 - \alpha_1\right) z^{-1}  } \right)^2 .

This leads to a straight forward and cheap implementation. The frequency response is however not very desirable. The crossover region is „smeared“ out, and the cutoff frequency (of the continuous-time transfer function) is shifted down to f_g = 0.6436 f_c, see also the table given in Chapter 4.1.

4.5. Second-order Butterworth, LPF6
For d = 1/\sqrt{2}, the most flat passband together with narrowest transition range is achieved. The continuous-time cutoff frequency matches the characteristic frequency, f^{(s)}_g = f_c. Commonly the filter is implemented using a biquad.

One may chose the filter constant such that it is represented exactly in integer arithmetic. For instance:

\dfrac{\alpha}{4} = \dfrac{1}{N}

where N is an integer with N>4. The filter is then represented by the equations

f_n = f_{n-1} - 4 f_{n-1}/N +  \left(1-\beta\right) \left( f_{n-1} - f_{n-2} \right)  +  x_n + 2 x_{n-1} + x_{n-2},

y_n = f_n / N.

The multiplication with (1-\beta) can be executed approximately but efficiently by choosing two integers X and M, such that M/X \approx (1-\beta), and thus

f_n = f_{n-1} - \left[ 4X f_{n-1} -  M N \left( f_{n-1} - f_{n-2} \right) \right]/(XN)  +  x_n + 2 x_{n-1} + x_{n-2}.

N X M T d f^{(s)}_g f_g
5 16 3 0.7016 \Delta t 0.7238 0.2269 f_s 0.1971 f_s
6 32 7 0.8149 \Delta t 0.7190 0.1953 f_s 0.1752 f_s
7 4 1 0.9186 \Delta t 0.7144 0.1733 f_s 0.1587 f_s
8 32 9 1.0155 \Delta t 0.7078 0.1567 f_s 0.1456 f_s
13 8 3 1.4087 \Delta t 0.7209 0.1130 f_s 0.1086 f_s
16 64 27 1.6105 \Delta t 0.7179 0.0988 f_s 0.0958 f_s
32 16 9 2.4450 \Delta t 0.7144 0.06497 f_s 0.06409 f_s

An interesting approximation is given by

N X M T d f^{(s)}_g f_g
8 4 1 1 \Delta t 0.75 0.1591 f_s 0.1476 f_s

which corresponds to the difference equations

f_n = f_{n-1} - \left( f_{n-1} + f_{n-2} \right)/4  +  x_n + 2 x_{n-1} + x_{n-2} ,

y_n = f_n / 8.

4.6. Biquad
For biquads see Chapter 1.2. For a low-pass filter, with unit zero-frequency gain and a zero at the Nyquist frequency, the biquad coefficients have to fulfill

a_0 = a_2 = a_1/2 ,
a_0 + a_1 + a_2 = 1 + b_1 + b_2 ,

yielding the parameterization in Chapter 4.3.


5. Second-order Inverse Chebyshev Low Pass Filter [1|2|3|4|5|6|7|A]

The inverse Chebyshev, or Chebyshev Type II, low-pass filter has some nice features. Especially, it provides a sharp transition region, which makes it e.g. a nice anti-aliasing filter, and at the same time retains a good phase behavior, comparable to a Butterworth filter. However, for even-order filters the passband attenuation stays finite at high frequencies, and odd-order filters are preferred from that perspective. In addition, for the discrete-time filter I could not find a single setting of interest which could be realized with integer-valued parameters, i.e., one needs to resort to a full biquad implementation, and pay attention to its pain-points when using fixed-point arithmetic.

5.1. Inverse Chebyshev Continuous Time
The amplitude frequency response of a second-order inverse Chebyshev low-pass filter is given by

|H(\omega)| = \dfrac{1}{ \sqrt{1 + \dfrac{1}{\epsilon^2 T_2^2(\omega_c/\omega)} } } ,

with the Chebyshev polynomial T_2(x) = 2 x^2 -1, the characteristic frequency f_c, and the ripple parameter \epsilon<1.

The -3 dB cut-off frequency f_g, the frequency location of zero amplitude f_0, and the gain in the stopband A are given by

f_g = \sqrt{\dfrac{2 \epsilon}{1 + \epsilon}}\, f_c ,    f_0 = \sqrt{2}\, f_c ,    A = \dfrac{ \epsilon}{\sqrt{1 + \epsilon^2}} .

The characteristic frequency f_c is generally considered as passband edge, and corresponds to the frequency where the amplitude falls below the stopband gain factor A the first time.

The transfer function H(s) is

H(s) = A \dfrac{ (s - \hat{z}_1)(s - \hat{z}_2) }{  (s - \hat{p}_1)(s - \hat{p}_2) } ,

with the zeros and poles located at

\hat{z}_{1,2} = \pm j \sqrt{2}\, \omega_c ,    \hat{p}_{1,2} = \left( -\sqrt{\dfrac{1-A}{2}} \pm j \sqrt{\dfrac{1+A}{2}} \right) \sqrt{2A}\, \omega_c ,

References I like:

5.2. Inverse Chebyshev Biquad, LPF7
The math involved with Chebyshev filters is somewhat lengthy. However, it can be bypassed by applying this set of conditions in the determination of the biquad coefficients:

  • (i) location of the zero-gain frequency f_0, implying H(\exp(-j\omega_0\Delta t)) = 0
  • (ii) gain A in the stopband, implying |H(-1)| = A
  • (iii) unity gain at zero frequency, implying H(1)=1
  • (iv) zeros are conjugated complex pairs on the unit circle
  • (v) poles are conjugated complex pairs within the unit circle

The filter coefficients are then uniquely determined by three parameters, the frequency of the zero f_0, the gain in the stopband A, and a third parameter which governs the behavior in the passband. The value of this third parameter would be fixed by requiring a Chebyshev characteristic. However, it is useful to keep it as a tuning parameter, which implies that the filter below is not strictly an inverse Chebyshev filter, but more general.

The discrete-time transfer function can be expressed as

H(z) = \dfrac{  \alpha \left( 1 - \hat{z}_0 z^{-1} \right)\left( 1 - \hat{z}_0^* z^{-1} \right) }{ \left( 1 - \hat{z}_p z^{-1} \right)\left( 1 - \hat{z}_p^* z^{-1} \right) } = \dfrac{  \alpha \left( 1 - 2 \cos\varphi_0 z^{-1} + z^{-2} \right) }{ 1 - 2 z_p \cos\varphi_p z^{-1} + z_p^2 z^{-2} } ,

with zeros at \hat{z}_0 = e^{j \varphi_0} and \hat{z}_0^*, and poles at \hat{z}_p = z_p e^{j \varphi_p} and \hat{z}_p^*, where z_p \le 1 . The frequency response is then

H(\omega) = \dfrac{  2 \alpha \left[ \cos(\omega \Delta t) - \cos(\varphi_0) \right] }{ (1 + z_p^2) \cos(\omega \Delta t) + j (1 - z_p^2) \sin(\omega \Delta t) - 2 z_p \cos(\varphi_p)  } .

The free parameters are chosen as the zero-gain frequency f_0, the stopband gain A, and the magnitude of the pole z_p. The remaining values are then given by

\varphi_0 = 2 \pi f_0 \Delta t ,

\alpha = \dfrac{ (1 + z_p^2) A }{ (1+A) + (1-A) \cos(\varphi_0) } ,

\cos(\varphi_p) = \dfrac{ (1-A) + (1+A) \cos\varphi_0  }{ (1+A) + (1-A) \cos\varphi_0  } \, \dfrac{ 1 + z_p^2 }{ 2 z_p } .

Inverse Chebyshev characteristic would be achieved by choosing z_p according to

|z_p|^2 = \dfrac{ (1+A) + (1-A) \cos\varphi_0 - \sqrt{2(1-A)A}\, \sin\varphi_0 }{ (1+A) + (1-A) \cos\varphi_0 + \sqrt{2(1-A)A}\, \sin\varphi_0 } .

A f_0 z_p z^{cheby}_p       
dB f_s
f_s Hz Hz
f_g = Hz    f_{max} = Hz    |H|_{max} = dB    f_1 = Hz

   log   phase
Inverse Chebyshev LPF: Discrete-time second-order inverse Chebyshev LPF based on a biquad with a generalized paramterization. Comment: The characteristic values are determined numerically.


6. Second-order Notch and Peaking Filter [1|2|3|4|5|6|7|A]

The notch filter is ideal for filtering out one specific frequency. While the notch filter has a zero at the center frequency f_0, a variant of it often refereed to as peaking filter allows us to set the gain (attenuation) at f_0. The notch filter is then recovered for infinite attenuation. The peaking filter is thus more versatile and considered here, but also often refered to as notch filter.

6.1. Continuous Time
The frequency response of a second-order peaking filter is given by

H(s) = \dfrac{1 + \dfrac{A}{Q} s + s^2}{1 + \dfrac{1}{Q} s + s^2}

The notch filter corresponds to A = 0. The following, equivalent formulations are also often found:

H(s) = \dfrac{1 + \dfrac{1}{Q'} s + s^2}{1 + \dfrac{1}{AQ'} s + s^2},    H(s) = \dfrac{1 + \dfrac{A'}{Q''} s + s^2}{1 + \dfrac{1}{A'Q''} s + s^2}

Note: The peak filter can be derived from a band-pass filter H_b(s) in several ways, e.g.: H(s) =  1 + V H_b(s) or H(s) = \left[1 - V H_b(s) \right]^{-1} or H(s) =  \left[ 1 + V H_b(s) \right] / \left[ 1 - V H_b(s) \right], corresponding to the three given forms. The meaning of Q changes thereby.

This gives rise to a multitude of equations which can be found on the web. In this article the first version is prefered, since here Q has the „usual“ meaning as it is known for the second-order low-pass, band-pass, and high-pass filters.

The definition of „bandwidth“ becomes somewhat fuzzy and allows for several possibilities, which further increases the number of equations found on the web. It is especially reminded that for the band-pass filter H_b(s) = (s/Q) /( 1 + s/Q + s^2), which underlies the peaking filter, it strictly holds

f_1 f_2 = f_0^2

whereas

\dfrac{f_2 - f_1}{f_0} = \dfrac{1}{Q}

holds only approximately, for large values of Q. Here, f_1 and f_2 are the lower and upper -3 dB frequencies, and the bandwidth would be defined as \delta f = f_2-f_1. For a peaking filter with gain this also would be an appropriate definition of bandwidth. However, for a peaking filter with attenuation the frequencies f_1 and f_2 do in general not correspond to the -3 dB points.

6.2. Biquad
The biquad coefficients can be derived by applying this set of conditions:

  • (i) unity gain at zero frequency, implying H(1)=1
  • (ii) unity gain at the Nyquist frequency, implying H(-1)=1
  • (iii) location of the center frequency f_0 and pass-through amplitude A, implying H(\exp(-j\omega_0\Delta t)) = A
  • (iv) bandwidth \delta f at the center frequency f_0

The filter coefficients are then uniquely determined by three parameters, the center frequency f_0, the pass-through amplitude A at the center frequency, and a third parameter which relates to the bandwidth.

The discrete-time transfer function can be expressed in many equivalent ways, which differ in how the biquad coefficients are parameterized. In fact, quite many different parameterizations can be found on the web, too many to be all discussed here. Most of these equations are obtained from discretizing the continous-time transfer function using the bilinear transformation.

A simple parameterization, which follows quite directly from imposing the above set of conditions, is

H(z) = \dfrac{  \left( 1 + \alpha A \right) - 2 \cos(\omega_0\Delta t) z^{-1} + \left( 1 - \alpha A \right) z^{-2} }{ \left( 1 + \alpha \right) - 2 \cos(\omega_0\Delta t) z^{-1} + \left( 1 - \alpha \right) z^{-2} }.

Here, the attenuation at the center frequency is directly specified by A. The parameter \alpha is related to the bandwidth (in a non-linear fashion). The equation is sometimes also written in the more „normalized“ fashion

H(z) = \left(\dfrac{1 + k \mu}{1 + k}\right) \dfrac{ 1 - \dfrac{2 \cos(\omega_0\Delta t)}{1 + k \mu} z^{-1} + \dfrac{1 - k \mu}{1 + k \mu} z^{-2} }{ 1 - \dfrac{2 \cos(\omega_0\Delta t)}{1 + k} z^{-1} + \dfrac{1 - k}{1 + k} z^{-2} }

with obvious relations between the \alpha, A and k, \mu parameters.

Another frequently found set of equations is obtained by pre-warping the bilinear transformation, i.e., by using s = \dfrac{1}{K} \dfrac{1-z^{-1}}{1+z^{-1}}, where K = \tan(\frac{1}{2}\omega_0 \Delta t). This results in

H(z) = \dfrac{  \left(1 + K^2 + \dfrac{K}{Q} A \right) + 2 \left( K^2 - 1 \right) z^{-1} + \left( 1 + K^2 - A \dfrac{K}{Q} \right) z^{-2} }{ \left( 1 + K^2 + \dfrac{K}{Q} \right) + 2 \left( K^2 - 1 \right) z^{-1} + \left( 1 + K^2 - \dfrac{K}{Q} \right) z^{-2} }

with obvious relations of its parameters to the above \alpha, A parameters. It should be noted however that the applied warping only matches the center-frequency of the discrete-time filter to that of the continuous-time filter, but not e.g. f_1 and f_2. Therefore, the frequency response of the discrete-time filter will in general be distorted, and its actual bandwidth be different from that of the continuous-time filter. In short, the „bandwidth“ used in this sort of equations is not a precisely defined aspect of the discrete-time filter’s frequency response but more of a figure-of-merit, which more or less represents the actual bandwidth.

Note: \dfrac{1-K^2}{1+K^2} = \cos(\omega_0 \Delta t), \dfrac{K}{1+K^2} = \dfrac{1}{2} \sin(\omega_0 \Delta t)

* http://www.eng.ucy.ac.cy/cpitris/courses/ECE623/presentations/DSP-LECT-20.pdf
* https://github.com/libaudioverse/libaudioverse/blob/master/audio%20eq%20cookbook.txt
* http://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
* https://pdfs.semanticscholar.org/56e6/1c03d338c293afa6cbb470b11d360195cf5f.pdf
* https://www.dummies.com/education/science/science-engineering/how-to-characterize-the-peaking-filter-for-an-audio-graphic-equalizer/
* https://ccrma.stanford.edu/~jos/fp/Peaking_Equalizers.html
* https://www.dsprelated.com/freebooks/filters/Peaking_Equalizers.html
* http://www.advsolned.com/allpass-peaking-bell-filter/
* bilinear pre-warped z transform: https://www.earlevel.com/main/2003/03/02/the-bilinear-z-transform/
* bandwidth in octaves: https://www.rane.com/note170.html
* http://dspguide.com/ch19/3.htm


7. Design Tools [1|2|3|4|5|6|7|A]

The combination of a one-stage CIC filter with a biquad can yield efficient and/or effective anti-aliasing filters, useful for instance for downsampling in real-time applications.

References I like:

7.1. CIC + Biquad Low Pass Filter

LPF order T d CIC order M
\Delta t
f_s Hz ms
f_g = Hz    f_{max} = Hz    |H|_{max} = dB    f_1 = Hz
\alpha = = 1/    \beta = = 1/

   log   phase
CIC + Biquad LPF: Discrete-time moving average plus first- or second-order LPF derived from bilinear transformation. Comment: The characteristic values are determined numerically.

Example settings:

M T d \alpha \beta f_g
12 3.20156 \Delta t 0.46852 1/12 1/4 0.04372 f_s
12 2.63629 \Delta t 0.30346 1/8 1/5 0.05581 f_s
3 0.59161 \Delta t 0.33806 5/4 1/2 0.02183 f_s

7.2. CIC + Inverse Chebyshev Biquad Low Pass Filter

A f_0 z_p CIC order M
dB f_s
f_s Hz Hz
f_g = Hz    f_{max} = Hz    |H|_{max} = dB    f_1 = Hz
\alpha = = 1/    a_1 = = 1/ = \alpha
b_1 = = 1/    b_2 = = 1/

   log   phase
CIC + Inverse Chebyshev LPF: Discrete-time moving average plus second-order inverse Chebyshev LPF based on a biquad usinga generalized parameterization. Comment: The characteristic values are determined numerically.


Appendix [1|2|3|4|5|6|7|A]

A.1. References

All given references, and some more:


A.2. Poles and Zeros

The continuous-time transfer function H(s) may be written in terms of its poles and zeros as

H(s) = A \dfrac{ \prod\limits_k \left( s - \hat{z}_k \right) }{ \prod\limits_k \left( s - \hat{p}_k \right) } ,

with a constant A, and the (complex valued) zeros \hat{z}_k and poles \hat{p}_k, written as

\hat{z}_k = z'_k + j z''_k ,    \hat{p}_k = p'_k + j p''_k .

The squared magnitude of the frequency-dependent transfer function H(\omega) = H_s(j\omega) then becomes

|H|^2(\omega) = H(\omega)H^*(\omega) = A^2 \dfrac{ \prod\limits_k \left[ \omega - (z''_k + j z'_k) \right]\left[ \omega - (z''_k - j z'_k) \right] }{ \prod\limits_k \left[ \omega - (p''_k + j p'_k) \right] \left[ \omega - (p''_k - j p'_k) \right] } ,

which has complex conjugated pairs of zeros at

\hat{Z}_k = z''_k + j z'_k = j \hat{z}_k^* ,    \hat{Z}_k^* = z''_k - j z'_k = -j \hat{z}_k ,

and complex conjugated pairs of poles at

\hat{P}_k = p''_k + j p'_k = j \hat{p}_k^* ,    \hat{P}_k^* = p''_k - j p'_k = -j \hat{p}_k .

From determining the zeros and poles of |H|^2(\omega) one thus can straightforwardly conclude the zeros and poles of H(s), selecting the poles in the left-half s plane, p'_k<0, for the stability of H(s).

Example: Second-order LPF
We consider a second-order LPF, with |H|^2(\omega) written as

|H|^2(\omega) = \dfrac{ \omega_0^4 }{ \omega^4 + 2(2d^2-1)\omega_0^2 \omega^2 +  \omega_0^4} .

Now d < 1/\sqrt{2} shall be assumed. The poles are then determined to

X^2 + 2(2d^2-1)\omega_0^2 X + \omega_0^4 = 0 ,

X_{1/2} = (1-2d^2) \omega_0^2 \pm 2d \omega_0^2 \sqrt{ d^2-1} = \left[ (1-2d^2)  \pm j 2 d \sqrt{ 1 - d^2} \right] \omega_0^2 ,

\omega_{1\pm} = \pm \sqrt{ X_1 } = \pm \sqrt{ (1-2d^2) + j 2 d \sqrt{1 - d^2} }\, \omega_0 = \pm \left( \sqrt{ 1-d^2 } + j d \right) \omega_0 ,

\omega_{2\pm} = \pm \sqrt{ X_2 } = \pm \sqrt{ (1-2d^2) - j 2 d \sqrt{1 - d^2} }\, \omega_0 = \pm \left( \sqrt{ 1-d^2 } - j d \right) \omega_0 ,

where in the last step the complex square roots were evaluated. This yields the poles

\hat{p} = \left( -d \pm j \sqrt{ 1-d^2 } \right) \omega_0 ,

with phases chosen such that the poles are in the left-half plane. One finally obtains

H(s) = \dfrac{ \omega_0^2 }{ (s - \hat{p})(s - \hat{p}^*) }   = \dfrac{ \omega_0^2 }{ s^2  + 2 d \omega_0 s + \omega_0^2 } .


A.3. Square Root of a Complex Number

\sqrt{\hat{z}} = \sqrt{x+jy} = \pm \left[ \sqrt{\dfrac{|\hat{z}| + x}{2}} + j \, \text{sign}(y) \sqrt{\dfrac{|\hat{z}| - x}{2}} \right] .


A.4. Bilinear Transformation

s = \dfrac{2}{\Delta t} \dfrac{z - 1}{z + 1} ,   z = \dfrac{ 1 + \frac{\Delta t}{2} s}{ 1 - \frac{\Delta t}{2} s} .

With s = x + j y one finds for z =|z| e^{j \varphi}:

z = \dfrac{ 1 - \left(\frac{\Delta t}{2}\right)^2 |s|^2 + j \Delta t y }        { 1 - \Delta t x + \left(\frac{\Delta t}{2}\right)^2 |s|^2 } = \dfrac{1 + |z|^2}{2}\cdot \dfrac{ 1 - \left(\frac{\Delta t}{2}\right)^2 |s|^2 + j \Delta t y}{ 1 + \left(\frac{\Delta t}{2}\right)^2 |s|^2 } .

\cos\varphi = \dfrac{1 + |z|^2}{2 |z|} \cdot \dfrac{ 1 - \left(\frac{\Delta t}{2}\right)^2 |s|^2}{ 1 + \left(\frac{\Delta t}{2}\right)^2 |s|^2 } = \dfrac{1 + |z|^2}{2 |z|}\, \cos\left[ 2 \arctan\left( \frac{\Delta t}{2} |s| \right) \right] .

|z|^2 = \dfrac{ 1 + \Delta t x + \left(\frac{\Delta t}{2}\right)^2 |s|^2 }         { 1 - \Delta t x + \left(\frac{\Delta t}{2}\right)^2 |s|^2 } ,    |z| = \cot\left[ \dfrac{1}{2} \arccos\left( \dfrac{ \Delta t x }{ 1 + \left(\frac{\Delta t}{2}\right)^2 |s|^2 }  \right) \right] .

Poles and zeros transform as expected:

H(s) = A \dfrac{ \prod\limits_k \left( s - \hat{z}_k \right) }{ \prod\limits_k \left( s - \hat{p}_k \right) } \quad\longrightarrow\quad H(z) = A' \dfrac{ \prod\limits_k \left( z - \bar{z}_k \right) }{ \prod\limits_k \left( z - \bar{p}_k \right) } ,

\hat{z}_k , \hat{p}_k \quad\longrightarrow\quad \bar{z}_k = \dfrac{ 1 + \frac{\Delta t}{2} \hat{z}_k}{ 1 - \frac{\Delta t}{2} \hat{z}_k} ,\; \bar{p}_k = \dfrac{ 1 + \frac{\Delta t}{2} \hat{p}_k}{ 1 - \frac{\Delta t}{2} \hat{p}_k} .

Example: Second-order Inverse Chebyshev LPF
Using the method of Chapter A.2., the zero and poles of H(s) are determined as

\hat{z}_{1,2} = \pm j \omega_0 ,    \hat{p}_{1,2} = \left( -\sqrt{\dfrac{1-A}{2}} \pm j \sqrt{\dfrac{1+A}{2}} \right) \sqrt{A}\, \omega_0 .

Using the notation of chapter 5.2, one finds the zeros of H(z) as

z^2_0 = 1,    \varphi_0 = 2 \arctan\left( \frac{\Delta t}{2} \omega_0 \right) .

Considering that \hat{z}_0 = e^{j\varphi_0} and that on the frequency axis z = e^{j \omega \Delta t}, one just has recovered the warping formula applied to \omega_0. For the poles one finds

z_p^2 = \dfrac{1 - \sqrt{2(1-A)A}\, X + A X^2}{1 + \sqrt{2(1-A)A}\, X + A X^2} ,    \cos\varphi_p = \dfrac{1 + z_p^2}{2 z_p} \cdot \dfrac{ 1 - A X^2}{ 1 + A X^2 },    with X = \dfrac{\Delta t}{2} \omega_0 .

Substituting \frac{\Delta t}{2} \omega_0 = \tan(\frac{\varphi_0}{2}) and exploiting the relation \tan^2(x/2) = (1-\cos x)/(1+\cos x) yields the result given in Chapter 5.2.


A.5. Frequency Response of a Biquad

The frequency response H(\omega) = H_z(e^{j \omega \Delta t}) of a biquad

H_z(z) = \dfrac{ a_0 + a_1 z^{-1} + a_2 z^{-2} }{ 1 + b_1 z^{-1} + b_2 z^{-2} }

is calculated as

H(\omega) = \dfrac{a_1 + \left(a_0+a_2\right) \cos(\omega \Delta t) + j \left(a_0-a_2\right) \sin(\omega \Delta t) }{ b_1 + \left(1+b_2\right) \cos(\omega \Delta t) + j \left(1-b_2\right) \sin(\omega \Delta t) }



Hinterlasse einen Kommentar