S Meunier

D Meneveaux

L Aveneau

D Ghazanfarpour

Cosine Lobes for Interactive Direct Lighting in Dynamic Scenes

Keywords: (according to ACM CCS): I.3.3 [Computer Graphics]: Bitmap and frame buffer operations I.3.7 [Computer Graphics]: Color, shading, shadowing, and texture, XLim laboratory

Many rendering systems rely on spherical harmonics or wavelets that provide a mean to solve the rendering equation using scalar products. In addition, these basis of functions may represent hemispherical lighting, visibility and reflectance with a small number of coefficients. However, their use requires the projection of each term of the rendering equation, which is computationally intensive at run-time (for instance with dynamic environments), and limits practically the number of basis functions used at the expense of precision. In this paper, we show how cosine lobes can also be used for representing each term of the rendering equation, without the drawbacks usually existing in the above basis of functions. Cosine lobe representations (such as Lambert, Phong or Lafortune models) have already been used intensively by many authors for their advantage of intuitively and compactly representing reflectance functions. Using this representation also for visibility and lighting, we explain how the rendering equation can be efficiently solved. As an application, we propose an interactive rendering system for direct lighting, naturally including soft shadows and spatially varying materials.

Introduction

Several methods in the literature are proposed for producing physically-based images with lighting simulations, complex BRDF (bidirectional reflectance distribution function), natural lighting, a high number of objects and so on [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF]. Unfortunately, producing realistic images at interactive frame rates in dynamic scenes remains a challenging task. Many authors have provided interesting solutions with the use of orthogonal basis functions such as spherical harmonics (SH) [START_REF] Kautz | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF], hemispherical harmonics (HSH) [START_REF] Gautron | Cache de luminance et cartes graphiques : une approche pour la simulation d'éclairage temps réel dans des scènes animées[END_REF] or Haar wavelets (HW) [START_REF] Ng R | All-frequency shadows using non-linear wavelet lighting approximation[END_REF]. This type of representation simplifies the rendering equation computation and replaces the continuous integration by products of coefficients. Nevertheless, the projection of visibility, BRDF and/or incoming radiance on the basis functions remains a costly process. This problem is often addressed offline, with pre-computed radiance transfer methods (PRT) [START_REF] Kautz | Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments[END_REF][START_REF] Ng R | All-frequency shadows using non-linear wavelet lighting approximation[END_REF], so as to provide real time rendering. However, PRT still is difficult to manage with fully dynamic scenes because the stored data increases rapidly, unless some parameters be fixed (e.g. viewpoint, objects, BRDF, etc.). Some applications compute only a few coefficients for ensuring interactivity rather than relying on pre-computed transfers [START_REF] Kautz | Hemispherical rasterization for self-shadowing of dynamic objects[END_REF]. These methods inhibit all frequency lighting environments with SH and produce blocking artifacts with HW.

In this paper we propose a method for overcoming the limitations of previous methods, using the flexibility of spherical radial basis functions (SRBF) [START_REF] Tsai Y.-T | All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation[END_REF]. They are defined on the sphere, not necessarily uniformly distributed, invariant under rotation and spatially localized; they only require choosing an axis and some shape parameters which is much meaningful; the representation is compact (a small number of coefficients accurately represent data). Our approach is based on cosine lobes, providing a mean for using directly Phong and Lafortune BRDF models, with potentially spatially varying, anisotropic, retro-reflection and/or off-specular reflectance data [START_REF] Mcallister | Efficient rendering of spatial bi-directional reflectance distribution functions[END_REF]. It provides practical solutions to represent many reflectance properties. Our contributions include:

• a method for solving rendering equation using a homogeneous representation based on cosine lobes; • an acceleration system based on table fetching and linear interpolations; • an application using cosine lobes for estimating shadows and direct lighting with various types of BRDF; • an interactive rendering system providing all-frequency shadows and realistic materials in dynamic scenes.

This paper shows how the rendering equation can be efficiently solved when each term is represented by cosine lobes. Our application does not address diffuse interreflections that should be handled using existing approximations proposed in [START_REF] Govindaraju | Image-based proxy accumulation for real-time soft global illumination[END_REF][START_REF] Guerrero | Real-time indirect illumination and soft shadows in dynamic scenes using spherical lights[END_REF].

The remaining of this paper is organized as follows: Section 2 reviews existing methods and highlights our contributions; Section 3 presents an overview of our approach. Section 4 details the use of cosine lobes for the rendering equation; Section 5 explains the representation of visibility, BRDF and incoming radiance using cosine lobe sums; Section 6 presents and discusses the performances provided by our method; Finally, Section 7 concludes and mentions future work.

Related Work

Let us consider the rendering equation :

Lo(x, -→ vo ) = Z !+ fr(x, -→ v , -→ vo )L i (x, -→ v )( -→ n • -→ v )d -→ v (1) = Z !+ fr(x, -→ v , -→ vo )L i (x, -→ v )d -→ v (2)
where Lo denotes the outgoing radiance from a point x in direction -→ vo ; !+ is the upper hemisphere with solid angles d -→ v around directions -→ v ; fr is the BRDF and fr the BRDF multiplied by ( -→ n • -→ v ); L i is the incoming radiance at x from direction -→ v and n the surface normal at x.

Let Ls be the environment lighting at a point x from any direction -→ v . When objects are located between x and the lighting environment, visibility has to be considered. Equation 3 can thus be rewritten :

Lo(x, -→ vo ) = Z !+ fr(x, -→ v , -→ vo )Ls(x, -→ v )V (x, -→ v )d -→ v (3)
V being the visibility term.

We aim at dealing with dynamic scenes (with potentially varying BRDF, moving lighting sources and objects), so that the triple product fr • Ls • V be solved at run-time. Several methods choose SH or HW basis of functions for representing fr, Ls and V , such as [KLA04, ZHL * 05, RWS * 06, SGNS07, GJW08].

Generally speaking, orthogonal basis functions such as SH and HW provide a useful representation of functions on the sphere. Once the projection is performed, the integration becomes a dot product, allowing real-time or interactive rendering for dynamic scenes. However, all the coefficients have to be estimated even though some (potentially many) of them only contribute negligibly to lighting computations [START_REF] Kautz | Hemispherical rasterization for self-shadowing of dynamic objects[END_REF].

Shadow fields [ZHL * 05] tabulate SH or HW visibility representation in the space surrounding each blocker, requiring large tables and many triple product computations. SH exponentiation [RWS * 06,SGNS07,GJW08] replaces the SH triple products by fast coefficients accumulation in log space. With bounding sphere sets, these methods are capable of handling more complex scenes and skinned objects. Nevertheless, SH limit the application to low frequency lighting environments. Kautz et al. [START_REF] Kautz | Hemispherical rasterization for self-shadowing of dynamic objects[END_REF] propose a hemispherical rasterization of blocker geometry for estimating a SH representation of the transfer function ( fr.V in this case). This method also uses SH with low memory requirements and no precomputation based on the scene geometry. [KLA04,ZHL * 05,RWS * 06,SGNS07,GJW08] are commonly denoted as blocker accumulation methods and imply visibility computations for each shaded point. [START_REF] Kozlowski | Is accurate occlusion of glossy reflections necessary? In APGV '07[END_REF] perceptually analyze possible simplifications for visibility. They propose a directional ambient occlusion method based on a piecewise constant approximation of V . SH are used to smooth the results and even with approximated shadows and highlights, images remain visually plausible. Green et al. [START_REF] Green | View-dependent precomputed light transport using nonlinear gaussian function approximations[END_REF][START_REF] Green | Efficient reflectance and visibility approximations for environment map rendering[END_REF] use a similar simplification for shadows, associated with isotropic gaussian functions (SRBF) for representing BRDF and/or visibility; incoming radiance is represented by environment maps, pre-filtered using gaussian functions, allowing real time rendering. These methods have been designed for existing environment maps and do not allow interactive dynamic lighting.

Kozlowski and Kautz

SRBF are well known in the computer graphics community. For instance the generalized cosine lobe [START_REF] Lafortune | Non-linear approximation of reflectance functions[END_REF] and the isotropic Gaussian kernel [START_REF] Ward | Measuring and modeling anisotropic reflection[END_REF] are used to model BRDF, leading to a compact, expressive and physically plausible representation. Tsai and Shih [START_REF] Tsai Y.-T | All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation[END_REF] adapt the Abel-Poisson kernel to PRT. The compactness of SRBF leads to a more efficient compression scheme than the clustered principal component analysis proposed in [START_REF] Hall | Clustered principal components for precomputed radiance transfer[END_REF] and provides faster rendering with more realistic results.

The attractiveness of SH/HW is practically limited by the projection cost for high-frequency environments or BRDF. For instance a small detail may require a lot of (sometimes null) coefficients with SH/HW when only one cosine lobe can be efficiently used. In addition, some useless computations have to be performed during the SH/HW projection, when null coefficients are produced. We propose to overcome this limitation with the use of an unfixed (contrary to [START_REF] Tsai Y.-T | All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation[END_REF]) set of cosine lobes for every term : visibility, lighting and reflectance.

Technical Report, XLim laboratory

Overview

We propose to represent each term of the triple product fr • Ls • V by cosine lobes in the rendering equation (Equation 3). The integration is performed based on the idea that the product of two lobes can be approximated by another lobe (Sections 4). We provide an interactive application dedicated to direct lighting (Section 5), where our approximation proves efficient with freely and dynamically distributed lobes, avoiding useless computations of null coefficients (Section 6).

We have implemented a blocker accumulation system that benefits from SRBF (with cosine lobes) advantages, while providing interactive rendering in dynamic scenes with allfrequency BRDF, visibility and incoming radiance. Therefore, scene geometry is simplified using sphere sets (blockers and light sources); from a given point, a sphere defines a cone, used for constructing a lobe.

The triple product fr • Ls • V cosine lobes are generated using:

• the BRDF (Lambert, Phong, Lafortune, etc.) representation for fr (see Section 5.1); • a set of spherical light sources for Ls, with one lobe for each sphere (see Section 5.2); • for V , a visibility mask is produced using the simplified objects geometry for constructing again spheres and thus cosine lobes (see section 5.3).

The Rendering Equation and Cosine Lobes

A cosine lobe function is written:

c( -→ v ) = s max( -→ a • -→ v , 0) e (4) = s "( -→ v ) ( 5
)
where s is the scaling factor, -→ a is the axis and e is the lobe thickness.

Fixing x, -→ vo and assuming fr, Ls and V approximated by any cosine lobe sums: fr(

-→ v ) = # i s i " i ( -→ v ) (6) Ls( -→ v ) = # j s j " j ( -→ v ) (7) V ( -→ v ) = # k s k " k ( -→ v ) (8)
Note that s i " i , s j " j and s k " k do generally not correspond to the same cosine lobes (even when i = j = k).

Equation 3 can be rewritten: However, the term

