Chapter 2

Additional Information About Filtering

Mathematically speaking, using the filter with the transfer function H(z) for filtering the sequence x(n) leads to a perfectly determined result. However, depending on the practical implementation of the filter, the results can vary in terms of precision, speed, etc. This chapter deals with the technical aspects of filtering. If you restrict yourself to a "simulation" approach, as we have up until now, the filter function is everything you will ever need. However, if this filtering operation has to be implemented, its effectiveness requires some additional knowledge that will be detailed in this chapter.

2.1 Filter implementation

Linear invariant filters as we consider in this section perform the convolution of the input sequence {x(n)} with a sequence {h(n)} called the impulse response of the filter:

(2.1)

When the series {h(n)} is of finite length, we speak of a finite impulse response filter or FIR filter; otherwise, we have an infinite impulse response filter or IIR filter.

The DTFT of {h(n)} is called the complex gain of the filter, and its modulus is the gain.

The ZT of {h(n)} is called the transfer function.

2.1.1 Examples of filter structures

In this section, we will study the implementation of the filtering function, in other words its programming. Figure 2.1 shows a particular implementation called the canonical direct form of a general recursive filter with the transfer function:

Choosing the same degree for both the numerator and the denominator does not restrict us in any way; you need only consider that some of the coefficients can be equal to zero.

Figure 2.1 - Processing architecture

This "implementation" first performs the calculation of:

t(n) = i(n) - a1t(n - 1) - . - apt(n - p)

then the calculation of:

o(n) = b0t(n) + b1t(n - 1) + . + bpt(n - p)

where {i(n)} and {o(n)} are the input and output sequences respectively.

For this algorithm, the vector x(n) [x0 (n) x1 (n) ... xp-1 (n)]T (p× 1) is called the filter state:

x(n) = [t(n) t(n - 1) ... t(n - p + 1)]T

Its components, referred to as state variables, are the input values of the "delay" cells denoted z-1 in Figure 2.1. By introducing the vectors:

a = [a1 a2 ... ap]T and b = [b1 b2 ... bp]T

we get the following expression for the algorithm:

(2.2)

The following filtrer function implements this algorithm:

Determining the initial state leading to a given input-output sequence is another problem altogether. Using the recursive equations (2.2) that lead to t(n) and o(n), we can also write:

This expression shows that the initial state T = [t(-1) ... t(-p)]T can be reconstructed so long as the matrix is invertible. The system theory demonstrates that the possibility of reconstruction is related to the concept of observability. A great number of observability criteria exist, based on the state representations associated with a system [22].

The filtric function detailed hereafter carries out the reconstruction of the state associated with the processing architecture implemented by filter:

The state reconstruction function is inseparably related to the filtering function, which is itself based on a particular processing architecture.

The reconstruction function is rarely used. The state vector xs, final state of the first block's processing, is transmitted as the initial state of the second block. The result yp is identical to the one obtained for the filtering of the entire block etot.

MATLAB®'s filtering function, filter, uses the Transpose-Form IIR structure [28], different from the previous one, represented in Figure 2.2. As in our example, filter transmits the state vector. The reconstruction function filtic.m is available as part of the Signal Toolbox. Exercise 2.1 is a study of this structure.

Exercise 2.1 (Filter architecture) (see p. 189)

Consider the Transpose-Form IIR structure (Figure 2.2) of a rational filter.

Figure 2.2 - Transpose-Form IIR structure

1. Determine the filter's transfer function. 2. By defining the state x = [x1(n) ... xp(n)] at the time n, determine the state representation and express it as follows: Use this to find the filtering program. It might be useful to notice that the matrix A is the transpose of the companion matrix (compan function) associated with the denominator polynomial [1 a1 a2 ... ap]. 3. Find the associated reconstruction function using only the filtering function. In order to do this, express xk (0) as the sum of an input filtering and an output filtering.

2.1.2 Distributing the calculation load in an FIR filter

We wish to distribute the calculation load for an FIR filtering algorithm among several processors. Only two methods will be presented. The first one consists of distributing the number of multiplication/accumulation operations (MAC operations) among M branches without changing the processing rate. The second one consists of organizing the calculation in different units, so as to reduce this speed, at the cost of a certain delay.

Paralleled calculations

Consider the filtering equation . For a given M, we define k = mM + r where r ? {0, ..., M - 1}. We obtain:

This expression shows y(n) as the sum of M terms:

The first term is the filtering of a sequence ... , x(n-M), x(n), x(n+M), ... by the filter with the impulse response h(0), h(M), ... The next terms correspond to translated sequences filtered by the filters hr (m) = {h(r), ..., h(r + mM), ... }. The filter hr (m) is called the r-th M-polyphase component of h(n).

Figure 2.3 illustrates this processing method.

Figure 2.3 - A representation of the paralleled process

Exercise 2.2 (Parallel implementation of the FIR filtering) (see p. 192)

Write a program designed to simulate the process described by Figure 2.3. Choose M = 4 and a low-pass, (-0.3, +0.3) band FIR filter with 25 coefficients. The result will be compared to the one obtained through direct filtering.

This method for paralleling does not reduce the processing speed in intermediate filters. Only the number of multiplications per filter is reduced.

2.1.3 FIR block filtering

Let us again consider the FIR filtering equation y(k) = h(0)x(k) + h(1)x(k - 1) + . + h(P)x(k - P). Let:

with:

By organizing the inputs modulo M , we obtain:

By restricting ourselves to the case M = N = 2, the previous expression can be written:

If we assume x0(n) = [x(2n), x(2n-2), ... ]T and x1(n) = [x(2n-1), x(2n-3), ... ]T , we can also write:

where h0 = [h(0), h(2), ... ]T and h1 = [h(1), h(3), ... ]T. If we develop y(2n) and y(2n - 1), we obtain:

Therefore, the calculation of the two terms y(2n) and y(2n - 1) requires the calculation of a total of four terms. However, one of them, x1(n)(h0 + h1), appears twice. If P refers to the length of the filter h, the lengths of h0 and h1 are at the most equal to P/2. Hence the three terms of the calculation of y(2n) and y(2n - 1) correspond to...