% fusionner max rejection a surface donnee v.s minimiser surface a rejection donnee
  % demontrer comment la quantification rejette du bruit vers les hautes frequences => 6 dB de

  %    rejection par bit et perte si moins de bits que rejection/6
  % developper programme lineaire en incluant le decalage de bits

  % insister que avant on etait synthetisable mais pas implementable, alors que maintenant on
  % implemente et on demontre que ca tourne

  %   gwen : pourquoi le FIR est desormais implementable et ne l'etait pas meme sur zedboard->new FIR ?
  % Gwen : peut-on faire un vrai banc de bruit de phase avec ce FIR, ie ajouter ADC, NCO et mixer
  %        (zedboard ou redpit)

  % label schema : verifier que "argumenter de la cascade de FIR" est fait

  \documentclass[a4paper,conference]{IEEEtran/IEEEtran}
  \usepackage{graphicx,color,hyperref}
  \usepackage{amsfonts}
  \usepackage{amsthm}
  \usepackage{amssymb}
  \usepackage{amsmath}
  \usepackage{algorithm2e}
  \usepackage{url,balance}
  \usepackage[normalem]{ulem}

  \usepackage{tikz}
  \usetikzlibrary{positioning,fit}
  \usepackage{multirow}
  \usepackage{scalefnt}

  % correct bad hyphenation here
  \hyphenation{op-tical net-works semi-conduc-tor}
  \textheight=26cm
  \setlength{\footskip}{30pt}
  \pagenumbering{gobble}
  \begin{document}
  \title{Filter optimization for real time digital processing of radiofrequency signals: application
  to oscillator metrology}
  
  \author{\IEEEauthorblockN{A. Hugeat\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, J. Bernard\IEEEauthorrefmark{2},
  G. Goavec-M\'erou\IEEEauthorrefmark{1},
  P.-Y. Bourgeois\IEEEauthorrefmark{1}, J.-M. Friedt\IEEEauthorrefmark{1}}
  \IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, Time \& Frequency department, Besan\c con, France }
  \IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Computer Science department DISC, Besan\c con, France \\
  Email: \{pyb2,jmfriedt\}@femto-st.fr}
  }
  \maketitle
  \thispagestyle{plain}
  \pagestyle{plain}
  
  ewtheorem{definition}{Definition}
  
  \begin{abstract}
  Software Defined Radio (SDR) provides stability, flexibility and reconfigurability to
  radiofrequency signal processing. Applied to oscillator characterization in the context
  of ultrastable clocks, stringent filtering requirements are defined by spurious signal or
  noise rejection needs. Since real time radiofrequency processing must be performed in a
  Field Programmable Array to meet timing constraints, we investigate optimization strategies
  to design filters meeting rejection characteristics while limiting the hardware resources

  required and keeping timing constraints within the targeted measurement bandwidths. The
  presented technique is applicable to scheduling any sequence of processing blocks characterized
  by a throughput, resource occupation and performance tabulated as a function of configuration
  characateristics, as is the case for filters with their coefficients and resolution yielding
  rejection and number of multipliers.

  \end{abstract}
  
  \begin{IEEEkeywords}
  Software Defined Radio, Mixed-Integer Linear Programming, Finite Impulse Response filter
  \end{IEEEkeywords}
  
  \section{Digital signal processing of ultrastable clock signals}
  
  Analog oscillator phase noise characteristics are classically performed by downconverting
  the radiofrequency signal using a saturated mixer to bring the radiofrequency signal to baseband,
  followed by a Fourier analysis of the beat signal to analyze phase fluctuations close to carrier. In
  a fully digital approach, the radiofrequency signal is digitized and numerically downconverted by
  multiplying the samples with a local numerically controlled oscillator (Fig. \ref{schema}) \cite{rsi}.
  
  \begin{figure}[h!tb]
  \begin{center}
  \includegraphics[width=.8\linewidth]{images/schema}
  \end{center}
  \caption{Fully digital oscillator phase noise characterization: the Device Under Test
  (DUT) signal is sampled by the radiofrequency grade Analog to Digital Converter (ADC) and
  downconverted by mixing with a Numerically Controlled Oscillator (NCO). Unwanted signals
  and noise aliases are rejected by a Low Pass Filter (LPF) implemented as a cascade of Finite
  Impulse Response (FIR) filters. The signal is then decimated before a Fourier analysis displays
  the spectral characteristics of the phase fluctuations.}
  \label{schema}
  \end{figure}
  
  As with the analog mixer,
  the non-linear behavior of the downconverter introduces noise or spurious signal aliasing as
  well as the generation of the frequency sum signal in addition to the frequency difference.
  These unwanted spectral characteristics must be rejected before decimating the data stream
  for the phase noise spectral characterization \cite{andrich2018high}. The characteristics introduced between the
  downconverter
  and the decimation processing blocks are core characteristics of an oscillator characterization
  system, and must reject out-of-band signals below the targeted phase noise -- typically in the
  sub -170~dBc/Hz for ultrastable oscillator we aim at characterizing. The filter blocks will
  use most resources of the Field Programmable Gate Array (FPGA) used to process the radiofrequency
  datastream: optimizing the performance of the filter while reducing the needed resources is
  hence tackled in a systematic approach using optimization techniques. Most significantly, we
  tackle the issue by attempting to cascade multiple Finite Impulse Response (FIR) filters with
  tunable number of coefficients and tunable number of bits representing the coefficients and the
  data being processed.
  
  \section{Finite impulse response filter}

  We select FIR filters for their unconditional stability and ease of design. A FIR filter is defined

  by a set of weights $b_k$ applied to the inputs $x_k$ through a convolution to generate the
  outputs $y_k$

  \begin{align}
      y_n=\sum_{k=0}^N b_k x_{n-k}
      \label{eq:fir_equation}
  \end{align}

  
  As opposed to an implementation on a general purpose processor in which word size is defined by the

  processor architecture, implementing such a filter on an FPGA offers more degrees of freedom since

  not only the coefficient values and number of taps must be defined, but also the number of bits
  defining the coefficients and the sample size. For this reason, and because we consider pipeline
  processing (as opposed to First-In, First-Out FIFO memory batch processing) of radiofrequency
  signals, High Level Synthesis (HLS) languages \cite{kasbah2008multigrid} are not considered but

  the problem is tackled at the Very-high-speed-integrated-circuit Hardware Description Language

  (VHDL) level.

  Since latency is not an issue in a openloop phase noise characterization instrument, the large
  numbre of taps in the FIR, as opposed to the shorter Infinite Impulse Response (IIR) filter,
  is not considered as an issue as would be in a closed loop system.
  
  The coefficients are classically expressed as floating point values. However, this binary
  number representation is not efficient for fast arithmetic computation by an FPGA. Instead,
  we select to quantify these floating point values into integer values. This quantization
  will result in some precision loss.

  \begin{figure}[h!tb]

  \includegraphics[width=\linewidth]{images/zero_values}

  \caption{Impact of the quantization resolution of the coefficients: the quantization is
  set to 6~bits -- with the horizontal black lines indicating $\pm$1 least significant bit -- setting
  the 30~first and 30~last coefficients out of the initial 128~band-pass
  filter coefficients to 0 (red dots).}
  \label{float_vs_int}
  \end{figure}
  
  The tradeoff between quantization resolution and number of coefficients when considering
  integer operations is not trivial. As an illustration of the issue related to the
  relation between number of fiter taps and quantization, Fig. \ref{float_vs_int} exhibits
  a 128-coefficient FIR bandpass filter designed using floating point numbers (blue). Upon
  quantization on 6~bit integers, 60 of the 128~coefficients in the beginning and end of the
  taps become null, making the large number of coefficients irrelevant and allowing to save
  processing resource by shrinking the filter length. This tradeoff aimed at minimizing resources
  to reach a given rejection level, or maximizing out of band rejection for a given computational
  resource, will drive the investigation on cascading filters designed with varying tap resolution
  and tap length, as will be shown in the next section. Indeed, our development strategy closely
  follows the skeleton approach \cite{crookes1998environment, crookes2000design, benkrid2002towards}
  in which basic blocks are defined and characterized before being assembled \cite{hide}
  in a complete processing chain. In our case, assembling the filter blocks is a simpler block
  combination process since we assume a single value to be processed and a single value to be
  generated at each clock cycle. The FIR filters will not be considered to decimate in the
  current implementation: the decimation is assumed to be located after the FIR cascade at the
  moment.

  \section{Methodology description}

  Our objective is to develop a new methodology applicable to any Digital Signal Processing (DSP)
  chain obtained by assembling basic processing blocks, with hardware and manufacturer independence.

  Achieving such a target requires defining an abstract model to represent some basic properties
  of DSP blocks such as perfomance (i.e. rejection or ripples in the bandpass for filters) and
  resource occupation. These abstract properties, not necessarily related to the detailed hardware
  implementation of a given platform, will feed a scheduler solver aimed at assembling the optimum
  target, whether in terms of maximizing performance for a given arbitrary resource occupation, or
  minimizing resource occupation for a given perfomance. In our approach, the solution of the
  solver is then synthesized using the dedicated tool provided by each platform manufacturer
  to assess the validity of our abstract resource occupation indicator, and the result of running
  the DSP chain on the FPGA allows for assessing the performance of the scheduler. We emphasize
  that all solutions found by the solver are synthesized and executed on hardware at the end
  of the analysis.
  
  In this demonstration , we focus on only two operations: filtering and shifting the number of
  bits needed to represent the data along the processing chain.

  We have chosen these basic operations because shifting and the filtering have already been studied

  in the literature \cite{lim_1996, lim_1988, young_1992, smith_1998} providing a framework for

  assessing our results. Furthermore, filtering is a core step in any radiofrequency frontend
  requiring pipelined processing at full bandwidth for the earliest steps, including for

  time and frequency transfer or characterization \cite{carolina1,carolina2,rsi}.
  
  Addressing only two operations allows for demonstrating the methodology but should not be
  considered as a limitation of the framework which can be extended to assembling any number
  of skeleton blocks as long as perfomance and resource occupation can be determined. Hence,
  in this paper we will apply our methodology on simple DSP chains: a white noise input signal

  is generated using a Pseudo-Random Number (PRN) generator or thanks at a radiofrequency-grade

  Analog to Digital Converter (ADC) loaded by a 50~$\Omega$ resistor. Once samples have been
  digitized at a rate of 125~MS/s, filtering is applied to qualify the processing block performance --
  practically meeting the radiofrequency frontend requirement of noise and bandwidth reduction
  by filtering and decimating. Finally, bursts of filtered samples are stored for post-processing,

  allowing to assess either filter rejection for a given resource usage, or validating the rejection

  when implementing a solution minimizing resource occupation.

  
  The first step of our approach is to model the DSP chain and since we just optimize
  the filtering, we have not modeling the PRN generator or the ADC. The filtering can be
  done by two ways. The first one we use only one FIR filter with lot of coefficients
  to rejection the noise, we called this approach a monolithic approach. And the second one
  we select different FIR filters with less coefficients the monolithic filter and we cascaded
  it to filtering the signal.
  
  After each filter we leave the possibility of shifting the filtered data to consume
  less resources. Hence in the case of cascaded filter, we define a stage as a filter
  and a shifter (the shift could be omitted if we do not need to divide the filtered data).
  
  \subsection{Model of a FIR filter}

  
  A cascade of filters is composed of $n$ FIR stages. In stage $i$ ($1 \leq i \leq n$)
  the FIR has $C_i$ coefficients and each coefficient is an integer value with $\pi^C_i$
  bits while the filtered data are shifted by $\pi^S_i$ bits. We define also $\pi^-_i$ as

  the size of input data and $\pi^+_i$ as the size of output data. The figure~\ref{fig:fir_stage}
  shows a filtering stage.
  
  \begin{figure}
    \centering
    \begin{tikzpicture}[node distance=2cm]
      
  ode[draw,minimum size=1.3cm] (FIR) { $C_i, \pi_i^C$ } ;
      
  ode[draw,minimum size=1.3cm] (Shift) [right of=FIR, ] { $\pi_i^S$ } ;
      
  ode (Start) [left of=FIR] { } ;
      
  ode (End) [right of=Shift] { } ;
  
      
  ode[draw,fit=(FIR) (Shift)] (Filter) { } ;
  
      \draw[->] (Start) edge node [above] { $\pi_i^-$ } (FIR) ;
      \draw[->] (FIR) -- (Shift) ;
      \draw[->] (Shift) edge node [above] { $\pi_i^+$ } (End) ;
    \end{tikzpicture}
    \caption{A single filter is composed of a FIR (on the left) and a Shifter (on the right)}
    \label{fig:fir_stage}
  \end{figure}

  FIR $i$ has been characterized through numerical simulation as able to reject $F(C_i, \pi_i^C)$ dB.
  This rejection has been computed using GNU Octave software FIR coefficient design functions
  (\texttt{firls} and \texttt{fir1}).

  For each configuration $(C_i, \pi_i^C)$, we first create a FIR with floating point coefficients and a given $C_i$ number of coefficients.
  Then, the floating point coefficients are discretized into integers. In order to ensure that the coefficients are coded on $\pi_i^C$~bits effectively,
  the coefficients are normalized by their absolute maximum before being scaled to integer coefficients.

  At least one coefficient is coded on $\pi_i^C$~bits, and in practice only $b_{C_i/2}$ is coded on $\pi_i^C$~bits while the others are coded on much fewer bits.

  With these coefficients, the \texttt{freqz} function is used to estimate the magnitude of the filter
  transfer function.
  Comparing the performance between FIRs requires however defining a unique criterion. As shown in figure~\ref{fig:fir_mag},

  the FIR magnitude exhibits two parts: we focus here on the transitions width and the rejection rather than on the

  bandpass ripples as emphasized in \cite{lim_1988,lim_1996}.

  
  \begin{figure}

  \begin{center}
  \scalebox{0.8}{

    \centering
    \begin{tikzpicture}[scale=0.3]
      \draw[<->] (0,15) -- (0,0) -- (21,0) ;
      \draw[thick] (0,12) -- (8,12) -- (20,0) ;
  
      \draw (0,14) node [left] { $P$ } ;
      \draw (20,0) node [below] { $f$ } ;
  
      \draw[>=latex,<->] (0,14) -- (8,14) ;
      \draw (4,14) node [above] { passband } node [below] { $40\%$ } ;
  
      \draw[>=latex,<->] (8,14) -- (12,14) ;
      \draw (10,14) node [above] { transition } node [below] { $20\%$ } ;
  
      \draw[>=latex,<->] (12,14) -- (20,14) ;
      \draw (16,14) node [above] { stopband } node [below] { $40\%$ } ;
  
      \draw[>=latex,<->] (16,12) -- (16,8) ;
      \draw (16,10) node [right] { rejection } ;
  
      \draw[dashed] (8,-1) -- (8,14) ;
      \draw[dashed] (12,-1) -- (12,14) ;
  
      \draw[dashed] (8,12) -- (16,12) ;
      \draw[dashed] (12,8) -- (16,8) ;
  
    \end{tikzpicture}

  }
  \end{center}

  \caption{Shape of the filter transmitted power $P$ as a function of frequency $f$:
  the passband is considered to occupy the initial 40\% of the Nyquist frequency range,
  the stopband the last 40\%, allowing 20\% transition width.}
  \label{fig:fir_mag}

  \end{figure}

  In the transition band, the behavior of the filter is left free, we only care about the passband and the stopband characteristics.

  Our initial criterion considered the mean value of the stopband rejection, as shown in figure~\ref{fig:mean_criterion}. This criterion

  yields unacceptable results since notches overestimate the rejection capability of the filter. Furthermore, the losses within
  the passband are not considered and might be excessive for excessively wide transitions widths introduced for filters with few coefficients.
  Such biases are compensated for by the second considered criterion which is based on computing the maximum rejection within the stopband minus the mean of the absolute value of passband rejection. With this criterion, the results are significantly improved as shown in figure~\ref{fig:custom_criterion} and meet the expected rejection capability of low pass filters.

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/colored_mean_criterion}

  \caption{Mean stopband rejection criterion comparison between monolithic filter and cascaded filters}

  \label{fig:mean_criterion}

  \end{figure}

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/colored_custom_criterion}

  \caption{Custom criterion (maximum rejection in the stopband minus the mean of the absolute value of the passband rejection)

  comparison between monolithic filter and cascaded filters}

  \label{fig:custom_criterion}
  \end{figure}

  Thanks to the latter criterion which will be used in the remainder of this paper, we are able to automatically generate multiple FIR taps
  and estimate their rejection. Figure~\ref{fig:rejection_pyramid} exhibits the
  rejection as a function of the number of coefficients and the number of bits representing these coefficients.
  The curve shaped as a pyramid exhibits optimum configurations sets at the vertex where both edges meet.

  Indeed for a given number of coefficients, increasing the number of bits over the edge will not improve the rejection.

  Conversely when setting the a given number of bits, increasing the number of coefficients will not improve
  the rejection. Hence the best coefficient set are on the vertex of the pyramid.

  
  \begin{figure}
  \centering
  \includegraphics[width=\linewidth]{images/rejection_pyramid}
  \caption{Rejection as a function of number of coefficients and number of bits}
  \label{fig:rejection_pyramid}
  \end{figure}

  Although we have an efficient criterion to estimate the rejection of one set of coefficients (taps),

  we have a problem when we cascade filters and estimate the criterion as a sum two or more individual criteria.

  If the FIR filter coefficients are the same between the stages, we have:

  $$F_{total} = F_1 + F_2$$

  But selecting two different sets of coefficient will yield a more complex situation in which
  the previous relation is no longer valid as illustrated on figure~\ref{fig:sum_rejection}. The red and blue curves
  are two different filters with maximums and notches not located at the same frequency offsets.

  Hence when summing the transfer functions, the resulting rejection shown as the dashed yellow line is improved

  with respect to a basic sum of the rejection criteria shown as a the dotted yellow line.

  Thus, estimating the rejection of filter cascades is more complex than takin the sum of all the rejection

  criteria of each filter. However since the this sum underestimates the rejection capability of the cascade,
  this upper bound is considered as a pessimistic and acceptable criterion for deciding on the suitability
  of the filter cascade to meet design criteria.

  
  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/cascaded_criterion}

  \caption{Rejection of two cascaded filters}
  \label{fig:sum_rejection}
  \end{figure}

  Based on this analysis, we address the estimate of resource consumption (called
  silicon area -- in the case of FPGAs meaning processing cells) as a function of
  filter characteristics. As a reminder, we do not aim at matching actual hardware
  configuration but consider an arbitrary silicon area occupied by each processing function,
  and will assess after synthesis the adequation of this arbitrary unit with actual
  hardware resources provided by FPGA manufacturers. The sum of individual processing
  unit areas is constrained by a total silicon area representative of FPGA global resources.
  Formally, variable $a_i$ is the area taken by filter~$i$

  (in arbitrary unit). Variable $r_i$ is the rejection of filter~$i$ (in dB).
  Constant $\mathcal{A}$ is the total available area. We model our problem as follows:

  \begin{align}
  \text{Maximize } & \sum_{i=1}^n r_i  
  otag \\
  \sum_{i=1}^n a_i & \leq \mathcal{A} & \label{eq:area} \\
  a_i & = C_i \times (\pi_i^C + \pi_i^-), & \forall i \in [1, n] \label{eq:areadef} \\
  r_i & = F(C_i, \pi_i^C), & \forall i \in [1, n] \label{eq:rejectiondef} \\
  \pi_i^+ & = \pi_i^- + \pi_i^C - \pi_i^S, & \forall i \in [1, n] \label{eq:bits} \\
  \pi_{i - 1}^+ & = \pi_i^-, & \forall i \in [2, n] \label{eq:inout} \\
  \pi_i^+ & \geq 1 + \sum_{k=1}^{i} \left(1 + \frac{r_j}{6}\right), & \forall i \in [1, n] \label{eq:maxshift} \\
  \pi_1^- &= \Pi^I \label{eq:init}
  \end{align}

  Equation~\ref{eq:area} states that the total area taken by the filters must be
  less than the available area. Equation~\ref{eq:areadef} gives the definition of

  the area used by a filter, considered as the area of the FIR since the Shifter is
  assumed not to require significant resources. We consider that the FIR needs $C_i$ registers of size

  $\pi_i^C + \pi_i^-$~bits to store the results of the multiplications of the

  input data with the coefficients. Equation~\ref{eq:rejectiondef} gives the
  definition of the rejection of the filter thanks to the tabulated function~$F$ that we defined
  previously. The Shifter does not introduce negative rejection as we will explain later,

  so the rejection only comes from the FIR. Equation~\ref{eq:bits} states the
  relation between $\pi_i^+$ and $\pi_i^-$. The multiplications in the FIR add
  $\pi_i^C$ bits as most coefficients are close to zero, and the Shifter removes
  $\pi_i^S$ bits. Equation~\ref{eq:inout} states that the output number of bits of
  a filter is the same as the input number of bits of the next filter.
  Equation~\ref{eq:maxshift} ensures that the Shifter does not introduce negative
  rejection. Indeed, the results of the FIR can be right shifted without compromising
  the quality of the rejection until a threshold. Each bit of the output data

  increases the maximum rejection level by 6~dB. We add one to take the sign bit

  into account. If equation~\ref{eq:maxshift} was not present, the Shifter could
  shift too much and introduce some noise in the output data. Each supplementary

  shift bit would cause an additional 6~dB rejection rise. A totally equivalent equation is:
  $\pi_i^S \leq \pi_i^- + \pi_i^C - 1 - \sum_{k=1}^{i} \left(1 + \frac{r_j}{6}\right)$.
  Finally, equation~\ref{eq:init} gives the number of bits of the global input.

  
  This model is non-linear and even non-quadratic, as $F$ does not have a known
  linear or quadratic expression. We introduce $p$ FIR configurations

  $(C_{ij}, \pi_{ij}^C), 1 \leq j \leq p$ that are constants. We define binary
  variable $\delta_{ij}$ that has value 1 if stage~$i$ is in configuration~$j$
  and 0 otherwise. The new equations are as follows:

  
  \begin{align}
  a_i & = \sum_{j=1}^p \delta_{ij} \times C_{ij} \times (\pi_{ij}^C + \pi_i^-), & \forall i \in [1, n] \label{eq:areadef2} \\
  r_i & = \sum_{j=1}^p \delta_{ij} \times F(C_{ij}, \pi_{ij}^C), & \forall i \in [1, n] \label{eq:rejectiondef2} \\
  \pi_i^+ & = \pi_i^- + \left(\sum_{j=1}^p \delta_{ij} \pi_{ij}^C\right) - \pi_i^S, & \forall i \in [1, n] \label{eq:bits2} \\
  \sum_{j=1}^p \delta_{ij} & \leq 1, & \forall i \in [1, n] \label{eq:config}
  \end{align}
  
  Equations \ref{eq:areadef2}, \ref{eq:rejectiondef2} and \ref{eq:bits2} replace
  respectively equations \ref{eq:areadef}, \ref{eq:rejectiondef} and \ref{eq:bits}.
  Equation~\ref{eq:config} states that for each stage, a single configuration is chosen at most.

  This modified model is quadratic, and it can be linearised if necessary. The Gurobi
  (\url{www.gurobi.com}) optimization software is used to solve this quadratic
  model, and since Gurobi is able to linearize, the model is left as is. This model
  has $O(np)$ variables and $O(n)$ constraints.

  Two problems will be addressed using the workflow described in the next section: on the one
  hand maximizing the rejection capability of a set of cascaded filters occupying a fixed arbitrary
  silcon area (section~\ref{sec:fixed_area}) and on the second hand the dual problem of minimizing the silicon area
  for a fixed rejection criterion (section~\ref{sec:fixed_rej}). In the latter case, the
  objective function is replaced with:

  \begin{align}
  \text{Minimize } & \sum_{i=1}^n a_i  
  otag
  \end{align}

  We adapt our constraints of quadratic program to replace equation \ref{eq:area}
  with equation \ref{eq:rejection_min} where $\mathcal{R}$ is the minimal

  rejection required.
  
  \begin{align}
  \sum_{i=1}^n r_i & \geq \mathcal{R} & \label{eq:rejection_min}
  \end{align}
  
  \section{Design workflow}
  \label{sec:workflow}

  In this section, we describe the workflow to compute all the results presented in sections~\ref{sec:fixed_area}

  and \ref{sec:fixed_rej}. Figure~\ref{fig:workflow} shows the global workflow and the different steps involved

  in the computation of the results.

  
  \begin{figure}
    \centering
    \begin{tikzpicture}[node distance=0.75cm and 2cm]
      
  ode[draw,minimum size=1cm] (Solver) { Filter Solver } ;
      
  ode (Start) [left= 3cm of Solver] { } ;
      
  ode[draw,minimum size=1cm] (TCL) [right= of Solver] { TCL Script } ;
      
  ode (Input) [above= of TCL] { } ;
      
  ode[draw,minimum size=1cm] (Deploy) [below= of Solver] { Deploy Script } ;
      
  ode[draw,minimum size=1cm] (Bitstream) [below= of TCL] { Bitstream } ;
      
  ode[draw,minimum size=1cm,rounded corners] (Board) [below right= of Deploy] { Board } ;
      
  ode[draw,minimum size=1cm] (Postproc) [below= of Deploy] { Post-Processing } ;
      
  ode (Results) [left= of Postproc] { } ;
  
      \draw[->] (Start) edge node [above] { $\mathcal{A}, n, \Pi^I$ } node [below] { $(C_{ij}, \pi_{ij}^C), F$ } (Solver) ;
      \draw[->] (Input) edge node [left] { ADC or PRN } (TCL) ;
      \draw[->] (Solver) edge node [below] { (1a) } (TCL) ;
      \draw[->] (Solver) edge node [right] { (1b) } (Deploy) ;
      \draw[->] (TCL) edge node [left] { (2) } (Bitstream) ;
      \draw[->,dashed] (Bitstream) -- (Deploy) ;
      \draw[->] (Deploy) to[out=-30,in=120] node [above] { (3) } (Board) ;
      \draw[->] (Board) to[out=150,in=-60] node [below] { (4) } (Deploy) ;
      \draw[->] (Deploy) edge node [left] { (5) } (Postproc) ;
      \draw[->] (Postproc) -- (Results) ;
    \end{tikzpicture}
    \caption{Design workflow from the input parameters to the results}
    \label{fig:workflow}
  \end{figure}
  
  The filter solver is a C++ program that takes as input the maximum area
  $\mathcal{A}$, the number of stages $n$, the size of the input signal $\Pi^I$,
  the FIR configurations $(C_{ij}, \pi_{ij}^C)$ and the function $F$. It creates

  the quadratic programs and uses the Gurobi solver to estimate the optimal results.

  Then it produces two scripts: a TCL script ((1a) on figure~\ref{fig:workflow})
  and a deploy script ((1b) on figure~\ref{fig:workflow}).
  
  The TCL script describes the whole digital processing chain from the beginning

  (the raw signal data) to the end (the filtered data) in a language compatible

  with proprietary synthesis software, namely Vivado for Xilinx and Quartus for

  Intel/Altera. The raw input data generated from a 20-bit Pseudo Random Number (PRN)

  generator inside the FPGA and $\Pi^I$ is fixed at 16~bits.
  Then the script builds each stage of the chain with a generic FIR task that
  comes from a skeleton library. The generic FIR is highly configurable
  with the number of coefficients and the size of the coefficients. The coefficients
  themselves are not stored in the script.

  As the signal is processed in real-time, the output signal is stored as
  consecutive bursts of data for post-processing, mainly assessing the consistency of the
  implemented FIR cascade transfer function with the design criteria and the expected
  transfer function.

  
  The TCL script is used by Vivado to produce the FPGA bitstream ((2) on figure~\ref{fig:workflow}).
  We use the 2018.2 version of Xilinx Vivado and we execute the synthesized
  bitstream on a Redpitaya board fitted with a Xilinx Zynq-7010 series

  FPGA (xc7z010clg400-1) and two LTC2145 14-bit 125~MS/s ADC, loaded with 50~$\Omega$ resistors to
  provide a broadband noise source.
  The board runs the Linux kernel and surrounding environment produced from the
  Buildroot framework available at \url{https://github.com/trabucayre/redpitaya/}: configuring

  the Zynq FPGA, feeding the FIR with the set of coefficients, executing the simulation and

  fetching the results is automated.

  
  The deploy script uploads the bitstream to the board ((3) on
  figure~\ref{fig:workflow}), flashes the FPGA, loads the different drivers,
  configures the coefficients of the FIR filters. It then waits for the results
  and retrieves the data to the main computer ((4) on figure~\ref{fig:workflow}).
  
  Finally, an Octave post-processing script computes the final results thanks to
  the output data ((5) on figure~\ref{fig:workflow}).
  The results are normalized so that the Power Spectrum Density (PSD) starts at zero
  and the different configurations can be compared.

  \section{Maximizing the rejection at fixed silicon area}

  \label{sec:fixed_area}
  This section presents the output of the filter solver {\em i.e.} the computed
  configurations for each stage, the computed rejection and the computed silicon area.

  Such results allow for understanding the choices made by the solver to compute its solutions.

  
  The experimental setup is composed of three cases. The raw input is generated
  by a Pseudo Random Number (PRN) generator, which fixes the input data size $\Pi^I$.
  Then the total silicon area $\mathcal{A}$ has been fixed to either 500, 1000 or 1500
  arbitrary units. Hence, the three cases have been named: MAX/500, MAX/1000, MAX/1500.
  The number of configurations $p$ is 1827, with $C_i$ ranging from 3 to 60 and $\pi^C$
  ranging from 2 to 22. In each case, the quadratic program has been able to give a
  result up to five stages ($n = 5$) in the cascaded filter.
  
  Table~\ref{tbl:gurobi_max_500} shows the results obtained by the filter solver for MAX/500.
  Table~\ref{tbl:gurobi_max_1000} shows the results obtained by the filter solver for MAX/1000.
  Table~\ref{tbl:gurobi_max_1500} shows the results obtained by the filter solver for MAX/1500.
  
  \renewcommand{\arraystretch}{1.4}
  
  \begin{table}
    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MAX/500}
    \label{tbl:gurobi_max_500}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area  \\
          \hline
              1 & (21, 7, 0)  & -           & -           & -           & -           & 32~dB           & 483   \\
              2 & (3, 3, 15)  & (31, 9, 0)  & -           & -           & -           & 58~dB           & 460   \\
              3 & (3, 3, 15)  & (27, 9, 0)  & (5, 3, 0)   & -           & -           & 66~dB           & 488   \\
              4 & (3, 3, 15)  & (19, 7, 0)  & (11, 5, 0)  & (3, 3, 0)   & -           & 74~dB           & 499   \\
              5 & (3, 3, 15)  & (23, 8, 0)  & (3, 3, 1)   & (3, 3, 0)   & (3, 3, 0)   & 78~dB           & 489   \\
          \hline
        \end{tabular}
      }
  \end{table}
  
  \begin{table}
    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MAX/1000}
    \label{tbl:gurobi_max_1000}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area \\
          \hline
              1 & (37, 11, 0) & -           & -           & -           & -           & 56~dB           & 999  \\
              2 & (3, 3, 15)  & (51, 14, 0) & -           & -           & -           & 87~dB           & 975  \\
              3 & (3, 3, 15)  & (35, 11, 0) & (19, 7, 0)  & -           & -           & 99~dB           & 1000 \\
              4 & (3, 4, 16)  & (27, 8, 0)  & (19, 7, 1)  & (11, 5, 0)  & -           & 103~dB          & 998  \\
              5 & (3, 3, 15)  & (31, 9, 0)  & (19, 7, 0)  & (3, 3, 1)   & (3, 3, 0)   & 111~dB          & 984  \\
          \hline
        \end{tabular}
      }
  \end{table}
  
  \begin{table}
    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MAX/1500}
    \label{tbl:gurobi_max_1500}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area  \\
          \hline
              1 & (47, 15, 0) & -           & -           & -           & -           & 71~dB           & 1457  \\
              2 & (19, 6, 15) & (51, 14, 0) & -           & -           & -           & 103~dB          & 1489  \\
              3 & (3, 3, 15)  & (35, 11, 0) & (35, 11, 0) & -           & -           & 122~dB          & 1492  \\
              4 & (3, 3, 15)  & (27, 8, 0)  & (19, 7, 0)  & (27, 9, 0)  & -           & 129~dB          & 1498  \\
              5 & (3, 3, 15)  & (23, 9, 2)  & (27, 9, 0)  & (19, 7, 0)  & (3, 3, 0)   & 136~dB          & 1499  \\
          \hline
        \end{tabular}
      }
  \end{table}
  
  \renewcommand{\arraystretch}{1}
  
  From these tables, we can first state that the more stages are used to define
  the cascaded FIR filters, the better the rejection. It was an expected result as it has
  been previously observed that many small filters are better than

  a single large filter \cite{lim_1988, lim_1996, young_1992}, despite such conclusions

  being hardly used in practice due to the lack of tools for identifying individual filter
  coefficients in the cascaded approach.
  
  Second, the larger the silicon area, the better the rejection. This was also an

  expected result as more area means a filter of better quality with more coefficients
  or more bits per coefficient.

  
  Then, we also observe that the first stage can have a larger shift than the other
  stages. This is explained by the fact that the solver tries to use just enough
  bits for the computed rejection after each stage. In the first stage, a
  balance between a strong rejection with a low number of bits is targeted. Equation~\ref{eq:maxshift}
  gives the relation between both values.
  
  Finally, we note that the solver consumes all the given silicon area.

  The following graphs present the rejection for real data on the FPGA. In all the following

  figures, the solid line represents the actual rejection of the filtered

  data on the FPGA as measured experimentally and the dashed line are the noise levels

  given by the quadratic solver. The configurations are those computed in the previous section.
  
  Figure~\ref{fig:max_500_result} shows the rejection of the different configurations in the case of MAX/500.
  Figure~\ref{fig:max_1000_result} shows the rejection of the different configurations in the case of MAX/1000.
  Figure~\ref{fig:max_1500_result} shows the rejection of the different configurations in the case of MAX/1500.

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/max_500}
  \caption{Signal spectrum for MAX/500}
  \label{fig:max_500_result}

  \end{figure}

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/max_1000}
  \caption{Signal spectrum for MAX/1000}
  \label{fig:max_1000_result}

  \end{figure}

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/max_1500}
  \caption{Signal spectrum for MAX/1500}
  \label{fig:max_1500_result}

  \end{figure}

  In all cases, we observe that the actual rejection is close to the rejection computed by the solver.
  
  We compare the actual silicon resources given by Vivado to the
  resources in arbitrary units.
  The goal is to check that our arbitrary units of silicon area models well enough
  the real resources on the FPGA. Especially we want to verify that, for a given
  number of arbitrary units, the actual silicon resources do not depend on the
  number of stages $n$. Most significantly, our approach aims
  at remaining far enough from the practical logic gate implementation used by
  various vendors to remain platform independent and be portable from one
  architecture to another.
  
  Table~\ref{tbl:resources_usage} shows the resources usage in the case of MAX/500, MAX/1000 and
  MAX/1500 \emph{i.e.} when the maximum allowed silicon area is fixed to 500, 1000
  and 1500 arbitrary units. We have taken care to extract solely the resources used by

  the FIR filters and remove additional processing blocks including FIFO and Programmable
  Logic (PL -- FPGA) to Processing System (PS -- general purpose processor) communication.

  \begin{table}[h!tb]

    \caption{Resource occupation. The last column refers to available resources on a Zynq-7010 as found on the Redpitaya.}
    \label{tbl:resources_usage}
    \centering
        \begin{tabular}{|c|c|ccc|c|}
          \hline
          $n$ &          & MAX/500  & MAX/1000 & MAX/1500 & \emph{Zynq 7010}         \\ \hline\hline
              & LUT      & 249      & 453      & 627      & \emph{17600}             \\
          1   & BRAM     & 1        & 1        & 1        & \emph{120}               \\
              & DSP      & 21       & 37       & 47       & \emph{80}                \\ \hline
              & LUT      & 2374     & 5494     & 691      & \emph{17600}             \\
          2   & BRAM     & 2        & 2        & 2        & \emph{120}               \\
              & DSP      & 0        & 0        & 70       & \emph{80}                \\ \hline
              & LUT      & 2443     & 3304     & 3521     & \emph{17600}             \\
          3   & BRAM     & 3        & 3        & 3        & \emph{120}               \\
              & DSP      & 0        & 19       & 35       & \emph{80}                \\ \hline
              & LUT      & 2634     & 3753     & 2557     & \emph{17600}             \\
          4   & BRAM     & 4        & 4        & 4        & \emph{120}               \\
              & DPS      & 0        & 19       & 46       & \emph{80}                \\ \hline
              & LUT      & 2423     & 3047     & 2847     & \emph{17600}             \\
          5   & BRAM     & 5        & 5        & 5        & \emph{120}               \\
              & DPS      & 0        & 22       & 46       & \emph{80}                \\ \hline
        \end{tabular}

  \end{table}

  In some cases, Vivado replaces the DSPs by Look Up Tables (LUTs). We assume that,

  when the filter coefficients are small enough, or when the input size is small
  enough, Vivado optimizes resource consumption by selecting multiplexers to

  implement the multiplications instead of a DSP. In this case, it is quite difficult
  to compare the whole silicon budget.

  However, a rough estimation can be made with a simple equivalence: looking at

  the first column (MAX/500), where the number of LUTs is quite stable for $n \geq 2$,
  we can deduce that a DSP is roughly equivalent to 100~LUTs in terms of silicon

  area use. With this equivalence, our 500 arbitraty units correspond to 2500 LUTs,
  1000 arbitrary units correspond to 5000 LUTs and 1500 arbitrary units correspond

  to 7300 LUTs. The conclusion is that the orders of magnitude of our arbitrary

  unit map well to actual hardware resources. The relatively small differences can probably be explained

  by the optimizations done by Vivado based on the detailed map of available processing resources.

  We now present the computation time needed to solve the quadratic problem.
  For each case, the filter solver software is executed on a Intel(R) Xeon(R) CPU E5606
  clocked at 2.13~GHz. The CPU has 8 cores that are used by Gurobi to solve
  the quadratic problem. Table~\ref{tbl:area_time} shows the time needed to solve the quadratic

  problem when the maximal area is fixed to 500, 1000 and 1500 arbitrary units.

  \begin{table}[h!tb]
  \caption{Time needed to solve the quadratic program with Gurobi}

  \label{tbl:area_time}

  \centering

  \begin{tabular}{|c|c|c|c|}\hline
  $n$ & Time (MAX/500)          & Time (MAX/1000)             & Time (MAX/1500)              \\\hline\hline
  1   & 0.1~s                   & 0.1~s                       & 0.3~s                        \\
  2   & 1.1~s                   & 2.2~s                       & 12~s                         \\
  3   & 17~s                    & 137~s  ($\approx$ 2~min)    & 275~s ($\approx$ 4~min)      \\
  4   & 52~s                    & 5448~s ($\approx$ 90~min)   & 5505~s ($\approx$ 17~h)      \\
  5   & 286~s ($\approx$ 4~min) & 4119~s ($\approx$ 68~min)   & 235479~s ($\approx$ 3~days)  \\\hline

  \end{tabular}

  \end{table}

  As expected, the computation time seems to rise exponentially with the number of stages. % TODO: exponentiel ?
  When the area is limited, the design exploration space is more limited and the solver is able to

  find an optimal solution faster.

  
  \subsection{Minimizing resource occupation at fixed rejection}\label{sec:fixed_rej}

  This section presents the results of the complementary quadratic program aimed at
  minimizing the area occupation for a targeted rejection level.

  
  The experimental setup is also composed of three cases. The raw input is the same

  as in the previous section, from a PRN generator, which fixes the input data size $\Pi^I$.

  Then the targeted rejection $\mathcal{R}$ has been fixed to either 40, 60 or 80~dB.
  Hence, the three cases have been named: MIN/40, MIN/60, MIN/80.
  The number of configurations $p$ is the same as previous section.
  
  Table~\ref{tbl:gurobi_min_40} shows the results obtained by the filter solver for MIN/40.
  Table~\ref{tbl:gurobi_min_60} shows the results obtained by the filter solver for MIN/60.
  Table~\ref{tbl:gurobi_min_80} shows the results obtained by the filter solver for MIN/80.
  
  \renewcommand{\arraystretch}{1.4}

  \begin{table}[h!tb]

    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MIN/40}
    \label{tbl:gurobi_min_40}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area  \\
          \hline
              1 & (27, 8, 0)  & -           & -           & -           & -           & 41~dB           & 648   \\
              2 & (3, 2, 14)  & (19, 7, 0)  & -           & -           & -           & 40~dB           & 263   \\
              3 & (3, 3, 15)  & (11, 5, 0)  & (3, 3, 0)   & -           & -           & 41~dB           & 192   \\
              4 & (3, 3, 15)  & (3, 3, 0)   & (3, 3, 0)   & (3, 3, 0)   & -           & 42~dB           & 147   \\
          \hline
        \end{tabular}
      }

  \end{table}

  \begin{table}[h!tb]

    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MIN/60}
    \label{tbl:gurobi_min_60}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area \\
          \hline
              1 & (39, 13, 0) & -           & -           & -           & -           & 60~dB           & 1131 \\
              2 & (3, 3, 15)  & (35, 10, 0) & -           & -           & -           & 60~dB           & 547  \\
              3 & (3, 3, 15)  & (27, 8, 0)  & (3, 3, 0)   & -           & -           & 62~dB           & 426  \\
              4 & (3, 2, 14)  & (11, 5, 1)  & (11, 5, 0)  & (3, 3, 0)   & -           & 60~dB           & 344  \\
              5 & (3, 2, 14)  & (3, 3, 1)   & (3, 3, 0)   & (3, 3, 0)   & (3, 3, 0)   & 60~dB           & 279  \\
          \hline
        \end{tabular}
      }

  \end{table}

  \begin{table}[h!tb]

    \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MIN/80}
    \label{tbl:gurobi_min_80}
    \centering
      {\scalefont{0.77}
        \begin{tabular}{|c|ccccc|c|c|}
          \hline
           $n$  & $i = 1$     & $i = 2$     & $i = 3$     & $i = 4$     & $i = 5$     & Rejection       & Area  \\
          \hline
              1 & (55, 16, 0) & -           & -           & -           & -           & 81~dB           & 1760  \\
              2 & (3, 3, 15)  & (47, 14, 0) & -           & -           & -           & 80~dB           & 903   \\
              3 & (3, 3, 15)  & (23, 9, 0)  & (19, 7, 0)  & -           & -           & 80~dB           & 698   \\
              4 & (3, 3, 15)  & (27, 9, 0)  & (7, 7, 4)   & (3, 3, 0)   & -           & 80~dB           & 605   \\
              5 & (3, 2, 14)  & (27, 8, 0)  & (3, 3, 1)   & (3, 3, 0)   & (3, 3, 0)   & 81~dB           & 534   \\
          \hline
        \end{tabular}
      }

  \end{table}

  \renewcommand{\arraystretch}{1}

  % JMF : je croyais que dans un cas le monolithique n'y arrivait juste pas : tu as retire' ce cas ?
  From these tables, we can first state that all configurations reach the targeted rejection
  level or even better thanks to our underestimate of the cascade rejection as the sum of the
  individual filter rejection
  % we have stages lesser is the area occupied in arbitrary unit. JMF : je ne comprends pas cette phrase
  Futhermore, the area of the monolithic filter is twice as big as the two cascaded filters
  (1131 and 1760  arbitrary units v.s 547 and 903 arbitrary units for 60 and 80~dB rejection
  respectively). More generally, the more filters are cascaded, the lower the occupied area.
  
  Like in previous section, the solver chooses always a little filter as first
  filter stage and the second one is often the biggest filter. This choice can be explained
  as in the previous section, with the solver using just enough bits not to degrade the input
  signal and in the second filter selecting a better filter to improve rejection without
  having too many bits in the output data.
  
  For the specific case of MIN/40 for $n = 5$ the solver has determined that the optimal
  number of filters is 4 so it did not chose any configuration for the last filter. Hence this

  solution is equivalent to the result for $n = 4$.

  The following graphs present the rejection for real data on the FPGA. In all the following

  figures, the solid line represents the actual rejection of the filtered

  data on the FPGA as measured experimentally and the dashed line is the noise level

  given by the quadratic solver.
  
  Figure~\ref{fig:min_40} shows the rejection of the different configurations in the case of MIN/40.
  Figure~\ref{fig:min_60} shows the rejection of the different configurations in the case of MIN/60.
  Figure~\ref{fig:min_80} shows the rejection of the different configurations in the case of MIN/80.

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/min_40}
  \caption{Signal spectrum for MIN/40}
  \label{fig:min_40}

  \end{figure}

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/min_60}
  \caption{Signal spectrum for MIN/60}
  \label{fig:min_60}

  \end{figure}

  \begin{figure}
  \centering

  \includegraphics[width=\linewidth]{images/min_80}
  \caption{Signal spectrum for MIN/80}
  \label{fig:min_80}

  \end{figure}

  We observe that all rejections given by the quadratic solver are close to the experimentally
  measured rejection. All curves prove that the constraint to reach the target rejection is
  respected with both monolithic or cascaded filters.

  Table~\ref{tbl:resources_usage} shows the resource usage in the case of MIN/40, MIN/60 and

  MIN/80 \emph{i.e.} when the target rejection is fixed to 40, 60 and 80~dB. We
  have taken care to extract solely the resources used by
  the FIR filters and remove additional processing blocks including FIFO and PL to
  PS communication.
  
  \begin{table}
    \caption{Resource occupation. The last column refers to available resources on a Zynq-7010 as found on the Redpitaya.}
    \label{tbl:resources_usage_comp}
    \centering
        \begin{tabular}{|c|c|ccc|c|}
          \hline
          $n$ &          & MIN/40   & MIN/60   & MIN/80   & \emph{Zynq 7010}         \\ \hline\hline
              & LUT      & 343      & 334      & 772      & \emph{17600}             \\
          1   & BRAM     & 1        & 1        & 1        & \emph{120}               \\
              & DSP      & 27       & 39       & 55       & \emph{80}                \\ \hline
              & LUT      & 1252     & 2862     & 5099     & \emph{17600}             \\
          2   & BRAM     & 2        & 2        & 2        & \emph{120}               \\
              & DSP      & 0        & 0        & 0        & \emph{80}                \\ \hline
              & LUT      & 891      & 2148     & 2023     & \emph{17600}             \\
          3   & BRAM     & 3        & 3        & 3        & \emph{120}               \\
              & DSP      & 0        & 0        & 19       & \emph{80}                \\ \hline
              & LUT      & 662      & 1729     & 2451     & \emph{17600}             \\
          4   & BRAM     & 4        & 4        & 4        & \emph{120}               \\
              & DPS      & 0        & 0        & 7        & \emph{80}                \\ \hline
              & LUT      & -        & 1259     & 2602     & \emph{17600}             \\
          5   & BRAM     & -        & 5        & 5        & \emph{120}               \\
              & DPS      & -        & 0        & 0        & \emph{80}                \\ \hline
        \end{tabular}
  \end{table}

  If we keep the previous estimation of cost of one DSP in terms of LUT (1 DSP $\approx$ 100 LUT)
  the real resource consumption decreases as a function of the number of stages in the cascaded
  filter according

  to the solution given by the quadratic solver. Indeed, we have always a decreasing
  consumption even if the difference between the monolithic and the two cascaded

  filters is less than expected.

  Finally, table~\ref{tbl:area_time_comp} shows the computation time to solve

  the quadratic program.

  \begin{table}[h!tb]

  \caption{Time to solve the quadratic program with Gurobi}
  \label{tbl:area_time_comp}
  \centering
  \begin{tabular}{|c|c|c|c|}\hline
  $n$ & Time (MIN/40)           & Time (MIN/60)               & Time (MIN/80)  \\\hline\hline
  1   & 0.07~s                  & 0.02~s                      & 0.01~s         \\
  2   & 7.8~s                   & 16~s                        & 14~s           \\
  3   & 4.7~s                   & 14~s                        & 28~s           \\
  4   & 39~s                    & 20~s                        & 193~s          \\
  5   & 126~s                   & 12~s                        & 170~s          \\\hline
  \end{tabular}
  \end{table}

  The time needed to solve this configuration is significantly shorter than the time
  needed in the previous section. Indeed the worst time in this case is only 3~minutes,
  compared to 3~days in the previous section: this problem is more easily solved than the
  previous one.

  \section{Conclusion}

  We have proposed a new approach to schedule a set of signal processing blocks whose performances
  and resource consumption has been tabulated, and applied this methodology to the practical
  case of implementing cascaded FIR filters inside a FPGA.
  This method aims to be hardware independent and focuses an a high-level of abstraction.
  We have modeled the FIR filter operation and the impact of data shift. Thanks to this model,
  we have created a quadratic program to select the optimal FIR taps to reach a targeted
  rejection. Individual filter taps have been identified using commonly available tools and the
  emphasis is on FIR assembly rather than individual FIR coefficient identification.

  
  Our experimental results are very promising in providing a rational approach to selecting
  the coefficients of each FIR filter in the context of a performance target for a chain of

  such filters. The FPGA design that is produced automatically by the proposed
  workflow is able to filter an input signal as expected, validating experimentally our model and our approach.
  The quadratic program can be adapted it to an other problem based on assembling skeleton blocks.

  
  A perspective is to model and add the decimators to the processing chain to have a classical

  FIR filter and decimator. The impact of the decimator is not trivial, especially in terms of silicon
  area usage for subsequent stages since some hardware optimization can be applied in

  this case.
  
  The software used to demonstrate the concepts developed in this paper is based on the
  CPU-FPGA co-design framework available at \url{https://github.com/oscimp/oscimpDigital}.

  \section*{Acknowledgement}
  
  This work is supported by the ANR Programme d'Investissement d'Avenir in
  progress at the Time and Frequency Departments of the FEMTO-ST Institute
  (Oscillator IMP, First-TF and Refimeve+), and by R\'egion de Franche-Comt\'e.

  The authors would like to thank E. Rubiola, F. Vernotte, and G. Cabodevila

  for support and fruitful discussions.
  
  \bibliographystyle{IEEEtran}
  \balance
  \bibliography{references,biblio}
  \end{document}