Is there a widely available implementation of an adaptive, recursive, numerically stable IIR filter?

Question

I need to identify the coefficients of a linear, causal, time-invariant physical system that can be described by a classical state-space formulation.

For the sake of the example, suppose that the system has two states, only one of which is observed, and one input. Suppose that this state-space is described by the following matrices:

  a = 
               x1          x2
   x1   -1.25e-05    1.25e-05
   x2   1.389e-05  -2.083e-05

  b = 
              u1
   x1          0
   x2  6.944e-06

  c = 
       x1  x2
   y1   1   0

  d = 
       u1
   y1   0

The discrete, equivalent, zeroth-order-hold transfer function is given by:

  3.481e-05 z + 3.446e-05
  -----------------------
   z^2 - 1.97 z + 0.9704

In filter terminology, this is clearly an IIR filter because of the higher-order coefficients in the denominator.

I'm looking for a suitable adaptive filter algorithm that will:

• be numerically stable
• allow for online updates, i.e. be recursive

From my research so far, I've read that if you treat the recorded signal values as inputs to this filter, then it's possible to treat this filter as a FIR filter. This is called the equation error approach. One can then apply, for example, the QR-RLS adaptive FIR filter to find the coefficients.

But in my trials on the real system I often find filter coefficients that are not physically possible, and I just don't know if my hand-crafted adaptive filter is correct or not. I would like to try with an "official" implementation of such an adaptive filter.

Therefore, my question is whether there exists (for MATLAB or in C) a freely available implementation of an adaptive IIR filter that has the same advantages of a QR-RLS FIR adaptive filter.

Answer 1

While I may not be the final authority, I have looked into this. In short, no.

From my own research, the problem is stability. When the filter is extended to have a wider bandwidth, which of course means that the y[n] values and lags 'move quickly', and then is quenched to a lower bandwidth, where the y[n] values and lags 'move slowly', the y[n] buffer may be loaded with fast oscillated entries that appear as initial conditions to the low bandwidth filter. The filter can easily go unstable. If the bandwidth is altered in an adiabatic manner this instability may be avoided (I'm supposing) but then there's an additional constraint on your adaptivity.

Now, you don't have to take my word for it. In "Adaptive Filters: Theory and Applications" (2nd) page 5 (Amazon), Boroujeny writes, "However, as we shall see in the later chapters, because of the many difficulties involved in adaptation of IIR filters, their application in the area of adaptive filters is rather limited." The author states that the adaptation process can place (digital) poles outside of the unit circle, even if they started within the unit circle. The rest of his book is mostly dedicated to FIR filters.

I hope this information helps.

Answer 2

The equation error approach should have worked (assuming your model was the same order as the unknown system). While it's not as efficient as QRD-RLS algorithms you can always use a Kalman filter to do recursive least squares.

There's an open source Kalman filtering toolbox available here. Below is an example of how to solve the parameter estimation problem for the example system you provided.

The basic idea is to model the output of an unknown system as a weighted combination of delayed outputs and inputs (also known as an ARMA model). The weights are unknown so they become the states that the Kalman filter tries to estimate, while the inputs and outputs form a time-varying observation matrix.

clear all
close all

addpath ekfukf

% unknown system
z = tf('z');
unknown_tf = (3.481e-5*z + 3.446e-5)/(z^2 - 1.97*z + 0.9704);
unknown_ss = ss(unknown_tf);
A = unknown_ss.a;
B = unknown_ss.b;
C = unknown_ss.c;
D = unknown_ss.d;
x = zeros(2, 1);

M = 1000; % simulation length
u = randn(M, 1); % vector of observed inputs
y = zeros(M, 1); % initialize output vector
e = zeros(M, 1); % output prediction error

% kalman filter parameters
dim = 4; % number of unknown parameters (2 numerator coeffs, 2 denominator)
F = eye(dim); % state transition matrix
R = 1e-6; % measurement noise variance  
Q = 0*eye(dim); % process noise covariance
P = 10*eye(dim); % state-error covariance
alpha = zeros(2, 1); % initial estimate of denominator coefficiens 
beta = zeros(2, 1); % initial estimate of numerator coefficients
state = [alpha; beta]; 
Hu = zeros(2, 1); % observation vector containing observed inputs
Hy = zeros(2, 1); % observation vector containing observed outputs

for k=2:M
    % compute kalman filter time-update
    [state, P] = kf_predict(state, P, F, Q);

    % update observation vectors
    Hy = circshift(Hy, 1);
    Hy(1) = y(k-1);
    Hu = circshift(Hu, 1);
    Hu(1) = u(k);

    % observe system output
    x = A*x + B*u(k);
    y(k) = C*x + D*u(k);

    % compute kalman filter measurement update
    [state, P, K, yfilt, S, likelihood] = kf_update(state, P, y(k), [Hy; Hu]', R); 
    e(k) = y(k) - yfilt;
end

% exctract numerator and denominator coefficients
alpha = state(1:2);
beta = state(3:4);

% construct numerator and denominator polynomial estimates
est_num = [0, beta'];
est_den = [1, -alpha'];

% plot estimated and true frequency response
W = pi*logspace(-3, 0, 300);
Ghat = freqz(est_num, est_den, W);
G = freqz(unknown_tf.num{1}, unknown_tf.den{1}, W);
semilogx(W./pi, 20*log10(abs(G)), 'k', 'linewidth', 2)
hold on
semilogx(W./pi, 20*log10(abs(Ghat)), 'g')
title('estimated (green) and true (black) magnitude response')
ylabel('magnitude (dB)');
xlabel('normalized frequency');

% plot prediction error 
figure
plot(e)
title('prediction error vs. time');
ylabel('error');
xlabel('sample index');

It works!

The prediction error (trying to predict the output of the unknown system) gives us an idea of how quickly the model converged. Convergence speed is highly dependent on the measurement noise variance.