% 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,journal]{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}
\usepackage{caption}
\usepackage{subcaption}

% 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}
\newtheorem{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: processing
resources
are hence saved 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 performance (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 performance. 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 performance 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 by sampling a wideband (125~MS/s)
14-bit 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. Since we aim at only optimizing
the filtering part of the signal processing chain, we have not included the PRN generator or the
ADC in the model: the input data size and rate are considered fixed and defined by the hardware.
The filtering can be done in two ways, either by considering a single monolithic FIR filter
requiring many coefficients to reach the targeted noise rejection ratio, or by
cascading multiple FIR filters, each with fewer coefficients than found in the monolithic filter.

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]
    \node[draw,minimum size=1.3cm] (FIR) { $C_i, \pi_i^C$ } ;
    \node[draw,minimum size=1.3cm] (Shift) [right of=FIR, ] { $\pi_i^S$ } ;
    \node (Start) [left of=FIR] { } ;
    \node (End) [right of=Shift] { } ;

    \node[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}. Throughout this demonstration,
we arbitrarily set a bandpass of 40\% of the Nyquist frequency and a bandstop from 60\%
of the Nyquist frequency to the end of the band, as would be typically selected to prevent
aliasing before decimating the dataflow by 2. The method is however generalized to any filter
shape as long as it is defined from the initial modeling steps: Fig. \ref{fig:rejection_pyramid}
as described below is indeed unique for each filter shape.

\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 define the passband and the stopband characteristics.
% r2.7
Initial considered criteria include the mean value of the stopband rejection which 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.
Our final criterion to compute the filter rejection considers
% r2.8 et r2.2 r2.3
the minimal rejection within the stopband, to which the sum of the absolute values
within the passband is subtracted to avoid filters with excessive ripples, normalized to the
bin width to remain consistent with the passband criterion (dBc/Hz units in all cases). With this
criterion, we meet the expected rejection capability of low pass filters as shown in figure~\ref{fig:custom_criterion}.

% \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 sum of the
absolute values of the passband rejection normalized to the bandwidth)
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{Filter rejection as a function of number of coefficients and number of bits
: this lookup table will be used to identify which filter parameters -- number of bits
representing coefficients and number of coefficients -- best match the targeted transfer function.}
\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.
% r2.9
Thus, estimating the rejection of filter cascades is more complex than taking the sum of all the rejection
criteria of each filter. However since the individual filter rejection sum underestimates the rejection capability of the cascade,
% r2.10
this upper bound is considered as a conservative 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{Transfer function of individual filters and after cascading the two filters,
demonstrating that the selected criterion of maximum rejection in the bandstop (horizontal
lines) is met. Notice that the cascaded filter has better rejection than summing the bandstop
maximum of each individual filter.
}
\label{fig:sum_rejection}
\end{figure}

Finally in our case, we consider that the input signal are fully known. The
resolution of the input data stream are fixed and still the same for all experiments
in this paper.

