# Matrix Filter Representations

This appendix introduces various*matrix representations*for digital filters, including the important

*state space formulation*. Additionally, elementary

*system identification*based on a matrix description is described.

##

Introduction

It is illuminating to look at *matrix representations*of digital filters.

^{F.1}Every

*linear*digital filter can be expressed as a

*constant matrix*multiplying the input signal (the

*input vector*) to produce the output signal (vector) ,

*i.e.*,

*finite-length*inputs (to avoid infinite matrices), and the output signal will also be length . Thus, the filter matrix is a square matrix, and the input/output signal vectors are column vectors. More generally, any finite-order

*linear operator*can be expressed as a matrix multiply. For example, the Discrete Fourier Transform (DFT) can be represented by the ``DFT matrix'' , where the column index and row index range from 0 to [84, p. 111].

^{F.2}Even infinite-order linear operators are often thought of as matrices having infinite extent. In summary, if a digital filter is

*linear*, it can be represented by a

*matrix*.

## General Causal Linear Filter Matrix

To be*causal*, the filter output at time cannot depend on the input at any times greater than . This implies that a causal filter matrix must be

*lower triangular*. That is, it must have zeros above the main diagonal. Thus, a causal linear filter matrix will have entries that satisfy for . For example, the general causal, linear, digital-filter matrix operating on three-sample sequences is

or, more explicitly,

While Eq.(F.2) covers the general case of linear, causal, digital filters operating on the space of three-sample sequences, it includes

*time varying*filters, in general. For example, the gain of the ``current input sample'' changes over time as .

## General LTI Filter Matrix

The general linear,*time-invariant*(LTI) matrix is

*Toeplitz*. A

*Toeplitz matrix*is

*constant along all its diagonals*. For example, the general LTI matrix is given by

*banded Toeplitz filter matrix*:

(F.3) |

We could add more rows to obtain more output samples, but the additional outputs would all be zero. In general, if a causal FIR filter is length , then its order is , so to avoid ``cutting off'' the output signal prematurely, we must append at least zeros to the input signal. Appending zeros in this way is often called

*zero padding*, and it is used extensively in spectrum analysis [84]. As a specific example, an order 5 causal FIR filter (length 6) requires 5 samples of zero-padding on the input signal to avoid output truncation. If the FIR filter is

*noncausal*, then zero-padding is needed

*before*the input signal in order not to ``cut off'' the ``pre-ring'' of the filter (the response before time ). To handle arbitrary-length input signals, keeping the filter length at 3 (an order 2 FIR filter), we may simply use a longer banded Toeplitz filter matrix:

*linear operators*[56]. Thus, we may say that every LTI filter corresponds to a

*Toeplitz linear operator*.

## Cyclic Convolution Matrix

An infinite Toeplitz matrix implements, in principle,*acyclic convolution*(which is what we normally mean when we just say ``convolution''). In practice, the convolution of a signal and an impulse response , in which both and are more than a hundred or so samples long, is typically implemented fastest using

*FFT convolution*(

*i.e.*, performing fast convolution using the Fast Fourier Transform (FFT) [84]

^{F.3}). However, the FFT computes

*cyclic convolution*unless sufficient zero padding is used [84]. The matrix representation of cyclic (or ``circular'') convolution is a

*circulant matrix*,

*e.g.*,

^{F.4}For example, the eigenvectors of an circulant matrix are the DFT sinusoids for a length DFT [84]. Similarly, the eigenvalues may be found by simply taking the DFT of the first row. The DFT eigenstructure of circulant matrices is directly related to the DFT convolution theorem [84]. The above circulant matrix , when multiplied times a length 6 vector , implements cyclic convolution of with . Using the DFT to perform the circular convolution can be expressed as

IDFTDFTDFT

where `
' denotes circular convolution. Let
denote
the matrix of sampled DFT sinusoids for a length DFT:
. Then
is the
*DFT matrix*, where `' denotes Hermitian transposition (transposition and complex-conjugation). The DFT of the length- vector can be written as , and the corresponding inverse DFT is . The DFT-eigenstructure of circulant matrices provides that a real circulant matrix having top row satisfies diag, where is the length DFT of , and diag denotes a diagonal matrix with the elements of along the diagonal. Therefore, diag. By the DFT convolution theorem,

