Clean file, add comments, link to doc, ...

This commit is contained in:
Thomas Dehaeze 2020-11-10 11:16:16 +01:00
parent 093e4fd1e9
commit 73d6d4b132
60 changed files with 244 additions and 2042 deletions

49
css/custom.css Normal file
View File

@ -0,0 +1,49 @@
.figure p{
text-align: center;
}
.figure img{
max-width:100%;
display: block;
margin: auto;
}
table {
margin-left: auto;
margin-right: auto;
}
.org-src-container > pre.src:before {
display: inline;
position: absolute;
color: #808080;
background-color: white;
top: -10px;
left: 10px;
padding: 0px 4px;
border: 1px solid #d0d0d0;
font-size: 80%;
}
.org-src-container > pre {
margin-top: 1.5em;
position: relative;
overflow: visible;
}
.org-src-container > pre > code.src:before {
display: inline;
position: absolute;
color: #808080;
background-color: white;
top: -10px;
left: 10px;
padding: 0px 4px;
border: 1px solid #d0d0d0;
font-size: 80%;
}
.org-src-container > pre.src-emacs-lisp:before { content: 'Emacs Lisp'; }
.org-src-container > pre.src-elisp:before { content: 'Emacs Lisp'; }
.org-src-container > pre.src-sh:before { content: 'shell'; }
.org-src-container > pre.src-bash:before { content: 'bash'; }
.org-src-container > pre.src-org:before { content: 'Org mode'; }
.org-src-container > pre.src-python:before { content: 'Python'; }
.org-src-container > pre.src-matlab:before { content: 'Matlab'; }

View File

