]> AND Private Git Repository - dmems12.git/blob - dmems12.tex
Logo AND Algorithmique Numérique Distribuée

Private GIT Repository
9ème :
[dmems12.git] / dmems12.tex
1
2 \documentclass[10pt, peerreview, compsocconf]{IEEEtran}
3 %\usepackage{latex8}
4 %\usepackage{times}
5 \usepackage[utf8]{inputenc}
6 %\usepackage[cyr]{aeguill}
7 %\usepackage{pstricks,pst-node,pst-text,pst-3d}
8 %\usepackage{babel}
9 \usepackage{amsmath}
10 \usepackage{url}
11 \usepackage{graphicx}
12 \usepackage{thumbpdf}
13 \usepackage{color}
14 \usepackage{moreverb}
15 \usepackage{commath}
16 \usepackage{subfigure}
17 %\input{psfig.sty}
18 \usepackage{fullpage}
19 \usepackage{fancybox}
20
21 \usepackage[ruled,lined,linesnumbered]{algorithm2e}
22
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
24 \newcommand{\noun}[1]{\textsc{#1}}
25
26 \newcommand{\tab}{\ \ \ }
27
28
29
30 \begin{document}
31
32
33 %% \author{\IEEEauthorblockN{Authors Name/s per 1st Affiliation (Author)}
34 %% \IEEEauthorblockA{line 1 (of Affiliation): dept. name of organization\\
35 %% line 2: name of organization, acronyms acceptable\\
36 %% line 3: City, Country\\
37 %% line 4: Email: name@xyz.com}
38 %% \and
39 %% \IEEEauthorblockN{Authors Name/s per 2nd Affiliation (Author)}
40 %% \IEEEauthorblockA{line 1 (of Affiliation): dept. name of organization\\
41 %% line 2: name of organization, acronyms acceptable\\
42 %% line 3: City, Country\\
43 %% line 4: Email: name@xyz.com}
44 %% }
45
46
47
48 \title{Using FPGAs for high speed and real time cantilever deflection estimation}
49 \author{\IEEEauthorblockN{Raphaël Couturier\IEEEauthorrefmark{1}, Stéphane Domas\IEEEauthorrefmark{1}, Gwenhaël Goavec-Merou\IEEEauthorrefmark{2} and Michel Lenczner\IEEEauthorrefmark{2}}
50 \IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, DISC, University of Franche-Comte, Belfort, France\\
51 \{raphael.couturier,stephane.domas\}@univ-fcomte.fr}
52 \IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Time-Frequency, University of Franche-Comte, Besançon, France\\
53 \{michel.lenczner@utbm.fr,gwenhael.goavec@trabucayre.com}
54 }
55
56
57
58
59
60
61 %\maketitle
62
63 \thispagestyle{empty}
64
65 \begin{abstract}
66
67   
68
69
70 \end{abstract}
71
72 \begin{IEEEkeywords}
73 FPGA, cantilever, interferometry.
74 \end{IEEEkeywords}
75
76
77 \IEEEpeerreviewmaketitle
78
79 \section{Introduction}
80
81 Cantilevers  are  used  inside  atomic  force  microscope (AFM) which  provides  high
82 resolution images of  surfaces.  Several technics have been  used to measure the
83 displacement  of cantilevers  in litterature.   For example,  it is  possible to
84 determine  accurately  the  deflection  with different  mechanisms. 
85 In~\cite{CantiPiezzo01},   authors  used   piezoresistor  integrated   into  the
86 cantilever.   Nevertheless this  approach  suffers from  the  complexity of  the
87 microfabrication  process needed  to  implement the  sensor  in the  cantilever.
88 In~\cite{CantiCapacitive03},  authors  have  presented an  cantilever  mechanism
89 based on  capacitive sensing. This kind  of technic also  involves to instrument
90 the cantiliver which result in a complex fabrication process.
91
92 In this  paper our attention is focused  on a method based  on interferometry to
93 measure cantilevers' displacements.  In  this method cantilevers are illuminated
94 by  an optic  source. The  interferometry produces  fringes on  each cantilevers
95 which enables to  compute the cantilever displacement.  In  order to analyze the
96 fringes a  high speed camera  is used. Images  need to be processed  quickly and
97 then  a estimation  method is  required to  determine the  displacement  of each
98 cantilever.  In~\cite{AFMCSEM11},  the authors have  used an algorithm  based on
99 spline to estimate the cantilevers' positions.
100
101    The overall  process gives
102 accurate results  but all the computation  are performed on  a standard computer
103 using labview.  Consequently,  the main drawback of this  implementation is that
104 the computer is a bootleneck in the overall process. In this paper we propose to
105 use a  method based on least  square and to  implement all the computation  on a
106 FGPA.
107
108 The remainder  of the paper  is organized as  follows. Section~\ref{sec:measure}
109 describes  more precisely  the measurement  process. Our  solution based  on the
110 least  square   method  and   the  implementation  on   FPGA  is   presented  in
111 Section~\ref{sec:solus}.       Experimentations      are       described      in
112 Section~\ref{sec:results}.  Finally  a  conclusion  and  some  perspectives  are
113 presented.
114
115
116
117 %% quelques ref commentées sur les calculs basés sur l'interférométrie
118
119 \section{Measurement principles}
120 \label{sec:measure}
121
122
123
124
125
126
127
128
129 \subsection{Architecture}
130 \label{sec:archi}
131 %% description de l'architecture générale de l'acquisition d'images
132 %% avec au milieu une unité de traitement dont on ne précise pas ce
133 %% qu'elle est.
134
135 In order to develop simple,  cost effective and user-friendly cantilever arrays,
136 authors   of    ~\cite{AFMCSEM11}   have   developped   a    system   based   of
137 interferometry. In opposition to other optical based systems, using a laser beam
138 deflection scheme and  sentitive to the angular displacement  of the cantilever,
139 interferometry  is sensitive  to  the  optical path  difference  induced by  the
140 vertical displacement of the cantilever.
141
142 The system build  by authors of~\cite{AFMCSEM11} has been  developped based on a
143 Linnick     interferomter~\cite{Sinclair:05}.    It     is     illustrated    in
144 Figure~\ref{fig:AFM}.  A  laser diode  is first split  (by the splitter)  into a
145 reference beam and a sample beam  that reachs the cantilever array.  In order to
146 be  able to  move  the cantilever  array, it  is  mounted on  a translation  and
147 rotational hexapod  stage with  five degrees of  freedom. The optical  system is
148 also fixed to the stage.  Thus,  the cantilever array is centered in the optical
149 system which  can be adjusted accurately.   The beam illuminates the  array by a
150 microscope objective  and the  light reflects on  the cantilevers.  Likewise the
151 reference beam  reflects on a  movable mirror.  A  CMOS camera chip  records the
152 reference and  sample beams which  are recombined in  the beam splitter  and the
153 interferogram.   At the  beginning of  each  experiment, the  movable mirror  is
154 fitted  manually in  order to  align the  interferometric  fringes approximately
155 parallel  to the cantilevers.   When cantilevers  move due  to the  surface, the
156 bending of  cantilevers produce  movements in the  fringes that can  be detected
157 with    the    CMOS    camera.     Finally    the    fringes    need    to    be
158 analyzed. In~\cite{AFMCSEM11}, the authors used a LabView program to compute the
159 cantilevers' movements from the fringes.
160
161 \begin{figure}    
162 \begin{center}
163 \includegraphics[width=\columnwidth]{AFM}
164 \end{center}
165 \caption{schema of the AFM}
166 \label{fig:AFM}   
167 \end{figure}
168
169
170 %% image tirée des expériences.
171
172 \subsection{Cantilever deflection estimation}
173 \label{sec:deflest}
174
175 As shown on image \ref{img:img-xp}, each cantilever is covered by
176 interferometric fringes. The fringes will distort when cantilevers are
177 deflected. Estimating the deflection is done by computing this
178 distortion. For that, (ref A. Meister + M Favre) proposed a method
179 based on computing the phase of the fringes, at the base of each
180 cantilever, near the tip, and on the base of the array. They assume
181 that a linear relation binds these phases, which can be use to
182 "unwrap" the phase at the tip and to determine the deflection.\\
183
184 More precisely, segment of pixels are extracted from images taken by a
185 high-speed camera. These segments are large enough to cover several
186 interferometric fringes and are placed at the base and near the tip of
187 the cantilevers. They are called base profile and tip profile in the
188 following. Furthermore, a reference profile is taken on the base of
189 the cantilever array.
190
191 The pixels intensity $I$ (in gray level) of each profile is modelized by :
192
193 \begin{equation}
194 \label{equ:profile}
195 I(x) = ax+b+A.cos(2\pi f.x + \theta)
196 \end{equation}
197
198 where $x$ is the position of a pixel in its associated segment.
199
200 The global method consists in two main sequences. The first one aims
201 to determin the frequency $f$ of each profile with an algorithm based
202 on spline interpolation (see section \ref{algo-spline}). It also
203 computes the coefficient used for unwrapping the phase. The second one
204 is the acquisition loop, while which images are taken at regular time
205 steps. For each image, the phase $\theta$ of all profiles is computed
206 to obtain, after unwrapping, the deflection of
207 cantilevers. Originally, this computation was also done with an
208 algorithm based on spline. This article proposes a new version based
209 on a least square method.
210
211 \subsection{Design goals}
212 \label{sec:goals}
213
214 The main goal is to implement a computing unit to estimate the
215 deflection of about $10\times10$ cantilevers, faster than the stream of
216 images coming from the camera. The accuracy of results must be close
217 to the maximum precision ever obtained experimentally on the
218 architecture, i.e. 0.3nm. Finally, the latency between an image
219 entering in the unit and the deflections must be as small as possible
220 (NB : future works plan to add some control on the cantilevers).\\
221
222 If we put aside some hardware issues like the speed of the link
223 between the camera and the computation unit, the time to deserialize
224 pixels and to store them in memory, ... the phase computation is
225 obviously the bottle-neck of the whole process. For example, if we
226 consider the camera actually in use, an exposition time of 2.5ms for
227 $1024\times 1204$ pixels seems the minimum that can be reached. For
228 100 cantilevers, if we neglect the time to extract pixels, it implies
229 that computing the deflection of a single
230 cantilever should take less than 25$\mu$s, thus 12.5$\mu$s by phase.\\
231
232 In fact, this timing is a very hard constraint. Let consider a very
233 small programm that initializes twenty million of doubles in memory
234 and then does 1000000 cumulated sums on 20 contiguous values
235 (experimental profiles have about this size). On an intel Core 2 Duo
236 E6650 at 2.33GHz, this program reaches an average of 155Mflops. 
237
238 %%Itimplies that the phase computation algorithm should not take more than
239 %%$155\times 12.5 = 1937$ floating operations. For integers, it gives $3000$ operations. 
240
241 Obviously, some cache effects and optimizations on
242 huge amount of computations can drastically increase these
243 performances : peak efficiency is about 2.5Gflops for the considered
244 CPU. But this is not the case for phase computation that used only few
245 tenth of values.\\
246
247 In order to evaluate the original algorithm, we translated it in C
248 language. As said further, for 20 pixels, it does about 1550
249 operations, thus an estimated execution time of $1550/155
250 =$10$\mu$s. For a more realistic evaluation, we constructed a file of
251 1Mo containing 200 profiles of 20 pixels, equally scattered. This file
252 is equivalent to an image stored in a device file representing the
253 camera. We obtained an average of 10.5$\mu$s by profile (including I/O
254 accesses). It is under are requirements but close to the limit. In
255 case of an occasional load of the system, it could be largely
256 overtaken. A solution would be to use a real-time operating system but
257 another one to search for a more efficient algorithm.
258
259 But the main drawback is the latency of such a solution : since each
260 profile must be treated one after another, the deflection of 100
261 cantilevers takes about $200\times 10.5 = 2.1$ms, which is inadequate
262 for an efficient control. An obvious solution is to parallelize the
263 computations, for example on a GPU. Nevertheless, the cost to transfer
264 profile in GPU memory and to take back results would be prohibitive
265 compared to computation time. It is certainly more efficient to
266 pipeline the computation. For example, supposing that 200 profiles of
267 20 pixels can be pushed sequentially in the pipelined unit cadenced at
268 a 100MHz (i.e. a pixel enters in the unit each 10ns), all profiles
269 would be treated in $200\times 20\times 10.10^{-9} =$ 40$\mu$s plus
270 the latency of the pipeline. This is about 500 times faster than
271 actual results.\\
272
273 For these reasons, an FPGA as the computation unit is the best choice
274 to achieve the required performance. Nevertheless, passing from
275 a C code to a pipelined version in VHDL is not obvious at all. As
276 explained in the next section, it can even be impossible because of
277 some hardware constraints specific to FPGAs.
278
279
280 \section{Proposed solution}
281 \label{sec:solus}
282
283 Project Oscar aims  to provide a hardware and  software architecture to estimate
284 and  control the  deflection of  cantilevers. The  hardware part  consists  in a
285 high-speed camera,  linked on an embedded  board hosting FPGAs. By  the way, the
286 camera output stream can be pushed  directly into the FPGA. The software part is
287 mostly the VHDL  code that deserializes the camera  stream, extracts profile and
288 computes  the deflection. Before  focusing on  our work  to implement  the phase
289 computation, we give some general information about FPGAs and the board we use.
290
291 \subsection{FPGAs}
292
293 A field-programmable gate array (FPGA) is an integrated circuit
294 designed to be configured by the customer. FGPAs are composed of
295 programmable logic components, called configurable logic blocks
296 (CLB). These blocks mainly contains look-up tables (LUT), flip/flops
297 (F/F) and latches, organized in one or more slices connected
298 together. Each CLB can be configured to perform simple (AND, XOR, ...)
299 or complex combinational functions. They are interconnected by
300 reconfigurable links. Modern FPGAs contain memory elements and
301 multipliers which enable to simplify the design and to increase the
302 performance. Nevertheless, all other complex operations, like
303 division, trigonometric functions, $\ldots$ are not available and must
304 be done by configuring a set of CLBs.
305
306 Since this configuration is not obvious at all, it can be done via a
307 framework that synthetize a design written in an hardware description
308 language (HDL), and after, that place and route 
309
310  is used to configure a FPGA.
311 FGPAs programming  is very different  from classic processors  programming. When
312 logic blocks are  programmed and linked to perform an  operation, they cannot be
313 reused anymore.  FPGAs are cadenced more slowly than classic processors but they
314 can perform pipeline  as well as parallel operations. A  pipeline provides a way
315 to  manipulate  data  quickly  since  at   each  clock  top  it  handles  a  new
316 data. However, using  a pipeline consumes more logics  and components since they
317 are not  reusable. Nevertheless it is  probably the most  efficient technique on
318 FPGA.   Parallel operations  can be  used in  order to  manipulate  several data
319 simultaneously. When  it is  possible, using  a pipeline is  a good  solution to
320 manipulate  new  data  at  each  clock  top  and  using  parallelism  to  handle
321 simultaneously several pipelines in order to handle multiple data streams.
322
323 %% parler du VHDL, synthèse et bitstream
324 \subsection{The board}
325
326 The board we use is designed by the Armadeus compagny, under the name
327 SP Vision. It consists in a development board hosting a i.MX27 ARM
328 processor (from Freescale). The board includes all classical
329 connectors : USB, Ethernet, ... A Flash memory contains a Linux kernel
330 that can be launched after booting the board via u-Boot.
331
332 The processor is directly connected to a Spartan3A FPGA (from Xilinx)
333 via its special interface called WEIM. The Spartan3A is itself
334 connected to a Spartan6 FPGA. Thus, it is possible to develop programs
335 that communicate between i.MX and Spartan6, using Spartan3 as a
336 tunnel. By default, the WEIM interface provides a clock signal at
337 100MHz that is connected to dedicated FPGA pins.
338
339 The Spartan6 is an LX100 version. It has 15822 slices, equivalent to
340 101261 logic cells. There are 268 internal block RAM of 18Kbits, and
341 180 dedicated multiply-adders (named DSP48), which is largely enough
342 for our project.
343
344 Some I/O pins of Spartan6 are connected to two $2\times 17$ headers
345 that can be used as user wants. For the project, they will be
346 connected to the interface card of the camera.
347
348 \subsection{Considered algorithms}
349
350 Two solutions have been studied to achieve phase computation. The
351 original one, proposed by A. Meister and M. Favre, is based on
352 interpolation by splines. It allows to compute frequency and
353 phase. The second one, detailed in this article, is based on a
354 classical least square method but suppose that frequency is already
355 known.
356
357 \subsubsection{Spline algorithm}
358 \label{sec:algo-spline}
359 Let consider a profile $P$, that is a segment of $M$ pixels with an
360 intensity in gray levels. Let call $I(x)$ the intensity of profile in $x
361 \in [0,M[$. 
362
363 At first, only $M$ values of $I$ are known, for $x = 0, 1,
364 \ldots,M-1$. A normalisation allows to scale known intensities into
365 $[-1,1]$. We compute splines that fit at best these normalised
366 intensities. Splines are used to interpolate $N = k\times M$ points
367 (typically $k=4$ is sufficient), within $[0,M[$. Let call $x^s$ the
368 coordinates of these $N$ points and $I^s$ their intensities.
369
370 In order to have the frequency, the mean line $a.x+b$ (see equation \ref{equ:profile}) of $I^s$ is
371 computed. Finding intersections of $I^s$ and this line allow to obtain
372 the period thus the frequency.
373
374 The phase is computed via the equation :
375 \begin{equation}
376 \theta = atan \left[ \frac{\sum_{i=0}^{N-1} sin(2\pi f x^s_i) \times I^s(x^s_i)}{\sum_{i=0}^{N-1} cos(2\pi f x^s_i) \times I^s(x^s_i)} \right]
377 \end{equation}
378
379 Two things can be noticed :
380 \begin{itemize}
381 \item the frequency could also be obtained using the derivates of
382   spline equations, which only implies to solve quadratic equations.
383 \item frequency of each profile is computed a single time, before the
384   acquisition loop. Thus, $sin(2\pi f x^s_i)$ and $cos(2\pi f x^s_i)$
385   could also be computed before the loop, which leads to a much faster
386   computation of $\theta$.
387 \end{itemize}
388
389 \subsubsection{Least square algorithm}
390
391 Assuming that we compute the phase during the acquisition loop,
392 equation \ref{equ:profile} has only 4 parameters :$a, b, A$, and
393 $\theta$, $f$ and $x$ being already known. Since $I$ is non-linear, a
394 least square method based an Gauss-newton algorithm must be used to
395 determine these four parameters. Since it is an iterative process
396 ending with a convergence criterion, it is obvious that it is not
397 particularly adapted to our design goals.
398
399 Fortunatly, it is quite simple to reduce the number of parameters to
400 only $\theta$. Let $x^p$ be the coordinates of pixels in a segment of
401 size $M$. Thus, $x^p = 0, 1, \ldots, M-1$. Let $I(x^p)$ be their
402 intensity. Firstly, we "remove" the slope by computing :
403
404 \[I^{corr}(x^p) = I(x^p) - a.x^p - b\]
405
406 Since linear equation coefficients are searched, a classical least
407 square method can be used to determine $a$ and $b$ :
408
409 \[a = \frac{covar(x^p,I(x^p))}{var(x^p)} \]
410
411 Assuming an overlined symbol means an average, then :
412
413 \[b = \overline{I(x^p)} - a.\overline{{x^p}}\]
414
415 Let $A$ be the amplitude of $I^{corr}$, i.e. 
416
417 \[A = \frac{max(I^{corr}) - min(I^{corr})}{2}\]
418
419 Then, the least square method to find $\theta$ is reduced to search the minimum of :
420
421 \[\sum_{i=0}^{M-1} \left[ cos(2\pi f.i + \theta) - \frac{I^{corr}(i)}{A} \right]^2\]
422
423 It is equivalent to derivate this expression and to solve the following equation :
424
425 \begin{eqnarray*}
426 2\left[ cos\theta \sum_{i=0}^{M-1} I^{corr}(i).sin(2\pi f.i) + sin\theta \sum_{i=0}^{M-1} I^{corr}(i).cos(2\pi f.i)\right] \\
427 - A\left[ cos2\theta \sum_{i=0}^{M-1} sin(4\pi f.i) + sin2\theta \sum_{i=0}^{M-1} cos(4\pi f.i)\right]   = 0
428 \end{eqnarray*}
429
430 Several points can be noticed :
431 \begin{itemize}
432 \item As in the spline method, some parts of this equation can be
433   computed before the acquisition loop. It is the case of sums that do
434   not depend on $\theta$ :
435
436 \[ \sum_{i=0}^{M-1} sin(4\pi f.i), \sum_{i=0}^{M-1} cos(4\pi f.i) \] 
437
438 \item Lookup tables for $sin(2\pi f.i)$ and $cos(2\pi f.i)$ can also be
439 computed.
440
441 \item The simplest method to find the good $\theta$ is to discretize
442   $[-\pi,\pi]$ in $nb_s$ steps, and to search which step leads to the
443   result closest to zero. By the way, three other lookup tables can
444   also be computed before the loop :
445
446 \[ sin \theta, cos \theta, \]
447
448 \[ \left[ cos 2\theta \sum_{i=0}^{M-1} sin(4\pi f.i) + sin 2\theta \sum_{i=0}^{M-1} cos(4\pi f.i)\right] \]
449
450 \item This search can be very fast using a dichotomous process in $log_2(nb_s)$ 
451
452 \end{itemize}
453
454 Finally, the whole summarizes in an algorithm (called LSQ in the following) in two parts, one before and one during the acquisition loop :
455 \begin{algorithm}[h]
456 \caption{LSQ algorithm - before acquisition loop.}
457 \label{alg:lsq-before}
458
459    $M \leftarrow $ number of pixels of the profile\\
460    I[] $\leftarrow $ intensities of pixels\\
461    $f \leftarrow $ frequency of the profile\\
462    $s4i \leftarrow \sum_{i=0}^{M-1} sin(4\pi f.i)$\\
463    $c4i \leftarrow \sum_{i=0}^{M-1} cos(4\pi f.i)$\\
464    $nb_s \leftarrow $ number of discretization steps of $[-\pi,\pi]$\\
465
466    \For{$i=0$ to $nb_s $}{
467      $\theta  \leftarrow -\pi + 2\pi\times \frac{i}{nb_s}$\\
468      lut$_s$[$i$] $\leftarrow sin \theta$\\
469      lut$_c$[$i$] $\leftarrow cos \theta$\\
470      lut$_A$[$i$] $\leftarrow cos 2 \theta \times s4i + sin 2 \theta \times c4i$\\
471      lut$_{sfi}$[$i$] $\leftarrow sin (2\pi f.i)$\\
472      lut$_{cfi}$[$i$] $\leftarrow cos (2\pi f.i)$\\
473    }
474 \end{algorithm}
475
476 \begin{algorithm}[ht]
477 \caption{LSQ algorithm - during acquisition loop.}
478 \label{alg:lsq-during}
479
480    $\bar{x} \leftarrow \frac{M-1}{2}$\\
481    $\bar{y} \leftarrow 0$, $x_{var} \leftarrow 0$, $xy_{covar} \leftarrow 0$\\
482    \For{$i=0$ to $M-1$}{
483      $\bar{y} \leftarrow \bar{y} + $ I[$i$]\\
484      $x_{var} \leftarrow x_{var} + (i-\bar{x})^2$\\
485    }
486    $\bar{y} \leftarrow \frac{\bar{y}}{M}$\\
487    \For{$i=0$ to $M-1$}{
488      $xy_{covar} \leftarrow xy_{covar} + (i-\bar{x}) \times (I[i]-\bar{y})$\\
489    }
490    $slope \leftarrow \frac{xy_{covar}}{x_{var}}$\\
491    $start \leftarrow y_{moy} - slope\times \bar{x}$\\
492    \For{$i=0$ to $M-1$}{
493      $I[i] \leftarrow I[i] - start - slope\times i$\\
494    }
495    
496    $I_{max} \leftarrow max_i(I[i])$, $I_{min} \leftarrow min_i(I[i])$\\
497    $amp \leftarrow \frac{I_{max}-I_{min}}{2}$\\
498
499    $Is \leftarrow 0$, $Ic \leftarrow 0$\\
500    \For{$i=0$ to $M-1$}{
501      $Is \leftarrow Is + I[i]\times $ lut$_{sfi}$[$i$]\\
502      $Ic \leftarrow Ic + I[i]\times $ lut$_{cfi}$[$i$]\\
503    }
504
505    $\delta \leftarrow \frac{nb_s}{2}$, $b_l \leftarrow 0$, $b_r \leftarrow \delta$\\
506    $v_l \leftarrow -2.I_s - amp.$lut$_A$[$b_l$]\\
507
508    \While{$\delta >= 1$}{
509
510      $v_r \leftarrow 2.[ Is.$lut$_c$[$b_r$]$ + Ic.$lut$_s$[$b_r$]$ ] - amp.$lut$_A$[$b_r$]\\
511
512      \If{$!(v_l < 0$ and $v_r >= 0)$}{
513        $v_l \leftarrow v_r$ \\
514        $b_l \leftarrow b_r$ \\
515      }
516      $\delta \leftarrow \frac{\delta}{2}$\\
517      $b_r \leftarrow b_l + \delta$\\
518    }
519    \uIf{$!(v_l < 0$ and $v_r >= 0)$}{
520      $v_l \leftarrow v_r$ \\
521      $b_l \leftarrow b_r$ \\
522      $b_r \leftarrow b_l + 1$\\
523      $v_r \leftarrow 2.[ Is.$lut$_c$[$b_r$]$ + Ic.$lut$_s$[$b_r$]$ ] - amp.$lut$_A$[$b_r$]\\
524    }
525    \Else {
526      $b_r \leftarrow b_l + 1$\\
527    }
528
529    \uIf{$ abs(v_l) < v_r$}{
530      $b_{\theta} \leftarrow b_l$ \\
531    }
532    \Else {
533      $b_{\theta} \leftarrow b_r$ \\
534    }
535    $\theta \leftarrow \pi\times \left[\frac{2.b_{ref}}{nb_s}-1\right]$\\
536
537 \end{algorithm}
538
539 \subsubsection{Comparison}
540
541 We compared the two algorithms on the base of three criterions :
542 \begin{itemize}
543 \item precision of results on a cosinus profile, distorted with noise,
544 \item number of operations,
545 \item complexity to implement an FPGA version.
546 \end{itemize}
547
548 For the first item, we produced a matlab version of each algorithm,
549 running with double precision values. The profile was generated for
550 about 34000 different values of period ($\in [3.1, 6.1]$, step = 0.1),
551 phase ($\in [-3.1 , 3.1]$, step = 0.062) and slope ($\in [-2 , 2]$,
552 step = 0.4). For LSQ, $nb_s = 1024$, which leads to a maximal error of
553 $\frac{\pi}{1024}$ on phase computation. Current A. Meister and
554 M. Favre experiments show a ratio of 50 between variation of phase and
555 the deflection of a lever. Thus, the maximal error due to
556 discretization correspond to an error of 0.15nm on the lever
557 deflection, which is smaller than the best precision they achieved,
558 i.e. 0.3nm.
559
560 For each test, we add some noise to the profile : each group of two
561 pixels has its intensity added to a random number picked in $[-N,N]$
562 (NB: it should be noticed that picking a new value for each pixel does
563 not distort enough the profile). The absolute error on the result is
564 evaluated by comparing the difference between the reference and
565 computed phase, out of $2\pi$, expressed in percents. That is : $err =
566 100\times \frac{|\theta_{ref} - \theta_{comp}|}{2\pi}$.
567
568 Table \ref{tab:algo_prec} gives the maximum and average error for the two algorithms and increasing values of $N$.
569
570 \begin{table}[ht]
571   \begin{center}
572     \begin{tabular}{|c|c|c|c|c|}
573       \hline
574   & \multicolumn{2}{c|}{SPL} & \multicolumn{2}{c|}{LSQ} \\ \cline{2-5}
575   noise & max. err. & aver. err. & max. err. & aver. err. \\ \hline
576   0 & 2.46 & 0.58 & 0.49 & 0.1 \\ \hline
577   2.5 & 2.75 & 0.62 & 1.16 & 0.22 \\ \hline
578   5 & 3.77 & 0.72 & 2.47 & 0.41 \\ \hline
579   7.5 & 4.72 & 0.86 & 3.33 & 0.62 \\ \hline
580   10 & 5.62 & 1.03 & 4.29 & 0.81 \\ \hline
581   15 & 7.96 & 1.38 & 6.35 & 1.21 \\ \hline
582   30 & 17.06 & 2.6 & 13.94 & 2.45 \\ \hline
583
584 \end{tabular}
585 \caption{Error (in \%) for cosinus profiles, with noise.}
586 \label{tab:algo_prec}
587 \end{center}
588 \end{table}
589
590 These results show that the two algorithms are very close, with a
591 slight advantage for LSQ. Furthemore, both behave very well against
592 noise. Assuming the experimental ratio of 50 (see above), an error of
593 1 percent on phase correspond to an error of 0.5nm on the lever
594 deflection, which is very close to the best precision.
595
596 Obviously, it is very hard to predict which level of noise will be
597 present in real experiments and how it will distort the
598 profiles. Nevertheless, we can see on figure \ref{fig:noise20} the
599 profile with $N=10$ that leads to the biggest error. It is a bit
600 distorted, with pikes and straight/rounded portions, and relatively
601 close to most of that come from experiments. Figure \ref{fig:noise60}
602 shows a sample of worst profile for $N=30$. It is completly distorted,
603 largely beyond the worst experimental ones. 
604
605 \begin{figure}[ht]
606 \begin{center}
607   \includegraphics[width=9cm]{intens-noise20}
608 \end{center}
609 \caption{Sample of worst profile for N=10}
610 \label{fig:noise20}
611 \end{figure}
612
613 \begin{figure}[ht]
614 \begin{center}
615   \includegraphics[width=9cm]{intens-noise60}
616 \end{center}
617 \caption{Sample of worst profile for N=30}
618 \label{fig:noise60}
619 \end{figure}
620
621 The second criterion is relatively easy to estimate for LSQ and harder
622 for SPL because of $atan$ operation. In both cases, it is proportional
623 to numbers of pixels $M$. For LSQ, it also depends on $nb_s$ and for
624 SPL on $N = k\times M$, i.e. the number of interpolated points. 
625
626 We assume that $M=20$, $nb_s=1024$, $k=4$, all possible parts are
627 already in lookup tables and a limited set of operations (+, -, *, /,
628 $<$, $>$) is taken account. Translating the two algorithms in C code, we
629 obtain about 430 operations for LSQ and 1550 (plus few tenth for
630 $atan$) for SPL. This result is largely in favor of LSQ. Nevertheless,
631 considering the total number of operations is not really pertinent for
632 an FPGA implementation : it mainly depends on the type of operations
633 and their
634 ordering. The final decision is thus driven by the third criterion.\\
635
636 The Spartan 6 used in our architecture has hard constraint : it has no
637 built-in floating point units. Obviously, it is possible to use some
638 existing "black-boxes" for double precision operations. But they have
639 a quite long latency. It is much simpler to exclusively use integers,
640 with a quantization of all double precision values. Obviously, this
641 quantization should not decrease too much the precision of
642 results. Furthermore, it should not lead to a design with a huge
643 latency because of operations that could not complete during a single
644 or few clock cycles. Divisions are in this case and, moreover, they
645 need an varying number of clock cycles to complete. Even
646 multiplications can be a problem : DSP48 take inputs of 18 bits
647 maximum. For larger multiplications, several DSP must be combined,
648 increasing the latency.
649
650 Nevertheless, the hardest constraint does not come from the FPGA
651 characteristics but from the algorithms. Their VHDL implentation will
652 be efficient only if they can be fully (or near) pipelined. By the
653 way, the choice is quickly done : only a small part of SPL can be.
654 Indeed, the computation of spline coefficients implies to solve a
655 tridiagonal system $A.m = b$. Values in $A$ and $b$ can be computed
656 from incoming pixels intensity but after, the back-solve starts with
657 the lastest values, which breaks the pipeline. Moreover, SPL relies on
658 interpolating far more points than profile size. Thus, the end
659 of SPL works on a larger amount of data than the beginning, which
660 also breaks the pipeline.
661
662 LSQ has not this problem : all parts except the dichotomial search
663 work on the same amount of data, i.e. the profile size. Furthermore,
664 LSQ needs less operations than SPL, implying a smaller output
665 latency. Consequently, it is the best candidate for phase
666 computation. Nevertheless, obtaining a fully pipelined version
667 supposes that operations of different parts complete in a single clock
668 cycle. It is the case for simulations but it completely fails when
669 mapping and routing the design on the Spartan6. By the way,
670 extra-latency is generated and there must be idle times between two
671 profiles entering into the pipeline.
672
673 %%Before obtaining the least bitstream, the crucial question is : how to
674 %%translate the C code the LSQ into VHDL ?
675
676
677 %\subsection{VHDL design paradigms}
678
679 \section{Experimental tests}
680
681 \subsection{VHDL implementation}
682
683 % - ecriture d'un code en C avec integer
684 % - calcul de la taille max en bit de chaque variable en fonction de la quantization.
685 % - tests de quantization : équilibre entre précision et contraintes FPGA
686 % - en parallèle : simulink et VHDL à la main
687 %
688 \subsection{Simulation}
689
690 % ghdl + gtkwave
691 % au mieux : une phase tous les 33 cycles, latence de 95 cycles.
692 % mais routage/placement impossible.
693 \subsection{Bitstream creation}
694
695 % pas fait mais prévision d'une sortie tous les 480ns avec une latence de 1120
696
697 \label{sec:results}
698
699
700
701
702 \section{Conclusion and perspectives}
703
704
705 \bibliographystyle{plain}
706 \bibliography{biblio}
707
708 \end{document}