Based on this analysis, we address the estimate of resource consumption (called
% r2.11
silicon area -- in the case of FPGAs this means 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  \notag \\
\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 since we multiply some variable with another variable
and it is even non-quadratic, as the cost function $F$ does not have a known
linear or quadratic expression. To linearize this problem, we introduce $p$ FIR configurations.
% AH: conflit merge
% This variable must be defined by the user, it represent the number of different
% set of coefficients generated (for memory, we use \texttt{firls} and \texttt{fir1}
% functions from GNU Octave). To choose this value, we consider a subset of the figure~\ref{fig:rejection_pyramid}
% to restrict the number of configurations. Indeed, it is useless to have too many coefficients or
% too many bits, hence we take the configurations close to edge of pyramid. Thank to theses
% configurations $C_{ij}$ and $\pi_{ij}^C$ ($1 \leq j \leq p$) become constant
% and the function $F$ can be estimate for each configurations
% thanks our rejection criterion. We also defined binary
This variable $p$ is defined by the user, and represents the number of different
set of coefficients generated (remember, we use \texttt{firls} and \texttt{fir1}
functions from GNU Octave) based on the targeted filter characteristics and implementation
assumptions (estimated number of bits defining the coefficients). Hence, $C_{ij}$ and
$\pi_{ij}^C$ become constants and
we define $1 \leq j \leq p$ so that the function $F$ can be estimated (Look Up Table)
for each configurations thanks to the rejection criterion. We also define the 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.

% JM: conflict merge
% However the problem remains quadratic at this stage since in the constraint~\ref{eq:areadef2}
% we multiply
% $\delta_{ij}$ and $\pi_i^-$. However, since $\delta_{ij}$ is a binary variable we can
% linearise this multiplication if we can bound $\pi_i^-$. As $\pi_i^-$ is the data size,
% we define $0 < \pi_i^- \leq 128$ which is the maximum data size whose estimation is
% assumed on hardware characteristics.
% The Gurobi (\url{www.gurobi.com}) optimization software used to solve this quadratic
% model is able to linearize the model provided as is. This model
% has $O(np)$ variables and $O(n)$ constraints.}
The problem remains quadratic at this stage since in the constraint~\ref{eq:areadef2}
we multiply
$\delta_{ij}$ and $\pi_i^-$. However, since $\delta_{ij}$ is a binary variable we can
linearize this multiplication. The following formula shows how to linearize
this situation in general case with $y$ a binary variable and $x$ a real variable ($0 \leq x \leq X^{max}$):
\begin{equation*}
  m = x \times y \implies
  \left \{
  \begin{split}
    m & \geq 0 \\
    m & \leq y \times X^{max} \\
    m & \leq x \\
    m & \geq x - (1 - y) \times X^{max} \\
  \end{split}
  \right .
\end{equation*}
So if we bound up $\pi_i^-$ by 128~bits which is the maximum data size whose estimation is
assumed on hardware characteristics,
the Gurobi (\url{www.gurobi.com}) optimization software will be able to linearize
for us the quadratic problem so the model is left as is. This model
has $O(np)$ variables and $O(n)$ constraints.

% 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.
% % r2.12
% This variable must be defined by the user, it represent the number of different
% set of coefficients generated (for memory, we use \texttt{firls} and \texttt{fir1}
% functions from GNU Octave).
% 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.
%
% % r2.13
% This modified model is quadratic since we multiply two variables in the
% equation~\ref{eq:areadef2} ($\delta_{ij}$ by $\pi_{ij}^-$) but 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
silicon 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  \notag
\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]
    \node[draw,minimum size=1cm] (Solver) { Filter Solver } ;
    \node (Start) [left= 3cm of Solver] { } ;
    \node[draw,minimum size=1cm] (TCL) [right= of Solver] { TCL Script } ;
    \node (Input) [above= of TCL] { } ;
    \node[draw,minimum size=1cm] (Deploy) [below= of Solver] { Deploy Script } ;
    \node[draw,minimum size=1cm] (Bitstream) [below= of TCL] { Bitstream } ;
    \node[draw,minimum size=1cm,rounded corners] (Board) [below right= of Deploy] { Board } ;
    \node[draw,minimum size=1cm] (Postproc) [below= of Deploy] { Post-Processing } ;
    \node (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 allowing for
a fully automated optimal solution search.}
  \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}

% r2.14 et r2.15 et r2.16
\begin{figure}
  \centering
  \begin{subfigure}{\linewidth}
    \includegraphics[width=\linewidth]{images/max_500}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MAX/500 problem of maximizing rejection for a given resource allocation (500~arbitrary units).}
    \label{fig:max_500_result}
  \end{subfigure}

  \begin{subfigure}{\linewidth}
    \includegraphics[width=\linewidth]{images/max_1000}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MAX/1000 problem of maximizing rejection for a given resource allocation (1000~arbitrary units).}
    \label{fig:max_1000_result}
  \end{subfigure}

  \begin{subfigure}{\linewidth}
    \includegraphics[width=\linewidth]{images/max_1500}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MAX/1500 problem of maximizing rejection for a given resource allocation (1500~arbitrary units).}
    \label{fig:max_1500_result}
  \end{subfigure}
  \caption{Solutions for the MAX/500, MAX/1000 and MAX/1500 problems of maximizing
rejection for a given resource allocation.
The filter shape constraint (bandpass and bandstop) is shown as thick
horizontal lines on each chart.}
\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 following synthesis of the solutions found for
the problem of maximizing rejection for a given resource allocation. 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 arbitrary 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 composed of four 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, 80 or 100~dB.
Hence, the three cases have been named: MIN/40, MIN/60, MIN/80 and MIN/100.
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.
Table~\ref{tbl:gurobi_min_100} shows the results obtained by the filter solver for MIN/100.

