Signals and Systems/Filter Design

The Laplace transform allows analyzing the frequency response of circuits based on the differential equations of their capacitive and inductive components. Filter design starts with finding the proper transfer function in order to amplify selected parts of a signal and to damp other ones as a function of their frequency.

Choosing the proper filter structure and deriving the coefficient values is a further topic prensented in the wikibook Signal Processing which deals with the application of signal and systems.

Brick-wall filters
Separating signal from noise or different signals in the same transmission channel basing on their frequency content is best done with a brick-wall filter which shows full transmission in the passband and complete attenuation in the nearby stopbands, with abrupt transitions.

This can be done with the help of the Fourier transform which provides complete information of the frequency content of a given signal. Having calculated a Fourier transform, one can zero out unwanted frequency contents and calculate the inverse Fourier Transform, in order to provide the signal filtered with a brick-wall gauge.

The Fourier transform being given by:
 * $$\mathcal{F}(f(t)) = F(j\omega) = \int_{-\infty}^\infty f(t)e^{-j\omega t}dt$$

one finds out that the Fourier transform integral, with its infinite bounds, would have to be calculated from the day of the creation of our universe and all the way up to the day of its decay before the integral could have been fully calculated. And only then can the ideal brick-wall filtered signal be delivered.

In more technical terms, the ideal brick-wall filter suffers from an infinite latency.

Analog filters
The analysis of analog circuits shows that their outputs are related to their input by a set of differential equations. The Laplace transform rewrites these differential equations as a set of linear equations of the complex variable $$s$$. With this, a polynomial function multiplying the Laplace transform of the input signal can be equated to another polynomial function multiplying the Laplace transform of the ouput signal:
 * $$(b_m s^{m} + b_{m-1} s^{m-1} + \ldots + b_1 s + b_0) \cdot X(s) = (a_n s^{n} + a_{n-1} s^{n-1} + \ldots + a_1 s + a_0) \cdot Y(s)$$

Thus, the transfer function of a realizable analog filter can be written as the ratio of two polynomial functions of $$s$$:
 * $$H(s) = \frac{Y(s)}{X(s)} = \frac{b_m s^{m} + b_{m-1} s^{m-1} + \ldots + b_1 s + b_0}{a_n s^{n} + a_{n-1} s^{n-1} + \ldots + a_1 s + a_0}$$

Hence, the problem of analog filter design is to find a pair of polynomial functions which, put together, best approximate the ideal but not realizable brick-wall transfer function.

In the early days of electric signal processing, scientists have come up with filter functions which are still largely used today. The functions they have devised are all of lowpass type. Frequency transformation techniques allow to find polynomials for other filter types such as highpass and bandpass.

The Complex Plane
The transfer function of an analog filter is the ratio of two polynomial functions of $$s$$:
 * $$H(s) = \frac{Y(s)}{X(s)} = \frac{num(s)}{den(s)}$$

The variable $$s$$ is a complex number which can be written as $$s = \sigma + i \cdot \omega$$. The complex plane is a plane with the imaginary axis vertical and the horizontal axis as the real part.

The roots of the transfer function numerator polynom are called the transfer function zeroes. The roots of the transfer function denominator polynom are called the transfer function poles.

The transfer function can be written as a function of its zeroes $$z_i$$, its poles $$p_i$$ and an additional gain factor $$k$$ in the form:
 * $$H(s) = k \cdot \frac{\prod_{i=1}^M (s-z_i)}{\prod_{i=1}^N (s-p_i)}$$

The poles and the zeroes of a transfer function can be drawn in the complex plane. Their position provide information about the frequency response of the system. Indeed, the frequency response is equal to the transfer function taken for $$s = i \omega$$, which is along the imaginary axis.
 * $$\frac{Y(i \omega)}{X(i \omega)} = H(s = i \omega)$$

Effect of Poles
A stable LTI system has all its poles on the left side half plane of $$s$$.

If a pole would be located on the imaginary axis, at $$p = i \omega_p$$, then the factor $$1/{(s-p)}$$ of the transfer function would be infinite at the point $$s = i \omega_p$$ and so would the global frequency response $$H(s = i \omega_p)$$. A special case of this is the integrator: it has a pole at $$p = 0$$, and indeed has a long-term infinite output for a constant, non-zero, input.

For poles close to the imaginary axis, the frequency response takes a large amplitude for frequencies close to them. In other words, poles close to the imaginary axis indicate the passband.

