Matlab System Identification Example 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.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  '', 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) [30, §6.2].6.14

 + 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 

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

[How to cite this work] [Order a printed hardcopy]