## 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,

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 [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:**

Markov Parameters

**Previous Section:**

State Space Realization