Effect of Zeros
If a zero is located on the imaginary axis, at $$z = i \omega_z$$, then the factor $$(s-z)$$ of the transfer function is zero at the point $$s = i \omega_z$$ and so is the global frequency response $$H(s = i \omega_z)$$.

Zeroes on, or close to the imaginary axis indicate the stopband.

Designing Filters
Devising the proper transfer function for a given filter function goes through the following steps:
 * selecting a normalized filter function
 * transforming and scaling the function for the particular needs

The coefficients of the numerator and denominator coefficients are finally used to calculate the element values of a selected filter circuit.



Example: Lowpass Filter
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 Function
As a first step, we have to choose a filter function.

Programs such as Octave or Matlab provide functions which allow to determine the minimal filter order required to fulfill a given specification. This is a good help when choosing from the possible functions.

Let's however here arbitrarily choose a Butterworth transfer function.

Normalized Filter Function
The following Octave script allows to plot the amplitudes of normalized Butterworth transfer functions from order 8 to 16.

#--- # fs = 40E3; fPass = 3000; rPass = 0.125; fStop1 = 4000; rStop1 = 14; fStop2 = 4600; rStop2 = 32; pointNb = 1000; AdbMin = 40; makeGraphics = 1; figureIndex = 0; # wLog = 2*pi*logspace(-1, 1, pointNb); fc = 0.87; Adb = []; for order = 8:16 [num, den] = butter(order, 2*pi, 's'); while ( length(num) < length(den) ) num = [0, num]; endwhile; Adb = [Adb; 20*log10(abs(freqs(num, den, wLog)))]; endfor Adb(Adb < -AdbMin) = -AdbMin; figureIndex = figureIndex+1; figure(figureIndex); semilogx(wLog/(2*pi), Adb); hold on; semilogx([wLog(1)/(2*pi), fc, fc], -[rPass, rPass, AdbMin], 'r'); semilogx([fStop1*fc/fPass, fStop1*fc/fPass, fStop2*fc/fPass, fStop2*fc/fPass, wLog(length(wLog))/(2*pi)], ...        -[0    , rStop1, rStop1, rStop2, rStop2            ], 'r'); hold off; axis([wLog(1)/(2*pi), wLog(length(wLog))/(2*pi), -AdbMin, 0]); grid; xlabel('frequency [Hz]'); ylabel('amplitude [dB]'); if (makeGraphics != 0) print -dsvg g712_butterworth_normalized.svg endif
 * 1) Specifications
 * 1) Normalized filter function
 * 1) Normalized filter function

The following figure shows the result: one needs at least a 13th order Butterworth filter to meet the specifications.



On the graph, one can note that all the amplitude responses go through the same point at -3 dB.

The specification frequencies have been scaled down to fit to the normalized cutoff frequency of 1 Hz. In the script, one might have noted an additional scaling factor of : this is due to the fact that the corner cutoff amplitude is -0.125 dB and not -3 dB. That value has been adjusted by hand for this example. Again, Octave or Matlab scripts automate this task.

Denormalized Filter Function
The frequency scaling of the normalized transfer function is done by replacing

$$s \to f_c \cdot s $$

The following Octave script does this by multiplying the numerator and denominator coefficients by the appropriate power of $$f_c$$.

#--- # order = 13; wLog = 2*pi*logspace(2, 5, pointNb); fc = 0.87; [num, den] = butter(order, 2*pi, 's'); while ( length(num) < length(den) ) num = [0, num]; endwhile; for index = 1:order+1 num(index) = num(index) * (fPass/fc)^(index-1); den(index) = den(index) * (fPass/fc)^(index-1); endfor Adb = 20*log10(abs(freqs(num, den, wLog))); Adb(Adb < -AdbMin) = -AdbMin; figureIndex = figureIndex+1; figure(figureIndex); semilogx(wLog/(2*pi), Adb); hold on; semilogx([wLog(1)/(2*pi), fPass, fPass], -[rPass, rPass, AdbMin], 'r'); semilogx([fStop1, fStop1, fStop2, fStop2, wLog(length(wLog))/(2*pi)], ...        -[0    , rStop1, rStop1, rStop2, rStop2            ], 'r'); hold off; axis([wLog(1)/(2*pi), wLog(length(wLog))/(2*pi), -AdbMin, 0]); grid; xlabel('frequency [Hz]'); ylabel('amplitude [dB]'); if (makeGraphics != 0) print -dsvg g712_butterworth.svg endif
 * 1) Denormalized filter function



The coefficients of the numerator and denominator coefficients are now ready to be used to calculate the element values of a selected filter circuit.