Signal Processing/Filter Design

Filter design
Designing a filter generally starts with the specification of its frequency response. From this, both a transfer function and a filter structure have to be chosen. Indeed, some filter structures such as chains of second order sections are adaptable to any filter function whilst the basic ladder filter structure with only inductances on the series branches and capacitances on the parallel branches merely implement all-pole lowpass transfer functions. Mapping the transfer function to the filter structure gives the element values of analog filters elements and the coefficients of digital filters.

From the ideal filter, further optimizations have to be done.

The lumped element or the coefficient values are ideal values. Lumped element values have to be mapped to the ones of existing components. Filter coefficients have to be rounded to values supported by the number representation of the filter circuit. For example, choosing a fixed-point or a floating point Digital Signal Processor (DSP) will result in coarser or finer coefficient values.

Then, the signal amplitudes inside the filter have to be considered. For active filters and digital filters, the internal signals have to be checked against saturation or overflow. For any filter, the signals should not become too small, because this would seriously affect the signal to noise ratio of the whole filter. So basically, the filter design process doesn't only analyse the transfer function from the input to the output, but also the transfer function from the input to the internal signals.

Filter representations
A filter can be represented in many equivalent ways:
 * obviously by its transfer function
 * by the zeroes and poles of the transfer function, together with a gain factor
 * by the state-space representation of the filter

Methods allow to switch from one representation to the other. The different representations cast a specific light on the filter's properties. The transfer function, plotted for $$s = j\omega$$ shows how frequencies are amplified or damped. The zero-pole plot gives a direct information about the system stability. The state-space representation gives more insight into the chosen filter structure and allows to analyse the amplitudes of some signals (the state variables) inside the filter.

Transfer Function
The most common way to choose a transfer function for simple filters such as lowpass, highpass, bandpass or bandstop is to select a usual filter function: Butterworth, Chebyshev, Elliptic, Bessel…

Realizable transfer functions are in the form of the quotient of two polynomials.
 * $$H(s) = \frac{num(s)}{den(s)}$$

The filter order is the maximal power of $$s$$ (or of $$z$$ for sampled systems) in the two polynomials.

Most filter types try to get as close as possible to a brick-wall shape, with different optimization goals in mind. For these, the higher the filter order, the closer the function comes to the brick-wall shape. Butterworth filters do this by having a maximally flat (ripple-free) transfer funtion. Bessel filters, on the other hand, aim to minimize group delay variation in the passband.

These transfer functions are defined for normalized lowpass filters. Filter transformation techniques allow to transform the lowpass function to lowpass with a desired cutoff frequency, highpass, bandpass or bandstop.

Also, the functions stand for time-continuous systems, understood with the help of the Laplace analysis. For sampled systems, such as switched-capacitor or digital filters, a further step is needed: the $$s$$-based transfer function has to be transformed in a $$z$$-based transfer function.

Zeroes, Poles and Gain
The zeroes are the roots of the transfer function numerator. The poles are the roots of the transfer function denominator. The gain is the DC value of the transfer function: the value for $$s = 0$$.

The transfer function of a system given by its zeroes $$z_i$$, poles $$p_i$$ and gain $$k$$ is:
 * $$H(s) = k \cdot \frac{\prod_{i=1}^M (s-z_i)}{\prod_{i=1}^N (s-p_i)}$$

A stable LTI system has all its poles on the left side half-plane of $$s$$ (i.e., will have negative real parts).

State Space Representation
The state space representation describes the filter as a set of differential equations:
 * $$\dot{\mathbf{x}}(t) = A \cdot \mathbf{x}(t) + B \cdot u(t)$$
 * $$y(t) = C \cdot \mathbf{x}(t) + D \cdot u(t)$$

Here, the input is $$u(t)$$, the internal filter signals are $$\mathbf{x}(t)$$ and the output is $$y(t)$$. The Laplace transform of these equations is given by:
 * $$s \cdot \mathbf{X}(s) = A \cdot \mathbf{X}(s) + B \cdot U(s)$$
 * $$Y(s) = C \cdot \mathbf{X}(s) + D \cdot U(s)$$

The transfer function is given by:
 * $$H(s) = \frac{Y(s)}{U(s)} = C(s\mathbf{I} - A)^{-1}B + D$$

The eigenvalues of $$A$$ are the poles of the system. The characteristic polynomial of $$A$$ is the denominator of the transfer function.

4th Order Normalized Butterworth Filter
 This example goes through the different representations of a 4th order normalized Butterworth filter.

Contrarily to other filter functions given by polynomials, the Butterworth filter is best described by the location of its poles. It is an all-pole function: it has no zeroes.



Poles
The normalized Butterworth filter has equally spaced poles on the left half plane part of the unit circle. For a 4th order filter, their values are:

