tdehaeze/dehaeze21_desig_compl_filte
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Rework all the Matlab computation

Browse Source
This commit is contained in:
Thomas Dehaeze 2020-10-26 21:35:47 +01:00
parent 7b870ae344
commit 5d6df8df40
54 changed files with 14350 additions and 1997 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
matlab/figs-tikz Symbolic link
View File

@ -0,0 +1 @@
../tikz/figs/

BIN
matlab/figs/comp_fir_ligo_hinf.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/comp_fir_ligo_hinf.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 130 KiB

BIN
matlab/figs/fir_filter_ligo.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/fir_filter_ligo.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 80 KiB

BIN
matlab/figs/hinf_filters_results.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/hinf_filters_results.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 134 KiB

BIN
matlab/figs/hinf_synthesis_ligo_results.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 75 KiB

BIN
matlab/figs/ligo_specifications.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/ligo_specifications.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 69 KiB

BIN
matlab/figs/ligo_weights.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 69 KiB

2111
matlab/figs/three_complementary_filters_results.pdf Normal file
View File

File diff suppressed because it is too large Load Diff

BIN
matlab/figs/three_complementary_filters_results.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 183 KiB

BIN
matlab/figs/three_weighting_functions.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/three_weighting_functions.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 101 KiB

BIN
matlab/figs/weights_W1_W2.pdf Normal file
View File

Binary file not shown.

BIN
matlab/figs/weights_W1_W2.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 76 KiB

802
matlab/index.html Normal file
View File

