Compare View

switch
from
...
to
 
Commits (2)

Diff

Showing 13 changed files Inline Diff

ifcs2018_abstract.aux 1 1 ifcs2018_abstract.aux
ifcs2018_abstract.bbl 2 2 ifcs2018_abstract.bbl
ifcs2018_abstract.blg 3 3 ifcs2018_abstract.blg
ifcs2018_abstract.log 4 4 ifcs2018_abstract.log
ifcs2018_abstract.out 5 5 ifcs2018_abstract.out
ifcs2018_abstract.pdf 6 6 ifcs2018_abstract.pdf
7 7
ifcs2018_proceeding.aux 8 8 ifcs2018_proceeding.aux
ifcs2018_proceeding.bbl 9 9 ifcs2018_proceeding.bbl
ifcs2018_proceeding.blg 10 10 ifcs2018_proceeding.blg
ifcs2018_proceeding.log 11 11 ifcs2018_proceeding.log
ifcs2018_proceeding.out 12 12 ifcs2018_proceeding.out
ifcs2018_proceeding.pdf 13 13 ifcs2018_proceeding.pdf
14 14
ifcs2018_poster.aux 15 15 ifcs2018_poster.aux
ifcs2018_poster.log 16 16 ifcs2018_poster.log
ifcs2018_poster.out 17 17 ifcs2018_poster.out
ifcs2018_poster.pdf 18 18 ifcs2018_poster.pdf
19 19
20 ifcs2018_journal.aux
21 ifcs2018_journal.bbl
22 ifcs2018_journal.blg
23 ifcs2018_journal.log
24 ifcs2018_journal.out
25 ifcs2018_journal.pdf
26
ifcs2018_journal.aux 20 27 *.bak
ifcs2018_journal.bbl 21 28
ifcs2018_journal.blg 22
ifcs2018_journal.log 23
# source: https://tex.stackexchange.com/questions/40738/how-to-properly-make-a-latex-project 1 1 # source: https://tex.stackexchange.com/questions/40738/how-to-properly-make-a-latex-project
2 2
TEX = pdflatex -shell-escape -interaction=nonstopmode -file-line-error 3 3 TEX = pdflatex -shell-escape -interaction=nonstopmode -file-line-error
BIB = bibtex 4 4 BIB = bibtex
TARGET = ifcs2018 5 5 TARGET = ifcs2018
6 6
all: $(TARGET)_abstract $(TARGET)_poster $(TARGET)_proceeding $(TARGET)_journal 7 7 all: $(TARGET)_abstract $(TARGET)_poster $(TARGET)_proceeding $(TARGET)_journal
8 8
view: $(TARGET) 9 9 view: $(TARGET)
evince $(TARGET).pdf 10 10 evince $(TARGET).pdf
11 11
view_poster: $(TARGET)_poster 12 12 view_poster: $(TARGET)_poster
evince $(TARGET)_poster.pdf 13 13 evince $(TARGET)_poster.pdf
14 14
$(TARGET)_abstract: $(TARGET)_abstract.tex references.bib biblio.bib 15 15 $(TARGET)_abstract: $(TARGET)_abstract.tex references.bib biblio.bib
$(TEX) $@.tex 16 16 $(TEX) $@.tex
$(BIB) $@ 17 17 $(BIB) $@
$(TEX) $@.tex 18 18 $(TEX) $@.tex
$(TEX) $@.tex 19 19 $(TEX) $@.tex
20 20
$(TARGET)_poster: 21 21 $(TARGET)_poster:
$(TEX) $@.tex 22 22 $(TEX) $@.tex
$(TEX) $@.tex 23 23 $(TEX) $@.tex
24 24
$(TARGET)_proceeding: $(TARGET)_proceeding.tex references.bib biblio.bib 25 25 $(TARGET)_proceeding: $(TARGET)_proceeding.tex references.bib biblio.bib
$(TEX) $@.tex 26 26 $(TEX) $@.tex
$(BIB) $@ 27 27 $(BIB) $@
$(TEX) $@.tex 28 28 $(TEX) $@.tex
$(TEX) $@.tex 29 29 $(TEX) $@.tex
30 30
31 $(TARGET)_journal: $(TARGET)_journal.tex references.bib biblio.bib
32 $(TEX) $@.tex
33 $(BIB) $@
34 $(TEX) $@.tex
35 $(TEX) $@.tex
36
$(TARGET)_journal: $(TARGET)_journal.tex references.bib biblio.bib 31 37 clean:
$(TEX) $@.tex 32 38 rm -f $(TARGET)_abstract.aux $(TARGET)_abstract.log $(TARGET)_abstract.out $(TARGET)_abstract.bbl $(TARGET)_abstract.blg
$(BIB) $@ 33
$(TEX) $@.tex 34 39 rm -f $(TARGET)_poster.aux $(TARGET)_poster.log $(TARGET)_poster.out
40 rm -f $(TARGET)_proceeding.aux $(TARGET)_proceeding.log $(TARGET)_proceeding.out $(TARGET)_proceeding.bbl $(TARGET)_proceeding.blg
41 rm -f $(TARGET)_journal.aux $(TARGET)_journal.log $(TARGET)_journal.out $(TARGET)_journal.bbl $(TARGET)_journal.blg
$(TEX) $@.tex 35 42
36 43 mrproper: clean
clean: 37 44 rm -f $(TARGET)_abstract.pdf $(TARGET)_proceeding.pdf $(TARGET)_poster.pdf $(TARGET)_journal.pdf
demo_critere_filtre.m
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
File was created 1 \documentclass[a4paper,conference]{IEEEtran/IEEEtran}
2 \usepackage{graphicx,color,hyperref}
3 \usepackage{amsfonts}
4 \usepackage{amsthm}
5 \usepackage{amssymb}
6 \usepackage{amsmath}
7 \usepackage{algorithm2e}
8 \usepackage{url,balance}
9 \usepackage[normalem]{ulem}
10 \usepackage{tikz}
11 \usetikzlibrary{positioning,fit}
12 \usepackage{multirow}
13 \usepackage{scalefnt}
14
15 % correct bad hyphenation here
16 \hyphenation{op-tical net-works semi-conduc-tor}
17 \textheight=26cm
18 \setlength{\footskip}{30pt}
19 \pagenumbering{gobble}
20 \begin{document}
21 \title{Filter optimization for real time digital processing of radiofrequency signals: application
22 to oscillator metrology}
23
24 \author{\IEEEauthorblockN{A. Hugeat\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, J. Bernard\IEEEauthorrefmark{2},
25 G. Goavec-M\'erou\IEEEauthorrefmark{1},
26 P.-Y. Bourgeois\IEEEauthorrefmark{1}, J.-M. Friedt\IEEEauthorrefmark{1}}
27 \IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, Time \& Frequency department, Besan\c con, France }
28 \IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Computer Science department DISC, Besan\c con, France \\
29 Email: \{pyb2,jmfriedt\}@femto-st.fr}
30 }
31 \maketitle
32 \thispagestyle{plain}
33 \pagestyle{plain}
34 \newtheorem{definition}{Definition}
35
36 \begin{abstract}
37 Software Defined Radio (SDR) provides stability, flexibility and reconfigurability to
38 radiofrequency signal processing. Applied to oscillator characterization in the context
39 of ultrastable clocks, stringent filtering requirements are defined by spurious signal or
40 noise rejection needs. Since real time radiofrequency processing must be performed in a
41 Field Programmable Array to meet timing constraints, we investigate optimization strategies
42 to design filters meeting rejection characteristics while limiting the hardware resources
43 required and keeping timing constraints within the targeted measurement bandwidths.
44 \end{abstract}
45
46 \begin{IEEEkeywords}
47 Software Defined Radio, Mixed-Integer Linear Programming, Finite Impulse Response filter
48 \end{IEEEkeywords}
49
50 \section{Digital signal processing of ultrastable clock signals}
51
52 Analog oscillator phase noise characteristics are classically performed by downconverting
53 the radiofrequency signal using a saturated mixer to bring the radiofrequency signal to baseband,
54 followed by a Fourier analysis of the beat signal to analyze phase fluctuations close to carrier. In
55 a fully digital approach, the radiofrequency signal is digitized and numerically downconverted by
56 multiplying the samples with a local numerically controlled oscillator (Fig. \ref{schema}) \cite{rsi}.
57
58 \begin{figure}[h!tb]
59 \begin{center}
60 \includegraphics[width=.8\linewidth]{images/schema}
61 \end{center}
62 \caption{Fully digital oscillator phase noise characterization: the Device Under Test
63 (DUT) signal is sampled by the radiofrequency grade Analog to Digital Converter (ADC) and
64 downconverted by mixing with a Numerically Controlled Oscillator (NCO). Unwanted signals
65 and noise aliases are rejected by a Low Pass Filter (LPF) implemented as a cascade of Finite
66 Impulse Response (FIR) filters. The signal is then decimated before a Fourier analysis displays
67 the spectral characteristics of the phase fluctuations.}
68 \label{schema}
69 \end{figure}
70
71 As with the analog mixer,
72 the non-linear behavior of the downconverter introduces noise or spurious signal aliasing as
73 well as the generation of the frequency sum signal in addition to the frequency difference.
74 These unwanted spectral characteristics must be rejected before decimating the data stream
75 for the phase noise spectral characterization \cite{andrich2018high}. The characteristics introduced between the
76 downconverter
77 and the decimation processing blocks are core characteristics of an oscillator characterization
78 system, and must reject out-of-band signals below the targeted phase noise -- typically in the
79 sub -170~dBc/Hz for ultrastable oscillator we aim at characterizing. The filter blocks will
80 use most resources of the Field Programmable Gate Array (FPGA) used to process the radiofrequency
81 datastream: optimizing the performance of the filter while reducing the needed resources is
82 hence tackled in a systematic approach using optimization techniques. Most significantly, we
83 tackle the issue by attempting to cascade multiple Finite Impulse Response (FIR) filters with
84 tunable number of coefficients and tunable number of bits representing the coefficients and the
85 data being processed.
86
87 \section{Finite impulse response filter}
88
89 We select FIR filter for their unconditional stability and ease of design. A FIR filter is defined
90 by a set of weights $b_k$ applied to the inputs $x_k$ through a convolution to generate the
91 outputs $y_k$
92 \begin{align}
93 y_n=\sum_{k=0}^N b_k x_{n-k}
94 \label{eq:fir_equation}
95 \end{align}
96
97 As opposed to an implementation on a general purpose processor in which word size is defined by the
98 processor architecture, implementing such a filter on an FPGA offer more degrees of freedom since
99 not only the coefficient values and number of taps must be defined, but also the number of bits
100 defining the coefficients and the sample size. For this reason, and because we consider pipeline
101 processing (as opposed to First-In, First-Out FIFO memory batch processing) of radiofrequency
102 signals, High Level Synthesis (HLS) languages \cite{kasbah2008multigrid} are not considered but
103 the problem is tackled at the Very-high-speed-integrated-circuit Hardware Description Language (VHDL) level.
104 Since latency is not an issue in a openloop phase noise characterization instrument, the large
105 numbre of taps in the FIR, as opposed to the shorter Infinite Impulse Response (IIR) filter,
106 is not considered as an issue as would be in a closed loop system.
107
108 The coefficients are classically expressed as floating point values. However, this binary
109 number representation is not efficient for fast arithmetic computation by an FPGA. Instead,
110 we select to quantify these floating point values into integer values. This quantization
111 will result in some precision loss.
112
113 \begin{figure}[h!tb]
114 \includegraphics[width=\linewidth]{images/demo_filtre}
115 \caption{Impact of the quantization resolution of the coefficients: the quantization is
116 set to 6~bits -- with the horizontal black lines indicating $\pm$1 least significant bit -- setting
117 the 30~first and 30~last coefficients out of the initial 128~band-pass
118 filter coefficients to 0 (red dots).}
119 \label{float_vs_int}
120 \end{figure}
121
122 The tradeoff between quantization resolution and number of coefficients when considering
123 integer operations is not trivial. As an illustration of the issue related to the
124 relation between number of fiter taps and quantization, Fig. \ref{float_vs_int} exhibits
125 a 128-coefficient FIR bandpass filter designed using floating point numbers (blue). Upon
126 quantization on 6~bit integers, 60 of the 128~coefficients in the beginning and end of the
127 taps become null, making the large number of coefficients irrelevant and allowing to save
128 processing resource by shrinking the filter length. This tradeoff aimed at minimizing resources
129 to reach a given rejection level, or maximizing out of band rejection for a given computational
130 resource, will drive the investigation on cascading filters designed with varying tap resolution
131 and tap length, as will be shown in the next section. Indeed, our development strategy closely
132 follows the skeleton approach \cite{crookes1998environment, crookes2000design, benkrid2002towards}
133 in which basic blocks are defined and characterized before being assembled \cite{hide}
134 in a complete processing chain. In our case, assembling the filter blocks is a simpler block
135 combination process since we assume a single value to be processed and a single value to be
136 generated at each clock cycle. The FIR filters will not be considered to decimate in the
137 current implementation: the decimation is assumed to be located after the FIR cascade at the
138 moment.
139
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.
144
145 For the moment, we are focused on only two operations: the filtering and the shift of data.
146 We have chosen this basic operation because the shifting and the filtering have already be studied in
147 lot of works {\color{red} mettre les nouvelles référence ici} hence it will be easier
148 to check and validate our results.
149
150 However having only two operations is insufficient to work with complex DSP but
151 in this paper we only want demonstrate the relevance and the efficiency of our approach.
152 In future work it will be possible to add more operations and we are able to
153 model any DSP chain.
154
155 We will apply our methodology on very simple DSP chain. We generate a digital signal
156 thanks at generator of Pseudo-Random Number (PRN) or thanks at an Analog to Digital
157 Converter (ADC). Once we have a digital signal, we filter it to decrease the noise level.
158 Finally we stored some burst of filtered samples before post-processing it.
159 % TODO: faire un schéma
160 In this particular case, we want optimize the filtering step to have the best noise
161 rejection for constrain number of resource or to have the minimal resources
162 consumption for a given rejection objective.
163
164 The first step of our approach is to model the DSP chain and since we just optimize
165 the filtering, we have not modeling the PRN generator or the ADC. The filtering can be
166 done by two ways. The first one we use only one FIR filter with lot of coefficients
167 to rejection the noise, we called this approach a monolithic approach. And the second one
168 we select different FIR filters with less coefficients the monolithic filter and we cascaded
169 it to filtering the signal.
170
171 After each filter we leave the possibility of shifting the filtered data to consume
172 less resources. Hence in the case of cascaded filter, we define a stage as a filter
173 and a shifter (the shift could be omitted if we do not need to divide the filtered data).
174
175 \subsection{Model of a FIR filter}
176 A cascade of filter are composed of $n$ stage. In stage $i$ ($1 \leq i \leq n$)
177 the FIR has $C_i$ coefficients and each coefficients are integer values with $\pi^C_i$
178 bits and the filtered data are shifted of $\pi^S_i$ bits. We define also $\pi^-_i$ as
179 the size of input data and $\pi^+_i$ as the size of output data. The figure~\ref{fig:fir_stage}
180 shows a filtering stage.
181
182 \begin{figure}
183 \centering
184 \begin{tikzpicture}[node distance=2cm]
185 \node[draw,minimum size=1.3cm] (FIR) { $C_i, \pi_i^C$ } ;
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] { } ;
189
190 \node[draw,fit=(FIR) (Shift)] (Filter) { } ;
191
192 \draw[->] (Start) edge node [above] { $\pi_i^-$ } (FIR) ;
193 \draw[->] (FIR) -- (Shift) ;
194 \draw[->] (Shift) edge node [above] { $\pi_i^+$ } (End) ;
195 \end{tikzpicture}
196 \caption{A single filter is composed of a FIR (on the left) and a Shifter (on the right)}
197 \label{fig:fir_stage}
198 \end{figure}
199
200 FIR $i$ can reject $F(C_i, \pi_i^C)$ dB. $F$ is determined numerically.
201 To measure this rejection, we use GNU Octave software to design FIR filter coefficients thanks to two
202 algorithms (\texttt{firls} and \texttt{fir1}).
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.
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.
207
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.
211
212 \begin{figure}
213 \centering
214 \begin{tikzpicture}[scale=0.3]
215 \draw[<->] (0,15) -- (0,0) -- (21,0) ;
216 \draw[thick] (0,12) -- (8,12) -- (20,0) ;
217
218 \draw (0,14) node [left] { $P$ } ;
219 \draw (20,0) node [below] { $f$ } ;
220
221 \draw[>=latex,<->] (0,14) -- (8,14) ;
222 \draw (4,14) node [above] { passband } node [below] { $40\%$ } ;
223
224 \draw[>=latex,<->] (8,14) -- (12,14) ;
225 \draw (10,14) node [above] { transition } node [below] { $20\%$ } ;
226
227 \draw[>=latex,<->] (12,14) -- (20,14) ;
228 \draw (16,14) node [above] { stopband } node [below] { $40\%$ } ;
229
230 \draw[>=latex,<->] (16,12) -- (16,8) ;
231 \draw (16,10) node [right] { rejection } ;
232
233 \draw[dashed] (8,-1) -- (8,14) ;
234 \draw[dashed] (12,-1) -- (12,14) ;
235
236 \draw[dashed] (8,12) -- (16,12) ;
237 \draw[dashed] (12,8) -- (16,8) ;
238
239 \end{tikzpicture}
240
241 % \includegraphics[width=.5\linewidth]{images/fir_magnitude}
242 \caption{Shape of the filter transmitted power $P$ as a function of frequency $f$:
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}
246 \end{figure}
247
248 In the transition band, the behavior of the filter is left free, we only care about the passband and the stopband.
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.
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}.
251
252 \begin{figure}
253 \centering
254 \includegraphics[width=\linewidth]{images/mean_criterion}
255 \caption{Mean criterion comparison between monolithic filter and cascade filters}
256 \label{fig:mean_criterion}
257 \end{figure}
258
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}
265
266 Although we have a efficient criterion to estimate the rejection of one set of coefficient
267 we have a problem when we sum two or more criterion. If the FIR filter coefficients are the same
268 between the stage, we have:
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.
278
279 \begin{figure}
280 \centering
281 \includegraphics[width=\linewidth]{images/sum_rejection}
282 \caption{Rejection of two cascaded filters}
283 \label{fig:sum_rejection}
284 \end{figure}
285
286 \section{Experiments with fixed area space}
287
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}
294
295 \begin{figure}
296 \centering
297 \includegraphics[width=\linewidth]{images/max_rejection/prn_1000}
298 \caption{Experimental results for design with PRN as data input and 1000 a.u. as max arbitrary space}
299 \label{fig:prn_1000}
300 \end{figure}
301
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}
308
309 \begin{table}
310 \centering
311 \begin{tabular}{|c|c|ccc|c|c|}
312 \hline
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}
328
329 \begin{table}
330 \centering
331 {\scalefont{0.85}
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}
356
357 \begin{table}
358 \centering
359 {\scalefont{0.85}
360 \begin{tabular}{|c|c|ccccc|c|c|}
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}
384
385 \begin{table}
386 \centering
387 \begin{tabular}{|c|c|c|c|c|}\hline
388 Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
389 & 1 & 0.02~s & $\approx$ 20 min & $\approx$ 1 min \\
390 PRN & 2 & 1.70~s & $\approx$ 20 min & $\approx$ 1 min \\
391 & 3 & 19~s & $\approx$ 20 min & $\approx$ 1 min \\\hline
392 \end{tabular}
393 \caption{Time to compute and deploy the designs for PRN 500}
394 \label{tbl:time_prn_500}
395 \end{table}
396
397 \begin{table}
398 \centering
399 \begin{tabular}{|c|c|c|c|c|}\hline
400 Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
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}
410
411 \begin{table}
412 \centering
413 \begin{tabular}{|c|c|c|c|c|}\hline
414 Input & Stages & Computation time & Vivado time & Redpitaya time \\\hline\hline
415 & 1 & 0.07~s & $\approx$ 20 min & $\approx$ 1 min \\
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}
424
425 \section{Experiments with fixed rejection target}
426
427 \begin{figure}
428 \centering
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}
433
434 \begin{figure}
435 \centering
436 \includegraphics[width=\linewidth]{images/min_area/prn_100}
437 \caption{Results for design with PRN as data input and 50 dB as aimed rejection level}
438 \label{fig:prn_100}
439 \end{figure}
440
441 \begin{figure}
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}
447
448 \section{Conclusion}
449
450 \section*{Acknowledgement}
451
452 This work is supported by the ANR Programme d'Investissement d'Avenir in
453 progress at the Time and Frequency Departments of the FEMTO-ST Institute
454 (Oscillator IMP, First-TF and Refimeve+), and by R\'egion de Franche-Comt\'e.
455 The authors would like to thank E. Rubiola, F. Vernotte, and G. Cabodevila
456 for support and fruitful discussions.
457
458 \bibliographystyle{IEEEtran}
459 \balance
460 \bibliography{references,biblio}
461 \end{document}
% JMF : revoir l'abstract : on y avait mis le Zynq7010 de la redpitaya en montrant 1 462
% 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
\usepackage{graphicx,color,hyperref} 7
\usepackage{amsfonts} 8
\usepackage{amsthm} 9
\usepackage{amssymb} 10
\usepackage{amsmath} 11
\usepackage{algorithm2e} 12
\usepackage{url,balance} 13
\usepackage[normalem]{ulem} 14
% correct bad hyphenation here 15
\hyphenation{op-tical net-works semi-conduc-tor} 16
\textheight=26cm 17
\setlength{\footskip}{30pt} 18
\pagenumbering{gobble} 19
\begin{document} 20
\title{Filter optimization for real time digital processing of radiofrequency signals: application 21
to oscillator metrology} 22
23
\author{\IEEEauthorblockN{A. Hugeat\IEEEauthorrefmark{1}\IEEEauthorrefmark{2}, J. Bernard\IEEEauthorrefmark{2}, 24
G. Goavec-M\'erou\IEEEauthorrefmark{1}, 25
P.-Y. Bourgeois\IEEEauthorrefmark{1}, J.-M. Friedt\IEEEauthorrefmark{1}} 26
\IEEEauthorblockA{\IEEEauthorrefmark{1}FEMTO-ST, Time \& Frequency department, Besan\c con, France } 27
\IEEEauthorblockA{\IEEEauthorrefmark{2}FEMTO-ST, Computer Science department DISC, Besan\c con, France \\ 28
Email: \{pyb2,jmfriedt\}@femto-st.fr} 29
} 30
\maketitle 31
\thispagestyle{plain} 32
\pagestyle{plain} 33
\newtheorem{definition}{Definition} 34
35
\begin{abstract} 36
Software Defined Radio (SDR) provides stability, flexibility and reconfigurability to 37
radiofrequency signal processing. Applied to oscillator characterization in the context 38
of ultrastable clocks, stringent filtering requirements are defined by spurious signal or 39
noise rejection needs. Since real time radiofrequency processing must be performed in a 40
Field Programmable Array to meet timing constraints, we investigate optimization strategies 41
to design filters meeting rejection characteristics while limiting the hardware resources 42
required and keeping timing constraints within the targeted measurement bandwidths. 43
\end{abstract} 44
45
\begin{IEEEkeywords} 46
Software Defined Radio, Mixed-Integer Linear Programming, Finite Impulse Response filter 47
\end{IEEEkeywords} 48
49
\section{Digital signal processing of ultrastable clock signals} 50
51
Analog oscillator phase noise characteristics are classically performed by downconverting 52
the radiofrequency signal using a saturated mixer to bring the radiofrequency signal to baseband, 53
followed by a Fourier analysis of the beat signal to analyze phase fluctuations close to carrier. In 54
a fully digital approach, the radiofrequency signal is digitized and numerically downconverted by 55
multiplying the samples with a local numerically controlled oscillator (Fig. \ref{schema}) \cite{rsi}. 56
57
\begin{figure}[h!tb] 58
\begin{center} 59
\includegraphics[width=.8\linewidth]{images/schema} 60
\end{center} 61
\caption{Fully digital oscillator phase noise characterization: the Device Under Test 62
(DUT) signal is sampled by the radiofrequency grade Analog to Digital Converter (ADC) and 63
downconverted by mixing with a Numerically Controlled Oscillator (NCO). Unwanted signals 64
and noise aliases are rejected by a Low Pass Filter (LPF) implemented as a cascade of Finite 65
Impulse Response (FIR) filters. The signal is then decimated before a Fourier analysis displays 66
the spectral characteristics of the phase fluctuations.} 67
\label{schema} 68
\end{figure} 69
70
As with the analog mixer, 71
the non-linear behavior of the downconverter introduces noise or spurious signal aliasing as 72
well as the generation of the frequency sum signal in addition to the frequency difference. 73
These unwanted spectral characteristics must be rejected before decimating the data stream 74
for the phase noise spectral characterization \cite{andrich2018high}. The characteristics introduced between the 75
downconverter 76
and the decimation processing blocks are core characteristics of an oscillator characterization 77
system, and must reject out-of-band signals below the targeted phase noise -- typically in the 78
sub -170~dBc/Hz for ultrastable oscillator we aim at characterizing. The filter blocks will 79
use most resources of the Field Programmable Gate Array (FPGA) used to process the radiofrequency 80
datastream: optimizing the performance of the filter while reducing the needed resources is 81
hence tackled in a systematic approach using optimization techniques. Most significantly, we 82
tackle the issue by attempting to cascade multiple Finite Impulse Response (FIR) filters with 83
tunable number of coefficients and tunable number of bits representing the coefficients and the 84
data being processed. 85
86
\section{Finite impulse response filter} 87
88
We select FIR filter for their unconditional stability and ease of design. A FIR filter is defined 89
by a set of weights $b_k$ applied to the inputs $x_k$ through a convolution to generate the 90
outputs $y_k$ 91
$$y_n=\sum_{k=0}^N b_k x_{n-k}$$ 92
93
As opposed to an implementation on a general purpose processor in which word size is defined by the 94
processor architecture, implementing such a filter on an FPGA offer more degrees of freedom since 95
not only the coefficient values and number of taps must be defined, but also the number of bits 96
defining the coefficients and the sample size. For this reason, and because we consider pipeline 97
processing (as opposed to First-In, First-Out FIFO memory batch processing) of radiofrequency 98
signals, High Level Synthesis (HLS) languages \cite{kasbah2008multigrid} are not considered but 99
the problem is tackled at the Very-high-speed-integrated-circuit Hardware Description Language (VHDL) level. 100
Since latency is not an issue in a openloop phase noise characterization instrument, the large 101
numbre of taps in the FIR, as opposed to the shorter Infinite Impulse Response (IIR) filter, 102
is not considered as an issue as would be in a closed loop system. 103
104
The coefficients are classically expressed as floating point values. However, this binary 105
number representation is not efficient for fast arithmetic computation by an FPGA. Instead, 106
we select to quantify these floating point values into integer values. This quantization 107
will result in some precision loss. 108
109
%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
\includegraphics[width=\linewidth]{images/demo_filtre} 125
\caption{Impact of the quantization resolution of the coefficients: the quantization is 126
set to 6~bits -- with the horizontal black lines indicating $\pm$1 least significant bit -- setting 127
the 30~first and 30~last coefficients out of the initial 128~band-pass 128
filter coefficients to 0 (red dots).} 129
\label{float_vs_int} 130
\end{figure} 131
132
The tradeoff between quantization resolution and number of coefficients when considering 133
integer operations is not trivial. As an illustration of the issue related to the 134
relation between number of fiter taps and quantization, Fig. \ref{float_vs_int} exhibits 135
a 128-coefficient FIR bandpass filter designed using floating point numbers (blue). Upon 136
quantization on 6~bit integers, 60 of the 128~coefficients in the beginning and end of the 137
taps become null, making the large number of coefficients irrelevant and allowing to save 138
processing resource by shrinking the filter length. This tradeoff aimed at minimizing resources 139
to reach a given rejection level, or maximizing out of band rejection for a given computational 140
resource, will drive the investigation on cascading filters designed with varying tap resolution 141
and tap length, as will be shown in the next section. Indeed, our development strategy closely 142
follows the skeleton approach \cite{crookes1998environment, crookes2000design, benkrid2002towards} 143
in which basic blocks are defined and characterized before being assembled \cite{hide} 144
in a complete processing chain. In our case, assembling the filter blocks is a simpler block 145
combination process since we assume a single value to be processed and a single value to be 146
generated at each clock cycle. The FIR filters will not be considered to decimate in the 147
current implementation: the decimation is assumed to be located after the FIR cascade at the 148
moment. 149
150
\section{Filter optimization} 151
152
A basic approach for implementing the FIR filter is to compute the transfer function of 153
a monolithic filter: this single filter defines all coefficients with the same resolution 154
(number of bits) and processes data represented with their own resolution. Meeting the 155
filter shape requires a large number of coefficients, limited by resources of the FPGA since 156
this filter must process data stream at the radiofrequency sampling rate after the mixer. 157
158
An optimization problem \cite{leung2004handbook} aims at improving one or many 159
performance criteria within a constrained resource environment. Amongst the tools 160
developed to meet this aim, Mixed-Integer Linear Programming (MILP) provides the framework to 161
formally define the stated problem and search for an optimal use of available 162
resources \cite{yu2007design, kodek1980design}. 163
164
First we need to ensure that our problem is a real optimization problem. When 165
designing a processing function in the FPGA, we aim at meeting some requirement such as 166
the throughput, the computation time or the noise rejection noise. However, due to limited 167
resources to design the process like BRAM (high performance RAM), DSP (Digital Signal Processor) 168
or LUT (Look Up Table), a tradeoff must be generally searched between performance and available 169
computational resources: optimizing some criteria within finite, limited 170
resources indeed matches the definition of a classical optimization problem. 171
172
Specifically the degrees of freedom when addressing the problem of replacing the single monolithic 173
FIR with a cascade of optimized filters are the number of coefficients $N_i$ of each filter $i$, 174
the number of bits $C_i$ representing the coefficients and the number of bits $D_i$ needed to represent 175
the data $x_k$ fed to each filter as provided by the acquisition or previous processing stage. 176
Because each FIR in the chain is fed the output of the previous stage, 177
the optimization of the complete processing chain within a constrained resource environment is not 178
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
\begin{figure}[h!tb] 196
\begin{center} 197
\includegraphics[width=.5\linewidth]{schema2} 198
\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
Following these considerations, the model is expressed as: 207
\begin{align} 208
\begin{cases} 209
\mathcal{R}_i &= \mathcal{F}(N_i, C_i)\\ 210
\mathcal{A}_i &= N_i \times C_i\\ 211
\Delta_i &= \Delta _{i-1} + \mathcal{P}_i 212
\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
\begin{figure}[h!tb] 229
\includegraphics[width=\linewidth]{images/noise-rejection.pdf} 230
\caption{Rejection as a function of number of coefficients and number of bits} 231
\label{noise-rejection} 232
\end{figure} 233
234
The objective function maximizes the noise rejection ($\max(\Delta_{i_{\max}})$) while keeping resource 235
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
Linear program formalism for solving the problem is well documented: an objective function is 250
defined which is linearly dependent on the parameters to be optimized. Constraints are expressed 251
as linear equations and solved using one of the available solvers, in our case GLPK\cite{glpk}. 252
With the notations used in the description of system \ref{model-FIR}, we have defined the linear problem as: 253
\paragraph{Variables} 254
\begin{align*} 255
x_{i,j} \in \lbrace 0,1 \rbrace & \text{ $i$ is a given filter} \\ 256
& \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
The constraint \ref{cstr_size} means the occupation for the current stage $j$ depends on 290
the previous occupation and the occupation of current selected filter (it is possible 291
that no filter is selected for this stage). And the second one \ref{cstr_rejection} 292
means the same thing but for the rejection, the rejection depends the previous rejection 293
plus the rejection of selected filter. 294
295
\subsection{Low bandpass ripple and maximum rejection criteria} 296
297
The MILP solver provides a solution to the problem by selecting a series of small FIR with 298
increasing number of bits representing data and coefficients as well as an increasing number 299
of coefficients, instead of a single monolithic filter. 300
301
\begin{figure}[h!tb] 302
% \includegraphics[width=\linewidth]{images/compare-fir.pdf} 303
\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
Fig. \ref{compare-fir} exhibits the 310
performance comparison between one solution and a monolithic FIR when selecting a cutoff 311
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
317
The resource occupation when synthesizing such FIR on a Xilinx FPGA is summarized as Tab. \ref{t1}. 318
We have considered a set of resources representative of the hardware platform we work on, 319
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
\begin{table}[h!tb] 331
\caption{Resource occupation on a Xilinx Zynq-7000 series FPGA when synthesizing the FIR cascade 332
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
\subsection{Alternate criteria}\label{median} 354
355
Fig. \ref{compare-fir} provides FIR solutions matching well the targeted transfer 356
function, namely low ripple in the bandpass defined as the first 40\% of the frequency 357
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
\begin{figure}[h!tb] 364
\includegraphics[width=\linewidth]{images/fir-mono-vs-fir-series-noise-fixe-mean-light.pdf} 365
\caption{Comparison of the rejection capability between a series of FIR and a monolithic FIR 366
with a cutoff frequency set at half the Nyquist frequency.} 367
\label{compare-mean} 368
\end{figure} 369
370
In the case of the mean value criterion (Fig. \ref{compare-mean}), the solution is not 371
acceptable since the notch at the end of the transition band compensates for some unacceptable 372
rise in the rejection close to the Nyquist frequency. Applying such a filter might yield excessive 373
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
\begin{figure}[h!tb] 382
\includegraphics[width=\linewidth]{images/fir-mono-vs-fir-series-noise-fixe-median-light.pdf} 383
\caption{Comparison of the rejection capability between a series of FIR and a monolithic FIR 384
with a cutoff frequency set at half the Nyquist frequency.} 385
\label{compare-median} 386
\end{figure} 387
388
\section{Filter coefficient selection} 389
390
The coefficients of a single monolithic filter are computed as the impulse response 391
of the filter transfer function, and practically approximated by a multitude of methods 392
including least square optimization (Matlab's {\tt firls} function), Hamming or Kaiser windowing 393
(Matlab's {\tt fir1} function). 394
395
\begin{figure}[h!tb] 396
\includegraphics[width=\linewidth]{images/fir1-vs-firls} 397
\caption{Evolution of the rejection capability of least-square optimized filters and Hamming 398
FIR filters as a function of the number of coefficients, for floating point numbers and 8-bit 399
encoded integers.} 400
\label{2} 401
\end{figure} 402
403
Cascading filters opens a new optimization opportunity by 404
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
\section{Conclusion} 411
412
We address the optimization problem of designing a low-pass filter chain in a Field Programmable Gate 413
Array for improved noise rejection within constrained resource occupation, as needed for 414
real time processing of radiofrequency signal when characterizing spectral phase noise 415
characteristics of stable oscillators. The flexibility of the digital approach makes the result 416
best suited for closing the loop and using the measurement output in a feedback loop for 417
controlling clocks, e.g. in a quartz-stabilized high performance clock whose long term behavior 418
is controlled by non-piezoelectric resonator (sapphire resonator, microwave or optical 419
atomic transition). 420
421
\section*{Acknowledgement} 422
423
This work is supported by the ANR Programme d'Investissement d'Avenir in 424
progress at the Time and Frequency Departments of the FEMTO-ST Institute 425
(Oscillator IMP, First-TF and Refimeve+), and by R\'egion de Franche-Comt\'e. 426
The authors would like to thank E. Rubiola, F. Vernotte, and G. Cabodevila 427
for support and fruitful discussions. 428
429
\bibliographystyle{IEEEtran} 430
\balance 431
\bibliography{references,biblio} 432
\end{document} 433
434
\section{Contexte d'ordonnancement} 435
Dans cette partie, nous donnerons des d\'efinitions de termes rattach\'es au domaine de l'ordonnancement 436
et nous verrons que le sujet trait\'e se rapproche beaucoup d'un problème d'ordonnancement. De ce fait 437
nous pourrons aller plus loin que les travaux vus pr\'ec\'edemment et nous tenterons des approches d'ordonnancement 438
et d'optimisation. 439
440
\subsection{D\'efinition du vocabulaire} 441
Avant tout, il faut d\'efinir ce qu'est un problème d'optimisation. Il y a deux d\'efinitions 442
importantes à donner. La première est propos\'ee par Legrand et Robert dans leur livre \cite{def1-ordo} : 443
\begin{definition} 444
\label{def-ordo1} 445
Un ordonnancement d'un système de t\^aches $G\ =\ (V,\ E,\ w)$ est une fonction $\sigma$ : 446
$V \rightarrow \mathbb{N}$ telle que $\sigma(u) + w(u) \leq \sigma(v)$ pour toute arête $(u,\ v) \in E$. 447
\end{definition} 448
449
Dit plus simplement, l'ensemble $V$ repr\'esente les t\^aches à ex\'ecuter, l'ensemble $E$ repr\'esente les d\'ependances 450
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
chacune des t\^aches. La d\'efinition dit que si une t\^ache $v$ d\'epend d'une t\^ache $u$ alors 452
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
temps d'ex\'ecution. 454
455
Une autre d\'efinition importante qui est propos\'ee par Leung et al. \cite{def2-ordo} est : 456
\begin{definition} 457
\label{def-ordo2} 458
L'ordonnancement traite de l'allocation de ressources rares à des activit\'es avec 459
l'objectif d'optimiser un ou plusieurs critères de performance. 460
\end{definition} 461
462
Cette d\'efinition est plus g\'en\'erique mais elle nous int\'eresse d'avantage que la d\'efinition \ref{def-ordo1}. 463
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
En revanche la d\'efinition \ref{def-ordo2} sera au c\oe{}ur du projet. Pour se convaincre de cela, 467
il nous faut d'abord d\'efinir quel est le type de problème d'ordonnancement qu'on traite et quelles 468
sont les m\'ethodes qu'on peut appliquer. 469
470
Les problèmes d'ordonnancement peuvent être class\'es en diff\'erentes cat\'egories : 471
\begin{itemize} 472
\item T\^aches ind\'ependantes : dans cette cat\'egorie de problèmes, les t\^aches sont complètement ind\'ependantes 473
les unes des autres. Dans notre cas, ce n'est pas le plus adapt\'e. 474
\item Graphe de t\^aches : la d\'efinition \ref{def-ordo1} d\'ecrit cette cat\'egorie. La plupart du temps, 475
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
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
Bien entendu, cette liste n'est pas exhaustive et il existe de nombreuses autres classifications et sous-classifications 485
de ces problèmes. Nous n'avons parl\'e ici que des cat\'egories les plus communes. 486
487
Un autre point à d\'efinir, est le critère d'optimisation. Il y a là encore un grand nombre de 488
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
En plus de cela, on peut avoir besoin de plusieurs critères d'optimisation. Il s'agit dans ce cas d'une optimisation 500
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
\subsection{Formalisation du problème} 505
\label{formalisation} 506
Maintenant que nous avons donn\'e le vocabulaire li\'e à l'ordonnancement, nous allons pouvoir essayer caract\'eriser 507
formellement notre problème. En effet, nous allons reprendre les contraintes \'enonc\'ees dans la sections \ref{def-contraintes} 508
et nous essayerons de les formaliser le plus finement possible. 509
510
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
que nous nommerons $\mathcal{P}_{i}$. Cet ensemble $\mathcal{P}_i$ est propre à chaque t\^ache et il variera d'une 512
t\^ache à l'autre. Nous reviendrons plus tard sur les paramètres qui peuvent composer cet ensemble. 513
514
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
images/custom_criterion.pdf
No preview for this file type
images/max_rejection/prn_1000.pdf
No preview for this file type
images/max_rejection/prn_2000.pdf
No preview for this file type
images/max_rejection/prn_500.pdf
No preview for this file type
images/mean_criterion.pdf
No preview for this file type
images/min_area/prn_100.pdf
No preview for this file type
images/min_area/prn_150.pdf
No preview for this file type
images/min_area/prn_50.pdf
No preview for this file type
images/sum_rejection.pdf
No preview for this file type