\begin{cases} p_1, p_1^* = e^{\pm i (\pi / 2 + \pi / 8)} \\ p_2, p_2^* = e^{\pm i (\pi / 2 + 3 \pi / 8)} \end{cases} or \begin{cases} p_1, p_1^* = -0.383 \pm 0.924 \cdot i \\ p_2, p_2^* = -0.924 \pm 0.383 \cdot i \end{cases} $$

Transfer Function
The transfer function is thus:
 * $$H(s) = \frac{1}{(s- p_1) \cdot (s- p_1^*) \cdot (s- p_2) \cdot (s- p_2^*) }$$
 * $$H(s) = \frac{1}{(s^2 - 2 \cdot re(p_1) \cdot s + 1) \cdot (s^2 - 2 \cdot re(p_2) \cdot s + 1) }$$
 * $$H(s) = \frac{1}{(s^2 + 0.765 s + 1) \cdot (s^2 + 1.848 s + 1) }$$
 * $$H(s) = \frac{1}{s^4 + 2.61 s^3 + 3.41 s^2 + 2.61 s +1}$$

State Space Representation
A possible state space representation is found from the transfer function:
 * $$Y(s) = U(s) \frac{1}{s^4 + 2.61 s^3 + 3.41 s^2 + 2.61 s +1}$$


 * $$s^4 \cdot Y(s) + 2.61 s^3 \cdot Y(s) + 3.41 s^2 \cdot Y(s) + 2.61 s \cdot Y(s) + Y(s) = U(s)$$
 * $$s \cdot s^3 \cdot Y(s) = - Y(s) - 2.61 s \cdot Y(s) - 3.41 s^2 \cdot Y(s) - 2.61 s^3 \cdot Y(s) + U(s)$$

Writing:
 * $$\begin{cases}

X_1(s) = Y(s) \\ X_2(s) = s \cdot Y(s) = s \cdot X_1(s) \\ X_3(s) = s^2 \cdot Y(s) = s \cdot X_2(s) \\ X_4(s) = s^3 \cdot Y(s) = s \cdot X_3(s) \\ \end{cases}$$

gives:
 * $$s \cdot X_4(s) = - X_1(s) - 2.61 \cdot X_2(s) - 3.41 \cdot X_3(s) - 2.61 \cdot X_4(s) + U(s) $$

This can be rewritten in matrix form as:
 * $$\begin{matrix}

s \cdot \begin{bmatrix} X_1 \\ X_2 \\ X_3 \\ X_4 \end{bmatrix} & = & \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ -1 & -2.61 & -3.41 & -2.61 \end{bmatrix} \begin{bmatrix} X_1 \\ X_2 \\ X_3 \\ X_4 \end{bmatrix} & + & \begin{bmatrix} 0 \\ 0 \\ 0 \\ 1 \end{bmatrix} U \\ Y & = & \begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} X_1 \\ X_2 \\ X_3 \\ X_4 \end{bmatrix} & + & \begin{bmatrix} 0 \end{bmatrix} U \end{matrix}$$

Octave Code
The following Octave code provides the same results and plots the amplitude response of the filter: > order = 4; > [z, p, g] = butter(order, 1, 's') z = [] (0x0) p = -0.38268 + 0.92388i -0.92388 + 0.38268i -0.92388 - 0.38268i -0.38268 - 0.92388i g = 1 > [b, a] = butter(order, 1, 's') b = 1 a = 1.0000  2.6131   3.4142   2.6131   1.0000 > w = logspace(-1, 1, 100); > while ( length(b) < length(a) ) >  b = [0, b]; > endwhile; > Adb = 20*log10(abs(freqs(b, a, w))); > figure(1); > semilogx(w, Adb); > grid; An amplitude response plot appears in a Gnuplot window > [A, B, C, D] = tf2ss(b, a)    A = 0.00000   1.00000   0.00000   0.00000 0.00000  0.00000   1.00000   0.00000         0.00000   0.00000   0.00000   1.00000        -1.00000  -2.61313  -3.41421  -2.61313     B = 0 0        0         1     C = 1   0   0   0 D = 0 > eig(A) ans = -0.38268 + 0.92388i -0.38268 - 0.92388i -0.92388 + 0.38268i -0.92388 - 0.38268i > poly(A) ans = 1.0000  2.6131   3.4142   2.6131   1.0000

Analog Lowpass Filter Design
 This example goes through the steps of designing an analog lowpass filter from given specifications.

Specification
A reduced version of CCITT G712 input filter specification, giving only the lowpass part, is shown in the plot on the side.

The passband goes up to 3 kHz and allows a maximal ripple of 0.125 dB. The stopband requires an attenuation of 14 dB at 4 kHz and an attenuation of 32 dB above 4.6 kHz.

Filter order
A first step is to evaluate the order needed for usual filter functions. As there is no specification on the phase, there is no need to consider a Bessel filter.

