Designing Complementary filters for sensor fusion using \(\mathcal{H}_\infty\) synthesis - Matlab Computation

1. H-Infinity synthesis of complementary filters
2. Design of complementary filters used in the Active Vibration Isolation System at the LIGO
3. “Closed-Loop” complementary filters
4. Synthesis of three complementary filters
5. Matlab Scripts

This report is also available as a pdf.

The present document is a companion file for the journal paper (Dehaeze, Vermat, and Collette 2021). All the Matlab (MATLAB 2020) scripts used in the paper are here shared and explained.

This document is divided into the following sections also corresponding to the paper sections:

Section 1: the shaping of complementary filters is written as an \(\mathcal{H}_\infty\) optimization problem using weighting functions. The weighting function design is discussed and the method is applied for the design of a set of simple complementary filters.
Section 2: the effectiveness of the proposed complementary filter synthesis strategy is demonstrated by designing complex complementary filters used in the first isolation stage at the LIGO
Section 3: complementary filters are designed using the classical feedback loop
Section 4: the proposed design method is generalized for the design of a set of three complementary filters
Section 5: complete Matlab scripts and functions developed are listed

1. H-Infinity synthesis of complementary filters

1.1. Synthesis Architecture

In order to generate two complementary filters with a wanted shape, the generalized plant of Figure 1 can be used. The included weights \(W_1(s)\) and \(W_2(s)\) are used to specify the upper bounds of the complementary filters being generated.

Figure 1: Generalized plant used for the \(\mathcal{H}_\infty\) synthesis of a set of two complementary fiters

Applied the standard \(\mathcal{H}_\infty\) synthesis on this generalized plant will give a transfer function \(H_2(s)\) (see Figure 2) such that the \(\mathcal{H}_\infty\) norm of the transfer function from \(w\) to \([z_1,\ z_2]\) is less than one \eqref{eq:h_inf_objective}.

Figure 2: Generalized plant with the synthesized filter obtained after the \(\mathcal{H}_\infty\) synthesis

\begin{equation} \left\| \begin{array}{c} (1 - H_2(s)) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \label{eq:h_inf_objective} \end{equation}

Thus, if the synthesis is successful and the above condition is verified, we can define \(H_1(s)\) to be the complementary of \(H_2(s)\) \eqref{eq:H1_complementary_of_H2} and we have condition \eqref{eq:shaping_comp_filters} verified.

\begin{equation} H_1(s) = 1 - H_2(s) \label{eq:H1_complementary_of_H2} \end{equation} \begin{equation} \left\| \begin{array}{c} H_1(s) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \quad \Longrightarrow \quad \left\{ \begin{array}{c} |H_1(j\omega)| < \frac{1}{|W_1(j\omega)|}, \quad \forall \omega \\ |H_2(j\omega)| < \frac{1}{|W_2(j\omega)|}, \quad \forall \omega \end{array} \right. \label{eq:shaping_comp_filters} \end{equation}

We then see that \(W_1(s)\) and \(W_2(s)\) can be used to set the wanted upper bounds of the magnitudes of both \(H_1(s)\) and \(H_2(s)\).

The presented synthesis method therefore allows to shape two filters \(H_1(s)\) and \(H_2(s)\) \eqref{eq:shaping_comp_filters} while ensuring their complementary property \eqref{eq:H1_complementary_of_H2}.

The complete Matlab script for this part is given in Section 5.1.

1.2. Design of Weighting Function - Proposed formula

A formula is proposed to help the design of the weighting functions:

\begin{equation} W(s) = \left( \frac{ \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{G_0}{G_c}\right)^{\frac{1}{n}} }{ \left(\frac{1}{G_\infty}\right)^{\frac{1}{n}} \frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{1}{G_c}\right)^{\frac{1}{n}} }\right)^n \label{eq:weighting_function_formula} \end{equation}

The parameters permits to specify:

the low frequency gain: \(G_0 = lim_{\omega \to 0} |W(j\omega)|\)
the high frequency gain: \(G_\infty = lim_{\omega \to \infty} |W(j\omega)|\)
the absolute gain at \(\omega_0\): \(G_c = |W(j\omega_0)|\)
the absolute slope between high and low frequency: \(n\)

The general shape of a weighting function generated using the formula is shown in figure 3.

Figure 3: Magnitude of the weighting function generated using formula \eqref{eq:weighting_function_formula}

1.3. Weighting functions for the design of two complementary filters

The weighting function formula \eqref{eq:weighting_function_formula} is used to generate the upper bounds of two complementary filters that we wish to design.

The matlab function generateWF is described in Section 5.5.

%% Design of the Weighting Functions
W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);

The inverse magnitude of these two weighting functions are shown in Figure 4.

Figure 4: Inverse magnitude of the design weighting functions

1.4. Synthesis of the complementary filters

The generalized plant of Figure 1 is defined as follows:

%% Generalized Plant
P = [W1 -W1;
     0   W2;
     1   0];

And the \(\mathcal{H}_\infty\) synthesis is performed using the hinfsyn command.

%% H-Infinity Synthesis
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');

[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');

  Test bounds:  0.3223 <=  gamma  <=  1000

    gamma        X>=0        Y>=0       rho(XY)<1    p/f
  1.795e+01     1.4e-07     0.0e+00     1.481e-16     p
  2.406e+00     1.4e-07     0.0e+00     3.604e-15     p
  8.806e-01    -3.1e+02 #  -1.4e-16     7.370e-19     f
  1.456e+00     1.4e-07     0.0e+00     1.499e-18     p
  1.132e+00     1.4e-07     0.0e+00     8.587e-15     p
  9.985e-01     1.4e-07     0.0e+00     2.331e-13     p
  9.377e-01    -7.7e+02 #  -6.6e-17     3.744e-14     f
  9.676e-01    -2.0e+03 #  -5.7e-17     1.046e-13     f
  9.829e-01    -6.6e+03 #  -1.1e-16     2.949e-13     f
  9.907e-01     1.4e-07     0.0e+00     2.374e-19     p
  9.868e-01    -1.6e+04 #  -6.4e-17     5.331e-14     f
  9.887e-01    -5.1e+04 #  -1.5e-17     2.703e-19     f
  9.897e-01     1.4e-07     0.0e+00     1.583e-11     p
  Limiting gains...
  9.897e-01     1.5e-07     0.0e+00     1.183e-12     p
  9.897e-01     6.9e-07     0.0e+00     1.365e-12     p

  Best performance (actual): 0.9897

As shown above, the obtained \(\mathcal{H}_\infty\) norm of the transfer function from \(w\) to \([z_1,\ z_2]\) is found to be less than one meaning the synthesis is successful.

We then define the filter \(H_1(s)\) to be the complementary of \(H_2(s)\) \eqref{eq:H1_complementary_of_H2}.

%% Define H1 to be the complementary of H2
H1 = 1 - H2;

The function generateCF can also be used to synthesize the complementary filters. This function is described in Section 5.6.

[H1, H2] = generateCF(W1, W2);

1.5. Obtained Complementary Filters

The obtained complementary filters are shown below and are found to be of order 5. Their bode plots are shown in figure 5 and compare with the defined upper bounds.

zpk(H1)
ans =

            (s+1.289e05) (s+153.6) (s+3.842)^3
  -------------------------------------------------------
  (s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)

zpk(H2)
ans =

         125.61 (s+3358)^2 (s^2 + 46.61s + 813.8)
  -------------------------------------------------------
  (s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)

Figure 5: Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis

2. Design of complementary filters used in the Active Vibration Isolation System at the LIGO

In this section, the proposed method for the design of complementary filters is validated for the design of a set of two complex complementary filters used for the first isolation stage at the LIGO (Hua, Adhikari, et al. 2004).

The complete Matlab script for this part is given in Section 5.2.

2.1. Specifications

The specifications for the set of complementary filters (\(L_1,H_1\)) used at the LIGO are summarized below (for further details, refer to (Hua, Debra, et al. 2004)):

From 0 to 0.008 Hz, the magnitude \(|L_1(j\omega)|\) should be less or equal to \(8 \times 10^{-4}\)
Between 0.008 Hz to 0.04 Hz, the filter \(L_1(s)\) should attenuate the input signal proportional to frequency cubed
Between 0.04 Hz to 0.1 Hz, the magnitude \(|L_1(j\omega)|\) should be less than \(3\)
Above 0.1 Hz, the magnitude \(|H_1(j\omega)|\) should be less than \(0.045\)

The specifications are translated into upper bounds of the complementary filters and are shown in Figure 6.

Figure 6: Specification for the LIGO complementary filters

2.2. FIR Filter

To replicated the complementary filters developed in (Hua, Adhikari, et al. 2004), the CVX Matlab toolbox (Grant and Boyd 2014) is used.

The CVX toolbox is initialized and the SeDuMi solver (Sturm 1999) is used.

%% Initialized CVX
cvx_startup;
cvx_solver sedumi;

We define the frequency vectors on which we will constrain the norm of the FIR filter.

%% Frequency vectors
w1 = 0:4.06e-4:0.008;
w2 = 0.008:4.06e-4:0.04;
w3 = 0.04:8.12e-4:0.1;
w4 = 0.1:8.12e-4:0.83;

The order \(n\) of the FIR filter is defined.

%% Filter order
n = 512;

%% Initialization of filter responses
A1 = [ones(length(w1),1),  cos(kron(w1'.*(2*pi),[1:n-1]))];
A2 = [ones(length(w2),1),  cos(kron(w2'.*(2*pi),[1:n-1]))];
A3 = [ones(length(w3),1),  cos(kron(w3'.*(2*pi),[1:n-1]))];
A4 = [ones(length(w4),1),  cos(kron(w4'.*(2*pi),[1:n-1]))];

B1 = [zeros(length(w1),1), sin(kron(w1'.*(2*pi),[1:n-1]))];
B2 = [zeros(length(w2),1), sin(kron(w2'.*(2*pi),[1:n-1]))];
B3 = [zeros(length(w3),1), sin(kron(w3'.*(2*pi),[1:n-1]))];
B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))];

And the convex optimization is run.

%% Convex optimization
cvx_begin

variable y(n+1,1)

% t
maximize(-y(1))

for i = 1:length(w1)
    norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3;
end

for  i = 1:length(w2)
    norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3*(2*pi*w2(i)/(0.008*2*pi))^3;
end

for i = 1:length(w3)
    norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3;
end

for i = 1:length(w4)
    norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1);
end

cvx_end

h = y(2:end);

cvx_begin
variable y(n+1,1)
% t
maximize(-y(1))
for i = 1:length(w1)
    norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3;
end
for  i = 1:length(w2)
    norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3*(2*pi*w2(i)/(0.008*2*pi))^3;
end
for i = 1:length(w3)
    norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3;
end
for i = 1:length(w4)
    norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1);
end
cvx_end

Calling SeDuMi 1.34: 4291 variables, 1586 equality constraints
   For improved efficiency, SeDuMi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.34 (beta) by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 1586, order n = 3220, dim = 4292, blocks = 1073
nnz(A) = 1100727 + 0, nnz(ADA) = 1364794, nnz(L) = 683190
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.11E+02 0.000
  1 :  -2.58E+00 1.25E+02 0.000 0.3049 0.9000 0.9000   4.87  1  1  3.0E+02
  2 :  -2.36E+00 3.90E+01 0.000 0.3118 0.9000 0.9000   1.83  1  1  6.6E+01
  3 :  -1.69E+00 1.31E+01 0.000 0.3354 0.9000 0.9000   1.76  1  1  1.5E+01
  4 :  -8.60E-01 7.10E+00 0.000 0.5424 0.9000 0.9000   2.48  1  1  4.8E+00
  5 :  -4.91E-01 5.44E+00 0.000 0.7661 0.9000 0.9000   3.12  1  1  2.5E+00
  6 :  -2.96E-01 3.88E+00 0.000 0.7140 0.9000 0.9000   2.62  1  1  1.4E+00
  7 :  -1.98E-01 2.82E+00 0.000 0.7271 0.9000 0.9000   2.14  1  1  8.5E-01
  8 :  -1.39E-01 2.00E+00 0.000 0.7092 0.9000 0.9000   1.78  1  1  5.4E-01
  9 :  -9.99E-02 1.30E+00 0.000 0.6494 0.9000 0.9000   1.51  1  1  3.3E-01
 10 :  -7.57E-02 8.03E-01 0.000 0.6175 0.9000 0.9000   1.31  1  1  2.0E-01
 11 :  -5.99E-02 4.22E-01 0.000 0.5257 0.9000 0.9000   1.17  1  1  1.0E-01
 12 :  -5.28E-02 2.45E-01 0.000 0.5808 0.9000 0.9000   1.08  1  1  5.9E-02
 13 :  -4.82E-02 1.28E-01 0.000 0.5218 0.9000 0.9000   1.05  1  1  3.1E-02
 14 :  -4.56E-02 5.65E-02 0.000 0.4417 0.9045 0.9000   1.02  1  1  1.4E-02
 15 :  -4.43E-02 2.41E-02 0.000 0.4265 0.9004 0.9000   1.01  1  1  6.0E-03
 16 :  -4.37E-02 8.90E-03 0.000 0.3690 0.9070 0.9000   1.00  1  1  2.3E-03
 17 :  -4.35E-02 3.24E-03 0.000 0.3641 0.9164 0.9000   1.00  1  1  9.5E-04
 18 :  -4.34E-02 1.55E-03 0.000 0.4788 0.9086 0.9000   1.00  1  1  4.7E-04
 19 :  -4.34E-02 8.77E-04 0.000 0.5653 0.9169 0.9000   1.00  1  1  2.8E-04
 20 :  -4.34E-02 5.05E-04 0.000 0.5754 0.9034 0.9000   1.00  1  1  1.6E-04
 21 :  -4.34E-02 2.94E-04 0.000 0.5829 0.9136 0.9000   1.00  1  1  9.9E-05
 22 :  -4.34E-02 1.63E-04 0.015 0.5548 0.9000 0.0000   1.00  1  1  6.6E-05
 23 :  -4.33E-02 9.42E-05 0.000 0.5774 0.9053 0.9000   1.00  1  1  3.9E-05
 24 :  -4.33E-02 6.27E-05 0.000 0.6658 0.9148 0.9000   1.00  1  1  2.6E-05
 25 :  -4.33E-02 3.75E-05 0.000 0.5972 0.9187 0.9000   1.00  1  1  1.6E-05
 26 :  -4.33E-02 1.89E-05 0.000 0.5041 0.9117 0.9000   1.00  1  1  8.6E-06
 27 :  -4.33E-02 9.72E-06 0.000 0.5149 0.9050 0.9000   1.00  1  1  4.5E-06
 28 :  -4.33E-02 2.94E-06 0.000 0.3021 0.9194 0.9000   1.00  1  1  1.5E-06
 29 :  -4.33E-02 9.73E-07 0.000 0.3312 0.9189 0.9000   1.00  2  2  5.3E-07
 30 :  -4.33E-02 2.82E-07 0.000 0.2895 0.9063 0.9000   1.00  2  2  1.6E-07
 31 :  -4.33E-02 8.05E-08 0.000 0.2859 0.9049 0.9000   1.00  2  2  4.7E-08
 32 :  -4.33E-02 1.43E-08 0.000 0.1772 0.9059 0.9000   1.00  2  2  8.8E-09

iter seconds digits       c*x               b*y
 32     49.4   6.8 -4.3334083581e-02 -4.3334090214e-02
|Ax-b| =   3.7e-09, [Ay-c]_+ =   1.1E-10, |x|=  1.0e+00, |y|=  2.6e+00

Detailed timing (sec)
   Pre          IPM          Post
3.902E+00    4.576E+01    1.035E-02
Max-norms: ||b||=1, ||c|| = 3,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 4.26267.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): -0.0433341
h = y(2:end);

Finally, the filter response is computed over the frequency vector defined and the result is shown on figure 7 which is very close to the filters obtain in (Hua, Adhikari, et al. 2004).

%% Combine the frequency vectors to form the obtained filter
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;

Figure 7: FIR Complementary filters obtain after convex optimization

2.3. Weighting function design

The weightings function that will be used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters are now designed.

These weights will determine the order of the obtained filters.

Here are the requirements on the filters:

reasonable order
to be as close as possible to the specified upper bounds
stable and minimum phase

The weighting function for the High Pass filter is defined as follows:

%% Design of the weight for the high pass filter
w1 = 2*pi*0.008; x1 = 0.35;
w2 = 2*pi*0.04;  x2 = 0.5;
w3 = 2*pi*0.05;  x3 = 0.5;

% Slope of +3 from w1
wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1);
% Little bump from w2 to w3
wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1);
% No Slope at high frequencies
wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1);
% Little bump between w2 and w3
w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1;
wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n;

wH = 1/wH;
wH = minreal(ss(wH));

And the weighting function for the Low pass filter is taken as a Chebyshev Type I filter.

%% Design of the weight for the low pass filter
n = 20; % Filter order
Rp = 1; % Peak to peak passband ripple
Wp = 2*pi*0.102; % Edge frequency

% Chebyshev Type I filter design
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);

wL = 1/wL;
wL = minreal(ss(wL));

The inverse magnitude of the weighting functions are shown in Figure 8.

Figure 8: Weights for the \(\mathcal{H}_\infty\) synthesis

2.4. Synthesis of the complementary filters

The generalized plant of figure 1 is defined.

%% Generalized plant for the H-infinity Synthesis
P = [0   wL;
     wH -wH;
     1   0];

And the standard \(\mathcal{H}_\infty\) synthesis using the hinfsyn command is performed.

%% Standard H-Infinity synthesis
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');

[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
Resetting value of Gamma min based on D_11, D_12, D_21 terms

Test bounds:      0.3276 <  gamma  <=      1.8063

  gamma    hamx_eig  xinf_eig  hamy_eig   yinf_eig   nrho_xy   p/f
    1.806   1.4e-02 -1.7e-16   3.6e-03   -4.8e-12    0.0000    p
    1.067   1.3e-02 -4.2e-14   3.6e-03   -1.9e-12    0.0000    p
    0.697   1.3e-02 -3.0e-01#  3.6e-03   -3.5e-11    0.0000    f
    0.882   1.3e-02 -9.5e-01#  3.6e-03   -1.2e-34    0.0000    f
    0.975   1.3e-02 -2.7e+00#  3.6e-03   -1.6e-12    0.0000    f
    1.021   1.3e-02 -8.7e+00#  3.6e-03   -4.5e-16    0.0000    f
    1.044   1.3e-02 -6.5e-14   3.6e-03   -3.0e-15    0.0000    p
    1.032   1.3e-02 -1.8e+01#  3.6e-03    0.0e+00    0.0000    f
    1.038   1.3e-02 -3.8e+01#  3.6e-03    0.0e+00    0.0000    f
    1.041   1.3e-02 -8.3e+01#  3.6e-03   -2.9e-33    0.0000    f
    1.042   1.3e-02 -1.9e+02#  3.6e-03   -3.4e-11    0.0000    f
    1.043   1.3e-02 -5.3e+02#  3.6e-03   -7.5e-13    0.0000    f

 Gamma value achieved:     1.0439

The obtained \(\mathcal{H}_\infty\) norm is found to be close than one meaning the synthesis is successful.

The high pass filter \(H_H(s)\) is defined to be the complementary of the synthesized low pass filter \(H_L(s)\):

\begin{equation} H_H(s) = 1 - H_L(s) \end{equation}

%% High pass filter as the complementary of the low pass filter
Hh = 1 - Hl;

The size of the filters is shown to be equal to the sum of the weighting functions orders.

size(Hh), size(Hl)
State-space model with 1 outputs, 1 inputs, and 27 states.
State-space model with 1 outputs, 1 inputs, and 27 states.

The magnitude of the obtained filters as well as the requirements are shown in Figure 9.

Figure 9: Obtained complementary filters using the \(\mathcal{H}_\infty\) synthesis

2.5. Comparison of the FIR filters and synthesized filters

Let’s now compare the FIR filters designed in (Hua, Adhikari, et al. 2004) with the with complementary filters obtained with the \(\mathcal{H}_\infty\) synthesis.

This is done in Figure 10, and both set of filters are found to be very close to each other.

Figure 10: Comparison between the FIR filters developped for LIGO and the \(\mathcal{H}_\infty\) complementary filters

3. “Closed-Loop” complementary filters

In this section, the classical feedback architecture shown in Figure 11 is used for the design of complementary filters.

Figure 11: “Closed-Loop” complementary filters

The complete Matlab script for this part is given in Section 5.3.

3.1. Weighting Function design

Weighting functions using the generateWF Matlab function are designed to specify the upper bounds of the complementary filters to be designed. These weighting functions are the same as the ones used in Section 1.3.

%% Design of the Weighting Functions
W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);

3.2. Generalized plant

The generalized plant of Figure 12 is defined below:

%% Generalized plant for "closed-loop" complementary filter synthesis
P = [ W1 0   1;
     -W1 W2 -1];

Figure 12: Generalized plant used for the \(\mathcal{H}_\infty\) synthesis of “closed-loop” complementary filters

3.3. Synthesis of the closed-loop complementary filters

And the standard \(\mathcal{H}_\infty\) synthesis is performed.

%% Standard H-Infinity Synthesis
[L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');

  Test bounds:  0.3191 <=  gamma  <=  1.669

    gamma        X>=0        Y>=0       rho(XY)<1    p/f
  7.299e-01    -1.5e-19    -2.4e+01 #   1.555e-18     f
  1.104e+00     0.0e+00     1.6e-07     2.037e-19     p
  8.976e-01    -3.2e-16    -1.4e+02 #   5.561e-16     f
  9.954e-01     0.0e+00     1.6e-07     1.041e-15     p
  9.452e-01    -1.1e-15    -3.8e+02 #   4.267e-15     f
  9.700e-01    -6.5e-16    -1.6e+03 #   9.876e-15     f
  9.826e-01     0.0e+00     1.6e-07     8.775e-39     p
  9.763e-01    -5.0e-16    -6.2e+03 #   3.519e-14     f
  9.795e-01     0.0e+00     1.6e-07     6.971e-20     p
  9.779e-01    -1.9e-31    -2.2e+04 #   5.600e-18     f
  9.787e-01     0.0e+00     1.6e-07     5.546e-19     p
  Limiting gains...
  9.789e-01     0.0e+00     1.6e-07     1.084e-13     p
  9.789e-01     0.0e+00     9.7e-07     1.137e-13     p

  Best performance (actual): 0.9789

3.4. Synthesized filters

The obtained filter \(L(s)\) can then be included in the feedback architecture shown in Figure 13.

The closed-loop transfer functions from \(\hat{x}_1\) to \(\hat{x}\) and from \(\hat{x}_2\) to \(\hat{x}\) corresponding respectively to the sensitivity and complementary sensitivity transfer functions are defined below:

%% Complementary filters
H1 = inv(1 + L);
H2 = 1 - H1;

zpk(H1) =
            (s+3.842)^3 (s+153.6) (s+1.289e05)
  -------------------------------------------------------
  (s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)

zpk(H2) =
         125.61 (s+3358)^2 (s^2 + 46.61s + 813.8)
  -------------------------------------------------------
  (s+1.29e05) (s^2 + 102.1s + 2733) (s^2 + 69.45s + 3272)

The bode plots of the synthesized complementary filters are compared with the upper bounds in Figure 13.

Figure 13: Bode plot of the obtained complementary filters

4. Synthesis of three complementary filters

In this section, the proposed synthesis method of complementary filters is generalized for the synthesis of a set of three complementary filters.

The complete Matlab script for this part is given in Section 5.4.

4.1. Synthesis Architecture

The synthesis objective is to shape three filters that are complementary. This corresponds to the conditions \eqref{eq:obj_three_cf} where \(W_1(s)\), \(W_2(s)\) and \(W_3(s)\) are weighting functions used to specify the maximum wanted magnitude of the three complementary filters.

\begin{equation} \begin{aligned} & |H_1(j\omega)| < \frac{1}{|W_1(j\omega)|}, \quad \forall\omega \\ & |H_2(j\omega)| < \frac{1}{|W_2(j\omega)|}, \quad \forall\omega \\ & |H_3(j\omega)| < \frac{1}{|W_3(j\omega)|}, \quad \forall\omega \\ & H_1(s) + H_2(s) + H_3(s) = 1 \end{aligned} \label{eq:obj_three_cf} \end{equation}

This synthesis can be done by performing the standard \(\mathcal{H}_\infty\) synthesis with on the generalized plant in Figure 14.

Figure 14: Generalized architecture for generating 3 complementary filters

After synthesis, filter \(H_2(s)\) and \(H_3(s)\) are obtained as shown in Figure 14. The last filter \(H_1(s)\) is defined as the complementary of the two others as in \eqref{eq:H1_complementary_of_H2_H3}.

\begin{equation} H_1(s) = 1 - H_2(s) - H_3(s) \label{eq:H1_complementary_of_H2_H3} \end{equation}

4.2. Weights

The three weighting functions are defined as shown below.

%% Design of the Weighting Functions
W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5);
W2 = 0.22*(1 + s/2/pi/1)^2/(sqrt(1e-4) + s/2/pi/1)^2*(1 + s/2/pi/10)^2/(1 + s/2/pi/1000)^2;
W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5);

Their inverse magnitudes are displayed in Figure 15.

Figure 15: Three weighting functions used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters

4.3. H-Infinity Synthesis

The generalized plant in Figure 14 containing the weighting functions is defined below.

%% Generalized plant for the synthesis of 3 complementary filters
P = [W1 -W1 -W1;
     0   W2  0 ;
     0   0   W3;
     1   0   0];

And the standard \(\mathcal{H}_\infty\) synthesis using the hinfsyn command is performed.

%% Standard H-Infinity Synthesis
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');

Resetting value of Gamma min based on D_11, D_12, D_21 terms

Test bounds:      0.1000 <  gamma  <=   1050.0000

  gamma    hamx_eig  xinf_eig  hamy_eig   yinf_eig   nrho_xy   p/f
1.050e+03   3.2e+00   4.5e-13   6.3e-02   -1.2e-11    0.0000    p
  525.050   3.2e+00   1.3e-13   6.3e-02    0.0e+00    0.0000    p
  262.575   3.2e+00   2.1e-12   6.3e-02   -1.5e-13    0.0000    p
  131.337   3.2e+00   1.1e-12   6.3e-02   -7.2e-29    0.0000    p
   65.719   3.2e+00   2.0e-12   6.3e-02    0.0e+00    0.0000    p
   32.909   3.2e+00   7.4e-13   6.3e-02   -5.9e-13    0.0000    p
   16.505   3.2e+00   1.4e-12   6.3e-02    0.0e+00    0.0000    p
    8.302   3.2e+00   1.6e-12   6.3e-02    0.0e+00    0.0000    p
    4.201   3.2e+00   1.6e-12   6.3e-02    0.0e+00    0.0000    p
    2.151   3.2e+00   1.6e-12   6.3e-02    0.0e+00    0.0000    p
    1.125   3.2e+00   2.8e-12   6.3e-02    0.0e+00    0.0000    p
    0.613   3.0e+00 -2.5e+03#  6.3e-02    0.0e+00    0.0000    f
    0.869   3.1e+00 -2.9e+01#  6.3e-02    0.0e+00    0.0000    f
    0.997   3.2e+00   1.9e-12   6.3e-02    0.0e+00    0.0000    p
    0.933   3.1e+00 -6.9e+02#  6.3e-02    0.0e+00    0.0000    f
    0.965   3.1e+00 -3.0e+03#  6.3e-02    0.0e+00    0.0000    f
    0.981   3.1e+00 -8.6e+03#  6.3e-02    0.0e+00    0.0000    f
    0.989   3.2e+00 -2.7e+04#  6.3e-02    0.0e+00    0.0000    f
    0.993   3.2e+00 -5.7e+05#  6.3e-02    0.0e+00    0.0000    f
    0.995   3.2e+00   2.2e-12   6.3e-02    0.0e+00    0.0000    p
    0.994   3.2e+00   1.6e-12   6.3e-02    0.0e+00    0.0000    p
    0.994   3.2e+00   1.0e-12   6.3e-02    0.0e+00    0.0000    p

 Gamma value achieved:     0.9936

The two synthesized filters \(H_2(s)\) and \(H_3(s)\) are defined below: And the third filter \(H_1(s)\) is defined using \eqref{eq:H1_complementary_of_H2_H3}.

%% H1 is defined as the complementary filter of H2 and H3
H1 = 1 - H2 - H3;

The bode plots of the three obtained complementary filters are shown in Figure 16.

Figure 16: The three complementary filters obtained after \(\mathcal{H}_\infty\) synthesis

5. Matlab Scripts

5.1. `1_synthesis_complementary_filters.m`

This scripts corresponds to section 3 of (Dehaeze, Vermat, and Collette 2021).

  1: %% Clear Workspace and Close figures
  2: clear; close all; clc;
  3: 
  4: %% Intialize Laplace variable
  5: s = zpk('s');
  6: 
  7: %% Initialize Frequency Vector
  8: freqs = logspace(-1, 3, 1000);
  9: 
 10: %% Weighting Function Design
 11: % Parameters
 12: n = 3; w0 = 2*pi*10; G0 = 1e-3; G1 = 1e1; Gc = 2;
 13: 
 14: % Formulas
 15: W = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
 16: 
 17: %% Magnitude of the weighting function with parameters
 18: figure;
 19: hold on;
 20: plot(freqs, abs(squeeze(freqresp(W, freqs, 'Hz'))), 'k-');
 21: 
 22: plot([1e-3 1e0], [G0 G0], 'k--', 'LineWidth', 1)
 23: text(1e0, G0, '$\quad G_0$')
 24: 
 25: plot([1e1 1e3], [G1 G1], 'k--', 'LineWidth', 1)
 26: text(1e1,G1,'$G_{\infty}\quad$','HorizontalAlignment', 'right')
 27: 
 28: plot([w0/2/pi w0/2/pi], [1 2*Gc], 'k--', 'LineWidth', 1)
 29: text(w0/2/pi,1,'$\omega_c$','VerticalAlignment', 'top', 'HorizontalAlignment', 'center')
 30: 
 31: plot([w0/2/pi/2  2*w0/2/pi], [Gc Gc], 'k--', 'LineWidth', 1)
 32: text(w0/2/pi/2, Gc, '$G_c \quad$','HorizontalAlignment', 'right')
 33: 
 34: text(w0/5/pi/2, abs(evalfr(W, j*w0/5)), 'Slope: $n \quad$', 'HorizontalAlignment', 'right')
 35: 
 36: text(w0/2/pi, abs(evalfr(W, j*w0)), '$\bullet$', 'HorizontalAlignment', 'center')
 37: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 38: xlabel('Frequency [Hz]'); ylabel('Magnitude');
 39: hold off;
 40: xlim([freqs(1), freqs(end)]);
 41: ylim([5e-4, 20]);
 42: yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]);
 43: 
 44: %% Design of the Weighting Functions
 45: W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
 46: W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);
 47: 
 48: %% Plot of the Weighting function magnitude
 49: figure;
 50: tiledlayout(1, 1, 'TileSpacing', 'None', 'Padding', 'None');
 51: ax1 = nexttile();
 52: hold on;
 53: set(gca,'ColorOrderIndex',1)
 54: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
 55: set(gca,'ColorOrderIndex',2)
 56: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
 57: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 58: xlabel('Frequency [Hz]', 'FontSize', 10); ylabel('Magnitude', 'FontSize', 10);
 59: hold off;
 60: xlim([freqs(1), freqs(end)]);
 61: xticks([0.1, 1, 10, 100, 1000]);
 62: ylim([8e-4, 20]);
 63: yticks([1e-3, 1e-2, 1e-1, 1, 1e1]);
 64: yticklabels({'', '$10^{-2}$', '', '$10^0$', ''});
 65: ax1.FontSize = 9;
 66: leg = legend('location', 'south', 'FontSize', 8);
 67: leg.ItemTokenSize(1) = 18;
 68: 
 69: %% Generalized Plant
 70: P = [W1 -W1;
 71:      0   W2;
 72:      1   0];
 73: 
 74: %% H-Infinity Synthesis
 75: [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
 76: 
 77: %% Define H1 to be the complementary of H2
 78: H1 = 1 - H2;
 79: 
 80: %% Bode plot of the complementary filters
 81: figure;
 82: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
 83: 
 84: % Magnitude
 85: ax1 = nexttile([2, 1]);
 86: hold on;
 87: set(gca,'ColorOrderIndex',1)
 88: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
 89: set(gca,'ColorOrderIndex',2)
 90: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
 91: 
 92: set(gca,'ColorOrderIndex',1)
 93: plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
 94: set(gca,'ColorOrderIndex',2)
 95: plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
 96: hold off;
 97: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 98: set(gca, 'XTickLabel',[]); ylabel('Magnitude');
 99: ylim([8e-4, 20]);
100: yticks([1e-3, 1e-2, 1e-1, 1, 1e1]);
101: yticklabels({'', '$10^{-2}$', '', '$10^0$', ''})
102: leg = legend('location', 'south', 'FontSize', 8, 'NumColumns', 2);
103: leg.ItemTokenSize(1) = 18;
104: 
105: % Phase
106: ax2 = nexttile;
107: hold on;
108: set(gca,'ColorOrderIndex',1)
109: plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
110: set(gca,'ColorOrderIndex',2)
111: plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
112: hold off;
113: set(gca, 'XScale', 'log');
114: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
115: yticks([-180:90:180]);
116: ylim([-180, 200])
117: yticklabels({'-180', '', '0', '', '180'})
118: 
119: linkaxes([ax1,ax2],'x');
120: xlim([freqs(1), freqs(end)]);

5.2. `2_ligo_complementary_filters.m`

This scripts corresponds to section 4 of (Dehaeze, Vermat, and Collette 2021).

  1: %% Clear Workspace and Close figures
  2: clear; close all; clc;
  3: 
  4: %% Intialize Laplace variable
  5: s = zpk('s');
  6: 
  7: %% Initialize Frequency Vector
  8: freqs = logspace(-3, 0, 1000);
  9: 
 10: %% Upper bounds for the complementary filters
 11: figure;
 12: hold on;
 13: set(gca,'ColorOrderIndex',1)
 14: plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
 15: set(gca,'ColorOrderIndex',1)
 16: plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
 17: set(gca,'ColorOrderIndex',1)
 18: plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
 19: set(gca,'ColorOrderIndex',2)
 20: plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
 21: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
 22: xlabel('Frequency [Hz]'); ylabel('Magnitude');
 23: hold off;
 24: xlim([freqs(1), freqs(end)]);
 25: ylim([1e-4, 10]);
 26: leg = legend('location', 'southeast', 'FontSize', 8);
 27: leg.ItemTokenSize(1) = 18;
 28: 
 29: %% Initialized CVX
 30: cvx_startup;
 31: cvx_solver sedumi;
 32: 
 33: %% Frequency vectors
 34: w1 = 0:4.06e-4:0.008;
 35: w2 = 0.008:4.06e-4:0.04;
 36: w3 = 0.04:8.12e-4:0.1;
 37: w4 = 0.1:8.12e-4:0.83;
 38: 
 39: %% Filter order
 40: n = 512;
 41: 
 42: %% Initialization of filter responses
 43: A1 = [ones(length(w1),1),  cos(kron(w1'.*(2*pi),[1:n-1]))];
 44: A2 = [ones(length(w2),1),  cos(kron(w2'.*(2*pi),[1:n-1]))];
 45: A3 = [ones(length(w3),1),  cos(kron(w3'.*(2*pi),[1:n-1]))];
 46: A4 = [ones(length(w4),1),  cos(kron(w4'.*(2*pi),[1:n-1]))];
 47: 
 48: B1 = [zeros(length(w1),1), sin(kron(w1'.*(2*pi),[1:n-1]))];
 49: B2 = [zeros(length(w2),1), sin(kron(w2'.*(2*pi),[1:n-1]))];
 50: B3 = [zeros(length(w3),1), sin(kron(w3'.*(2*pi),[1:n-1]))];
 51: B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))];
 52: 
 53: %% Convex optimization
 54: cvx_begin
 55: 
 56: variable y(n+1,1)
 57: 
 58: % t
 59: maximize(-y(1))
 60: 
 61: for i = 1:length(w1)
 62:     norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3;
 63: end
 64: 
 65: for  i = 1:length(w2)
 66:     norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3*(2*pi*w2(i)/(0.008*2*pi))^3;
 67: end
 68: 
 69: for i = 1:length(w3)
 70:     norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3;
 71: end
 72: 
 73: for i = 1:length(w4)
 74:     norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1);
 75: end
 76: 
 77: cvx_end
 78: 
 79: h = y(2:end);
 80: 
 81: %% Combine the frequency vectors to form the obtained filter
 82: w = [w1 w2 w3 w4];
 83: H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
 84: 
 85: %% Bode plot of the obtained complementary filters
 86: figure;
 87: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
 88: 
 89: % Magnitude
 90: ax1 = nexttile([2, 1]);
 91: hold on;
 92: set(gca,'ColorOrderIndex',1)
 93: plot(w, abs(1-H), '-', 'DisplayName', '$L_1$');
 94: plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
 95: 
 96: set(gca,'ColorOrderIndex',2)
 97: plot(w, abs(H), '-', 'DisplayName', '$H_1$');
 98: plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
 99: plot([0.008 0.04], [8e-3, 1],       'k--', 'HandleVisibility', 'off');
100: plot([0.04 0.1], [3, 3],            'k--', 'HandleVisibility', 'off');
101: 
102: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
103: set(gca, 'XTickLabel',[]); ylabel('Magnitude');
104: hold off;
105: ylim([5e-3, 10]);
106: leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
107: leg.ItemTokenSize(1) = 16;
108: 
109: % Phase
110: ax2 = nexttile;
111: hold on;
112: plot(w, 180/pi*unwrap(angle(1-H)), '-');
113: plot(w, 180/pi*unwrap(angle(H)), '-');
114: hold off;
115: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
116: set(gca, 'XScale', 'log');
117: yticks([-360:180:180]); ylim([-380, 200]);
118: 
119: linkaxes([ax1,ax2],'x');
120: xlim([1e-3, 1]);
121: 
122: %% Design of the weight for the high pass filter
123: w1 = 2*pi*0.008; x1 = 0.35;
124: w2 = 2*pi*0.04;  x2 = 0.5;
125: w3 = 2*pi*0.05;  x3 = 0.5;
126: 
127: % Slope of +3 from w1
128: wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1);
129: % Little bump from w2 to w3
130: wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1);
131: % No Slope at high frequencies
132: wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1);
133: % Little bump between w2 and w3
134: w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1;
135: wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n;
136: 
137: wH = 1/wH;
138: wH = minreal(ss(wH));
139: 
140: %% Design of the weight for the low pass filter
141: n = 20; % Filter order
142: Rp = 1; % Peak to peak passband ripple
143: Wp = 2*pi*0.102; % Edge frequency
144: 
145: % Chebyshev Type I filter design
146: [b,a] = cheby1(n, Rp, Wp, 'high', 's');
147: wL = 0.04*tf(a, b);
148: 
149: wL = 1/wL;
150: wL = minreal(ss(wL));
151: 
152: %% Magnitude of the designed Weights and initial specifications
153: figure;
154: hold on;
155: set(gca,'ColorOrderIndex',1);
156: plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|W_L|^{-1}$');
157: plot([0.1, 10], [0.045, 0.045],     'k:',  'DisplayName', 'Spec. on $L_1$');
158: 
159: set(gca,'ColorOrderIndex',2);
160: plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|W_H|^{-1}$');
161: plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
162: plot([0.008 0.04], [8e-3, 1],       'k--', 'HandleVisibility', 'off');
163: plot([0.04 0.1], [3, 3],            'k--', 'HandleVisibility', 'off');
164: 
165: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
166: xlabel('Frequency [Hz]'); ylabel('Magnitude');
167: hold off;
168: xlim([freqs(1), freqs(end)]);
169: ylim([5e-3, 10]);
170: leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
171: leg.ItemTokenSize(1) = 16;
172: 
173: %% Generalized plant for the H-infinity Synthesis
174: P = [0   wL;
175:      wH -wH;
176:      1   0];
177: 
178: %% Standard H-Infinity synthesis
179: [Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
180: 
181: %% High pass filter as the complementary of the low pass filter
182: Hh = 1 - Hl;
183: 
184: %% Minimum realization of the filters
185: Hh = minreal(Hh);
186: Hl = minreal(Hl);
187: 
188: %% Bode plot of the obtained filters and comparison with the upper bounds
189: figure;
190: hold on;
191: set(gca,'ColorOrderIndex',1);
192: plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$L_1^\prime$');
193: plot([0.1, 10], [0.045, 0.045],     'k:',  'DisplayName', 'Spec. on $L_1$');
194: 
195: set(gca,'ColorOrderIndex',2);
196: plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_1^\prime$');
197: plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
198: plot([0.008 0.04], [8e-3, 1],       'k--', 'HandleVisibility', 'off');
199: plot([0.04 0.1], [3, 3],            'k--', 'HandleVisibility', 'off');
200: 
201: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
202: xlabel('Frequency [Hz]'); ylabel('Magnitude');
203: hold off;
204: xlim([freqs(1), freqs(end)]);
205: ylim([5e-3, 10]);
206: leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
207: leg.ItemTokenSize(1) = 16;
208: 
209: %% Comparison of the complementary filters obtained with H-infinity and with CVX
210: figure;
211: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
212: 
213: % Magnitude
214: ax1 = nexttile([2, 1]);
215: hold on;
216: set(gca,'ColorOrderIndex',1);
217: plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ...
218:      'DisplayName', '$L_1(s)$ - $\mathcal{H}_\infty$');
219: set(gca,'ColorOrderIndex',2);
220: plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ...
221:      'DisplayName', '$H_1(s)$ - $\mathcal{H}_\infty$');
222: 
223: set(gca,'ColorOrderIndex',1);
224: plot(w, abs(1-H), '--', ...
225:      'DisplayName', '$L_1(s)$ - FIR');
226: set(gca,'ColorOrderIndex',2);
227: plot(w, abs(H), '--', ...
228:      'DisplayName', '$H_1(s)$ - FIR');
229: hold off;
230: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
231: ylabel('Magnitude');
232: set(gca, 'XTickLabel',[]);
233: ylim([5e-3, 10]);
234: leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
235: leg.ItemTokenSize(1) = 16;
236: 
237: % Phase
238: ax2 = nexttile;
239: hold on;
240: set(gca,'ColorOrderIndex',1);
241: plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hl, freqs, 'Hz')))), '-');
242: set(gca,'ColorOrderIndex',2);
243: plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
244: 
245: set(gca,'ColorOrderIndex',1);
246: plot(w, 180/pi*unwrap(angle(1-H)), '--');
247: set(gca,'ColorOrderIndex',2);
248: plot(w, 180/pi*unwrap(angle(H)), '--');
249: set(gca, 'XScale', 'log');
250: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
251: hold off;
252: yticks([-360:180:180]); ylim([-380, 200]);
253: 
254: linkaxes([ax1,ax2],'x');
255: xlim([freqs(1), freqs(end)]);

5.3. `3_closed_loop_complementary_filters.m`

This scripts corresponds to section 5.1 of (Dehaeze, Vermat, and Collette 2021).

 1: %% Clear Workspace and Close figures
 2: clear; close all; clc;
 3: 
 4: %% Intialize Laplace variable
 5: s = zpk('s');
 6: 
 7: %% Initialize Frequency Vector
 8: freqs = logspace(-1, 3, 1000);
 9: 
10: %% Design of the Weighting Functions
11: W1 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.45);
12: W2 = generateWF('n', 2, 'w0', 2*pi*10, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.45);
13: 
14: %% Generalized plant for "closed-loop" complementary filter synthesis
15: P = [ W1 0   1;
16:      -W1 W2 -1];
17: 
18: %% Standard H-Infinity Synthesis
19: [L, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
20: 
21: %% Complementary filters
22: H1 = inv(1 + L);
23: H2 = 1 - H1;
24: 
25: %% Bode plot of the obtained Complementary filters with upper-bounds
26: freqs = logspace(-1, 3, 1000);
27: figure;
28: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
29: 
30: % Magnitude
31: ax1 = nexttile([2, 1]);
32: hold on;
33: set(gca,'ColorOrderIndex',1)
34: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
35: set(gca,'ColorOrderIndex',2)
36: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
37: 
38: set(gca,'ColorOrderIndex',1)
39: plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
40: set(gca,'ColorOrderIndex',2)
41: plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
42: 
43: plot(freqs, abs(squeeze(freqresp(L, freqs, 'Hz'))), 'k--', 'DisplayName', '$|L|$');
44: hold off;
45: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
46: set(gca, 'XTickLabel',[]); ylabel('Magnitude');
47: ylim([1e-3, 1e3]);
48: yticks([1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3]);
49: yticklabels({'', '$10^{-2}$', '', '$10^0$', '', '$10^2$', ''});
50: leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 3);
51: leg.ItemTokenSize(1) = 18;
52: 
53: % Phase
54: ax2 = nexttile;
55: hold on;
56: set(gca,'ColorOrderIndex',1)
57: plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
58: set(gca,'ColorOrderIndex',2)
59: plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
60: hold off;
61: set(gca, 'XScale', 'log');
62: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
63: yticks([-180:90:180]);
64: ylim([-180, 200])
65: yticklabels({'-180', '', '0', '', '180'})
66: 
67: linkaxes([ax1,ax2],'x');
68: xlim([freqs(1), freqs(end)]);

5.4. `4_three_complementary_filters.m`

This scripts corresponds to section 5.2 of (Dehaeze, Vermat, and Collette 2021).

 1: %% Clear Workspace and Close figures
 2: clear; close all; clc;
 3: 
 4: %% Intialize Laplace variable
 5: s = zpk('s');
 6: 
 7: freqs = logspace(-2, 3, 1000);
 8: 
 9: %% Design of the Weighting Functions
10: W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'Ginf', 1000, 'Gc', 0.5);
11: W2 = 0.22*(1 + s/2/pi/1)^2/(sqrt(1e-4) + s/2/pi/1)^2*(1 + s/2/pi/10)^2/(1 + s/2/pi/1000)^2;
12: W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'Ginf', 1/10, 'Gc', 0.5);
13: 
14: %% Inverse magnitude of the weighting functions
15: figure;
16: hold on;
17: set(gca,'ColorOrderIndex',1)
18: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
19: set(gca,'ColorOrderIndex',2)
20: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
21: set(gca,'ColorOrderIndex',3)
22: plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
23: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
24: xlabel('Frequency [Hz]'); ylabel('Magnitude');
25: hold off;
26: xlim([freqs(1), freqs(end)]); ylim([2e-4, 1.3e1])
27: leg = legend('location', 'northeast', 'FontSize', 8);
28: leg.ItemTokenSize(1) = 18;
29: 
30: %% Generalized plant for the synthesis of 3 complementary filters
31: P = [W1 -W1 -W1;
32:      0   W2  0 ;
33:      0   0   W3;
34:      1   0   0];
35: 
36: %% Standard H-Infinity Synthesis
37: [H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
38: 
39: %% Synthesized H2 and H3 filters
40: H2 = tf(H(1));
41: H3 = tf(H(2));
42: 
43: %% H1 is defined as the complementary filter of H2 and H3
44: H1 = 1 - H2 - H3;
45: 
46: %% Bode plot of the obtained complementary filters
47: figure;
48: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
49: 
50: % Magnitude
51: ax1 = nexttile([2, 1]);
52: hold on;
53: set(gca,'ColorOrderIndex',1)
54: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
55: set(gca,'ColorOrderIndex',2)
56: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
57: set(gca,'ColorOrderIndex',3)
58: plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
59: set(gca,'ColorOrderIndex',1)
60: plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
61: set(gca,'ColorOrderIndex',2)
62: plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
63: set(gca,'ColorOrderIndex',3)
64: plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$');
65: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
66: hold off;
67: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
68: ylabel('Magnitude');
69: set(gca, 'XTickLabel',[]);
70: ylim([1e-4, 20]);
71: leg = legend('location', 'northeast', 'FontSize', 8);
72: leg.ItemTokenSize(1) = 18;
73: 
74: % Phase
75: ax2 = nexttile;
76: hold on;
77: set(gca,'ColorOrderIndex',1)
78: plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))));
79: set(gca,'ColorOrderIndex',2)
80: plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))));
81: set(gca,'ColorOrderIndex',3)
82: plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz'))));
83: hold off;
84: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
85: set(gca, 'XScale', 'log');
86: yticks([-180:90:180]); ylim([-220, 220]);
87: 
88: linkaxes([ax1,ax2],'x');
89: xlim([freqs(1), freqs(end)]);

5.5. `generateWF`: Generate Weighting Functions

This function is used to easily generate weighting functions from classical requirements.

function [W] = generateWF(args)
% createWeight -
%
% Syntax: [W] = generateWeight(args)
%
% Inputs:
%    - n  - Weight Order (integer)
%    - G0 - Low frequency Gain
%    - G1 - High frequency Gain
%    - Gc - Gain of the weight at frequency w0
%    - w0 - Frequency at which |W(j w0)| = Gc [rad/s]
%
% Outputs:
%    - W - Generated Weighting Function

%% Argument validation
arguments
    args.n    (1,1) double {mustBeInteger, mustBePositive} = 1
    args.G0   (1,1) double {mustBeNumeric, mustBePositive} = 0.1
    args.Ginf (1,1) double {mustBeNumeric, mustBePositive} = 10
    args.Gc   (1,1) double {mustBeNumeric, mustBePositive} = 1
    args.w0   (1,1) double {mustBeNumeric, mustBePositive} = 1
end

% Verification of correct relation between G0, Gc and Ginf
mustBeBetween(args.G0, args.Gc, args.Ginf);

%% Initialize the Laplace variable
s = zpk('s');

%% Create the weighting function according to formula
W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ...
      (args.G0/args.Gc)^(1/args.n))/...
     ((1/args.Ginf)^(1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ...
      (1/args.Gc)^(1/args.n)))^args.n;

%% Custom validation function
function mustBeBetween(a,b,c)
    if ~((a > b && b > c) || (c > b && b > a))
        eid = 'createWeight:inputError';
        msg = 'Gc should be between G0 and Ginf.';
        throwAsCaller(MException(eid,msg))
    end

5.6. `generateCF`: Generate Complementary Filters

This function is used to easily synthesize a set of two complementary filters using the \(\mathcal{H}_\infty\) synthesis.

function [H1, H2] = generateCF(W1, W2, args)
% createWeight -
%
% Syntax: [H1, H2] = generateCF(W1, W2, args)
%
% Inputs:
%    - W1 - Weighting Function for H1
%    - W2 - Weighting Function for H2
%    - args:
%      - method  - H-Infinity solver ('lmi' or 'ric')
%      - display - Display synthesis results ('on' or 'off')
%
% Outputs:
%    - H1 - Generated H1 Filter
%    - H2 - Generated H2 Filter

%% Argument validation
arguments
    W1
    W2
    args.method  char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric'
    args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on'
end

%% The generalized plant is defined
P = [W1 -W1;
     0   W2;
     1   0];

%% The standard H-infinity synthesis is performed
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display);

%% H1 is defined as the complementary of H2
H1 = 1 - H2;

Bibliography

Dehaeze, Thomas, Mohit Vermat, and Christophe Collette. 2021. “A New Method of Designing Complementary Filters for Sensor Fusion Using the H-Infinity Synthesis.” Mechanical Systems and Signal Processing, November.

Grant, Michael, and Stephen Boyd. 2014. “CVX: Matlab Software for Disciplined Convex Programming, Version 2.1.”

Hua, Wensheng, R Adhikari, Daniel B DeBra, Joseph A Giaime, Giles Dominic Hammond, C Hardham, Mike Hennessy, Jonathan P How, et al. 2004. “Low-Frequency Active Vibration Isolation for Advanced LIGO.” In Gravitational Wave and Particle Astrophysics Detectors, 5500:194–205. International Society for Optics and Photonics.

Hua, Wensheng, B. Debra, T. Hardham, T. Lantz, and A. Giaime. 2004. “Polyphase FIR Complementary Filters for Control Systems.” In Proceedings of ASPE Spring Topical Meeting on Control of Precision Systems, 109–14.

MATLAB. 2020. Version 9.9.0 (R2020b). Natick, Massachusetts: The MathWorks Inc.

Sturm, Jos F. 1999. “Using SeDuMi 1.02, a Matlab Toolbox for Optimization over Symmetric Cones.” Optimization Methods and Software 11 (1-4). Taylor & Francis:625–53. https://doi.org/10.1080/10556789908805766.