Signal Processing/Filter Design

From Wikibooks, open books for an open world
Jump to: navigation, search

Filter design[edit]

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 sould 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[edit]

A filter can be represented in many equivalent ways:

  • obviously by its transfer function
  • by the zeroes and poles ot 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[edit]

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 polynoms.

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 polynoms.

A 4th order transfer function is in the form:

H(s) = \frac{b_4 s^{4} + b_3 s^{3} + b_2 s^{2} + b_1 s + b_0}{a_4 s^{4} + a_3 s^{3} + a_2 s^{2} + a_1 s + a_0}

A 4th order all-pole transfer function is in the form:

H(s) = \frac{1}{a_4 s^{4} + a_3 s^{3} + a_2 s^{2} + a_1 s + a_0}

A 4th order FIR transfer function is in the form:

H(z) = \frac{b_4 z^{4} + b_3 z^{3} + b_2 z^{2} + b_1 z + b_0}{1} = b_4 z^{4} + b_3 z^{3} + b_2 z^{2} + b_1 z + b_0

The different filter functions try to get as close as possible to a brick-wall shape, with different optimization goals in mind. The higher the filter order, the closest the function comes to the brick-wall shape. Butterworth filter do this by having a maximally flat (ripple-free) transfer funtion. Bessel functions put a special weight on keeping a constant group delay 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[edit]

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.

State Space Representation[edit]

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.

Examples[edit]

4th Order Normalized Butterworth Filter[edit]

Example
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 of a 4th order Butterworth filter

Poles[edit]

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[edit]

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[edit]

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[edit]

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[edit]

Example
This example goes through the steps of designing an analog lowpass filter from given specifications.

Specification[edit]

CCITT G712 input lowpass filter 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[edit]

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 check for a Bessel function.

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[edit]

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:

G712-Chebyshev1.svg

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[edit]

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.