Signals and Systems/Filter Design

From Wikibooks, open books for an open world
< Signals and Systems
Jump to: navigation, search

The Laplace transform has shown to allow to analyse 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 ampify 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[edit]

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

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

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 complex plane of 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[edit]

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

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

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

Devising the proper transfer function for a given filter function goes through the following steps:

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

Example: Lowpass Filter[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 Function[edit]

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

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

#-------------------------------------------------------------------------------
# Specifications
#
fs = 40E3;
fPass = 3000;
rPass = 0.125;
fStop1 = 4000;
rStop1 = 14;
fStop2 = 4600;
rStop2 = 32;

pointNb = 1000;
AdbMin = 40;

makeGraphics = 1;
figureIndex = 0;

#===============================================================================
# Normalized filter function
#
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

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

G712 butterworth normalized

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 fc = 0.87: 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[edit]

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.

#-------------------------------------------------------------------------------
# Denormalized filter function
#
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

G712 butterworth

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