## Inverse Filters

Note that the filter matrix is often*invertible*[58]. In that case, we can effectively run the filter

*backwards*:

*not*necessarily correspond to a

*stable*inverse-filter when the lengths of the input and output vectors are allowed to grow larger. For example, the inverted filter matrix may contain truncated

*growing*exponentials, as illustrated in the following

`matlab`example:

> h = toeplitz([1,2,0,0,0],[1,0,0,0,0]) h = 1 0 0 0 0 2 1 0 0 0 0 2 1 0 0 0 0 2 1 0 0 0 0 2 1 > inv(h) ans = 1 0 0 0 0 -2 1 0 0 0 4 -2 1 0 0 -8 4 -2 1 0 16 -8 4 -2 1The inverse of the FIR filter is in fact unstable, having impulse response , , which grows to with . Another point to notice is that the inverse of a banded Toeplitz matrix is not banded (although the inverse of lower-triangular [causal] matrix remains lower triangular). This corresponds to the fact that the inverse of an FIR filter is an IIR filter.

## State Space Realization

Above, we used a matrix multiply to represent convolution of the filter input signal with the filter's impulse response. This only works for FIR filters since an IIR filter would require an infinite impulse-response matrix. IIR filters have an extensively used matrix representation called*state space form*(or ``state space realizations''). They are especially convenient for representing filters with

*multiple inputs*and

*multiple outputs*(MIMO filters). An order digital filter with inputs and outputs can be written in state-space form as follows:

where is the length

*state vector*at discrete time , is a vector of inputs, and the output vector. is the

*state transition matrix*, and it determines the

*dynamics*of the system (its

*poles*, or resonant

*modes*). State-space models are described further in Appendix G. Here, we will only give an illustrative example and make a few observations:

### State Space Filter Realization Example

The digital filter having difference equation^{F.5}

Thus, is the vector of state variables at time , is the state-input gain vector, is the vector of state-gains for the output, and the direct-path gain is . This example is repeated using matlab in §G.7.8 (after we have covered

*transfer functions*). A general procedure for converting any difference equation to state-space form is described in §G.7. The particular state-space model shown in Eq.(F.5) happens to be called

*controller canonical form*, for reasons discussed in Appendix G. The set of all state-space realizations of this filter is given by exploring the set of all

*similarity transformations*applied to any particular realization, such as the control-canonical form in Eq.(F.5). Similarity transformations are discussed in §G.8, and in books on linear algebra [58]. Note that the state-space model replaces an

*th-order difference equation*by a

*vector first-order difference equation*. This provides elegant simplifications in the theory and analysis of digital filters. For example, consider the case , and , so that Eq.(F.4) reduces to

where is the transition matrix, and both and are signal vectors. (This filter has inputs and outputs.) This vector first-order difference equation is analogous to the following scalar first-order difference equation:

*zero-input response*of the filter,

*i.e.*, .) Similarly, setting to in Eq.(F.6) yields

## Time Domain Filter Estimation

*System identification*is the subject of identifying filter coefficients given measurements of the input and output signals [46,78]. For example, one application is

*amplifier modeling*, in which we measure (1) the normal output of an electric guitar (provided by the pick-ups), and (2) the output of a microphone placed in front of the amplifier we wish to model. The guitar may be played in a variety of ways to create a collection of input/output data to use in identifying a model of the amplifier's ``sound.'' There are many commercial products which offer ``virtual amplifier'' presets developed partly in such a way.

^{F.6}One can similarly model electric guitars themselves by measuring the pick signal delivered to the string (as the input) and the normal pick-up-mix output signal. A separate identification is needed for each switch and tone-control position. After identifying a sampling of models, ways can be found to interpolate among the sampled settings, thereby providing ``virtual'' tone-control knobs that respond like the real ones [101]. In the notation of the §F.1, assume we know and and wish to solve for the filter impulse response . We now outline a simple yet practical method for doing this, which follows readily from the discussion of the previous section. Recall that convolution is

*commutative*. In terms of the matrix representation of §F.3, this implies that the input signal and the filter can switch places to give

