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:
Example State Space Filter Transfer Function
Previous Section:
Effect of Measurement Noise