[dmems12.git] / dmems12.tex

\documentclass[10pt, peerreview, compsocconf]{IEEEtran}
%\usepackage{latex8}
%\usepackage{times}
\usepackage[utf8]{inputenc}
%\usepackage[cyr]{aeguill}
%\usepackage{pstricks,pst-node,pst-text,pst-3d}
%\usepackage{babel}
\usepackage{amsmath}
\usepackage{url}
\usepackage{graphicx}
\usepackage{thumbpdf}
\usepackage{color}
\usepackage{moreverb}
\usepackage{commath}
\usepackage{subfigure}
%\input{psfig.sty}
\usepackage{fullpage}
\usepackage{fancybox}

\usepackage[ruled,lined,linesnumbered]{algorithm2e}

%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
\newcommand{\noun}[1]{\textsc{#1}}

\newcommand{\tab}{\ \ \ }


\begin{document}


%% \author{\IEEEauthorblockN{Authors Name/s per 1st Affiliation (Author)}
%% \IEEEauthorblockA{line 1 (of Affiliation): dept. name of organization\\
%% line 2: name of organization, acronyms acceptable\\
%% line 3: City, Country\\
%% line 4: Email: name@xyz.com}
%% \and
%% \IEEEauthorblockN{Authors Name/s per 2nd Affiliation (Author)}
%% \IEEEauthorblockA{line 1 (of Affiliation): dept. name of organization\\
%% line 2: name of organization, acronyms acceptable\\
%% line 3: City, Country\\
%% line 4: Email: name@xyz.com}
%% }



\title{A new approach based on least square methods to estimate in real time cantilevers deflection with a FPGA}
\author{\IEEEauthorblockN{Raphaël Couturier\IEEEauthorrefmark{1}, Stéphane Domas\IEEEauthorrefmark{1}, Gwenhaël Goavec-Merou\IEEEauthorrefmark{2} and Michel Lenczner\IEEEauthorrefmark{2}}
\IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, DISC, University of Franche-Comte, Belfort, France\\
\{raphael.couturier,stephane.domas\}@univ-fcomte.fr}
\IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Time-Frequency, University of Franche-Comte, Besançon, France\\
\{michel.lenczner@utbm.fr,gwenhael.goavec@trabucayre.com}
}






%\maketitle

\thispagestyle{empty}

\begin{abstract}

Atomics force  microscope (AFM) provide  high resolution images of  surfaces. In
this paper, we  focus our attention on an interferometry  method to estimate the
cantilevers deflection.   The initial method  was based on splines  to determine
the phase  of interference fringes,  and thus the deflection.  Computations were
performed on a  PC with LabView.  Here,  we propose a new approach  based on the
least square  methods and its implementation  that we have developed  on a FPGA,
using the pipelining  technique. Simulations and real tests  showed us that this
implementation is  very efficient and should  allow us to  control a cantilevers
array in real time.


\end{abstract}

\begin{IEEEkeywords}
FPGA, cantilever, interferometry.
\end{IEEEkeywords}


\IEEEpeerreviewmaketitle

\section{Introduction}

Cantilevers are  used inside atomic  force microscopes (AFM) which  provide high
resolution images of surfaces.  Several techniques have been used to measure the
displacement  of cantilevers  in literature.   For  example, it  is possible  to
determine    accurately    the    deflection    with    different    mechanisms.
In~\cite{CantiPiezzo01},   authors  used   piezoresistor  integrated   into  the
cantilever.   Nevertheless this  approach  suffers from  the  complexity of  the
microfabrication  process needed  to  implement the  sensor  in the  cantilever.
In~\cite{CantiCapacitive03},  authors  have  presented a  cantilever  mechanism
based on capacitive sensing. This  kind of technique also involves to instrument
the cantilever which results in a complex fabrication process.

In this  paper our attention is focused  on a method based  on interferometry to
measure cantilevers' displacements.  In  this method cantilevers are illuminated
by an optic source. The interferometry produces fringes on each cantilever which
enables  us to compute  the cantilever  displacement.  In  order to  analyze the
fringes a  high speed camera  is used. Images  need to be processed  quickly and
then  an estimation method  is required  to determine  the displacement  of each
cantilever.  In~\cite{AFMCSEM11}, authors have used an algorithm based on spline
to estimate the cantilevers' positions.

The  overall  process  gives  accurate  results but  all  the  computations  are
performed on a standard computer using LabView.  Consequently, the main drawback
of this  implementation is that the computer  is a bottleneck. In  this paper we
propose  to use  a  method  based on  least  squares and  to  implement all  the
computation on a FPGA.

The remainder  of the paper  is organized as  follows. Section~\ref{sec:measure}
describes  more precisely  the measurement  process. Our  solution based  on the
least  square   method  and   the  implementation  on   FPGA  is   presented  in
Section~\ref{sec:solus}.       Experimentations      are       described      in
Section~\ref{sec:results}.  Finally  a  conclusion  and  some  perspectives  are
presented.





\section{Measurement principles}
\label{sec:measure}

\subsection{Architecture}
\label{sec:archi}


In order to develop simple,  cost effective and user-friendly cantilever arrays,
authors   of    ~\cite{AFMCSEM11}   have   developed   a    system   based   on
interferometry. In opposition to other optical based systems, using a laser beam
deflection scheme and  sensitive to the angular displacement  of the cantilever,
interferometry  is sensitive  to  the  optical path  difference  induced by  the
vertical displacement of the cantilever.

The    system   built    by   these    authors   is    based   on    a   Linnick
interferometer~\cite{Sinclair:05}.   It is illustrated  in Figure~\ref{fig:AFM}.
A laser  diode is  first split  (by the splitter)  into a  reference beam  and a
sample beam that reach the cantilever array.   In order to be able to move the
cantilever array,  it is mounted on  a translation and  rotational hexapod stage
with five  degrees of freedom.  The optical system  is also fixed to  the stage.
Thus,  the cantilever  array is  centered  in the  optical system  which can  be
adjusted accurately.  The  beam illuminates the array by  a microscope objective
and the light reflects on the cantilevers.  Likewise the reference beam reflects
on a movable mirror.  A CMOS  camera chip records the reference and sample beams
which  are  recombined in  the  beam splitter  and  the  interferogram.  At  the
beginning of each experiment, the movable  mirror is fitted manually in order to
align  the interferometric  fringes approximately  parallel to  the cantilevers.
When cantilevers  move due  to the surface,  the bending of  cantilevers produce
movements in the fringes that can be detected with the CMOS camera.  Finally the
fringes need to be analyzed. In~\cite{AFMCSEM11}, authors used a LabView program
to compute the cantilevers' deflections from the fringes.

\begin{figure}    
\begin{center}
\includegraphics[width=\columnwidth]{AFM}
\end{center}
\caption{schema of the AFM}
\label{fig:AFM}   
\end{figure}


%% image tirée des expériences.

\subsection{Cantilever deflection estimation}
\label{sec:deflest}

\begin{figure}    
\begin{center}
\includegraphics[width=\columnwidth]{lever-xp}
\end{center}
\caption{Portion of an image picked by the camera}
\label{fig:img-xp}   
\end{figure}

As shown on image \ref{fig:img-xp}, each cantilever is covered by
several interferometric fringes. The fringes will distort when
cantilevers are deflected. Estimating the deflection is done by
computing this distortion. For that, authors of \cite{AFMCSEM11}
proposed a method based on computing the phase of the fringes, at the
base of each cantilever, near the tip, and on the base of the
array. They assume that a linear relation binds these phases, which
can be used to "unwrap" the phase at the tip and to determine the deflection.\\

More precisely, segments of pixels are extracted from images taken by
the camera. These segments are large enough to cover several
interferometric fringes. As said above, they are placed at the base
and near the tip of the cantilevers. They are called base profile and
tip profile in the following. Furthermore, a reference profile is
taken on the base of the cantilever array.

The pixels intensity $I$ (in gray level) of each profile is modelized by:

\begin{equation}
\label{equ:profile}
I(x) = ax+b+A.cos(2\pi f.x + \theta)
\end{equation}

where $x$ is the position of a pixel in its associated segment.

The global method consists in two main sequences. The first one aims
to determine the frequency $f$ of each profile with an algorithm based
on spline interpolation (see section \ref{sec:algo-spline}). It also
computes the coefficient used for unwrapping the phase. The second one
is the acquisition loop, during which images are taken at regular time
steps. For each image, the phase $\theta$ of all profiles is computed
to obtain, after unwrapping, the deflection of
cantilevers. Originally, this computation was also done with an
algorithm based on spline. This article proposes a new version based
on a least square method.

\subsection{Design goals}
\label{sec:goals}

The main goal is to implement a computing unit to estimate the
deflection of about $10\times10$ cantilevers, faster than the stream of
images coming from the camera. The accuracy of results must be close
to the maximum precision ever obtained experimentally on the
architecture, i.e. 0.3nm. Finally, the latency between an image
entering in the unit and the deflections must be as small as possible
(NB: future works plan to add some control on the cantilevers).\\

If we put aside some hardware issues like the speed of the link
between the camera and the computation unit, the time to deserialize
pixels and to store them in memory, ... the phase computation is
obviously the bottle-neck of the whole process. For example, if we
consider the camera actually in use, an exposition time of 2.5ms for
$1024\times 1204$ pixels seems the minimum that can be reached. For
100 cantilevers, if we neglect the time to extract pixels, it implies
that computing the deflection of a single
cantilever should take less than 25$\mu$s, thus 12.5$\mu$s by phase.\\

In fact, this timing is a very hard constraint. Let us consider a very
small program that initializes twenty million of doubles in memory
and then does 1,000,000 cumulated sums on 20 contiguous values
(experimental profiles have about this size). On an intel Core 2 Duo
E6650 at 2.33GHz, this program reaches an average of 155Mflops. 

%%Itimplies that the phase computation algorithm should not take more than
%%$155\times 12.5 = 1937$ floating operations. For integers, it gives $3000$ operations. 

Obviously, some cache effects and optimizations on
huge amount of computations can drastically increase these
performances: peak efficiency is about 2.5Gflops for the considered
CPU. But this is not the case for phase computation that used only a few
tenth of values.\\

In order to evaluate the original algorithm, we translated it in C
language. As stated before, for 20 pixels, it does about 1,550
operations, thus an estimated execution time of $1,550/155
=$10$\mu$s. For a more realistic evaluation, we constructed a file of
1Mo containing 200 profiles of 20 pixels, equally scattered. This file
is equivalent to an image stored in a device file representing the
camera. We obtained an average of 10.5$\mu$s by profile (including I/O
accesses). It is under our requirements but close to the limit. In
case of an occasional load of the system, it could be largely
overtaken. A solution would be to use a real-time operating system but
another one to search for a more efficient algorithm.

But the main drawback is the latency of such a solution: since each
profile must be treated one after another, the deflection of 100
cantilevers takes about $200\times 10.5 = 2.1$ms, which is inadequate
for an efficient control. An obvious solution is to parallelize the
computations, for example on a GPU. Nevertheless, the cost of transferring
profile in GPU memory and of taking back results would be prohibitive
compared to computation time. It is certainly more efficient to
pipeline the computation. For example, supposing that 200 profiles of
20 pixels can be pushed sequentially in the pipelined unit cadenced at
a 100MHz (i.e. a pixel enters in the unit each 10ns), all profiles
would be treated in $200\times 20\times 10.10^{-9} =$ 40$\mu$s plus
the latency of the pipeline. This is about 500 times faster than
actual results.\\

For these reasons, a FPGA as the computation unit is the best choice
to achieve the required performance. Nevertheless, passing from
a C code to a pipelined version in VHDL is not obvious at all. As
explained in the next section, it can even be impossible because of
some hardware constraints specific to FPGAs.


\section{Proposed solution}
\label{sec:solus}

Project Oscar aims  to provide a hardware and  software architecture to estimate
and  control the  deflection of  cantilevers. The  hardware part  consists  in a
high-speed camera, linked  on an embedded board hosting FPGAs.  In this way, the
camera output stream can be pushed  directly into the FPGA. The software part is
mostly the VHDL  code that deserializes the camera  stream, extracts profile and
computes  the deflection. Before  focusing on  our work  to implement  the phase
computation, we give some general information about FPGAs and the board we use.

\subsection{FPGAs}

A field-programmable gate  array (FPGA) is an integrated  circuit designed to be
configured by the customer. FGPAs are composed of programmable logic components,
called  configurable logic  blocks (CLB).  These blocks  mainly  contain look-up
tables  (LUT), flip/flops (F/F)  and latches,  organized in  one or  more slices
connected together. Each CLB can be configured to perform simple (AND, XOR, ...)
or complex  combinational functions.  They are  interconnected by reconfigurable
links.  Modern  FPGAs contain  memory elements and  multipliers which  enable to
simplify the  design and  to increase the  performance. Nevertheless,  all other
complex  operations, like  division, trigonometric  functions, $\ldots$  are not
available  and  must  be  done  by  configuring  a  set  of  CLBs.   Since  this
configuration  is not  obvious at  all, it  can be  done via  a  framework, like
ISE~\cite{ISE}. Such  a software can synthetize  a design written  in a hardware
description language  (HDL), map it onto  CLBs, place/route them  for a specific
FPGA, and finally produce a bitstream  that is used to configure the FPGA. Thus,
from  the developer's  point of  view, the  main difficulty  is to  translate an
algorithm in HDL  code, taking into account FPGA  resources and constraints like
clock signals and I/O values that drive the FPGA.

Indeed, HDL programming is very different from classic languages like
C. A program can be seen as a state-machine, manipulating signals that
evolve from state to state. Moreover, HDL instructions can executed
concurrently. Basic logic operations are used to aggregate signals to
produce new states and assign it to another signal. States are mainly
expressed as arrays of bits. Fortunately, libraries propose some
higher levels representations like signed integers, and arithmetic
operations.

Furthermore, even if FPGAs are cadenced more slowly than classic
processors, they can perform pipeline as well as parallel
operations. A pipeline consists in cutting a process in a sequence of
small tasks, taking the same execution time. It accepts a new data at
each clock top, thus, after a known latency, it also provides a result
at each clock top. However, using a pipeline consumes more logics
since the components of a task are not reusable by another
one. Nevertheless it is probably the most efficient technique on
FPGA. Because of its architecture, it is also very easy to process
several data concurrently. Whenever possible, the best performance
is reached using parallelism to handle simultaneously several
pipelines in order to handle multiple data streams.

\subsection{The board}

The board we use is designed by the Armadeus company, under the name
SP Vision. It consists in a development board hosting a i.MX27 ARM
processor (from Freescale). The board includes all classical
connectors: USB, Ethernet, ... A Flash memory contains a Linux kernel
that can be launched after booting the board via u-Boot.

The processor is directly connected to a Spartan3A FPGA (from Xilinx)
via its special interface called WEIM. The Spartan3A is itself
connected to a Spartan6 FPGA. Thus, it is possible to develop programs
that communicate between i.MX and Spartan6, using Spartan3 as a
tunnel. By default, the WEIM interface provides a clock signal at
100MHz that is connected to dedicated FPGA pins.

The Spartan6 is an LX100 version. It has 15822 slices, each slice
containing 4 LUTs and 8 flip/flops. It is equivalent to 101261 logic
cells. There are 268 internal block RAM of 18Kbits, and 180 dedicated
multiply-adders (named DSP48), which is largely enough for our
project.

Some I/O pins of Spartan6 are connected to two $2\times 17$ headers
that can be used as user wants. For the project, they will be
connected to the interface card of the camera.

\subsection{Considered algorithms}

Two solutions have been studied to achieve phase computation. The
original one, proposed by A. Meister and M. Favre, is based on
interpolation by splines. It allows to compute frequency and
phase. The second one, detailed in this article, is based on a
classical least square method which supposes that the frequency is already
known.

\subsubsection{Spline algorithm (SPL)}
\label{sec:algo-spline}
Let us consider a profile $P$, that is a segment of $M$ pixels with an
intensity in gray levels. Let us call $I(x)$ the intensity of profile in $x
\in [0,M[$. 

At first, only $M$ values of $I$ are known, for $x = 0, 1,
\ldots,M-1$. A normalization allows to scale known intensities into
$[-1,1]$. We compute splines that fit at best these normalized
intensities. Splines are used to interpolate $N = k\times M$ points
(typically $k=4$ is sufficient), within $[0,M[$. Let  $x^s$ be the
coordinates of these $N$ points and $I^s$ their intensities.

In order to have the frequency, the mean line $a.x+b$ (see equation \ref{equ:profile}) of $I^s$ is
computed. Finding intersections of $I^s$ and this line allows us to obtain
the period and thus the frequency.

The phase is computed via the equation:
\begin{equation}
\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]
\end{equation}

Two things can be noticed:
\begin{itemize}
\item the frequency could also be obtained using the derivates of
  spline equations, which only implies to solve quadratic equations.
\item frequency of each profile is computed only once, before the
  acquisition loop. Thus, $sin(2\pi f x^s_i)$ and $cos(2\pi f x^s_i)$
  could also be computed before the loop, which would lead to a much faster
  computation of $\theta$.
\end{itemize}

\subsubsection{Least square algorithm (LSQ)}

Assuming that we compute the phase during the acquisition loop,
equation \ref{equ:profile} has only 4 parameters: $a, b, A$, and
$\theta$, $f$ and $x$ being already known. Since $I$ is non-linear, a
least square method based on a Gauss-newton algorithm can be used to
determine these four parameters. Since it is an iterative process
ending with a convergence criterion, it is obvious that it is not
particularly adapted to our design goals.

Fortunately, it is quite simple to reduce the number of parameters to
only $\theta$. Let $x^p$ be the coordinates of pixels in a segment of
size $M$. Thus, $x^p = 0, 1, \ldots, M-1$. Let $I(x^p)$ be their
intensity. Firstly, we "remove" the slope by computing:

\[I^{corr}(x^p) = I(x^p) - a.x^p - b\]

Since linear equation coefficients are searched, a classical least
square method can be used to determine $a$ and $b$:

\[a = \frac{covar(x^p,I(x^p))}{var(x^p)} \]

Assuming an overlined symbol means an average, then:

\[b = \overline{I(x^p)} - a.\overline{{x^p}}\]

Let $A$ be the amplitude of $I^{corr}$, i.e. 

\[A = \frac{max(I^{corr}) - min(I^{corr})}{2}\]

Then, the least square method to find $\theta$ is reduced to search the minimum of:

\[\sum_{i=0}^{M-1} \left[ cos(2\pi f.i + \theta) - \frac{I^{corr}(i)}{A} \right]^2\]

It is equivalent to derivating this expression and to solving the following equation:


%\begin{eqnarray*}{l}
$$2\left[ cos\theta \sum_{i=0}^{M-1} I^{corr}(i).sin(2\pi f.i) \right.$$
$$\left. + sin\theta \sum_{i=0}^{M-1} I^{corr}(i).cos(2\pi f.i)\right]- $$
$$ 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$$
%\end{eqnarray*}


Several points can be noticed:
\begin{itemize}
\item As in the spline method, some parts of this equation can be
  computed before the acquisition loop. It is the case of sums that do
  not depend on $\theta$:

\[ \sum_{i=0}^{M-1} sin(4\pi f.i), \sum_{i=0}^{M-1} cos(4\pi f.i) \] 

\item Lookup tables for $sin(2\pi f.i)$ and $cos(2\pi f.i)$ can also be
computed.

\item The simplest method to find the good $\theta$ is to discretize
  $[-\pi,\pi]$ in $nb_s$ steps, and to search which step leads to the
  result closest to zero. Hence, three other lookup tables can
  also be computed before the loop:

\[ 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] \]

\item This search can be very fast using a dichotomous process in $log_2(nb_s)$ 

\end{itemize}

Finally, this is synthetized in an algorithm (called LSQ in the following) in two parts, one before and one during the acquisition loop:
\begin{algorithm}[htbp]
\caption{LSQ algorithm - before acquisition loop.}
\label{alg:lsq-before}

   $M \leftarrow $ number of pixels of the profile\\
   I[] $\leftarrow $ intensity of pixels\\
   $f \leftarrow $ frequency of the profile\\
   $s4i \leftarrow \sum_{i=0}^{M-1} sin(4\pi f.i)$\\
   $c4i \leftarrow \sum_{i=0}^{M-1} cos(4\pi f.i)$\\
   $nb_s \leftarrow $ number of discretization steps of $[-\pi,\pi]$\\

   \For{$i=0$ to $nb_s $}{
     $\theta  \leftarrow -\pi + 2\pi\times \frac{i}{nb_s}$\\
     lut$_s$[$i$] $\leftarrow sin \theta$\\
     lut$_c$[$i$] $\leftarrow cos \theta$\\
     lut$_A$[$i$] $\leftarrow cos 2 \theta \times s4i + sin 2 \theta \times c4i$\\
     lut$_{sfi}$[$i$] $\leftarrow sin (2\pi f.i)$\\
     lut$_{cfi}$[$i$] $\leftarrow cos (2\pi f.i)$\\
   }
\end{algorithm}

\begin{algorithm}[htbp]
\caption{LSQ algorithm - during acquisition loop.}
\label{alg:lsq-during}

   $\bar{x} \leftarrow \frac{M-1}{2}$\\
   $\bar{y} \leftarrow 0$, $x_{var} \leftarrow 0$, $xy_{covar} \leftarrow 0$\\
   \For{$i=0$ to $M-1$}{
     $\bar{y} \leftarrow \bar{y} + $ I[$i$]\\
     $x_{var} \leftarrow x_{var} + (i-\bar{x})^2$\\
   }
   $\bar{y} \leftarrow \frac{\bar{y}}{M}$\\
   \For{$i=0$ to $M-1$}{
     $xy_{covar} \leftarrow xy_{covar} + (i-\bar{x}) \times (I[i]-\bar{y})$\\
   }
   $slope \leftarrow \frac{xy_{covar}}{x_{var}}$\\
   $start \leftarrow y_{moy} - slope\times \bar{x}$\\
   \For{$i=0$ to $M-1$}{
     $I[i] \leftarrow I[i] - start - slope\times i$\\
   }
   
   $I_{max} \leftarrow max_i(I[i])$, $I_{min} \leftarrow min_i(I[i])$\\
   $amp \leftarrow \frac{I_{max}-I_{min}}{2}$\\

   $Is \leftarrow 0$, $Ic \leftarrow 0$\\
   \For{$i=0$ to $M-1$}{
     $Is \leftarrow Is + I[i]\times $ lut$_{sfi}$[$i$]\\
     $Ic \leftarrow Ic + I[i]\times $ lut$_{cfi}$[$i$]\\
   }

   $\delta \leftarrow \frac{nb_s}{2}$, $b_l \leftarrow 0$, $b_r \leftarrow \delta$\\
   $v_l \leftarrow -2.I_s - amp.$lut$_A$[$b_l$]\\

   \While{$\delta >= 1$}{

     $v_r \leftarrow 2.[ Is.$lut$_c$[$b_r$]$ + Ic.$lut$_s$[$b_r$]$ ] - amp.$lut$_A$[$b_r$]\\

     \If{$!(v_l < 0$ and $v_r >= 0)$}{
       $v_l \leftarrow v_r$ \\
       $b_l \leftarrow b_r$ \\
     }
     $\delta \leftarrow \frac{\delta}{2}$\\
     $b_r \leftarrow b_l + \delta$\\
   }
   \uIf{$!(v_l < 0$ and $v_r >= 0)$}{
     $v_l \leftarrow v_r$ \\
     $b_l \leftarrow b_r$ \\
     $b_r \leftarrow b_l + 1$\\
     $v_r \leftarrow 2.[ Is.$lut$_c$[$b_r$]$ + Ic.$lut$_s$[$b_r$]$ ] - amp.$lut$_A$[$b_r$]\\
   }
   \Else {
     $b_r \leftarrow b_l + 1$\\
   }

   \uIf{$ abs(v_l) < v_r$}{
     $b_{\theta} \leftarrow b_l$ \\
   }
   \Else {
     $b_{\theta} \leftarrow b_r$ \\
   }
   $\theta \leftarrow \pi\times \left[\frac{2.b_{ref}}{nb_s}-1\right]$\\

\end{algorithm}

\subsubsection{Comparison}

We compared the two algorithms on the base of three criteria:
\begin{itemize}
\item precision of results on a cosines profile, distorted by noise,
\item number of operations,
\item complexity of implementating an FPGA version.
\end{itemize}

For the first item, we produced a matlab version of each algorithm,
running with double precision values. The profile was generated for
about 34000 different values of period ($\in [3.1, 6.1]$, step = 0.1),
phase ($\in [-3.1 , 3.1]$, step = 0.062) and slope ($\in [-2 , 2]$,
step = 0.4). For LSQ, $nb_s = 1024$, which leads to a maximal error of
$\frac{\pi}{1024}$ on phase computation. Current A. Meister and
M. Favre's experiments show a ratio of 50 between the variation of a phase and
the deflection of a lever. Thus, the maximal error due to
discretization corresponds to an error of 0.15nm on the lever
deflection, which is smaller than the best precision they achieved,
i.e. 0.3nm.

For each test, we add some noise to the profile: each group of two
pixels has its intensity added to a random number picked in $[-N,N]$
(NB: it should be noticed that picking a new value for each pixel does
not distort enough the profile). The absolute error on the result is
evaluated by comparing the difference between the reference and
computed phase, out of $2\pi$, expressed in percentage. That is: $err =
100\times \frac{|\theta_{ref} - \theta_{comp}|}{2\pi}$.

Table  \ref{tab:algo_prec}  gives  the   maximum  and  average  error  for  both
algorithms and increasing values of $N$.

\begin{table}[ht]
  \begin{center}
    \begin{tabular}{|c|c|c|c|c|}
      \hline
  & \multicolumn{2}{c|}{SPL} & \multicolumn{2}{c|}{LSQ} \\ \cline{2-5}
  noise & max. err. & aver. err. & max. err. & aver. err. \\ \hline
  0 & 2.46 & 0.58 & 0.49 & 0.1 \\ \hline
  2.5 & 2.75 & 0.62 & 1.16 & 0.22 \\ \hline
  5 & 3.77 & 0.72 & 2.47 & 0.41 \\ \hline
  7.5 & 4.72 & 0.86 & 3.33 & 0.62 \\ \hline
  10 & 5.62 & 1.03 & 4.29 & 0.81 \\ \hline
  15 & 7.96 & 1.38 & 6.35 & 1.21 \\ \hline
  30 & 17.06 & 2.6 & 13.94 & 2.45 \\ \hline

\end{tabular}
\caption{Error (in \%) for cosines profiles, with noise.}
\label{tab:algo_prec}
\end{center}
\end{table}

These results show that the two algorithms are very close, with a
slight advantage for LSQ. Furthermore, both behave very well against
noise. Assuming the experimental ratio of 50 (see above), an error of
1 percent on the phase corresponds to an error of 0.5nm on the lever
deflection, which is very close to the best precision.

Obviously, it is very hard to predict which level of noise will be
present in real experiments and how it will distort the
profiles. Nevertheless, we can see on figure \ref{fig:noise20} the
profile with $N=10$ that leads to the biggest error. It is a bit
distorted, with pikes and straight/rounded portions, and relatively
close to experiments. Figure \ref{fig:noise60}
shows a sample of worst profile for $N=30$. It is completely distorted,
largely beyond the worst experimental ones. 

\begin{figure}[ht]
\begin{center}
  \includegraphics[width=\columnwidth]{intens-noise20}
\end{center}
\caption{Sample of worst profile for N=10}
\label{fig:noise20}
\end{figure}

\begin{figure}[ht]
\begin{center}
  \includegraphics[width=\columnwidth]{intens-noise60}
\end{center}
\caption{Sample of worst profile for N=30}
\label{fig:noise60}
\end{figure}

The second criterion is relatively easy to estimate for LSQ and harder
for SPL because of $atan$ operation. In both cases, it is proportional
to the numbers of pixels $M$. For LSQ, it also depends on $nb_s$ and for
SPL on $N = k\times M$, i.e. the number of interpolated points. 

We assume that $M=20$, $nb_s=1024$, $k=4$, all possible parts are
already in lookup tables and a limited set of operations (+, -, *, /,
$<$, $>$) is taken into account. Translating both algorithms in C code, we
obtain about 430 operations for LSQ and 1,550 (plus a few tenth for
$atan$) for SPL. This result is largely in favor of LSQ. Nevertheless,
considering the total number of operations is not really pertinent for
an FPGA implementation: it mainly depends on the type of operations
and their
ordering. The final decision is thus driven by the third criterion.\\

The Spartan 6 used in our architecture has a hard constraint: it has no built-in
floating  point  units.   Obviously,  it   is  possible  to  use  some  existing
"black-boxes"  for double  precision operations.   But  they have  quite a  long
latency. It is much simpler to  exclusively use integers, with a quantization of
all double  precision values. Obviously,  this quantization should  not decrease
too much the  precision of results. Furthermore, it should not  lead to a design
with  a huge  latency because  of operations  that could  not complete  during a
single or  few clock  cycles. Divisions fall  into that category  and, moreover,
they need a varying number of clock cycles to complete. Even multiplications can
be a problem: a DSP48 takes inputs of 18 bits maximum.  For larger multiplications,
several DSP must be combined, increasing the latency.

Nevertheless, the hardest constraint does not come from the FPGA characteristics
but from  the algorithms.  Their VHDL implementation  will be efficient  only if
they can be fully (or near) pipelined.  Thus, the choice is quickly made: only a
small  part  of  SPL  can  be  pipelined.  Indeed,  the  computation  of  spline
coefficients implies to solve a tridiagonal  system $A.m = b$. Values in $A$ and
$b$ can  be computed  from incoming pixels  intensity but after,  the back-solve
starts with the latest values,  which breaks the pipeline.  Moreover, SPL relies
on interpolating far  more points than profile size. Thus, the  end of SPL works
on a larger amount of data than at the beginning, which also breaks the pipeline.

LSQ has  not this problem: all parts  except the dichotomial search  work on the
same  amount  of  data, i.e.  the  profile  size.  Furthermore, LSQ  needs  less
operations than SPL, implying a  smaller output latency. Consequently, it is the
best candidate for phase  computation. Nevertheless, obtaining a fully pipelined
version supposes that  operations of different parts complete  in a single clock
cycle. It is  the case for simulations but it completely  fails when mapping and
routing the design  on the Spartan6. Thus,  extra-latency is generated and
there must be idle times between two profiles entering into the pipeline.

%%Before obtaining the least bitstream, the crucial question is: how to
%%translate the C code the LSQ into VHDL ?


%\subsection{VHDL design paradigms}

\section{Experimental tests}

%In this section we explain what  we have done yet. Until now, we could not perform
%real experiments  since we just have  received the FGPA  board. Nevertheless, we
%will include real experiments in the final version of this paper.

\subsection{VHDL implementation}

From the LSQ algorithm, we have written a C program that uses only
integer values. We use a very simple quantization by multiplying
double precision values by a power of two, keeping the integer
part. For example, all values stored in lut$_s$, lut$_c$, $\ldots$ are
scaled by 1,024. Since LSQ also computes average, variance, ... to
remove the slope, the result of implied Euclidean divisions may be
relatively wrong. To avoid that, we also scale the pixel intensities
by a power of two. Furthermore, assuming $nb_s$ is fixed, these
divisions have a known denominator. Thus, they can be replaced by
their multiplication/shift counterpart. Finally, all other
multiplications or divisions by a power of two have been replaced by
left or right bit shifts. Thus, the code only contains
additions, subtractions and multiplications of signed integers, which
are perfectly adapted to FGPAs.

As  mentioned above, hardware  constraints have  a great  influence on  the VHDL
implementation. Consequently, we searched the  maximum value of each variable as
a function of the different scale  factors and the size of profiles, which gives
their maximum size in bits. That  size determines the maximum scale factors that
allow to  use the  least possible  RAMs and DSPs.  Actually, we  implemented our
algorithm  with  this  maximum  size  but  current works  study  the  impact  of
quantization on  the results precision  and design complexity. We  have compared
the result of  the LSQ version using integers and doubles  and observed that the
precision of both were similar.

Then we built two versions of VHDL codes: one directly by hand coding
and the other with Matlab using the Simulink HDL coder
feature~\cite{HDLCoder}. Although the approach is completely different
we obtained VHDL codes that are quite comparable. Each approach has
advantages and drawbacks.  Roughly speaking, hand coding provides
beautiful and much better structured code while Simulink enables us to
produce a code faster.  In terms of throughput and latency,
simulations show that the two approaches are close with a slight
advantage for hand coding.  We hope that real experiments will confirm
that.

\subsection{Simulation}

Before experimental tests on the board, we simulated our two VHDL
codes with GHDL and GTKWave (two free tools with linux). For that, we
built a testbench based on profiles taken from experimentations and
compared the results to values given by the SPL algorithm. Both
versions lead to correct results.

Our first codes were highly optimized : the pipeline could compute a
new phase each 33 cycles and its latency was equal to 95 cycles. Since
the Spartan6 is clocked at 100MHz, it implies that estimating the
deflection of 100 cantilevers would take about $(95 + 200\times 33).10
= 66.95\mu$s, i.e. nearly 15,000 estimations by second.

\subsection{Bitstream creation}

In order to test our code on the SP Vision board, the design was
extended with a component that keeps profiles in RAM, flushes them in
the phase computation component and stores its output in another
RAM. We also added a wishbone : a component that can "drive" signals
to communicate between i.MX and other components. It is mainly used
to start to flush profiles and to retrieve the computed phases in RAM.

Unfortunately, the first  designs could not be placed and route  with ISE on the
Spartan6 with  a 100MHz clock. The  main problems came from  routing values from
RAMs to DSPs and obtaining a result  under 10ns. So, we needed to decompose some
parts of  the pipeline, which  adds some cycles.  For example, some  delays have
been introduced between  RAMs output and DSPs. Finally,  we obtained a bitstream
that has a latency  of 112 cycles and computes a new  phase every 40 cycles. For
100 cantilevers, it takes $(112 + 200\times 40).10 = 81.12\mu$s to compute their
deflection.

This bitstream has been successfully tested on the board.



\label{sec:results}




\section{Conclusion and perspectives}
In this paper we have presented a new method to estimate the
cantilevers deflection in an AFM.  This method is based on least
square methods.  We have used quantization to produce an algorithm
based exclusively on integer values, which is adapted to a FPGA
implementation. We obtained a precision on results similar to the
initial version based on splines.  Our solution has been implemented
with a pipeline technique.  Consequently, it enables to handle a new
profile image very quickly. Currently we have performed simulations
and real tests on a Spartan6 FPGA.

In future work,  we plan to study  the quantization. Then we want  to couple our
algorithm with a high speed camera and we plan to control the whole AFM system.

\bibliographystyle{plain}
\bibliography{biblio}

\end{document}