2 \documentclass[10pt, conference, compsocconf]{IEEEtran}
5 \usepackage[utf8]{inputenc}
6 %\usepackage[cyr]{aeguill}
7 %\usepackage{pstricks,pst-node,pst-text,pst-3d}
16 \usepackage{subfigure}
21 \usepackage[ruled,lined,linesnumbered]{algorithm2e}
23 %%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
24 \newcommand{\noun}[1]{\textsc{#1}}
26 \newcommand{\tab}{\ \ \ }
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}
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}
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}
69 {\it keywords}: FPGA, cantilever, interferometry.
72 \section{Introduction}
74 Cantilevers are used inside atomic force microscope (AFM) which provides high
75 resolution images of surfaces. Several technics have been used to measure the
76 displacement of cantilevers in litterature. For example, it is possible to
77 determine accurately the deflection with different mechanisms.
78 In~\cite{CantiPiezzo01}, authors used piezoresistor integrated into the
79 cantilever. Nevertheless this approach suffers from the complexity of the
80 microfabrication process needed to implement the sensor in the cantilever.
81 In~\cite{CantiCapacitive03}, authors have presented an cantilever mechanism
82 based on capacitive sensing. This kind of technic also involves to instrument
83 the cantiliver which result in a complex fabrication process.
85 In this paper our attention is focused on a method based on interferometry to
86 measure cantilevers' displacements. In this method cantilevers are illuminated
87 by an optic source. The interferometry produces fringes on each cantilevers
88 which enables to compute the cantilever displacement. In order to analyze the
89 fringes a high speed camera is used. Images need to be processed quickly and
90 then a estimation method is required to determine the displacement of each
91 cantilever. In~\cite{AFMCSEM11}, the authors have used an algorithm based on
92 spline to estimate the cantilevers' positions.
94 The overall process gives
95 accurate results but all the computation are performed on a standard computer
96 using labview. Consequently, the main drawback of this implementation is that
97 the computer is a bootleneck in the overall process. In this paper we propose to
98 use a method based on least square and to implement all the computation on a
101 The remainder of the paper is organized as follows. Section~\ref{sec:measure}
102 describes more precisely the measurement process. Our solution based on the
103 least square method and the implementation on FPGA is presented in
104 Section~\ref{sec:solus}. Experimentations are described in
105 Section~\ref{sec:results}. Finally a conclusion and some perspectives are
110 %% quelques ref commentées sur les calculs basés sur l'interférométrie
112 \section{Measurement principles}
122 \subsection{Architecture}
124 %% description de l'architecture générale de l'acquisition d'images
125 %% avec au milieu une unité de traitement dont on ne précise pas ce
128 In order to develop simple, cost effective and user-friendly cantilever arrays,
129 authors of ~\cite{AFMCSEM11} have developped a system based of
130 interferometry. In opposition to other optical based systems, using a laser beam
131 deflection scheme and sentitive to the angular displacement of the cantilever,
132 interferometry is sensitive to the optical path difference induced by the
133 vertical displacement of the cantilever.
135 The system build by authors of~\cite{AFMCSEM11} has been developped based on a
136 Linnick interferomter~\cite{Sinclair:05}. It is illustrated in Figure~\ref{fig:AFM}. A
137 laser beam is first split (by the splitter) into a reference beam and a sample
138 beam that reachs the cantilever array. In order to be able to move the
139 cantilever array, it is mounted on a translation and rotational stage with five
140 degrees of freedom. The optical system is also fixed to the stage. Thus, the
141 cantilever array is centered in the optical system which can be adjusted
142 accurately. The beam illuminates the array by a microscope objective and the
143 light reflects on the cantilevers. Likewise the reference beam reflects on a
144 movable mirror. A CMOS camera chip records the reference and sample beams which
145 are recombined in the beam splitter and the interferogram. At the beginning of
146 each experiment, the movable mirror is fitted manually in order to align the
147 interferometric fringes approximately parallel to the cantilevers. When
148 cantilevers move due to the surface, the bending of cantilevers produce
149 movements in the fringes that can be detected with the CMOS camera. Finally the
150 fringes need to be analyzed. In~\cite{AFMCSEM11}, the authors used a LabView
151 program to compute the cantilevers' movements from the fringes.
155 \includegraphics[width=\columnwidth]{AFM}
157 \caption{schema of the AFM}
162 %% image tirée des expériences.
164 \subsection{Cantilever deflection estimation}
167 As shown on image \ref{img:img-xp}, each cantilever is covered by
168 interferometric fringes. The fringes will distort when cantilevers are
169 deflected. Estimating the deflection is done by computing this
170 distortion. For that, (ref A. Meister + M Favre) proposed a method
171 based on computing the phase of the fringes, at the base of each
172 cantilever, near the tip, and on the base of the array. They assume
173 that a linear relation binds these phases, which can be use to
174 "unwrap" the phase at the tip and to determine the deflection.\\
176 More precisely, segment of pixels are extracted from images taken by a
177 high-speed camera. These segments are large enough to cover several
178 interferometric fringes and are placed at the base and near the tip of
179 the cantilevers. They are called base profile and tip profile in the
180 following. Furthermore, a reference profile is taken on the base of
181 the cantilever array.
183 The pixels intensity $I$ (in gray level) of each profile is modelized by :
187 I(x) = ax+b+A.cos(2\pi f.x + \theta)
190 where $x$ is the position of a pixel in its associated segment.
192 The global method consists in two main sequences. The first one aims
193 to determin the frequency $f$ of each profile with an algorithm based
194 on spline interpolation (see section \ref{algo-spline}). It also
195 computes the coefficient used for unwrapping the phase. The second one
196 is the acquisition loop, while which images are taken at regular time
197 steps. For each image, the phase $\theta$ of all profiles is computed
198 to obtain, after unwrapping, the deflection of cantilevers.
200 \subsection{Design goals}
203 If we put aside some hardware issues like the speed of the link
204 between the camera and the computation unit, the time to deserialize
205 pixels and to store them in memory, ... the phase computation is
206 obviously the bottle-neck of the whole process. For example, if we
207 consider the camera actually in use, an exposition time of 2.5ms for
208 $1024\times 1204$ pixels seems the minimum that can be reached. For a
209 $10\times 10$ cantilever array, if we neglect the time to extract
210 pixels, it implies that computing the deflection of a single
211 cantilever should take less than 25$\mu$s, thus 12.5$\mu$s by phase.\\
213 In fact, this timing is a very hard constraint. Let consider a very
214 small programm that initializes twenty million of doubles in memory
215 and then does 1000000 cumulated sums on 20 contiguous values
216 (experimental profiles have about this size). On an intel Core 2 Duo
217 E6650 at 2.33GHz, this program reaches an average of 155Mflops. It
218 implies that the phase computation algorithm should not take more than
219 $240\times 12.5 = 1937$ floating operations. For integers, it gives
222 %% to be continued ...
224 %% � faire : timing de l'algo spline en C avec atan et tout le bordel.
229 \section{Proposed solution}
233 \subsection{FPGA constraints}
235 A field-programmable gate array (FPGA) is an integrated circuit designed to be
236 configured by the customer. A hardware description language (HDL) is used to
237 configure a FPGA. FGPAs are composed of programmable logic components, called
238 logic blocks. These blocks can be configured to perform simple (AND, XOR, ...)
239 or complex combinational functions. Logic blocks are interconnected by
240 reconfigurable links. Modern FPGAs contains memory elements and multipliers
241 which enables to simplify the design and increase the speed. As the most complex
242 operation operation on FGPAs is the multiplier, design of FGPAs should not used
243 complex operations. For example, a divider is not an available operation and it
244 should be programmed using simple components.
246 FGPAs programming is very different from classic processors programming. When
247 logic block are programmed and linked to performed an operation, they cannot be
248 reused anymore. FPGA are cadenced slowly than classic processors but they can
249 performed pipelined as well as pipelined operations. A pipeline provides a way
250 manipulate data quickly since at each clock top to handle a new data. However,
251 using a pipeline consomes more logics and components since they are not
252 reusable, nevertheless it is probably the most efficient technique on FPGA.
253 Parallel operations can be used in order to manipulate several data
254 simultaneously. When it is possible, using a pipeline is a good solution to
255 manipulate new data at each clock top and using parallelism to handle
256 simultaneously several data streams.
258 %% contraintes imposées par le FPGA : algo pipeline/parallele, pas d'op math complexe, ...
261 \subsection{Considered algorithms}
263 Two solutions have been studied to achieve phase computation. The
264 original one, proposed by A. Meister and M. Favre, is based on
265 interpolation by splines. It allows to compute frequency and
266 phase. The second one, detailed in this article, is based on a
267 classical least square method but suppose that frequency is already
270 \subsubsection{Spline algorithm}
271 \label{sec:algo-spline}
272 Let consider a profile $P$, that is a segment of $M$ pixels with an
273 intensity in gray levels. Let call $I(x)$ the intensity of profile in $x
276 At first, only $M$ values of $I$ are known, for $x = 0, 1,
277 \ldots,M-1$. A normalisation allows to scale known intensities into
278 $[-1,1]$. We compute splines that fit at best these normalised
279 intensities. Splines are used to interpolate $N = k\times M$ points
280 (typically $k=3$ is sufficient), within $[0,M[$. Let call $x^s$ the
281 coordinates of these $N$ points and $I^s$ their intensities.
283 In order to have the frequency, the mean line $a.x+b$ (see equation \ref{equ:profile}) of $I^s$ is
284 computed. Finding intersections of $I^s$ and this line allow to obtain
285 the period thus the frequency.
287 The phase is computed via the equation :
289 \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]
292 Two things can be noticed. Firstly, the frequency could also be
293 obtained using the derivates of spline equations, which only implies
294 to solve quadratic equations. Secondly, frequency of each profile is
295 computed a single time, before the acquisition loop. Thus, $sin(2\pi f
296 x^s_i)$ and $cos(2\pi f x^s_i)$ could also be computed before the loop, which leads to a
297 much faster computation of $\theta$.
299 \subsubsection{Least square algorithm}
301 Assuming that we compute the phase during the acquisition loop,
302 equation \ref{equ:profile} has only 4 parameters :$a, b, A$, and
303 $\theta$, $f$ and $x$ being already known. Since $I$ is non-linear, a
304 least square method based an Gauss-newton algorithm must be used to
305 determine these four parameters. Since it is an iterative process
306 ending with a convergence criterion, it is obvious that it is not
307 particularly adapted to our design goals.
309 Fortunatly, it is quite simple to reduce the number of parameters to
310 only $\theta$. Let $x^p$ be the coordinates of pixels in a segment of
311 size $M$. Thus, $x^p = 0, 1, \ldots, M-1$. Let $I(x^p)$ be their
312 intensity. Firstly, we "remove" the slope by computing :
314 \[I^{corr}(x^p) = I(x^p) - a.x^p - b\]
316 Since linear equation coefficients are searched, a classical least
317 square method can be used to determine $a$ and $b$ :
319 \[a = \frac{covar(x^p,I(x^p))}{var(x^p)} \]
321 Assuming an overlined symbol means an average, then :
323 \[b = \overline{I(x^p)} - a.\overline{{x^p}}\]
325 Let $A$ be the amplitude of $I^{corr}$, i.e.
327 \[A = \frac{max(I^{corr}) - min(I^{corr})}{2}\]
329 Then, the least square method to find $\theta$ is reduced to search the minimum of :
331 \[\sum_{i=0}^{M-1} \left[ cos(2\pi f.i + \theta) - \frac{I^{corr}(i)}{A} \right]^2\]
333 It is equivalent to derivate this expression and to solve the following equation :
336 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] \\
337 - 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
340 Several points can be noticed :
342 \item As in the spline method, some parts of this equation can be
343 computed before the acquisition loop. It is the case of sums that do
344 not depend on $\theta$ :
346 \[ \sum_{i=0}^{M-1} sin(4\pi f.i), \sum_{i=0}^{M-1} cos(4\pi f.i) \]
348 \item Lookup tables for $sin(2\pi f.i)$ and $cos(2\pi f.i)$ can also be
351 \item The simplest method to find the good $\theta$ is to discretize
352 $[-\pi,\pi]$ in $N$ steps, and to search which step leads to the
353 result closest to zero. By the way, three other lookup tables can
354 also be computed before the loop :
356 \[ sin \theta, cos \theta, \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] \]
358 \item This search can be very fast using a dichotomous process in $log_2(N)$
362 Finally, the whole summarizes in an algorithm (called LSQ in the following) in two parts, one before and one during the acquisition loop :
364 \caption{LSQ algorithm - before acquisition loop.}
365 \label{alg:lsq-before}
367 $M \leftarrow $ number of pixels of the profile\\
368 I[] $\leftarrow $ intensities of pixels\\
369 $f \leftarrow $ frequency of the profile\\
370 $s4i \leftarrow \sum_{i=0}^{M-1} sin(4\pi f.i)$\\
371 $c4i \leftarrow \sum_{i=0}^{M-1} cos(4\pi f.i)$\\
372 $nb_s \leftarrow $ number of discretization steps of $[-\pi,\pi]$\\
374 \For{$i=0$ to $nb_s $}{
375 $\theta \leftarrow -\pi + 2\pi\times \frac{i}{nb_s}$\\
376 lut\_sin[$i$] $\leftarrow sin \theta$\\
377 lut\_cos[$i$] $\leftarrow cos \theta$\\
378 lut\_A[$i$] $\leftarrow cos 2 \theta \times s4i + sin 2 \theta \times c4i$\\
379 lut\_sinfi[$i$] $\leftarrow sin (2\pi f.i)$\\
380 lut\_cosfi[$i$] $\leftarrow cos (2\pi f.i)$\\
385 \caption{LSQ algorithm - during acquisition loop.}
386 \label{alg:lsq-during}
388 $\bar{x} \leftarrow \frac{M-1}{2}$\\
389 $\bar{y} \leftarrow 0$, $x_{var} \leftarrow 0$, $xy_{covar} \leftarrow 0$\\
390 \For{$i=0$ to $M-1$}{
391 $\bar{y} \leftarrow \bar{y} + $ I[$i$]\\
392 $x_{var} \leftarrow x_{var} + (i-\bar{x})^2$\\
394 $\bar{y} \leftarrow \frac{\bar{y}}{M}$\\
395 \For{$i=0$ to $M-1$}{
396 $xy_{covar} \leftarrow xy_{covar} + (i-\bar{x}) \times (I[i]-\bar{y})$\\
398 $slope \leftarrow \frac{xy_{covar}}{x_{var}}$\\
399 $start \leftarrow y_{moy} - slope\times \bar{x}$\\
400 \For{$i=0$ to $M-1$}{
401 $I[i] \leftarrow I[i] - start - slope\times i$\tcc*[f]{slope removal}\\
404 $I_{max} \leftarrow max_i(I[i])$, $I_{min} \leftarrow min_i(I[i])$\\
405 $amp \leftarrow \frac{I_{max}-I_{min}}{2}$\\
407 $Is \leftarrow 0$, $Ic \leftarrow 0$\\
408 \For{$i=0$ to $M-1$}{
409 $Is \leftarrow Is + I[i]\times $ lut\_sinfi[$i$]\\
410 $Ic \leftarrow Ic + I[i]\times $ lut\_cosfi[$i$]\\
413 $\theta \leftarrow -\pi$\\
414 $val_1 \leftarrow 2\times \left[ Is.\cos(\theta) + Ic.\sin(\theta) \right] - amp\times \left[ c4i.\sin(2\theta) + s4i.\cos(2\theta) \right]$\\
415 \For{$i=1-n_s$ to $n_s$}{
416 $\theta \leftarrow \frac{i.\pi}{n_s}$\\
417 $val_2 \leftarrow 2\times \left[ Is.\cos(\theta) + Ic.\sin(\theta) \right] - amp\times \left[ c4i.\sin(2\theta) + s4i.\cos(2\theta) \right]$\\
419 \lIf{$val_1 < 0$ et $val_2 >= 0$}{
420 $\theta_s \leftarrow \theta - \left[ \frac{val_2}{val_2-val_1}\times \frac{\pi}{n_s} \right]$\\
422 $val_1 \leftarrow val_2$\\
428 \subsubsection{Comparison}
430 \subsection{VHDL design paradigms}
432 \subsection{VHDL implementation}
434 \section{Experimental results}
440 \section{Conclusion and perspectives}
443 \bibliographystyle{plain}
444 \bibliography{biblio}