Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

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 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) [29, §6.2].F.8

Figure F.1: 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


Next  |  Prev  |  Up  |  Top  |  Index  |  JOS Index  |  JOS Pubs  |  JOS Home  |  Search

[How to cite this work]  [Order a printed hardcopy]  [Comment on this page via email]

``Introduction to Digital Filters with Audio Applications'', by Julius O. Smith III, (September 2007 Edition).
Copyright © 2015-04-22 by Julius O. Smith III
Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
CCRMA