\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}

\begin{table}[h!tb]
  \caption{Configurations $(C_i, \pi_i^C, \pi_i^S)$, rejections and areas (in arbitrary units) for MIN/100}
  \label{tbl:gurobi_min_100}
  \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 & -           & -           & -           & -           & -           & -               & -     \\
            2 & (15, 7, 17) & (51, 14, 0) & -           & -           & -           & 100~dB          & 1365  \\
            3 & (3, 3, 15)  & (27, 9, 0)  & (27, 9, 0)  & -           & -           & 100~dB          & 1002  \\
            4 & (3, 3, 15)  & (31, 9, 0)  & (19, 7, 0)  & (3, 3, 0)   & -           & 101~dB          & 909   \\
            5 & (3, 3, 15)  & (23, 8, 1)  & (19, 7, 0)  & (3, 3, 0)   & (3, 3, 0)   & 101~dB          & 810   \\
        \hline
      \end{tabular}
    }
\end{table}
\renewcommand{\arraystretch}{1}

From these tables, we can first state that almost 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. The only exception is for the monolithic case ($n = 1$) in
MIN/100: no solution is found for a single monolithic filter reach a 100~dB rejection.
Furthermore, 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.
Figure~\ref{fig:min_100} shows the rejection of the different configurations in the case of MIN/100.

% \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}
%
% \begin{figure}
% \centering
% \includegraphics[width=\linewidth]{images/min_100}
% \caption{Signal spectrum for MIN/100}
% \label{fig:min_100}
% \end{figure}

% r2.14 et r2.15 et r2.16
\begin{figure}
  \centering
  \begin{subfigure}{\linewidth}
    \includegraphics[width=.91\linewidth]{images/min_40}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MIN/40 problem of minimizing resource allocation for reaching a 40~dB rejection.}
    \label{fig:min_40}
  \end{subfigure}

  \begin{subfigure}{\linewidth}
    \includegraphics[width=.91\linewidth]{images/min_60}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MIN/60 problem of minimizing resource allocation for reaching a 60~dB rejection.}
    \label{fig:min_60}
  \end{subfigure}

  \begin{subfigure}{\linewidth}
    \includegraphics[width=.91\linewidth]{images/min_80}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MIN/80 problem of minimizing resource allocation for reaching a 80~dB rejection.}
    \label{fig:min_80}
  \end{subfigure}

  \begin{subfigure}{\linewidth}
    \includegraphics[width=.91\linewidth]{images/min_100}
    \caption{Filter transfer functions for varying number of cascaded filters solving
the MIN/100 problem of minimizing resource allocation for reaching a 100~dB rejection.}
    \label{fig:min_100}
  \end{subfigure}
  \caption{Solutions for the MIN/40, MIN/60, MIN/80 and MIN/100 problems of reaching a
given rejection while minimizing resource allocation. The filter shape constraint (bandpass and
bandstop) is shown as thick
horizontal lines on each chart.}
\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 (except in MIN/100 which has no monolithic solution) or cascaded filters.

Table~\ref{tbl:resources_usage} shows the resource usage in the case of MIN/40, MIN/60;
MIN/80 and MIN/100 \emph{i.e.} when the target rejection is fixed to 40, 60, 80 and 100~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.