Lo = Z ! # i s i " i ( -→ v ) # j s j " j ( -→ v ) # k s k " k ( -→ v ) d -→ v (9) = # i # j # k s i s j s k Z ! " i ( -→ v )" j ( -→ v )" k ( -→ v ) d -→ v (10) Y X Z Y X Z X Y Z X Y Z X Z X Z
R ! " i ( -→ v )" j ( -→ v )" k ( -→ v ) d -→ v
cannot be precomputed because cosine lobes are not fixed. Thus, starting from equation 9:

Lo = Z ! # i, j s i " i ( -→ v )s j " j ( -→ v ) # k s k " k ( -→ v ) d -→ v (11) we propose to approximate s i " i ( -→ v )s j " j ( -→ v ) by a single cosine lobe s l " l ( -→ v ): Lo ≈ Z ! # l s l " l ( -→ v ) # k s k " k ( -→ v ) d -→ v (12) = Z ! # l,k s l " l ( -→ v )s k " k ( -→ v ) d -→ v (13)
The rendering equation finally becomes a sum of cosine lobe integrations:

Lo ≈ Z ! # m sm"m( -→ v ) d -→ v = # m sm Z ! "m( -→ v ) d -→ v (14)
Sections 4.1 and 4.2 describe how the product of two lobes is managed and how the cosine lobe integration can be precomputed.

Product

We propose to approximate the product of two cosine lobes c 1 and c 2 by another single cosine lobe cr:

c 1 ( -→ v )c 2 ( -→ v ) ≈ cr( -→ v ) (15) ≈ s 1 s 2 pr max( -→ ar • -→ v , 0) er (16)
pr is called the partial scaling factor. The parameters pr, -→ ar and er depend on e 1 , e 2 and $ (the angle between -→ a 1 and -→ a 2 ). pr, -→ ar and er are precomputed for a wide range of (e 1 , e 2 , $) triplets and store the result in a table. For maintaining directionality and shape for the approximated cosine lobe (see Figure 1), we choose to minimize the

L 2 distance between c 1 ( -→ v ) • c 2 ( -→ v ) and cr( -→ v ): {pr, -→ ar , er} = arg min {pr, - → ar ,er} Z ! c 1 ( -→ v )c 2 ( -→ v ) -cr( -→ v ) 2 d -→ v (17)
Since solving this optimization problem for the parameters altogether is time-consuming and prone to numerical errors, we propose an approach in two steps: (i) finding the axis -→ ar and (ii) fitting the partial scaling factor pr and the exponent er.

First step:

{ar} = arg max { - → ar } c 1 ( -→ ar )c 2 ( -→ ar ) (18)
is done with the Nelder-Mead method. The algorithm is initialized with

- → a1 + - → a2 - → a1 + - → a2 since -→ ar is necessary between -→ a 1 and -→ a 2 .

Second step:

{pr, er} = arg min {pr ,er}

Z ! c 1 ( -→ v )c 2 ( -→ v ) -cr( -→ v ) 2 d -→ v (19)
is done with the Levenberg-Marquardt algorithm. er is initialized with e 1 + e 2 which is the solution for -→ a 1 = -→ a 2 and pr with 1.

In practice e 1 and e 2 are sampled in [0, 100000] so that any scene type be handled. $ is ranged in [0, %]. Our sampling strategy is based on the cosine lobe integration variation (see Figure 2). We pick 16 linearly spaced samples in each interval [0, 1], [1, 10]. . . , [10000, 100000], which is empirically an acceptable trade-off between storage, computing time and accuracy. Such an approximation produces precise results compared to the actual lobe product integration: we have estimated the integration difference value between (i) the approximated lobe product in the table and (ii) the corresponding lobe product. The integration is processed using an adaptive Simpson quadrature with a tolerance set to 10 -7 ; the mean error is equal to 7.1379 × 10 -4 , with a variance equal to 1.2091 × 10 -8 ; the median error value is 0. Note that the maximum error (Figure 3) occurs with very low exponents e 1 and e 2 , that are practically rarely used, even with Lafortune model. In addition, as explained in Section 4.3, cosine lobe exponents produced during rendering always increase.

Integration

Note that the cosine lobe integration value only depends on the exponent:

Z ! c( -→ v ) d -→ v = Z ! s max( -→ a • -→ v , 0) e d -→ v (20) = s Z ! max( -→ X • -→ v , 0) e d -→ v ( 21
)
where c is an arbitrary lobe and -→ X is any axis.

Since this integration is time-consuming, we choose to pre-compute a set of values for e in the range [0, 100000] using the common adaptive Gauss quadrature. For limiting memory requirements, the table sampling strategy also relies on the cosine lobe integration variation (see Figure 2): the first samples are densely distributed while the last ones can be sparser; 128 linearly spaced samples are estimated in each interval [0, 1], [1, 10]. . . , [10000, 100000]. At run-time, the rendering process uses linear interpolation; it is fast and the error remains low with this sampling strategy (see section 6).

Rendering

The shading of a 3D point seen from a given viewpoint is performed using only cosine lobe sums for BRDF, visibility and lighting. Equation 3 can be rewritten:

Lo = Z !+ I # i c i ( -→ v ) J # j c j ( -→ v ) K # k c k ( -→ v ) c dot ( -→ v )d -→ v (22) where c dot ( -→ v ) is the cosine lobe -→ n • -→ v .

Using approximation 16, equation 22 becomes

Lo = I•J•K # m Z !+ cm( -→ v )d -→ v (23)
Many cosine lobe products only negligibly contribute to the outgoing radiance Lo, for two main reasons: (i) when a product c 1 • c 2 is approximated as a new cosine lobe cr (see Equation 16), the exponent er is greater than e 1 and e 2 (see Figure 6); (ii) integration values of cosine lobes decrease when exponents increase. Thus, a negligible cosine lobe multiplied with another lobe is approximated by a negligible lobe.

When the product of two cosine lobe sums is expanded (see Equations 11 and 13), we propose to check the integration value of each approximated cosine lobe product. If the integration value is smaller than a given threshold t i , the new Technical Report, XLim laboratory Figure 4 presents a statistical distribution of lobe products integration values (10 million of randomly generated products). This Figure shows that even a small threshold avoids many lobes computations.

In our application, a single threshold t i is set manually by the user for all shaded points (see Figure 5); other strategies with varying thresholds depending on time constraints or the scene configuration could also be investigated.

In addition, given a product of several lobes, when the first one can be neglected, all the remaining computations can obviously be avoided. Consequently, the ordering of computations affects the number of products actually performed. Section 5.4 describes the ordering strategy we propose for our application.

Cosine Lobes for Direct Lighting

This section describes an application that uses our representation for direct lighting. We discuss cosine lobes construction for each term of the integral. Once lobes are defined, the computations are performed using the precomputed products and integration described above.

Review of Cosine Lobe-Based BRDF Models

Phong and Lafortune models [Pho75, LFTG97] rely on cosine lobes for compactly representing BRDF. This section expresses those models according to our formalism.

Phong Model

fr( -→ v ) = k d ( -→ n • -→ v ) + ks max( -→ r • -→ v , 0) es (24) = c d ( -→ v ) + cs( -→ v ) ( 25
)
t i = 5 × 10 -2 t i = 1 × 10 -3 t i = 0 Figure 5:
The scene is rendered using several values of t i . Visibility is represented with a high number of thin cosinus lobes (small contributions) in penumbra, while umbra require a small number of thick lobes. Consequently, when t i is too high, shadows become aliased. where

c d ( -→ v ) = k d max( -→ n • -→ v , 0) 1 (26) and cs( -→ v ) = ks max( -→ r • -→ v , 0) es (27)
k d defines the diffuse cosine lobe and -→ n is the surface normal. es is the specular cosine lobe exponent, and ks the scale factor. -→ r is the reflection of -→ vo by -→ n .

Note that the original Phong model is not physically correct, but we can as well express the modified Phong model [START_REF] Lewis R | Making shaders more physically plausible[END_REF] with our representation.

Lafortune Model

fr( -→ v ) = # i & i max( -→ v T M i -→ vo ) ei (28) = # i c i ( -→ v ) (29)
where, according to [START_REF] Lafortune | Non-linear approximation of reflectance functions[END_REF], and using Helmoltz reciprocity :

c i ( -→ v ) = k i M i -→ vo ei max( M i -→ vo M i -→ vo • -→ v , 0) ei ( 30
)
k i is the scale factor of c i , M i represents the lobe axis in a local frame depending on material anisotropy and e i the exponent.

Incoming Radiance

Given a spherical light source (r, p, v) defined by a radius r, a position p and a radiance value v, the incident radiance on the hemisphere centered around a given point x is delimited by a spherical cap (a cone); the cosine lobe corresponding to the incident radiance coming from this light source is:

c light ( -→ v ) = v max p -x p -x • -→ v , 0 ecap (31)
ecap parameter corresponds the the spherical cap angle $cap (see Figure 7). ecap is estimated solving the following opti- mization problem

{ecap} = arg min {ecap} Z ! fcap( -→ v ) -ccap( -→ v ) 2 d -→ v ( 32
)
where fcap is a function representing a spherical cap:

fcap( -→ v ) = 1 if -→ v • -→ a ≥ cos($cap) 0 else (33) Technical Report, XLim laboratory and ccap( -→ v ) = 1 max( -→ a • -→ v , 0) ecap (34)
In practice, ecap is precomputed and stored in a table containing 1000 values, and linearly interpolated at run-time.

When using several light sources, the corresponding cosine lobes may overlap. However, our visibility processing (described in the next Section) naturally handles this case.

Visibility

In the off-line process, objects are simplified by sphere sets using the methods proposed by Wang et al. [WZS * 06] or by Bradshaw and O'Sullivan [START_REF] Bradshaw | Spheretree construction using dynamic medial axis approximation[END_REF]; they have been already successfully used for visibility computations by many authors [RWS * 06, SGNS07, GJW08]. These methods closely approximate the object surface with only a few spheres (for instance only 64 spheres for the Stanford bunny provide satisfactory shadow results in our application). The computed sphere sets strongly simplify objects geometry and are still deformable with articulated skeletons systems.

As for light sources, we aim at defining cosine lobes using spherical caps (cones):

V ( -→ v ) = 1 -focc( -→ v ) (35)
where focc( -→ v ) is (see Equation 33 and Figure 7):

focc( -→ v ) = 1 if -→ v • --→ aocc ≥ cos($occ) 0 else (36)
With N spheres, visibility becomes:

V ( -→ v ) = N ' j=1 1 -f j ( -→ v ) (37) = 1 - 2 N -1 # k=1 g k ( -→ v ) ( 38
)
where g k is a term corresponding to a product of spherical caps.

Unfortunately, this expression leads to the computation of 2 N terms. Replacing each spherical by only one lobe still requires 2 N lobes and introduces imprecisions that cumulate during the computation of the above product. The main reason is that one cosine lobe does not accurately approximate a spherical cap, and finally caps overlapping (that is correctly handled for visibility in Equation 38) produces errors with this representation.

To tackle these problems, we propose to hierarchically build a visibility mask, to control lobes generation (including overlapping problems) and to limit their number. One visibility mask is built for each (spherical) light source so as to precisely identify the incoming light directions and the effective blockers. This process is repeated for each point at each frame.

