Code covered by the BSD License  

Highlights from
H-infinity sub-optimal causal FIR inverse filter via LMI

from H-infinity sub-optimal causal FIR inverse filter via LMI by Masaaki Nagahara
This computes the optimal inverse filter minimizing the H-infinity norm of the error system via LMI.

dfdesign_w_lmi(phi, w, d, n); 
function [psi gopt] = dfdesign_w_lmi(phi, w, d, n);

% [psi gopt] = dfdesign(phi, d, n);
% DFDESIGN computes the H-infinity optimal inverse FIR filter of phi(z).
% The resulting filter minimizes the H-infinity norm of the error system
%     E_w(z) = [z^(-d) - psi(z)phi(z)]w(z).
% The LMI method based on KYP lemma is used in this function.
%
% [INPUT]
% phi: target system (discrete-time)
% w: weighting function
% d: delay
% n: length of the FIR filter
%
% [OUTPUT]
% psi: resulting inverse FIR filter (in state-space representation)
% gopt: optimal value


Ts = phi.Ts;
delayden = zeros(1, d+1);
delayden(1) = 1;
delay = tf(1, delayden, Ts);

[Aphi Bphi Cphi Dphi] = ssdata(ss(phi, 'min')*w);
[Ad Bd Cd Dd] = ssdata(ss(delay, 'min')*w);

Apsi11 = zeros(n-1,1);
Apsi21 = 0;
Apsi12 = eye(n-1);
Apsi22 = zeros(1,n-1);
Apsi = [Apsi11 Apsi12; Apsi21 Apsi22];
Bpsi = zeros(n, 1);
Bpsi(n) = 1;
Kpsi1 = eye(n);
Kpsi2 = zeros(1,n);
Kpsi = [Kpsi1; Kpsi2];
Mpsi = zeros(n+1,1);
Mpsi(size(Mpsi,1)) = 1;

A13 = zeros(size(Apsi, 1), size(Ad, 2));
A21 = zeros(size(Aphi, 1), size(Apsi, 2));
A23 = zeros(size(Aphi, 1), size(Ad, 2));
A31 = zeros(size(Ad, 1), size(Apsi, 2));
A32 = zeros(size(Ad, 1), size(Aphi, 2));
A = [Apsi Bpsi*Cphi A13; A21 Aphi A23; A31 A32 Ad];
B = [Bpsi*Dphi; Bphi; -Bd];
Cp13 = zeros(size(Kpsi, 1), size(Cd, 2));
Cp = [Kpsi Mpsi*Cphi Cp13];
Cq1 = zeros(1, size(Kpsi, 2));
Cq2 = zeros(1, size(Mpsi*Cphi, 2));
Cq = [Cq1 Cq2 Cd];
Dp = Mpsi*Dphi;
I22 = eye(size(B, 2), size(B, 2));
I33 = eye(1, 1);
Pdim = size(A, 2);

  setlmis([]);
  a=lmivar(2,[1 n+1]);
  P=lmivar(1,[Pdim 1]);
  x=lmivar(2,[1 1]);
  
  lmiterm([1 1 1 P],A',A);                        % LMI #1: A'*P*A
  lmiterm([1 1 1 P],1,-1);                        % LMI #1: -P
  lmiterm([1 2 1 P],B',A);                        % LMI #1: B'*P*A
  lmiterm([1 2 2 x],.5*1,-I22,'s');               % LMI #1: -x*I22 (NON SYMMETRIC?)
  lmiterm([1 2 2 P],B',B);                        % LMI #1: B'*P*B
  lmiterm([1 3 1 a],1,Cp);                        % LMI #1: a*Cp
  lmiterm([1 3 1 0],Cq);                          % LMI #1: Cq
  lmiterm([1 3 2 a],1,Dp);                        % LMI #1: a*Dp
  lmiterm([1 3 3 x],.5*1,-I33,'s');               % LMI #1: -x*I33 (NON SYMMETRIC?)
   
  lmiterm([-2 1 1 0],P);                          % LMI #2: P
   
  lmisys=getlmis;

% mincx
lmin = decnbr(lmisys);
lmic = zeros(lmin, 1);
for j=1:lmin,
    [xj] = defcx(lmisys,j,x);
    lmic(j) = xj;
end  
[gopt, xopt] = mincx(lmisys, lmic);

aopt = dec2mat(lmisys, xopt, 1);
Cpsi = aopt * Kpsi;
Dpsi = aopt * Mpsi;
psi = ss(Apsi, Bpsi, Cpsi, Dpsi, Ts);

Contact us