% 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)

% ajouter pyramide "juste"
% 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}
\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.
\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 filter 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 offer 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}
We want create a new methodology to develop any Digital Signal Processing (DSP) chain
and for any hardware platform (Altera, Xilinx...). To do this we have defined an
abstract model to represent some basic operations of DSP.

For the moment, we are focused on only two operations: the filtering and the shifting of data.
We have chosen this basic operation because the shifting and the filtering have already be studied in
lot of works \cite{lim_1996, lim_1988, young_1992, smith_1998} hence it will be easier
to check and validate our results.

However having only two operations is insufficient to work with complex DSP but
in this paper we only want demonstrate the relevance and the efficiency of our approach.
In future work it will be possible to add more operations and we are able to
model any DSP chain.

We will apply our methodology on very simple DSP chain. We generate a digital signal
thanks at generator of Pseudo-Random Number (PRN) or thanks at an Analog to Digital
Converter (ADC). Once we have a digital signal, we filter it to decrease the noise level.
Finally we stored some burst of filtered samples before post-processing it.
In this particular case, we want optimize the filtering step to have the best noise
rejection for constrain number of resource or to have the minimal resources
consumption for a given rejection objective.

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 filter are composed of $n$ stage. In stage $i$ ($1 \leq i \leq n$)
the FIR has $C_i$ coefficients and each coefficients are integer values with $\pi^C_i$
bits and the filtered data are shifted of $\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$ can reject $F(C_i, \pi_i^C)$ dB. $F$ is determined numerically.
To measure this rejection, we use GNU Octave software to design FIR filter coefficients thanks to two
algorithms (\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 other are coded on very fewer bits.

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

\begin{figure}
  \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}
\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.
Our first criterion considers the mean value of the stopband rejection, as shown in figure~\ref{fig:mean_criterion}. This criterion does not work because we do not consider the shape of the passband.
A second criterion considers 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}.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/colored_mean_criterion}
\caption{Mean criterion comparison between monolithic filter and cascade filters}
\label{fig:mean_criterion}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/colored_custom_criterion}
\caption{Custom criterion comparison between monolithic filter and cascade filters}
\label{fig:custom_criterion}
\end{figure}

Thanks to this criterion we are able to automatically generate lot of fir coefficients
and estimate their rejection. The figure~\ref{fig:rejection_pyramid} exhibits the
rejection in function of the number of coefficients and their number of bits.
We can observe it looks like a pyramid so the edge represents the best
coefficient set. Indeed if we choose a number of coefficients, increasing the number
of bits over the edge will not improve the rejection. Conversely when we choose
a number of bits, too much increase the number of coefficients will not improve
the rejection. Hence the best coefficient set are on the edge of 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 a efficient criterion to estimate the rejection of one set of coefficient
we have a problem when we sum two or more criterion. If the FIR filter coefficients are the same
between the stage, we have:
$$F_{total} = F_1 + F_2$$
But when we choose two different set of coefficient, the previous equality are not
true. The figure~\ref{fig:sum_rejection} illustrates the problem. The red and blue curves
are two different filter coefficient and we can see that their maximum on the stopband
are not at the same frequency. So when we sum the rejection criteria (the dotted yellow line)
we do not meet the dashed yellow line. Define the rejection of cascaded filters
is more difficult than just take the summation between all the rejection criteria of each filter.
However this summation gives us an upper bound for rejection although in fact we obtain
better rejection than expected.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{images/cascaded_criterion}
\caption{Rejection of two cascaded filters}
\label{fig:sum_rejection}
\end{figure}

The first problem we address is to maximize the rejection under bounded silicon area
and feasibility constraints. 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:

Finally we can describe our abstract model with following expressions :
\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 for a filter. More precisely, it is the area of the FIR as the Shifter
does not need any circuitry. 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 and the coefficients. Equation~\ref{eq:rejectiondef} gives the
definition of the rejection of the filter thanks to function~$F$ that we defined
previously. The Shifter does not introduce negative rejection as we 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 of 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 6~dB of noise. 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 global input's number of bits.

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.