@ -0,0 +1,802 @@
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2020-10-26 lun. 18:59 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<title>Complementary Filters Shaping Using \(\mathcal{H}_\infty\) Synthesis - Matlab Computation</title>
<meta name="generator" content="Org mode" />
<meta name="author" content="Thomas Dehaeze" />
<link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
<link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
<script src="../js/jquery.min.js"></script>
<script src="../js/bootstrap.min.js"></script>
<script src="../js/jquery.stickytableheaders.min.js"></script>
<script src="../js/readtheorg.js"></script>
<script>MathJax = {
tex: {
tags: 'ams',
macros: {bm: ["\\boldsymbol{#1}",1],}
}
};
</script>
<script type="text/javascript" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
</head>
<body>
<div id="org-div-home-and-up">
<a accesskey="h" href="../index.html"> UP </a>
|
<a accesskey="H" href="../index.html"> HOME </a>
</div><div id="content">
<h1 class="title">Complementary Filters Shaping Using \(\mathcal{H}_\infty\) Synthesis - Matlab Computation</h1>
<div id="table-of-contents">
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#orgcd571e1">1. H-Infinity synthesis of complementary filters</a>
<ul>
<li><a href="#orgea65f42">1.1. Synthesis Architecture</a></li>
<li><a href="#org8a8a446">1.2. Design of Weighting Function</a></li>
<li><a href="#orgfe22503">1.3. H-Infinity Synthesis</a></li>
<li><a href="#orgf9471db">1.4. Obtained Complementary Filters</a></li>
</ul>
</li>
<li><a href="#org0e59118">2. Generating 3 complementary filters</a>
<ul>
<li><a href="#orgfbd127e">2.1. Theory</a></li>
<li><a href="#orgb689a85">2.2. Weights</a></li>
<li><a href="#org24d1ff1">2.3. H-Infinity Synthesis</a></li>
<li><a href="#orgef81b1e">2.4. Obtained Complementary Filters</a></li>
</ul>
</li>
<li><a href="#org2c0b99c">3. Implement complementary filters for LIGO</a>
<ul>
<li><a href="#org6ea0071">3.1. Specifications</a></li>
<li><a href="#org8f811fb">3.2. FIR Filter</a></li>
<li><a href="#org9a34a27">3.3. Weights</a></li>
<li><a href="#orgc9f5d68">3.4. H-Infinity Synthesis</a></li>
<li><a href="#orgd83a581">3.5. Compare FIR and H-Infinity Filters</a></li>
</ul>
</li>
</ul>
</div>
</div>
<p>
In this document, the design of complementary filters is studied.
</p>
<p>
One use of complementary filter is described below:
</p>
<blockquote>
<p>
The basic idea of a complementary filter involves taking two or more sensors, filtering out unreliable frequencies for each sensor, and combining the filtered outputs to get a better estimate throughout the entire bandwidth of the system.
To achieve this, the sensors included in the filter should complement one another by performing better over specific parts of the system bandwidth.
</p>
</blockquote>
<p>
This document is divided into several sections:
</p>
<ul class="org-ul">
<li>in section <a href="#org954a80e">1</a>, the \(\mathcal{H}_\infty\) synthesis is used for generating two complementary filters</li>
<li>in section <a href="#orga319f0c">2</a>, a method using the \(\mathcal{H}_\infty\) synthesis is proposed to shape three of more complementary filters</li>
<li>in section <a href="#orgb69dc24">3</a>, the \(\mathcal{H}_\infty\) synthesis is used and compared with FIR complementary filters used for LIGO</li>
</ul>
<div class="note" id="org3c0f990">
<p>
Add the Matlab code use to obtain the results presented in the paper are accessible <a href="matlab.zip">here</a> and presented below.
</p>
</div>
<div id="outline-container-orgcd571e1" class="outline-2">
<h2 id="orgcd571e1"><span class="section-number-2">1</span> H-Infinity synthesis of complementary filters</h2>
<div class="outline-text-2" id="text-1">
<p>
<a id="org954a80e"></a>
</p>
<div class="note" id="orgca3b1ba">
<p>
The Matlab file corresponding to this section is accessible <a href="matlab/h_inf_synthesis_complementary_filters.m">here</a>.
</p>
</div>
</div>
<div id="outline-container-orgea65f42" class="outline-3">
<h3 id="orgea65f42"><span class="section-number-3">1.1</span> Synthesis Architecture</h3>
<div class="outline-text-3" id="text-1-1">
<p>
We here synthesize two complementary filters using the \(\mathcal{H}_\infty\) synthesis.
The goal is to specify upper bounds on the norms of the two complementary filters \(H_1(s)\) and \(H_2(s)\) while ensuring their complementary property (\(H_1(s) + H_2(s) = 1\)).
</p>
<p>
In order to do so, we use the generalized plant shown on figure <a href="#org46317e2">1</a> where \(W_1(s)\) and \(W_2(s)\) are weighting transfer functions that will be used to shape \(H_1(s)\) and \(H_2(s)\) respectively.
</p>
<div id="org46317e2" class="figure">
<p><img src="figs-tikz/h_infinity_robust_fusion.png" alt="h_infinity_robust_fusion.png" />
</p>
<p><span class="figure-number">Figure 1: </span>\(\mathcal{H}_\infty\) synthesis of the complementary filters</p>
</div>
<p>
The \(\mathcal{H}_\infty\) synthesis applied on this generalized plant will give a transfer function \(H_2\) (figure <a href="#org46317e2">1</a>) such that the \(\mathcal{H}_\infty\) norm of the transfer function from \(w\) to \([z_1,\ z_2]\) is less than one:
\[ \left\| \begin{array}{c} (1 - H_2(s)) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \]
</p>
<p>
Thus, if the above condition is verified, we can define \(H_1(s) = 1 - H_2(s)\) and we have that:
\[ \left\| \begin{array}{c} H_1(s) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \]
Which is almost (with an maximum error of \(\sqrt{2}\)) equivalent to:
</p>
\begin{align*}
|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{align*}
<p>
We then see that \(W_1(s)\) and \(W_2(s)\) can be used to shape both \(H_1(s)\) and \(H_2(s)\) while ensuring their complementary property by the definition of \(H_1(s) = 1 - H_2(s)\).
</p>
</div>
</div>
<div id="outline-container-org8a8a446" class="outline-3">
<h3 id="org8a8a446"><span class="section-number-3">1.2</span> Design of Weighting Function</h3>
<div class="outline-text-3" id="text-1-2">
<p>
A formula is proposed to help the design of the weighting functions:
</p>
\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
\end{equation}
<p>
The parameters permits to specify:
</p>
<ul class="org-ul">
<li>the low frequency gain: \(G_0 = lim_{\omega \to 0} |W(j\omega)|\)</li>
<li>the high frequency gain: \(G_\infty = lim_{\omega \to \infty} |W(j\omega)|\)</li>
<li>the absolute gain at \(\omega_0\): \(G_c = |W(j\omega_0)|\)</li>
<li>the absolute slope between high and low frequency: \(n\)</li>
</ul>
<p>
The general shape of a weighting function generated using the formula is shown in figure <a href="#org4ee6999">2</a>.
</p>
<div id="org4ee6999" class="figure">
<p><img src="figs-tikz/weight_formula.png" alt="weight_formula.png" />
</p>
<p><span class="figure-number">Figure 2: </span>Amplitude of the proposed formula for the weighting functions</p>
</div>
<div class="org-src-container">
<pre class="src src-matlab">n = 2; w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>11; G0 = 1<span class="org-type">/</span>10; G1 = 1000; Gc = 1<span class="org-type">/</span>2;
W1 = (((1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n))<span class="org-type">/</span>((1<span class="org-type">/</span>G1)<span class="org-type">^</span>(1<span class="org-type">/</span>n)<span class="org-type">*</span>(1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (1<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n)))<span class="org-type">^</span>n;
n = 3; w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>10; G0 = 1000; G1 = 0.1; Gc = 1<span class="org-type">/</span>2;
W2 = (((1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n))<span class="org-type">/</span>((1<span class="org-type">/</span>G1)<span class="org-type">^</span>(1<span class="org-type">/</span>n)<span class="org-type">*</span>(1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (1<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n)))<span class="org-type">^</span>n;
</pre>
</div>
<div id="org79474f8" class="figure">
<p><img src="figs/weights_W1_W2.png" alt="weights_W1_W2.png" />
</p>
<p><span class="figure-number">Figure 3: </span>Weights on the complementary filters \(W_1\) and \(W_2\) and the associated performance weights</p>
</div>
</div>
</div>
<div id="outline-container-orgfe22503" class="outline-3">
<h3 id="orgfe22503"><span class="section-number-3">1.3</span> H-Infinity Synthesis</h3>
<div class="outline-text-3" id="text-1-3">
<p>
We define the generalized plant \(P\) on matlab.
</p>
<div class="org-src-container">
<pre class="src src-matlab">P = [W1 <span class="org-type">-</span>W1;
0 W2;
1 0];
</pre>
</div>
<p>
And we do the \(\mathcal{H}_\infty\) synthesis using the <code>hinfsyn</code> command.
</p>
<div class="org-src-container">
<pre class="src src-matlab">[H2, <span class="org-type">~</span>, gamma, <span class="org-type">~</span>] = hinfsyn(P, 1, 1,<span class="org-string">'TOLGAM'</span>, 0.001, <span class="org-string">'METHOD'</span>, <span class="org-string">'ric'</span>, <span class="org-string">'DISPLAY'</span>, <span class="org-string">'on'</span>);
</pre>
</div>
<pre class="example" id="org9fcd456">
[H2, ~, 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.1000 &lt; gamma &lt;= 1050.0000
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.050e+03 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
525.050 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
262.575 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
131.337 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
65.719 2.8e+01 2.4e-07 4.1e+00 -9.5e-14 0.0000 p
32.909 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
16.505 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
8.302 2.8e+01 2.4e-07 4.1e+00 -7.2e-14 0.0000 p
4.201 2.8e+01 2.4e-07 4.1e+00 -2.5e-25 0.0000 p
2.151 2.7e+01 2.4e-07 4.1e+00 -3.8e-14 0.0000 p
1.125 2.6e+01 2.4e-07 4.1e+00 -5.4e-24 0.0000 p
0.613 2.3e+01 -3.7e+01# 4.1e+00 0.0e+00 0.0000 f
0.869 2.6e+01 -3.7e+02# 4.1e+00 0.0e+00 0.0000 f
0.997 2.6e+01 -1.1e+04# 4.1e+00 0.0e+00 0.0000 f
1.061 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.029 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.013 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.005 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.001 2.6e+01 -3.1e+04# 4.1e+00 -3.8e-14 0.0000 f
1.003 2.6e+01 -2.8e+05# 4.1e+00 0.0e+00 0.0000 f
1.004 2.6e+01 2.4e-07 4.1e+00 -5.8e-24 0.0000 p
1.004 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
Gamma value achieved: 1.0036
</pre>
<p>
We then define the high pass filter \(H_1 = 1 - H_2\). The bode plot of both \(H_1\) and \(H_2\) is shown on figure <a href="#orgbc23afe">4</a>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">H1 = 1 <span class="org-type">-</span> H2;
</pre>
</div>
</div>
</div>
<div id="outline-container-orgf9471db" class="outline-3">
<h3 id="orgf9471db"><span class="section-number-3">1.4</span> Obtained Complementary Filters</h3>
<div class="outline-text-3" id="text-1-4">
<p>
The obtained complementary filters are shown on figure <a href="#orgbc23afe">4</a>.
</p>
<div id="orgbc23afe" class="figure">
<p><img src="figs/hinf_filters_results.png" alt="hinf_filters_results.png" />
</p>
<p><span class="figure-number">Figure 4: </span>Obtained complementary filters using \(\mathcal{H}_\infty\) synthesis</p>
</div>
</div>
</div>
</div>
<div id="outline-container-org0e59118" class="outline-2">
<h2 id="org0e59118"><span class="section-number-2">2</span> Generating 3 complementary filters</h2>
<div class="outline-text-2" id="text-2">
<p>
<a id="orga319f0c"></a>
</p>
<div class="note" id="org2ba2674">
<p>
The Matlab file corresponding to this section is accessible <a href="matlab/three_comp_filters.m">here</a>.
</p>
</div>
</div>
<div id="outline-container-orgfbd127e" class="outline-3">
<h3 id="orgfbd127e"><span class="section-number-3">2.1</span> Theory</h3>
<div class="outline-text-3" id="text-2-1">
<p>
We want:
</p>
\begin{align*}
& |H_1(j\omega)| < 1/|W_1(j\omega)|, \quad \forall\omega\\
& |H_2(j\omega)| < 1/|W_2(j\omega)|, \quad \forall\omega\\
& |H_3(j\omega)| < 1/|W_3(j\omega)|, \quad \forall\omega\\
& H_1(s) + H_2(s) + H_3(s) = 1
\end{align*}
<p>
For that, we use the \(\mathcal{H}_\infty\) synthesis with the architecture shown on figure <a href="#org353bdb6">5</a>.
</p>
<div id="org353bdb6" class="figure">
<p><img src="figs-tikz/comp_filter_three_hinf.png" alt="comp_filter_three_hinf.png" />
</p>
<p><span class="figure-number">Figure 5: </span>Generalized architecture for generating 3 complementary filters</p>
</div>
<p>
The \(\mathcal{H}_\infty\) objective is:
</p>
\begin{align*}
& |(1 - H_2(j\omega) - H_3(j\omega)) W_1(j\omega)| < 1, \quad \forall\omega\\
& |H_2(j\omega) W_2(j\omega)| < 1, \quad \forall\omega\\
& |H_3(j\omega) W_3(j\omega)| < 1, \quad \forall\omega\\
\end{align*}
<p>
And thus if we choose \(H_1 = 1 - H_2 - H_3\) we have solved the problem.
</p>
</div>
</div>
<div id="outline-container-orgb689a85" class="outline-3">
<h3 id="orgb689a85"><span class="section-number-3">2.2</span> Weights</h3>
<div class="outline-text-3" id="text-2-2">
<p>
First we define the weights.
</p>
<div class="org-src-container">
<pre class="src src-matlab">n = 2; w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>1; G0 = 1<span class="org-type">/</span>10; G1 = 1000; Gc = 1<span class="org-type">/</span>2;
W1 = (((1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n))<span class="org-type">/</span>((1<span class="org-type">/</span>G1)<span class="org-type">^</span>(1<span class="org-type">/</span>n)<span class="org-type">*</span>(1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (1<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n)))<span class="org-type">^</span>n;
W2 = 0.22<span class="org-type">*</span>(1 <span class="org-type">+</span> s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>1)<span class="org-type">^</span>2<span class="org-type">/</span>(sqrt(1e<span class="org-type">-</span>4) <span class="org-type">+</span> s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>1)<span class="org-type">^</span>2<span class="org-type">*</span>(1 <span class="org-type">+</span> s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>10)<span class="org-type">^</span>2<span class="org-type">/</span>(1 <span class="org-type">+</span> s<span class="org-type">/</span>2<span class="org-type">/</span><span class="org-constant">pi</span><span class="org-type">/</span>1000)<span class="org-type">^</span>2;
n = 3; w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>10; G0 = 1000; G1 = 0.1; Gc = 1<span class="org-type">/</span>2;
W3 = (((1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n))<span class="org-type">/</span>((1<span class="org-type">/</span>G1)<span class="org-type">^</span>(1<span class="org-type">/</span>n)<span class="org-type">*</span>(1<span class="org-type">/</span>w0)<span class="org-type">*</span>sqrt((1<span class="org-type">-</span>(G0<span class="org-type">/</span>Gc)<span class="org-type">^</span>(2<span class="org-type">/</span>n))<span class="org-type">/</span>(1<span class="org-type">-</span>(Gc<span class="org-type">/</span>G1)<span class="org-type">^</span>(2<span class="org-type">/</span>n)))<span class="org-type">*</span>s <span class="org-type">+</span> (1<span class="org-type">/</span>Gc)<span class="org-type">^</span>(1<span class="org-type">/</span>n)))<span class="org-type">^</span>n;
</pre>
</div>
<div id="org743463d" class="figure">
<p><img src="figs/three_weighting_functions.png" alt="three_weighting_functions.png" />
</p>
<p><span class="figure-number">Figure 6: </span>Three weighting functions used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters</p>
</div>
</div>
</div>
<div id="outline-container-org24d1ff1" class="outline-3">
<h3 id="org24d1ff1"><span class="section-number-3">2.3</span> H-Infinity Synthesis</h3>
<div class="outline-text-3" id="text-2-3">
<p>
Then we create the generalized plant <code>P</code>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">P = [W1 <span class="org-type">-</span>W1 <span class="org-type">-</span>W1;
0 W2 0 ;
0 0 W3;
1 0 0];
</pre>
</div>
<p>
And we do the \(\mathcal{H}_\infty\) synthesis.
</p>
<div class="org-src-container">
<pre class="src src-matlab">[H, <span class="org-type">~</span>, gamma, <span class="org-type">~</span>] = hinfsyn(P, 1, 2,<span class="org-string">'TOLGAM'</span>, 0.001, <span class="org-string">'METHOD'</span>, <span class="org-string">'ric'</span>, <span class="org-string">'DISPLAY'</span>, <span class="org-string">'on'</span>);
</pre>
</div>
<pre class="example" id="org4dc4918">
[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 &lt; gamma &lt;= 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
</pre>
</div>
</div>
<div id="outline-container-orgef81b1e" class="outline-3">
<h3 id="orgef81b1e"><span class="section-number-3">2.4</span> Obtained Complementary Filters</h3>
<div class="outline-text-3" id="text-2-4">
<p>
The obtained filters are:
</p>
<div class="org-src-container">
<pre class="src src-matlab">H2 = tf(H(1));
H3 = tf(H(2));
H1 = 1 <span class="org-type">-</span> H2 <span class="org-type">-</span> H3;
</pre>
</div>
<div id="org7b51da2" class="figure">
<p><img src="figs/three_complementary_filters_results.png" alt="three_complementary_filters_results.png" />
</p>
<p><span class="figure-number">Figure 7: </span>The three complementary filters obtained after \(\mathcal{H}_\infty\) synthesis</p>
</div>
</div>
</div>
</div>
<div id="outline-container-org2c0b99c" class="outline-2">
<h2 id="org2c0b99c"><span class="section-number-2">3</span> Implement complementary filters for LIGO</h2>
<div class="outline-text-2" id="text-3">
<p>
<a id="orgb69dc24"></a>
</p>
<div class="note" id="orgce8d730">
<p>
The Matlab file corresponding to this section is accessible <a href="matlab/comp_filters_ligo.m">here</a>.
</p>
</div>
<p>
Let&rsquo;s try to design complementary filters that are corresponding to the complementary filters design for the LIGO and described in (<a href="#citeproc_bib_item_1">Hua 2005</a>).
</p>
<p>
The FIR complementary filters designed in (<a href="#citeproc_bib_item_1">Hua 2005</a>) are of order 512.
</p>
</div>
<div id="outline-container-org6ea0071" class="outline-3">
<h3 id="org6ea0071"><span class="section-number-3">3.1</span> Specifications</h3>
<div class="outline-text-3" id="text-3-1">
<p>
The specifications for the filters are:
</p>
<ol class="org-ol">
<li>From \(0\) to \(0.008\text{ Hz}\),the magnitude of the filter’s transfer function should be less than or equal to \(8 \times 10^{-3}\)</li>
<li>From \(0.008\text{ Hz}\) to \(0.04\text{ Hz}\), it attenuates the input signal proportional to frequency cubed</li>
<li>Between \(0.04\text{ Hz}\) and \(0.1\text{ Hz}\), the magnitude of the transfer function should be less than 3</li>
<li>Above \(0.1\text{ Hz}\), the maximum of the magnitude of the complement filter should be as close to zero as possible. In our system, we would like to have the magnitude of the complementary filter to be less than \(0.1\). As the filters obtained in (<a href="#citeproc_bib_item_1">Hua 2005</a>) have a magnitude of \(0.045\), we will set that as our requirement</li>
</ol>
<p>
The specifications are translated in upper bounds of the complementary filters are shown on figure <a href="#orgcfe60ed">8</a>.
</p>
<div id="orgcfe60ed" class="figure">
<p><img src="figs/ligo_specifications.png" alt="ligo_specifications.png" />
</p>
<p><span class="figure-number">Figure 8: </span>Specification for the LIGO complementary filters</p>
</div>
</div>
</div>
<div id="outline-container-org8f811fb" class="outline-3">
<h3 id="org8f811fb"><span class="section-number-3">3.2</span> FIR Filter</h3>
<div class="outline-text-3" id="text-3-2">
<p>
We here try to implement the FIR complementary filter synthesis as explained in (<a href="#citeproc_bib_item_1">Hua 2005</a>).
For that, we use the <a href="http://cvxr.com/cvx/">CVX matlab Toolbox</a>.
</p>
<p>
We setup the CVX toolbox and use the <code>SeDuMi</code> solver.
</p>
<div class="org-src-container">
<pre class="src src-matlab">cvx_startup;
cvx_solver sedumi;
</pre>
</div>
<p>
We define the frequency vectors on which we will constrain the norm of the FIR filter.
</p>
<div class="org-src-container">
<pre class="src src-matlab">w1 = 0<span class="org-type">:</span>4.06e<span class="org-type">-</span>4<span class="org-type">:</span>0.008;
w2 = 0.008<span class="org-type">:</span>4.06e<span class="org-type">-</span>4<span class="org-type">:</span>0.04;
w3 = 0.04<span class="org-type">:</span>8.12e<span class="org-type">-</span>4<span class="org-type">:</span>0.1;
w4 = 0.1<span class="org-type">:</span>8.12e<span class="org-type">-</span>4<span class="org-type">:</span>0.83;
</pre>
</div>
<p>
We then define the order of the FIR filter.
</p>
<div class="org-src-container">
<pre class="src src-matlab">n = 512;
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab">A1 = [ones(length(w1),1), cos(kron(w1<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
A2 = [ones(length(w2),1), cos(kron(w2<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
A3 = [ones(length(w3),1), cos(kron(w3<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
A4 = [ones(length(w4),1), cos(kron(w4<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
B1 = [zeros(length(w1),1), sin(kron(w1<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
B2 = [zeros(length(w2),1), sin(kron(w2<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
B3 = [zeros(length(w3),1), sin(kron(w3<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
B4 = [zeros(length(w4),1), sin(kron(w4<span class="org-type">'.*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span>),[1<span class="org-type">:</span>n<span class="org-type">-</span>1]))];
</pre>
</div>
<p>
We run the convex optimization.
</p>
<div class="org-src-container">
<pre class="src src-matlab">cvx_begin
variable y(n<span class="org-type">+</span>1,1)
<span class="org-comment">% t</span>
maximize(<span class="org-type">-</span>y(1))
<span class="org-keyword">for</span> <span class="org-variable-name"><span class="org-constant">i</span></span> = <span class="org-constant">1:length(w1)</span>
norm([0 A1(<span class="org-constant">i</span>,<span class="org-type">:</span>); 0 B1(<span class="org-constant">i</span>,<span class="org-type">:</span>)]<span class="org-type">*</span>y) <span class="org-type">&lt;=</span> 8e<span class="org-type">-</span>3;
<span class="org-keyword">end</span>
<span class="org-keyword">for</span> <span class="org-variable-name"><span class="org-constant">i</span></span> = <span class="org-constant">1:length(w2)</span>
norm([0 A2(<span class="org-constant">i</span>,<span class="org-type">:</span>); 0 B2(<span class="org-constant">i</span>,<span class="org-type">:</span>)]<span class="org-type">*</span>y) <span class="org-type">&lt;=</span> 8e<span class="org-type">-</span>3<span class="org-type">*</span>(2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>w2(<span class="org-constant">i</span>)<span class="org-type">/</span>(0.008<span class="org-type">*</span>2<span class="org-type">*</span><span class="org-constant">pi</span>))<span class="org-type">^</span>3;
<span class="org-keyword">end</span>
<span class="org-keyword">for</span> <span class="org-variable-name"><span class="org-constant">i</span></span> = <span class="org-constant">1:length(w3)</span>
norm([0 A3(<span class="org-constant">i</span>,<span class="org-type">:</span>); 0 B3(<span class="org-constant">i</span>,<span class="org-type">:</span>)]<span class="org-type">*</span>y) <span class="org-type">&lt;=</span> 3;
<span class="org-keyword">end</span>
<span class="org-keyword">for</span> <span class="org-variable-name"><span class="org-constant">i</span></span> = <span class="org-constant">1:length(w4)</span>
norm([[1 0]<span class="org-type">'-</span> [0 A4(<span class="org-constant">i</span>,<span class="org-type">:</span>); 0 B4(<span class="org-constant">i</span>,<span class="org-type">:</span>)]<span class="org-type">*</span>y]) <span class="org-type">&lt;=</span> y(1);
<span class="org-keyword">end</span>
cvx_end
h = y(2<span class="org-type">:</span>end);
</pre>
</div>
<pre class="example" id="orgd562951">
cvx_begin
variable y(n+1,1)
% t
maximize(-y(1))
for i = 1:length(w1)
norm([0 A1(i,:); 0 B1(i,:)]*y) &lt;= 8e-3;
end
for i = 1:length(w2)
norm([0 A2(i,:); 0 B2(i,:)]*y) &lt;= 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) &lt;= 3;
end
for i = 1:length(w4)
norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) &lt;= 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);
</pre>
<p>
Finally, we compute the filter response over the frequency vector defined and the result is shown on figure <a href="#org0d25e01">9</a> which is very close to the filters obtain in (<a href="#citeproc_bib_item_1">Hua 2005</a>).
</p>
<div class="org-src-container">
<pre class="src src-matlab">w = [w1 w2 w3 w4];
H = [exp(<span class="org-type">-</span><span class="org-constant">j</span><span class="org-type">*</span>kron(w<span class="org-type">'.*</span>2<span class="org-type">*</span><span class="org-constant">pi</span>,[0<span class="org-type">:</span>n<span class="org-type">-</span>1]))]<span class="org-type">*</span>h;
</pre>
</div>
<div id="org0d25e01" class="figure">
<p><img src="figs/fir_filter_ligo.png" alt="fir_filter_ligo.png" />
</p>
<p><span class="figure-number">Figure 9: </span>FIR Complementary filters obtain after convex optimization (<a href="./figs/fir_filter_ligo.png">png</a>, <a href="./figs/fir_filter_ligo.pdf">pdf</a>)</p>
</div>
</div>
</div>
<div id="outline-container-org9a34a27" class="outline-3">
<h3 id="org9a34a27"><span class="section-number-3">3.3</span> Weights</h3>
<div class="outline-text-3" id="text-3-3">
<p>
We design weights that will be used for the \(\mathcal{H}_\infty\) synthesis of the complementary filters.
These weights will determine the order of the obtained filters.
Here are the requirements on the filters:
</p>
<ul class="org-ul">
<li>reasonable order</li>
<li>to be as close as possible to the specified upper bounds</li>
<li>stable minimum phase</li>
</ul>
<p>
The bode plot of the weights is shown on figure <a href="#orge27ef5e">10</a>.
</p>
<div id="orge27ef5e" class="figure">
<p><img src="figs/ligo_weights.png" alt="ligo_weights.png" />
</p>
<p><span class="figure-number">Figure 10: </span>Weights for the \(\mathcal{H}_\infty\) synthesis</p>
</div>
</div>
</div>
<div id="outline-container-orgc9f5d68" class="outline-3">
<h3 id="orgc9f5d68"><span class="section-number-3">3.4</span> H-Infinity Synthesis</h3>
<div class="outline-text-3" id="text-3-4">
<p>
We define the generalized plant as shown on figure <a href="#org46317e2">1</a>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">P = [0 wL;
wH <span class="org-type">-</span>wH;
1 0];
</pre>
</div>
<p>
And we do the \(\mathcal{H}_\infty\) synthesis using the <code>hinfsyn</code> command.
</p>
<div class="org-src-container">
<pre class="src src-matlab">[Hl, <span class="org-type">~</span>, gamma, <span class="org-type">~</span>] = hinfsyn(P, 1, 1,<span class="org-string">'TOLGAM'</span>, 0.001, <span class="org-string">'METHOD'</span>, <span class="org-string">'ric'</span>, <span class="org-string">'DISPLAY'</span>, <span class="org-string">'on'</span>);
</pre>
</div>
<pre class="example" id="org51b2320">
[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 &lt; gamma &lt;= 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
</pre>
<p>
The high pass filter is defined as \(H_H = 1 - H_L\).
</p>
<div class="org-src-container">
<pre class="src src-matlab">Hh = 1 <span class="org-type">-</span> Hl;
</pre>
</div>
<p>
The size of the filters is shown below.
</p>
<pre class="example" id="orgede384d">
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.
</pre>
<p>
The bode plot of the obtained filters as shown on figure <a href="#org941eede">11</a>.
</p>
<div id="org941eede" class="figure">
<p><img src="figs/hinf_synthesis_ligo_results.png" alt="hinf_synthesis_ligo_results.png" />
</p>
<p><span class="figure-number">Figure 11: </span>Obtained complementary filters using the \(\mathcal{H}_\infty\) synthesis</p>
</div>
</div>
</div>
<div id="outline-container-orgd83a581" class="outline-3">
<h3 id="orgd83a581"><span class="section-number-3">3.5</span> Compare FIR and H-Infinity Filters</h3>
<div class="outline-text-3" id="text-3-5">
<p>
Let&rsquo;s now compare the FIR filters designed in (<a href="#citeproc_bib_item_1">Hua 2005</a>) and the one obtained with the \(\mathcal{H}_\infty\) synthesis on figure <a href="#org4f70e8a">12</a>.
</p>
<div id="org4f70e8a" class="figure">
<p><img src="figs/comp_fir_ligo_hinf.png" alt="comp_fir_ligo_hinf.png" />
</p>
<p><span class="figure-number">Figure 12: </span>Comparison between the FIR filters developped for LIGO and the \(\mathcal{H}_\infty\) complementary filters</p>
</div>
</div>
</div>
</div>
<p>
</p>
<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><h2 class='citeproc-org-bib-h2'>Bibliography</h2>
<div class="csl-bib-body">
<div class="csl-entry"><a name="citeproc_bib_item_1"></a>Hua, Wensheng. 2005. “Low Frequency Vibration Isolation and Alignment System for Advanced LIGO.” stanford university.</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Thomas Dehaeze</p>
<p class="date">Created: 2020-10-26 lun. 18:59</p>
</div>
</body>
</html>

931
matlab/index.org Normal file
View File

@ -0,0 +1,931 @@
#+TITLE: Complementary Filters Shaping Using $\mathcal{H}_\infty$ Synthesis - Matlab Computation
:DRAWER:
#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP: ../index.html
#+LATEX_CLASS: cleanreport
#+LATEX_CLASS_OPTIONS: [tocnp, secbreak, minted]
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="../css/readtheorg.css"/>
#+HTML_HEAD: <script src="../js/jquery.min.js"></script>
#+HTML_HEAD: <script src="../js/bootstrap.min.js"></script>
#+HTML_HEAD: <script src="../js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script src="../js/readtheorg.js"></script>
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :tangle matlab/comp_filters_design.m
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:
* Introduction :ignore:
In this document, the design of complementary filters is studied.
One use of complementary filter is described below:
#+begin_quote
The basic idea of a complementary filter involves taking two or more sensors, filtering out unreliable frequencies for each sensor, and combining the filtered outputs to get a better estimate throughout the entire bandwidth of the system.
To achieve this, the sensors included in the filter should complement one another by performing better over specific parts of the system bandwidth.
#+end_quote
This document is divided into several sections:
- in section [[sec:h_inf_synthesis_complementary_filters]], the $\mathcal{H}_\infty$ synthesis is used for generating two complementary filters
- in section [[sec:three_comp_filters]], a method using the $\mathcal{H}_\infty$ synthesis is proposed to shape three of more complementary filters
- in section [[sec:comp_filters_ligo]], the $\mathcal{H}_\infty$ synthesis is used and compared with FIR complementary filters used for LIGO
#+begin_note
Add the Matlab code use to obtain the results presented in the paper are accessible [[file:matlab.zip][here]] and presented below.
#+end_note
* H-Infinity synthesis of complementary filters
:PROPERTIES:
:header-args:matlab+: :tangle matlab/h_inf_synthesis_complementary_filters.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:h_inf_synthesis_complementary_filters>>
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/h_inf_synthesis_complementary_filters.m][here]].
#+end_note
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab
freqs = logspace(-1, 3, 1000);
#+end_src
** Synthesis Architecture
We here synthesize two complementary filters using the $\mathcal{H}_\infty$ synthesis.
The goal is to specify upper bounds on the norms of the two complementary filters $H_1(s)$ and $H_2(s)$ while ensuring their complementary property ($H_1(s) + H_2(s) = 1$).
In order to do so, we use the generalized plant shown on figure [[fig:h_infinity_robst_fusion]] where $W_1(s)$ and $W_2(s)$ are weighting transfer functions that will be used to shape $H_1(s)$ and $H_2(s)$ respectively.
#+name: fig:h_infinity_robst_fusion
#+caption: $\mathcal{H}_\infty$ synthesis of the complementary filters
[[file:figs-tikz/h_infinity_robust_fusion.png]]
The $\mathcal{H}_\infty$ synthesis applied on this generalized plant will give a transfer function $H_2$ (figure [[fig:h_infinity_robst_fusion]]) such that the $\mathcal{H}_\infty$ norm of the transfer function from $w$ to $[z_1,\ z_2]$ is less than one:
\[ \left\| \begin{array}{c} (1 - H_2(s)) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \]
Thus, if the above condition is verified, we can define $H_1(s) = 1 - H_2(s)$ and we have that:
\[ \left\| \begin{array}{c} H_1(s) W_1(s) \\ H_2(s) W_2(s) \end{array} \right\|_\infty < 1 \]
Which is almost (with an maximum error of $\sqrt{2}$) equivalent to:
\begin{align*}
|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{align*}
We then see that $W_1(s)$ and $W_2(s)$ can be used to shape both $H_1(s)$ and $H_2(s)$ while ensuring their complementary property by the definition of $H_1(s) = 1 - H_2(s)$.
** Design of Weighting Function
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
\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 [[fig:weight_formula]].
#+name: fig:weight_formula
#+caption: Amplitude of the proposed formula for the weighting functions
[[file:figs-tikz/weight_formula.png]]
#+begin_src matlab
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-4, 20]);
xticks([0.1, 1, 10, 100, 1000]);
legend('location', 'northeast', 'FontSize', 8);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/weights_W1_W2.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:weights_W1_W2
#+caption: Weights on the complementary filters $W_1$ and $W_2$ and the associated performance weights
#+RESULTS:
[[file:figs/weights_W1_W2.png]]
** H-Infinity Synthesis
We define the generalized plant $P$ on matlab.
#+begin_src matlab
P = [W1 -W1;
0 W2;
1 0];
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
#+end_src
#+RESULTS:
#+begin_example
[H2, ~, 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.1000 < gamma <= 1050.0000
gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
1.050e+03 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
525.050 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
262.575 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
131.337 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
65.719 2.8e+01 2.4e-07 4.1e+00 -9.5e-14 0.0000 p
32.909 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
16.505 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
8.302 2.8e+01 2.4e-07 4.1e+00 -7.2e-14 0.0000 p
4.201 2.8e+01 2.4e-07 4.1e+00 -2.5e-25 0.0000 p
2.151 2.7e+01 2.4e-07 4.1e+00 -3.8e-14 0.0000 p
1.125 2.6e+01 2.4e-07 4.1e+00 -5.4e-24 0.0000 p
0.613 2.3e+01 -3.7e+01# 4.1e+00 0.0e+00 0.0000 f
0.869 2.6e+01 -3.7e+02# 4.1e+00 0.0e+00 0.0000 f
0.997 2.6e+01 -1.1e+04# 4.1e+00 0.0e+00 0.0000 f
1.061 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.029 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.013 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.005 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
1.001 2.6e+01 -3.1e+04# 4.1e+00 -3.8e-14 0.0000 f
1.003 2.6e+01 -2.8e+05# 4.1e+00 0.0e+00 0.0000 f
1.004 2.6e+01 2.4e-07 4.1e+00 -5.8e-24 0.0000 p
1.004 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
Gamma value achieved: 1.0036
#+end_example
We then define the high pass filter $H_1 = 1 - H_2$. The bode plot of both $H_1$ and $H_2$ is shown on figure [[fig:hinf_filters_results]].
#+begin_src matlab
H1 = 1 - H2;
#+end_src
** Obtained Complementary Filters
The obtained complementary filters are shown on figure [[fig:hinf_filters_results]].
#+begin_src matlab :exports none
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$');
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([1e-4, 20]);
yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',2)
plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/hinf_filters_results.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:hinf_filters_results
#+caption: Obtained complementary filters using $\mathcal{H}_\infty$ synthesis
#+RESULTS:
[[file:figs/hinf_filters_results.png]]
* Generating 3 complementary filters
:PROPERTIES:
:header-args:matlab+: :tangle matlab/three_comp_filters.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:three_comp_filters>>
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/three_comp_filters.m][here]].
#+end_note
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab
freqs = logspace(-2, 4, 1000);
#+end_src
** Theory
We want:
\begin{align*}
& |H_1(j\omega)| < 1/|W_1(j\omega)|, \quad \forall\omega\\
& |H_2(j\omega)| < 1/|W_2(j\omega)|, \quad \forall\omega\\
& |H_3(j\omega)| < 1/|W_3(j\omega)|, \quad \forall\omega\\
& H_1(s) + H_2(s) + H_3(s) = 1
\end{align*}
For that, we use the $\mathcal{H}_\infty$ synthesis with the architecture shown on figure [[fig:comp_filter_three_hinf]].
#+name: fig:comp_filter_three_hinf
#+caption: Generalized architecture for generating 3 complementary filters
[[file:figs-tikz/comp_filter_three_hinf.png]]
The $\mathcal{H}_\infty$ objective is:
\begin{align*}
& |(1 - H_2(j\omega) - H_3(j\omega)) W_1(j\omega)| < 1, \quad \forall\omega\\
& |H_2(j\omega) W_2(j\omega)| < 1, \quad \forall\omega\\
& |H_3(j\omega) W_3(j\omega)| < 1, \quad \forall\omega\\
\end{align*}
And thus if we choose $H_1 = 1 - H_2 - H_3$ we have solved the problem.
** Weights
First we define the weights.
#+begin_src matlab
n = 2; w0 = 2*pi*1; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W3 = (((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;
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
legend('location', 'northeast', 'FontSize', 8);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/three_weighting_functions.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:three_weighting_functions
#+caption: Three weighting functions used for the $\mathcal{H}_\infty$ synthesis of the complementary filters
#+RESULTS:
[[file:figs/three_weighting_functions.png]]
** H-Infinity Synthesis
Then we create the generalized plant =P=.
#+begin_src matlab
P = [W1 -W1 -W1;
0 W2 0 ;
0 0 W3;
1 0 0];
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis.
#+begin_src matlab :results output replace :exports both
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
#+end_src
#+RESULTS:
#+begin_example
[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
#+end_example
** Obtained Complementary Filters
The obtained filters are:
#+begin_src matlab
H2 = tf(H(1));
H3 = tf(H(2));
H1 = 1 - H2 - H3;
#+end_src
#+begin_src matlab :exports none
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca,'ColorOrderIndex',3)
plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([1e-4, 20]);
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))));
set(gca,'ColorOrderIndex',2)
plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))));
set(gca,'ColorOrderIndex',3)
plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz'))));
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:90:360]); ylim([-270, 270]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/three_complementary_filters_results.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:three_complementary_filters_results
#+caption: The three complementary filters obtained after $\mathcal{H}_\infty$ synthesis
#+RESULTS:
[[file:figs/three_complementary_filters_results.png]]
* Implement complementary filters for LIGO
:PROPERTIES:
:header-args:matlab+: :tangle matlab/comp_filters_ligo.m
:header-args:matlab+: :comments org :mkdirp yes
:END:
<<sec:comp_filters_ligo>>
** Introduction :ignore:
#+begin_note
The Matlab file corresponding to this section is accessible [[file:matlab/comp_filters_ligo.m][here]].
#+end_note
Let's try to design complementary filters that are corresponding to the complementary filters design for the LIGO and described in cite:hua05_low_ligo.
The FIR complementary filters designed in cite:hua05_low_ligo are of order 512.
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
<<matlab-dir>>
#+end_src
#+begin_src matlab :exports none :results silent :noweb yes
<<matlab-init>>
#+end_src
#+begin_src matlab
freqs = logspace(-3, 0, 1000);
#+end_src
** Specifications
The specifications for the filters are:
1. From $0$ to $0.008\text{ Hz}$,the magnitude of the filter’s transfer function should be less than or equal to $8 \times 10^{-3}$
2. From $0.008\text{ Hz}$ to $0.04\text{ Hz}$, it attenuates the input signal proportional to frequency cubed
3. Between $0.04\text{ Hz}$ and $0.1\text{ Hz}$, the magnitude of the transfer function should be less than 3
4. Above $0.1\text{ Hz}$, the maximum of the magnitude of the complement filter should be as close to zero as possible. In our system, we would like to have the magnitude of the complementary filter to be less than $0.1$. As the filters obtained in cite:hua05_low_ligo have a magnitude of $0.045$, we will set that as our requirement
The specifications are translated in upper bounds of the complementary filters are shown on figure [[fig:ligo_specifications]].
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1)
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1)
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2)
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-4, 10]);
legend('location', 'southeast', 'FontSize', 8);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/ligo_specifications.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:ligo_specifications
#+caption: Specification for the LIGO complementary filters
#+RESULTS:
[[file:figs/ligo_specifications.png]]
** TODO FIR Filter
We here try to implement the FIR complementary filter synthesis as explained in cite:hua05_low_ligo.
For that, we use the [[http://cvxr.com/cvx/][CVX matlab Toolbox]].
We setup the CVX toolbox and use the =SeDuMi= solver.
#+begin_src matlab
cvx_startup;
cvx_solver sedumi;
#+end_src
We define the frequency vectors on which we will constrain the norm of the FIR filter.
#+begin_src matlab
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;
#+end_src
We then define the order of the FIR filter.
#+begin_src matlab
n = 512;
#+end_src
#+begin_src matlab
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]))];
#+end_src
We run the convex optimization.
#+begin_src matlab :results output replace :wrap example
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);
#+end_src
#+RESULTS:
#+begin_example
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);
#+end_example
Finally, we compute the filter response over the frequency vector defined and the result is shown on figure [[fig:fir_filter_ligo]] which is very close to the filters obtain in cite:hua05_low_ligo.
#+begin_src matlab
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
#+end_src
#+begin_src matlab :exports none
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
plot(w, abs(H), 'k-');
plot(w, abs(1-H), 'k--');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 5]);
% Phase
ax2 = nexttile;
hold on;
plot(w, 180/pi*angle(H), 'k-');
plot(w, 180/pi*angle(1-H), 'k--');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-180:90:180]); ylim([-180, 180]);
linkaxes([ax1,ax2],'x');
xlim([1e-3, 1]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:fir_filter_ligo
#+caption: FIR Complementary filters obtain after convex optimization
#+RESULTS:
[[file:figs/fir_filter_ligo.png]]
** Weights
We design weights that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters.
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 minimum phase
The bode plot of the weights is shown on figure [[fig:ligo_weights]].
#+begin_src matlab :exports none
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));
#+end_src
#+begin_src matlab :exports none
n = 20; Rp = 1; Wp = 2*pi*0.102;
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
wL = minreal(ss(wL));
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|w_H|^{-1}$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|w_L|^{-1}$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec.');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
plot([0.1, 10], [0.045, 0.045], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast', 'FontSize', 8);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/ligo_weights.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:ligo_weights
#+caption: Weights for the $\mathcal{H}_\infty$ synthesis
#+RESULTS:
[[file:figs/ligo_weights.png]]
** H-Infinity Synthesis
We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion]].
#+begin_src matlab
P = [0 wL;
wH -wH;
1 0];
#+end_src
And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
#+begin_src matlab :results output replace :exports both :wrap example
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
#+end_src
#+RESULTS:
#+begin_example
[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
#+end_example
The high pass filter is defined as $H_H = 1 - H_L$.
#+begin_src matlab
Hh = 1 - Hl;
#+end_src
#+begin_src matlab :exports none
Hh = minreal(Hh);
Hl = minreal(Hl);
#+end_src
The size of the filters is shown below.
#+begin_src matlab :exports results :results output replace :wrap example
size(Hh), size(Hl)
#+end_src
#+RESULTS:
#+begin_example
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.
#+end_example
The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_ligo_results]].
#+begin_src matlab :exports none
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1);
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1);
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2);
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast', 'FontSize', 8);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/hinf_synthesis_ligo_results.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:hinf_synthesis_ligo_results
#+caption: Obtained complementary filters using the $\mathcal{H}_\infty$ synthesis
#+RESULTS:
[[file:figs/hinf_synthesis_ligo_results.png]]
** TODO Compare FIR and H-Infinity Filters
Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one obtained with the $\mathcal{H}_\infty$ synthesis on figure [[fig:comp_fir_ligo_hinf]].
#+begin_src matlab :exports none
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_H(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',1);
plot(w, abs(H), '--', ...
'DisplayName', '$H_H(s)$ - FIR');
set(gca,'ColorOrderIndex',2);
plot(w, abs(1-H), '--', ...
'DisplayName', '$H_L(s)$ - FIR');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 10]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hh, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hl, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*angle(H), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*angle(1-H), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-180:90:180]); ylim([-180, 180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/comp_fir_ligo_hinf.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:comp_fir_ligo_hinf
#+caption: Comparison between the FIR filters developped for LIGO and the $\mathcal{H}_\infty$ complementary filters
#+RESULTS:
[[file:figs/comp_fir_ligo_hinf.png]]
* Bibliography :ignore:
bibliographystyle:unsrt
bibliography:ref.bib

1074
matlab/mat/comp_ligo_fir.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/comp_ligo_hinf.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/hinf_filters_results.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/hinf_three_results.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/hinf_three_weights.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/hinf_weights.csv Normal file
View File

File diff suppressed because it is too large Load Diff

1001
matlab/mat/ligo_weights.csv Normal file
View File

File diff suppressed because it is too large Load Diff

501
matlab/mat/weight_formula.csv Normal file
View File

@ -0,0 +1,501 @@
freqs,ampl
0.1,0.00124984376952881
0.101393945767529,0.00125685750941772
0.102807322383087,0.00126406813733846
0.104240400702156,0.00127148117998881
0.105693455355799,0.00127910231918553
0.107166764803286,0.00128693739621714
0.10866061138546,0.00129499241631883
0.110175281378839,0.00130327355327284
0.111711065050482,0.00131178715413789
0.113268256713615,0.00132053974411124
0.114847154784029,0.00132953803152714
0.116448061837269,0.00133878891299542
0.118071284666619,0.00134829947868417
0.119717134341897,0.00135807701775069
0.121385926269063,0.00136812902392454
0.123077980250667,0.00137846320124731
0.124793620547131,0.00138908746997326
0.126533175938894,0.00140000997263542
0.128296979789415,0.00141123908028183
0.130085370109057,0.0014227833988865
0.131898689619867,0.00143465177594028
0.13373728582125,0.00144685330722638
0.135601511056563,0.00145939734378589
0.137491722580642,0.00147229349907857
0.139408282628258,0.0014855516563444
0.14135155848354,0.00149918197617139
0.143321922550357,0.00151319490427561
0.14531975242369,0.00152760117949929
0.147345430961984,0.001542411842033
0.149399346360526,0.00155763824186835
0.151481892225835,0.00157329204748753
0.15359346765109,0.00158938525479641
0.155734477292613,0.00160593019630794
0.157905331447418,0.00162293955058283
0.160106446131832,0.00164042635193486
0.162338243161228,0.00165840400040801
0.164601150230855,0.00167688627203319
0.166895600997802,0.00169588732937223
0.169222035164104,0.00171542173235734
0.171580898561001,0.0017355044494341
0.17397264323438,0.00175615086901664
0.176397727531404,0.00177737681126359
0.178856616188346,0.00179919854018389
0.18134978041965,0.00182163277608158
0.183877698008233,0.00184469670834908
0.186440853397049,0.00186840800861865
0.189039737781922,0.00189278484428217
0.191674849205682,0.00191784589238927
0.194346692653602,0.00194361035393467
0.19705578015018,0.00197009796854536
0.199802630857255,0.00199732902957889
0.202587771173502,0.0020253243996442
0.205411734835306,0.00205410552655685
0.208275063019052,0.00208369445974066
0.211178304444824,0.00211411386708823
0.214122015481573,0.00214538705229329
0.217106760253727,0.00217753797266776
0.220133110749303,0.0022105912574573
0.223201646929523,0.00224457222666906
0.226312956839953,0.00227950691042593
0.229467636723194,0.00231542206886212
0.232666291133146,0.00235234521257481
0.235909533050864,0.00239030462364776
0.239197984002024,0.00242932937726244
0.242532274176035,0.00246944936391327
0.245913042546805,0.00251069531224364
0.249340936995188,0.00255309881252008
0.25281661443315,0.00259669234076224
0.256340740929651,0.00264150928354692
0.259913991838293,0.00268758396350498
0.263537051926739,0.00273495166553027
0.267210615507943,0.00278364866372035
0.270935386573205,0.00283371224906948
0.274712078927081,0.00288518075793453
0.278541416324177,0.00293809360129545
0.282424132607843,0.00299249129483226
0.286360971850812,0.0030484154898412
0.290352688497781,0.0031059090050133
0.294400047510004,0.00316501585909936
0.298503824511873,0.0032257813044857
0.302664805939569,0.00328825186170615
0.306883789191764,0.00335247535491597
0.311161582782436,0.00341850094835445
0.315499006495809,0.00348637918382354
0.319896891543454,0.00355616201921048
0.324356080723581,0.00362790286808339
0.328877428582551,0.00370165664038944
0.333461801578636,0.00377747978428603
0.338110078248068,0.00385543032913618
0.342823149373397,0.0039355679297004
0.347601918154198,0.00401795391155783
0.352447300380159,0.00410265131779068
0.357360224606579,0.00418972495696669
0.362341632332315,0.00427924145245532
0.367392478180207,0.00437126929311439
0.372513730080021,0.0044658788853849
0.377706369453937,0.00456314260683253
0.382971391404628,0.00466313486117581
0.388309804905961,0.00476593213484159
0.393722632996348,0.0048716130550896
0.399210912974805,0.00498025844974942
0.404775696599732,0.00509195140861367
0.410418050290471,0.00520677734653297
0.416139055331671,0.0053248240682591
0.421939808080503,0.00544618183508416
0.427821420176762,0.0055709434333248
0.433785018755899,0.00569920424470177
0.439831746665023,0.0058310623186666
0.445962762681909,0.00596661844672819
0.45217924173707,0.00610597623883409
0.45848237513891,0.00624924220186199
0.464873370802026,0.00639652582027892
0.471353453478691,0.00654793963902696
0.47792386499356,0.00670359934869576
0.48458586448165,0.00686362387304391
0.491340728629636,0.00702813545893254
0.498189751920516,0.00719725976873649
0.505134246881677,0.00737112597529983
0.512175544336424,0.00754986685950441
0.519314993659021,0.0077336189105217
0.526553963033276,0.00792252242882021
0.533893839714735,0.00811672163200241
0.541336030296538,0.0083163647635472
0.548881960978967,0.00852160420453555
0.556533077842765,0.00873259658843936
0.564290847126254,0.00894950291905528
0.572156755506325,0.00917248869166733
0.580132310383338,0.00940172401752446
0.588219040169996,0.00963738375172108
0.596418494584246,0.00987964762457099
0.604732244946265,0.0101287003765674
0.613161884479578,0.0103847318970236
0.621709028616383,0.0106479373664921
0.630375315307128,0.0109185174030614
0.6391624053344,0.011196678212632
0.648071982631197,0.011482631743278
0.657105754603627,0.0117765958437987
0.666265452458115,0.0120787944265724
0.675552831533164,0.0123894576348226
0.684969671635745,0.0127088220144133
0.69451777738237,0.0130371306902897
0.704198978544929,0.0133746335476852
0.71401513040134,0.0137215874182187
0.723968114091088,0.0140782562710072
0.734059836975721,0.0144449114089228
0.744292233004376,0.0148218316701279
0.754667263084389,0.0152093036350209
0.765186915457082,0.0156076218387322
0.775853206078784,0.0160170889893121
0.786668179007158,0.0164380161917532
0.797633906792928,0.0168707231779958
0.808752490877045,0.017315538543068
0.820026061993413,0.0177727999875122
0.831456780577207,0.0182428545662572
0.843046837178897,0.0187260589440964
0.854798452884041,0.0192227796579347
0.866713879738923,0.0197333933859734
0.878795401182132,0.0202582872240027
0.891045332482152,0.0207978589689765
0.903466021181053,0.0213525174100481
0.916059847544371,0.0219226826272462
0.92882922501725,0.0225087862979786
0.941776600686952,0.0231112720115502
0.954904455751808,0.0237305955918877
0.968215305996708,0.0243672254286671
0.981711702275219,0.0250216428170423
0.995396230998424,0.0256943423061784
1.00927151463057,0.026385832056794
1.02334021219164,0.0270966342079253
1.03760501976691,0.0278272852531217
1.05206867102362,0.0285783364262923
1.06673393773486,0.0293503540974215
1.08160363031071,0.0301439201783789
1.09668059833687,0.0309596325390481
1.1119677311207,0.0317981054340067
1.12746795824495,0.0326599699399889
1.14318425012915,0.0335458744043682
1.15911961859889,0.0344564849048982
1.17527711746295,0.0353924857209554
1.19165984309856,0.036354579816527
1.20827093504478,0.037343489335192
1.22511357660412,0.0383599561073468
1.24219099545262,0.0394047421699245
1.25950646425836,0.0404786302988653
1.27706330130864,0.0415824245545925
1.2948648711459,0.0427169508407526
1.31291458521247,0.043883057476478
1.33121590250431,0.0450816157824338
1.34977233023394,0.0463135206809095
1.36858742450252,0.0475796913102165
1.38766479098131,0.0488810716536558
1.40700808560268,0.0502186311833147
1.42662101526074,0.0515933655189563
1.44650733852165,0.0530062971022605
1.46667086634397,0.0544584758866766
1.48711546280895,0.0559509800431428
1.50784504586105,0.0574849166819288
1.52886358805873,0.0590614225908512
1.55017511733577,0.0606816649901108
1.57178371777316,0.0623468423039935
1.59369353038178,0.0640581849496757
1.61590875389592,0.0658169561433637
1.63843364557798,0.0676244527239964
1.66127252203429,0.0694820059947267
1.68442976004232,0.071390982582396
1.70790979738943,0.0733527853151986
1.73171713372335,0.0753688541187309
1.75585633141447,0.0774406669306033
1.78033201643011,0.0795697406337812
1.80514887922111,0.0817576320088108
1.83031167562061,0.0840059387050673
1.85582522775552,0.086316300231147
1.88169442497056,0.0886903989645098
1.90792422476527,0.0911299611804565
1.93451965374405,0.0936367581005048
1.96148580857943,0.0962126069602082
1.98882785698881,0.0988593720964329
2.01655103872475,0.101578966054087
2.04466066657912,0.104373350712264
2.07316212740123,0.107244538429732
2.10206088313016,0.110194593209675
2.13136247184144,0.113225631883545
2.16107250880838,0.116339825313863
2.19119668757815,0.119539399615746
2.22174078106288,0.12282663739692
2.25271064264598,0.12620387901591
2.28411220730382,0.129673523858072
2.31595149274315,0.133238031629055
2.34823460055428,0.136899923665257
2.38096771738036,0.140661784260755
2.41415711610302,0.144526262010144
2.44780915704444,0.148496071166648
2.48193028918626,0.152573993014796
2.51652705140539,0.156762877256901
2.55160607372718,0.161065643412464
2.58717407859592,0.1654852822296
2.62323788216312,0.170024857107441
2.65980439559376,0.174687505528422
2.69688062639069,0.179476440499244
2.73447367973758,0.184394951999205
2.77259075986048,0.189446408434499
2.81123917140846,0.194634258096962
2.85042632085343,0.199962030625633
2.89015971790951,0.205433338469367
2.93044697697214,0.211051878348636
2.97129581857733,0.216821432714486
3.01271407088116,0.222745871202514
3.05470967115997,0.228829152079549
3.09729066733141,0.235075323680605
3.14046521949675,0.241488525833486
3.18424160150465,0.24807299126828
3.22862820253673,0.254833047008791
3.27363352871525,0.261773115742795
3.31926620473319,0.268897717167807
3.36553497550709,0.276211469308868
3.41244870785289,0.283719089804643
3.46001639218511,0.291425397157943
3.50824714423979,0.299335311946534
3.55715020682139,0.307453857989918
3.60673495157403,0.315786163467503
3.65701088077749,0.324337461983377
3.70798762916817,0.333113093572643
3.75967496578547,0.342118505644029
3.81208279584389,0.351359253853236
3.86522116263126,0.36084100290123
3.9191002494334,0.370569527251405
3.97373038148561,0.380550711759293
4.02912202795135,0.390790552208207
4.08528580392856,0.401295155743936
4.1422324724839,0.412070741201301
4.19997294671531,0.423123639315149
4.25851829184341,0.434460292807991
4.31787972733202,0.446087256346283
4.37806862903817,0.45801119635701
4.43909653139217,0.470238890695935
4.50097512960805,0.482777228158622
4.56371628192476,0.495633207825031
4.62733201187869,0.508813938228191
4.6918345106078,0.522326636337199
4.75723613918789,0.53617862634452
4.82354943100147,0.550377338247268
4.89078709413959,0.564930306211964
4.95896201383721,0.579845166711924
5.02808725494248,0.595129656426328
5.09817606442043,0.610791609889683
5.1692418738916,0.626838956880323
5.24129830220606,0.643279719536317
5.31435915805324,0.660122009187079
5.38843844260822,0.677374022888832
5.46355035221488,0.695044039651974
5.5397092811064,0.713140416348375
5.61692982416381,0.731671583286571
5.69522677971282,0.750646039442876
5.77461515235982,0.770072347336476
5.85511015586724,0.789959127536677
5.93672721606913,0.810315052790655
6.01948197382727,0.831148841760246
6.10339028802862,0.852469252356602
6.18846823862439,0.874285074661867
6.27473212971158,0.896605123427405
6.36219849265749,0.919438230138643
6.45088408926769,0.94279323463703
6.54080591499826,0.966678976290355
6.63198120221267,0.991104284703253
6.72442742348425,1.01607796996062
6.81816229494448,1.04160881239745
6.91320377967813,1.06770555188962
7.00957009116562,1.09437687666126
7.10727969677342,1.12163141160536
7.20635132129305,1.14947770611576
7.3068039505295,1.17792422142984
7.40865683493956,1.20697931748277
7.51192949332097,1.23665123927596
7.61664171655289,1.26694810276361
7.72281357138865,1.29787788026358
7.83046540430119,1.3294483854004
7.93961784538228,1.36166725759032
8.05029181229598,1.3945419460805
8.16250851428723,1.42807969355667
8.27628945624635,1.46228751933592
8.39165644283016,1.49717220216369
8.50863158264058,1.53274026263664
8.62723729246145,1.5689979452756
8.74749630155442,1.60595120027567
8.86943165601471,1.643605664963
8.99306672318762,1.68196664499107
9.11842519614657,1.72103909531164
9.24553109823358,1.76082760095885
9.374408787663,1.80133635768771
9.50508296218951,1.84256915251122
9.63757866384109,1.88452934418316
9.771921283718,1.92721984367676
9.90813656685867,1.97064309471195
10.0462506171734,2.01480105438701
10.1862899024469,2.05969517397301
10.3282812594103,2.10532637993194
10.4722518988843,2.15169505522199
10.6182294109938,2.19880102095579
10.7662417704549,2.24664351847938
10.9163173419361,2.29522119194176
11.0684848854941,2.34453207142641
11.2227735620851,2.39457355671776
11.3792129391532,2.44534240177657
11.5378329962966,2.49683469999911
11.6986641310131,2.54904587033564
11.8617371645248,2.60197064434355
12.0270833476851,2.6556030542507
12.1947343669674,2.70993642210359
12.3647223505371,2.76496335007402
12.5370798744092,2.82067571199669
12.7118399686903,2.87706464620789
12.8890361239089,2.93412054975349
13.0687022974335,2.9918330740313
13.2508729199795,3.05019112192989
13.4355829022083,3.10918284652194
13.6228676414165,3.16879565136619
13.8127630283201,3.22901619246727
14.0053054539322,3.28983038193766
14.2005318165368,3.35122339340032
14.3984795287601,3.41317966916477
14.5991865247398,3.47568292920294
14.8026912673951,3.53871618194413
15.0090327557974,3.60226173690188
15.2182505326439,3.66630121913765
15.4303846918356,3.73081558555893
15.645475886161,3.79578514304128
15.8635653350859,3.86118956835606
16.0846948326536,3.92700792987717
16.3089067554933,3.99321871103223
16.5362440709418,4.059799835455
16.766750345277,4.12672869378805
17.0004697520672,4.19398217207637
17.2374470806362,4.26153668168462
17.4777277446468,4.32936819066333
17.7213577908036,4.3974522564814
17.9683839076772,4.4657640600358
18.2188534346517,4.53427844084184
18.4728143709963,4.60296993330202
18.7303153850644,4.67181280394469
18.9914058236194,4.74078108951949
19.256135721292,4.80984863583088
19.5245558101686,4.87898913718794
19.7967175295133,4.94817617634452
20.0726730356257,5.01738326480138
20.3524752118358,5.08658388334002
20.6361776786386,5.15575152265655
20.9238348039698,5.22485972396316
21.2155017136245,5.29388211942536
21.5112343018217,5.36279247230368
21.8110892419152,5.43156471667044
22.1151239972549,5.50017299657441
22.4233968321984,5.56859170452924
22.7359668232772,5.63679551920541
23.0528938705171,5.70475944220955
23.3742387089181,5.77245883384028
23.7000629200933,5.83986944771488
24.0304289440697,5.90696746416749
24.3654000912547,5.97372952232568
24.7050405545683,6.04013275077966
25.049415421745,6.10615479676495
25.3985906878073,6.17177385378773
25.752633267712,6.23696868762925
26.1116110091746,6.30171866067413
26.4755927056707,6.36600375451527
26.8446481096197,6.42980459079649
27.2188479457518,6.49310245026205
27.5982639246618,6.55587928999065
27.9829687565512,6.61811775879929
28.3730361651621,6.67980121081061
28.7685409019059,6.74091371718493
29.169558760188,6.80144007602565
29.5761665899325,6.8613658204741
29.9884423123103,6.92067722501673
30.4064649346707,6.97936131003399
30.8303145656828,7.03740584462663
31.260072430687,7.09479934776086
31.6958208872612,7.15153108777898
32.1376434410028,7.20759108032728
32.5856247615323,7.26297008475719
33.0398506987186,7.31765959905996
33.5004082991313,7.37165185339816
33.9673858227221,7.42493980230096
34.4408727597382,7.47751711559223
34.9209598478727,7.52937816812287
35.4077390896527,7.5805180283804
35.9013037700707,7.63093244605004
36.4017484744614,7.68061783860244
36.9091691066278,7.72957127698342
37.4236629072198,7.77779047048138
37.9453284723694,7.8252737508475
38.4742657725851,7.87202005574321
39.0105761719099,7.91802891158854
39.554362447347,7.96330041588347
40.105728808555,8.00783521907301
40.6647809178186,8.05163450602489
41.2316259102975,8.09469997718672
41.8063724145576,8.13703382948712
42.3891305733878,8.17863873704334
42.980012064908,8.21951783173482
43.5791301239703,8.25967468369993
44.1865995638594,8.29911328181031
44.8025367982949,8.33783801417395
45.4270598637405,8.37585364871619
46.060288442024,8.41316531388383
46.7023438832734,8.44977847951546
47.3533492291712,8.48569893791774
48.0134292365345,8.52093278518449
48.6827104012228,8.55548640279275
49.3613209823792,8.58936643950677
50.0493910270095,8.62257979361857
50.7470523949047,8.65513359555045
51.4544387839093,8.68703519084252
52.1716857555435,8.71829212354567
52.8989307609815,8.74891212003788
53.6363131673924,8.77890307327924
54.383974284648,8.80827302751912
55.142057392403,8.83703016346628
55.9107077675529,8.86518278393096
56.6900727120743,8.89273929994598
57.4803015812536,8.91970821737171
58.2815458123085,8.94609812398846
59.0939589534097,8.97191767707788
59.9176966931062,8.99717559149351
60.7529168901607,9.02188062821919
61.5997796038017,9.04604158341276
62.4584471243962,9.06966727793122
63.3290840045511,9.09276654733233
64.2118570906476,9.11534823234681
65.1069355548146,9.13742116981416
66.0144909273489,9.15899418407432
66.9346971295866,9.18007607880696
67.8677305072329,9.20067562930889
68.8137698641567,9.2208015752002
69.7729964966554,9.24046261354877
70.7455942281988,9.25966739240248
71.7317494446561,9.27842450471818
72.7316511300146,9.296742482676
73.7454909025955,9.31462979236761
74.7734630517759,9.33209482884641
75.8157645752211,9.34914591152809
76.8725952166374,9.36579127992929
77.9441575040495,9.38203908973263
79.0306567886135,9.3978974091658
80.1323012839689,9.41337421568309
81.2493021061405,9.42847739293711
82.3818733139961,9.4432147280292
83.5302319502678,9.45759390902663
84.6945980831458,9.47162252273536
85.8751948484517,9.48530805271663
87.0722484923992,9.4986578775368
88.2859884149515,9.51167926923896
89.516647213783,9.52437939202607
90.7644607288536,9.53676530114503
92.0296680876042,9.54884394196147
93.3125117507825,9.5606221492155
94.6132375589077,9.57210664644862
95.9320947793824,9.58330404559261
97.2693361542618,9.59422084671128
98.6252179486878,9.60486343788619
100,9.6152380952381
1 freqs ampl
2 0.1 0.00124984376952881
3 0.101393945767529 0.00125685750941772
4 0.102807322383087 0.00126406813733846
5 0.104240400702156 0.00127148117998881
6 0.105693455355799 0.00127910231918553
7 0.107166764803286 0.00128693739621714
8 0.10866061138546 0.00129499241631883
9 0.110175281378839 0.00130327355327284
10 0.111711065050482 0.00131178715413789
11 0.113268256713615 0.00132053974411124
12 0.114847154784029 0.00132953803152714
13 0.116448061837269 0.00133878891299542
14 0.118071284666619 0.00134829947868417
15 0.119717134341897 0.00135807701775069
16 0.121385926269063 0.00136812902392454
17 0.123077980250667 0.00137846320124731
18 0.124793620547131 0.00138908746997326
19 0.126533175938894 0.00140000997263542
20 0.128296979789415 0.00141123908028183
21 0.130085370109057 0.0014227833988865
22 0.131898689619867 0.00143465177594028
23 0.13373728582125 0.00144685330722638
24 0.135601511056563 0.00145939734378589
25 0.137491722580642 0.00147229349907857
26 0.139408282628258 0.0014855516563444
27 0.14135155848354 0.00149918197617139
28 0.143321922550357 0.00151319490427561
29 0.14531975242369 0.00152760117949929
30 0.147345430961984 0.001542411842033
31 0.149399346360526 0.00155763824186835
32 0.151481892225835 0.00157329204748753
33 0.15359346765109 0.00158938525479641
34 0.155734477292613 0.00160593019630794
35 0.157905331447418 0.00162293955058283
36 0.160106446131832 0.00164042635193486
37 0.162338243161228 0.00165840400040801
38 0.164601150230855 0.00167688627203319
39 0.166895600997802 0.00169588732937223
40 0.169222035164104 0.00171542173235734
41 0.171580898561001 0.0017355044494341
42 0.17397264323438 0.00175615086901664
43 0.176397727531404 0.00177737681126359
44 0.178856616188346 0.00179919854018389
45 0.18134978041965 0.00182163277608158
46 0.183877698008233 0.00184469670834908
47 0.186440853397049 0.00186840800861865
48 0.189039737781922 0.00189278484428217
49 0.191674849205682 0.00191784589238927
50 0.194346692653602 0.00194361035393467
51 0.19705578015018 0.00197009796854536
52 0.199802630857255 0.00199732902957889
53 0.202587771173502 0.0020253243996442
54 0.205411734835306 0.00205410552655685
55 0.208275063019052 0.00208369445974066
56 0.211178304444824 0.00211411386708823
57 0.214122015481573 0.00214538705229329
58 0.217106760253727 0.00217753797266776
59 0.220133110749303 0.0022105912574573
60 0.223201646929523 0.00224457222666906
61 0.226312956839953 0.00227950691042593
62 0.229467636723194 0.00231542206886212
63 0.232666291133146 0.00235234521257481
64 0.235909533050864 0.00239030462364776
65 0.239197984002024 0.00242932937726244
66 0.242532274176035 0.00246944936391327
67 0.245913042546805 0.00251069531224364
68 0.249340936995188 0.00255309881252008
69 0.25281661443315 0.00259669234076224
70 0.256340740929651 0.00264150928354692
71 0.259913991838293 0.00268758396350498
72 0.263537051926739 0.00273495166553027
73 0.267210615507943 0.00278364866372035
74 0.270935386573205 0.00283371224906948
75 0.274712078927081 0.00288518075793453
76 0.278541416324177 0.00293809360129545
77 0.282424132607843 0.00299249129483226
78 0.286360971850812 0.0030484154898412
79 0.290352688497781 0.0031059090050133
80 0.294400047510004 0.00316501585909936
81 0.298503824511873 0.0032257813044857
82 0.302664805939569 0.00328825186170615
83 0.306883789191764 0.00335247535491597
84 0.311161582782436 0.00341850094835445
85 0.315499006495809 0.00348637918382354
86 0.319896891543454 0.00355616201921048
87 0.324356080723581 0.00362790286808339
88 0.328877428582551 0.00370165664038944
89 0.333461801578636 0.00377747978428603
90 0.338110078248068 0.00385543032913618
91 0.342823149373397 0.0039355679297004
92 0.347601918154198 0.00401795391155783
93 0.352447300380159 0.00410265131779068
94 0.357360224606579 0.00418972495696669
95 0.362341632332315 0.00427924145245532
96 0.367392478180207 0.00437126929311439
97 0.372513730080021 0.0044658788853849
98 0.377706369453937 0.00456314260683253
99 0.382971391404628 0.00466313486117581
100 0.388309804905961 0.00476593213484159
101 0.393722632996348 0.0048716130550896
102 0.399210912974805 0.00498025844974942
103 0.404775696599732 0.00509195140861367
104 0.410418050290471 0.00520677734653297
105 0.416139055331671 0.0053248240682591
106 0.421939808080503 0.00544618183508416
107 0.427821420176762 0.0055709434333248
108 0.433785018755899 0.00569920424470177
109 0.439831746665023 0.0058310623186666
110 0.445962762681909 0.00596661844672819
111 0.45217924173707 0.00610597623883409
112 0.45848237513891 0.00624924220186199
113 0.464873370802026 0.00639652582027892
114 0.471353453478691 0.00654793963902696
115 0.47792386499356 0.00670359934869576
116 0.48458586448165 0.00686362387304391
117 0.491340728629636 0.00702813545893254
118 0.498189751920516 0.00719725976873649
119 0.505134246881677 0.00737112597529983
120 0.512175544336424 0.00754986685950441
121 0.519314993659021 0.0077336189105217
122 0.526553963033276 0.00792252242882021
123 0.533893839714735 0.00811672163200241
124 0.541336030296538 0.0083163647635472
125 0.548881960978967 0.00852160420453555
126 0.556533077842765 0.00873259658843936
127 0.564290847126254 0.00894950291905528
128 0.572156755506325 0.00917248869166733
129 0.580132310383338 0.00940172401752446
130 0.588219040169996 0.00963738375172108
131 0.596418494584246 0.00987964762457099
132 0.604732244946265 0.0101287003765674
133 0.613161884479578 0.0103847318970236
134 0.621709028616383 0.0106479373664921
135 0.630375315307128 0.0109185174030614
136 0.6391624053344 0.011196678212632
137 0.648071982631197 0.011482631743278
138 0.657105754603627 0.0117765958437987
139 0.666265452458115 0.0120787944265724
140 0.675552831533164 0.0123894576348226
141 0.684969671635745 0.0127088220144133
142 0.69451777738237 0.0130371306902897
143 0.704198978544929 0.0133746335476852
144 0.71401513040134 0.0137215874182187
145 0.723968114091088 0.0140782562710072
146 0.734059836975721 0.0144449114089228
147 0.744292233004376 0.0148218316701279
148 0.754667263084389 0.0152093036350209
149 0.765186915457082 0.0156076218387322
150 0.775853206078784 0.0160170889893121
151 0.786668179007158 0.0164380161917532
152 0.797633906792928 0.0168707231779958
153 0.808752490877045 0.017315538543068
154 0.820026061993413 0.0177727999875122
155 0.831456780577207 0.0182428545662572
156 0.843046837178897 0.0187260589440964
157 0.854798452884041 0.0192227796579347
158 0.866713879738923 0.0197333933859734
159 0.878795401182132 0.0202582872240027
160 0.891045332482152 0.0207978589689765
161 0.903466021181053 0.0213525174100481
162 0.916059847544371 0.0219226826272462
163 0.92882922501725 0.0225087862979786
164 0.941776600686952 0.0231112720115502
165 0.954904455751808 0.0237305955918877
166 0.968215305996708 0.0243672254286671
167 0.981711702275219 0.0250216428170423
168 0.995396230998424 0.0256943423061784
169 1.00927151463057 0.026385832056794
170 1.02334021219164 0.0270966342079253
171 1.03760501976691 0.0278272852531217
172 1.05206867102362 0.0285783364262923
173 1.06673393773486 0.0293503540974215
174 1.08160363031071 0.0301439201783789
175 1.09668059833687 0.0309596325390481
176 1.1119677311207 0.0317981054340067
177 1.12746795824495 0.0326599699399889
178 1.14318425012915 0.0335458744043682
179 1.15911961859889 0.0344564849048982
180 1.17527711746295 0.0353924857209554
181 1.19165984309856 0.036354579816527
182 1.20827093504478 0.037343489335192
183 1.22511357660412 0.0383599561073468
184 1.24219099545262 0.0394047421699245
185 1.25950646425836 0.0404786302988653
186 1.27706330130864 0.0415824245545925
187 1.2948648711459 0.0427169508407526
188 1.31291458521247 0.043883057476478
189 1.33121590250431 0.0450816157824338
190 1.34977233023394 0.0463135206809095
191 1.36858742450252 0.0475796913102165
192 1.38766479098131 0.0488810716536558
193 1.40700808560268 0.0502186311833147
194 1.42662101526074 0.0515933655189563
195 1.44650733852165 0.0530062971022605
196 1.46667086634397 0.0544584758866766
197 1.48711546280895 0.0559509800431428
198 1.50784504586105 0.0574849166819288
199 1.52886358805873 0.0590614225908512
200 1.55017511733577 0.0606816649901108
201 1.57178371777316 0.0623468423039935
202 1.59369353038178 0.0640581849496757
203 1.61590875389592 0.0658169561433637
204 1.63843364557798 0.0676244527239964
205 1.66127252203429 0.0694820059947267
206 1.68442976004232 0.071390982582396
207 1.70790979738943 0.0733527853151986
208 1.73171713372335 0.0753688541187309
209 1.75585633141447 0.0774406669306033
210 1.78033201643011 0.0795697406337812
211 1.80514887922111 0.0817576320088108
212 1.83031167562061 0.0840059387050673
213 1.85582522775552 0.086316300231147
214 1.88169442497056 0.0886903989645098
215 1.90792422476527 0.0911299611804565
216 1.93451965374405 0.0936367581005048
217 1.96148580857943 0.0962126069602082
218 1.98882785698881 0.0988593720964329
219 2.01655103872475 0.101578966054087
220 2.04466066657912 0.104373350712264
221 2.07316212740123 0.107244538429732
222 2.10206088313016 0.110194593209675
223 2.13136247184144 0.113225631883545
224 2.16107250880838 0.116339825313863
225 2.19119668757815 0.119539399615746
226 2.22174078106288 0.12282663739692
227 2.25271064264598 0.12620387901591
228 2.28411220730382 0.129673523858072
229 2.31595149274315 0.133238031629055
230 2.34823460055428 0.136899923665257
231 2.38096771738036 0.140661784260755
232 2.41415711610302 0.144526262010144
233 2.44780915704444 0.148496071166648
234 2.48193028918626 0.152573993014796
235 2.51652705140539 0.156762877256901
236 2.55160607372718 0.161065643412464
237 2.58717407859592 0.1654852822296
238 2.62323788216312 0.170024857107441
239 2.65980439559376 0.174687505528422
240 2.69688062639069 0.179476440499244
241 2.73447367973758 0.184394951999205
242 2.77259075986048 0.189446408434499
243 2.81123917140846 0.194634258096962
244 2.85042632085343 0.199962030625633
245 2.89015971790951 0.205433338469367
246 2.93044697697214 0.211051878348636
247 2.97129581857733 0.216821432714486
248 3.01271407088116 0.222745871202514
249 3.05470967115997 0.228829152079549
250 3.09729066733141 0.235075323680605
251 3.14046521949675 0.241488525833486
252 3.18424160150465 0.24807299126828
253 3.22862820253673 0.254833047008791
254 3.27363352871525 0.261773115742795
255 3.31926620473319 0.268897717167807
256 3.36553497550709 0.276211469308868
257 3.41244870785289 0.283719089804643
258 3.46001639218511 0.291425397157943
259 3.50824714423979 0.299335311946534
260 3.55715020682139 0.307453857989918
261 3.60673495157403 0.315786163467503
262 3.65701088077749 0.324337461983377
263 3.70798762916817 0.333113093572643
264 3.75967496578547 0.342118505644029
265 3.81208279584389 0.351359253853236
266 3.86522116263126 0.36084100290123
267 3.9191002494334 0.370569527251405
268 3.97373038148561 0.380550711759293
269 4.02912202795135 0.390790552208207
270 4.08528580392856 0.401295155743936
271 4.1422324724839 0.412070741201301
272 4.19997294671531 0.423123639315149
273 4.25851829184341 0.434460292807991
274 4.31787972733202 0.446087256346283
275 4.37806862903817 0.45801119635701
276 4.43909653139217 0.470238890695935
277 4.50097512960805 0.482777228158622
278 4.56371628192476 0.495633207825031
279 4.62733201187869 0.508813938228191
280 4.6918345106078 0.522326636337199
281 4.75723613918789 0.53617862634452
282 4.82354943100147 0.550377338247268
283 4.89078709413959 0.564930306211964
284 4.95896201383721 0.579845166711924
285 5.02808725494248 0.595129656426328
286 5.09817606442043 0.610791609889683
287 5.1692418738916 0.626838956880323
288 5.24129830220606 0.643279719536317
289 5.31435915805324 0.660122009187079
290 5.38843844260822 0.677374022888832
291 5.46355035221488 0.695044039651974
292 5.5397092811064 0.713140416348375
293 5.61692982416381 0.731671583286571
294 5.69522677971282 0.750646039442876
295 5.77461515235982 0.770072347336476
296 5.85511015586724 0.789959127536677
297 5.93672721606913 0.810315052790655
298 6.01948197382727 0.831148841760246
299 6.10339028802862 0.852469252356602
300 6.18846823862439 0.874285074661867
301 6.27473212971158 0.896605123427405
302 6.36219849265749 0.919438230138643
303 6.45088408926769 0.94279323463703
304 6.54080591499826 0.966678976290355
305 6.63198120221267 0.991104284703253
306 6.72442742348425 1.01607796996062
307 6.81816229494448 1.04160881239745
308 6.91320377967813 1.06770555188962
309 7.00957009116562 1.09437687666126
310 7.10727969677342 1.12163141160536
311 7.20635132129305 1.14947770611576
312 7.3068039505295 1.17792422142984
313 7.40865683493956 1.20697931748277
314 7.51192949332097 1.23665123927596
315 7.61664171655289 1.26694810276361
316 7.72281357138865 1.29787788026358
317 7.83046540430119 1.3294483854004
318 7.93961784538228 1.36166725759032
319 8.05029181229598 1.3945419460805
320 8.16250851428723 1.42807969355667
321 8.27628945624635 1.46228751933592
322 8.39165644283016 1.49717220216369
323 8.50863158264058 1.53274026263664
324 8.62723729246145 1.5689979452756
325 8.74749630155442 1.60595120027567
326 8.86943165601471 1.643605664963
327 8.99306672318762 1.68196664499107
328 9.11842519614657 1.72103909531164
329 9.24553109823358 1.76082760095885
330 9.374408787663 1.80133635768771
331 9.50508296218951 1.84256915251122
332 9.63757866384109 1.88452934418316
333 9.771921283718 1.92721984367676
334 9.90813656685867 1.97064309471195
335 10.0462506171734 2.01480105438701
336 10.1862899024469 2.05969517397301
337 10.3282812594103 2.10532637993194
338 10.4722518988843 2.15169505522199
339 10.6182294109938 2.19880102095579
340 10.7662417704549 2.24664351847938
341 10.9163173419361 2.29522119194176
342 11.0684848854941 2.34453207142641
343 11.2227735620851 2.39457355671776
344 11.3792129391532 2.44534240177657
345 11.5378329962966 2.49683469999911
346 11.6986641310131 2.54904587033564
347 11.8617371645248 2.60197064434355
348 12.0270833476851 2.6556030542507
349 12.1947343669674 2.70993642210359
350 12.3647223505371 2.76496335007402
351 12.5370798744092 2.82067571199669
352 12.7118399686903 2.87706464620789
353 12.8890361239089 2.93412054975349
354 13.0687022974335 2.9918330740313
355 13.2508729199795 3.05019112192989
356 13.4355829022083 3.10918284652194
357 13.6228676414165 3.16879565136619
358 13.8127630283201 3.22901619246727
359 14.0053054539322 3.28983038193766
360 14.2005318165368 3.35122339340032
361 14.3984795287601 3.41317966916477
362 14.5991865247398 3.47568292920294
363 14.8026912673951 3.53871618194413
364 15.0090327557974 3.60226173690188
365 15.2182505326439 3.66630121913765
366 15.4303846918356 3.73081558555893
367 15.645475886161 3.79578514304128
368 15.8635653350859 3.86118956835606
369 16.0846948326536 3.92700792987717
370 16.3089067554933 3.99321871103223
371 16.5362440709418 4.059799835455
372 16.766750345277 4.12672869378805
373 17.0004697520672 4.19398217207637
374 17.2374470806362 4.26153668168462
375 17.4777277446468 4.32936819066333
376 17.7213577908036 4.3974522564814
377 17.9683839076772 4.4657640600358
378 18.2188534346517 4.53427844084184
379 18.4728143709963 4.60296993330202
380 18.7303153850644 4.67181280394469
381 18.9914058236194 4.74078108951949
382 19.256135721292 4.80984863583088
383 19.5245558101686 4.87898913718794
384 19.7967175295133 4.94817617634452
385 20.0726730356257 5.01738326480138
386 20.3524752118358 5.08658388334002
387 20.6361776786386 5.15575152265655
388 20.9238348039698 5.22485972396316
389 21.2155017136245 5.29388211942536
390 21.5112343018217 5.36279247230368
391 21.8110892419152 5.43156471667044
392 22.1151239972549 5.50017299657441
393 22.4233968321984 5.56859170452924
394 22.7359668232772 5.63679551920541
395 23.0528938705171 5.70475944220955
396 23.3742387089181 5.77245883384028
397 23.7000629200933 5.83986944771488
398 24.0304289440697 5.90696746416749
399 24.3654000912547 5.97372952232568
400 24.7050405545683 6.04013275077966
401 25.049415421745 6.10615479676495
402 25.3985906878073 6.17177385378773
403 25.752633267712 6.23696868762925
404 26.1116110091746 6.30171866067413
405 26.4755927056707 6.36600375451527
406 26.8446481096197 6.42980459079649
407 27.2188479457518 6.49310245026205
408 27.5982639246618 6.55587928999065
409 27.9829687565512 6.61811775879929
410 28.3730361651621 6.67980121081061
411 28.7685409019059 6.74091371718493
412 29.169558760188 6.80144007602565
413 29.5761665899325 6.8613658204741
414 29.9884423123103 6.92067722501673
415 30.4064649346707 6.97936131003399
416 30.8303145656828 7.03740584462663
417 31.260072430687 7.09479934776086
418 31.6958208872612 7.15153108777898
419 32.1376434410028 7.20759108032728
420 32.5856247615323 7.26297008475719
421 33.0398506987186 7.31765959905996
422 33.5004082991313 7.37165185339816
423 33.9673858227221 7.42493980230096
424 34.4408727597382 7.47751711559223
425 34.9209598478727 7.52937816812287
426 35.4077390896527 7.5805180283804
427 35.9013037700707 7.63093244605004
428 36.4017484744614 7.68061783860244
429 36.9091691066278 7.72957127698342
430 37.4236629072198 7.77779047048138
431 37.9453284723694 7.8252737508475
432 38.4742657725851 7.87202005574321
433 39.0105761719099 7.91802891158854
434 39.554362447347 7.96330041588347
435 40.105728808555 8.00783521907301
436 40.6647809178186 8.05163450602489
437 41.2316259102975 8.09469997718672
438 41.8063724145576 8.13703382948712
439 42.3891305733878 8.17863873704334
440 42.980012064908 8.21951783173482
441 43.5791301239703 8.25967468369993
442 44.1865995638594 8.29911328181031
443 44.8025367982949 8.33783801417395
444 45.4270598637405 8.37585364871619
445 46.060288442024 8.41316531388383
446 46.7023438832734 8.44977847951546
447 47.3533492291712 8.48569893791774
448 48.0134292365345 8.52093278518449
449 48.6827104012228 8.55548640279275
450 49.3613209823792 8.58936643950677
451 50.0493910270095 8.62257979361857
452 50.7470523949047 8.65513359555045
453 51.4544387839093 8.68703519084252
454 52.1716857555435 8.71829212354567
455 52.8989307609815 8.74891212003788
456 53.6363131673924 8.77890307327924
457 54.383974284648 8.80827302751912
458 55.142057392403 8.83703016346628
459 55.9107077675529 8.86518278393096
460 56.6900727120743 8.89273929994598
461 57.4803015812536 8.91970821737171
462 58.2815458123085 8.94609812398846
463 59.0939589534097 8.97191767707788
464 59.9176966931062 8.99717559149351
465 60.7529168901607 9.02188062821919
466 61.5997796038017 9.04604158341276
467 62.4584471243962 9.06966727793122
468 63.3290840045511 9.09276654733233
469 64.2118570906476 9.11534823234681
470 65.1069355548146 9.13742116981416
471 66.0144909273489 9.15899418407432
472 66.9346971295866 9.18007607880696
473 67.8677305072329 9.20067562930889
474 68.8137698641567 9.2208015752002
475 69.7729964966554 9.24046261354877
476 70.7455942281988 9.25966739240248
477 71.7317494446561 9.27842450471818
478 72.7316511300146 9.296742482676
479 73.7454909025955 9.31462979236761
480 74.7734630517759 9.33209482884641
481 75.8157645752211 9.34914591152809
482 76.8725952166374 9.36579127992929
483 77.9441575040495 9.38203908973263
484 79.0306567886135 9.3978974091658
485 80.1323012839689 9.41337421568309
486 81.2493021061405 9.42847739293711
487 82.3818733139961 9.4432147280292
488 83.5302319502678 9.45759390902663
489 84.6945980831458 9.47162252273536
490 85.8751948484517 9.48530805271663
491 87.0722484923992 9.4986578775368
492 88.2859884149515 9.51167926923896
493 89.516647213783 9.52437939202607
494 90.7644607288536 9.53676530114503
495 92.0296680876042 9.54884394196147
496 93.3125117507825 9.5606221492155
497 94.6132375589077 9.57210664644862
498 95.9320947793824 9.58330404559261
499 97.2693361542618 9.59422084671128
500 98.6252179486878 9.60486343788619
501 100 9.6152380952381

BIN
matlab/matlab.zip Normal file
View File

Binary file not shown.

395
matlab/matlab/comp_filters_ligo.m Normal file
View File

@ -0,0 +1,395 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
freqs = logspace(-3, 0, 1000);
% Specifications
% The specifications for the filters are:
% 1. From $0$ to $0.008\text{ Hz}$,the magnitude of the filter’s transfer function should be less than or equal to $8 \times 10^{-3}$
% 2. From $0.008\text{ Hz}$ to $0.04\text{ Hz}$, it attenuates the input signal proportional to frequency cubed
% 3. Between $0.04\text{ Hz}$ and $0.1\text{ Hz}$, the magnitude of the transfer function should be less than 3
% 4. Above $0.1\text{ Hz}$, the maximum of the magnitude of the complement filter should be as close to zero as possible. In our system, we would like to have the magnitude of the complementary filter to be less than $0.1$. As the filters obtained in cite:hua05_low_ligo have a magnitude of $0.045$, we will set that as our requirement
% The specifications are translated in upper bounds of the complementary filters are shown on figure [[fig:ligo_specifications]].
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1)
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1)
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2)
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-4, 10]);
legend('location', 'northeast');
% FIR Filter
% We here try to implement the FIR complementary filter synthesis as explained in cite:hua05_low_ligo.
% For that, we use the [[http://cvxr.com/cvx/][CVX matlab Toolbox]].
% We setup the CVX toolbox and use the =SeDuMi= solver.
cvx_startup;
cvx_solver sedumi;
% We define the frequency vectors on which we will constrain the norm of the FIR filter.
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;
% We then define the order of the FIR filter.
n = 512;
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]))];
% We run the 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);
% #+RESULTS:
% #+begin_example
% 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);
% #+end_example
% Finally, we compute the filter response over the frequency vector defined and the result is shown on figure [[fig:fir_filter_ligo]] which is very close to the filters obtain in cite:hua05_low_ligo.
w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
figure;
ax1 = subplot(2,1,1);
hold on;
plot(w, abs(H), 'k-');
plot(w, abs(1-H), 'k--');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 5]);
ax2 = subplot(2,1,2);
hold on;
plot(w, 180/pi*angle(H), 'k-');
plot(w, 180/pi*angle(1-H), 'k--');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-540:90:360]);
linkaxes([ax1,ax2],'x');
xlim([1e-3, 1]);
xticks([0.01, 0.1, 1, 10, 100, 1000]);
% Weights
% We design weights that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters.
% 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 minimum phase
% The bode plot of the weights is shown on figure [[fig:ligo_weights]].
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));
n = 20; Rp = 1; Wp = 2*pi*0.102;
[b,a] = cheby1(n, Rp, Wp, 'high', 's');
wL = 0.04*tf(a, b);
wL = 1/wL;
wL = minreal(ss(wL));
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|w_H|^{-1}$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|w_L|^{-1}$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec.');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
plot([0.1, 10], [0.045, 0.045], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast');
% H-Infinity Synthesis
% We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion]].
P = [0 wL;
wH -wH;
1 0];
% And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
% #+RESULTS:
% #+begin_example
% [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
% #+end_example
% The high pass filter is defined as $H_H = 1 - H_L$.
Hh = 1 - Hl;
Hh = minreal(Hh);
Hl = minreal(Hl);
% The size of the filters is shown below.
size(Hh), size(Hl)
% #+RESULTS:
% #+begin_example
% 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.
% #+end_example
% The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_ligo_results]].
figure;
hold on;
set(gca,'ColorOrderIndex',1);
plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$');
set(gca,'ColorOrderIndex',1);
plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1);
plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',2);
plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$');
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast');
% Compare FIR and H-Infinity Filters
% Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one obtained with the $\mathcal{H}_\infty$ synthesis on figure [[fig:comp_fir_ligo_hinf]].
figure;
ax1 = subplot(2,1,1);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',1);
plot(w, abs(H), '--');
set(gca,'ColorOrderIndex',2);
plot(w, abs(1-H), '--');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([1e-3, 10]);
ax2 = subplot(2,1,2);
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$\mathcal{H}_\infty$ filters');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'HandleVisibility', 'off');
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*angle(H), '--', 'DisplayName', 'FIR filters');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*angle(1-H), '--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-540:90:360]);
legend('location', 'northeast');
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.001, 0.01, 0.1, 1]);