Given a point x, for each light source represented by a sphere (r, p), a unit square S unit surrounding the corresponding cone (see Figure 8(a)) is placed at a distance dist, deduced from the cone angle ($cap). Blockers are clipped and the remaining spheres are splatted onto S unit . Note that S unit may be placed behind the light source. Splatted spheres are used to build a quad-tree (see Figure 9). The quad-tree depth is practically limited to 4 or 5. We have chosen an aggressive approach where partially covered leaves are ignored so as to reduce the number of cosine lobes. However, the conservative alternative can be easily implemented, keeping the partially covered leaves (see Figure 10). A cosine lobe is built for each shadow leaf using a bounding sphere (see Figure 11).

Nearby cosine lobes mutually overlap and the error must be corrected. Therefore, we propose to associate a weight to each lobe:

V ( -→ v ) = c unit ( -→ v ) - N-1 # i w i c i ( -→ v ) (39)
where and w i depends on the distance dist and on the leaves depth.

c unit ( -→ v ) = 1. max( -→ n • -→ v , 0) 0 = 1 (40) c unit ( -→ v )ca( -→ v ) = ca( -→ v ) (41)
When dist varies, the cosine lobes overlap varies as well (see Figure 12). w i values are fixed experimentally for various values of dist ∈ {0.25, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} and for each depth. Weights can be adjusted manually once with a single scene (one plane, one object and one light source is enough) so as to obtain smooth shadows. First, weight values are initialized to 0; the first level weight is then chosen to obtain a perfect dark in the umbra area. The weight associated with each successive level is then fixed so as to avoid shadows sharpness with the previous ones. This process is performed only once and the resulting weights can be used with any scene. At run-time, weights are linearly interpolated according to dist.

Ordering Computations

This section describes our strategy for ordering lobe products in Equation 22, so as to reduce computing time. Let us recall that the cosine lobe cr resulting from the product of two lobes c 1 and c 2 depends on e 1 , e 2 and $ (see Section 4.1). Firstly, when $ increases, cr integration value decreases. Secondly, with visibility masks, all the products between visibility lobes and incoming radiance lobes are actually useful (cannot be neglected). Thirdly, in our application (direct lighting), incoming radiance requires less cosine lobes than visibility. Thus, firstly multiplying reflectance and lighting (before visibility) produces fast results. This strategy produces results remaining close to Monte-Carlo integration with interactive rendering (Figure 13), though it does not ensure the best precision.

Results and Discussion

This section presents some results obtained with our application, implemented in C++. All the results have been produced with an Intel Core 2 Duo T7500 (2.2GHz) CPU and an NVidia Geforce 8400 GS. Constructing our 3 pre-computed tables (product, integration and spherical caps) requires approximately 4 hours. Let us recall that these tables do not depend on the scene characteristics: they are computed once and for all. The memory required for their storage is less than 10 MB. For the sake of clarity, the following results only corresponds to cosine lobes representation, without any optimization: for each frame, the whole image is completely computed (no caching technique is employed, nor spatial or temporal coherency). Our application consists in shading a selected set of 3D points (either triangle vertices or estimated per pixel, employing a G-buffer). Shading is performed on the CPU and rasterization uses OpenGL.

Due to the flexibility and adaptivity of cosine lobes representation (number of lobes, bandwith, free distribution, etc.), computing efforts can be concentrated on regions where lighting needs to be detailed. For instance, penumbra areas requires a higher number of lobes for representing visibility than umbra areas or fully-lit regions, where only a few computations can produce precise results. Consequently, our application performances vary according to chosen materials (Figure 14), number of shaded points, maximum depth of visibility masks and chosen threshold of cosine lobe products (Figures 16 and5) or relative position between occluders and lights. Frame rates are given in the caption of each figure.

Presently, many approaches dedicated to dynamic scenes are based on SH. Generally speaking such orthogonal basis functions are adapted to integrations. However, the projection on these basis functions is costly. For instance, on Figure 16, umbra and lit regions would have required the same amount of computation and are restricted to low-frequency effects, contrary to our cosine-lobe based approach.

Once cosine lobe sums are defined, the approach given in section 4 approximates lobe products and uses tables for reducing computing time. Figure 15 shows that these approximations lead to a small difference, compared to Monte-Carlo integration. For the Monte-Carlo reference image, the maximum radiance value is 0.701961, the mean radiance value is 0.127906; with our approach, the maximum difference is 0.054902 while the mean difference is 0.003884. Note that our method introduces a bias: regions located in shadows, where more cosine lobes products are computed, exhibit more negative difference values.

We also compare the images produced by our application with path tracing for direct lighting and soft shadows. As shown in Figure 13, the results are hardly distinguishable. sue for realistic rendering; their processing is often handled specifically. Figures 13 and17 show that our method naturally produces realistic direct lighting (including hard and soft shadows), without introducing specific processing. Our application can be used with any type of material, using a spatially varying Lafortune BRDF model (see Figure 14). We do not introduce any occluder in this scene for highlighting the low cost for representing reflectance and lighting with cosine lobe sums. Figure 18 illustrates a difficult case for the construction of the visibility mask. Self-shadowing is robustly supported with adapted sphere sets.

Conclusion

In this paper, we present a new method for solving the rendering equation using approximated cosine lobe products.

We propose an interactive application for direct lighting including all-frequency lighting / materials, soft shadows and dynamic scenes. As shown in the results, cosine lobes as basis functions are useful for properly targeting relevant hemispherical directions (flexibility) in terms of BRDF, incoming radiance or visibility. In addition, their representation only requires a few parameters (compactness).

Shading is only performed on the CPU and performances can be further improved using a full GPU implementation, adapting for instance the streaming reduction algorithm proposed in [START_REF] Assarsson | Efficient stream reduction on the gpu[END_REF] for avoiding useless computations concerning negligible lobe products. Furthermore, our application demonstrates the feasibility of rapidly computing direct lighting, without using any spatial and temporal coherency. Visibility masks are often identical for close neighborhoods of shaded points; this coherency can be used for reducing the computing time.

Global illumination in dynamic scenes has been addressed using spherical proxies and spherical harmonics [START_REF] Govindaraju | Image-based proxy accumulation for real-time soft global illumination[END_REF][START_REF] Guerrero | Real-time indirect illumination and soft shadows in dynamic scenes using spherical lights[END_REF]. Similarly to these authors, our direct lighting application makes use of blocker accumulation. In the future, Technical Report, XLim laboratory Figure 18: Two chocolate-flavoured bunnies in a green box. These pictures take 6.25s for rendering. 480000 points are shaded, using 1024 spherical occluders, 2 cosine lobes for the reflectance models and 1 spherical light source. We choose 5 for the maximum depth of the visibility mask.

we aim at adapting their interesting simplifications to cosine lobes for overcoming the usual SH drawbacks.

Finally, cosine lobes can also be used in various SRBF applications, such as the representation of high frequency environment maps [START_REF] Tsai Y.-T | All-frequency precomputed radiance transfer using spherical radial basis functions and clustered tensor approximation[END_REF], anisotropic materials [START_REF] Lafortune | Non-linear approximation of reflectance functions[END_REF], or real-time environment map rendering applications, such as Green et al. [START_REF] Green | Efficient reflectance and visibility approximations for environment map rendering[END_REF] method.

Figure 1 :

 1 Figure 1: Example of cosine lobe product. Cosine lobes in the left box are multiplied. The top-right box provides the resulting surface while the bottom-right box illustrates the (single) lobe corresponding to our approximation: the global shape and direction are preserved.

Figure 2 :

 2 Figure 2: Integration value according to cosine lobe exponents (logarithmic scale), with a scaling factor s set to 1. The low variation for high exponents allows a sparser sampling than with lower exponents.

Figure 3 :

 3 Figure3: Lobe product approximation (Section 4.1) error: e 1 and e 2 correspond to the lobes exponents and $ is the angle between the lobes axis. The error is defined as the absolute difference between the integration of c 1 • c 2 and cr. The error value is represented by a color; the white regions correspond to a low error; the black regions to a high error. An important error (≈ 50%) is noticeable for e 1 and e 2 between 0 and 1 approximatively and an angle $ around %/2. Left: full representation for e 1 and e 2 in range [0..100000]; Right: zoom on the region corresponding to the highest error.

Figure 4 :

 4 Figure 4: Percentage of lobe products integration (ordinate) smaller than or equal to a given value. The dashed curve represents integrated products in a low frequency environment (exponents ∈ [0, 10]); the dotted curve corresponds to a higher frequency environment (exponents ∈ [0, 100]) and the last curve to exponents ∈ [0, 1000]. In practice, exponents often grow up to 15000.

Figure 6 :

 6 Figure6: Representation of er according to e 1 and e 2 for the approximation of the product c 1 .c 2 for a fixed $. Note that er is greater or equal than e 1 and e 2 (the corresponding surface is not a plane).

Figure 7 :

 7 Figure 7: $cap represents the spherical cap angle corresponding to sphere (p, r).

Figure 8 :Figure 9 :

 89 Figure 8: Occluding spheres splatting: (a) a unit square is placed according to the spherical light source; (b) useless blockers are clipped; (c) remaining spheres are splatted onto the unit square.

Figure 10 :

 10 Figure 10: Shadows artifacts when depth of visibility mask is not high enough. Top row: with aggressive visibility; bottom row: with conservative visibility. From left to right: visibility depth is equal to 2, 3 and 5.

Figure 11 :

 11 Figure 11: Once visibility mask is built (left), one bounding sphere is constructed per shadow leaf and the associated cosine lobes are produced (right).

Figure 12 :

 12 Figure 12: The cosine lobes overlapping vary according to the distance dist, even with the same visibility mask.

Figure 13 :

 13 Figure 13: Shadows with direct lighting for a scene composed of 5 and a plane: (a) reference path tracing image and (b) image computed using our application. Maximal depth of the visibility mask is set to 5.

Figure 14 :

 14 Figure 14: Spatially varying BRDF are possible using Lafortune model. This scene is composed of 65536 shaded points, 1 spherical light source and no occluder, 13.73 frame per second.

Figure 15 :

 15 Figure 15: Comparison between Monte-Carlo integration (top left) and our approach (top right) for the cosine lobe sums integration. Lobe sums are generated using the methods described in Section 5. The bottom images show the difference. The bottom left image represents positive values and the bottom right the negatives.

Figure 17 :

 17 Figure 17: Shadow smoothing according to the distance between plane and the Bunny. Frame rates decrease with respect to the (increasing) shadowed area. Penumbra implies complex visibility masks and more cosine lobes. 16384 points are shaded considering a strictly diffuse material, 1 spherical light source, 64 spherical occluders, fixing the maximal depth of the visibility mask to 5 and the threshold to 0.001 (Section 4.3).

Table 6

 6 presents the computing time ratio corresponding to each rendering part, for different types of scenes, with hard and soft shadows, complex or simple material and various numbers of occluders. The last row corresponds to a scene that does not contain any occluder, implying no visibility computations.The main part of computing time is dedicated to visibility mask construction, as for any other blocker accumulation method.

	Shadows are considered by many authors as a key is-