The following Octave code: fs = 40E3; fPass = 3000; rPass = 0.125; fStop1 = 4000; rStop1 = 14; fStop2 = 4600; rStop2 = 32; order1 = buttord(fPass/(fs/2), fStop1/(fs/2), rPass, rStop1); order2 = buttord(fPass/(fs/2), fStop2/(fs/2), rPass, rStop2); buttOrder = max(order1, order2) order1 = cheb1ord(fPass/(fs/2), fStop1/(fs/2), rPass, rStop1); order2 = cheb1ord(fPass/(fs/2), fStop2/(fs/2), rPass, rStop2); cheb1Order = max(order1, order2) order1 = cheb2ord(fPass/(fs/2), fStop1/(fs/2), rPass, rStop1); order2 = cheb2ord(fPass/(fs/2), fStop2/(fs/2), rPass, rStop2); cheb2Order = max(order1, order2) order1 = ellipord(fPass/(fs/2), fStop1/(fs/2), rPass, rStop1); order2 = ellipord(fPass/(fs/2), fStop2/(fs/2), rPass, rStop2); ellipticOrder = max(order1, order2) disp('Required filter orders:'); printf(" Butterworth      : %2i\n", buttOrder) printf(" Chebyshev type 1 : %2i\n", cheb1Order) printf(" Chebyshev type 2 : %2i\n", cheb2Order) printf(" Elliptic (Cauer) : %2i\n", ellipticOrder)

gives us the results: Required filter orders: Butterworth     : 13 Chebyshev type 1 : 6 Chebyshev type 2 : 6 Elliptic (Cauer) : 4

The Butterworth filter needs roughly twice as much hardware as the Chebyshev type 1. The Chebyshev type 2 and the elliptic filter have zeroes which the Chebyshev type 1 filter doesn't, as it is an all-pole function. So they too will need more hardware than the Chebyshev type 1.

Chebyshev type 1 will be the candidate chosen for this filter.

Transfer function
Octave does the frequency transforms automatically. The amplitude response is obtained with: pointNb = 100; wLog = logspace(2, 5, pointNb); AdbMin = 40; [order, wc] = cheb1ord(fPass/(fs/2), fStop2/(fs/2), rPass, rStop2); order = order + 1; wc = 1.03 * wc; [num, den] = cheby1(order, rPass, (fs/2)*wc, 's'); while ( length(num) < length(den) ) num = [0, num]; endwhile; Adb = 20*log10(abs(freqs(num, den, wLog))); Adb(Adb < -AdbMin) = -AdbMin; figure(1); semilogx(wLog, Adb); grid; hold on; semilogx([wLog(1), fPass, fPass], -[rPass, rPass, AdbMin], 'r'); semilogx([fStop1, fStop1, fStop2, fStop2, wLog(length(wLog))], ...        -[0    , rStop1, rStop1, rStop2, rStop2            ], 'r'); hold off; xlabel('frequency [Hz]'); ylabel('amplitude [dB]');

This gives the following plot:

As one might have noticed, the filter order has been pushed up to 7 and the cutoff frequency has been pushed a little higher in order to balance the distance to the specification borders on the passband and the stopband side.

Filter coefficients
Obviously, the filter coefficients depend on the chosen filter structure.

For a chain of second order sections, the higher order transfer function has to be split into the product of second order transfer functions: [zeroes, poles, gain] = cheby1(order, rPass, (fs/2)*wc, 's'); [SOS, G] = zp2sos(zeroes, poles, gain)

This gives: SOS = 1.0000e+00  0.0000e+00   0.0000e+00   1.0000e+00   2.0033e+03   3.0335e+06 1.0000e+00  0.0000e+00   0.0000e+00   1.0000e+00   1.3863e+03   7.0724e+06 1.0000e+00  0.0000e+00   0.0000e+00   1.0000e+00   4.9478e+02   1.0311e+07 1.0000e+00  0.0000e+00   0.0000e+00   1.0000e+00   1.1118e+03  -0.0000e+00 G = 2.4594e+23 - 1.3998e+07i

The "1, 0, 0" beginning of the SOS matrix lines shows that they are all-pole functions.

The last SOS matrix line with values 0 at the third and at the last positions shows that it corresponds to a first order section. Indeed, the 7th order function is made out of 3 second order sections and one first order section.

The imaginary part of the gain is to be ignored: it is certainly an error of fixed point calculus. The real part of the gain is equal to the product of the last (non-zero) elements of the SOS matrix lines.

From these coefficients, one can find the analog element values of state variable filters or Sallen-Key circuits which both are of all-pole kind. The coefficients can also used almost as such in the canonical biquads, but with $$b_1 = 0$$ and $$b_2 = 0$$, which simplifies the circuits.