144
matlab/matlab/h_inf_synthesis_complementary_filters.m Normal file
View File

@ -0,0 +1,144 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
freqs = logspace(-1, 3, 1000);
% Design of Weighting Function
% 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
% \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 [[fig:weight_formula]].
% #+name: fig:weight_formula
% #+caption: Amplitude of the proposed formula for the weighting functions
% [[file:figs-tikz/weight_formula.png]]
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W2 = (((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;
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-4, 20]);
xticks([0.1, 1, 10, 100, 1000]);
legend('location', 'northeast');
% H-Infinity Synthesis
% We define the generalized plant $P$ on matlab.
P = [W1 -W1;
0 W2;
1 0];
% And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command.
[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
% #+RESULTS:
% #+begin_example
% [H2, ~, 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.1000 < gamma <= 1050.0000
% gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
% 1.050e+03 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 525.050 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 262.575 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 131.337 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
% 65.719 2.8e+01 2.4e-07 4.1e+00 -9.5e-14 0.0000 p
% 32.909 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 16.505 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
% 8.302 2.8e+01 2.4e-07 4.1e+00 -7.2e-14 0.0000 p
% 4.201 2.8e+01 2.4e-07 4.1e+00 -2.5e-25 0.0000 p
% 2.151 2.7e+01 2.4e-07 4.1e+00 -3.8e-14 0.0000 p
% 1.125 2.6e+01 2.4e-07 4.1e+00 -5.4e-24 0.0000 p
% 0.613 2.3e+01 -3.7e+01# 4.1e+00 0.0e+00 0.0000 f
% 0.869 2.6e+01 -3.7e+02# 4.1e+00 0.0e+00 0.0000 f
% 0.997 2.6e+01 -1.1e+04# 4.1e+00 0.0e+00 0.0000 f
% 1.061 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.029 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.013 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.005 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.001 2.6e+01 -3.1e+04# 4.1e+00 -3.8e-14 0.0000 f
% 1.003 2.6e+01 -2.8e+05# 4.1e+00 0.0e+00 0.0000 f
% 1.004 2.6e+01 2.4e-07 4.1e+00 -5.8e-24 0.0000 p
% 1.004 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% Gamma value achieved: 1.0036
% #+end_example
% We then define the high pass filter $H_1 = 1 - H_2$. The bode plot of both $H_1$ and $H_2$ is shown on figure [[fig:hinf_filters_results]].
H1 = 1 - H2;
% Obtained Complementary Filters
% The obtained complementary filters are shown on figure [[fig:hinf_filters_results]].
figure;
ax1 = subplot(2,1,1);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$');
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-4, 20]);
legend('location', 'northeast');
ax2 = subplot(2,1,2);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
set(gca,'ColorOrderIndex',2)
plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:90:360]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.1, 1, 10, 100, 1000]);