\renewcommand{\arraystretch}{1.2}
\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
  {\scalefont{0.90}
      \begin{tabular}{|c|c|cccc|c|}
        \hline
        $n$ &          & MIN/40   & MIN/60   & MIN/80   & MIN/100  & \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     & 640      & \emph{17600}             \\
        2   & BRAM     & 2        & 2        & 2        & 2        & \emph{120}               \\
            & DSP      & 0        & 0        & 0        & 66       & \emph{80}                \\ \hline
            & LUT      & 891      & 2148     & 2023     & 2448     & \emph{17600}             \\
        3   & BRAM     & 3        & 3        & 3        & 3        & \emph{120}               \\
            & DSP      & 0        & 0        & 19       & 27       & \emph{80}                \\ \hline
            & LUT      & 662      & 1729     & 2451     & 2893     & \emph{17600}             \\
        4   & BRAM     & 4        & 4        & 4        & 4        & \emph{120}               \\
            & DPS      & 0        & 0        & 7        & 19       & \emph{80}                \\ \hline
            & LUT      & -        & 1259     & 2602     & 2505     & \emph{17600}             \\
        5   & BRAM     & -        & 5        & 5        & 5        & \emph{120}               \\
            & DPS      & -        & 0        & 0        & 19       & \emph{80}                \\ \hline
      \end{tabular}
  }
\end{table}
\renewcommand{\arraystretch}{1}

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.

\renewcommand{\arraystretch}{1.2}
\begin{table}[h!tb]
\caption{Time to solve the quadratic program with Gurobi}
\label{tbl:area_time_comp}
\centering
{\scalefont{0.90}
\begin{tabular}{|c|c|c|c|c|}\hline
$n$ & Time (MIN/40)           & Time (MIN/60)               & Time (MIN/80) & Time (MIN/100)               \\\hline\hline
1   & 0.07~s                  & 0.02~s                      & 0.01~s        & -                            \\
2   & 7.8~s                   & 16~s                        & 14~s          & 1.8~s                        \\
3   & 4.7~s                   & 14~s                        & 28~s          & 39~s                         \\
4   & 39~s                    & 20~s                        & 193~s         & 522~s  ($\approx$ 9~min)     \\
5   & -                       & 12~s                        & 170~s         & 1048~s ($\approx$ 17~min)    \\\hline
\end{tabular}
}
\end{table}
\renewcommand{\arraystretch}{1}

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 17~minutes,
compared to 3~days in the previous section: this problem is more easily solved than the
previous one.

To conclude, we compare our monolithic filters with the FIR Compiler provided by
Xilinx in the Vivado software suite (v.2018.2). For each experiment we use the
same coefficient set and we compare the resource consumption, having checked that
the transfer functions are indeed the same with both implementations.
Table~\ref{tbl:xilinx_resources} exhibits the results.
The FIR Compiler never uses BRAM while our filter implementation uses one block. This difference
is explained be our wish to have a dynamically reconfigurable FIR filter whose
coefficients can be updated from the processing system without having to update the FPGA design.
With the FIR compiler, the coefficients are defined during the FPGA design so that
changing coefficients required generating a new design. The difference with the LUT consumption
is also attributed to the reconfigurability logic. However the DSP consumption, the scarcest
resource, is the same between the Xilinx FIR Compiler end
our FIR block: we hence conclude that our solutions are as good as the Xilinx implementation.

\renewcommand{\arraystretch}{1.2}
\begin{table}
\centering
\caption{Resource consumption compared between the FIR Compiler from Xilinx and our FIR block}
\label{tbl:xilinx_resources}
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
\multirow{2}{*}{} & \multicolumn{3}{c|}{Xilinx} & \multicolumn{3}{c|}{Our FIR block} \\ \cline{2-7}
                  & LUT     & BRAM     & DSP    & LUT       & BRAM       & DSP       \\ \hline
MAX/500           & 177     & 0        & 21     & 249       & 1          & 21        \\ \hline
MAX/1000          & 306     & 0        & 37     & 453       & 1          & 37        \\ \hline
MAX/1500          & 418     & 0        & 47     & 627       & 1          & 47        \\ \hline
MIN/40            & 225     & 0        & 27     & 347       & 1          & 27        \\ \hline
MIN/60            & 322     & 0        & 39     & 334       & 1          & 39        \\ \hline
MIN/80            & 482     & 0        & 55     & 772       & 1          & 55        \\ \hline
\end{tabular}
\end{table}
\renewcommand{\arraystretch}{1}

\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}