Matlab System Identification Example

Chapter Contents:

Search Introduction to Digital Filters

Would you like to be notified by email when Julius Orion Smith III publishes a new entry into his blog?

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

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

Order a Hardcopy of Introduction to Digital Filters

Previous: Effect of Measurement Noise
Next: State Space Filters

written by Julius Orion Smith III
Julius Smith's background is in electrical engineering (BS Rice 1975, PhD Stanford 1983). He is presently Professor of Music and Associate Professor (by courtesy) of Electrical Engineering at Stanford's Center for Computer Research in Music and Acoustics (CCRMA), teaching courses and pursuing research related to signal processing applied to music and audio systems. See http://ccrma.stanford.edu/~jos/ for details.

Sign in

Search Online Books

Free Online Books

Ads

Chapters

Introduction to Digital Filters
    Matrix Filter Representations
       Time Domain Filter Estimation
          Matlab System Identification Example