The section~\ref{sec:fixed_area} shows the results for the first version of quadratic program but the section~\ref{sec:fixed_rej}
presents the results for the complementary problem. In this case we want
minimize the occupied area for a targeted rejection level. Hence we have replace
the objective function with:
\begin{align}
\text{Minimize } & \sum_{i=1}^n a_i  \notag
\end{align}
We adapt our constraints of quadratic program to replace the equation \ref{eq:area}
by the 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 section~\ref{sec:fixed_area}.
Figure~\ref{fig:workflow} shows the global workflow and the different steps involved in the computations 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}
  \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 get 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).
The raw input data generated from a 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.
Whereas the signal is processed in real-time, the output signal is stored as
consecutive bursts of data.

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 125~MS/s ADC.
The board works with a Buildroot Linux image. We have developed some tools and
drivers to flash and communicate with the FPGA. They are used to automatize all
the workflow inside the board: load the filter coefficients and retrieve the
computed data.

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.

The workflow used to compute the results in section~\ref{sec:fixed_rej}, we
have just adapted the quadratic program but the rest of the workflow is unchanged.

\section{Experiments with fixed area space}
\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.
This is interesting to understand 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 conclusion
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 (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 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 level
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 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}
  \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 filters coefficients are small enough, or when the input size is small
enough, Vivado optimized 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 corresponds to 2500 LUTs,
1000 arbitrary units corresponds to 5000 LUTs and 1500 arbitrary units corresponds
to 7300 LUTs. The conclusion is that the orders of magnitude of our arbitrary
unit are quite good. The relatively small differences can probably be explained
by the optimizations done by Vivado based on the detailed map of available processing resources.

We present the computation time to solve the quadratic problem.
For each case, the filter solver software are executed with a Intel(R) Xeon(R) CPU E5606
cadenced 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}
\caption{Time 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. On the contrary, in the case of MAX/1500 with
5~stages, we were not able to obtain a result after 40~hours of computation so we decided to stop.

\section{Experiments with fixed rejection target}
\label{sec:fixed_rej}
This section presents the results of complementary quadratic program which we
minimize the area occupation for a targeted noise level.

The experimental setup is also composed of three cases. The raw input is the same
as previous section, 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}
  \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}
  \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}
  \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}

From these tables, we can first state that all configuration reach the target rejection
level and more we have stages lesser is the area occupied in arbitrary unit.
Futhermore, the area of the monolithic filter is twice bigger than the two cascaded.
More generally, more there is filters lower is the occupied area.

Like in previous section, the solver choose always a little filter as first
filter stage and the second one is often the biggest filter. this choice can be explain
as the previous section. The solver uses just enough bits to not degrade the input
signal and in second filter it can choose a better filter to improve rejection without
have too bits in the output data.

For the specific case in MIN/40 for $n = 5$ the solver has determined that the optimal
number of filter is 4 so it not chose any configuration in 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 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 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 real
rejection. All curves prove that the constraint to reach the target rejection is
respected both monolithic filter or cascaded filters.

Table~\ref{tbl:resources_usage} shows the resources 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 term of LUT (1 DSP $\approx$ 100 LUT)
the real resource consumption decrease in function of number of stage 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 lesser than expected.

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

\begin{table}
\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 are substantially faster than time
needed in the previous section. Indeed the worst time in this case is only 3~minutes
in balance of 3~days on previous section. We are able to solve more easily this
problem than the previous one.

\section{Conclusion}

In this paper, we have proposed a new approach to work with a cascade of FIR filter inside a FPGA.
This method aims to be hardware independent and focus an high-level of abstraction.
We have modeled the FIR filter operation and the data shift impact. With this model
we have created a quadratic program to select the optimal FIR coefficient set to reject a
maximum of noise. In our experiments we have chosen deliberately some common tools
to design the filter coefficients but we can use any other method.

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 our
workflow is able to filter an input signal as expected which validates our model and our approach.
We can easily change the quadratic program to adapt it to an other problem.

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 so trivial, especially in terms of silicon
area for the 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}