Figure 16: Performance gains obtained according to the threshold t i (Section 4.3). The scene is composed of 4096 shaded point, one spherical light source and 64 spherical occluders (the bunny). The threshold is chosen so as to keep indistinguishable the difference between resulting images (left and right columns). The false-color images illustrate the number of cosine lobe products computed (scale is given at the bottom).

The specular highlight properly interacts with the shadows, our method naturally allows to avoid many computations due to the threshold. Note that frame rates increase when the highlight moves away from the shadow. We can see a reduction of about 50% of the number of product computations (corresponding to 2 frames per second in practice, because of visibility mask construction cost).

Technical Report, XLim laboratory


Memb, IEEE Blaise Ravelo
email: blaise.ravelo@yahoo.fr

Member, IEEE Mathieu Guerin
email: mathieu.guerin@im2np.fr

Member, IEEE Wenceslas Rahajandrai

Membe, IEEE Valentin Gies
email: valentin.gies@im2np.fr

Member, IEEE Lala Rajaoarisoa
email: lala.rajaoarisoa@imt-lille-douai.fr

Member, IEEE Sébastien Lalléchère
email: sebastien.lallechere@uca.fr

Wenceslas Rahajandraibe
email: wenceslas.rahajandraibe@im2np.fr

ON INDUSTRIAL ELECTRONICS 1 Low-Pass NGD Numerical Function and STM32 MCU Emulation Test 1 2

Keywords: Circuit design, emulation, microcontroller 25 unit

This article introduces an original microcon-6 troller unit (MCU) design of numerical low-pass (LP) neg-7 ative group delay (NGD) function. The innovative theory 8 of the numerical LP-NGD function is developed based on 9 the first-order analog transfer function discretization. The 10 infinite impulse response (IIR) LP-NGD is fundamentally 11 formulated in function of the desired NGD value, cutoff 12 frequency, gain, and the MCU sampling frequency. A STM32 13 MCU proof-of-concept (POC) is tested to implement the 14 IIR LP-NGD function. Different real-time tests with visual-15 ization of input and output analog signals from the MCU 16 LP-NGD POC were performed. As expected, time-advance 17 demonstration tests with milli-second short-and several 18 hour long-duration time-scale with arbitrary waveform sig-19 nals from temperature and humidity sensors. The signal 20 time-advance is not in contradiction with the causality. The 21 proposed digital MCU function opens a potential future in-22 dustrial application of LP-NGD function via sensored signal Q1 23 anticipation.

system. The quality of service (QoS) of remote sensing systems forecasted Industry 4.0 factories to operate correctly depends undeniably on the delay effects. The networked modern industrial sensing systems present an open risk of anomaly detection because of unintentional delays. So far, no rigorous solution is available to overcome the time delay effects in almost all areas of engineering [2]- [9].

One can cite different challenges to be overcome to face up the delay impacts. The signal delay degrades the digital communication QoS due to the wireless propagation channel dispersion and data desynchronization [3]- [4]. In the area of automation and control engineering, different industrial application examples as optimal vehicle activities [5], networked control system, and power system damping control [6] present performance limitations because of the delay effect. Time-domain analyses confirm that the stability of dynamic systems depend unavoidably to the feedback delays [7].

To take into account the delay effects on the system performance, analytical, and modeling methods for the enhancement of printed circuit board (PCB) signal integrity were proposed [8]- [11]. Delay estimation techniques were developed for estimating communication system performance [3]. Diverse estimation techniques based on complex analytical signals [12] have been developed. In addition to the noise effect issues, the group delay (GD) limits undesirably the communication performances [13]- [14].

In the area of biomedical engineering, under the telemedicine revolution, new generation of wearable biosensors is currently developed [15]. The delay effects play also a vital role on the performance of the wearable sensor signal monitoring. An innovative photoplethysmography signal measurements of respiration, oxygen, and blood pressure rate measurement were proposed [16]- [18]. A technique allowing to use an optimal number and positioning of multiple accelerometers for respiratory waveform estimation is developed in [19]. To improve the medical sensor performances, an advanced signal detection technique using negative GD (NGD) function is developed in [20]. Under the same challenge, an NGD predictor based on anticipatory relax dynamics using basically a delayed induced feedback function is proposed in [21] and [22].

In electronics engineering, NGD microwave monolithic integrated circuits have been employed to equalize the undesirable signal degradations as regularly found in the PCBs [23]. This NGD equalization technique was initiated in [24] for RC-effect The technique was introduced in [25] to target industrial applications by overcoming multisensor time-delays. Acting as an unfamiliar function, it is worth to define what is the NGD function meaning. To answer to this curiosity, a brief state of the art about the NGD function will be described in the following paragraph.

I E E E P r o o f

The NGD phenomenon was initially validated with lowfrequency (LF) electronic circuits in 1990s [26], [27]. Despite the counterintuitive aspect, the NGD circuits can be designed and implemented like classical and familiar electrical circuits as filters, amplifiers, phase shifters, oscillators. The LF NGD topologies proposed initially in [26]- [29] were built with classical electronic passive (R, L, and C) and operational amplifier components. The NGD phenomenon meaning was illustrated with time-domain experimentations showing the possibility to propagate the NGD output signal in time-advance with the input one [26]- [29]. This demonstration was realized with PCB showing an NGD output lamp lightened on before the input one [29]. This counterintuitive NGD phenomenon is not in contradiction with the causality principle [26]- [29]. For further understanding about the NGD concept, an investigation on the similitude between the NGD and linear filter gain behaviors is introduced in [30]. The concept of low-pass (LP) [31], [32] NGD function was initiated. Despite the research work progress in the area of electronics engineering, the NGD function remains an unmatured concept for most of industrial engineers. Open questions are regularly asked about the industrial potential applications of the NGD function.

For this reason, a completely original research work on the design of numerical LP NGD circuit is developed in the present paper. The LP-NGD function is implemented on microcontroller unit (MCU) operating in different time scales. The NGD function can be potentially used as a pioneer work of anticipator for the real-time predetection of actuator and sensor faults and intelligent system diagnoses [33], [34]. The research work originality relies on the MCU numerical LP-NGD function design and test. For the better understanding, this paper is organized in three main sections. Section II will describe the numerical LP-NGD transfer function (TF) theory. Section III will present the LP-NGD MCU demonstrator as proof-of-concept (POC).

Finally, Section IV concludes this article.

129

A. Input Signal and MCU Specifications

130

A LP-NGD function is dedicated to operate with base band 131 signal, which can be noticed by v in . The ideal specifications of 132 its frequency spectrum magnitude in dB of V in (jω) are illustrated 133 in Fig. 2 by denoting the angular frequency variable ω. In this 134 diagram, the signal spectrum magnitude is assumed V in = 0 in 135 the forbidden zone represented by the shadowed area. Signal 136 v in can present any arbitrary waveform but should be limited 137 bandwidth denoted ω max . The signal minimal and maximal 138 strengths, denoted by V min and V max , need to be defined.

139

The MCU architecture is assumed powered by V dd = V 0 . It is 140 mainly constituted by the block of computing core interfaced by 141 analog digital converters (ADC) and DAC as input and output 142 interfaces, respectively. The MCU operation is cadenced by 143 its internal clock generating the sampling clock supposed with 144 period denoted T s = 1/f s (f s is the sampling frequency).

v in (k) = v in (kT s ) v out (k) = v out (kT s ) . ( 1
)

B. Theoretical Investigation of Discrete

151

LP-NGD Function 152

This section describes the main specifications and character-153 ization of numerical LP-NGD function. 1) The NGD value at very low frequencies

τ n = τ (ω ≈ 0) < 0. (4)
2) The NGD cutoff angular frequency ω n = 2πf n , which is a root of the equation τ (ω) = 0.

Fig. 3(c) highlights the LP-NGD magnitude response specifications. For the case of first order linear NGD system (which will be explored in the following paragraph), it associates three main characteristics.

3) The minimal magnitude N dBmin = N 0dB = N dB (ω ≈ 0), which should correspond to the magnitude at very LFs.

4) The maximal magnitude

N dB max = N dB (ω ≈ ∞),
which should be reached at very high frequencies.

5) And the magnitude

N ndB = N dB (ω = ω n ) at cutoff frequency.

2) Canonical Form of First Order LP-NGD TF:

With the real parameters defined previously N n , τ n , and ω n , the canonical form of LP-NGD analog TF is defined as [31] N

(s) = N n (1 + τ a s)/(1 + τ b s) (5)
where

τ a = τ 2 n + 4/ω 2 n -τ n /2 (6)
τ b = τ 2 n + 4/ω 2 n + τ n /2. ( 7
)
The integrity between NGD input and output signals, v in and v out can be quantified with the quantity M n = N (ω n ). The LP-NGD bandwidth magnitude flatness and maximal flatness factor can be assessed by the ratio, respectively, as follows:

ξ( ) = M n ( )/N n ( 8
)
ξ max ( ) = N max ( )/N n ( 9
)
where the NGD value-bandwidth product

= -ω n τ n . ( 10
)
By considering the LP-NGD TF expressed in (5), the previous magnitude flatness factors are rewritten as ξ( )

ξ max ( ) = ( 2 + 4 + )/( 2 + 4 -). ( 12
)
3) Characterization of LP-NGD TF Flatness: The magni-190 tude flatness constitutes one of the main parameters affecting 191 the LP-NGD function performance. It is particularly important 192 to investigate the LP-NGD TF magnitude variation range in the 193 NGD frequency band. Following this point, we considered the 194 general plot of the flatness defined by (11)

v out (t) -N n v in (t) = N n τ a dv in (t) dt -τ b dv out (t) dt . ( 14
)
The NGD function operation can be rigorously realized under

214 v out (0) = v in (1) = v in (0) = 0. (15)
The

Δv out (kT s )/Δt = [v out (k + 1) -v out (k)] /T s . ( 17
)
Substituting previous (16) and ( 17) into (14), the numerical 219 LP-NGD difference equation can be fundamentally formulated

220 by 221 v out (k + 1) = ζ 2 v in (k + 1) + ζ 1 v in (k) + ζ 0 v out (k) (18)
where

222 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ζ 0 = √ 4+ω 2 n τ 2 n +ω n τ n √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) ζ 1 = N n (ω n τ n - √ 4+ω 2 n τ 2 n ) √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) ζ 2 = N n √ 4+ω 2 n τ 2 n +ω n (2T s -τ n ) √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) . ( 19
)
The LP-NGD infinite impulse response (IIR) coefficient vari-

223

ations are semi-logarithmically plotted in Fig. 5

k 0 = Ent [τ n /T s ] (20)
Fig. 6. STM32 MCU integrated circuit top view [37].

with Ent[x] represents the superior integer part of real x.

As a validation, the innovative LP-NGD numerical function, original experimental results with STM32 MCU will be explored in the following section.

III. LP-NGD NUMERICAL FUNCTION EXPERIMENTATION WITH STM32 MCU