95
matlab/matlab/three_comp_filters.m Normal file
View File

@ -0,0 +1,95 @@
%% Clear Workspace and Close figures
clear; close all; clc;
%% Intialize Laplace variable
s = zpk('s');
freqs = logspace(-2, 4, 1000);
% Weights
% First we define the weights.
n = 2; w0 = 2*pi*1; G0 = 1/10; G1 = 1000; Gc = 1/2;
W1 = (((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;
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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2;
W3 = (((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;
figure;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
xticks([0.01, 0.1, 1, 10, 100, 1000]);
legend('location', 'northeast');
% H-Infinity Synthesis
% Then we create the generalized plant =P=.
P = [W1 -W1 -W1;
0 W2 0 ;
0 0 W3;
1 0 0];
% And we do the $\mathcal{H}_\infty$ synthesis.
[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
% Obtained Complementary Filters
% The obtained filters are:
H2 = tf(H(1));
H3 = tf(H(2));
H1 = 1 - H2 - H3;
figure;
ax1 = subplot(2,1,1);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
set(gca,'ColorOrderIndex',2)
plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$');
set(gca,'ColorOrderIndex',3)
plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$');
set(gca,'ColorOrderIndex',1)
plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$');
set(gca,'ColorOrderIndex',2)
plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$');
set(gca,'ColorOrderIndex',3)
plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-4, 20]);
legend('location', 'northeast');
ax2 = subplot(2,1,2);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))));
set(gca,'ColorOrderIndex',2)
plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))));
set(gca,'ColorOrderIndex',3)
plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz'))));
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:90:360]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.1, 1, 10, 100, 1000]);