@ -513,7 +513,7 @@ legend{
padding:0;
white-space:normal}
.fa:before,#content .admonition-title:before,#content h1 .headerlink:before,#content h2 .headerlink:before,#content h3 .headerlink:before,#content h4 .headerlink:before,#content h5 .headerlink:before,#content h6 .headerlink:before,#content dl dt .headerlink:before,.icon:before,.wy-dropdown .caret:before,.wy-inline-validate.wy-inline-validate-success .wy-input-context:before,.wy-inline-validate.wy-inline-validate-danger .wy-input-context:before,.wy-inline-validate.wy-inline-validate-warning .wy-input-context:before,.wy-inline-validate.wy-inline-validate-info .wy-input-context:before,.wy-alert,#content .note,#content .attention,#content .caution,#content .danger,#content .error,#content .hint,#content .important,#content .tip,#content .warning,#content .seealso,#content .admonitiontodo,.btn,input[type="text"],input[type="password"],input[type="email"],input[type="url"],input[type="date"],input[type="month"],input[type="time"],input[type="datetime"],input[type="datetime-local"],input[type="week"],input[type="number"],input[type="search"],input[type="tel"],input[type="color"],select,textarea,#table-of-contents li.on a,#table-of-contents li.current>a,.wy-side-nav-search>a,.wy-side-nav-search .wy-dropdown>a,.wy-nav-top a{
.fa:before,#content .admonition-title:before,#content h1 .headerlink:before,#content h2 .headerlink:before,#content h3 .headerlink:before,#content h4 .headerlink:before,#content h5 .headerlink:before,#content h6 .headerlink:before,#content dl dt .headerlink:before,.icon:before,.wy-dropdown .caret:before,.wy-inline-validate.wy-inline-validate-success .wy-input-context:before,.wy-inline-validate.wy-inline-validate-danger .wy-input-context:before,.wy-inline-validate.wy-inline-validate-warning .wy-input-context:before,.wy-inline-validate.wy-inline-validate-info .wy-input-context:before,.wy-alert,#content .note,#content .attention,#content .caution,#content .danger,#content .error,#content .summary,#content .hint,#content .important,#content .tip,#content .warning,#content .question,#content .seealso,#content .admonitiontodo,.btn,input[type="text"],input[type="password"],input[type="email"],input[type="url"],input[type="date"],input[type="month"],input[type="time"],input[type="datetime"],input[type="datetime-local"],input[type="week"],input[type="number"],input[type="search"],input[type="tel"],input[type="color"],select,textarea,#table-of-contents li.on a,#table-of-contents li.current>a,.wy-side-nav-search>a,.wy-side-nav-search .wy-dropdown>a,.wy-nav-top a{
-webkit-font-smoothing:antialiased}
/*!
@ -576,7 +576,7 @@ a .fa,a #content .admonition-title,#content a .admonition-title{
.nav #content .admonition-title,#content .nav .admonition-title,.nav .icon{
display:inline}
.wy-alert,#content .note,#content .attention,#content .caution,#content .danger,#content .error,#content .hint,#content .important,#content .tip,#content .warning,#content .seealso,#content .admonitiontodo{
.wy-alert,#content .note,#content .attention,#content .caution,#content .danger,#content .error,#content .summary,#content .hint,#content .important,#content .tip,#content .warning,#content .question,#content .seealso,#content .admonitiontodo{
padding:12px;
line-height:24px;
margin-bottom:24px;
@ -596,32 +596,45 @@ a .fa,a #content .admonition-title,#content a .admonition-title{
#content .danger,#content .error{
background:#fdf3f2}
.wy-alert.wy-alert-warning,#content .wy-alert-warning.note,#content .attention,#content .caution,#content .wy-alert-warning.danger,#content .wy-alert-warning.error,#content .wy-alert-warning.hint,#content .wy-alert-warning.important,#content .wy-alert-warning.tip,#content .warning,#content .wy-alert-warning.seealso,#content .admonitiontodo{
.wy-alert.wy-alert-warning,#content .wy-alert-warning.note,#content .attention,#content .caution,#content .wy-alert-warning.danger,#content .wy-alert-warning.error,#content .wy-alert-warning.summary,#content .wy-alert-warning.hint,#content .wy-alert-warning.important,#content .wy-alert-warning.tip,#content .warning,#content .wy-alert-warning.seealso,#content .admonitiontodo{
background:#ffedcc}
#content .admonition-title.note:before, #content .admonition-title.seealso:before,
#content .admonition-title.warning:before, #content .admonition-title.caution:before,
#content .admonition-title.warning:before,
#content .admonition-title.caution:before,
#content .admonition-title.attention:before,
#content .admonition-title.tip:before, #content .admonition-title.hint:before,
#content .admonition-title.important:before,
#content .admonition-title.error:before, #content .admonition-title.danger:before{
#content .admonition-title.error:before,
#content .admonition-title.danger:before{
font-family:FontAwesome;
content: "";}
#content .note,#content .seealso{
#content .admonition-title.question:before{
font-family:FontAwesome;
content: "";}
#content .admonition-title.note:before,
#content .admonition-title.seealso:before,
#content .admonition-title.tip:before,
#content .admonition-title.summary:before,
#content .admonition-title.hint:before{
font-family:FontAwesome;
content: "";}
#content .note,#content .question,#content .seealso{
background:#e7f2fa}
.wy-alert p:last-child,#content .note p:last-child,#content .attention p:last-child,#content .caution p:last-child,#content .danger p:last-child,#content .error p:last-child,#content .hint p:last-child,#content .important p:last-child,#content .tip p:last-child,#content .warning p:last-child,#content .seealso p:last-child,#content .admonitiontodo p:last-child{
.wy-alert p:last-child,#content .note p:last-child,#content .attention p:last-child,#content .caution p:last-child,#content .danger p:last-child,#content .error p:last-child,#content .summary p:last-child,#content .hint p:last-child,#content .important p:last-child,#content .tip p:last-child,#content .warning p:last-child,#content .question p:last-child,#content .seealso p:last-child,#content .admonitiontodo p:last-child{
margin-bottom:0}
#content .admonition-title.tip,#content .admonition-title.important,#content .admonition-title.hint{
#content .admonition-title.tip,#content .admonition-title.important,#content .admonition-title.summary,#content .admonition-title.hint{
line-height: 1;
background:#1abc9c}
#content .important,#content .tip,#content .hint{
#content .important,#content .tip,#content .summary,#content .hint{
background:#dbfaf4}
#content .admonition-title.note,#content .admonition-title.seealso{
#content .admonition-title.note,#content .admonition-title.question,#content .admonition-title.seealso{
line-height: 1;
background:#6ab0de}
@ -938,7 +951,7 @@ footer p{
font-style:italic;
}
#content .note .last,#content .attention .last,#content .caution .last,#content .danger .last,#content .error .last,#content .hint .last,#content .important .last,#content .tip .last,#content .warning .last,#content .seealso .last,#content .admonitiontodo .last{
#content .note .last,#content .attention .last,#content .caution .last,#content .danger .last,#content .error .last,#content .hint .summary,#content .hint .last,#content .important .last,#content .tip .last,#content .warning .last,#content .question .last,#content .seealso .last,#content .admonitiontodo .last{
margin-bottom:0}
#content .admonition-title:before{

BIN
doc/APA95ML.pdf Normal file

Binary file not shown.

BIN
doc/IDS3010.pdf Normal file

Binary file not shown.

Binary file not shown.

Binary file not shown.

BIN
doc/LA75B.pdf Normal file

Binary file not shown.

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 34 KiB

Binary file not shown.

Binary file not shown.

After

Width:  |  Height:  |  Size: 52 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 37 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 37 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 39 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 41 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 39 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 40 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 41 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 132 KiB

After

Width:  |  Height:  |  Size: 136 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 50 KiB

After

Width:  |  Height:  |  Size: 57 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 35 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 48 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 25 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 36 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 47 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 27 KiB

Binary file not shown.

Binary file not shown.

Before

Width:  |  Height:  |  Size: 134 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 92 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 56 KiB

View File

@ -3,26 +3,17 @@
"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-11-04 mer. 20:38 -->
<!-- 2020-11-10 mar. 10:42 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<title>Encoder - Test Bench</title>
<meta name="generator" content="Org mode" />
<meta name="author" content="Dehaeze Thomas" />
<link rel="stylesheet" type="text/css" href="./css/htmlize.css"/>
<link rel="stylesheet" type="text/css" href="./css/readtheorg.css"/>
<link rel="stylesheet" type="text/css" href="./css/zenburn.css"/>
<link rel="stylesheet" type="text/css" href="./css/custom.css"/>
<script type="text/javascript" src="./js/jquery.min.js"></script>
<script type="text/javascript" src="./js/bootstrap.min.js"></script>
<script type="text/javascript" src="./js/jquery.stickytableheaders.min.js"></script>
<script type="text/javascript" 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">
@ -35,73 +26,69 @@
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org2b2f049">1. Experimental Setup</a></li>
<li><a href="#org9ae614e">2. Huddle Test</a>
<li><a href="#org8cea100">1. Experimental Setup</a></li>
<li><a href="#orgd6a888d">2. Noise Spectral Density of the Encoder</a>
<ul>
<li><a href="#orgc034802">2.1. Load Data</a></li>
<li><a href="#orgd8ae4ae">2.2. Time Domain Results</a></li>
<li><a href="#org864c60c">2.3. Frequency Domain Noise</a></li>
<li><a href="#org9e03e60">2.1. Load Data</a></li>
<li><a href="#orgb74e6ae">2.2. Time Domain Results</a></li>
<li><a href="#org9257580">2.3. Frequency Domain Noise</a></li>
</ul>
</li>
<li><a href="#org96539c5">3. Comparison Interferometer / Encoder</a>
<li><a href="#org09de461">3. Dynamics from Actuator to Encoder</a>
<ul>
<li><a href="#orgc6453b2">3.1. Load Data</a></li>
<li><a href="#orgc362641">3.2. Time Domain Results</a></li>
<li><a href="#orgb50503c">3.3. Difference between Encoder and Interferometer as a function of time</a></li>
<li><a href="#org5ec03aa">3.4. Difference between Encoder and Interferometer as a function of position</a></li>
</ul>
</li>
<li><a href="#org6b73303">4. Identification</a>
<ul>
<li><a href="#orga15e9e2">4.1. Load Data</a></li>
<li><a href="#org6974653">4.2. Identification</a></li>
</ul>
</li>
<li><a href="#org5be1a0f">5. Change of Stiffness due to Sensors stack being open/closed circuit</a>
<ul>
<li><a href="#org00fc76c">5.1. Load Data</a></li>
<li><a href="#org757917f">5.2. Transfer Functions</a></li>
</ul>
</li>
<li><a href="#org8dab16f">6. Generated Number of Charge / Voltage</a>
<ul>
<li><a href="#org28149b5">6.1. Steps</a></li>
<li><a href="#org900387e">6.2. Add Parallel Resistor</a></li>
<li><a href="#org7cfe5a8">6.3. Sinus</a></li>
<li><a href="#orgff48429">3.1. Load Data</a></li>
<li><a href="#org43af2dd">3.2. Identification</a></li>
</ul>
</li>
</ul>
</div>
</div>
<div id="outline-container-org2b2f049" class="outline-2">
<h2 id="org2b2f049"><span class="section-number-2">1</span> Experimental Setup</h2>
<p>
In this document, we wish to study the use of an encoder in parallel with an Amplified Piezoelectric Actuator.
</p>
<p>
The document is divided into the following Sections:
</p>
<ul class="org-ul">
<li>Section <a href="#org6bad4bb">1</a>: the test-bench used is described</li>
<li>Section <a href="#org4c3ff10">2</a>: the noise spectral density of the encoder is estimated</li>
<li>Section <a href="#org5214466">3</a>: the dynamics from the amplified piezoelectric actuator to the encoder measured displacement is identified</li>
</ul>
<div id="outline-container-org8cea100" class="outline-2">
<h2 id="org8cea100"><span class="section-number-2">1</span> Experimental Setup</h2>
<div class="outline-text-2" id="text-1">
<p>
The experimental Setup is schematically represented in Figure <a href="#orgff06236">1</a>.
<a id="org6bad4bb"></a>
</p>
<p>
The experimental Setup is schematically represented in Figure <a href="#org34db538">1</a>.
</p>
<p>
The mass can be vertically moved using the amplified piezoelectric actuator.
The displacement of the mass (relative to the mechanical frame) is measured both by the interferometer and by the encoder.
The displacement of the mass (relative to the mechanical frame) is measured both by the interferometer and by an encoder.
</p>
<div id="orgff06236" class="figure">
<div id="org34db538" class="figure">
<p><img src="figs/exp_setup_schematic.png" alt="exp_setup_schematic.png" />
</p>
<p><span class="figure-number">Figure 1: </span>Schematic of the Experiment</p>
</div>
<div id="orgc6b5ad5" class="figure">
<div id="org4486dea" class="figure">
<p><img src="figs/IMG_20201023_153905.jpg" alt="IMG_20201023_153905.jpg" />
</p>
<p><span class="figure-number">Figure 2: </span>Side View of the encoder</p>
</div>
<div id="orgcf538be" class="figure">
<div id="org8d81753" class="figure">
<p><img src="figs/IMG_20201023_153914.jpg" alt="IMG_20201023_153914.jpg" />
</p>
<p><span class="figure-number">Figure 3: </span>Front View of the encoder</p>
@ -109,10 +96,13 @@ The displacement of the mass (relative to the mechanical frame) is measured both
</div>
</div>
<div id="outline-container-org9ae614e" class="outline-2">
<h2 id="org9ae614e"><span class="section-number-2">2</span> Huddle Test</h2>
<div id="outline-container-orgd6a888d" class="outline-2">
<h2 id="orgd6a888d"><span class="section-number-2">2</span> Noise Spectral Density of the Encoder</h2>
<div class="outline-text-2" id="text-2">
<p>
<a id="org4c3ff10"></a>
</p>
<p>
The goal in this section is the estimate the noise of both the encoder and the intereferometer.
</p>
@ -122,8 +112,8 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
</p>
</div>
<div id="outline-container-orgc034802" class="outline-3">
<h3 id="orgc034802"><span class="section-number-3">2.1</span> Load Data</h3>
<div id="outline-container-org9e03e60" class="outline-3">
<h3 id="org9e03e60"><span class="section-number-3">2.1</span> Load Data</h3>
<div class="outline-text-3" id="text-2-1">
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'mat/int_enc_huddle_test.mat'</span>, <span class="org-string">'interferometer'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'t'</span>);
@ -138,11 +128,11 @@ encoder = detrend(encoder, 0);
</div>
</div>
<div id="outline-container-orgd8ae4ae" class="outline-3">
<h3 id="orgd8ae4ae"><span class="section-number-3">2.2</span> Time Domain Results</h3>
<div id="outline-container-orgb74e6ae" class="outline-3">
<h3 id="orgb74e6ae"><span class="section-number-3">2.2</span> Time Domain Results</h3>
<div class="outline-text-3" id="text-2-2">
<div id="orgff47510" class="figure">
<div id="org170b58a" class="figure">
<p><img src="figs/huddle_test_time_domain.png" alt="huddle_test_time_domain.png" />
</p>
<p><span class="figure-number">Figure 4: </span>Huddle test - Time domain signals</p>
@ -154,7 +144,7 @@ encoder = detrend(encoder, 0);
</div>
<div id="org21436f3" class="figure">
<div id="orgca372e2" class="figure">
<p><img src="figs/huddle_test_time_domain_filtered.png" alt="huddle_test_time_domain_filtered.png" />
</p>
<p><span class="figure-number">Figure 5: </span>Huddle test - Time domain signals filtered with a LPF at 10Hz</p>
@ -162,8 +152,8 @@ encoder = detrend(encoder, 0);
</div>
</div>
<div id="outline-container-org864c60c" class="outline-3">
<h3 id="org864c60c"><span class="section-number-3">2.3</span> Frequency Domain Noise</h3>
<div id="outline-container-org9257580" class="outline-3">
<h3 id="org9257580"><span class="section-number-3">2.3</span> Frequency Domain Noise</h3>
<div class="outline-text-3" id="text-2-3">
<div class="org-src-container">
<pre class="src src-matlab">Ts = 1e<span class="org-type">-</span>4;
@ -175,7 +165,7 @@ win = hann(ceil(10<span class="org-type">/</span>Ts));
</div>
<div id="org63bcf43" class="figure">
<div id="org58bfefc" class="figure">
<p><img src="figs/huddle_test_asd.png" alt="huddle_test_asd.png" />
</p>
<p><span class="figure-number">Figure 6: </span>Amplitude Spectral Density of the signals during the Huddle test</p>
@ -184,139 +174,18 @@ win = hann(ceil(10<span class="org-type">/</span>Ts));
</div>
</div>
<div id="outline-container-org96539c5" class="outline-2">
<h2 id="org96539c5"><span class="section-number-2">3</span> Comparison Interferometer / Encoder</h2>
<div id="outline-container-org09de461" class="outline-2">
<h2 id="org09de461"><span class="section-number-2">3</span> Dynamics from Actuator to Encoder</h2>
<div class="outline-text-2" id="text-3">
<p>
The goal here is to make sure that the interferometer and encoder measurements are coherent.
We may see non-linearity in the interferometric measurement.
<a id="org5214466"></a>
</p>
</div>
<div id="outline-container-orgc6453b2" class="outline-3">
<h3 id="orgc6453b2"><span class="section-number-3">3.1</span> Load Data</h3>
<div id="outline-container-orgff48429" class="outline-3">
<h3 id="orgff48429"><span class="section-number-3">3.1</span> Load Data</h3>
<div class="outline-text-3" id="text-3-1">
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'mat/int_enc_comp.mat'</span>, <span class="org-string">'interferometer'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>, <span class="org-string">'t'</span>);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab">interferometer = detrend(interferometer, 0);
encoder = detrend(encoder, 0);
u = detrend(u, 0);
</pre>
</div>
</div>
</div>
<div id="outline-container-orgc362641" class="outline-3">
<h3 id="orgc362641"><span class="section-number-3">3.2</span> Time Domain Results</h3>
<div class="outline-text-3" id="text-3-2">
<div id="orgfd33423" class="figure">
<p><img src="figs/int_enc_one_cycle.png" alt="int_enc_one_cycle.png" />
</p>
<p><span class="figure-number">Figure 7: </span>One cycle measurement</p>
</div>
<div id="orgf82fd3e" class="figure">
<p><img src="figs/int_enc_one_cycle_error.png" alt="int_enc_one_cycle_error.png" />
</p>
<p><span class="figure-number">Figure 8: </span>Difference between the Encoder and the interferometer during one cycle</p>
</div>
</div>
</div>
<div id="outline-container-orgb50503c" class="outline-3">
<h3 id="orgb50503c"><span class="section-number-3">3.3</span> Difference between Encoder and Interferometer as a function of time</h3>
<div class="outline-text-3" id="text-3-3">
<div class="org-src-container">
<pre class="src src-matlab">Ts = 1e<span class="org-type">-</span>4;
d_i_mean = reshape(interferometer, [2<span class="org-type">/</span>Ts floor(Ts<span class="org-type">/</span>2<span class="org-type">*</span>length(interferometer))]);
d_e_mean = reshape(encoder, [2<span class="org-type">/</span>Ts floor(Ts<span class="org-type">/</span>2<span class="org-type">*</span>length(encoder))]);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab">w0 = 2<span class="org-type">*</span><span class="org-constant">pi</span><span class="org-type">*</span>5; <span class="org-comment">% [rad/s]</span>
xi = 0.7;
G_lpf = 1<span class="org-type">/</span>(1 <span class="org-type">+</span> 2<span class="org-type">*</span>xi<span class="org-type">/</span>w0<span class="org-type">*</span>s <span class="org-type">+</span> s<span class="org-type">^</span>2<span class="org-type">/</span>w0<span class="org-type">^</span>2);
d_err_mean = reshape(lsim(G_lpf, encoder <span class="org-type">-</span> interferometer, t), [2<span class="org-type">/</span>Ts floor(Ts<span class="org-type">/</span>2<span class="org-type">*</span>length(encoder))]);
d_err_mean = d_err_mean <span class="org-type">-</span> mean(d_err_mean);
</pre>
</div>
<div id="org6eb78a0" class="figure">
<p><img src="figs/int_enc_error_mean_time.png" alt="int_enc_error_mean_time.png" />
</p>
<p><span class="figure-number">Figure 9: </span>Difference between the two measurement in the time domain, averaged for all the cycles</p>
</div>
</div>
</div>
<div id="outline-container-org5ec03aa" class="outline-3">
<h3 id="org5ec03aa"><span class="section-number-3">3.4</span> Difference between Encoder and Interferometer as a function of position</h3>
<div class="outline-text-3" id="text-3-4">
<p>
Compute the mean of the interferometer measurement corresponding to each of the encoder measurement.
</p>
<div class="org-src-container">
<pre class="src src-matlab">[e_sorted, <span class="org-type">~</span>, e_ind] = unique(encoder);
i_mean = zeros(length(e_sorted), 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(e_sorted)</span>
i_mean(<span class="org-constant">i</span>) = mean(interferometer(e_ind <span class="org-type">==</span> <span class="org-constant">i</span>));
<span class="org-keyword">end</span>
i_mean_error = (i_mean <span class="org-type">-</span> e_sorted);
</pre>
</div>
<div id="org2aecb58" class="figure">
<p><img src="figs/int_enc_error_mean_position.png" alt="int_enc_error_mean_position.png" />
</p>
<p><span class="figure-number">Figure 10: </span>Difference between the two measurement as a function of the measured position by the encoder, averaged for all the cycles</p>
</div>
<p>
The period of the non-linearity seems to be \(1.53 \mu m\) which corresponds to the wavelength of the Laser.
</p>
<div class="org-src-container">
<pre class="src src-matlab">win_length = 1530; <span class="org-comment">% length of the windows (corresponds to 1.53 um)</span>
num_avg = floor(length(e_sorted)<span class="org-type">/</span>win_length); <span class="org-comment">% number of averaging</span>
i_init = ceil((length(e_sorted) <span class="org-type">-</span> win_length<span class="org-type">*</span>num_avg)<span class="org-type">/</span>2); <span class="org-comment">% does not start at the extremity</span>
e_sorted_mean_over_period = mean(reshape(i_mean_error(i_init<span class="org-type">:</span>i_init<span class="org-type">+</span>win_length<span class="org-type">*</span>num_avg<span class="org-type">-</span>1), [win_length num_avg]), 2);
</pre>
</div>
<div id="orgd0cea4e" class="figure">
<p><img src="figs/int_non_linearity_period_wavelength.png" alt="int_non_linearity_period_wavelength.png" />
</p>
<p><span class="figure-number">Figure 11: </span>Non-Linearity of the Interferometer over the period of the wavelength</p>
</div>
</div>
</div>
</div>
<div id="outline-container-org6b73303" class="outline-2">
<h2 id="org6b73303"><span class="section-number-2">4</span> Identification</h2>
<div class="outline-text-2" id="text-4">
</div>
<div id="outline-container-orga15e9e2" class="outline-3">
<h3 id="orga15e9e2"><span class="section-number-3">4.1</span> Load Data</h3>
<div class="outline-text-3" id="text-4-1">
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'mat/int_enc_id_noise_bis.mat'</span>, <span class="org-string">'interferometer'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>, <span class="org-string">'t'</span>);
</pre>
</div>
@ -330,9 +199,9 @@ u = detrend(u, 0);
</div>
</div>
<div id="outline-container-org6974653" class="outline-3">
<h3 id="org6974653"><span class="section-number-3">4.2</span> Identification</h3>
<div class="outline-text-3" id="text-4-2">
<div id="outline-container-org43af2dd" class="outline-3">
<h3 id="org43af2dd"><span class="section-number-3">3.2</span> Identification</h3>
<div class="outline-text-3" id="text-3-2">
<div class="org-src-container">
<pre class="src src-matlab">Ts = 1e<span class="org-type">-</span>4; <span class="org-comment">% Sampling Time [s]</span>
win = hann(ceil(10<span class="org-type">/</span>Ts));
@ -349,579 +218,24 @@ win = hann(ceil(10<span class="org-type">/</span>Ts));
</div>
<div id="org37fcd05" class="figure">
<div id="orgb8b5d9d" class="figure">
<p><img src="figs/identification_dynamics_coherence.png" alt="identification_dynamics_coherence.png" />
</p>
</div>
<div id="orge65c21c" class="figure">
<div id="orgb95b979" class="figure">
<p><img src="figs/identification_dynamics_bode.png" alt="identification_dynamics_bode.png" />
</p>
</div>
</div>
</div>
</div>
<div id="outline-container-org5be1a0f" class="outline-2">
<h2 id="org5be1a0f"><span class="section-number-2">5</span> Change of Stiffness due to Sensors stack being open/closed circuit</h2>
<div class="outline-text-2" id="text-5">
</div>
<div id="outline-container-org00fc76c" class="outline-3">
<h3 id="org00fc76c"><span class="section-number-3">5.1</span> Load Data</h3>
<div class="outline-text-3" id="text-5-1">
<div class="org-src-container">
<pre class="src src-matlab">oc = load(<span class="org-string">'./mat/identification_open_circuit.mat'</span>, <span class="org-string">'t'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>);
sc = load(<span class="org-string">'./mat/identification_short_circuit.mat'</span>, <span class="org-string">'t'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>);
</pre>
</div>
</div>
</div>
<div id="outline-container-org757917f" class="outline-3">
<h3 id="org757917f"><span class="section-number-3">5.2</span> Transfer Functions</h3>
<div class="outline-text-3" id="text-5-2">
<div class="org-src-container">
<pre class="src src-matlab">Ts = 1e<span class="org-type">-</span>4; <span class="org-comment">% Sampling Time [s]</span>
win = hann(ceil(10<span class="org-type">/</span>Ts));
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab">[tf_oc_est, f] = tfestimate(oc.u, oc.encoder, win, [], [], 1<span class="org-type">/</span>Ts);
[co_oc_est, <span class="org-type">~</span>] = mscohere( oc.u, oc.encoder, win, [], [], 1<span class="org-type">/</span>Ts);
[tf_sc_est, <span class="org-type">~</span>] = tfestimate(sc.u, sc.encoder, win, [], [], 1<span class="org-type">/</span>Ts);
[co_sc_est, <span class="org-type">~</span>] = mscohere( sc.u, sc.encoder, win, [], [], 1<span class="org-type">/</span>Ts);
</pre>
</div>
<div id="orgcba7c3a" class="figure">
<p><img src="figs/stiffness_force_sensor_coherence.png" alt="stiffness_force_sensor_coherence.png" />
</p>
</div>
<div id="orgfdef3b4" class="figure">
<p><img src="figs/stiffness_force_sensor_bode.png" alt="stiffness_force_sensor_bode.png" />
</p>
</div>
<div id="org655038d" class="figure">
<p><img src="figs/stiffness_force_sensor_bode_zoom.png" alt="stiffness_force_sensor_bode_zoom.png" />
</p>
<p><span class="figure-number">Figure 16: </span>Zoom on the change of resonance</p>
</div>
<div class="important" id="orga5e301b">
<p>
The change of resonance frequency / stiffness is very small and is not important here.
</p>
</div>
</div>
</div>
</div>
<div id="outline-container-org8dab16f" class="outline-2">
<h2 id="org8dab16f"><span class="section-number-2">6</span> Generated Number of Charge / Voltage</h2>
<div class="outline-text-2" id="text-6">
<p>
Two stacks are used as actuator (in parallel) and one stack is used as sensor.
</p>
<p>
The amplifier gain is 20V/V (Cedrat LA75B).
</p>
</div>
<div id="outline-container-org28149b5" class="outline-3">
<h3 id="org28149b5"><span class="section-number-3">6.1</span> Steps</h3>
<div class="outline-text-3" id="text-6-1">
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'./mat/force_sensor_steps.mat'</span>, <span class="org-string">'t'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>, <span class="org-string">'v'</span>);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-type">figure</span>;
tiledlayout(2, 1, <span class="org-string">'TileSpacing'</span>, <span class="org-string">'None'</span>, <span class="org-string">'Padding'</span>, <span class="org-string">'None'</span>);
nexttile;
plot(t, v);
xlabel(<span class="org-string">'Time [s]'</span>); ylabel(<span class="org-string">'Measured voltage [V]'</span>);
nexttile;
plot(t, u);
xlabel(<span class="org-string">'Time [s]'</span>); ylabel(<span class="org-string">'Actuator Voltage [V]'</span>);
</pre>
</div>
<div id="org0c7c950" class="figure">
<p><img src="figs/force_sen_steps_time_domain.png" alt="force_sen_steps_time_domain.png" />
</p>
<p><span class="figure-number">Figure 17: </span>Time domain signal during the 3 actuator voltage steps</p>
</div>
<p>
Three steps are performed at the following time intervals:
</p>
<div class="org-src-container">
<pre class="src src-matlab">t_s = [ 2.5, 23;
23.8, 35;
35.8, 50];
</pre>
</div>
<p>
Fit function:
</p>
<div class="org-src-container">
<pre class="src src-matlab">f = @(b,x) b(1)<span class="org-type">.*</span>exp(b(2)<span class="org-type">.*</span>x) <span class="org-type">+</span> b(3);
</pre>
</div>
<p>
We are interested by the <code>b(2)</code> term, which is the time constant of the exponential.
</p>
<div class="org-src-container">
<pre class="src src-matlab">tau = zeros(size(t_s, 1),1);
V0 = zeros(size(t_s, 1),1);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-keyword">for</span> <span class="org-variable-name">t_i</span> = <span class="org-constant">1:size(t_s, 1)</span>
t_cur = t(t_s(t_i, 1) <span class="org-type">&lt;</span> t <span class="org-type">&amp;</span> t <span class="org-type">&lt;</span> t_s(t_i, 2));
t_cur = t_cur <span class="org-type">-</span> t_cur(1);
y_cur = v(t_s(t_i, 1) <span class="org-type">&lt;</span> t <span class="org-type">&amp;</span> t <span class="org-type">&lt;</span> t_s(t_i, 2));
nrmrsd = @(b) norm(y_cur <span class="org-type">-</span> f(b,t_cur)); <span class="org-comment">% Residual Norm Cost Function</span>
B0 = [0.5, <span class="org-type">-</span>0.15, 2.2]; <span class="org-comment">% Choose Appropriate Initial Estimates</span>
[B,rnrm] = fminsearch(nrmrsd, B0); <span class="org-comment">% Estimate Parameters &#8216;B&#8217;</span>
tau(t_i) = 1<span class="org-type">/</span>B(2);
V0(t_i) = B(3);
<span class="org-keyword">end</span>
</pre>
</div>
<table border="2" cellspacing="0" cellpadding="6" rules="groups" frame="hsides">
<colgroup>
<col class="org-right" />
<col class="org-right" />
</colgroup>
<thead>
<tr>
<th scope="col" class="org-right">\(tau\) [s]</th>
<th scope="col" class="org-right">\(V_0\) [V]</th>
</tr>
</thead>
<tbody>
<tr>
<td class="org-right">6.47</td>
<td class="org-right">2.26</td>
</tr>
<tr>
<td class="org-right">6.76</td>
<td class="org-right">2.26</td>
</tr>
<tr>
<td class="org-right">6.49</td>
<td class="org-right">2.25</td>
</tr>
</tbody>
</table>
<p>
With the capacitance being \(C = 4.4 \mu F\), the internal impedance of the Speedgoat ADC can be computed as follows:
</p>
<div class="org-src-container">
<pre class="src src-matlab">Cp = 4.4e<span class="org-type">-</span>6; <span class="org-comment">% [F]</span>
Rin = abs(mean(tau))<span class="org-type">/</span>Cp;
</pre>
</div>
<pre class="example">
1494100.0
</pre>
<p>
The input impedance of the Speedgoat&rsquo;s ADC should then be close to \(1.5\,M\Omega\) (specified at \(1\,M\Omega\)).
</p>
<div class="important" id="org60ccf75">
<p>
How can we explain the voltage offset?
</p>
</div>
<p>
As shown in Figure <a href="#org39e1694">18</a> (taken from (<a href="#citeproc_bib_item_1">Reza and Andrew 2006</a>)), an input voltage offset is due to the input bias current \(i_n\).
</p>
<div id="org39e1694" class="figure">
<p><img src="figs/piezo_sensor_model_instrumentation.png" alt="piezo_sensor_model_instrumentation.png" />
</p>
<p><span class="figure-number">Figure 18: </span>Model of a piezoelectric transducer (left) and instrumentation amplifier (right)</p>
</div>
<p>
The estimated input bias current is then:
</p>
<div class="org-src-container">
<pre class="src src-matlab">in = mean(V0)<span class="org-type">/</span>Rin;
</pre>
</div>
<pre class="example">
1.5119e-06
</pre>
<p>
An additional resistor in parallel with \(R_{in}\) would have two effects:
</p>
<ul class="org-ul">
<li>reduce the input voltage offset
\[ V_{off} = \frac{R_a R_{in}}{R_a + R_{in}} i_n \]</li>
<li>increase the high pass corner frequency \(f_c\)
\[ C_p \frac{R_{in}R_a}{R_{in} + R_a} = \tau_c = \frac{1}{f_c} \]
\[ R_a = \frac{R_i}{f_c C_p R_i - 1} \]</li>
</ul>
<p>
If we allow the high pass corner frequency to be equals to 3Hz:
</p>
<div class="org-src-container">
<pre class="src src-matlab">fc = 3;
Ra = Rin<span class="org-type">/</span>(fc<span class="org-type">*</span>Cp<span class="org-type">*</span>Rin <span class="org-type">-</span> 1);
</pre>
</div>
<pre class="example">
79804
</pre>
<p>
With this parallel resistance value, the voltage offset would be:
</p>
<div class="org-src-container">
<pre class="src src-matlab">V_offset = Ra<span class="org-type">*</span>Rin<span class="org-type">/</span>(Ra <span class="org-type">+</span> Rin) <span class="org-type">*</span> in;
</pre>
</div>
<pre class="example">
0.11454
</pre>
<p>
Which is much more acceptable.
</p>
</div>
</div>
<div id="outline-container-org900387e" class="outline-3">
<h3 id="org900387e"><span class="section-number-3">6.2</span> Add Parallel Resistor</h3>
<div class="outline-text-3" id="text-6-2">
<p>
A resistor of \(\approx 100\,k\Omega\) is added in parallel with the force sensor and the same kin.
</p>
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'./mat/force_sensor_steps_R_82k7.mat'</span>, <span class="org-string">'t'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>, <span class="org-string">'v'</span>);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-type">figure</span>;
tiledlayout(2, 1, <span class="org-string">'TileSpacing'</span>, <span class="org-string">'None'</span>, <span class="org-string">'Padding'</span>, <span class="org-string">'None'</span>);
nexttile;
plot(t, v);
xlabel(<span class="org-string">'Time [s]'</span>); ylabel(<span class="org-string">'Measured voltage [V]'</span>);
nexttile;
plot(t, u);
xlabel(<span class="org-string">'Time [s]'</span>); ylabel(<span class="org-string">'Actuator Voltage [V]'</span>);
</pre>
</div>
<div id="orgaee44e2" class="figure">
<p><img src="figs/force_sen_steps_time_domain_par_R.png" alt="force_sen_steps_time_domain_par_R.png" />
</p>
<p><span class="figure-number">Figure 19: </span>Time domain signal during the actuator voltage steps</p>
</div>
<p>
Three steps are performed at the following time intervals:
</p>
<div class="org-src-container">
<pre class="src src-matlab">t_s = [1.9, 6;
8.5, 13;
15.5, 21;
22.6, 26;
30.0, 36;
37.5, 41;
46.2, 49.5]
</pre>
</div>
<p>
Fit function:
</p>
<div class="org-src-container">
<pre class="src src-matlab">f = @(b,x) b(1)<span class="org-type">.*</span>exp(b(2)<span class="org-type">.*</span>x) <span class="org-type">+</span> b(3);
</pre>
</div>
<p>
We are interested by the <code>b(2)</code> term, which is the time constant of the exponential.
</p>
<div class="org-src-container">
<pre class="src src-matlab">tau = zeros(size(t_s, 1),1);
V0 = zeros(size(t_s, 1),1);
</pre>
</div>
<div class="org-src-container">
<pre class="src src-matlab"><span class="org-keyword">for</span> <span class="org-variable-name">t_i</span> = <span class="org-constant">1:size(t_s, 1)</span>
t_cur = t(t_s(t_i, 1) <span class="org-type">&lt;</span> t <span class="org-type">&amp;</span> t <span class="org-type">&lt;</span> t_s(t_i, 2));
t_cur = t_cur <span class="org-type">-</span> t_cur(1);
y_cur = v(t_s(t_i, 1) <span class="org-type">&lt;</span> t <span class="org-type">&amp;</span> t <span class="org-type">&lt;</span> t_s(t_i, 2));
nrmrsd = @(b) norm(y_cur <span class="org-type">-</span> f(b,t_cur)); <span class="org-comment">% Residual Norm Cost Function</span>
B0 = [0.5, <span class="org-type">-</span>0.2, 0.2]; <span class="org-comment">% Choose Appropriate Initial Estimates</span>
[B,rnrm] = fminsearch(nrmrsd, B0); <span class="org-comment">% Estimate Parameters &#8216;B&#8217;</span>
tau(t_i) = 1<span class="org-type">/</span>B(2);
V0(t_i) = B(3);
<span class="org-keyword">end</span>
</pre>
</div>
<p>
And indeed, we obtain a much smaller offset voltage and a much faster time constant.
</p>
<table border="2" cellspacing="0" cellpadding="6" rules="groups" frame="hsides">
<colgroup>
<col class="org-right" />
<col class="org-right" />
</colgroup>
<thead>
<tr>
<th scope="col" class="org-right">\(tau\) [s]</th>
<th scope="col" class="org-right">\(V_0\) [V]</th>
</tr>
</thead>
<tbody>
<tr>
<td class="org-right">0.43</td>
<td class="org-right">0.15</td>
</tr>
<tr>
<td class="org-right">0.45</td>
<td class="org-right">0.16</td>
</tr>
<tr>
<td class="org-right">0.43</td>
<td class="org-right">0.15</td>
</tr>
<tr>
<td class="org-right">0.43</td>
<td class="org-right">0.15</td>
</tr>
<tr>
<td class="org-right">0.45</td>
<td class="org-right">0.15</td>
</tr>
<tr>
<td class="org-right">0.46</td>
<td class="org-right">0.16</td>
</tr>
<tr>
<td class="org-right">0.48</td>
<td class="org-right">0.16</td>
</tr>
</tbody>
</table>
<p>
Knowing the capacitance value, we can estimate the value of the added resistor (neglecting the input impedance of \(\approx 1\,M\Omega\)):
</p>
<div class="org-src-container">
<pre class="src src-matlab">Cp = 4.4e<span class="org-type">-</span>6; <span class="org-comment">% [F]</span>
Rin = abs(mean(tau))<span class="org-type">/</span>Cp;
</pre>
</div>
<pre class="example">
101200.0
</pre>
<p>
And we can verify that the bias current estimation stays the same:
</p>
<div class="org-src-container">
<pre class="src src-matlab">in = mean(V0)<span class="org-type">/</span>Rin;
</pre>
</div>
<pre class="example">
1.5305e-06
</pre>
<p>
This validates the model of the ADC and the effectiveness of the added resistor.
</p>
</div>
</div>
<div id="outline-container-org7cfe5a8" class="outline-3">
<h3 id="org7cfe5a8"><span class="section-number-3">6.3</span> Sinus</h3>
<div class="outline-text-3" id="text-6-3">
<div class="org-src-container">
<pre class="src src-matlab">load(<span class="org-string">'./mat/force_sensor_sin.mat'</span>, <span class="org-string">'t'</span>, <span class="org-string">'encoder'</span>, <span class="org-string">'u'</span>, <span class="org-string">'v'</span>);
u = u(t<span class="org-type">&gt;</span>25);
v = v(t<span class="org-type">&gt;</span>25);
encoder = encoder(t<span class="org-type">&gt;</span>25) <span class="org-type">-</span> mean(encoder(t<span class="org-type">&gt;</span>25));
t = t(t<span class="org-type">&gt;</span>25);
</pre>
</div>
<p>
The driving voltage is a sinus at 0.5Hz centered on 3V and with an amplitude of 3V (Figure <a href="#orga42b4f8">20</a>).
</p>
<div id="orga42b4f8" class="figure">
<p><img src="figs/force_sensor_sin_u.png" alt="force_sensor_sin_u.png" />
</p>
<p><span class="figure-number">Figure 20: </span>Driving Voltage</p>
</div>
<p>
The full stroke as measured by the encoder is:
</p>
<div class="org-src-container">
<pre class="src src-matlab">max(encoder)<span class="org-type">-</span>min(encoder)
</pre>
</div>
<pre class="example">
5.005e-05
</pre>
<p>
Its signal is shown in Figure <a href="#orgcda82ea">21</a>.
</p>
<div id="orgcda82ea" class="figure">
<p><img src="figs/force_sensor_sin_encoder.png" alt="force_sensor_sin_encoder.png" />
</p>
<p><span class="figure-number">Figure 21: </span>Encoder measurement</p>
</div>
<p>
The generated voltage by the stack is shown in Figure
</p>
<div id="orgdd5fd9b" class="figure">
<p><img src="figs/force_sensor_sin_stack.png" alt="force_sensor_sin_stack.png" />
</p>
<p><span class="figure-number">Figure 22: </span>Voltage measured on the stack used as a sensor</p>
</div>
<p>
The capacitance of the stack is
</p>
<div class="org-src-container">
<pre class="src src-matlab">Cp = 4.4e<span class="org-type">-</span>6; <span class="org-comment">% [F]</span>
</pre>
</div>
<p>
The corresponding generated charge is then shown in Figure <a href="#org15cf433">23</a>.
</p>
<div id="org15cf433" class="figure">
<p><img src="figs/force_sensor_sin_charge.png" alt="force_sensor_sin_charge.png" />
</p>
<p><span class="figure-number">Figure 23: </span>Generated Charge</p>
</div>
<p>
The relation between the generated voltage and the measured displacement is almost linear as shown in Figure <a href="#orge276186">24</a>.
</p>
<div class="org-src-container">
<pre class="src src-matlab">b1 = encoder<span class="org-type">\</span>(v<span class="org-type">-</span>mean(v));
</pre>
</div>
<div id="orge276186" class="figure">
<p><img src="figs/force_sensor_linear_relation.png" alt="force_sensor_linear_relation.png" />
</p>
<p><span class="figure-number">Figure 24: </span>Almost linear relation between the relative displacement and the generated voltage</p>
</div>
<p>
With a 16bits ADC, the resolution will then be equals to (in [nm]):
</p>
<div class="org-src-container">
<pre class="src src-matlab">abs((20<span class="org-type">/</span>2<span class="org-type">^</span>16)<span class="org-type">/</span>(b1<span class="org-type">/</span>1e9))
</pre>
</div>
<pre class="example">
3.9838
</pre>
<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>Reza, Moheimani, and Fleming Andrew. 2006. <i>Piezoelectric Transducers for Vibration Control and Damping</i>. London: Springer.</div>
</div>
</div>
</div>
</div>
</div>
<div id="postamble" class="status">
<p class="author">Author: Dehaeze Thomas</p>
<p class="date">Created: 2020-11-04 mer. 20:38</p>
<p class="date">Created: 2020-11-10 mar. 10:42</p>
</div>
</body>
</html>

775
index.org
View File

@ -9,10 +9,9 @@
#+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: <link rel="stylesheet" type="text/css" href="./css/zenburn.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/custom.css"/>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/readtheorg.js"></script>
#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
@ -27,7 +26,7 @@
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
#+PROPERTY: header-args:matlab :session *MATLAB*
#+PROPERTY: header-args:matlab+ :tangle script.m
#+PROPERTY: header-args:matlab+ :tangle matlab/script.m
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :results none
@ -37,9 +36,28 @@
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:
* Introduction :ignore:
In this document, we wish to study the use of an encoder in parallel with an Amplified Piezoelectric Actuator.
The document is divided into the following Sections:
- Section [[sec:experimental_setup]]: the test-bench used is described
- Section [[sec:encoder_noise]]: the noise spectral density of the encoder is estimated
- Section [[sec:dynamics_encoder]]: the dynamics from the amplified piezoelectric actuator to the encoder measured displacement is identified
* Experimental Setup
<<sec:experimental_setup>>
The experimental Setup is schematically represented in Figure [[fig:exp_setup_schematic]].
#+begin_note
Here are the equipment used in the test bench:
- Renishaw Resolution Encoder with 1nm resolution ([[file:doc/L-9517-9448-05-B_Data_sheet_RESOLUTE_BiSS_en.pdf][doc]])
- Attocube interferometer ([[file:doc/IDS3010.pdf][doc]])
- Cedrat Amplified Piezoelectric Actuator APA95ML ([[file:doc/APA95ML.pdf][doc]])
- Voltage Amplifier LA75B ([[file:doc/LA75B.pdf][doc]])
- Speedgoat IO131 with 16bits ADC and DAC ([[file:doc/IO130 IO131 OEM Datasheet.pdf][doc]])
#+end_note
The mass can be vertically moved using the amplified piezoelectric actuator.
The displacement of the mass (relative to the mechanical frame) is measured both by the interferometer and by the encoder.
@ -48,6 +66,7 @@ The displacement of the mass (relative to the mechanical frame) is measured both
[[file:figs/exp_setup_schematic.png]]
#+name: fig:encoder_side_view
#+ATTR_ORG: :width 300
#+caption: Side View of the encoder
[[file:figs/IMG_20201023_153905.jpg]]
@ -55,11 +74,13 @@ The displacement of the mass (relative to the mechanical frame) is measured both
#+caption: Front View of the encoder
[[file:figs/IMG_20201023_153914.jpg]]
* Huddle Test
* Noise Spectral Density of the Encoder
<<sec:encoder_noise>>
** Introduction :ignore:
The goal in this section is the estimate the noise of both the encoder and the intereferometer.
Nothing is then to the actuator such that the relative motion between the mass and the frame is as small as possible.
The actuator is not excited, thus the relative motion between the mass and the frame is as small as possible.
Ideally, a mechanical part would clamp the two together, we here suppose that the APA is still enough to clamp the two together.
** Matlab Init :noexport:ignore:
@ -72,6 +93,8 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
#+end_src
** Load Data
The measurement data are loaded and the offset are removed using the =detrend= command.
#+begin_src matlab
load('mat/int_enc_huddle_test.mat', 'interferometer', 'encoder', 't');
#+end_src
@ -82,6 +105,8 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
#+end_src
** Time Domain Results
The measurement of both the encoder and interferometer are shown in Figure [[fig:huddle_test_time_domain]].
#+begin_src matlab :exports none
figure;
hold on;
@ -101,6 +126,7 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
#+RESULTS:
[[file:figs/huddle_test_time_domain.png]]
The raw signals are filtered with a Low Pass filter (defined below) such that we can see the low frequency motion (Figure [[fig:huddle_test_time_domain_filtered]]).
#+begin_src matlab
G_lpf = 1/(1 + s/2/pi/10);
#+end_src
@ -125,6 +151,8 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
[[file:figs/huddle_test_time_domain_filtered.png]]
** Frequency Domain Noise
The noise of the measurement (supposing there is no motion) is now translated in the frequency domain by computed the Amplitude Spectral Density.
#+begin_src matlab
Ts = 1e-4;
win = hann(ceil(10/Ts));
@ -133,6 +161,10 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
[p_e, ~] = pwelch(encoder, win, [], [], 1/Ts);
#+end_src
The comparison of the ASD of the encoder and interferometer are shown in Figure [[fig:huddle_test_asd]].
It is clear that although the encoder exhibit higher frequency noise, is it more stable at low frequency as the length of the beam path in the air is much smaller and thus changed of temperature/pressure/humity of the air has much smaller effect on the measured displacement.
#+begin_src matlab :exports none
figure;
hold on;
@ -154,10 +186,11 @@ Ideally, a mechanical part would clamp the two together, we here suppose that th
#+RESULTS:
[[file:figs/huddle_test_asd.png]]
* Comparison Interferometer / Encoder
* Dynamics from Actuator to Encoder
<<sec:dynamics_encoder>>
** Introduction :ignore:
The goal here is to make sure that the interferometer and encoder measurements are coherent.
We may see non-linearity in the interferometric measurement.
Now the dynamics from the force actuator to the measurement by the encoder is identified.
** Matlab Init :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@ -169,181 +202,72 @@ We may see non-linearity in the interferometric measurement.
#+end_src
** Load Data
#+begin_src matlab
load('mat/int_enc_comp.mat', 'interferometer', 'encoder', 'u', 't');
#+end_src
#+begin_src matlab
interferometer = detrend(interferometer, 0);
encoder = detrend(encoder, 0);
u = detrend(u, 0);
#+end_src
** Time Domain Results
#+begin_src matlab :exports none
figure;
hold on;
plot(t, encoder, '-', 'DisplayName', 'Encoder')
plot(t, interferometer, '--', 'DisplayName', 'Interferometer')
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
legend('location', 'northeast');
xlim([50, 52])
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/int_enc_one_cycle.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:int_enc_one_cycle
#+caption: One cycle measurement
#+RESULTS:
[[file:figs/int_enc_one_cycle.png]]
#+begin_src matlab :exports none
figure;
hold on;
plot(t, encoder - interferometer, 'DisplayName', 'Difference')
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
legend('location', 'northeast');
xlim([50, 52])
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/int_enc_one_cycle_error.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:int_enc_one_cycle_error
#+caption: Difference between the Encoder and the interferometer during one cycle
#+RESULTS:
[[file:figs/int_enc_one_cycle_error.png]]
** Difference between Encoder and Interferometer as a function of time
#+begin_src matlab
Ts = 1e-4;
d_i_mean = reshape(interferometer, [2/Ts floor(Ts/2*length(interferometer))]);
d_e_mean = reshape(encoder, [2/Ts floor(Ts/2*length(encoder))]);
#+end_src
#+begin_src matlab
w0 = 2*pi*5; % [rad/s]
xi = 0.7;
G_lpf = 1/(1 + 2*xi/w0*s + s^2/w0^2);
d_err_mean = reshape(lsim(G_lpf, encoder - interferometer, t), [2/Ts floor(Ts/2*length(encoder))]);
d_err_mean = d_err_mean - mean(d_err_mean);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
for i_i = 1:size(d_err_mean, 2)
plot(t(1:size(d_err_mean, 1)), d_err_mean(:, i_i), 'k-')
end
plot(t(1:size(d_err_mean, 1)), mean(d_err_mean, 2), 'r-')
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/int_enc_error_mean_time.pdf', 'width', 'wide', 'height', 'normal', 'pdf', false);
#+end_src
#+name: fig:int_enc_error_mean_time
#+caption: Difference between the two measurement in the time domain, averaged for all the cycles
#+RESULTS:
[[file:figs/int_enc_error_mean_time.png]]
** Difference between Encoder and Interferometer as a function of position
Compute the mean of the interferometer measurement corresponding to each of the encoder measurement.
#+begin_src matlab
[e_sorted, ~, e_ind] = unique(encoder);
i_mean = zeros(length(e_sorted), 1);
for i = 1:length(e_sorted)
i_mean(i) = mean(interferometer(e_ind == i));
end
i_mean_error = (i_mean - e_sorted);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
% plot(encoder, interferometer - encoder, 'k.', 'DisplayName', 'Difference')
plot(1e6*(e_sorted), 1e9*(i_mean_error))
hold off;
xlabel('Encoder Measurement [$\mu m$]'); ylabel('Measrement Error [nm]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/int_enc_error_mean_position.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:int_enc_error_mean_position
#+caption: Difference between the two measurement as a function of the measured position by the encoder, averaged for all the cycles
#+RESULTS:
[[file:figs/int_enc_error_mean_position.png]]
The period of the non-linearity seems to be $1.53 \mu m$ which corresponds to the wavelength of the Laser.
#+begin_src matlab
win_length = 1530; % length of the windows (corresponds to 1.53 um)
num_avg = floor(length(e_sorted)/win_length); % number of averaging
i_init = ceil((length(e_sorted) - win_length*num_avg)/2); % does not start at the extremity
e_sorted_mean_over_period = mean(reshape(i_mean_error(i_init:i_init+win_length*num_avg-1), [win_length num_avg]), 2);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(1e-3*(0:win_length-1), 1e9*(e_sorted_mean_over_period))
hold off;
xlabel('Displacement [$\mu m$]'); ylabel('Measurement Non-Linearity [nm]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/int_non_linearity_period_wavelength.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:int_non_linearity_period_wavelength
#+caption: Non-Linearity of the Interferometer over the period of the wavelength
#+RESULTS:
[[file:figs/int_non_linearity_period_wavelength.png]]
* Identification
** 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
** Load Data
As usual, the measurement data are loaded.
#+begin_src matlab
load('mat/int_enc_id_noise_bis.mat', 'interferometer', 'encoder', 'u', 't');
#+end_src
The first 0.1 seconds are removed as it corresponds to transient behavior.
#+begin_src matlab
interferometer = interferometer(t>0.1);
encoder = encoder(t>0.1);
u = u(t>0.1);
t = t(t>0.1);
#+end_src
Finally the offset are removed using the =detrend= command.
#+begin_src matlab
interferometer = detrend(interferometer, 0);
encoder = detrend(encoder, 0);
u = detrend(u, 0);
#+end_src
** Excitation and Measured Signals
The excitation signal is a white noise filtered by a low pass filter to not excite too much the high frequency modes.
The excitation signal is shown in Figure [[fig:encoder_identification_excitation_time]].
#+begin_src matlab :exports none
figure;
plot(t, u);
xlabel('Time [s]'); ylabel('Voltage [V]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/encoder_identification_excitation_time.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:encoder_identification_excitation_time
#+caption:
#+RESULTS:
[[file:figs/encoder_identification_excitation_time.png]]
The measured motion by the interferometer and encoder is shown in Figure
#+begin_src matlab :exports none
figure;
hold on;
plot(t, interferometer, 'DisplayName', 'Interferometer');
plot(t, encoder, 'DisplayName', 'Encoder');
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
legend('location', 'southeast');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/encoder_identification_motion.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:encoder_identification_motion
#+caption:
#+RESULTS:
[[file:figs/encoder_identification_motion.png]]
** Identification
Now the dynamics from the voltage sent to the voltage amplitude driving the APA95ML to the measured displacement by both the encoder and interferometer are computed.
#+begin_src matlab
Ts = 1e-4; % Sampling Time [s]
win = hann(ceil(10/Ts));
#+end_src
#+begin_src matlab
[tf_i_est, f] = tfestimate(u, interferometer, win, [], [], 1/Ts);
[co_i_est, ~] = mscohere(u, interferometer, win, [], [], 1/Ts);
@ -351,19 +275,23 @@ The period of the non-linearity seems to be $1.53 \mu m$ which corresponds to th
[co_e_est, ~] = mscohere(u, encoder, win, [], [], 1/Ts);
#+end_src
The obtained coherence is shown in Figure [[fig:identification_dynamics_coherence]].
It is shown that the identification is good until 500Hz for the interferometer and until 1kHz for the encoder.
#+begin_src matlab :exports none
figure;
hold on;
plot(f, co_i_est, '-')
plot(f, co_e_est, '-')
plot(f, co_i_est, '-', 'DisplayName', 'Interferometer')
plot(f, co_e_est, '-', 'DisplayName', 'Encoder')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([0.5, 5e3]);
legend('location', 'southwest');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/identification_dynamics_coherence.pdf', 'width', 'normal', 'height', 'normal');
exportFig('figs/identification_dynamics_coherence.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:identification_dynamics_coherence
@ -371,27 +299,29 @@ The period of the non-linearity seems to be $1.53 \mu m$ which corresponds to th
#+RESULTS:
[[file:figs/identification_dynamics_coherence.png]]
The compared dynamics as measured by the intereferometer and encoder are shown in Figure [[fig:identification_dynamics_bode]].
#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = nexttile;
ax1 = nexttile([2, 1]);
hold on;
plot(f, abs(tf_i_est), '-', 'DisplayName', 'Int')
plot(f, abs(tf_e_est), '-', 'DisplayName', 'Enc')
plot(f, abs(tf_i_est), '-', 'DisplayName', 'Interferometer')
plot(f, abs(tf_e_est), '-', 'DisplayName', 'Encoder')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-7, 3e-4]);
legend('location', 'southwest');
ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(tf_i_est), '-')
plot(f, 180/pi*angle(tf_e_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
hold off;
yticks(-360:90:360);
axis padded 'auto x'
@ -408,494 +338,5 @@ The period of the non-linearity seems to be $1.53 \mu m$ which corresponds to th
#+RESULTS:
[[file:figs/identification_dynamics_bode.png]]
* Change of Stiffness due to Sensors stack being open/closed circuit
** 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
** Load Data
#+begin_src matlab
oc = load('./mat/identification_open_circuit.mat', 't', 'encoder', 'u');
sc = load('./mat/identification_short_circuit.mat', 't', 'encoder', 'u');
#+end_src
** Transfer Functions
#+begin_src matlab
Ts = 1e-4; % Sampling Time [s]
win = hann(ceil(10/Ts));
#+end_src
#+begin_src matlab
[tf_oc_est, f] = tfestimate(oc.u, oc.encoder, win, [], [], 1/Ts);
[co_oc_est, ~] = mscohere( oc.u, oc.encoder, win, [], [], 1/Ts);
[tf_sc_est, ~] = tfestimate(sc.u, sc.encoder, win, [], [], 1/Ts);
[co_sc_est, ~] = mscohere( sc.u, sc.encoder, win, [], [], 1/Ts);
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(f, co_oc_est, '-')
plot(f, co_sc_est, '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([0.5, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/stiffness_force_sensor_coherence.pdf', 'width', 'wide', 'height', 'normal');
#+end_src
#+name: fig:stiffness_force_sensor_coherence
#+caption:
#+RESULTS:
[[file:figs/stiffness_force_sensor_coherence.png]]
#+begin_src matlab :exports none
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = nexttile;
hold on;
plot(f, abs(tf_oc_est), '-', 'DisplayName', 'Open-Circuit')
plot(f, abs(tf_sc_est), '-', 'DisplayName', 'Short-Circuit')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
hold off;
ylim([1e-7, 3e-4]);
legend('location', 'southwest');
ax2 = nexttile;
hold on;
plot(f, 180/pi*angle(tf_oc_est), '-')
plot(f, 180/pi*angle(tf_sc_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
yticks(-360:90:360);
axis padded 'auto x'
linkaxes([ax1,ax2], 'x');
xlim([0.5, 5e3]);
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/stiffness_force_sensor_bode.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:stiffness_force_sensor_bode
#+caption:
#+RESULTS:
[[file:figs/stiffness_force_sensor_bode.png]]
#+begin_src matlab :tangle no :exports results :results file replace
xlim([180, 280]);
exportFig('figs/stiffness_force_sensor_bode_zoom.pdf', 'width', 'small', 'height', 'tall');
#+end_src
#+name: fig:stiffness_force_sensor_bode_zoom
#+caption: Zoom on the change of resonance
#+RESULTS:
[[file:figs/stiffness_force_sensor_bode_zoom.png]]
#+begin_important
The change of resonance frequency / stiffness is very small and is not important here.
#+end_important
* Generated Number of Charge / Voltage
** Introduction :ignore:
Two stacks are used as actuator (in parallel) and one stack is used as sensor.
The amplifier gain is 20V/V (Cedrat LA75B).
** 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
** Steps
#+begin_src matlab
load('./mat/force_sensor_steps.mat', 't', 'encoder', 'u', 'v');
#+end_src
#+begin_src matlab
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
nexttile;
plot(t, v);
xlabel('Time [s]'); ylabel('Measured voltage [V]');
nexttile;
plot(t, u);
xlabel('Time [s]'); ylabel('Actuator Voltage [V]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sen_steps_time_domain.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:force_sen_steps_time_domain
#+caption: Time domain signal during the 3 actuator voltage steps
#+RESULTS:
[[file:figs/force_sen_steps_time_domain.png]]
Three steps are performed at the following time intervals:
#+begin_src matlab
t_s = [ 2.5, 23;
23.8, 35;
35.8, 50];
#+end_src
Fit function:
#+begin_src matlab
f = @(b,x) b(1).*exp(b(2).*x) + b(3);
#+end_src
We are interested by the =b(2)= term, which is the time constant of the exponential.
#+begin_src matlab
tau = zeros(size(t_s, 1),1);
V0 = zeros(size(t_s, 1),1);
#+end_src
#+begin_src matlab
for t_i = 1:size(t_s, 1)
t_cur = t(t_s(t_i, 1) < t & t < t_s(t_i, 2));
t_cur = t_cur - t_cur(1);
y_cur = v(t_s(t_i, 1) < t & t < t_s(t_i, 2));
nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
B0 = [0.5, -0.15, 2.2]; % Choose Appropriate Initial Estimates
[B,rnrm] = fminsearch(nrmrsd, B0); % Estimate Parameters B
tau(t_i) = 1/B(2);
V0(t_i) = B(3);
end
#+end_src
#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
data2orgtable([abs(tau), V0], {}, {'$tau$ [s]', '$V_0$ [V]'}, ' %.2f ');
#+end_src
#+RESULTS:
| $tau$ [s] | $V_0$ [V] |
|-----------+-----------|
| 6.47 | 2.26 |
| 6.76 | 2.26 |
| 6.49 | 2.25 |
With the capacitance being $C = 4.4 \mu F$, the internal impedance of the Speedgoat ADC can be computed as follows:
#+begin_src matlab
Cp = 4.4e-6; % [F]
Rin = abs(mean(tau))/Cp;
#+end_src
#+begin_src matlab :results value replace :exports results
ans = Rin
#+end_src
#+RESULTS:
: 1494100.0
The input impedance of the Speedgoat's ADC should then be close to $1.5\,M\Omega$ (specified at $1\,M\Omega$).
#+begin_important
How can we explain the voltage offset?
#+end_important
As shown in Figure [[fig:piezo_sensor_model_instrumentation]] (taken from cite:reza06_piezoel_trans_vibrat_contr_dampin), an input voltage offset is due to the input bias current $i_n$.
#+name: fig:piezo_sensor_model_instrumentation
#+caption: Model of a piezoelectric transducer (left) and instrumentation amplifier (right)
[[file:figs/piezo_sensor_model_instrumentation.png]]
The estimated input bias current is then:
#+begin_src matlab
in = mean(V0)/Rin;
#+end_src
#+begin_src matlab :results value replace :exports results
ans = in
#+end_src
#+RESULTS:
: 1.5119e-06
An additional resistor in parallel with $R_{in}$ would have two effects:
- reduce the input voltage offset
\[ V_{off} = \frac{R_a R_{in}}{R_a + R_{in}} i_n \]
- increase the high pass corner frequency $f_c$
\[ C_p \frac{R_{in}R_a}{R_{in} + R_a} = \tau_c = \frac{1}{f_c} \]
\[ R_a = \frac{R_i}{f_c C_p R_i - 1} \]
If we allow the high pass corner frequency to be equals to 3Hz:
#+begin_src matlab
fc = 3;
Ra = Rin/(fc*Cp*Rin - 1);
#+end_src
#+begin_src matlab :results value replace :exports results
ans = Ra
#+end_src
#+RESULTS:
: 79804
With this parallel resistance value, the voltage offset would be:
#+begin_src matlab
V_offset = Ra*Rin/(Ra + Rin) * in;
#+end_src
#+begin_src matlab :results value replace :exports results
ans = V_offset
#+end_src
#+RESULTS:
: 0.11454
Which is much more acceptable.
** Add Parallel Resistor
A resistor of $\approx 100\,k\Omega$ is added in parallel with the force sensor and the same kin.
#+begin_src matlab
load('./mat/force_sensor_steps_R_82k7.mat', 't', 'encoder', 'u', 'v');
#+end_src
#+begin_src matlab
figure;
tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');
nexttile;
plot(t, v);
xlabel('Time [s]'); ylabel('Measured voltage [V]');
nexttile;
plot(t, u);
xlabel('Time [s]'); ylabel('Actuator Voltage [V]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sen_steps_time_domain_par_R.pdf', 'width', 'wide', 'height', 'tall');
#+end_src
#+name: fig:force_sen_steps_time_domain_par_R
#+caption: Time domain signal during the actuator voltage steps
#+RESULTS:
[[file:figs/force_sen_steps_time_domain_par_R.png]]
Three steps are performed at the following time intervals:
#+begin_src matlab
t_s = [1.9, 6;
8.5, 13;
15.5, 21;
22.6, 26;
30.0, 36;
37.5, 41;
46.2, 49.5]
#+end_src
Fit function:
#+begin_src matlab
f = @(b,x) b(1).*exp(b(2).*x) + b(3);
#+end_src
We are interested by the =b(2)= term, which is the time constant of the exponential.
#+begin_src matlab
tau = zeros(size(t_s, 1),1);
V0 = zeros(size(t_s, 1),1);
#+end_src
#+begin_src matlab
for t_i = 1:size(t_s, 1)
t_cur = t(t_s(t_i, 1) < t & t < t_s(t_i, 2));
t_cur = t_cur - t_cur(1);
y_cur = v(t_s(t_i, 1) < t & t < t_s(t_i, 2));
nrmrsd = @(b) norm(y_cur - f(b,t_cur)); % Residual Norm Cost Function
B0 = [0.5, -0.2, 0.2]; % Choose Appropriate Initial Estimates
[B,rnrm] = fminsearch(nrmrsd, B0); % Estimate Parameters B
tau(t_i) = 1/B(2);
V0(t_i) = B(3);
end
#+end_src
And indeed, we obtain a much smaller offset voltage and a much faster time constant.
#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
data2orgtable([abs(tau), V0], {}, {'$tau$ [s]', '$V_0$ [V]'}, ' %.2f ');
#+end_src
#+RESULTS:
| $tau$ [s] | $V_0$ [V] |
|-----------+-----------|
| 0.43 | 0.15 |
| 0.45 | 0.16 |
| 0.43 | 0.15 |
| 0.43 | 0.15 |
| 0.45 | 0.15 |
| 0.46 | 0.16 |
| 0.48 | 0.16 |
Knowing the capacitance value, we can estimate the value of the added resistor (neglecting the input impedance of $\approx 1\,M\Omega$):
#+begin_src matlab
Cp = 4.4e-6; % [F]
Rin = abs(mean(tau))/Cp;
#+end_src
#+begin_src matlab :results value replace :exports results
ans = Rin
#+end_src
#+RESULTS:
: 101200.0
And we can verify that the bias current estimation stays the same:
#+begin_src matlab
in = mean(V0)/Rin;
#+end_src
#+begin_src matlab :results value replace :exports results
ans = in
#+end_src
#+RESULTS:
: 1.5305e-06
This validates the model of the ADC and the effectiveness of the added resistor.
** Sinus
#+begin_src matlab
load('./mat/force_sensor_sin.mat', 't', 'encoder', 'u', 'v');
u = u(t>25);
v = v(t>25);
encoder = encoder(t>25) - mean(encoder(t>25));
t = t(t>25);
#+end_src
The driving voltage is a sinus at 0.5Hz centered on 3V and with an amplitude of 3V (Figure [[fig:force_sensor_sin_u]]).
#+begin_src matlab :exports none
figure;
plot(t, u)
xlabel('Time [s]'); ylabel('Control Voltage [V]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sensor_sin_u.pdf', 'width', 'normal', 'height', 'small');
#+end_src
#+name: fig:force_sensor_sin_u
#+caption: Driving Voltage
#+RESULTS:
[[file:figs/force_sensor_sin_u.png]]
The full stroke as measured by the encoder is:
#+begin_src matlab :results value replace
max(encoder)-min(encoder)
#+end_src
#+RESULTS:
: 5.005e-05
Its signal is shown in Figure [[fig:force_sensor_sin_encoder]].
#+begin_src matlab :exports none
figure;
plot(t, encoder)
xlabel('Time [s]'); ylabel('Encoder [m]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sensor_sin_encoder.pdf', 'width', 'normal', 'height', 'small');
#+end_src
#+name: fig:force_sensor_sin_encoder
#+caption: Encoder measurement
#+RESULTS:
[[file:figs/force_sensor_sin_encoder.png]]
The generated voltage by the stack is shown in Figure
#+begin_src matlab :exports none
figure;
plot(t, v)
xlabel('Time [s]'); ylabel('Force Sensor Output [V]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sensor_sin_stack.pdf', 'width', 'normal', 'height', 'small');
#+end_src
#+name: fig:force_sensor_sin_stack
#+caption: Voltage measured on the stack used as a sensor
#+RESULTS:
[[file:figs/force_sensor_sin_stack.png]]
The capacitance of the stack is
#+begin_src matlab
Cp = 4.4e-6; % [F]
#+end_src
The corresponding generated charge is then shown in Figure [[fig:force_sensor_sin_charge]].
#+begin_src matlab :exports none
figure;
plot(t, 1e6*Cp*(v-mean(v)))
xlabel('Time [s]'); ylabel('Generated Charge [$\mu C$]');
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sensor_sin_charge.pdf', 'width', 'normal', 'height', 'small');
#+end_src
#+name: fig:force_sensor_sin_charge
#+caption: Generated Charge
#+RESULTS:
[[file:figs/force_sensor_sin_charge.png]]
The relation between the generated voltage and the measured displacement is almost linear as shown in Figure [[fig:force_sensor_linear_relation]].
#+begin_src matlab
b1 = encoder\(v-mean(v));
#+end_src
#+begin_src matlab :exports none
figure;
hold on;
plot(encoder, v-mean(v), 'DisplayName', 'Measured Voltage');
plot(encoder, encoder*b1, 'DisplayName', sprintf('Linear Fit: $U_s \\approx %.3f [V/\\mu m] \\cdot d$', 1e-6*abs(b1)));
hold off;
xlabel('Measured Displacement [m]'); ylabel('Generated Voltage [V]');
legend();
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
exportFig('figs/force_sensor_linear_relation.pdf', 'width', 'normal', 'height', 'small');
#+end_src
#+name: fig:force_sensor_linear_relation
#+caption: Almost linear relation between the relative displacement and the generated voltage
#+RESULTS:
[[file:figs/force_sensor_linear_relation.png]]
With a 16bits ADC, the resolution will then be equals to (in [nm]):
#+begin_src matlab :results value replace
abs((20/2^16)/(b1/1e9))
#+end_src
#+RESULTS:
: 3.9838
The second resonance at around 900Hz most likely corresponds to the resonance of either the ruler support or the head support.

View File

@ -1 +0,0 @@
!function(a,b){"use strict";function c(c,g){var h=this;h.$el=a(c),h.el=c,h.id=e++,h.$window=a(b),h.$document=a(document),h.$el.bind("destroyed",a.proxy(h.teardown,h)),h.$clonedHeader=null,h.$originalHeader=null,h.isSticky=!1,h.hasBeenSticky=!1,h.leftOffset=null,h.topOffset=null,h.init=function(){h.$el.each(function(){var b=a(this);b.css("padding",0),h.$originalHeader=a("thead:first",this),h.$clonedHeader=h.$originalHeader.clone(),b.trigger("clonedHeader."+d,[h.$clonedHeader]),h.$clonedHeader.addClass("tableFloatingHeader"),h.$clonedHeader.css("display","none"),h.$originalHeader.addClass("tableFloatingHeaderOriginal"),h.$originalHeader.after(h.$clonedHeader),h.$printStyle=a('<style type="text/css" media="print">.tableFloatingHeader{display:none !important;}.tableFloatingHeaderOriginal{position:static !important;}</style>'),a("head").append(h.$printStyle)}),h.setOptions(g),h.updateWidth(),h.toggleHeaders(),h.bind()},h.destroy=function(){h.$el.unbind("destroyed",h.teardown),h.teardown()},h.teardown=function(){h.isSticky&&h.$originalHeader.css("position","static"),a.removeData(h.el,"plugin_"+d),h.unbind(),h.$clonedHeader.remove(),h.$originalHeader.removeClass("tableFloatingHeaderOriginal"),h.$originalHeader.css("visibility","visible"),h.$printStyle.remove(),h.el=null,h.$el=null},h.bind=function(){h.$scrollableArea.on("scroll."+d,h.toggleHeaders),h.isWindowScrolling||(h.$window.on("scroll."+d+h.id,h.setPositionValues),h.$window.on("resize."+d+h.id,h.toggleHeaders)),h.$scrollableArea.on("resize."+d,h.toggleHeaders),h.$scrollableArea.on("resize."+d,h.updateWidth)},h.unbind=function(){h.$scrollableArea.off("."+d,h.toggleHeaders),h.isWindowScrolling||(h.$window.off("."+d+h.id,h.setPositionValues),h.$window.off("."+d+h.id,h.toggleHeaders)),h.$scrollableArea.off("."+d,h.updateWidth)},h.toggleHeaders=function(){h.$el&&h.$el.each(function(){var b,c=a(this),d=h.isWindowScrolling?isNaN(h.options.fixedOffset)?h.options.fixedOffset.outerHeight():h.options.fixedOffset:h.$scrollableArea.offset().top+(isNaN(h.options.fixedOffset)?0:h.options.fixedOffset),e=c.offset(),f=h.$scrollableArea.scrollTop()+d,g=h.$scrollableArea.scrollLeft(),i=h.isWindowScrolling?f>e.top:d>e.top,j=(h.isWindowScrolling?f:0)<e.top+c.height()-h.$clonedHeader.height()-(h.isWindowScrolling?0:d);i&&j?(b=e.left-g+h.options.leftOffset,h.$originalHeader.css({position:"fixed","margin-top":h.options.marginTop,left:b,"z-index":3}),h.leftOffset=b,h.topOffset=d,h.$clonedHeader.css("display",""),h.isSticky||(h.isSticky=!0,h.updateWidth()),h.setPositionValues()):h.isSticky&&(h.$originalHeader.css("position","static"),h.$clonedHeader.css("display","none"),h.isSticky=!1,h.resetWidth(a("td,th",h.$clonedHeader),a("td,th",h.$originalHeader)))})},h.setPositionValues=function(){var a=h.$window.scrollTop(),b=h.$window.scrollLeft();!h.isSticky||0>a||a+h.$window.height()>h.$document.height()||0>b||b+h.$window.width()>h.$document.width()||h.$originalHeader.css({top:h.topOffset-(h.isWindowScrolling?0:a),left:h.leftOffset-(h.isWindowScrolling?0:b)})},h.updateWidth=function(){if(h.isSticky){h.$originalHeaderCells||(h.$originalHeaderCells=a("th,td",h.$originalHeader)),h.$clonedHeaderCells||(h.$clonedHeaderCells=a("th,td",h.$clonedHeader));var b=h.getWidth(h.$clonedHeaderCells);h.setWidth(b,h.$clonedHeaderCells,h.$originalHeaderCells),h.$originalHeader.css("width",h.$clonedHeader.width())}},h.getWidth=function(c){var d=[];return c.each(function(c){var e,f=a(this);if("border-box"===f.css("box-sizing"))e=f[0].getBoundingClientRect().width;else{var g=a("th",h.$originalHeader);if("collapse"===g.css("border-collapse"))if(b.getComputedStyle)e=parseFloat(b.getComputedStyle(this,null).width);else{var i=parseFloat(f.css("padding-left")),j=parseFloat(f.css("padding-right")),k=parseFloat(f.css("border-width"));e=f.outerWidth()-i-j-k}else e=f.width()}d[c]=e}),d},h.setWidth=function(a,b,c){b.each(function(b){var d=a[b];c.eq(b).css({"min-width":d,"max-width":d})})},h.resetWidth=function(b,c){b.each(function(b){var d=a(this);c.eq(b).css({"min-width":d.css("min-width"),"max-width":d.css("max-width")})})},h.setOptions=function(c){h.options=a.extend({},f,c),h.$scrollableArea=a(h.options.scrollableArea),h.isWindowScrolling=h.$scrollableArea[0]===b},h.updateOptions=function(a){h.setOptions(a),h.unbind(),h.bind(),h.updateWidth(),h.toggleHeaders()},h.init()}var d="stickyTableHeaders",e=0,f={fixedOffset:0,leftOffset:0,marginTop:0,scrollableArea:b};a.fn[d]=function(b){return this.each(function(){var e=a.data(this,"plugin_"+d);e?"string"==typeof b?e[b].apply(e):e.updateOptions(b):"destroy"!==b&&a.data(this,"plugin_"+d,new c(this,b))})}}(jQuery,window);

View File

@ -9,6 +9,8 @@ $(function() {
$('.hint').before("<p class='admonition-title hint'>Hint</p>");
$('.error').before("<p class='admonition-title error'>Error</p>");
$('.danger').before("<p class='admonition-title danger'>Danger</p>");
$('.question').before("<p class='admonition-title question'>Question</p>");
$('.summary').before("<p class='admonition-title hint'>Summary</p>");
});
$( document ).ready(function() {

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

Binary file not shown.

View File

@ -1,616 +0,0 @@
%% Huddle Test
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
% Detrend Data
ht.d = detrend(ht.d, 0);
ht.acc_1 = detrend(ht.acc_1, 0);
ht.acc_2 = detrend(ht.acc_2, 0);
ht.geo_1 = detrend(ht.geo_1, 0);
ht.geo_2 = detrend(ht.geo_2, 0);
ht.f_meas = detrend(ht.f_meas, 0);
% Compute PSD
run setup;
win = hann(ceil(10/Ts));
[p_d, f] = pwelch(ht.d, win, [], [], 1/Ts);
[p_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
[p_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
[p_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
[p_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
[p_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts);
% Plot PSD
figure;
hold on;
plot(f, p_acc1);
plot(f, p_acc2);
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Accelerometers')
figure;
hold on;
plot(f, p_geo1);
plot(f, p_geo2);
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Geophones')
figure;
hold on;
plot(f, p_d);
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$m^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Interferometers')
figure;
hold on;
plot(f, p_fmeas);
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Force Sensor')
%% Accelerometer and Geophone Models
% Accelerometer used: https://www.pcb.com/products?model=393B05
% Geophone used: L22 https://www.sercel.com/products/Lists/ProductSpecification/Geophones_brochure_Sercel_EN.pdf
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_geo = 120*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)]
% PSD of intertial sensors in [m^2/Hz]
figure;
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, sqrt(p_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
'DisplayName', 'Accelerometer');
set(gca, 'ColorOrderIndex', 1);
plot(f, sqrt(p_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
'HandleVisibility', 'off');
set(gca, 'ColorOrderIndex', 2);
plot(f, sqrt(p_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
'DisplayName', 'Geophone');
set(gca, 'ColorOrderIndex', 2);
plot(f, sqrt(p_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
'HandleVisibility', 'off');
set(gca, 'ColorOrderIndex', 3);
plot(f, sqrt(p_d), 'DisplayName', 'Interferometer');
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
title('Huddle Test')
legend();
%% Compare Theoretical model with identified one
id_ol = load('./mat/identification_noise_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
% Detrend Data
id_ol.d = detrend(id_ol.d, 0);
id_ol.acc_1 = detrend(id_ol.acc_1, 0);
id_ol.acc_2 = detrend(id_ol.acc_2, 0);
id_ol.geo_1 = detrend(id_ol.geo_1, 0);
id_ol.geo_2 = detrend(id_ol.geo_2, 0);
id_ol.f_meas = detrend(id_ol.f_meas, 0);
id_ol.u = detrend(id_ol.u, 0);
% Identification Parameters
run setup;
win = hann(ceil(10/Ts));
% IFF Plant
[tf_fmeas_est, f] = tfestimate(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts); % [V/m]
[co_fmeas_est, ~] = mscohere(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_fmeas_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_fmeas_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([40, 400]);
% Geophones
[tf_geo1_est, ~] = tfestimate(id_ol.d, id_ol.geo_1, win, [], [], 1/Ts); % [V/m]
[co_geo1_est, ~] = mscohere(id_ol.d, id_ol.geo_1, win, [], [], 1/Ts);
[tf_geo2_est, ~] = tfestimate(id_ol.d, id_ol.geo_2, win, [], [], 1/Ts); % [V/m]
[co_geo2_est, ~] = mscohere(id_ol.d, id_ol.geo_2, win, [], [], 1/Ts);
figure;
ax1 = subplot(2, 1, 1);
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, abs(tf_geo1_est), '.')
set(gca, 'ColorOrderIndex', 1);
plot(f, abs(tf_geo2_est), '.')
plot(f, abs(squeeze(freqresp(G_geo, f, 'Hz')).*(1i*2*pi*f)), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, 180/pi*angle(tf_geo1_est), '.')
set(gca, 'ColorOrderIndex', 1);
plot(f, 180/pi*angle(tf_geo2_est), '.')
plot(f, 180/pi*angle(-squeeze(freqresp(G_geo, f, 'Hz')).*(1i*2*pi*f)), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([40, 400]);
% Accelerometers
[tf_acc1_est, ~] = tfestimate(id_ol.d, id_ol.acc_1, win, [], [], 1/Ts); % [V/m]
[co_acc1_est, ~] = mscohere(id_ol.d, id_ol.acc_1, win, [], [], 1/Ts);
[tf_acc2_est, ~] = tfestimate(id_ol.d, id_ol.acc_2, win, [], [], 1/Ts); % [V/m]
[co_acc2_est, ~] = mscohere(id_ol.d, id_ol.acc_2, win, [], [], 1/Ts);
figure;
ax1 = subplot(2, 1, 1);
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, abs(tf_acc1_est), '.')
set(gca, 'ColorOrderIndex', 1);
plot(f, abs(tf_acc2_est), '.')
plot(f, abs(squeeze(freqresp(G_acc, f, 'Hz')).*(1i*2*pi*f).^2), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, 180/pi*angle(tf_acc1_est), '.')
set(gca, 'ColorOrderIndex', 1);
plot(f, 180/pi*angle(tf_acc2_est), '.')
plot(f, 180/pi*angle(squeeze(freqresp(G_acc, f, 'Hz')).*(1i*2*pi*f).^2), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([40, 400]);
%% IFF development
[tf_fmeas_est, f] = tfestimate(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts); % [V/m]
[co_fmeas_est, ~] = mscohere(id_ol.u, id_ol.f_meas, win, [], [], 1/Ts);
% Model
wz = 2*pi*103;
xi_z = 0.01;
wp = 2*pi*238;
xi_p = 0.015;
Giff = 20*(s^2 + 2*xi_z*s*wz + wz^2)/(s^2 + 2*xi_p*s*wp + wp^2)*(s/3/pi/(1 + s/3/pi));
% Comparison model and identification
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_fmeas_est), '.')
plot(f, abs(squeeze(freqresp(Giff, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_fmeas_est), '.')
plot(f, 180/pi*angle(squeeze(freqresp(Giff, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([40, 400]);
% Root Locus
gains = logspace(0, 5, 1000);
figure;
hold on;
plot(real(pole(Giff)), imag(pole(Giff)), 'kx');
plot(real(tzero(Giff)), imag(tzero(Giff)), 'ko');
for i = 1:length(gains)
cl_poles = pole(feedback(Giff, gains(i)/(s + 2*pi*2)));
plot(real(cl_poles), imag(cl_poles), 'k.');
end
ylim([0, 1800]);
xlim([-1600,200]);
xlabel('Real Part')
ylabel('Imaginary Part')
axis square
% Optimal Controller
Kiff_opt = 110/(s + 2*pi*2);
%% New identification
id_ol = load('./mat/identification_chirp_40_400.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
id_cl = load('./mat/identification_chirp_40_400_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
% Used controller
Kiff = -110/(s + 2*pi*2);
[tf_G_ol_est, f] = tfestimate(id_ol.u, id_ol.d, win, [], [], 1/Ts);
[co_G_ol_est, ~] = mscohere(id_ol.u, id_ol.d, win, [], [], 1/Ts);
[tf_G_cl_est, ~] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts);
[co_G_cl_est, ~] = mscohere(id_cl.u, id_cl.d, win, [], [], 1/Ts);
figure;
hold on;
plot(f, co_G_ol_est, '-')
plot(f, co_G_cl_est, '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([40, 400]); ylim([0, 1])
% Comparison model and identification
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_G_ol_est), '-')
plot(f, abs(tf_G_cl_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_G_ol_est), '-')
plot(f, 180/pi*angle(tf_G_cl_est), '-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([40, 400]);
%% Excitation Signal
run setup;
% Get trasnfer function from input [V] to output displacement [m]
id_cl = load('./mat/identification_noise_iff_bis.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
win = hann(ceil(10/Ts));
[tf_G_cl_est, f] = tfestimate(id_cl.u, id_cl.d, win, [], [], 1/Ts);
[co_G_cl_est, ~] = mscohere(id_cl.u, id_cl.d, win, [], [], 1/Ts);
G_d_est = -5e-6*(2*pi*230)^2/(s^2 + 2*0.3*2*pi*240*s + (2*pi*240)^2);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_G_cl_est), '-')
plot(f, abs(squeeze(freqresp(G_d_est, f, 'Hz'))), '--')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude [m/V]'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_G_cl_est), '-')
plot(f, 180/pi*angle(squeeze(freqresp(G_d_est, f, 'Hz'))), '--')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([10, 1000]);
%
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
ht.d = detrend(ht.d, 0);
ht.acc_1 = detrend(ht.acc_1, 0);
ht.acc_2 = detrend(ht.acc_2, 0);
ht.geo_1 = detrend(ht.geo_1, 0);
ht.geo_2 = detrend(ht.geo_2, 0);
win = hann(ceil(10/Ts));
[p_d, f] = pwelch(ht.d, win, [], [], 1/Ts);
[p_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
[p_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
[p_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
[p_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
% Generate Time domain signal with wanted PSD
Fs = 1/Ts; % Sampling Frequency [Hz]
t = 0:Ts:180; % Time Vector [s]
u = sqrt(Fs/2)*randn(length(t), 1); % Signal with an ASD equal to one
G_exc = 0.2e-6/(1 + s/2/pi/2)/(1 + s/2/pi/50);
y_d = lsim(G_exc, u, t);
[pxx, ~] = pwelch(y_d, win, 0, [], Fs);
figure;
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, sqrt(p_acc1)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
'DisplayName', 'Accelerometer');
set(gca, 'ColorOrderIndex', 1);
plot(f, sqrt(p_acc2)./abs(squeeze(freqresp(G_acc*s^2, f, 'Hz'))), ...
'HandleVisibility', 'off');
set(gca, 'ColorOrderIndex', 2);
plot(f, sqrt(p_geo1)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
'DisplayName', 'Geophone');
set(gca, 'ColorOrderIndex', 2);
plot(f, sqrt(p_geo2)./abs(squeeze(freqresp(G_geo*s, f, 'Hz'))), ...
'HandleVisibility', 'off');
plot(f, sqrt(pxx), 'k-', ...
'DisplayName', 'Excitation');
set(gca, 'ColorOrderIndex', 3);
plot(f, sqrt(p_d), 'DisplayName', 'Interferometer');
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('ASD [$m/\sqrt{Hz}$]'); xlabel('Frequency [Hz]');
title('Huddle Test')
legend();
% From displacement to Voltage
y_v = lsim(G_exc*(1 + s/2/pi/50)/G_d_est/(1 + s/2/pi/5e3), u, t);
figure; plot(t, y_v)
figure; plot(t, lsim(G_pf, y_v, t))
%% Transfer function of inertial sensors
load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
%% Estimation of the inertial sensor transfer functions
id = load('./mat/identification_noise_opt_iff.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
ht = load('./mat/huddle_test.mat', 'd', 'acc_1', 'acc_2', 'geo_1', 'geo_2', 'f_meas', 'u', 't');
ht.d = detrend(ht.d, 0);
ht.acc_1 = detrend(ht.acc_1, 0);
ht.acc_2 = detrend(ht.acc_2, 0);
ht.geo_1 = detrend(ht.geo_1, 0);
ht.geo_2 = detrend(ht.geo_2, 0);
ht.f_meas = detrend(ht.f_meas, 0);
id.d = detrend(id.d, 0);
id.acc_1 = detrend(id.acc_1, 0);
id.acc_2 = detrend(id.acc_2, 0);
id.geo_1 = detrend(id.geo_1, 0);
id.geo_2 = detrend(id.geo_2, 0);
id.f_meas = detrend(id.f_meas, 0);
% Compare PSD
run setup;
win = hann(ceil(10/Ts));
[p_id_d, f] = pwelch(id.d, win, [], [], 1/Ts);
[p_id_acc1, ~] = pwelch(id.acc_1, win, [], [], 1/Ts);
[p_id_acc2, ~] = pwelch(id.acc_2, win, [], [], 1/Ts);
[p_id_geo1, ~] = pwelch(id.geo_1, win, [], [], 1/Ts);
[p_id_geo2, ~] = pwelch(id.geo_2, win, [], [], 1/Ts);
[p_id_fmeas, ~] = pwelch(id.f_meas, win, [], [], 1/Ts);
[p_ht_d, ~] = pwelch(ht.d, win, [], [], 1/Ts);
[p_ht_acc1, ~] = pwelch(ht.acc_1, win, [], [], 1/Ts);
[p_ht_acc2, ~] = pwelch(ht.acc_2, win, [], [], 1/Ts);
[p_ht_geo1, ~] = pwelch(ht.geo_1, win, [], [], 1/Ts);
[p_ht_geo2, ~] = pwelch(ht.geo_2, win, [], [], 1/Ts);
[p_ht_fmeas, ~] = pwelch(ht.f_meas, win, [], [], 1/Ts);
figure;
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, p_ht_acc1, 'DisplayName', 'Huddle Test');
set(gca, 'ColorOrderIndex', 1);
plot(f, p_ht_acc2, 'HandleVisibility', 'off');
set(gca, 'ColorOrderIndex', 2);
plot(f, p_id_acc1, 'DisplayName', 'Identification Test');
set(gca, 'ColorOrderIndex', 2);
plot(f, p_id_acc2, 'HandleVisibility', 'off');
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Accelerometers')
legend();
figure;
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, p_ht_geo1, 'DisplayName', 'Huddle Test');
set(gca, 'ColorOrderIndex', 1);
plot(f, p_ht_geo2, 'HandleVisibility', 'off');
set(gca, 'ColorOrderIndex', 2);
plot(f, p_id_geo1, 'DisplayName', 'Identification Test');
set(gca, 'ColorOrderIndex', 2);
plot(f, p_id_geo2, 'HandleVisibility', 'off');
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$V^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Geophones')
legend();
figure;
hold on;
plot(f, p_ht_d, 'DisplayName', 'Huddle Test');
plot(f, p_id_d, 'DisplayName', 'Identification Test');
hold off;
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('PSD [$m^2/Hz$]'); xlabel('Frequency [Hz]');
title('Huddle Test - Interferometers')
legend();
% tf and coh computation
[tf_acc1_est, f] = tfestimate(id.d, id.acc_1, win, [], [], 1/Ts);
[co_acc1_est, ~] = mscohere(id.d, id.acc_1, win, [], [], 1/Ts);
[tf_acc2_est, ~] = tfestimate(id.d, id.acc_2, win, [], [], 1/Ts);
[co_acc2_est, ~] = mscohere(id.d, id.acc_2, win, [], [], 1/Ts);
[tf_geo1_est, ~] = tfestimate(id.d, id.geo_1, win, [], [], 1/Ts);
[co_geo1_est, ~] = mscohere(id.d, id.geo_1, win, [], [], 1/Ts);
[tf_geo2_est, ~] = tfestimate(id.d, id.geo_2, win, [], [], 1/Ts);
[co_geo2_est, ~] = mscohere(id.d, id.geo_2, win, [], [], 1/Ts);
% Coherence
figure;
hold on;
set(gca, 'ColorOrderIndex', 1);
plot(f, co_acc1_est, '-', 'DisplayName', 'Accelerometer')
set(gca, 'ColorOrderIndex', 1);
plot(f, co_acc2_est, '-', 'HandleVisibility', 'off')
set(gca, 'ColorOrderIndex', 2);
plot(f, co_geo1_est, '-', 'DisplayName', 'Geophone')
set(gca, 'ColorOrderIndex', 2);
plot(f, co_geo2_est, '-', 'HandleVisibility', 'off')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Coherence'); xlabel('Frequency [Hz]');
hold off;
xlim([2, 2e3]); ylim([0, 1])
legend();
% Models
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [[V/(m/s)]
% Transfer Functions
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_acc1_est./(1i*2*pi*f).^2), '-')
plot(f, abs(tf_acc2_est./(1i*2*pi*f).^2), '-')
plot(f, abs(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude [V/(m/s^2)]'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_acc1_est./(1i*2*pi*f).^2), '-')
plot(f, 180/pi*angle(tf_acc2_est./(1i*2*pi*f).^2), '-')
plot(f, 180/pi*angle(squeeze(freqresp(G_acc, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([2, 2e3]);
figure;
ax1 = subplot(2, 1, 1);
hold on;
plot(f, abs(tf_geo1_est./(1i*2*pi*f)), '-')
plot(f, abs(tf_geo2_est./(1i*2*pi*f)), '-')
plot(f, abs(squeeze(freqresp(G_geo, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'log');
ylabel('Amplitude[V/(m/s)]'); xlabel('Frequency [Hz]');
hold off;
ax2 = subplot(2, 1, 2);
hold on;
plot(f, 180/pi*angle(tf_geo1_est./(1i*2*pi*f)), '-')
plot(f, 180/pi*angle(tf_geo2_est./(1i*2*pi*f)), '-')
plot(f, 180/pi*angle(squeeze(freqresp(G_geo, f, 'Hz'))), 'k-')
set(gca, 'Xscale', 'log'); set(gca, 'Yscale', 'lin');
ylabel('Phase'); xlabel('Frequency [Hz]');
hold off;
linkaxes([ax1,ax2], 'x');
xlim([0.5, 2e3]);
%% Compare signal
id.acc_1 = detrend(id.acc_1, 0);
id.acc_2 = detrend(id.acc_2, 0);
id.geo_1 = detrend(id.geo_1, 0);
id.geo_2 = detrend(id.geo_2, 0);
id.d = detrend(id.d, 0);
G_acc = 1/(1 + s/2/pi/2500); % [V/(m/s2)]
G_geo = -1200*s^2/(s^2 + 2*0.7*2*pi*2*s + (2*pi*2)^2); % [V/(m/s)]
G_hpf = (s/2/pi/2)/(1 + s/2/pi/2);
acc1_d = lsim(G_hpf*1/G_acc/(s + 2*pi)^2, id.acc_1, id.t);
acc2_d = lsim(G_hpf*1/G_acc/(s + 2*pi)^2, id.acc_2, id.t);
geo1_d = lsim(G_hpf*1/G_geo/(s + 2*pi), id.geo_1, id.t);
geo2_d = lsim(G_hpf*1/G_geo/(s + 2*pi), id.geo_2, id.t);
figure;
hold on;
plot(id.t, id.d);
plot(id.t, acc1_d);
plot(id.t, acc2_d);
plot(id.t, geo1_d);
plot(id.t, geo2_d);
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');
% Fusion
wc = 2*pi*200;
G_hpf = (s/wc)/(1 + s/wc);
G_lpf = 1/(1 + s/wc);
ss_d = lsim(G_hpf, acc1_d, id.t) + lsim(G_lpf, geo1_d, id.t);
figure;
hold on;
plot(id.t, id.d);
plot(id.t, ss_d);
hold off;
xlabel('Time [s]'); ylabel('Displacement [m]');