To verify the established numerical LP-NGD theory feasibility, a POC is designed via C-program code emulation driving an STM32 MCU board. This section presents the first NGD numerical experimental results based on the calculation and measurement comparison.

A. Description of the MCU-Based LP-NGD IIR Function

A Nucleo L476RG development board using an STM32L476RG MCU was used for this article [37]. The MCU package is 64-pin integrated circuit. Fig. 6 presents the top view of the employed STM32 MCU manufactured by STMicroelectronics [37] and the active pins needed to implement the LP-NGD TF. The MCU board is fed by V 0 = 5 V dc power from USB connector. The microprocessor operating frequency is fixed at 80 MHz ADC are used in 12-bit precision.

The dynamic of operating input analog voltage can vary in the range between V min = 0 V to V max = 3.3 V with quantum resolution equal to ΔV = 800 μV under sampling frequency, f s = 5.33 MHz. Then, the IIR NGD TF described in (18) is coded in C-language program, as indicated in Fig. 7.

During program execution, analog input voltage applied on pin PA0 is converted in a digital signal and processed. Output analog signal is delivered from pin PA5 through the DAC. During the operation test, the MCU configuration including the input and output signal interactions with the microprocessor is carried out by using the STM32CubeIDE software installed on the driver PC as described by the experimental test of the following section. During the test, the oscilloscope input impedance is configured as R in = 1 GΩ. The employed test equipment specifications are addressed in Table I. Then, the tested signals are recorded in csv format and replotted in MATLAB environment. It is noteworthy that the MCU latency is about 21 μs.

B. Numerical LP-NGD Experimental Setup

C. Demonstration of Time<0 Signal Prediction

This section discusses the demonstration results with three types of sensors. Prediction of arbitrary waveform signals in different time-scales is presented. The Gaussian pulse and sinc wave pulse signals were chosen for the LP-NGD time-domain validation because of the following reason. 1) They present a limited bandwidth, which can be assumed as the LP-NGD circuit cutoff frequency.

2) They enable to calculate the pulse duration in function of the LP-NGD circuit cutoff frequency.

3) To show the LP-NGD signature with the time-advance of leading-and tailing-edge between the input and output pulses.

4) To demonstrate the signal integrity by comparing the pulsewidth and amplitudes.

5) And to measure the magnitude effect by comparing the maximum peaks of input and output pulse signals.

1) LP-NGD Function Frequency

316 v in (t) = V max e -(t-t 0 ) 2 /τ 2 a ( 21
)
with maximal peak V max = v in (t 0 ), at instant time t 0 , and 317 spectrum bandwidth depending to the time parameter [38] 318

τ a = a dB ln(10)/5/(πf a ) (22)
where a dB is the Gaussian spectrum attenuation at frequency

τ f = t out f - 322 t in f , with v in (t in r,f < t 0 ) = max[v in (t)]/2 and v out (t out r,f < t 0 ) = 323 max[v out (t)]/2
, combined to voltage gain between the input and 324 output signals

325 a = max [v out (t)] / max [v in (t)] . ( 23
)
And also, the input and output relative cross-correlation can 326 be written as

327 (v in , v out ) = k max k=1 [v in (k) -vin ] [v out (k) -vout ] k max k=1 [v in (k) -vin ] 2 k max k=1 [v out (k) -vout ] 2 (24)
where the input and output average values, vin =

328 k max k=1 v in (k)/k max and vout = k max k=1 v out (k)/k max .

329

The transient calculation was carried out with Gaussian signal 330 with t max = 20 ms and different parameters indicated by Table II. Moreover, Table II summarizes the LP-NGD transient characterization parameters. It is noteworthy that the numerical LP-NGD MCU input and output signal cross correlations are better than 92%. The generality of the sensational time-advance was validated with two different arbitrary waveform and time-scale signals as explored in the following paragraphs.

3) Numerical LP-NGD Sinc Waveform Signal Response:

Another transient test of the LP-NGD numerical circuit was performed with 10-ms pulse duration and 1-V peak-to-peak amplitude sinc waveform. The considered sinc signal presents a frequency spectrum belonging into the NGD frequency band. Fig. 11(a) displays the comparison between the calculated v out C and measured v out M . Once again, the output signal time-advance was observed, as seen in the screenshot of Fig. 11(b). During the calculation, the time step was fixed to 0.1 ms. Once again, based on this result showing about -1 ms time advance average of output sinc leading and trailing edges, the STM32 MCU is able to operate as a LP-NGD numerical circuit. In this case, based on formula (24), measured and calculated output signals have a relative cross-correlation of about 89.3% and 87.7%. In addition to the signal generator input signal test, computational investigation with more industrial context of sensored arbitrary waveform signals plotted in several hours and very long-time duration is examined in the following paragraph. originally designed and emulated on STM32 MCU board [37].

I E E E P r o o f

TABLE III TEMPERATURE AND HUMIDITY WIRELESS SENSOR CHARACTERISTICS

4) Hour Time-Scale Arbitrary Sensored Signal Predic-

QUERIES

Q1. Author: Please provide the expansion for the acronym DAC. Q2. Author: Please provide the ORCID for the corresponding author Mathieu Guerin. Q3. Author: Please update the Ref. [37] Q4. Author: Please provide the educational details (degree/university/location/subject) for author Blaise Ravelo. Q5. Author: Please check whether the biography of author Mathieu Guerin is correct as set. Q6. Author: Please provide the subject in which the author Mathieu Guerin received the doctorate degree. Q7. Author: Please provide the year in which the author Wenceslas Rahajandraibe received the Ph.D. degree.

I E E E P r o o f

Low-Pass NGD Numerical Function and STM32 MCU Emulation Test Blaise Ravelo is with the Nanjing University of Information Science and Technology, Nanjing 210044, China (e-mail: blaise.ravelo@yahoo.fr).

Mathieu Guerin, Wenceslas Rahajandraibe, and Valentin Gies are with the CNRS, University of Toulon, 83130 Marseille, France, and also with the Aix-Marseille University, 3007 Marseille, France (e-mail: mathieu.guerin@im2np.fr; wenceslas.rahajandraibe@im2np.fr; valentin.gies@im2np.fr).

Lala Rajaoarisoa is with the IMT Lille Douai, Institut Mines-Télécom, Centre for Digital Systems, University Lille, F-59000 Lille, France (e-mail: lala.rajaoarisoa@imt-lille-douai.fr).

Sébastien Lalléchère is with the Institute Pascal, SIGMA Clermont, Université Clermont Auvergne, 63001 Clermont Ferrand, France (e-mail: sebastien.lallechere@uca.fr).

Color as optimal vehicle activities [5], networked control system, and power system damping control [6] present performance limitations because of the delay effect. Time-domain analyses confirm that the stability of dynamic systems depend unavoidably to the feedback delays [7].

To take into account the delay effects on the system performance, analytical, and modeling methods for the enhancement of printed circuit board (PCB) signal integrity were proposed [8]- [11]. Delay estimation techniques were developed for estimating communication system performance [3]. Diverse estimation techniques based on complex analytical signals [12] have been developed. In addition to the noise effect issues, the group delay (GD) limits undesirably the communication performances [13]- [14].

In the area of biomedical engineering, under the telemedicine revolution, new generation of wearable biosensors is currently developed [15]. The delay effects play also a vital role on the performance of the wearable sensor signal monitoring. An innovative photoplethysmography signal measurements of respiration, oxygen, and blood pressure rate measurement were proposed [16]- [18]. A technique allowing to use an optimal number and positioning of multiple accelerometers for respiratory waveform estimation is developed in [19]. To improve the medical sensor performances, an advanced signal detection technique using negative GD (NGD) function is developed in [20]. Under the same challenge, an NGD predictor based on anticipatory relax dynamics using basically a delayed induced feedback function is proposed in [21] and [22].

In electronics engineering, NGD microwave monolithic integrated circuits have been employed to equalize the undesirable signal degradations as regularly found in the PCBs [23]. This NGD equalization technique was initiated in [24] for RC-effect The NGD phenomenon was initially validated with lowfrequency (LF) electronic circuits in 1990s [26], [27]. Despite the counterintuitive aspect, the NGD circuits can be designed for the real-time predetection of actuator and sensor faults and intelligent system diagnoses [33], [34]. The research work originality relies on the MCU numerical LP-NGD function design and test. For the better understanding, this paper is organized in three main sections. Section II will describe the numerical LP-NGD transfer function (TF) theory. Section III will present the LP-NGD MCU demonstrator as proof-of-concept (POC).

Finally, Section IV concludes this article.

A. Input Signal and MCU Specifications

130

A LP-NGD function is dedicated to operate with base band 131 signal, which can be noticed by v in . The ideal specifications of 132 its frequency spectrum magnitude in dB of V in (jω) are illustrated 133 in Fig. 2 by denoting the angular frequency variable ω. In this 134 diagram, the signal spectrum magnitude is assumed V in = 0 in 135 the forbidden zone represented by the shadowed area. Signal 136 v in can present any arbitrary waveform but should be limited 137 bandwidth denoted ω max . The signal minimal and maximal 138 strengths, denoted by V min and V max , need to be defined.

139

The MCU architecture is assumed powered by V dd = V 0 . It is 140 mainly constituted by the block of computing core interfaced by 141 analog digital converters (ADC) and DAC as input and output 142 interfaces, respectively. The MCU operation is cadenced by 143 its internal clock generating the sampling clock supposed with 144 period denoted T s = 1/f s (f s is the sampling frequency). Given 145 the case of input signal duration t max , the MCU must operate 146 with time step following the integer index k = 01,2, …, k max = 147 Ent[t max /T s ] by denoting Ent[x] the superior integer part of real 148 positive x. The discrete input and output signals at the instant 149 time t = kT s are represented by 1) The NGD value at very low frequencies

150 v in (k) = v in (kT s ) v out (k) = v out (kT s ) . ( 1
)

B. Theoretical Investigation of Discrete

τ n = τ (ω ≈ 0) < 0. (4)
2) The NGD cutoff angular frequency ω n = 2πf n , which is a root of the equation τ (ω) = 0.

Fig. 3(c) highlights the LP-NGD magnitude response specifications. For the case of first order linear NGD system (which will be explored in the following paragraph), it associates three main characteristics.

3) The minimal magnitude N dBmin = N 0dB = N dB (ω ≈ 0), which should correspond to the magnitude at very LFs.

4) The maximal magnitude

N dB max = N dB (ω ≈ ∞),
which should be reached at very high frequencies.

5) And the magnitude N ndB = N dB (ω = ω n ) at cutoff frequency.

2) Canonical Form of First Order LP-NGD TF:

With the real parameters defined previously N n , τ n , and ω n , the canonical form of LP-NGD analog TF is defined as [31] N

(s) = N n (1 + τ a s)/(1 + τ b s) (5)
where

τ a = τ 2 n + 4/ω 2 n -τ n /2 (6)
τ b = τ 2 n + 4/ω 2 n + τ n /2. ( 7
)
The integrity between NGD input and output signals, v in and v out can be quantified with the quantity M n = N (ω n ). The LP-NGD bandwidth magnitude flatness and maximal flatness factor can be assessed by the ratio, respectively, as follows:

ξ( ) = M n ( )/N n ( 8
)
ξ max ( ) = N max ( )/N n ( 9
)
where the NGD value-bandwidth product

= -ω n τ n . ( 10
)
By considering the LP-NGD TF expressed in (5), the previous magnitude flatness factors are rewritten as ξ( )

= 4 + ( 2 + 4 + ) 2 / 4 + ( 2 + 4 -) 2 (11)
ξ max ( ) = ( 2 + 4 + )/( 2 + 4 -). (12)
3) Characterization of LP-NGD TF Flatness: The magni-190 tude flatness constitutes one of the main parameters affecting 191 the LP-NGD function performance. It is particularly important 192 to investigate the LP-NGD TF magnitude variation range in the 193 NGD frequency band. Following this point, we considered the 194 general plot of the flatness defined by (11)

v out (t) -N n v in (t) = N n τ a dv in (t) dt -τ b dv out (t) dt . ( 14
)
The NGD function operation can be rigorously realized under

v out (0) = v in (1) = v in (0) = 0. (15)
The

Δv out (kT s )/Δt = [v out (k + 1) -v out (k)] /T s . ( 17
)
Substituting previous ( 16) and ( 17) into ( 14), the numerical 219 LP-NGD difference equation can be fundamentally formulated

220 by 221 v out (k + 1) = ζ 2 v in (k + 1) + ζ 1 v in (k) + ζ 0 v out (k) (18)
where

222 ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ζ 0 = √ 4+ω 2 n τ 2 n +ω n τ n √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) ζ 1 = N n (ω n τ n - √ 4+ω 2 n τ 2 n ) √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) ζ 2 = N n √ 4+ω 2 n τ 2 n +ω n (2T s -τ n ) √ 4+ω 2 n τ 2 n +ω n (τ n +2T s ) . ( 19
)
The LP-NGD infinite impulse response (IIR) coefficient vari-

k 0 = Ent [τ n /T s ] (20)
Fig. 6. STM32 MCU integrated circuit top view [37].

with Ent[x] represents the superior integer part of real x.

As a validation, the innovative LP-NGD numerical function, original experimental results with STM32 MCU will be explored in the following section.

III. LP-NGD NUMERICAL FUNCTION EXPERIMENTATION WITH STM32 MCU

To verify the established numerical LP-NGD theory feasibility, a POC is designed via C-program code emulation driving an STM32 MCU board. This section presents the first NGD numerical experimental results based on the calculation and measurement comparison.

A. Description of the MCU-Based LP-NGD IIR Function

A Nucleo L476RG development board using an STM32L476RG MCU was used for this article [37]. The MCU package is 64-pin integrated circuit. Fig. 6 presents the top view of the employed STM32 MCU manufactured by STMicroelectronics [37] and the active pins needed to implement the LP-NGD TF. The MCU board is fed by V 0 = 5 V dc power from USB connector. The microprocessor operating frequency is fixed at 80 MHz ADC are used in 12-bit precision.

The dynamic of operating input analog voltage can vary in the range between V min = 0 V to V max = 3.3 V with quantum resolution equal to ΔV = 800 μV under sampling frequency, f s = 5.33 MHz. Then, the IIR NGD TF described in (18) is coded in C-language program, as indicated in Fig. 7.

During program execution, analog input voltage applied on pin PA0 is converted in a digital signal and processed. Output analog signal is delivered from pin PA5 through the DAC. During the operation test, the MCU configuration including the input and output signal interactions with the microprocessor is carried out by using the STM32CubeIDE software installed on the driver PC as described by the experimental test of the following section. During the test, the oscilloscope input impedance is configured as R in = 1 GΩ. The employed test equipment specifications are addressed in Table I. Then, the tested signals are recorded in csv format and replotted in MATLAB environment. It is noteworthy that the MCU latency is about 21 μs.

B. Numerical LP-NGD Experimental Setup

C. Demonstration of Time<0 Signal Prediction

This section discusses the demonstration results with three types of sensors. Prediction of arbitrary waveform signals in different time-scales is presented. The Gaussian pulse and sinc wave pulse signals were chosen for the LP-NGD time-domain validation because of the following reason. 1) They present a limited bandwidth, which can be assumed as the LP-NGD circuit cutoff frequency.

2) They enable to calculate the pulse duration in function of the LP-NGD circuit cutoff frequency.

3) To show the LP-NGD signature with the time-advance of leading-and tailing-edge between the input and output pulses.

4) To demonstrate the signal integrity by comparing the pulsewidth and amplitudes.

5) And to measure the magnitude effect by comparing the maximum peaks of input and output pulse signals.

v in (t) = V max e -(t-t 0 ) 2 /τ 2 a ( 21
)
with maximal peak V max = v in (t 0 ), at instant time t 0 , and 317 spectrum bandwidth depending to the time parameter [38] 318 τ a = a dB ln(10)/5/(πf a )

where a dB is the Gaussian spectrum attenuation at frequency

a = max [v out (t)] / max [v in (t)] . ( 23
)
And also, the input and output relative cross-correlation can 326 be written as

327 (v in , v out ) = k max k=1 [v in (k) -vin ] [v out (k) -vout ] k max k=1 [v in (k) -vin ] 2 k max k=1 [v out (k) -vout ] 2 (24)
where the input and output average values, vin =

328 k max k=1 v in (k)/k max and vout = k max k=1 v out (k)/k max . 329
The transient calculation was carried out with Gaussian signal 330 with t max = 20 ms and different parameters indicated by Table II. Moreover, Table II summarizes the LP-NGD transient characterization parameters. It is noteworthy that the numerical LP-NGD MCU input and output signal cross correlations are better than 92%. The generality of the sensational time-advance was validated with two different arbitrary waveform and time-scale signals as explored in the following paragraphs.

3) Numerical LP-NGD Sinc Waveform Signal Response:

Another transient test of the LP-NGD numerical circuit was performed with 10-ms pulse duration and 1-V peak-to-peak amplitude sinc waveform. The considered sinc signal presents a frequency spectrum belonging into the NGD frequency band. Fig. 11(a) displays the comparison between the calculated v out C and measured v out M . Once again, the output signal time-advance was observed, as seen in the screenshot of Fig. 11(b). During the calculation, the time step was fixed to 0.1 ms. Once again, based on this result showing about -1 ms time advance average of output sinc leading and trailing edges, the STM32 MCU is able to operate as a LP-NGD numerical circuit. In this case, based on formula (24), measured and calculated output signals have a relative cross-correlation of about 89.3% and 87.7%. In addition to the signal generator input signal test, computational investigation with more industrial context of sensored arbitrary waveform signals plotted in several hours and very long-time duration is examined in the following paragraph.

Fig. 1 .

 1 Fig. 1. Illustration scenario of NGD application. (a) Delay penalizing system. (b) NGD corrected system. (b) Input and output signal chronogram.

Fig. 2 .

 2 Fig. 2. Specific characterization diagram of numerical LP-NGD input signal.II. MCU DESIGN PRINCIPLE OF DISCRETE

154 1 )

 1 NGD Specifications: The NGD TF can be represented as 155 a single-input single-output block diagram, as shown in Fig.3(a). 156 With the Laplace variable s, the associated analog TF can be 157 formulated by N (s) = V out (s)/V in (s). The NGD analysis in the 158 frequency domain is based on the consideration of magnitude 159 response, N (ω) = |N (jω)| and the GD defined by 160 τ (ω) = -∂ϕ(ω)/∂ω (2) with 161 ϕ(ω) = arctan {Im [N (jω)] /Re [N (jω)]} .(3)

Fig. 3 .

 3 Fig. 3. (a) Block, (b) GD, and (c) magnitude response diagrams of LP-NGD function.

Fig. 3 (

 3 Fig. 3(b) represents the GD diagram specific to the ideal LP-NGD function, which is defined by two fundamental characteristics as follows.

Fig. 4 .

 4 Fig. 4. LP-NGD TF magnitude flatness versus product ω n τ n in (a) semilog x-axis, (b) log-log plots, and (c) ratio from (9).

Fig. 4 Fig. 5 .

 45 Fig. 5. Variations of LP-NGD IIR coefficients versus product ω n T s .

  212 condition |τ n | = 4T s . The discrete TF is designed by taking into 213 account to the initial condition

Fig. 8 (

 8 Fig. 8(a) highlights the real-time experimental setup of the STM32 MCU based IIR LP-NGD function. As indicated by

Fig. 7 .

 7 Fig. 7. Numerical LP-NGD circuit code.

Fig. 8 .

 8 Fig. 8. (a) Illustrative diagram and (b) photograph of the LP-NGD MCU-based experimental setup.

Fig. 8 (Fig. 9 .

 89 Fig. 8(a), the input pulse analog voltage provided by the signal268

  Responses: Before the transient test, the LP-NGD circuit were synthesized under the desired specifications, static gain N 0 = 0 dB, cutoff frequency f n = 1 kHz, and NGD value τ n = -1 ms. Fig. 9 plots the magnitude and GD frequency responses from MATLAB computations of TF (5) in from dc to 1.7 kHz. These results confirm the LP-NGD behavior varying from -1 ms to 0 ms in the LF interval [0 Hz,1 kHz] of the synthesized circuit. It is noteworthy that the gain at the cutoff frequency is approximately 17 dB. A limited bandwidth smoothed enough signal must be considered to realize the I E E E P r o o f

311 2 )

 2 Time<0 Demonstration With Numerical LP-NGD 312 Gaussian Response: The LP-NGD transient analysis must be 313 performed with smoothed signal presenting spectrum belonging 314 to the NGD frequency band. Accordingly, the Gaussian wave-315 form pulse as input test signal is analytically expressed as

  319 f a . Knowing the maximal instant time t 0 , the LP-NGD cir-320 cuit transient characterization can be performed by rise-time 321 advance τ r = t out r t in r , and fall-time advance

331 333 ζ 2 = 8 .Fig. 10 .

 3332810 Fig. 10. (a) Calculated and measured LP-NGD MCU transient Gaussian pulse responses. (b) Corresponding oscilloscope screenshot of tested signals.

Fig. 11 .

 11 Fig. 11. (a) Calculated and measured LP-NGD MCU transient sinc waveform signal responses. (b) Corresponding oscilloscope screenshot of tested signals.

372tion:

  Fig. 12. LP-NGD responses of (a) temperature and (b) humidity sensors.

406 1 )

 1 An innovative theory of LP-NGD IIR TF was established. The 407 main specifications of the input signal susceptible to operate 408 with the LP-NGD function were defined in function of the MCU 409 specifications. The design method and implementation of the 410 numerical LP-NGD function was described. Then, the canonical 411 fundamental discrete function of LP-NGD numerical circuit was 412 analytically introduced. The LP-NGD analysis and synthesis 413 methodology was originally formulated in function of the NGD 414 specifications. 415 To validate the numerical LP-NGD function, an original 416 experimental setup with STM32 MCU board was consid-417 ered as a POC. The numerical LP-NGD circuit was exper-418 imentally validated by millisecond short-duration and com-419 putationally validated by hour long-duration arbitrary wave-420 form signal testing. Results showed a very good agree-421 ment between the analytical computations and experimenta-422 tions were introduced and discussed. The numerical LP-NGD 423 function industrial application was promisingly highlighted 424 by arbitrary waveform sensored signals representing temper-425 ature and humidity 12.5-day duration data from embedded 426 sensors. 427 Compared with the traditional LP-NGD analog circuit, the 428 investigated numerical NGD circuit presents. 429 Benefits on i) the possibility to operate with long period 430 signals having a duration of more than ten seconds (which 431 cannot be treated with analog NGD circuits), ii) to inter-432 face the NGD circuit directly to digital signals, iii) to 433 increase the order of LP-NGD circuit in order to widen 434 the NGD bandwidth without increasing the physical size 435 of the PCB prototype, and 4) to adjust negative delay pro-436 grammatically according to the needs, making it versatile 437 for a wide range of applications.