Here we have indicated the general case for a length causal FIR filter, with input and output signals that go on forever. While is not invertible because it is not square, we can solve for under general conditions by taking the

*pseudoinverse*of . Doing this provides a

*least-squares*

*system identification*method [46]. The

*Moore-Penrose pseudoinverse*is easy to derive.

^{F.7}First multiply Eq.(F.7) on the left by the transpose of in order to obtain a ``square'' system of equations:

Thus, is the

*Moore-Penrose pseudoinverse*of . If the input signal is an

*impulse*(a 1 at time zero and 0 at all other times), then is simply the identity matrix, which is its own inverse, and we obtain . We expect this by definition of the impulse response. More generally, is invertible whenever the input signal is ``sufficiently exciting'' at all frequencies. An LTI filter frequency response can be identified only at frequencies that are excited by the input, and the accuracy of the estimate at any given frequency can be improved by increasing the input signal power at that frequency. [46].

### Effect of Measurement Noise

In practice, measurements are never perfect. Let denote the measured output signal, where is a vector of ``measurement noise'' samples. Then we have*orthogonality principle*[38], the least-squares estimate of is obtained by orthogonally projecting onto the space spanned by the columns of . Geometrically speaking, choosing to minimize the Euclidean distance between and is the same thing as choosing it to minimize the sum of squared estimated measurement errors . The distance from to is minimized when the

*projection error*is orthogonal to every column of , which is true if and only if [84]. Thus, we have, applying the orthogonality principle,

*weighting function*in the least-squares estimate for by replacing in the derivations above by , where is any positive definite matrix (often taken to be diagonal and positive). In the present time-domain formulation, it is difficult to choose a weighting function that corresponds well to

*audio perception*. Therefore, in audio applications, frequency-domain formulations are generally more powerful for linear-time-invariant system identification. A practical example is the frequency-domain equation-error method described in §I.4.4 [78].

### Matlab System Identification Example

The Octave output for the following small matlab example is listed in Fig.F.1:delete('sid.log'); diary('sid.log'); % Log session echo('on'); % Show commands as well as responses N = 4; % Input signal length %x = rand(N,1) % Random input signal - snapshot: x = [0.056961, 0.081938, 0.063272, 0.672761]' h = [1 2 3]'; % FIR filter y = filter(h,1,x) % Filter output xb = toeplitz(x,[x(1),zeros(1,N-1)]) % Input matrix hhat = inv(xb' * xb) * xb' * y % Least squares estimate % hhat = pinv(xb) * y % Numerically robust pseudoinverse hhat2 = xb\y % Numerically superior (and faster) estimate diary('off'); % Close log fileOne fine point is the use of the syntax `` '', which has been a matlab language feature from the very beginning [82]. It is usually more accurate (and faster) than multiplying by the explicit pseudoinverse. It uses the QR decomposition to convert the system of linear equations into upper-triangular form (typically using Householder reflections), determine the effective rank of , and backsolve the reduced triangular system (starting at the bottom, which goes very fast) [29, §6.2].

^{F.8}

+ echo('on'); % Show commands as well as responses + N = 4; % Input signal length + %x = rand(N,1) % Random input signal - snapshot: + x = [0.056961, 0.081938, 0.063272, 0.672761]' x = 0.056961 0.081938 0.063272 0.672761 + h = [1 2 3]'; % FIR filter + y = filter(h,1,x) % Filter output y = 0.056961 0.195860 0.398031 1.045119 + xb = toeplitz(x,[x(1),zeros(1,N-1)]) % Input matrix xb = 0.05696 0.00000 0.00000 0.00000 0.08194 0.05696 0.00000 0.00000 0.06327 0.08194 0.05696 0.00000 0.67276 0.06327 0.08194 0.05696 + hhat = inv(xb' * xb) * xb' * y % Least squares estimate hhat = 1.0000 2.0000 3.0000 3.7060e-13 + % hhat = pinv(xb) * y % Numerically robust pseudoinverse + hhat2 = xb\y % Numerically superior (and faster) estimate hhat2 = 1.0000 2.0000 3.0000 3.6492e-16 |

**Next Section:**

State Space Filters

**Previous Section:**

Analog Filters