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?
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.
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];
For your own implementation see: