# A new method of designing complementary filters for sensor fusion using the $\mathcal{H}_\infty$ synthesis - Matlab Computation
Author: Dehaeze Thomas

UP | HOME

# A new method of designing complementary filters for sensor fusion using the \(\mathcal{H}_\infty\) synthesis - Matlab Computation

## Table of Contents

- 1. H-Infinity synthesis of complementary filters
  - 1.1. Synthesis Architecture
  - 1.2. Design of Weighting Function - Proposed formula
  - 1.3. Weighting functions for the design of two complementary filters
  - 1.4. Synthesis of the complementary filters
  - 1.5. Obtained Complementary Filters
- 2. Design of complementary filters used in the Active Vibration Isolation System at the LIGO
  - 2.1. Specifications
  - 2.2. FIR Filter
  - 2.3. Weighting function design
  - 2.4. Synthesis of the complementary filters
  - 2.5. Comparison of the FIR filters and synthesized filters
- 3. “Closed-Loop” complementary filters
  - 3.1. Weighting Function design
  - 3.2. Generalized plant
  - 3.3. Synthesis of the closed-loop complementary filters
  - 3.4. Synthesized filters
- 4. Synthesis of three complementary filters
  - 4.1. Synthesis Architecture
  - 4.2. Weights
  - 4.3. H-Infinity Synthesis
- 5. Matlab Scripts
  - 5.1. `1_synthesis_complementary_filters.m`
  - 5.2. `2_ligo_complementary_filters.m`
  - 5.3. `3_closed_loop_complementary_filters.m`
  - 5.4. `4_three_complementary_filters.m`
  - 5.5. `generateWF`: Generate Weighting Functions
  - 5.6. `generateCF`: Generate Complementary Filters

This report is also available as a pdf.

The present document is a companion file for the journal paper (Dehaeze, Vermat, and Collette 2004). 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 2004).

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

### 5.2. `2_ligo_complementary_filters.m`

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

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

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

### 5.4. `4_three_complementary_filters.m`

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

      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: addpath('./src');
     10:
     11: % Weights
     12: % First we define the weights.
     13:
     14: %% Design of the Weighting Functions
     15: W1 = generateWF('n', 2, 'w0', 2*pi*1, 'G0', 1/10, 'G1', 1000, 'Gc', 0.5);
     16: 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;
     17: W3 = generateWF('n', 3, 'w0', 2*pi*10, 'G0', 1000, 'G1', 1/10, 'Gc', 0.5);
     18:
     19: %% Inverse magnitude of the weighting functions
     20: figure;
     21: hold on;
     22: set(gca,'ColorOrderIndex',1)
     23: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
     24: set(gca,'ColorOrderIndex',2)
     25: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
     26: set(gca,'ColorOrderIndex',3)
     27: plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
     28: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
     29: xlabel('Frequency [Hz]'); ylabel('Magnitude');
     30: hold off;
     31: xlim([freqs(1), freqs(end)]); ylim([2e-4, 1.3e1])
     32: leg = legend('location', 'northeast', 'FontSize', 8);
     33: leg.ItemTokenSize(1) = 18;
     34:
     35: % H-Infinity Synthesis
     36: % Then we create the generalized plant =P=.
     37:
     38: %% Generalized plant for the synthesis of 3 complementary filters
     39: P = [W1 -W1 -W1;
     40:      0   W2  0 ;
     41:      0   0   W3;
     42:      1   0   0];
     43:
     44:
     45:
     46: % And we do the $\mathcal{H}_\infty$ synthesis.
     47:
     48: %% Standard H-Infinity Synthesis
     49: [H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
     50:
     51: % Obtained Complementary Filters
     52: % The obtained filters are:
     53:
     54: %%
     55: H2 = tf(H(1));
     56: H3 = tf(H(2));
     57: H1 = 1 - H2 - H3;
     58:
     59: %% Bode plot of the obtained complementary filters
     60: figure;
     61: tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
     62:
     63: % Magnitude
     64: ax1 = nexttile([2, 1]);
     65: hold on;
     66: set(gca,'ColorOrderIndex',1)
     67: plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
     68: set(gca,'ColorOrderIndex',2)
     69: plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
     70: set(gca,'ColorOrderIndex',3)
     71: plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
     72: set(gca,'ColorOrderIndex',1)
     73: plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
     74: set(gca,'ColorOrderIndex',2)
     75: plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
     76: set(gca,'ColorOrderIndex',3)
     77: plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$');
     78: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
     79: hold off;
     80: set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
     81: ylabel('Magnitude');
     82: set(gca, 'XTickLabel',[]);
     83: ylim([1e-4, 20]);
     84: leg = legend('location', 'northeast', 'FontSize', 8);
     85: leg.ItemTokenSize(1) = 18;
     86:
     87: % Phase
     88: ax2 = nexttile;
     89: hold on;
     90: set(gca,'ColorOrderIndex',1)
     91: plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))));
     92: set(gca,'ColorOrderIndex',2)
     93: plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))));
     94: set(gca,'ColorOrderIndex',3)
     95: plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz'))));
     96: hold off;
     97: xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
     98: set(gca, 'XScale', 'log');
     99: yticks([-180:90:180]); ylim([-220, 220]);
    100:
    101: linkaxes([ax1,ax2],'x');
    102: 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. 2004. “A New Method of Designing Complementary Fil.” In Gravitational Wave and Particle Astrophysics Detectors, 5500:194–205. International Society for Optics and Photonics.
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.

Author: Dehaeze Thomas

Created: 2021-09-01 mer. 14:58