438 2 )

 2 Disadvantages on i) the latency of the digital computa-439 tions and ii) the delay introduced by ADC and DAC im-440 pacting the NGD response. However, especially for long 441 NGD targeted by this digital implementation, conversion 442 delays are usually negligible. 443 With this article, we envisage a revolution NGD technology 444 for the future smart factories based on the sensored signal 445 anticipation. To do this, the NGD MCU concept can be integrated 446 in real-time simulation and deep learning interfaces[35],[36]. 447 I E E E P r o o f Blaise Ravelo (Member, IEEE) is currently a University Full Professor with the Nanjing University of Information Science and Technology, Nanjing, China. His research interests include on multiphysics and electronics engineering. Prof. Ravelo is a Pioneer of the Negative Group Delay (NGD) concept about t<0 signal travelling physical space. He was research director of 11 Ph.D. students (eight defended), postdocs, research engineers, and Master internships. With US, Chinese, Indian, European, and African partners, he is actively involved and contributes on several international research projects (ANR, FUI, FP7, INTERREG, H2020, Euripides², Eurostars …). He is Member of IET Electronics Letters Editorial Board as Circuit and System Subject Editor. He is ranked in Top 2% world's scientists based on years 2019 by Stanford University, US. He has Google scholar h-index (2021) = 23 and i10-index(2021) = 68. He is member of research groups: IEEE, URSI, GDR Ondes, Radio Society and (co-) authors of more than 360 scientific research papers in new technologies published in International Conference and Journals. He is Lecturer on Circuit and System Theory, STEM (science, technology, engineering, and maths) and applied physics. He is regularly invited to review papers submitted for publication to international journals (IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS, IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS, IEEE ACCESS, IET Circuits, Devices, and Systems, and IET Microwaves, Antennas and Propagation) and books (Wiley, Intech Science). Q4 Mathieu Guerin (Member, IEEE) received the 596 Engineering degree in microelectronics and 597 telecommunications from Polytech Marseille, 598 Marseille, France, and the Research Master de-599 gree in integrated circuits design from the Uni-600 versity of Aix-Marseille, Marseille, France, both 601 in 2010, and the doctorate degree from the Uni-602 versity of Aix-Marseille, in 2010. Q5 Q6 603 He was a Technical Leader of the Analog 604 and Radio-Frequency Design Team of IDEMIA-605 StarChip for five years and designed chips em-606 bedded in SIM cards and contactless bank cards with biometric recog-607 nition. He joined Aix-Marseille University as an Assistant Professor in 608 2020 and joined the CCSI Team of the IM2NP Laboratory. He is working 609 on methods of modeling and characterizing circuits in analog electron-610 ics. His research interest includes the design and synthesis of circuits in 611 digital electronics.

  Member, IEEE) re-614 ceived the B.Sc. degree in electrical engineer-615 ing from Nice Sophia-Antipolis University, Nice, 616 France, in 1996, the M.Sc. degree (with distinc-617 tion) in electrical engineering from the Science 618 Department, University of Montpellier, Montpel-619 lier, France, in 1998, and the Ph.D. degree in mi-620 croelectronics form the University of Montpellier. Q7 621 He is currently a Full Professor with the Uni-622 versity of Aix-Marseille, Marseille, France. Since 623 1998, he has been with the Microelectronics 624 Department of Informatics, Robotics and Microelectronics, Laboratory of 625 Montpellier (LIRMM). Since 2003, he has been with the Microelectronic 626 Department of Materials, Microelectronics and Nanoscience, Laboratory 627 of Provence (IM2NP), Marseille, France, where he was an Associate 628 Professor. Since 2014, he has been a Professor with Aix Marseille 629 University where he heads the Integrated Circuit Design group of the 630 IM2NP Laboratory. He is regularly involved to participate and to lead 631 national and international research projects (ANR, H2020, FP7 KIC-632 InnoEnergy …). He directed and cosupervised 18 Ph.D. and 15 Master 633 students. He has authored or coauthored 11 patents and more than 150 634 papers published in refereed journals and conferences. His research in-635 terests include AMS and RF circuit design from transistor to architectural 636 level. His current research interests include ultralow power circuit design 637 for smart sensor interface and embedded electronic in bioelectronic and 638 e-health applications, wireless systems, design technique, and architec-639 ture for multistandard transceiver. 640 Prof. Rahajandraibe is an Expert for the ANR, the French Agency 641 for Research. He has served on program committees of IEEE INTER-642 NATIONAL NEW CIRCUITS AND SYSTEMS CONFERENCE and IEEE INTER-643 NATIONAL CONFERENCE ON ELECTRONICS, CIRCUITS, AND SYSTEMS. He 644 has been and is a Reviewer of contributions submitted to several IEEE 645 conferences and journals such as IEEE INTERNATIONAL SYMPOSIUM ON 646 CIRCUITS AND SYSTEMS, IEEE INTERNATIONAL NEW CIRCUITS AND SYS-647 TEMS CONFERENCE, IEEE INTERNATIONAL MIDWEST SYMPOSIUM ON CIR-648 CUITS AND SYSTEMS, ESSCIRC, ESSDERC, RFIC, IEEE TRANSACTIONS 649 ON CIRCUITS AND SYSTEMS I AND II, IET Electronics Letters.

r

  Member, IEEE) received the 652 Graduate degree in applied physics from Ecole 653 Normale Supérieure, Paris, France, in 2001, 654 and the Ph.D. degree in electronics from Ecole 655 Nationale Supérieure de Techniques Avancées, 656 Paris XI Orsay University, Orsay, France, in 657 2005. 658 His Ph.D. was focused on both circuits and 659 algorithms for artificial Retinas. He is Lecturer 660 and Researcher with Toulon University since 661 2007, and an Associate Professor in Robotics, 662 Embedded Electronics and IoT at SeaTech, ISEN Toulon, and ENSTA 663 ParisTech. He was with the IM2NP CNRS Laboratory in 2017, in the 664 Circuits Design Team. His current research interests include embedded 665 algorithms and circuits for ultralow power systems. 666 Prof. Gies is Scientific Advisor of several start-ups in IoT and Head of 667 Scientific Microsystems for Internet of Things, Toulon University. 668 I E E E P r o o f 10 IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS Lala Rajaoarisoa (Member, IEEE) received the M.Sc. and Ph.D. degree in automatic and computer sciences from the University of Aix-Marseille, Marseille, France, in 2005 and 2009, respectively. He is currently an Assistant Professor with the Institut Mines-Télécom Lille Douai. Develop predictive models and controllers to assess system behavior and optimize its performance. This development includes the analysis of intrinsic properties such as stability, observability, identifiability, and controllability. He is involved in research activities dedicated to the optimization of energy efficiency of building systems and the control and management of hydraulic systems with more than 80 papers published in refereed journals and conferences. He regularly participates and contributes on several international projects (ANR, FUI, and INTERREG) and was the supervisor of more than of 15 Ph.D. students, postdocs, research engineers, and Master internships. His research interests include the development of data-driven tools and methods for the observation and control of large-scale distributed systems. Sébastien Lalléchère (Member, IEEE) was 691 born in Nevers, France, in 1979. He received the 692 M.Sc. and Ph.D. degrees in computational mod-693 eling and electronics/electromagnetism from 694 Polytech Clermont and Universite Blaise Pascal, 695 Clermont-Ferrand, France, in 2002 and 2006, 696 respectively. 697 He was Research Engineer with LASMEA, 698 Clermont-Ferrand, France, in 2007 focusing on 699 intensive computational methods for electro-700 magnetics. He is currently an Associate Profes-701 sor with the Institut Pascal and Universite Clermont Auvergne, Clermont-702 Ferrand, France. His research interests include electromagnetic compat-703 ibility including antennas and propagation, complex and reverberating 704 electromagnetic environments, electromagnetic coupling, computational 705 electromagnetics, stochastic modeling, and sensitivity analysis in elec-706 trical engineering. Authors: Carefully check the page proofs (and coordinate with all authors); additional changes or updates WILL NOT be accepted after the article is published online/print in its final form. Please check author names and affiliations, funding, as well as the overall article for any errors prior to sending in your author proof corrections. r Authors: We cannot accept new source files as corrections for your article. If possible, please annotate the PDF proof we have sent you with your corrections and upload it via the Author Gateway. Alternatively, you may send us your corrections in list format. You may also upload revised graphics via the Author Gateway. r Authors: Unless invited or otherwise informed, there is a mandatory Excessive Article Length charge of $250 per page ($200 for IES members) in excess of eight (8) pages. If you have any questions regarding overlength page charges, need an invoice, or have any other billing questions, please contact apcinquiries@ieee.org as they handle these billing requests.

28 I

 28 of the desired NGD value, cutoff 12 frequency, gain, and the MCU sampling frequency. A STM32 13 MCU proof-of-concept (POC) is tested to implement the 14 IIR LP-NGD function. Different real-time tests with visual-15 ization of input and output analog signals from the MCU 16 LP-NGD POC were performed. As expected, time-advance 17 demonstration tests with milli-second short-and several 18 hour long-duration time-scale with arbitrary waveform sig-19 nals from temperature and humidity sensors. The signal 20 time-advance is not in contradiction with the causality. The 21 proposed digital MCU function opens a potential future in-N EARLY age of Industry 4.0 revolution, the role of re- 29 mote communication is one of the design research engineer 30 breakthrough [1]. The delay effect plays a major role on the com-31 munication system general performance [2]-[4]. For example, 32 these effects are critical for reliable operation of remote-control 33 Manuscript received March 29, 2021; revised June 17, 2021 and July 3, 2021; accepted August 23, 2021. (Corresponding author: Mathieu Guerin.)

  system. The quality of service (QoS) of remote sensing systems forecasted Industry 4.0 factories to operate correctly depends undeniably on the delay effects. The networked modern industrial sensing systems present an open risk of anomaly detection because of unintentional delays. So far, no rigorous solution is available to overcome the time delay effects in almost all areas of engineering [2]-[9]. One can cite different challenges to be overcome to face up the delay impacts. The signal delay degrades the digital communication QoS due to the wireless propagation channel dispersion and data desynchronization [3]-[4]. In the area of automation and control engineering, different industrial application examples

