Matlab System Identification Example

The Octave output for the following small matlab example is listed in Fig.5.9:

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 file

One fine point is the use of the syntax `` $\underline{h}= \mathbf{x}\backslash\underline{y}$'', 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 $\mathbf{x}$, and backsolve the reduced triangular system (starting at the bottom, which goes very fast) [30, §6.2].^6.14

Figure 5.9: Time-domain system-identification matlab example.

 

+ 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

+ diary('off'); % Close log file