Identify the system from input and output

Question

I am using Excel. If I have an input and an output, and I know that a recursive band-pass filter produced the output, how can I can find the coefficients of that filter?

Answer 1

You can use an equation error method and get an estimation of the coefficients by solving an overdetermined system of linear equations. The transfer function of an IIR filter is given by

$$H(z)=\frac{B(z)}{A(z)}=\frac{\sum_{m=0}^{M}b[m]z^{-m}}{1+\sum_{n=1}^{N}a[n]z^{-n}}\tag{1}$$

From (1), for given input and output sequences $x[n]$ and $y[n]$, respectively, the following equation must be satisfied:

$$(a*y)[n]=(x*b)[n]\tag{2}$$

where $*$ denotes convolution. Equation (2) can be written in matrix/vector form:

$$\mathbf{y}=\mathbf{Xb}-\mathbf{Ya}\tag{3}$$

where $\mathbf{y}$ is the vector of the output signal, $\mathbf{b}$ and $\mathbf{a}$ are the unknown coefficient vectors (excluding $a[0]=1$), and $\mathbf{X}$ and $\mathbf{Y}$ are the convolution matrices consisting of the values of the input and output sequences, respectively. Equation (3) can be rewritten as

$$\mathbf{y}=\mathbf{[-Y,X]\cdot[a,b]^T}\tag{4}$$

where $^T$ denotes transposition. The overdetermined system (4) can be solved in a least squares sense using standard methods, such as the one implemented in Matlab or Octave.

Note that in any case you need to estimate/choose the orders $M$ and $N$ of the numerator and denominator polynomials in (1). For most standard filters you have $M=N$.

This little Matlab/Octave script shows how it works:

x=randn(100,1);       % some input signal
a=[1 .3 -.4 .2 .1];   % denominator coefficients
b=[2 3 4 3 2];        % numerator coefficients
y=filter(b,a,x);      % filter output
M=3;                  % guessed numerator order
N=3;                  % guessed denominator order

% build matrices and solve overdetermined linear system
A = flipud(hankel(flipud(y(:)))); A1 = A(:,1); A = A(:,2:N+1);
B = flipud(hankel(flipud(x(:)))); B = B(:,1:M+1);
C = [-A,B];
c = C\A1;
a1 = [1;c(1:N)];
b1 = c(N+1:M+N+1);

% compare outputs of original and estimated filter
y1=filter(b1,a1,x);
n=1:100;
plot(n,y,n,y1)