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.
It is illuminating to look at matrix representations of digital filters.F.1Every 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.,
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.2Even 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.
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
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 . 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:
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) F.3). However, the FFT computes cyclic convolution unless sufficient zero padding is used . The matrix representation of cyclic (or ``circular'') convolution is a circulant matrix, e.g.,
The DFT eigenstructure of circulant matrices is directly related to the DFT convolution theorem . 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
Premultiplying by the IDFT matrix yields
> 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 Filter Realization Example
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 .
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 .
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:
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 .
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.
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 .
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. .
It is also straightforward to introduce a 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 .
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 . 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
State Space Filters