171
matlab/ref.bib Normal file
View File

@ -0,0 +1,171 @@
@article{collette15_sensor_fusion_method_high_perfor,
author = {C. Collette and F. Matichard},
title = {Sensor Fusion Methods for High Performance Active Vibration Isolation Systems},
journal = {Journal of Sound and Vibration},
volume = {342},
number = {nil},
pages = {1-21},
year = {2015},
doi = {10.1016/j.jsv.2015.01.006},
url = {https://doi.org/10.1016/j.jsv.2015.01.006},
keywords = {},
}
@phdthesis{hua05_low_ligo,
author = {Hua, Wensheng},
school = {stanford university},
title = {Low frequency vibration isolation and alignment system for
advanced LIGO},
year = 2005,
}
@inproceedings{hua04_polyp_fir_compl_filter_contr_system,
author = {Wensheng Hua and Dan B. Debra and Corwin T. Hardham and Brian T. Lantz and Joseph A. Giaime},
title = {Polyphase FIR Complementary Filters for Control Systems},
booktitle = {Proceedings of ASPE Spring Topical Meeting on Control of Precision Systems},
year = {2004},
pages = {109--114},
}
@article{matichard15_seism_isolat_advan_ligo,
author = {Matichard, F and Lantz, B and Mittleman, R and Mason, K and Kissel, J and Abbott, B and Biscans, S and McIver, J and Abbott, R and Abbott, S and others},
title = {Seismic Isolation of Advanced Ligo: Review of Strategy, Instrumentation and Performance},
journal = {Classical and Quantum Gravity},
volume = {32},
number = {18},
pages = {185003},
year = {2015},
publisher = {IOP Publishing},
}
@article{min15_compl_filter_desig_angle_estim,
author = {Min, Hyung Gi and Jeung, Eun Tae},
title = {Complementary Filter Design for Angle Estimation Using Mems Accelerometer and Gyroscope},
journal = {Department of Control and Instrumentation, Changwon National University, Changwon, Korea},
pages = {641--773},
year = {2015},
}
@article{corke04_inert_visual_sensin_system_small_auton_helic,
author = {Peter Corke},
title = {An Inertial and Visual Sensing System for a Small Autonomous Helicopter},
journal = {Journal of Robotic Systems},
volume = {21},
number = {2},
pages = {43-51},
year = {2004},
doi = {10.1002/rob.10127},
url = {https://doi.org/10.1002/rob.10127},
}
@inproceedings{jensen13_basic_uas,
author = {Austin Jensen and Cal Coopmans and YangQuan Chen},
title = {Basics and guidelines of complementary filters for small UAS
navigation},
booktitle = {2013 International Conference on Unmanned Aircraft Systems
(ICUAS)},
year = 2013,
pages = {nil},
doi = {10.1109/icuas.2013.6564726},
url = {https://doi.org/10.1109/icuas.2013.6564726},
month = 5,
}
@inproceedings{pascoal99_navig_system_desig_using_time,
author = {A. Pascoal and I. Kaminer and P. Oliveira},
title = {Navigation System Design Using Time-Varying Complementary Filters},
booktitle = {Guidance, Navigation, and Control Conference and Exhibit},
year = {1999},
pages = {nil},
doi = {10.2514/6.1999-4290},
url = {https://doi.org/10.2514/6.1999-4290},
}
@article{zimmermann92_high_bandw_orien_measur_contr,
author = {M. Zimmermann and W. Sulzer},
title = {High Bandwidth Orientation Measurement and Control Based on Complementary Filtering},
journal = {Robot Control 1991},
pages = {525-530},
year = {1992},
doi = {10.1016/b978-0-08-041276-4.50093-5},
url = {https://doi.org/10.1016/b978-0-08-041276-4.50093-5},
publisher = {Elsevier},
series = {Robot Control 1991},
}
@inproceedings{baerveldt97_low_cost_low_weigh_attit,
author = {A.-J. Baerveldt and R. Klang},
title = {A Low-Cost and Low-Weight Attitude Estimation System for an Autonomous Helicopter},
booktitle = {Proceedings of IEEE International Conference on Intelligent Engineering Systems},
year = {1997},
pages = {nil},
doi = {10.1109/ines.1997.632450},
url = {https://doi.org/10.1109/ines.1997.632450},
month = {-},
}
@article{shaw90_bandw_enhan_posit_measur_using_measur_accel,
author = {F.R. Shaw and K. Srinivasan},
title = {Bandwidth Enhancement of Position Measurements Using Measured
Acceleration},
journal = {Mechanical Systems and Signal Processing},
volume = 4,
number = 1,
pages = {23-38},
year = 1990,
doi = {10.1016/0888-3270(90)90038-m},
url = {https://doi.org/10.1016/0888-3270(90)90038-m},
}
@article{brown72_integ_navig_system_kalman_filter,
author = {R. G. Brown},
title = {Integrated Navigation Systems and Kalman Filtering: a
Perspective},
journal = {Navigation},
volume = 19,
number = 4,
pages = {355-362},
year = 1972,
doi = {10.1002/j.2161-4296.1972.tb01706.x},
url = {https://doi.org/10.1002/j.2161-4296.1972.tb01706.x},
}
@article{matichard15_seism_isolat_advan_ligo,
author = {Matichard, F and Lantz, B and Mittleman, R and Mason, K and
Kissel, J and Abbott, B and Biscans, S and McIver, J and
Abbott, R and Abbott, S and others},
title = {Seismic Isolation of Advanced Ligo: Review of Strategy,
Instrumentation and Performance},
journal = {Classical and Quantum Gravity},
volume = 32,
number = 18,
pages = 185003,
year = 2015,
publisher = {IOP Publishing},
}
@article{mahony08_nonlin_compl_filter_special_orthog_group,
author = {Robert Mahony and Tarek Hamel and Jean-Michel Pflimlin},
title = {Nonlinear Complementary Filters on the Special Orthogonal
Group},
journal = {IEEE Transactions on Automatic Control},
volume = 53,
number = 5,
pages = {1203-1218},
year = 2008,
doi = {10.1109/tac.2008.923738},
url = {https://doi.org/10.1109/tac.2008.923738},
}

16
paper/paper.org
View File

@ -320,7 +320,7 @@ Let's validate the proposed design method of complementary filters with a simple
- the gain of both filters is equal to $10^{-3}$ away from the merging frequency - the gain of both filters is equal to $10^{-3}$ away from the merging frequency
The weighting functions $W_1(s)$ and $W_2(s)$ are designed using eqref:eq:weight_formula. The weighting functions $W_1(s)$ and $W_2(s)$ are designed using eqref:eq:weight_formula.
The parameters used are summarized in table ref:tab:weights_params and the magnitude of the weighting functions is shown in Fig. ref:fig:hinf_synthesis_results. The parameters used are summarized in table ref:tab:weights_params and the magnitude of the weighting functions is shown in Fig. ref:fig:hinf_filters_results.
#+name: tab:weights_params #+name: tab:weights_params
#+caption: Parameters used for $W_1(s)$ and $W_2(s)$ #+caption: Parameters used for $W_1(s)$ and $W_2(s)$
@ -334,16 +334,16 @@ The parameters used are summarized in table ref:tab:weights_params and the magni
| $G_c$ | $0.5$ | $0.5$ | | $G_c$ | $0.5$ | $0.5$ |
| $n$ | $2$ | $3$ | | $n$ | $2$ | $3$ |
The bode plots of the obtained complementary filters are shown in Fig. ref:fig:hinf_synthesis_results and their transfer functions in the Laplace domain are given below. The bode plots of the obtained complementary filters are shown in Fig. ref:fig:hinf_filters_results and their transfer functions in the Laplace domain are given below.
\begin{align*} \begin{align*}
H_1(s) &= \frac{10^{-8} (s+6.6e^9) (s+3450)^2 (s^2 + 49s + 895)}{(s+6.6e^4) (s^2 + 106 s + 3e^3) (s^2 + 72s + 3580)}\\ H_1(s) &= \frac{10^{-8} (s+6.6e^9) (s+3450)^2 (s^2 + 49s + 895)}{(s+6.6e^4) (s^2 + 106 s + 3e^3) (s^2 + 72s + 3580)}\\
H_2(s) &= \frac{(s+6.6e^4) (s+160) (s+4)^3}{(s+6.6e^4) (s^2 + 106 s + 3e^3) (s^2 + 72s + 3580)} H_2(s) &= \frac{(s+6.6e^4) (s+160) (s+4)^3}{(s+6.6e^4) (s^2 + 106 s + 3e^3) (s^2 + 72s + 3580)}
\end{align*} \end{align*}
#+name: fig:hinf_synthesis_results #+name: fig:hinf_filters_results
#+caption: Frequency response of the weighting functions and complementary filters obtained using $\mathcal{H}_\infty$ synthesis #+caption: Frequency response of the weighting functions and complementary filters obtained using $\mathcal{H}_\infty$ synthesis
#+attr_latex: :scale 1 #+attr_latex: :scale 1
[[file:figs/hinf_synthesis_results.pdf]] [[file:figs/hinf_filters_results.pdf]]
** Synthesis of Three Complementary Filters ** Synthesis of Three Complementary Filters
<<sec:hinf_three_comp_filters>> <<sec:hinf_three_comp_filters>>
@ -380,13 +380,13 @@ By choosing $H_1(s) \triangleq 1 - H_2(s) - H_3(s)$, the proposed $\mathcal{H}_\
*** Example of generated complementary filters :ignore: *** Example of generated complementary filters :ignore:
An example is given to validate the method where three sensors are used in different frequency bands (up to $\SI{1}{Hz}$, from $1$ to $\SI{10}{Hz}$ and above $\SI{10}{Hz}$ respectively). An example is given to validate the method where three sensors are used in different frequency bands (up to $\SI{1}{Hz}$, from $1$ to $\SI{10}{Hz}$ and above $\SI{10}{Hz}$ respectively).
Three weighting functions are designed using eqref:eq:weight_formula and shown by dashed curves in Fig. ref:fig:hinf_three_synthesis_results. Three weighting functions are designed using eqref:eq:weight_formula and shown by dashed curves in Fig. ref:fig:three_complementary_filters_results.
The bode plots of the obtained complementary filters are shown in Fig. ref:fig:hinf_three_synthesis_results. The bode plots of the obtained complementary filters are shown in Fig. ref:fig:three_complementary_filters_results.
#+name: fig:hinf_three_synthesis_results #+name: fig:three_complementary_filters_results
#+caption: Frequency response of the weighting functions and three complementary filters obtained using $\mathcal{H}_\infty$ synthesis #+caption: Frequency response of the weighting functions and three complementary filters obtained using $\mathcal{H}_\infty$ synthesis
#+attr_latex: :scale 1 #+attr_latex: :scale 1
[[file:figs/hinf_three_synthesis_results.pdf]] [[file:figs/three_complementary_filters_results.pdf]]
* Application: Design of Complementary Filters used in the Active Vibration Isolation System at the LIGO * Application: Design of Complementary Filters used in the Active Vibration Isolation System at the LIGO
<<sec:application_ligo>> <<sec:application_ligo>>

BIN
paper/paper.pdf
View File

Binary file not shown.

BIN
tikz/figs/comp_fir_ligo_hinf.pdf
View File

Binary file not shown.

BIN
tikz/figs/comp_fir_ligo_hinf.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 48 KiB

487
tikz/figs/comp_fir_ligo_hinf.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 206 KiB

66
tikz/figs/comp_fir_ligo_hinf.tex
View File

@ -1,66 +0,0 @@
\setlength\fwidth{6.5cm}
\setlength\fheight{7cm}
\begin{tikzpicture}
\begin{axis}[%
width=1.0\fwidth,
height=0.60\fheight,
at={(0.0\fwidth, 0.35\fheight)},
scale only axis,
xmode=log,
xmin=0.001,
xmax=1,
xtick={0.001,0.01,0.1,1},
xticklabels={{}},
xminorticks=true,
ymode=log,
ymin=0.002,
ymax=5,
ytick={0.001, 0.01, 0.1, 1, 10},
yminorticks=true,
ylabel={Magnitude},
xminorgrids,
yminorgrids,
legend style={at={(1,0)}, outer sep=2pt, anchor=south east, legend cell align=left, align=left, draw=black, nodes={scale=0.7, transform shape}}
]
\addplot [color=mycolor1, line width=1.5pt]
table [x=freqs, y=Hhm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_hinf.csv};
\addlegendentry{$H_H(s)$ - $\mathcal{H}_\infty$}
\addplot [color=mycolor1, dashed, line width=1.5pt]
table [x=freqs, y=Hhm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_fir.csv};
\addlegendentry{$H_H(s)$ - FIR}
\addplot [color=mycolor2, line width=1.5pt]
table [x=freqs, y=Hlm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_hinf.csv};
\addlegendentry{$H_L(s)$ - $\mathcal{H}_\infty$}
\addplot [color=mycolor2, dashed, line width=1.5pt]
table [x=freqs, y=Hlm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_fir.csv};
\addlegendentry{$H_L(s)$ - FIR}
\end{axis}
\begin{axis}[%
width=1.0\fwidth,
height=0.3\fheight,
at={(0.0\fwidth, 0.0\fheight)},
scale only axis,
xmode=log,
xmin=0.001,
xmax=1,
xtick={0.001, 0.01, 0.1, 1},
xminorticks=true,
xlabel={Frequency [Hz]},
ymin=-180,
ymax=180,
ytick={-180, -90, 0, 90, 180},
ylabel={Phase [deg]},
xminorgrids,
]
\addplot [color=mycolor1, line width=1.5pt, forget plot]
table [x=freqs, y=Hhp, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_hinf.csv};
\addplot [color=mycolor1, dashed, line width=1.5pt, forget plot]
table [x=freqs, y=Hhp, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_fir.csv};
\addplot [color=mycolor2, line width=1.5pt, forget plot]
table [x=freqs, y=Hlp, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_hinf.csv};
\addplot [color=mycolor2, dashed, line width=1.5pt, forget plot]
table [x=freqs, y=Hlp, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matcomp_ligo_fir.csv};
\end{axis}
\end{tikzpicture}

BIN
tikz/figs/fir_filter_ligo.pdf Normal file
View File

Binary file not shown.

BIN
tikz/figs/hinf_filters_results.pdf Normal file
View File

Binary file not shown.

BIN
tikz/figs/hinf_synthesis_ligo_results.pdf Normal file
View File

Binary file not shown.

BIN
tikz/figs/hinf_synthesis_results.pdf
View File

Binary file not shown.

BIN
tikz/figs/hinf_synthesis_results.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 36 KiB

408
tikz/figs/hinf_synthesis_results.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 128 KiB

67
tikz/figs/hinf_synthesis_results.tex
View File

@ -1,67 +0,0 @@
\setlength\fwidth{6.5cm}
\setlength\fheight{6cm}
\begin{tikzpicture}
\begin{axis}[%
width=1.0\fwidth,
height=0.5\fheight,
at={(0.0\fwidth, 0.5\fheight)},
scale only axis,
xmode=log,
xmin=0.1,
xmax=1000,
xtick={0.1, 1, 10, 100, 1000},
xticklabels={{}},
xminorticks=true,
ymode=log,
ymin=0.0005,
ymax=20,
ytick={0.001, 0.01, 0.1, 1, 10},
yminorticks=true,
ylabel={Magnitude},
xminorgrids,
yminorgrids,
]
\addplot [color=mycolor1, line width=1.5pt, forget plot]
table [x=freqs, y=H1, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_filters_results.csv};
\addplot [color=mycolor2, line width=1.5pt, forget plot]
table [x=freqs, y=H2, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_filters_results.csv};
\addplot [color=mycolor1, dashed, line width=1.5pt, forget plot]
table [x=freqs, y=W1, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_weights.csv};
\addplot [color=mycolor2, dashed, line width=1.5pt, forget plot]
table [x=freqs, y=W2, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_weights.csv};
\end{axis}
\begin{axis}[%
width=1.0\fwidth,
height=0.45\fheight,
at={(0.0\fwidth, 0.0\fheight)},
scale only axis,
xmode=log,
xmin=0.1,
xmax=1000,
xtick={0.1, 1, 10, 100, 1000},
xminorticks=true,
xlabel={Frequency [Hz]},
ymin=-200,
ymax=200,
ytick={-180, -90, 0, 90, 180},
ylabel={Phase [deg]},
xminorgrids,
legend style={at={(1,1.1)}, outer sep=2pt , anchor=north east, legend cell align=left, align=left, draw=black, nodes={scale=0.7, transform shape}},
]
\addlegendimage{color=mycolor1, dashed, line width=1.5pt}
\addlegendentry{$W_1^{-1}$};
\addlegendimage{color=mycolor2, dashed, line width=1.5pt}
\addlegendentry{$W_2^{-1}$};
\addplot [color=mycolor1, line width=1.5pt]
table [x=freqs, y=H1p, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_filters_results.csv};
\addlegendentry{$H_1$};
\addplot [color=mycolor2, line width=1.5pt]
table [x=freqs, y=H2p, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_filters_results.csv};
\addlegendentry{$H_2$};
\end{axis}
\end{tikzpicture}

BIN
tikz/figs/hinf_three_synthesis_results.pdf
View File

Binary file not shown.

BIN
tikz/figs/hinf_three_synthesis_results.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 52 KiB

442
tikz/figs/hinf_three_synthesis_results.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 196 KiB

80
tikz/figs/hinf_three_synthesis_results.tex
View File

@ -1,80 +0,0 @@
\setlength\fwidth{6.5cm}
\setlength\fheight{6cm}
\begin{tikzpicture}
\begin{axis}[%
width=1.0\fwidth,
height=0.55\fheight,
at={(0.0\fwidth, 0.45\fheight)},
scale only axis,
xmode=log,
xmin=0.1,
xmax=100,
xticklabels={{}},
xminorticks=true,
ymode=log,
ymin=0.0005,
ymax=20,
ytick={0.001, 0.01, 0.1, 1, 10},
yminorticks=true,
ylabel={Magnitude},
xminorgrids,
yminorgrids,
legend columns=2,
legend style={
/tikz/column 2/.style={
column sep=5pt,
},
at={(1,0)}, outer sep=2pt , anchor=south east, legend cell align=left, align=left, draw=black, nodes={scale=0.7, transform shape}
},
]
\addplot [color=mycolor1, dashed, line width=1.5pt]
table [x=freqs, y=W1, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_weights.csv};
\addlegendentry{${W_1}^{-1}$};
\addplot [color=mycolor1, line width=1.5pt]
table [x=freqs, y=H1, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\addlegendentry{$H_1$};
\addplot [color=mycolor2, dashed, line width=1.5pt]
table [x=freqs, y=W2, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_weights.csv};
\addlegendentry{${W_2}^{-1}$};
\addplot [color=mycolor2, line width=1.5pt]
table [x=freqs, y=H2, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\addlegendentry{$H_2$};
\addplot [color=mycolor3, dashed, line width=1.5pt]
table [x=freqs, y=W3, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_weights.csv};
\addlegendentry{${W_3}^{-1}$};
\addplot [color=mycolor3, line width=1.5pt]
table [x=freqs, y=H3, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\addlegendentry{$H_3$};
\end{axis}
\begin{axis}[%
width=1.0\fwidth,
height=0.4\fheight,
at={(0.0\fwidth, 0.0\fheight)},
scale only axis,
xmode=log,
xmin=0.1,
xmax=100,
xminorticks=true,
xlabel={Frequency [Hz]},
ymin=-240,
ymax=240,
ytick={-180, -90, 0, 90, 180},
ylabel={Phase [deg]},
xminorgrids,
]
\addplot [color=mycolor1, line width=1.5pt]
table [x=freqs, y=H1p, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\addplot [color=mycolor2, line width=1.5pt]
table [x=freqs, y=H2p, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\addplot [color=mycolor3, line width=1.5pt]
table [x=freqs, y=H3p, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/mathinf_three_results.csv};
\end{axis}
\end{tikzpicture}

BIN
tikz/figs/ligo_weights.pdf
View File

Binary file not shown.

BIN
tikz/figs/ligo_weights.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 21 KiB

373
tikz/figs/ligo_weights.svg
View File

File diff suppressed because one or more lines are too long
Side by Side

Before

Width:  |  Height:  |  Size: 94 KiB

66
tikz/figs/ligo_weights.tex
View File

@ -1,66 +0,0 @@
\setlength\fwidth{6.5cm}
\setlength\fheight{4cm}
\begin{tikzpicture}
\begin{axis}[%
width=1.0\fwidth,
height=1.0\fheight,
at={(0.0\fwidth, 0.0\fheight)},
scale only axis,
separate axis lines,
every outer x axis line/.append style={black},
every x tick label/.append style={font=\color{black}},
every x tick/.append style={black},
xmode=log,
xmin=0.001,
xmax=1,
xminorticks=true,
xlabel={Frequency [Hz]},
every outer y axis line/.append style={black},
every y tick label/.append style={font=\color{black}},
every y tick/.append style={black},
ymode=log,
ymin=0.005,
ymax=20,
yminorticks=true,
ylabel={Magnitude},
axis background/.style={fill=white},
xmajorgrids,
xminorgrids,
ymajorgrids,
yminorgrids,
legend style={at={(0,1)}, outer sep=2pt, anchor=north west, legend cell align=left, align=left, draw=black, nodes={scale=0.7, transform shape}}
]
\addplot [color=mycolor1, line width=1.5pt]
table [x=freqs, y=wHm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matligo_weights.csv};
\addlegendentry{$|w_H|^{-1}$}
\addplot [color=mycolor2, line width=1.5pt]
table [x=freqs, y=wLm, col sep=comma] {/home/thomas/Cloud/thesis/papers/dehaeze19_desig_compl_filte/matlab/matligo_weights.csv};
\addlegendentry{$|w_L|^{-1}$}
\addplot [color=black, dotted, line width=1.5pt]
table[row sep=crcr]{%
0.0005 0.008\\
0.008 0.008\\
};
\addlegendentry{Specifications}
\addplot [color=black, dotted, line width=1.5pt, forget plot]
table[row sep=crcr]{%
0.008 0.008\\
0.04 1\\
};
\addplot [color=black, dotted, line width=1.5pt, forget plot]
table[row sep=crcr]{%
0.04 3\\
0.1 3\\
};
\addplot [color=black, dotted, line width=1.5pt]
table[row sep=crcr]{%
0.1 0.045\\
2 0.045\\
};
\end{axis}
\end{tikzpicture}

2111
tikz/figs/three_complementary_filters_results.pdf Normal file
View File

File diff suppressed because it is too large Load Diff