Fig. 1 .

 1 Fig. 1. Illustration scenario of NGD application. (a) Delay penalizing system. (b) NGD corrected system. (b) Input and output signal chronogram.

  and implemented like classical and familiar electrical circuits as filters, amplifiers, phase shifters, oscillators. The LF NGD topologies proposed initially in[26]-[29] were built with classical electronic passive (R, L, and C) and operational amplifier components. The NGD phenomenon meaning was illustrated with time-domain experimentations showing the possibility to propagate the NGD output signal in time-advance with the input one[26]-[29]. This demonstration was realized with PCB showing an NGD output lamp lightened on before the input one[29]. This counterintuitive NGD phenomenon is not in contradiction with the causality principle[26]-[29]. For further understanding about the NGD concept, an investigation on the similitude between the NGD and linear filter gain behaviors is introduced in[30]. The concept of low-pass (LP)[31],[32] NGD function was initiated. Despite the research work progress in the area of electronics engineering, the NGD function remains an unmatured concept for most of industrial engineers. Open questions are regularly asked about the industrial potential applications of the NGD function. For this reason, a completely original research work on the design of numerical LP NGD circuit is developed in the present paper. The LP-NGD function is implemented on microcontroller unit (MCU) operating in different time scales. The NGD function can be potentially used as a pioneer work of anticipator

Fig. 2 .

 2 Fig. 2. Specific characterization diagram of numerical LP-NGD input signal.II. MCU DESIGN PRINCIPLE OF DISCRETE

  129

151LP-NGD Function 152

 152 This section describes the main specifications and character-153 ization of numerical LP-NGD function.

154 1 )

 1 NGD Specifications: The NGD TF can be represented as 155 a single-input single-output block diagram, as shown in Fig.3(a). 156 With the Laplace variable s, the associated analog TF can be 157 formulated by N (s) = V out (s)/V in (s). The NGD analysis in the 158 frequency domain is based on the consideration of magnitude 159 response, N (ω) = |N (jω)| and the GD defined by 160 τ (ω) = -∂ϕ(ω)/∂ω (2) with 161 ϕ(ω) = arctan {Im [N (jω)] /Re [N (jω)]} .

Fig. 3 .

 3 Fig. 3. (a) Block, (b) GD, and (c) magnitude response diagrams of LP-NGD function.

Fig. 3 (

 3 Fig. 3(b) represents the GD diagram specific to the ideal LP-NGD function, which is defined by two fundamental characteristics as follows.

Fig. 4 .

 4 Fig. 4. LP-NGD TF magnitude flatness versus product ω n τ n in (a) semilog x-axis, (b) log-log plots, and (c) ratio from (9).

Fig. 4 Fig. 5 .

 45 Fig. 5. Variations of LP-NGD IIR coefficients versus product ω n T s .

  212 condition |τ n | = 4T s . The discrete TF is designed by taking into 213 account to the initial condition214

223

  ations are semi-logarithmically plotted in Fig.5for the particular 224 case of NGD cutoff frequency-value product ω n τ n = 2π or f n 225 = 1/τ n . As the present study is aimed to the NGD analysis with 226 LF signal and high sampling rate, we choose the small range of 227 NGD cutoff frequency-sampling rate period product ω n T s from 228 10 -6 to 1. Coefficient ζ 0 is decreasing and the two other ones ζ 1 229 and ζ 2 increases with product ω n T s . 230 It is worth noting that the initial value of the discrete variable 231 allowing to realize a reasonable prediction must be 232

Fig. 8 (

 8 Fig. 8(a) highlights the real-time experimental setup of the STM32 MCU based IIR LP-NGD function. As indicated by

Fig. 7 .

 7 Fig. 7. Numerical LP-NGD circuit code.

Fig. 8 .

 8 Fig. 8. (a) Illustrative diagram and (b) photograph of the LP-NGD MCU-based experimental setup.

Fig. 8 (Fig. 9 .

 89 Fig. 8(a), the input pulse analog voltage provided by the signal268

1 )

 1 LP-NGD Function Frequency Responses: Before the transient test, the LP-NGD circuit were synthesized under the desired specifications, static gain N 0 = 0 dB, cutoff frequency f n = 1 kHz, and NGD value τ n = -1 ms. Fig. 9 plots the magnitude and GD frequency responses from MATLAB computations of TF (5) in from dc to 1.7 kHz. These results confirm the LP-NGD behavior varying from -1 ms to 0 ms in the LF interval [0 Hz,1 kHz] of the synthesized circuit. It is noteworthy that the gain at the cutoff frequency is approximately 17 dB. A limited bandwidth smoothed enough signal must be considered to realize the

311 2 )

 2 Time<0 Demonstration With Numerical LP-NGD 312 Gaussian Response: The LP-NGD transient analysis must be 313 performed with smoothed signal presenting spectrum belonging 314 to the NGD frequency band. Accordingly, the Gaussian wave-315 form pulse as input test signal is analytically expressed as 316

  319 f a . Knowing the maximal instant time t 0 , the LP-NGD cir-320 cuit transient characterization can be performed by rise-time 321 advance τ r = t out r t in r , and fall-time advance τ f = t out f -322 t in f , with v in (t in r,f < t 0 ) = max[v in (t)]/2 and v out (t out r,f < t 0 ) = 323 max[v out (t)]/2, combined to voltage gain between the input and 324 output signals 325

331 333 ζ 2 = 8 .Fig. 10 .

 3332810 Fig. 10. (a) Calculated and measured LP-NGD MCU transient Gaussian pulse responses. (b) Corresponding oscilloscope screenshot of tested signals.

Fig. 11 .

 11 Fig. 11. (a) Calculated and measured LP-NGD MCU transient sinc waveform signal responses. (b) Corresponding oscilloscope screenshot of tested signals.

4 )

 4 Hour Time-Scale Arbitrary Sensored Signal Predic-372 tion: An illustrative industrial application case of LP-NGD 373 function is described in this paragraph by using realm data 374 recorded from temperature and humidity sensors. The realm data 375 used in this paragraph were collected by sensors deployed on a 376 student residence located in Douai city, in the north of France 377 [39]. The adopted technology is wireless sensor network plat-378 forms based on a universal solution of the CLEODE company. 379 Each smart-meter/device is set with a sampling time of T s = 2.4 380 min, and its characteristics are given in Table III. Thus, these 381 configurations allow us to monitor the external environment 382 of buildings. Fig. 12 illustrate the evolution of the outdoor 383 climate during the experimentation (from Feb. 9, 2019 to Feb.

384 21 ,Fig. 13 . 1 )

 21131 Fig. 12. LP-NGD responses of (a) temperature and (b) humidity sensors.

438 2 )

 2 Disadvantages on i) the latency of the digital computa-439 tions and ii) the delay introduced by ADC and DAC im-440 pacting the NGD response. However, especially for long 441 NGD targeted by this digital implementation, conversion 442 delays are usually negligible. 443 With this article, we envisage a revolution NGD technology 444 for the future smart factories based on the sensored signal 445 anticipation. To do this, the NGD MCU concept can be integrated 446 in real-time simulation and deep learning interfaces [35], [36]. 447 I E E E P r o o f

  Given 145 the case of input signal duration t max , the MCU must operate 146 with time step following the integer index k = 01,2, …, k max = 147 Ent[t max /T s ] by denoting Ent[x] the superior integer part of real 148 positive x. The discrete input and output signals at the instant 149 time t = kT s are represented by

	150

  LP-NGD function discretization is established from the 215 differential infinitesimal elements, dv in (t)/dt and dv out (t)/dt.

	216
	217	By denoting the integer k = 01,2, …,k max , we have the discrete
	218	quantities expressed as
		Δv in (kT s )/Δt = [v in (k + 1) -v in (k)] /T s	(16)

TABLE II SPECIFICATIONS

 II OF THE INPUT GAUSSIAN SIGNALS, AND COMPARISON OF SIMULATED AND CALCULATED LP-NGD TRANSIENT SPECIFICATIONS time-domain effect of the LP-NGD function. For this reason,

307 the preliminary frequency domain analysis is necessary before 308 the LP-NGD in the time-domain. To complete validation of the 309 numerical LP-NGD function, time-domain measurement results 310 will be examined in the following paragraph.

TABLE IV

 IV

		8	IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS
	395	and (b). The cross correlation between the input-output tem-
	396	peratures and humidity's are 81.8% and 92.4%, as tabulated in
		NGD FUNCTION SPECIFICATIONS, AND TEMPERATURE AND HUMIDITY CROSS-CORRELATION P r o o f
	I E E E
			Fig. 13. Zoom in plots of LP-NGD responses of (a) [50 h, 100 h] and
			(b) [150 h, 180 h] time-interval temperature, and (c) [50 h, 100 h] and
			(d) [150 h, 180 h] time-interval humidity sensors.

in Table

IV

, τ n = -3600 s = -1 h and cutoff frequency,

f n = 278 μHz.

The feasibility study is based 300 h (h) or 12.5 days duration temperature and humidity measured sensored signals. The numerical LP-NGD computed responses are presented in Fig.

12(a)

Table IV ,

 IV respectively. As illustrated by the time-domain results

	397
	398	zoomed in time-intervals [50 h, 100 h] [see Fig. 13(a) and (c)],
	399	and [150 h, 180 h] [see Fig. 13(b) and (d)], an exciting result
	400	showing prediction of signal waveform is observed with output
	401	advance compared to the input under medium time-advance of
	402	about 1 h.
	403	IV. CONCLUSION
	404	The first-time investigation on design and test of LP-NGD
	405	numerical circuit was developed. The LP-NGD function was

  Wenceslas Rahajandraibe , Member, IEEE, Valentin Gies, Member, IEEE, Lala Rajaoarisoa , Member, IEEE, and Sébastien Lalléchère , Member, IEEE

1 2 Blaise Ravelo , Member, IEEE, Mathieu Guerin, Member, IEEE,

  LP-NGD function discretization is established from the 215 differential infinitesimal elements, dv in (t)/dt and dv out (t)/dt. By denoting the integer k = 01,2, …,k max , we have the discrete

	217
	218	quantities expressed as
		Δv in (kT s )/Δt = [v in (k + 1) -v in (k)] /T s	(16)

216

TABLE II SPECIFICATIONS

 II OF THE INPUT GAUSSIAN SIGNALS, AND COMPARISON OF SIMULATED AND CALCULATED LP-NGD TRANSIENT SPECIFICATIONS time-domain effect of the LP-NGD function. For this reason,

	307
	308	the preliminary frequency domain analysis is necessary before
	309	the LP-NGD in the time-domain. To complete validation of the
		numerical LP-NGD function, time-domain measurement results

310

will be examined in the following paragraph.

TABLE III TEMPERATURE

 III AND HUMIDITY WIRELESS SENSOR CHARACTERISTICS

V out (s) -N n V in (s) = s [N n τ a V in (s)τ b V out (s)] .(13)

0278-0046 © 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See https://www.ieee.org/publications/rights/index.html for more information.