System identification

Estimate an LTI system from input-output data using least squares.

Ivan Selesnick
selesi@poly.edu

Contents

Start

clear
close all

Load data

load data.txt;          % First column: input.   Second column: output

whos

x = data(:, 1);         % input signal
y = data(:, 2);         % output signal

N = length(y);
n = 0:N-1;
  Name        Size            Bytes  Class     Attributes

  data      100x2              1600  double              

Display data

figure(1)
clf
subplot(2, 1, 1)
plot(n, x)
title('Input signal')
YL1 = [-2 2];
ylim(YL1)

subplot(2, 1, 2)
plot(n, y)
title('Output signal')
ylim(YL1)

Least squares solution (length 5)

M = 5;                                      % M : length of impulse response

X = toeplitz(x, [x(1) zeros(1, M-1)]);      % X : convolution matrix

h = X \ y                                   % h : impulse response estimate
h =

    0.1086
    0.2395
    0.5281
    0.5150
    0.3517

figure(1)

subplot(2, 1, 1)
plot(0:M-1, h)
title(sprintf('Estimated impulse response (length %d)', M))

r = y - X * h;                               % r : residual
RMSE = sqrt( sum(r.^2) );

subplot(2, 1, 2)
plot(n, r)
title(sprintf('Residual (RMSE = %.2f)', RMSE))
ylim(YL1)

Least squares solution (length 10)

Increasing the impulse response leads to a smaller residual

M = 10;                                     % M : length of impulse response

X = toeplitz(x, [x(1) zeros(1, M-1)]);      % X : convolution matrix

h = X \ y                                   % h : impulse response estimate
h =

    0.1021
    0.2076
    0.4910
    0.5052
    0.3843
    0.3010
    0.0996
    0.1115
    0.0024
    0.0083

figure(1)
clf

subplot(2, 1, 1)
plot(0:M-1, h)
title(sprintf('Estimated impulse response (length %d)', M))

r = y - X * h;                               % r : residual
RMSE = sqrt( sum(r.^2) );

subplot(2, 1, 2)
plot(n, r)
title(sprintf('Residual (RMSE = %.2f)', RMSE))
ylim(YL1)

RMSE versus impulse response length

M = 20;
for m = 1:M                                     % m : impulse response length
    X = toeplitz(x, [x(1) zeros(1, m-1)]);
    h = X \ y;
    err(m) = sqrt(sum( (X*h - y).^2 ));
end

figure(1)
clf
plot(1:M, err, '.-')
xlabel('Length of impulse response')
ylabel('RMSE')
title('RMSE vs impulse response length');