4
$\begingroup$

I designed an IIR filter and want to do zero, step and projection initialization. There is mention in different reference papers but no explicit examples for Matlab. How do you do IIR initialization, for the previously mentioned types, in Matlab?

$\endgroup$
  • 1
    $\begingroup$ Is this question really about How to do it in MATLAB? or do you really want to understand about IIR initialization conceptually? $\endgroup$ – Dipan Mehta Apr 17 '12 at 8:29
4
$\begingroup$

Not easy. The Matlab filter() function implements a general purpose IIR filter and allows you to pass in an initial state. However, they don't publish the exact definition of the state and the algorithm deployed. It seems to be close to a transposed form II filter but with some subtle numerical differences, so the state doesn't quite match.

In any case it would be better to break this down into second order sections using sosfilter(). However that doesn't allow you to set an initial state. You may have to write your own filter function for that.

YOu still can use filter() to some extent. The default is zero initialization. You can also do step initialization by simply running a constant input for a while and collecting the state as an output argument from the filter, like this:

[b a] = butter(4,1000*2/44100);
% create state for a constant input of "1"
[~,state1] = filter(b,a,ones(1000,1));
% now filter for a constant of 2 but with a state for "1" input
ystep = filter(b,a,2*ones(100,1),state1);
plot(ystep);

Not sure about projection initialization, though.

$\endgroup$
1
$\begingroup$

Here is a Matlab implementation of projection initialization for a 2nd order IIR. For details, see Initialization for Improved IIR Filter Performance by Chornoboy.

% x is your input signal.
% y is the output signal.
% fs is the sample rate.

% This example uses a 2nd order Butterworth high-pass filter with fc = 0.3 Hz.
fc = 0.3;
[bh, ah] = butter(2, fc/fs, 'high');

A = [bh(2)-bh(1)*ah(2); bh(3)-bh(1)*ah(3)];
B = [-ah(2), -ah(3); 1 0];
C = [1;0];

Np = 1.0 * fs; % Set Np to the length of data you want to train on.
F = zeros(Np, 2);
G = diag(bh(1) * ones(1, Np));
for ix = 1:Np
   F(ix, :) = A'*B^(ix-1);
   if ix > 1
      for cix = (ix-1):-1:1
         G(ix, cix) = A'*B^(ix-cix-1)*C;
      end
   end
end

% Filter the data using state variables.
Mm1 = -inv(F'*F)*F'*G * x(1:Np); % M(-1)
y = zeros(size(x));
for ix = 1:length(y)
   Mn = B*Mm1 + C*x(ix);
   y(ix) = A'*Mm1 + bh(1)*x(ix);
   Mm1 = Mn;
end

Note that step initialization is done by setting Mm1 = 1/(1 + ah(2) + ah(3)) * [1; 1];

And zero initialization is done by setting Mm1 = [0;0];

$\endgroup$
-4
$\begingroup$

For your own implementation see:

http://www.emt.tugraz.at/publications/diplomarbeiten/da_hoebenreich/node21.html#fig:StateSpaceDiagram

$\endgroup$
  • $\begingroup$ Welcome to DSP.SE. This site, and the whole StackExchange universe frowns on link-only answers. Can you please summarize what that page says, and edit this answer to add it? $\endgroup$ – Peter K. Jan 21 '16 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.