jfriedt / IFCS2018 article

Browse Code »

Commit 842e804be4a295da293fe60a2399ce6987d0a453

Authored by Arthur HUGEAT 2019-01-14 16:51:29 +0100

1 parent 27f5f41088

Exists in master

Permier pas vers l'article journal.

Showing 11 changed files with 382 additions and 493 deletions Inline Diff

demo_critere_filtre.m
ifcs2018_journal.tex
images/custom_criterion.pdf
images/max_rejection/prn_1000.pdf
images/max_rejection/prn_2000.pdf
images/max_rejection/prn_500.pdf
images/mean_criterion.pdf
images/min_area/prn_100.pdf
images/min_area/prn_150.pdf
images/min_area/prn_50.pdf
images/sum_rejection.pdf

demo_critere_filtre.m

Diff comments View file @ 842e804

File was created	1	clear all;
	2
	3	global N = 2048;
	4
	5	function rejection = rejection_criteria(log_data, fc)
	6	N = length(log_data);
	7
	8	% Index of the first point in the tail fir
	9	index_tail = round((fc + 0.1) * N) + 1;
	10
	11	% Index of th last point on the band
	12	index_band = round((fc - 0.1) * N);
	13
	14	% Get the worst rejection in stopband
	15	worst_rejection = -max(log_data(index_tail:end));
	16
	17	% Get the total of deviation in passband
	18	worst_band = mean(-1 * abs(log_data(1:index_band)));
	19
	20	% Compute the rejection
	21	passband_malus = 10; % weighted value to penalize the deviation in passband
	22	rejection = worst_band * passband_malus + worst_rejection;
	23	endfunction
	24
	25	function [h, log_curve, rejection] = compute_freqz(filename)
	26	global N;
	27
	28	b = load(filename);
	29	[h, w] = freqz(b, 1, N/2);
	30	mag = abs(h);
	31	mag = mag ./ mag(1);
	32	log_curve = 20 * log10(mag);
	33
	34	rejection = rejection_criteria(log_curve, 0.5);
	35	endfunction
	36
	37	# Stages
	38	hTotal = ones(N/2, 1);
	39	% [h, curve1, c1] = compute_freqz("filters/fir1/fir1_033_int08"); % 1) -8dB
	40	% [h, curve1, c1] = compute_freqz("filters/fir1/fir1_037_int08"); % 2) -9dB
	41	[h, curve1, c1] = compute_freqz("filters/fir1/fir1_037_int08");
	42	hTotal = hTotal .* h;
	43	% [h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10"); % 1) -8dB
	44	% [h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10"); % 2) -9dB
	45	[h, curve2, c2] = compute_freqz("filters/fir1/fir1_033_int10");
	46	hTotal = hTotal .* h;
	47	% [h, curve3, c3] = compute_freqz("filters/fir1/fir1_033_int08");
	48	% hTotal = hTotal .* h;
	49	% [h, curve4, c4] = compute_freqz("filters/fir1/fir1_033_int10");
	50	% hTotal = hTotal .* h;
	51	% [h, curve5, c5] = compute_freqz("filters/fir1/fir1_015_int11");
	52	% hTotal = hTotal .* h;
	53
	54	# Log total
	55	mag = abs(hTotal);
	56	mag = mag ./ mag(1);
	57	log_freqz = 20 * log10(mag);
	58	cTotal = rejection_criteria(log_freqz, 0.5);
	59
	60	[ c1+c2 cTotal ]
	61	% [ c1+c2+c3+c4+c5 cTotal ]
	62
	63	clf;
	64	f_axe = [1:N/2] * 2/N;
	65	hold on;
	66	color = [0/255 114/255 189/255];
	67	plot(f_axe, curve1, "linewidth", 1.5, "color", color);
	68	plot([0 1], [-c1 -c1], "--", "linewidth", 1.5, "color", color);
	69
	70	color = [217/255 83/255 25/255];
	71	plot(f_axe, curve2, "linewidth", 1.5, "color", color);

ifcs2018_journal.tex

Diff comments View file @ 842e804

% JMF : revoir l'abstract : on y avait mis le Zynq7010 de la redpitaya en montrant	1
% comment optimiser les perfs a surface finie. Ici aussi on tombait dans le cas ou`	2
% la solution a 1 seul FIR n'etait simplement pas synthetisable => fusionner les deux	3
% contributions pour le papier TUFFC	4
	5
\documentclass[a4paper,conference]{IEEEtran/IEEEtran}	6	1	\documentclass[a4paper,conference]{IEEEtran/IEEEtran}
\usepackage{graphicx,color,hyperref}	7	2	\usepackage{graphicx,color,hyperref}
\usepackage{amsfonts}	8	3	\usepackage{amsfonts}
\usepackage{amsthm}	9	4	\usepackage{amsthm}
\usepackage{amssymb}	10	5	\usepackage{amssymb}
\usepackage{amsmath}	11	6	\usepackage{amsmath}
\usepackage{algorithm2e}	12	7	\usepackage{algorithm2e}
\usepackage{url,balance}	13	8	\usepackage{url,balance}
\usepackage[normalem]{ulem}	14	9	\usepackage[normalem]{ulem}
		10	\usepackage{tikz}
		11	\usetikzlibrary{positioning,fit}
		12	\usepackage{multirow}
		13	\usepackage{scalefnt}
		14
% correct bad hyphenation here	15	15	% correct bad hyphenation here
\hyphenation{op-tical net-works semi-conduc-tor}	16	16	\hyphenation{op-tical net-works semi-conduc-tor}
\textheight=26cm	17	17	\textheight=26cm
\setlength{\footskip}{30pt}	18	18	\setlength{\footskip}{30pt}
\pagenumbering{gobble}	19	19	\pagenumbering{gobble}
\begin{document}	20	20	\begin{document}
\title{Filter optimization for real time digital processing of radiofrequency signals: application	21	21	\title{Filter optimization for real time digital processing of radiofrequency signals: application
to oscillator metrology}	22	22	to oscillator metrology}
	23	23
\author{\IEEEauthorblockN{A. Hugeat\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, J. Bernard\IEEEauthorrefmark{2},	24	24	\author{\IEEEauthorblockN{A. Hugeat\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, J. Bernard\IEEEauthorrefmark{2},
G. Goavec-M\'erou\IEEEauthorrefmark{1},	25	25	G. Goavec-M\'erou\IEEEauthorrefmark{1},
P.-Y. Bourgeois\IEEEauthorrefmark{1}, J.-M. Friedt\IEEEauthorrefmark{1}}	26	26	P.-Y. Bourgeois\IEEEauthorrefmark{1}, J.-M. Friedt\IEEEauthorrefmark{1}}
\IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, Time \& Frequency department, Besan\c con, France }	27	27	\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 \\	28	28	\IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Computer Science department DISC, Besan\c con, France \\
Email: \{pyb2,jmfriedt\}@femto-st.fr}	29	29	Email: \{pyb2,jmfriedt\}@femto-st.fr}
}	30	30	}
\maketitle	31	31	\maketitle
\thispagestyle{plain}	32	32	\thispagestyle{plain}
\pagestyle{plain}	33	33	\pagestyle{plain}
\newtheorem{definition}{Definition}	34	34	\newtheorem{definition}{Definition}
	35	35
\begin{abstract}	36	36	\begin{abstract}
Software Defined Radio (SDR) provides stability, flexibility and reconfigurability to	37	37	Software Defined Radio (SDR) provides stability, flexibility and reconfigurability to
radiofrequency signal processing. Applied to oscillator characterization in the context	38	38	radiofrequency signal processing. Applied to oscillator characterization in the context
of ultrastable clocks, stringent filtering requirements are defined by spurious signal or	39	39	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	40	40	noise rejection needs. Since real time radiofrequency processing must be performed in a
Field Programmable Array to meet timing constraints, we investigate optimization strategies	41	41	Field Programmable Array to meet timing constraints, we investigate optimization strategies
to design filters meeting rejection characteristics while limiting the hardware resources	42	42	to design filters meeting rejection characteristics while limiting the hardware resources
required and keeping timing constraints within the targeted measurement bandwidths.	43	43	required and keeping timing constraints within the targeted measurement bandwidths.
\end{abstract}	44	44	\end{abstract}
	45	45
\begin{IEEEkeywords}	46	46	\begin{IEEEkeywords}
Software Defined Radio, Mixed-Integer Linear Programming, Finite Impulse Response filter	47	47	Software Defined Radio, Mixed-Integer Linear Programming, Finite Impulse Response filter
\end{IEEEkeywords}	48	48	\end{IEEEkeywords}
	49	49
\section{Digital signal processing of ultrastable clock signals}	50	50	\section{Digital signal processing of ultrastable clock signals}
	51	51
Analog oscillator phase noise characteristics are classically performed by downconverting	52	52	Analog oscillator phase noise characteristics are classically performed by downconverting
the radiofrequency signal using a saturated mixer to bring the radiofrequency signal to baseband,	53	53	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	54	54	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	55	55	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}.	56	56	multiplying the samples with a local numerically controlled oscillator (Fig. \ref{schema}) \cite{rsi}.
	57	57
\begin{figure}[h!tb]	58	58	\begin{figure}[h!tb]
\begin{center}	59	59	\begin{center}
\includegraphics[width=.8\linewidth]{images/schema}	60	60	\includegraphics[width=.8\linewidth]{images/schema}
\end{center}	61	61	\end{center}
\caption{Fully digital oscillator phase noise characterization: the Device Under Test	62	62	\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	63	63	(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	64	64	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	65	65	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	66	66	Impulse Response (FIR) filters. The signal is then decimated before a Fourier analysis displays
the spectral characteristics of the phase fluctuations.}	67	67	the spectral characteristics of the phase fluctuations.}
\label{schema}	68	68	\label{schema}
\end{figure}	69	69	\end{figure}
	70	70
As with the analog mixer,	71	71	As with the analog mixer,
the non-linear behavior of the downconverter introduces noise or spurious signal aliasing as	72	72	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.	73	73	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	74	74	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	75	75	for the phase noise spectral characterization \cite{andrich2018high}. The characteristics introduced between the
downconverter	76	76	downconverter
and the decimation processing blocks are core characteristics of an oscillator characterization	77	77	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	78	78	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	79	79	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	80	80	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	81	81	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	82	82	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	83	83	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	84	84	tunable number of coefficients and tunable number of bits representing the coefficients and the
data being processed.	85	85	data being processed.
	86	86
\section{Finite impulse response filter}	87	87	\section{Finite impulse response filter}
	88	88
We select FIR filter for their unconditional stability and ease of design. A FIR filter is defined	89	89	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	90	90	by a set of weights $b_k$ applied to the inputs $x_k$ through a convolution to generate the
outputs $y_k$	91	91	outputs $y_k$
$$y_n=\sum_{k=0}^N b_k x_{n-k}$$	92	92	\begin{align}
		93	y_n=\sum_{k=0}^N b_k x_{n-k}
		94	\label{eq:fir_equation}
		95	\end{align}
	93	96
As opposed to an implementation on a general purpose processor in which word size is defined by the	94	97	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	95	98	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	96	99	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	97	100	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	98	101	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	99	102	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.	100	103	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	101	104	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,	102	105	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.	103	106	is not considered as an issue as would be in a closed loop system.
	104	107
The coefficients are classically expressed as floating point values. However, this binary	105	108	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,	106	109	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	107	110	we select to quantify these floating point values into integer values. This quantization
will result in some precision loss.	108	111	will result in some precision loss.
	109	112
%As illustrated in Fig. \ref{float_vs_int}, we see that we aren't	110
%need too coefficients or too sample size. If we have lot of coefficients but a small sample size,	111
%the first and last are equal to zero. But if we have too sample size for few coefficients that not improve the quality.	112
	113
% JMF je ne comprends pas la derniere phrase ci-dessus ni la figure ci dessous	114
% AH en gros je voulais dire que prendre trop peu de bit avec trop de coeff, ça induit ta figure (bien mieux faite que moi)	115
% et que l'inverse trop de bit sur pas assez de coeff on ne gagne rien, je vais essayer de la reformuler	116
	117
%\begin{figure}[h!tb]	118
%\includegraphics[width=\linewidth]{images/float-vs-integer.pdf}	119
%\caption{Impact of the quantization resolution of the coefficients}	120
%\label{float_vs_int}	121
%\end{figure}	122
	123
\begin{figure}[h!tb]	124	113	\begin{figure}[h!tb]
\includegraphics[width=\linewidth]{images/demo_filtre}	125	114	\includegraphics[width=\linewidth]{images/demo_filtre}
\caption{Impact of the quantization resolution of the coefficients: the quantization is	126	115	\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	127	116	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	128	117	the 30~first and 30~last coefficients out of the initial 128~band-pass
filter coefficients to 0 (red dots).}	129	118	filter coefficients to 0 (red dots).}
\label{float_vs_int}	130	119	\label{float_vs_int}
\end{figure}	131	120	\end{figure}
	132	121
The tradeoff between quantization resolution and number of coefficients when considering	133	122	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	134	123	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	135	124	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	136	125	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	137	126	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	138	127	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	139	128	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	140	129	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	141	130	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	142	131	and tap length, as will be shown in the next section. Indeed, our development strategy closely
follows the skeleton approach \cite{crookes1998environment, crookes2000design, benkrid2002towards}	143	132	follows the skeleton approach \cite{crookes1998environment, crookes2000design, benkrid2002towards}
in which basic blocks are defined and characterized before being assembled \cite{hide}	144	133	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	145	134	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	146	135	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	147	136	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	148	137	current implementation: the decimation is assumed to be located after the FIR cascade at the
moment.	149	138	moment.
	150	139
\section{Filter optimization}	151	140	\section{Methodology description}
		141	We want create a new methodology to develop any Digital Signal Processing (DSP) chain
		142	and for any hardware platform (Altera, Xilinx...). To do this we have defined an
		143	abstract model to represent some basic operations of DSP.
	152	144
A basic approach for implementing the FIR filter is to compute the transfer function of	153	145	For the moment, we are focused on only two operations: the filtering and the shift of data.
a monolithic filter: this single filter defines all coefficients with the same resolution	154	146	We have chosen this basic operation because the shifting and the filtering have already be studied in
(number of bits) and processes data represented with their own resolution. Meeting the	155	147	lot of works {\color{red} mettre les nouvelles référence ici} hence it will be easier
filter shape requires a large number of coefficients, limited by resources of the FPGA since	156	148	to check and validate our results.
this filter must process data stream at the radiofrequency sampling rate after the mixer.	157
	158	149
An optimization problem \cite{leung2004handbook} aims at improving one or many	159	150	However having only two operations is insufficient to work with complex DSP but
performance criteria within a constrained resource environment. Amongst the tools	160	151	in this paper we only want demonstrate the relevance and the efficiency of our approach.
developed to meet this aim, Mixed-Integer Linear Programming (MILP) provides the framework to	161	152	In future work it will be possible to add more operations and we are able to
formally define the stated problem and search for an optimal use of available	162	153	model any DSP chain.
resources \cite{yu2007design, kodek1980design}.	163
	164	154
First we need to ensure that our problem is a real optimization problem. When	165	155	We will apply our methodology on very simple DSP chain. We generate a digital signal
designing a processing function in the FPGA, we aim at meeting some requirement such as	166	156	thanks at generator of Pseudo-Random Number (PRN) or thanks at an Analog to Digital
the throughput, the computation time or the noise rejection noise. However, due to limited	167	157	Converter (ADC). Once we have a digital signal, we filter it to decrease the noise level.
resources to design the process like BRAM (high performance RAM), DSP (Digital Signal Processor)	168	158	Finally we stored some burst of filtered samples before post-processing it.
or LUT (Look Up Table), a tradeoff must be generally searched between performance and available	169	159	% TODO: faire un schéma
computational resources: optimizing some criteria within finite, limited	170	160	In this particular case, we want optimize the filtering step to have the best noise
resources indeed matches the definition of a classical optimization problem.	171	161	rejection for constrain number of resource or to have the minimal resources
		162	consumption for a given rejection objective.
	172	163
Specifically the degrees of freedom when addressing the problem of replacing the single monolithic	173	164	The first step of our approach is to model the DSP chain and since we just optimize
FIR with a cascade of optimized filters are the number of coefficients $N_i$ of each filter $i$,	174	165	the filtering, we have not modeling the PRN generator or the ADC. The filtering can be
the number of bits $C_i$ representing the coefficients and the number of bits $D_i$ needed to represent	175	166	done by two ways. The first one we use only one FIR filter with lot of coefficients
the data $x_k$ fed to each filter as provided by the acquisition or previous processing stage.	176	167	to rejection the noise, we called this approach a monolithic approach. And the second one
Because each FIR in the chain is fed the output of the previous stage,	177	168	we select different FIR filters with less coefficients the monolithic filter and we cascaded
the optimization of the complete processing chain within a constrained resource environment is not	178	169	it to filtering the signal.
trivial. The resource occupation of a FIR filter is considered as $C_i \times N_i$ which aims	179
at approximating the number of bits needed in a worst case condition to represent the output of the	180
FIR. Indeed, the number of bits generated by the $i$th FIR is $(C_i+D_i)\times\log_2(N_i)$, but the	181
$\log$ function is avoided for its incompatibility with a linear programming description, and	182
the simple product is approximated as the number of gates needed to perform the calculation. Such an	183
occupied area estimate assumes that the number of gates scales as the number of bits and the number	184
of coefficients, but does not account for the detailed implementation of the hardware. Indeed,	185
various FPGA implementations will provide different hardware functionalities, and we shall consider	186
at the end of the design a synthesis step using vendor software to assess the validity of the solution	187
found. As an example of the limitation linked to the lack of detailed hardware consideration, Block Random	188
Access Memory (BRAM) used to store filter coefficients are not shared amongst filters, and multiplications	189
are most efficiently implemented by using DSP blocks whose input word	190
size is finite. DSPs are a scarce resource to be saved in a practical implementation. Keeping a high	191
abstraction on the resource occupation is nevertheless selected in the following discussion in order	192
to leave enough degrees of freedom in the problem to try and find original solutions: too many	193
constraints in the initial statement of the problem leave little room for finding an optimal solution.	194
	195	170
\begin{figure}[h!tb]	196	171	After each filter we leave the possibility of shifting the filtered data to consume
\begin{center}	197	172	less resources. Hence in the case of cascaded filter, we define a stage as a filter
\includegraphics[width=.5\linewidth]{schema2}	198	173	and a shifter (the shift could be omitted if we do not need to divide the filtered data).
\caption{Shape of the filter transmitted power $P$ as a function of frequency:	199
the bandpass BP is considered to occupy the initial	200
40\% of the Nyquist frequency range, the stopband the last 40\%, allowing 20\% transition	201
width.}	202
\label{rejection-shape}	203
\end{center}	204
\end{figure}	205
	206	174
Following these considerations, the model is expressed as:	207	175	\subsection{Model of a FIR filter}
\begin{align}	208	176	A cascade of filter are composed of $n$ stage. In stage $i$ ($1 \leq i \leq n$)
\begin{cases}	209	177	the FIR has $C_i$ coefficients and each coefficients are integer values with $\pi^C_i$
\mathcal{R}_i &= \mathcal{F}(N_i, C_i)\\	210	178	bits and the filtered data are shifted of $\pi^S_i$ bits. We define also $\pi^-_i$ as
\mathcal{A}_i &= N_i \times C_i\\	211	179	the size of input data and $\pi^+_i$ as the size of output data. The figure~\ref{fig:fir_stage}
\Delta_i &= \Delta _{i-1} + \mathcal{P}_i	212	180	shows a filtering stage.
\end{cases}	213
\label{model-FIR}	214
\end{align}	215
To explain the system \ref{model-FIR}, $\mathcal{R}_i$ represents the stopband rejection dependence with $N_i$ and $C_i$, $\mathcal{A}_i$	216
is a theoretical area occupation of the processing block on the FPGA as discussed earlier, and $\Delta_i$ is the total rejection for the current stage $i$.	217
Since the function $\mathcal{F}$ cannot be explictly expressed, we run simulations to determine the rejection depending	218
on $N_i$ and $C_i$. However, selecting the right filter requires a clear definition of the rejection criterion. Selecting an	219
incorrect criterion will lead the linear program solver to produce a solution which might not meet the user requirements.	220
Hence, amongst various criteria including the mean or median value of the FIR response in the stopband as will	221
be illustrated lated (section \ref{median}), we have designed	222
a criterion aimed at avoiding ripples in the passband and considering the maximum of the FIR spectral response in the stopband	223
(Fig. \ref{rejection-shape}). The bandpass criterion is defined as the sum of the absolute values of the spectral response	224
in the bandpass, reminiscent of a standard deviation of the spectral response: this criterion must be minimized to avoid	225
ripples in the passband. The stopband transfer function maximum must also be minimized in order to improve the filter	226
rejection capability. Weighing these two criteria allows designing the linear program to be solved.	227
	228	181
\begin{figure}[h!tb]	229	182	\begin{figure}
\includegraphics[width=\linewidth]{images/noise-rejection.pdf}	230	183	\centering
\caption{Rejection as a function of number of coefficients and number of bits}	231	184	\begin{tikzpicture}[node distance=2cm]
\label{noise-rejection}	232	185	\node[draw,minimum size=1.3cm] (FIR) { $C_i, \pi_i^C$ } ;
\end{figure}	233	186	\node[draw,minimum size=1.3cm] (Shift) [right of=FIR, ] { $\pi_i^S$ } ;
		187	\node (Start) [left of=FIR] { } ;
		188	\node (End) [right of=Shift] { } ;
	234	189
The objective function maximizes the noise rejection ($\max(\Delta_{i_{\max}})$) while keeping resource	235	190	\node[draw,fit=(FIR) (Shift)] (Filter) { } ;
occupation below a user-defined threshold, or as will be discussed here, aims at minimizing the area	236
needed to reach a given rejection ($\min(S_q)$ in the forthcoming discussion, Eqs. \ref{cstr_size}	237
and \ref{cstr_rejection}). The MILP solver is allowed to choose the number of successive	238
filters, within an upper bound. The last problem is to model the noise rejection. Since filter	239
noise rejection capability is not modeled with linear equations, a look-up-table is generated	240
for multiple filter configurations in which the $C_i$, $D_i$ and $N_i$ parameters are varied: for each	241
one of these conditions, the low-pass filter rejection is stored as computed by the frequency response	242
of the digital filter (Fig. \ref{noise-rejection}). Various rejection criteria have been investigated,	243
including mean value of the stopband response, median value of the stopband response, or as finally	244
selected, maximum value in the stopband. An intuitive analysis of the chart of Fig. \ref{noise-rejection}	245
hints at an optimum	246
set of tap length and number of bit for representing the coefficients along the line of the pyramidal	247
shaped rejection capability function.	248
	249	191
Linear program formalism for solving the problem is well documented: an objective function is	250	192	\draw[->] (Start) edge node [above] { $\pi_i^-$ } (FIR) ;
defined which is linearly dependent on the parameters to be optimized. Constraints are expressed	251	193	\draw[->] (FIR) -- (Shift) ;
as linear equations and solved using one of the available solvers, in our case GLPK\cite{glpk}.	252	194	\draw[->] (Shift) edge node [above] { $\pi_i^+$ } (End) ;
With the notations used in the description of system \ref{model-FIR}, we have defined the linear problem as:	253	195	\end{tikzpicture}
\paragraph{Variables}	254	196	\caption{A single filter is composed of a FIR (on the left) and a Shifter (on the right)}
\begin{align*}	255	197	\label{fig:fir_stage}
x_{i,j} \in \lbrace 0,1 \rbrace & \text{ $i$ is a given filter} \\	256	198	\end{figure}
& \text{ $j$ is the stage} \\	257
& \text{ If $x_{i,j}$ is equal to 1, the filter is selected} \\	258
\end{align*}	259
\paragraph{Constants}	260
\begin{align*}	261
\mathcal{F} = \lbrace F_1 ... F_p \rbrace & \text{ All possible filters}\\	262
& \text{ $p$ is the number of different filters} \\	263
% N(i) & \text{ % Constant to let the	264
% number of coefficients %} \\ & \text{	265
% for filter $i$}\\	266
% C(i) & \text{ % Constant to let the	267
% number of bits of %}\\ & \text{	268
% each coefficient for filter $i$}\\	269
\mathcal{S}_{\max} & \text{ Total space available inside the FPGA}	270
\end{align*}	271
\paragraph{Constraints}	272
\begin{align}	273
1 \leq i \leq p & \nonumber\\	274
1 \leq j \leq q & \text{ $q$ is the max of filter stage} \nonumber \\	275
\forall j, \mathlarger{\sum_{i}} x_{i,j} = 1 & \text{ At most one filter by stage} \nonumber\\	276
\mathcal{S}_0 = 0 & \text{ initial occupation} \nonumber\\	277
\forall j, \mathcal{S}_j = \mathcal{S}_{j-1} + \mathlarger{\sum_i (x_{i,j} \times \mathcal{A}_i)} \label{cstr_size} \\	278
\mathcal{S}_j \leq \mathcal{S}_{\max}\nonumber \\	279
\mathcal{N}_0 = 0 & \text{ initial rejection}\nonumber\\	280
\forall j, \mathcal{N}_j = \mathcal{N}_{j-1} + \mathlarger{\sum_i (x_{i,j} \times \mathcal{R}_i)} \label{cstr_rejection} \\	281
\mathcal{N}_q \geqslant 160 & \text{ an user defined bound}\nonumber\\	282
& \text{ (e.g. 160~dB here)}\nonumber\\\nonumber	283
\end{align}	284
\paragraph{Goal}	285
\begin{align*}	286
\min \mathcal{S}_q	287
\end{align*}	288
	289	199
The constraint \ref{cstr_size} means the occupation for the current stage $j$ depends on	290	200	FIR $i$ can reject $F(C_i, \pi_i^C)$ dB. $F$ is determined numerically.
the previous occupation and the occupation of current selected filter (it is possible	291	201	To measure this rejection, we use GNU Octave software to design FIR filter coefficients thanks to two
that no filter is selected for this stage). And the second one \ref{cstr_rejection}	292	202	algorithms (\texttt{firls} and \texttt{fir1}).
means the same thing but for the rejection, the rejection depends the previous rejection	293	203	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.
plus the rejection of selected filter.	294	204	Then, the floating point coefficients are discretized into integers. In order to ensure that the coefficients are coded on $\pi_i^C$~bits effectively,
		205	the coefficients are normalized by their absolute maximum before being scaled to integer coefficients.
		206	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.
	295	207
\subsection{Low bandpass ripple and maximum rejection criteria}	296	208	With these coefficients, the \texttt{freqz} function is used to estimate the magnitude of the filter.
		209	Comparing the performance between FIRs requires however a unique criterion. As shown in figure~\ref{fig:fir_mag},
		210	the FIR magnitude exhibits two parts.
	297	211
The MILP solver provides a solution to the problem by selecting a series of small FIR with	298	212	\begin{figure}
increasing number of bits representing data and coefficients as well as an increasing number	299	213	\centering
of coefficients, instead of a single monolithic filter.	300	214	\begin{tikzpicture}[scale=0.3]
		215	\draw[<->] (0,15) -- (0,0) -- (21,0) ;
		216	\draw[thick] (0,12) -- (8,12) -- (20,0) ;
	301	217
\begin{figure}[h!tb]	302	218	\draw (0,14) node [left] { $P$ } ;
% \includegraphics[width=\linewidth]{images/compare-fir.pdf}	303	219	\draw (20,0) node [below] { $f$ } ;
\includegraphics[width=\linewidth]{images/fir-mono-vs-fir-series-noise-fixe-jmf-light.pdf}	304
\caption{Comparison of the rejection capability between a series of FIR and a monolithic FIR	305
with a cutoff frequency set at half the Nyquist frequency.}	306
\label{compare-fir}	307
\end{figure}	308
	309	220
Fig. \ref{compare-fir} exhibits the	310	221	\draw[>=latex,<->] (0,14) -- (8,14) ;
performance comparison between one solution and a monolithic FIR when selecting a cutoff	311	222	\draw (4,14) node [above] { passband } node [below] { $40\%$ } ;
frequency of half the Nyquist frequency: a series of 5 FIR and a series of 10 FIR with the	312
same space usage are provided as selected by the MILP solver. The FIR cascade provides improved	313
rejection than the monolithic FIR at the expense of a lower cutoff frequency which remains to	314
be tuned or compensated for.	315
	316	223
		224	\draw[>=latex,<->] (8,14) -- (12,14) ;
		225	\draw (10,14) node [above] { transition } node [below] { $20\%$ } ;
	317	226
The resource occupation when synthesizing such FIR on a Xilinx FPGA is summarized as Tab. \ref{t1}.	318	227	\draw[>=latex,<->] (12,14) -- (20,14) ;
We have considered a set of resources representative of the hardware platform we work on,	319	228	\draw (16,14) node [above] { stopband } node [below] { $40\%$ } ;
Avnet's Zedboard featuring a Xilinx XC7Z020-CLG484-1 Zynq System on Chip (SoC). The results reported in	320
Tab. \ref{t1} emphasize that implementing the monolithic single FIR is impossible due to	321
the insufficient hardware resources (exhausted LUT resources), while the FIR cascading 5 or 10	322
filters fit in the available resources. However, in all cases the DSP resources are fully	323
used: while the design can be synthesized using Xilinx proprietary Vivado 2016.2 software,	324
implementing the design fails due to the excessive resource usage preventing routing the signals	325
on the FPGA. Such results emphasize on the one hand the improvement prospect of the optimization	326
procedure by finding non-trivial solutions matching resource constraints, but on the other	327
hand also illustrates the limitation of a model with an abstraction layer that does not account	328
for the detailed architecture of the hardware.	329
	330	229
\begin{table}[h!tb]	331	230	\draw[>=latex,<->] (16,12) -- (16,8) ;
\caption{Resource occupation on a Xilinx Zynq-7000 series FPGA when synthesizing the FIR cascade	332	231	\draw (16,10) node [right] { rejection } ;
identified as optimal by the MILP solver within a finite resource criterion. The last line refers	333
to available resources on a Zynq-7020 as found on the Zedboard.}	334
\begin{center}	335
\begin{tabular}{\|c\|cccc\|}\hline	336
FIR & BlockRAM & LookUpTables & DSP & rejection (dB)\\\hline\hline	337
1 (monolithic) & 1 & 76183 & 220 & -162 \\	338
5 & 5 & 18597 & 220 & -160 \\	339
10 & 8 & 24729 & 220 & -161 \\\hline\hline	340
\textbf{Zynq 7020} & \textbf{420} & \textbf{53200} & \textbf{220} & \\\hline	341
%\begin{tabular}{\|c\|ccccc\|}\hline	342
%FIR & BRAM36 & BRAM18 & LUT & DSP & rejection (dB)\\\hline\hline	343
%1 (monolithic) & 1 & 0 & {\color{Red}76183} & 220 & -162 \\	344
%5 & 0 & 5 & {\color{Green}18597} & 220 & -160 \\	345
%10 & 0 & 8 & {\color{Green}24729} & 220 & -161 \\\hline\hline	346
%\textbf{Zynq 7020} & \textbf{140} & \textbf{280} & \textbf{53200} & \textbf{220} & \\\hline	347
\end{tabular}	348
\end{center}	349
%\vspace{-0.7cm}	350
\label{t1}	351
\end{table}	352
	353	232
\subsection{Alternate criteria}\label{median}	354	233	\draw[dashed] (8,-1) -- (8,14) ;
		234	\draw[dashed] (12,-1) -- (12,14) ;
	355	235
Fig. \ref{compare-fir} provides FIR solutions matching well the targeted transfer	356	236	\draw[dashed] (8,12) -- (16,12) ;
function, namely low ripple in the bandpass defined as the first 40\% of the frequency	357	237	\draw[dashed] (12,8) -- (16,8) ;
range and maximum rejection of 160~dB in the last 40\% stopband. We illustrate now, for	358
demonstrating the need to properly select the optimization criterion, two cases of poor	359
filter shapes obtained by selecting the mean value and median value of the rejection,	360
with no consideration for the ripples in the bandpass. The results of the optimizations,	361
in these cases, are shown in Figs. \ref{compare-mean} and \ref{compare-median}.	362
	363	238
\begin{figure}[h!tb]	364	239	\end{tikzpicture}
\includegraphics[width=\linewidth]{images/fir-mono-vs-fir-series-noise-fixe-mean-light.pdf}	365	240
\caption{Comparison of the rejection capability between a series of FIR and a monolithic FIR	366	241	% \includegraphics[width=.5\linewidth]{images/fir_magnitude}
with a cutoff frequency set at half the Nyquist frequency.}	367	242	\caption{Shape of the filter transmitted power $P$ as a function of frequency $f$:
\label{compare-mean}	368	243	the passband is considered to occupy the initial 40\% of the Nyquist frequency range,
		244	the stopband the last 40\%, allowing 20\% transition width.}
		245	\label{fig:fir_mag}
\end{figure}	369	246	\end{figure}
	370	247
In the case of the mean value criterion (Fig. \ref{compare-mean}), the solution is not	371	248	In the transition band, the behavior of the filter is left free, we only care about the passband and the stopband.
acceptable since the notch at the end of the transition band compensates for some unacceptable	372	249	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.
rise in the rejection close to the Nyquist frequency. Applying such a filter might yield excessive	373	250	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}.
high frequency spurious components to be aliased at low frequency when decimating the signal.	374
Similarly, the lack of criterion on the bandpass shape induces a shape with poor flatness and	375
and slowly decaying transfer function starting to attenuate spectral components well before the	376
transition band starts. Such issues are partly aleviated by replacing a mean rejection value with	377
a median rejection value (Fig. \ref{compare-median}) but solutions remain unacceptable for	378
the reasons stated previously and much poorer than those found with the maximum rejection criterion	379
selected earlier (Fig. \ref{compare-fir}).	380
	381	251
\begin{figure}[h!tb]	382	252	\begin{figure}
\includegraphics[width=\linewidth]{images/fir-mono-vs-fir-series-noise-fixe-median-light.pdf}	383	253	\centering
\caption{Comparison of the rejection capability between a series of FIR and a monolithic FIR	384	254	\includegraphics[width=\linewidth]{images/mean_criterion}
with a cutoff frequency set at half the Nyquist frequency.}	385	255	\caption{Mean criterion comparison between monolithic filter and cascade filters}
\label{compare-median}	386	256	\label{fig:mean_criterion}
\end{figure}	387	257	\end{figure}
	388	258
\section{Filter coefficient selection}	389	259	\begin{figure}
		260	\centering
		261	\includegraphics[width=\linewidth]{images/custom_criterion}
		262	\caption{Custom criterion comparison between monolithic filter and cascade filters}
		263	\label{fig:custom_criterion}
		264	\end{figure}
	390	265
The coefficients of a single monolithic filter are computed as the impulse response	391	266	Although we have a efficient criterion to estimate the rejection of one set of coefficient
of the filter transfer function, and practically approximated by a multitude of methods	392	267	we have a problem when we sum two or more criterion. If the FIR filter coefficients are the same
including least square optimization (Matlab's {\tt firls} function), Hamming or Kaiser windowing	393	268	between the stage, we have:
(Matlab's {\tt fir1} function).	394	269	$$F_{total} = F_1 + F_2$$
		270	But when we choose two different set of coefficient, the previous equality are not
		271	true. The figure~\ref{fig:sum_rejection} illustrates the problem. The red and blue curves
		272	are two different filter coefficient and we can see that their maximum on the stopband
		273	are not at the same frequency. So when we sum the rejection criteria (the dotted yellow line)
		274	we do not meet the dashed yellow line. Define the rejection of cascaded filters
		275	is more difficult than just take the summation between all the rejection criteria of each filter.
		276	However this summation gives us an upper bound for rejection although in fact we obtain
		277	better rejection than expected.
	395	278
\begin{figure}[h!tb]	396	279	\begin{figure}
\includegraphics[width=\linewidth]{images/fir1-vs-firls}	397	280	\centering
\caption{Evolution of the rejection capability of least-square optimized filters and Hamming	398	281	\includegraphics[width=\linewidth]{images/sum_rejection}
FIR filters as a function of the number of coefficients, for floating point numbers and 8-bit	399	282	\caption{Rejection of two cascaded filters}
encoded integers.}	400	283	\label{fig:sum_rejection}
\label{2}	401
\end{figure}	402	284	\end{figure}
	403	285
Cascading filters opens a new optimization opportunity by	404	286	\section{Experiments with fixed area space}
selecting various coefficient sets depending on the number of coefficients. Fig. \ref{2}	405
illustrates that for a number of coefficients ranging from 8 to 47, {\tt fir1} provides a better	406
rejection than {\tt firls}: since the linear solver increases the number of coefficients along	407
the processing chain, the type of selected filter also changes depending on the number of coefficients	408
and evolves along the processing chain.	409
	410	287
\section{Conclusion}	411	288	\begin{figure}
		289	\centering
		290	\includegraphics[width=\linewidth]{images/max_rejection/prn_500}
		291	\caption{Experimental results for design with PRN as data input and 500 a.u. as max arbitrary space}
		292	\label{fig:prn_500}
		293	\end{figure}
	412	294
We address the optimization problem of designing a low-pass filter chain in a Field Programmable Gate	413	295	\begin{figure}
Array for improved noise rejection within constrained resource occupation, as needed for	414	296	\centering
real time processing of radiofrequency signal when characterizing spectral phase noise	415	297	\includegraphics[width=\linewidth]{images/max_rejection/prn_1000}
characteristics of stable oscillators. The flexibility of the digital approach makes the result	416	298	\caption{Experimental results for design with PRN as data input and 1000 a.u. as max arbitrary space}
best suited for closing the loop and using the measurement output in a feedback loop for	417	299	\label{fig:prn_1000}
controlling clocks, e.g. in a quartz-stabilized high performance clock whose long term behavior	418	300	\end{figure}
is controlled by non-piezoelectric resonator (sapphire resonator, microwave or optical	419
atomic transition).	420
	421	301
\section*{Acknowledgement}	422	302	\begin{figure}
		303	\centering
		304	\includegraphics[width=\linewidth]{images/max_rejection/prn_2000}
		305	\caption{Experimental results for design with PRN as data input and 2000 a.u. as max arbitrary space}
		306	\label{fig:prn_2000}
		307	\end{figure}
	423	308
This work is supported by the ANR Programme d'Investissement d'Avenir in	424	309	\begin{table}
progress at the Time and Frequency Departments of the FEMTO-ST Institute	425	310	\centering
(Oscillator IMP, First-TF and Refimeve+), and by R\'egion de Franche-Comt\'e.	426	311	\begin{tabular}{\|c\|c\|ccc\|c\|c\|}
The authors would like to thank E. Rubiola, F. Vernotte, and G. Cabodevila	427	312	\hline
for support and fruitful discussions.	428	313	\multicolumn{2}{\|c\|}{\multirow{2}{}{Stage}} & \multicolumn{3}{c\|}{Stage} & \multirow{2}{}{Rejection} & \multirow{2}{*}{Area} \\ \cline{3-5}
		314	\multicolumn{2}{\|c\|}{} & i = 1 & i = 2 & i = 3 & & \\ \hline
		315	& C & 19 & - & - & & \\
		316	n = 1 & $pi^C$ & 7 & - & - & 33 dB & 437 a.u. \\
		317	& $pi^S$ & 0 & - & - & & \\ \hline
		318	& C & 11 & 19 & - & & \\
		319	n = 2 & $pi^C$ & 5 & 7 & - & 53 dB & 478 a.u. \\
		320	& $pi^S$ & 16 & 0 & - & & \\ \hline
		321	& C & 9 & 15 & 11 & & \\
		322	n = 3 & $pi^C$ & 4 & 6 & 5 & 57 dB & 499 a.u. \\
		323	& $pi^S$ & 16 & 3 & 0 & & \\ \hline
		324	\end{tabular}
		325	\caption{Solver results for design with PRN as data input and 500 a.u. as max arbitrary space}
		326	\label{tbl:prn_500}
		327	\end{table}
	429	328
\bibliographystyle{IEEEtran}	430	329	\begin{table}
\balance	431	330	\centering
\bibliography{references,biblio}	432	331	{\scalefont{0.85}
\end{document}	433	332	\begin{tabular}{\|c\|c\|ccccc\|c\|c\|}
		333	\hline
		334	\multicolumn{2}{\|c\|}{\multirow{2}{}{Stage}} & \multicolumn{5}{c\|}{Stage} & \multirow{2}{}{Rejection} & \multirow{2}{*}{Area} \\ \cline{3-7}
		335	\multicolumn{2}{\|c\|}{} & i = 1 & i = 2 & i = 3 & i = 4 & i = 5 & & \\ \hline
		336	& C & 37 & - & - & - & - & & \\
		337	n = 1 & $pi^C$ & 11 & - & - & - & - & 56 dB & 999 a.u. \\
		338	& $pi^S$ & 0 & - & - & - & - & & \\ \hline
		339	& C & 11 & 39 & - & - & - & & \\
		340	n = 2 & $pi^C$ & 5 & 13 & - & - & - & 82 dB & 972 a.u. \\
		341	& $pi^S$ & 16 & 0 & - & - & - & & \\ \hline
		342	& C & 9 & 31 & 19 & - & - & & \\
		343	n = 3 & $pi^C$ & 7 & 8 & 7 & - & - & 93 dB & 990 a.u. \\
		344	& $pi^S$ & 19 & 2 & 0 & - & - & & \\ \hline
		345	& C & 9 & 19 & 17 & 11 & - & & \\
		346	n = 4 & $pi^C$ & 4 & 7 & 7 & 5 & - & 99 dB & 992 a.u. \\
		347	& $pi^S$ & 16 & 3 & 3 & 0 & - & & \\ \hline
		348	& C & 9 & 15 & 11 & 11 & 11 & & \\
		349	n = 5 & $pi^C$ & 4 & 7 & 5 & 5 & 5 & 99 dB & 998 a.u. \\
		350	& $pi^S$ & 16 & 3 & 2 & 1 & 1 & & \\ \hline
		351	\end{tabular}
		352	}
		353	\caption{Solver results for design with PRN as data input and 1000 a.u. as max arbitrary space}
		354	\label{tbl:prn_1000}
		355	\end{table}
	434	356
\section{Contexte d'ordonnancement}	435	357	\begin{table}
Dans cette partie, nous donnerons des d\'efinitions de termes rattach\'es au domaine de l'ordonnancement	436	358	\centering
et nous verrons que le sujet trait\'e se rapproche beaucoup d'un problème d'ordonnancement. De ce fait	437	359	{\scalefont{0.85}
nous pourrons aller plus loin que les travaux vus pr\'ec\'edemment et nous tenterons des approches d'ordonnancement	438	360	\begin{tabular}{\|c\|c\|ccccc\|c\|c\|}
et d'optimisation.	439	361	\hline
		362	\multicolumn{2}{\|c\|}{\multirow{2}{}{Stage}} & \multicolumn{5}{c\|}{Stage} & \multirow{2}{}{Rejection} & \multirow{2}{*}{Area} \\ \cline{3-7}
		363	\multicolumn{2}{\|c\|}{} & i = 1 & i = 2 & i = 3 & i = 4 & i = 5 & & \\ \hline
		364	& C & 39 & - & - & - & - & & \\
		365	n = 1 & $pi^C$ & 13 & - & - & - & - & 61 dB & 1131 a.u. \\
		366	& $pi^S$ & 0 & - & - & - & - & & \\ \hline
		367	& C & 37 & 39 & - & - & - & & \\
		368	n = 2 & $pi^C$ & 11 & 13 & - & - & - & 117 dB & 1974 a.u. \\
		369	& $pi^S$ & 17 & 0 & - & - & - & & \\ \hline
		370	& C & 15 & 35 & 35 & - & - & & \\
		371	n = 3 & $pi^C$ & 9 & 11 & 11 & - & - & 138 dB & 1985 a.u. \\
		372	& $pi^S$ & 19 & 3 & 0 & - & - & & \\ \hline
		373	& C & 11 & 27 & 27 & 23 & - & & \\
		374	n = 4 & $pi^C$ & 5 & 9 & 9 & 9 & - & 148 dB & 1993 a.u. \\
		375	& $pi^S$ & 16 & 3 & 2 & 0 & - & & \\ \hline
		376	& C & 11 & 27 & 31 & 11 & 11 & & \\
		377	n = 5 & $pi^C$ & 5 & 9 & 8 & 5 & 5 & 153 dB & 2000 a.u. \\
		378	& $pi^S$ & 16 & 3 & 1 & 0 & 1 & & \\ \hline
		379	\end{tabular}
		380	}
		381	\caption{Solver results for design with PRN as data input and 2000 a.u. as max arbitrary space}
		382	\label{tbl:prn_2000}
		383	\end{table}
	440	384
\subsection{D\'efinition du vocabulaire}	441	385	\begin{table}
Avant tout, il faut d\'efinir ce qu'est un problème d'optimisation. Il y a deux d\'efinitions	442	386	\centering
importantes à donner. La première est propos\'ee par Legrand et Robert dans leur livre \cite{def1-ordo} :	443	387	\begin{tabular}{\|c\|c\|c\|c\|c\|}\hline
\begin{definition}	444	388	Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
\label{def-ordo1}	445	389	& 1 & 0.02~s & $\approx$ 20 min & $\approx$ 1 min \\
Un ordonnancement d'un système de t\^aches $G\ =\ (V,\ E,\ w)$ est une fonction $\sigma$ :	446	390	PRN & 2 & 1.70~s & $\approx$ 20 min & $\approx$ 1 min \\
$V \rightarrow \mathbb{N}$ telle que $\sigma(u) + w(u) \leq \sigma(v)$ pour toute arête $(u,\ v) \in E$.	447	391	& 3 & 19~s & $\approx$ 20 min & $\approx$ 1 min \\\hline
\end{definition}	448	392	\end{tabular}
		393	\caption{Time to compute and deploy the designs for PRN 500}
		394	\label{tbl:time_prn_500}
		395	\end{table}
	449	396
Dit plus simplement, l'ensemble $V$ repr\'esente les t\^aches à ex\'ecuter, l'ensemble $E$ repr\'esente les d\'ependances	450	397	\begin{table}
des t\^aches et $w$ les temps d'ex\'ecution de la t\^ache. La fonction $\sigma$ donne donc l'heure de d\'ebut de	451	398	\centering
chacune des t\^aches. La d\'efinition dit que si une t\^ache $v$ d\'epend d'une t\^ache $u$ alors	452	399	\begin{tabular}{\|c\|c\|c\|c\|c\|}\hline
la date de d\'ebut de $v$ sera plus grande ou \'egale au d\'ebut de l'ex\'ecution de la t\^ache $u$ plus son	453	400	Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
temps d'ex\'ecution.	454	401	& 1 & 0.07~s & $\approx$ 20 min & $\approx$ 1 min \\
		402	& 2 & 1.31~s & $\approx$ 20 min & $\approx$ 1 min \\
		403	PRN & 3 & 119~s ($\approx$ 2~min) & $\approx$ 20 min & $\approx$ 1 min \\
		404	& 4 & 270~s ($\approx$ 5~min) & $\approx$ 20 min & $\approx$ 1 min \\
		405	& 5 & 5998~s ($\approx$ 2~h) & $\approx$ 20 min & $\approx$ 1 min \\\hline
		406	\end{tabular}
		407	\caption{Time to compute and deploy the designs for PRN 1000}
		408	\label{tbl:time_prn_1000}
		409	\end{table}
	455	410
Une autre d\'efinition importante qui est propos\'ee par Leung et al. \cite{def2-ordo} est :	456	411	\begin{table}
\begin{definition}	457	412	\centering
\label{def-ordo2}	458	413	\begin{tabular}{\|c\|c\|c\|c\|c\|}\hline
L'ordonnancement traite de l'allocation de ressources rares à des activit\'es avec	459	414	Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
l'objectif d'optimiser un ou plusieurs critères de performance.	460	415	& 1 & 0.07~s & $\approx$ 20 min & $\approx$ 1 min \\
\end{definition}	461	416	& 2 & 0.75~s & $\approx$ 20 min & $\approx$ 1 min \\
		417	PRN & 3 & 36~s & - & - \\
		418	& 4 & 14500~s ($\approx$ 4~h) & $\approx$ 20 min & $\approx$ 1 min \\
		419	& 5 & 74237~s ($\approx$ 20~h) & $\approx$ 20 min & $\approx$ 1 min \\\hline
		420	\end{tabular}
		421	\caption{Time to compute and deploy the designs for PRN 2000}
		422	\label{tbl:time_prn_2000}
		423	\end{table}
	462	424
Cette d\'efinition est plus g\'en\'erique mais elle nous int\'eresse d'avantage que la d\'efinition \ref{def-ordo1}.	463	425	\section{Experiments with fixed rejection target}
En effet, la partie qui nous int\'eresse dans cette première d\'efinition est le respect de la pr\'ec\'edance des t\^aches.	464
Dans les faits les dates de d\'ebut ne nous int\'eressent pas r\'eellement.	465
	466	426
En revanche la d\'efinition \ref{def-ordo2} sera au c\oe{}ur du projet. Pour se convaincre de cela,	467	427	\begin{figure}
il nous faut d'abord d\'efinir quel est le type de problème d'ordonnancement qu'on traite et quelles	468	428	\centering
sont les m\'ethodes qu'on peut appliquer.	469	429	\includegraphics[width=\linewidth]{images/min_area/prn_50}
		430	\caption{Results for design with PRN as data input and 50 dB as aimed rejection level}
		431	\label{fig:prn_500}
		432	\end{figure}
	470	433
Les problèmes d'ordonnancement peuvent être class\'es en diff\'erentes cat\'egories :	471	434	\begin{figure}
\begin{itemize}	472	435	\centering
\item T\^aches ind\'ependantes : dans cette cat\'egorie de problèmes, les t\^aches sont complètement ind\'ependantes	473	436	\includegraphics[width=\linewidth]{images/min_area/prn_100}
les unes des autres. Dans notre cas, ce n'est pas le plus adapt\'e.	474	437	\caption{Results for design with PRN as data input and 50 dB as aimed rejection level}
\item Graphe de t\^aches : la d\'efinition \ref{def-ordo1} d\'ecrit cette cat\'egorie. La plupart du temps,	475	438	\label{fig:prn_100}
les t\^aches sont repr\'esent\'ees par une DAG. Cette cat\'egorie est très proche de notre cas puisque nous devons \'egalement ex\'ecuter	476	439	\end{figure}
des t\^aches qui ont un certain nombre de d\'ependances. On pourra même dire que dans certain cas,	477
on a des anti-arbres, c'est à dire que nous avons une multitude de t\^aches d'entr\'ees qui convergent vers une	478
t\^ache de fin.	479
\item Workflow : cette cat\'egorie est une sous cat\'egorie des graphes de t\^aches dans le sens où	480
il s'agit d'un graphe de t\^aches r\'ep\'et\'e de nombreuses de fois. C'est exactement ce type de problème	481
que nous traitons ici.	482
\end{itemize}	483
	484	440
Bien entendu, cette liste n'est pas exhaustive et il existe de nombreuses autres classifications et sous-classifications	485	441	\begin{figure}
de ces problèmes. Nous n'avons parl\'e ici que des cat\'egories les plus communes.	486	442	\centering
		443	\includegraphics[width=\linewidth]{images/min_area/prn_150}
		444	\caption{Results for design with PRN as data input and 2000 a.u. as max arbitrary space}
		445	\label{fig:prn_150}
		446	\end{figure}
	487	447
Un autre point à d\'efinir, est le critère d'optimisation. Il y a là encore un grand nombre de	488	448	\section{Conclusion}
critères possibles. Nous allons donc parler des principaux :	489
\begin{itemize}	490
\item Temps de compl\'etion total (ou Makespan en anglais) : ce critère est l'un des critères d'optimisation	491
les plus courant. Il s'agit donc de minimiser la date de fin de la dernière t\^ache de l'ensemble des	492
t\^aches à ex\'ecuter. L'enjeu de cette optimisation est donc de trouver l'ordonnancement optimal permettant	493
la fin d'ex\'ecution au plus tôt.	494
\item Somme des temps d'ex\'ecution (Flowtime en anglais) : il s'agit de faire la somme des temps d'ex\'ecution de toutes les t\^aches	495
et d'optimiser ce r\'esultat.	496
\item Le d\'ebit : ce critère quant à lui, vise à augmenter au maximum le d\'ebit de traitement des donn\'ees.	497
\end{itemize}	498
	499	449
En plus de cela, on peut avoir besoin de plusieurs critères d'optimisation. Il s'agit dans ce cas d'une optimisation	500	450	\section*{Acknowledgement}
multi-critères. Bien entendu, cela complexifie d'autant plus le problème car la solution la plus optimale pour un	501
des critères peut être très mauvaise pour un autre critère. De ce cas, il s'agira de trouver une solution qui permet	502
de faire le meilleur compromis entre tous les critères.	503
	504	451
\subsection{Formalisation du problème}	505	452	This work is supported by the ANR Programme d'Investissement d'Avenir in
\label{formalisation}	506	453	progress at the Time and Frequency Departments of the FEMTO-ST Institute
Maintenant que nous avons donn\'e le vocabulaire li\'e à l'ordonnancement, nous allons pouvoir essayer caract\'eriser	507	454	(Oscillator IMP, First-TF and Refimeve+), and by R\'egion de Franche-Comt\'e.
formellement notre problème. En effet, nous allons reprendre les contraintes \'enonc\'ees dans la sections \ref{def-contraintes}	508	455	The authors would like to thank E. Rubiola, F. Vernotte, and G. Cabodevila
et nous essayerons de les formaliser le plus finement possible.	509	456	for support and fruitful discussions.
	510	457
Comme nous l'avons dit, une t\^ache est un bloc de traitement. Chaque t\^ache $i$ dispose d'un ensemble de paramètres	511	458	\bibliographystyle{IEEEtran}
que nous nommerons $\mathcal{P}_{i}$. Cet ensemble $\mathcal{P}_i$ est propre à chaque t\^ache et il variera d'une	512	459	\balance
t\^ache à l'autre. Nous reviendrons plus tard sur les paramètres qui peuvent composer cet ensemble.	513	460	\bibliography{references,biblio}
	514	461	\end{document}
Outre cet ensemble $\mathcal{P}_i$, chaque t\^ache dispose de paramètres communs :	515
\begin{itemize}	516
\item Dur\'ee de la t\^ache : Comme nous l'avons dit auparavant, dans le cadre d'un FPGA le temps est compt\'e en nombre de coup d'horloge.	517
En outre, les blocs sont toujours sollicit\'es, certains même sont capables de lire et de renvoyer une r\'esultat à chaque coups d'horloge.	518
Donc la dur\'ee d'une t\^ache ne peut être le laps de temps entre l'entr\'ee d'une donn\'ee et la sortie d'une autre. Nous d\'efinirons la	519
dur\'ee comme le temps de traitement d'une donn\'ee, c'est à dire la diff\'erence de temps entre la date de sortie d'une donn\'ee	520
et de sa date d'entr\'ee. Nous nommerons cette dur\'ee $\delta_i$. % Je devrais la nomm\'ee w comme dans la def2	521
\item La pr\'ecision : La pr\'ecision d'une donn\'ee est le nombre de bits significatifs qu'elle compte. En effet, au fil des traitements	522
les pr\'ecisions peuvent varier. On nomme donc la pr\'ecision d'entr\'ee d'une t\^ache $i$ comme $\pi_i^-$ et la pr\'ecision en sortie $\pi_i^+$.	523
\item La fr\'equence du flux en entr\'ee (ou sortie) : Cette fr\'equence repr\'esente la fr\'equence des donn\'ees qui arrivent (resp. sortent).	524
Selon les t\^aches, les fr\'equences varieront. En effet, certains blocs ralentissent le flux c'est pourquoi on distingue la fr\'equence du	525
flux en entr\'ee et la fr\'equence en sortie. Nous nommerons donc la fr\'equence du flux en entr\'ee $f_i^-$ et la fr\'equence en sortie $f_i^+$.	526
\item La quantit\'e de donn\'ees en entr\'ee (ou en sortie) : Il s'agit de la quantit\'e de donn\'ees que le bloc s'attend à traiter (resp.	527
est capable de produire). Les t\^aches peuvent avoir à traiter des gros volumes de donn\'ees et n'en ressortir qu'une partie. Cette	528
fois encore, il nous faut donc diff\'erencier l'entr\'ee et la sortie. Nous nommerons donc la quantit\'e de donn\'ees entrantes $q_i^-$	529
et la quantit\'e de donn\'ees sortantes $q_i^+$ pour une t\^ache $i$.	530
\item Le d\'ebit d'entr\'ee (ou de sortie) : Ce paramètre correspond au d\'ebit de donn\'ees que la t\^ache est capable de traiter ou qu'elle	531
fournit en sortie. Il s'agit simplement de l'expression des deux pr\'ec\'edents paramètres. Nous d\'efinirons donc la d\'ebit entrant de la	532
t\^ache $i$ comme $d_i^-\ =\ q_i^-\ \ f_i^-$ et le d\'ebit sortant comme $d_i^+\ =\ q_i^+\ \ f_i^+$.	533
\item La taille de la t\^ache : La taille dans les FPGA \'etant limit\'ee, ce paramètre exprime donc la place qu'occupe la t\^ache au sein du bloc.	534
Nous nommerons $\mathcal{A}_i$ cette taille.	535
\item Les pr\'ed\'ecesseurs et successeurs d'une t\^ache : cela nous permet de connaître les t\^aches requises pour pouvoir traiter	536
la t\^ache $i$ ainsi que les t\^aches qui en d\'ependent. Ces ensemble sont not\'es $\Gamma _i ^-$ et $ \Gamma _i ^+$ \\	537
%TODO Est-ce vraiment un paramètre ?	538
\end{itemize}	539
	540
Ces diff\'erents paramètres communs sont fortement li\'es aux \'el\'ements de $\mathcal{P}_i$. Voici quelques exemples de relations	541
que nous avons identifi\'ees :	542
\begin{itemize}	543
\item $ \delta _i ^+ \ = \ \mathcal{F}_{\delta}(\pi_i^-,\ \pi_i^+,\ d_i^-,\ d_i^+,\ \mathcal{P}_i) $ donne le temps d'ex\'ecution	544
de la t\^ache en fonction de la pr\'ecision voulue, du d\'ebit et des paramètres internes.	545
\item $ \pi _i ^+ \ = \ \mathcal{F}_{p}(\pi_i^-,\ \mathcal{P}_i) $, la fonction $F_p$ donne la pr\'ecision en sortie selon la pr\'ecision de d\'epart	546
et les paramètres internes de la t\^ache.	547
\item $d_i^+\ =\ \mathcal{F}_d(d_i^-, \mathcal{P}_i)$, la fonction $F_d$ donne le d\'ebit sortant de la t\^ache en fonction du d\'ebit	548
sortant et des variables internes de la t\^ache.	549
\item $A_i^+\ =\ \mathcal{F}_A(\pi_i^-,\ \pi_i^+,\ d_i^-,\ d_i^+, \mathcal{P}_i)$	550
\end{itemize}	551
Pour le moment, nous ne sommes pas capables de donner une d\'efinition g\'en\'erale de ces fonctions. Mais en revanche,	552
sur quelques exemples simples (cf. \ref{def-contraintes}), nous parvenons à donner une \'evaluation de ces fonctions.	553
	554
Maintenant que nous avons donn\'e toutes les notations utiles, nous allons \'enoncer des contraintes relatives à notre problème. Soit	555
un DGA $G(V,\ E)$, on a pour toutes arêtes $(i, j)\ \in\ E$ les in\'equations suivantes :	556
	557
\paragraph{Contrainte de pr\'ecision :}	558
Cette in\'equation traduit la contrainte de pr\'ecision d'une t\^ache à l'autre :	559
\begin{align*}	560
\pi _i ^+ \geq \pi _j ^-	561
\end{align*}	562
	563
\paragraph{Contrainte de d\'ebit :}	564
Cette in\'equation traduit la contrainte de d\'ebit d'une t\^ache à l'autre :	565
\begin{align*}	566
d _i ^+ = q _j ^- * (f_i + (1 / s_j) ) & \text{ où } s_j \text{ est une valeur positive de temporisation de la t\^ache}	567
\end{align*}	568
	569
\paragraph{Contrainte de synchronisation :}	570
Il s'agit de la contrainte qui impose que si à un moment du traitement, le DAG se s\'epare en plusieurs branches parallèles	571
et qu'elles se rejoignent plus tard, la somme des latences sur chacune des branches soit la même.	572
Plus formellement, s'il existe plusieurs chemins disjoints, partant de la t\^ache $s$ et allant à la t\^ache de $f$ alors :	573
\begin{align*}	574
\forall \text{ chemin } \mathcal{C}1(s, .., f),	575
\forall \text{ chemin } \mathcal{C}2(s, .., f)	576
\text{ tel que } \mathcal{C}1 \neq \mathcal{C}2	577
\Rightarrow	578
\sum _{i} ^{i \in \mathcal{C}1} \delta_i = \sum _{i} ^{i \in \mathcal{C}2} \delta_i	579
\end{align*}	580
	581
\paragraph{Contrainte de place :}	582
Cette in\'equation traduit la contrainte de place dans le FPGA. La taille max de la puce FPGA est nomm\'e $\mathcal{A}_{FPGA}$ :	583
\begin{align*}	584
\sum ^{\text{t\^ache } i} \mathcal{A}_i \leq \mathcal{A}_{FPGA}	585
\end{align*}	586
	587
\subsection{Exemples de mod\'elisation}	588
\label{exemples-modeles}	589
Nous allons maintenant prendre quelques blocs de traitement simples afin d'illustrer au mieux notre modèle.	590
Pour tous nos exemple, nous prendrons un d\'ebit en entr\'ee de 200 Mo/s avec une pr\'ecision de 16 bit.	591
	592
Prenons tout d'abord l'exemple d'un bloc de d\'ecimation. Le but de ce bloc est de ralentir le flux en ne gardant	593
que certaines donn\'ees à intervalle r\'egulier. Cet intervalle est appel\'e le facteur de d\'ecimation, on le notera $N$.	594
	595
Donc d'après notre mod\'elisation :	596
\begin{itemize}	597
\item $N \in \mathcal{P}_i$	598
%TODO N ou 1 ?	599
\item $\delta _i = N\ c.h.$ (coup d'horloge)	600
\item $\pi _i ^+ = \pi _i ^- = 16 bits$	601
\item $f _i ^+ = f _i ^-$	602
\item $q _i ^+ = q _i ^- / N$	603
\item $d _i ^+ = q _i ^- / N / f _i ^-$	604
\item $\Gamma _i ^+ = \Gamma _i ^- = 1$\\	605
%TODO Je ne sais pas trouver la taille...	606
\end{itemize}	607
	608
Un autre exemple int\'eressant que l'on peut donner, c'est le cas des spliters. Il s'agit la aussi d'un bloc très	609
simple qui permet de dupliquer un flux. On peut donc donner un nombre de sorties à cr\'eer, on note ce paramètre	610
%TODO pas très inspir\'e...	611
$X$. Voici ce que donne notre mod\'elisation :	612
\begin{itemize}	613
\item $X \in \mathcal{P}_i$	614
\item $\delta _i = 1\ c.h.$	615
\item $\pi _i ^+ = \pi _i ^- = 16 bits$	616
\item $f _i ^+ = f _i ^-$	617
\item $q _i ^+ = q _i ^-$	618
\item $d _i ^+ = d _i ^-$	619
\item $\Gamma _i ^- = 1$	620
\item $\Gamma _i ^+ = X$\\	621
\end{itemize}	622
	623
L'exemple suivant traite du cas du shifter. Il s'agit d'un bloc qui a pour but de diminuer le nombre de bits des	624
donn\'ees afin d'acc\'el\'erer les traitement sur les blocs suivants. On peut donc donner le nombre de bits à shifter,	625
on note ce paramètre $S$. Voici ce que donne notre mod\'elisation :	626
\begin{itemize}	627
\item $S \in \mathcal{P}_i$	628
\item $\delta _i = 1\ c.h.$	629
\item $\pi _i ^+ = \pi _i ^- - S$	630
\item $f _i ^+ = f _i ^-$	631
\item $q _i ^+ = q _i ^-$	632
\item $d _i ^+ = d _i ^-$	633
\item $\Gamma _i ^+ = \Gamma _i ^- = 1$\\	634
\end{itemize}	635
	636
Nous allons traiter un dernier exemple un peu plus complexe, le cas d'un filtre d\'ecimateur (ou FIR). Ce bloc	637
est compos\'e de beaucoup de paramètres internes. On peut d\'efinir un nombre d'\'etages $E$, qui repr\'esente le nombre	638
d'it\'erations à faire avant d'arrêter le traitement. Afin d'effectuer son filtrage, on doit donner au bloc un ensemble	639
de coefficients $C$ et par cons\'equent ces coefficients ont leur propre pr\'ecision $\pi _C$. Pour finir, le dernier	640
paramètre à donner est le facteur de d\'ecimation $N$. Si on applique notre mod\'elisation, on peut obtenir cela :	641
\begin{itemize}	642
\item $E \in \mathcal{P}_i$	643
\item $C \in \mathcal{P}_i$	644
\item $\pi _C \in \mathcal{P}_i$	645
\item $N \in \mathcal{P}_i$	646
\item $\delta _i = E * \|C\| * q_i^-\ c.h.$ %Trop simpliste	647
\item $\pi _i ^+ = \pi _i ^- * \pi _C$	648
\item $f _i ^+ = f _i ^-$	649
\item $q _i ^+ = q _i ^- / N$	650
\item $d _i ^+ = q _i ^- / N / f _i ^-$	651
\item $\Gamma _i ^+ = \Gamma _i ^- = 1$\\	652
\end{itemize}	653
	654
Ces exemples ne sont que des modèles provisoires; pour s'assurer de leur performance, il faudra les	655
confronter à des simulations.	656
	657
	658
Bien que les articles sur les skeletons, \cite{gwen-cogen}, \cite{skeleton} et \cite{hide}, nous aient donn\'e des indices sur une possible	659
mod\'elisation, ils \'etaient encore trop focalis\'es sur l'optimisation spatiale des blocs. Nous nous sommes donc inspir\'es de ces travaux	660
pour proposer notre modèle, en faisant abstraction des optimisations bas niveau.	661
	662	462