Georges Gielen

Wim Dehaene

Philips Phillip Christie

Dieter Draxelmayr

Edmond Janssens

Karen Maex

Ted Vucurevich

Analog and Digital Circuit Design in 65 nm CMOS: End of the Road?

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

With the evolution towards ultra-deep-submicron and nanometer CMOS technologies [START_REF]International Technology Roadmap for Semiconductors[END_REF] the design of complex Systems on a Chip (SoC) is emerging in consumer-market applic ations such as telecom and multimedia. These integrated systems are increasingly mixed-signal designs, embedding high-performance analog or mixed-signal blocks and possibly sensitive RF frontends together with complex digital circuitry (multiple processors, some logic blocks, and several large memory blocks) on the same chip. In addition, the growth of wireless services and other telecom applications increases the need for low-cost highly integrated solutions with very demanding performance specifications.

The use of nanometer CMOS technologies (below 90nm) however also brings along significant challenges for circuit design (both analog and digital). Some of these challenges were never encountered before, while others existed before but have become even stronger limitations today. As an example, in this paper the 65 nm technology node is used but the observations will only be more true if technology scaling is taken to even smaller dimensions.

Technology scaling down to 0.18 µm was based on the following reasoning. Reducing the feature size in the technology front-end (i.e. of the devices) and in the technology back-end (i.e. of the interconnect) combined with the addition of ever more interconnect layers, drastically increased the density of the digital circuits while reducing the intrinsic gate switching delay. The main disadvantage of this is that the power supply voltage must be decreased too, potentially leading to an increase in gate switching delay. However, the threshold voltage (VT) of the MOS transistors is lowered to compensate for this. Assume a full scaling scenario. This means that the scaling factor 1/S is the same for all geometry parameters and all voltage parameters. This ideal case leads for standard CMOS gates to a density increase of S 2 , a decrease of the intrinsic delay of 1/S, a decrease in power consumption of 1/S 2 at a constant power density. The noise margin of the logic is decreasing but remains acceptable. This ideal scenario cannot be realized in practice but it shows that for digital circuitry scaling is (was) advantageous: it brings faster, denser logic for less power and an acceptable robustness penalty [START_REF] Rabaey | Digital integrated circuits (2 nd edition)[END_REF]. For the analog circuits, scaling did not bring much area reduction, but mainly provided the transistors with extra speed, allowing the silicon implementation of RF circuits and high-speed analog blocks like data converters.

However, as technology scaling goes on to 90 nm and below (65nm, 45nm) , physical and quantummechanical effects that were previously not relevant become influential or even dominant. An example of this is the advent of no longer negligible leakage currents. Secondly as all geometrical dimensions become smaller, variability of technological parameters also gains in importance. The same absolute tolerance becomes relatively more important when the absolute value of a parameter decreases. This holds for VT, doping levels, widths and lengths, … A V T variation of 50 mV has more effect on the circuits when V T is 200 mV compared to a V T of 700 mV for instance. The central question in this discussion is whether the scaling advantages are still realized under the presence of these new effects. In other words: will the design of circuits in nanometer technologies of 65 nm and beyond still bring the benefits we are used to from the past, or will the design bottlenecks become so large that the design road ends at the 65nm technology node ? This paper will give an introductory overview of problems encountered in nanometer technologies.

The paper is organized as follows. Section 2 gives a more concrete description of the effects that nanometer scaling has on the circuit behavior in general. The 65nm node is used as an example. Section 3 then presents in detail the challenges in digital circuit design, while section 4 focuses on challenges in analog circuit design. Finally, section 5 presents some conclusions.

General consequences of 65 nm CMOS at the circuit level

Scaling technology to nanometer sizes has brought effects to the surface that were negligible before. The most important ones are discussed here: subthreshold leakage current, gate leakage current, interconnect delay and increased process variability. Their effect on circuit behavior will be described in the next subsections.

Subthreshold leakage current

The current conducted by a MOS transistor when it is off (V GS =0) is given by :

) exp( / 0 q nkT V ld subthresho T I I - = (1)
The subthreshold current thus increases exp onentially with decreasing V T . Moreover I 0 in this equation is invers ely proportional to L. Nanometer scaling thus leads to a dramatic increase in leakage current up to the point where it can no longer be ignored in the power consumption nor the functionality of the circuits.

For very short transistors the subthreshold leakage not only increases because L scaling implies V T scaling. When the drain region of a MOS transistor is close to the source region, the depletion region around the source and drain regions interact with each other. This lowers the potential barrier at the source side. This is usually modeled with an equivalent, V ds -dependent, V T decrease, as illustrated in Fig. 1.

DIBL ⇒ ∆V T <0

Fig. 1. Subthreshold current dependency on V ds and V gs .

Gate leakage current

Before the nanometer era, gate currents were considered as pure dynamic. MOS gates were capacitors and did not conduct DC currents. However, since the advent of oxide thicknesses of only a few nm, current can tunnel through the gate and have a static "leakage" component. The gate leakage current is described by :

) / exp( 2 _ gb ox ox gb leakage gate V t t V KW I α -         = ( 2
)
where K and α are fit factors. Remark that gate leakage is only present when there is a voltage across the gate, i.e.

when the transistor is on. Subthreshold leakage on the contrary occurs when a transistor is off. Gate leakage currents affect the power consumption and functionality at the circuit level. Using new technologies not only means scaling of the size. New materials are also introduced. A lot of research attention goes to the introduction of high-k dielectrics for the gate. This would lead to reduced gate leakage because thicker "oxides" can be used without compromising gate capacitance and related to that VT. A similar story holds for the (re -)introduction of metal gates.

Interconnect delay and scaling

A simplified, first-order approximation of interconnect delay is given by :

2 2 2 λ ρκ       = = L rcL t wire (3)
where r and c are the resistance and capacitance per unit length, ρ and κ are the same per unit area, and λ is the technology-related wire pitch. This shows that the interconnect delay is inversely proportional to the square of the wire pitch ? of the technology. This pitch scales down with the technology. This means that if the interconnect length (l) and interconnect pitch (?) scale identically, the wire delay will remain constant with technology scaling. As stated in the introduction, the intrinsic gate delay of a technology decreases with technology scaling. This means that the interconnect delay will gain in relative importance when the technologies are scaled down.

Interconnect length and pitch however scale identically for longer interconnects, e.g. busses. The length of these lines has the tendency to remain constant even as technology scales. For that case the interconnect delay becomes dominant even faster.

A similar reasoning can be built up for the power that is consumed due to the interconnect capacitance. Also the relative importance of this contribution to the power consumption is increasing. Also in this case the introduction of new materials could bring some relieve: low-k materials for the isolation between the interconnects for instance, or better conductor materials. [START_REF] Frank | Parameter variations in sub-100nm technology[END_REF].

Fig. 2. Number of dopant atoms vs. channel length

Process variability

The tolerance on some technological parameters will increase when the technology scales. There are numerous examples of this , each with a different physical explanation. Two examples will be given here.

The first example is the discrete nature of doping levels. The statistical variation of the number of dopants N varies with N 1/2 . This increases the uncertainty on V T for small N. Remark that for decreasing channel dimensions, the number of dopants that realizes a given doping level also decreases (see Fig. 2). The random placement of the dopants creating the source and drain also causes an uncertainty on the effective channel length (see Fig. 3). This effect is also enforced as the number of dopants goes down.

Fig. 3. Representation of the source/drain dopants in

a MOS transistor [START_REF] Frank | Simulation of Stochastic Doping Effects in Si MOSFETs[END_REF].

A second example is line edge roughness. As the L of a transistor scales down , the same roughness on the edge of a transistor becomes relatively more important, leading to a larger variation on Leff and thus on the transistor current.

Process variability is split in inter-die and intra-die variations. The latter are also referred to as device matching. Both types of variability are affected by the nanometer scaling. However the methods that can be used at the circuit level to deal with both types are different.

Digital circuit design challenges in 65 nm

In this section some examples will be given to illustrate the influence of nanometer technology scaling on digital circuit design.

The influence of VT variability on the energy-delay trade-off

The influence of V T variations on the delay of a logic gate increases as the technology becomes smaller, as is illustrated in Fig. 4. This is due to the fact that the overdrive voltages (V dd -V T ) that are used also become smaller. Therefore the relative importance of ? V T on V dd -V T grows.

The effect of this V T variation is compensated for by designing the circuit such that the delay requirements are still fulfilled in the worst case (largest V T in this case). This however has a penalty on the energy consumption.

The sizing of the circuit, and thus the parasitic capacitance, will be larger than it needs to be when the V T is not at its largest. However the dynamic energy consumption is determined by CV dd 2 , independent of V T . The effect of sizing on a circuit's energy consumption is larger than the effect on the circuit's delay. Therefore, the effect of "worst-case oversized" design on the energy consumption of circuits will be significant. [START_REF] Kobayashi | Self-adjusting thresholdvoltage scheme (SATS) for low-voltage high-speed operation[END_REF].

Fig. 4. Influence of V T variations on gate delay

Leakage-aware digital circuits

As explained above, nanometer scaling increases the leakage currents of MOS transistors up to a level where they can no longer be ignored. To deal with this, a range of techniques has been developed and are still under development. A first class of techniques called MTCMOS (multi threshold CMOS) increase the threshold voltage of transistors in the circuit that are not on a critical timing path. This reduces their subthreshold leakage without affecting the speed of the circuit. To apply this technique multiple V T technologie s are required. VTCMOS is usually combined with supply and/or ground switches that switch off the leaky circuits when they are inactive.

A second class of techniques is called VTCMOS (variable threshold CMOS). In this technique the V T is influenced via the bias of the transistors body. There are a number of strategies to adapt the body bias (see e.g. [START_REF] Sakurai | Adaptive Circuit Techniques for Managing Variations[END_REF]). There is however one problem with this technique: as technology scales down, the bulk factor becomes smaller. This means that the effect of body bias on V T reduces, thus limiting the effectiveness of the VTCMOS technique.

The examples given above suggest that there is a point where further scaling of the intrinsic MOS device is not really meaningful anymore. The delay still decreases but the penalty in terms of energy consumption due to leakage and worst case design will become unacceptable. Furthermore, the problems posed to the circuit designer, mainly due to intra-die matching, are hard to deal with. It is more than likely that the circuits needed to solve these problems will also consume power and silicon area…

Increasing interconnect delay and synchronous design

As the relative importance of interconnect delay tends to grow, it becomes harder and harder to maintain the synchronism between different parts of a digital design. In a typical 100 nm technology the max length of a wire is around 2 mm to keep the skew below 20 % of a 1 GHz clock (see Fig. 5). With decreasing interconnect pitches and line widths, this distance will also decrease. This calls for architectural solutions: localization of storage and computation, complemented with a communication strategy that is locally synchronous but globally asynchronous. This will again lead to power and silicon area overhead along with an increased design complexity.

Analog circuit design challenges in 65 nm

In this section some challenges will be discussed that are posed on analog circuit design in nanometer technology scaling.

The influence of reducing supply voltages

For analog circuits technology scaling does not bring large area reductions, because the active area (width times length) of key analog transistors in a circuits is determined by noise or mismatch constraints, both of which limit the dynamic range or accuracy levels that can be achieved. As today's applications tend to increase rather than relax the accuracy requirements, this means that the area of most analog blocks does not really become significantly smaller with technology scaling. In general, both the kT/C thermal noise as the mismatchinduced offset result in a relationship between achievable speed, dynamic range and power consumption of a circuit according to the following formula [START_REF] Kinget | Impact of transistor mismatch on the speed-accuracy-power trade-off of analog CMOS circuits[END_REF] :

const techn Power Accuracy Speed = * 2 (4)
In the case of thermal noise, the right-hand side constant depends only on temperature; in the case of mismatch it depends on the amount of mismatch in the technology process used. These relationships are plotted for some real technology process in Fig. 6. Clearly, for untrimmed or uncalibrated circuits, the mismatch limit is determining the minimum required power consumption for given speed and dynamic range requirements. The red squares mark real A/D converter designs. When the technology scales, then the transistor mismatch improves slightly. Hence, if one wants to exploit the higher speed offered by the scaled technology, then this will come at an increased power consumption for the same level of dynamic range. On the other hand, for fixed speed and fixed accuracy, the power would decrease due to the improved matching. However, this is only true when ignoring the reduction in supply voltage which comes along with employing nanometer technologies. This reduced supply voltage also reduces the input signal range, hence making the constraint on thermal noise and offset even more stringent. As a result, even for fixed speed and fixed accuracy, the power consumption will no longer decrease with scaling, but will remain the same or may even increase slightly (see Fig. 7), according to :

where P 1 , P 2 are the power consumption in technology 1 and 2 respectively, t ox1 and t ox2 the oxide thickness in both technologies, and m the ratio of the supply voltages. Hence, there is no real benefit for using nanometer technologies in this case. In addition, the reduced supply voltage also makes it difficult to implement analog circuits, since many circuit techniques like cascading and device stacking become no longer possible.

Fig. 7. Power consumption with technology scaling for analog circuits with fixed speed and accuracy requirements (red curve).

Also the increased process parameter variability has an impact on analog circuit performance, but analog designers have always had to cope with process tolerances and mismatches, and have been using statistical methods already a long time ago [START_REF] Director | Statistical integrated circuit design[END_REF].

Boosting analog design productivity

Due to the knowledge-intensive nature of analog design, most analog designs today are still handcrafted manually by analog expert designers, with only a SPICElike simulation shell and an interactive layout environment (with parameterized procedural device generators) as supporting facilities. This makes the design cycle for analog circuits long and error-prone. Therefore, although analog circuits typically occupy only a small fraction of the total area of mixed-signal ICs and SoCs , their design is often the bottleneck in mixed-signal systems, both in design time and effort as well as test cost, and they are often responsible for design errors and expensive reruns. This handcrafting is also increasingly at odds with the shortening time -to-market constraints of current consumer market products. This explains the growing need observed in industry today for analog CAD tools that boost analog design productivity by assisting designers with fast and first-time-correct design of analog circuits, or even by automating certain tasks or the entire circuit design process where possible.

While the basic level of design abstraction for analog circuits is mainly still the transistor level, commercial CAD tool support for analog cell-level circuit and layout synthesis is emerging. There has been remarkable progress at research level over the past decade, and in recent years several comme rcial offerings have appeared on the market. [START_REF] Gielen | Computer-aided design of analog and mixed-signal integrated circuits[END_REF] offers a fairly complete survey of the domain. Most of the basic techniques in both circuit and layout synthesis today rely on powerful numerical optimization engines coupled to "evaluation engines" that qualify the merit of some evolving analog circuit or layout candidate. These tools have been shown to produce designs comparable or better than manual designs in much shorter design times. As an example, Fig. 8 shows a particle/radiation detector frontend that has been generated with the AMGIE/LAYLA analog synthesis tools [START_REF] Van Der Plas | AMGIE -A synthesis environment for CMOS analog integrated circuits[END_REF].

Mixed-signal signal integrity

A difficult problem in mixed-signal designs, where sensitive analog and RF circuits are integrated on the same die with large digital circuitry, is signal integrity analysis, i.e. the verification of all unwanted signal interactions through crosstalk or couplings at the system level t hat can cause parametric malfunctioning of the chip. Parasitic signals are generated (e.g. digital switching noise) and coupled into the signal of interest, degrading or even destroying the performance of the analog/RF circuitry. These interactions can come from capacitive or (at higher frequencies) inductive crosstalk, from supply line or substrate couplings, from thermal interactions, from coupling through the package, from electromagnetic interference (EMC/EMI), etc.

Especially the analysis of digital switching noise that propagates through the substrate shared by the analog and digital circuits has received much attention in recent years [START_REF] Donnay | Analysis and reduction techniques for substrate noise coupling in mixed-signal integrated circuits[END_REF]. At the instants of switching, digital circuitry can inject spiky signals into the substrate, which then will propagate to and be picked up by the sensitive analog/RF circuits. As exa mple, consider a VCO at 2.3 GHz and a digital circuit block (250kgates) running at 13 MHz. As shown on the measurement plot of Fig. 9, the digital clock is visible as FM modulation around the VCO frequency and may cause conflicts with out-of-band emission requirements [START_REF] Leenaerts | CAD solutions and outstanding challenges for mixedsignal and RF IC design[END_REF]. Fig. 9. Measured FM modulation due to substrate switching noise coupling [START_REF] Leenaerts | CAD solutions and outstanding challenges for mixedsignal and RF IC design[END_REF].

In recent years research has been going on to find efficient yet accurate techniques to analyze these problems. For the noise propagation through the substrate, typically finite difference methods or boundary element methods are used to solve for the substrate potential distribution due to injected noise sources [START_REF] Donnay | Analysis and reduction techniques for substrate noise coupling in mixed-signal integrated circuits[END_REF], allowing to simulate the propagation of digital switching noise injected in the substrate to sensitive analog nodes elsewhere in the same substrate. This propagation analysis then has to be combined with an analysis of the (signaldependent) digital switching activity to know t he actual (time -varying) injected signals, and with an analysis of the impact of the local substrate voltage variations on the analog/RF circuit performance (e.g. the reduction of the effective number of bits of an embedded analog-to-digital converter) in order to cover the entire problem.

As an example, the SWAN methodology determines the switching noise that is generated by the digital circuits in a system, by a-priori characterizing every cell in a digital standard cell library with a macromodel that includes the current injected in the substrate due to an input transition, and then calculating the total injection of a complex system by combining the contributions of all switching cells over time depending on the event information obtained from a VHDL simulation of the system [START_REF] Van Heijningen | Substrate noise generation in complex digital systems: efficient modeling and simulaton methodology and experimental verification[END_REF]. Fig. 10 for example shows the comparison between time -domain SWAN simulations and measurements on a large experimental WLAN SoC with 220k gates, that contains a scalable OFDM -WLAN baseband modem and a low-IF digital IQ (de)modulator, fabricated in a 3.3 V 0.35 µm CMOS 2P5M process on an EPI-type substrate. Compared to the measurements, the simulated substrate-noise voltage from zero to 100 ns is within an error of 20% in its RMS value and within an error of 4% in its peak-to-peak value, which is a very good result for a difficult crosstalk effect like substrate noise couplings. Techniques to analyse the impact of this on the performance of the embedded analog blocks are also being developed, but still require further work [START_REF] Donnay | Analysis and reduction techniques for substrate noise coupling in mixed-signal integrated circuits[END_REF]. In future nanometer technologies, however, also other signal integrity problems will show up that need to be analysed and modeled. These include electromagnetic interactions, EMC/EMI, etc.

Conclusions

This paper has presented an overview of the challenges at hand when designing integrated electronic systems in nanometer-scale CMOS technologies. Key problems that affect circuit design with scaling technologies include increased leakage and variability. Next, the impact of this on digital circuit design has been discussed. Finally, problems bothering embedded analog circuits have been presented, such as reducing supply voltages, poor design productivity and signal integrity troubles. Whether we manage to address these problems or not will determine whether the design road ends at CMOS technology marker "65 nm" or not.

Fig. 5 .

 5 Fig. 5. Max interconnect length for 20% clock skew as a function of the clock frequency for a typical M1,M2 interconnect in a 100 nm technology.

Fig. 6 .

 6 Fig.[START_REF] Sakurai | Adaptive Circuit Techniques for Managing Variations[END_REF]. Thermal noise and mismatch limit in the powerspeed-accuracy trade-off governing analog circuits.

Fig. 8 .

 8 Fig. 8. Layout of a particle/radiation detector frontend generated with the AMGIE/LAYLA analog synthesis tools [10].

Fig. 10 .

 10 Fig.[START_REF] Van Der Plas | AMGIE -A synthesis environment for CMOS analog integrated circuits[END_REF]. Measured and SWAN[START_REF] Van Heijningen | Substrate noise generation in complex digital systems: efficient modeling and simulaton methodology and experimental verification[END_REF] simulated substrate noise in an experimental 220k-gates WLAN SoC .

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE


Alexandre Bône
email: alexandre.bone@icm-institute.org

Olivier Colliot
email: olivier.colliot@icm-institute.org

Stanley Durrleman
email: stanley.durrleman@icm-institute.org

Learning distributions of shape trajectories from longitudinal datasets: a hierarchical model on a manifold of diffeomorphisms

We propose a method to learn a distribution of shape trajectories from longitudinal data, i.e. the collection of individual objects repeatedly observed at multiple time-points. The method allows to compute an average spatiotemporal trajectory of shape changes at the group level, and the individual variations of this trajectory both in terms of geometry and time dynamics. First, we formulate a non-linear mixedeffects statistical model as the combination of a generic statistical model for manifold-valued longitudinal data, a deformation model defining shape trajectories via the action of a finite-dimensional set of diffeomorphisms with a manifold structure, and an efficient numerical scheme to compute parallel transport on this manifold. Second, we introduce a MCMC-SAEM algorithm with a specific approach to shape sampling, an adaptive scheme for proposal variances, and a log-likelihood tempering strategy to estimate our model. Third, we validate our algorithm on 2D simulated data, and then estimate a scenario of alteration of the shape of the hippocampus 3D brain structure during the course of Alzheimer's disease. The method shows for instance that hippocampal atrophy progresses more quickly in female subjects, and occurs earlier in APOE4 mutation carriers. We finally illustrate the potential of our method for classifying pathological trajectories versus normal ageing.

Introduction 1.Motivation

At the interface of geometry, statistics, and computer science, statistical shape analysis meets a growing number of applications in computer vision and medical image analysis. This research field has addressed two main statistical questions: atlas construction for cross-sectional shape datasets, and shape regression for shape time series. The former is the classical extension of a mean-variance analy-sis, which aims to estimate a mean shape and a covariance structure from observations of several individual instances of the same object or organ. The latter extends the concept of regression by estimating a spatiotemporal trajectory of shape changes from a series of observations of the same individual object at different time-points. The emergence of longitudinal shape datasets, which consist in the collection of individual objects repeatedly observed at multiple timepoints, has raised the need for a combined approach. One needs statistical methods to estimate normative spatiotemporal models from series of individual observations which differ in shape and dynamics of shape changes across individuals. Such model should capture and disentangle the inter-subject variability in shape at each time-point and the temporal variability due to shifts in time or scalings in pace of shape changes. Considering individual series as samples along a trajectory of shape changes, this approach amounts to estimate a spatiotemporal distribution of trajectories, and has potential applications in various fields including silhouette tracking in videos, analysis of growth patterns in biology or modelling disease progression in medicine.

Related work

The central difficulty in shape analysis is that shape spaces are either defined by invariance properties [21,40,[START_REF] Su | Rateinvariant analysis of trajectories on riemannian manifolds with application in visual speech recognition[END_REF] or by the conservation of topological properties [5,[START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF][START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF][START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF]34], and therefore have intrinsically a structure of infinite dimensional Riemannian manifolds or Lie groups. Statistical Shape Models [9] are linear but require consistent points labelling across observations and have no topology preservation guarantees. A now usual approach is to use the action of a group of diffeomorphisms to define a metric on a shape space [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF]43]. This approach has been used to compute a "Fréchet mean" together with a covariance matrix in the tangent-space of the mean [START_REF] Allassonnière | Bayesian mixed effect atlas estimation with a diffeomorphic deformation model[END_REF][START_REF] Fletcher | Geodesic regression and the theory of least squares on riemannian manifolds[END_REF][START_REF] Gori | A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes[END_REF]33,44] from a cross-sectional dataset, and regression from time series of shape data [4, 15,[START_REF] Fletcher | Geodesic regression and the theory of least squares on riemannian manifolds[END_REF][START_REF] Hinkle | Polynomial Regression on Riemannian Manifolds[END_REF][START_REF] Lorenzi | 4D registration of serial brains MR images: a robust measure of changes applied to Alzheimer's disease[END_REF][START_REF] Niethammer | Geodesic regression for image time-series[END_REF]. In [START_REF] Allassonnière | Construction of bayesian deformable models via a stochastic approximation algorithm: a convergence study[END_REF][START_REF] Durrleman | Spatiotemporal atlas estimation for developmental delay detection in longitudinal datasets[END_REF]31], these tools have been used to estimate an average trajectory of shape changes from a longitudinal dataset using the convenient assumption that the parameters encoding inter-individual variability are independent of time. The work in [27] introduced the idea to use the parallel transport to translate the spatiotemporal patterns seen in one individual into the geometry of another one. The co-adjoint transport is used in [START_REF] Singh | Hierarchical geodesic models in diffeomorphisms[END_REF] for the same purpose. Both estimate a group average trajectory from individual trajectories. The proposed models do not account for inter-individual variability in the time dynamics, which is of key importance in the absence of temporal markers of the progression to align the sequences. The same remarks can be applied to [START_REF] Chakraborty | Statistics on the space of trajectories for longitudinal data analysis[END_REF], which introduces a nice theoretical setting for spaces of trajectories, in the case of a fixed number of temporal observations across subjects. The need for temporal alignment in longitudinal data analysis is highlighted for instance in [START_REF] Durrleman | Toward a comprehensive framework for the spatiotemporal statistical analysis of longitudinal shape data[END_REF] with a diffeomorphismbased morphometry approach, or in [40,[START_REF] Su | Rateinvariant analysis of trajectories on riemannian manifolds with application in visual speech recognition[END_REF] with quotient manifolds. In [START_REF] Koval | Statistical learning of spatiotemporal patterns from longitudinal manifold-valued networks[END_REF][START_REF] Schiratti | Learning spatiotemporal trajectories from manifoldvalued longitudinal data[END_REF], a generative mixed-effects model for the statistical analysis of manifold-valued longitudinal data is introduced for the analysis of feature vectors. This model describes both the variability in the direction of the individual trajectories by introducing the concept of "expparallelization" which relies on parallel transport, and the pace at which those trajectories are followed using "timewarp" functions. Similar time-warps are used by the authors of [START_REF] Kim | Riemannian nonlinear mixed effects models: Analyzing longitudinal deformations in neuroimaging[END_REF] to refine their linear modeling approach [22].

Contributions

In this paper, we propose to extend the approach of [START_REF] Schiratti | Learning spatiotemporal trajectories from manifoldvalued longitudinal data[END_REF] from low-dimensional feature vectors to shape data. Using an approach designed for manifold-valued data for shape spaces defined by the action of a group of diffeomorphisms raises several theoretical and computational difficulties. Are notably needed: a finite-dimensional set of diffeomorphisms with a Riemannian manifold structure; stable and efficient numerical schemes to compute Riemannian exponential and parallel transport operators on this manifold as no closed-form expressions are available; accelerated algorithmic convergence to cope with an hundreds of times larger dimensionality. To this end, we propose here:

• to formulate a generative non-linear mixed-effects model for a finite-dimensional set of diffeomorphisms defined by control points, to show that this set is stable under exp-parallelization, and to use an efficient numerical scheme for parallel transport; • to introduce an adapted MCMC-SAEM algorithm with an adaptive block sampling of the latent variables, a specific sampling strategy for shape parameters based on random local displacements of the shape contours, and a vanishing tempering of the target log-likelihood; • to validate our method on 2D simulated data and a large dataset of 3D brain structures in the context of Alzheimer's disease progression, and to illustrate the potential of our method for classifying spatiotemporal patterns, e.g. to discriminate pathological versus normal trajectories of ageing. All in one, the proposed method estimates an average spatiotemporal trajectory of shape changes from a longitudinal dataset, together with distributions of space-shifts, time-shifts and acceleration factors describing the variability in shape, onset and pace of shape changes respectively.

Deformation model 2.1. The manifold of diffeomorphisms D c0

We follow the approach taken in [START_REF] Allassonnière | Construction of bayesian deformable models via a stochastic approximation algorithm: a convergence study[END_REF] built on the principles of the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework [START_REF] Miller | Group actions, homeomorphisms, and matching: A general framework[END_REF]. We note d ∈ {2, 3} the dimension of the ambient space. We choose k a Gaussian kernel of width σ ∈ R * + and c a set of n cp ∈ N "control" points c = (c 1 , ..., c ncp ) of the ambient space R d . For any set of "momentum" vectors m = (m 1 , ..., m ncn ), we define the "velocity" vector field v : R d → R d as the convolution v(x) = ncp k=1 k(c k , x) m k for any point x of the ambient space R d . From initial sets of n cp control points c 0 and corresponding momenta m 0 , we obtain the trajectories t → (c t , m t ) by integrating the Hamiltonian equations:

ċ = K c m ; ṁ = - 1 2 ∇ c m T K c m (1)
where K c is the n cp ×n cp "kernel" matrix [k(c i , c j )] ij , ∇(.) the gradient operator, and (.) T the matrix transposition.

Those trajectories prescribe the trajectory t → v t in the space of velocity fields. The integration along such a path from the identity generates a flow of diffeomorphisms t → φ t of the ambient space [5]. We can now define:

D c0 = φ 1 ; ∂ t φ t = v t • φ t , φ 0 = Id, v t = Conv(c t , m t ) ( ċt , ṁt ) = Ham(c t , m t ), m 0 ∈ R dncp (2)
where Conv(., .) and Ham(., .) are compact notations for the convolution operator and the Hamiltonian equations (1) respectively. D c0 has the structure of a manifold of finite dimension, where the metric at the tangent space T φ1 D c0 is given by K -1 c1 . It is shown in [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF] that the proposed paths t → φ t are the paths of minimum deformation energy, and are therefore the geodesics of D c0 . These geodesics are fully determined by an initial set of momenta m 0 .

Then, any point x ∈ R d of the ambient space follows the trajectory t → φ t (x). Such trajectories are used to deform any point cloud or mesh embedded in the ambient space, defining a diffeomorphic deformation of the shape. Formally, this operation defines a shape space S c0,y0 as the orbit of a reference shape y 0 under the action of D c0 . The manifold of diffeomorphisms D c0 is used as a proxy to manipulate shapes: all computations are performed in D c0 or more concretely on a finite set of control points and momentum vectors, and applied back to the template shape y 0 to obtain a result in S c0,y0 .

Riemmanian exponentials on D c0

For any set of control points c 0 , we define the exponential operator Exp c0 :

m 0 ∈ R dncp → φ 1 ∈ D c0 . Note that D c0 = Exp c0 (m 0 ) | m 0 ∈ R dncp .
The following proposition ensures the stability of D c0 by the exponential operator, i.e. that the control points obtained by applying successive compatible exponential maps with arbitrary momenta are reachable by an unique integration of the Hamiltonian equations from c 0 : Proposition. Let c 0 be a set of control points.

∀φ 1 ∈ D c0 , ∀w momenta, we have Exp φ1(c0) (w) ∈ D c0 . Proof. We note φ 1 = Exp φ1(c0) (w) ∈ D φ1(c0) and c 1 = φ 1 • φ 1 (c 0 ). By construction, there exist two paths t → ϕ t in D c0 and s → ϕ s in D φ1(c0) such that ϕ 1 • ϕ 1 (c 0 ) = c 1 .
Therefore there exist a diffeomorphic path u → ψ u such that ψ(c 0 ) = c 1 . Concluding with [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF], the path u → ψ u of minimum energy exists, and is written u → Exp c0 (u.m 0 ) for some m 0 ∈ R dncp .

As a physical interpretation might be given to the integration time t when building a statistical model, we introduce the notation Exp c0,t0,t :

m 0 ∈ R dncp → φ t ∈ D c0
where φ t is obtained by integrating from t = t 0 . Note that Exp c0 = Exp c0,0,1 .

On the considered manifold D c0 , computing exponentials -i.e. geodesics -therefore consists in integrating ordinary differential equations. This operation is direct and computationally tractable. The top line on Figure 1 plots a geodesic γ : t → φ t applied to the top-left shape y 0 .

Parallel transport and exp-parallels on D c0

In [START_REF] Schiratti | Learning spatiotemporal trajectories from manifoldvalued longitudinal data[END_REF] is introduced the exp-parallelism, which is a generalization of Euclidian parallelism to geodesicallycomplete manifolds. It relies on the Riemmanian parallel transport operator, which we propose to compute using the fanning scheme [START_REF] Louis | Parallel transport in shape analysis: a scalable numerical scheme[END_REF]. This numerical scheme only requires the exponential operator to approximate the parallel transport along a geodesic, with proved convergence.

We note P c0,m0,t0,t : R dncp → R dncp the parallel transport operator, which transports any momenta w along the geodesic γ : t → φ t = Exp c0,t0,t (m 0 ) from φ t0 to φ t . For any c 0 , m 0 , w and t 0 , we can now define the curve:

t → η c0,m0,t0,t (w) = Exp γ(t)(c0) P c0,m0,t0,t (w) . (3)
This curve, that we will call exp-parallel to γ, is welldefined on the manifold D c0 , according to the proposition of Section 2.2. Figure 1 illustrates the whole procedure. From the top-left shape, the computational scheme is as follow: integrate the Hamiltonian equations to obtain the control points c(t) (red crosses) and momenta m(t) (bold blue arrows); compute the associated velocity fields v t (light blue arrows); compute the flow γ : t → φ t (shape progression); transport the momenta w along γ (red arrows); compute the exp-parallel curve η by repeating the three first steps along the transported momenta.

Statistical model

For each individual 1 ≤ i ≤ N are available n i longitudinal shape measurements y = (y i,j ) 1≤j≤ni and associated times (t i,j ) 1≤j≤ni .

The generative statistical model

Let c 0 a set of control points and m 0 associated momenta. We call γ the geodesic t → Exp c0,t0,t (m 0 ) of D c0 . Let y 0 be a template mesh shape embedded in the ambient space. For a subject i, the observed longitudinal shape measurements y i,0 , ..., y i,ni are modeled as sample points at times ψ i (t i,j ) of an exp-parallel curve t → η c0,m0,t0,t (w i ) to this geodesic γ, plus additional noise i,j :

y i,j = η c0,m0,t0,ψi(ti,j ) (w i ) • y 0 + i,j . (4)
The time warp function ψ i and the space-shift momenta w i respectively encode for the individual time and space variability. The time-warp is defined as an affine reparametrization of the reference time t: ψ i (t) = α i (tt 0τ i ) + t 0 where the individual time-shift τ i ∈ R allows an interindividual variability in the stage of evolution, and the individual acceleration factor α i ∈ R * + a variability in the pace of evolution. For convenience, we write α i = exp(ξ i ).

In the spirit of Independent Component Analysis [START_REF] Cootes | Active appearance models[END_REF], the space-shift momenta w i is modeled as the linear combination of n s sources, gathered in the n cp × n s matrix A:

w i = A m ⊥ 0 s i .
Before computing this superposition, each column c l (A) of A has been projected on the hyperplane m ⊥ 0 for the metric K c0 , ensuring the orthogonality between m 0 and w i . As argued in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifoldvalued longitudinal data[END_REF], this orthogonality is fundamental for the identifiability of the model. Without this constraint, the projection of the space shifts (w i ) i on m 0 could be confounded with the acceleration factors (α i ) i .

Mixed-effects formulation

We use either the current [START_REF] Vaillant | Surface matching via currents[END_REF] or varifold [7] noise model for the residuals i,j , allowing our method to work with input meshes without any point correspondence. In this setting, we note i,j iid ∼ N (0, σ 2 ). The other previously introduced variables are modeled as random effects z, with:

y 0 ∼ N (y 0 , σ 2 y ), c 0 ∼ N (c 0 , σ 2 c ), m 0 ∼ N (m 0 , σ 2 m ), A ∼ N (A, σ 2 A ), τ i iid ∼ N (0, σ 2 τ ), ξ i iid ∼ N (0, σ 2 ξ ), s i iid ∼ N (0, 1).
We define θ = (y 0 , c 0 , m 0 , A, t 0 , σ 2 τ , σ 2 ξ , σ 2 ) the fixed effects i.e. the parameters of the model. The remaining variance parameters σ 2 y , σ 2 c , σ 2 m and σ 2 A can be chosen arbitrarily small. Standard conjugate distributions are chosen as Bayesian priors on the model parameters:

y 0 ∼ N (y 0 , ς 2 O ), c 0 ∼ N (c 0 , ς 2 c ), m 0 ∼ N (m 0 , ς 2 m ), A ∼ N (A, ς 2 A ), t 0 ∼ N (t 0 , ς 2 t ), σ 2 τ ∼ IG(m τ , σ 2 τ,0 ), σ 2 ξ ∼ IG(m ξ , σ 2 ξ,0 ), and σ 2 ∼ IG(m , σ 2
,0 ) with inverse-gamma distributions on variance parameters. Those priors ensure the existence of the maximum a posteriori (MAP) estimator. In practice, they regularize and guide the estimation procedure.

The proposed model belongs to the curved exponential family (see supplementary material, which gives the complete log-likelihood). In this setting, the algorithm introduced in the following section has a proved convergence.

We have defined a distribution of trajectories that could be noted t → y(t) = f θ,t (z) where z is a random variable following a normal distribution. We call t → f θ,t (E[z]) the average trajectory, which may not be equal to the expected trajectory t → E[f θ,t (z)] in the general non-linear case.

Estimation

The MCMC-SAEM algorithm

The Expectation Maximization (EM) algorithm [11] allows to estimate the parameters of a mixed-effects model with latent variables, here the random effects z. It alternates between an expectation (E) step and a maximization (M) one. The E step is intractable in our case, due to the non-linearity of the model. In [10] is introduced and proved a stochastic approximation of the EM algorithm, where the E step is replaced by a simulation (S) step followed by an approximation (A) one. The S step requires to sample q(z|y, θ k ), which is also intractable in our case. In the case of curved exponential models, the authors in [2] show that the convergence holds if the S step is replaced by a single transition of an ergodic Monte-Carlo Markov Chain (MCMC) whose stationary distribution is q(z|y, θ k ). This global algorithm is called the Monte-Carlo Markov Chain Stochastic Approximation Expectation-Maximization (MCMC-SAEM), and is exploited in this paper to compute the MAP estimate of the model parameters θ map = max θ q(y, z|θ) dz.

The adaptative block sampler

We use a block formulation of the Metropolis-Hasting within Gibbs (MHwG) sampler in the S-MCMC step. The latent variables z are decomposed into n b natural blocks:

z = y 0 , c 0 , m 0 , [c l (A)] l , [τ i , ξ i , s i ] i .
Those blocks have highly heterogeneous sizes, e.g. a single scalar for τ i versus possibly thousands for y 0 , for which we introduce a specific proposal distribution in Section 4.3.

For all the other blocks, we use a symmetric random walk MHwG sampler with normal proposal distributions of the form N (0, σ 2 b Id) to perturb the current block state z k b . In order to achieve reasonable acceptance rates ar i.e. around ar = 30% [36], the proposal standard deviations σ b are dynamically adapted every n adapt iterations by measuring the mean acceptance rates ar over the last n detect iterations, and applying, for any b:

σ b ← σ b + 1 k δ ar -ar (1 -ar )1 ar≥ar + ar 1 ar<ar (5)
with δ > 0.5. Inspired by [3], this dynamic adaptation is performed with a geometrically decreasing step-size k -δ , ensuring the vanishing property of the adaptation scheme and the convergence of the whole algorithm [2, 3]. It proved very efficient in practice with n adapt = n detect = 10 and δ = 0.51, for any kind of data.

Efficient sampling of smooth template shapes

The first block z 1 = y 0 i.e. the coordinates of the points of the template mesh, is of very high dimension: naively sampling over each scalar value of its numerical description would result both in unnatural distorted shapes and a daunting computational burden.

We propose to take advantage of the geometrical nature of y 0 and leverage the framework introduced in Section 2 by perturbing the current block state z k 1 with a small displacement field v, obtained by the convolution of random momenta on a pre-selected set of control points. This proposal distribution can be seen as a normal distribution N (0, σ 2 1 D T D) where σ 2 1 is the variance associated with the random momenta, and D the convolution matrix. In practice, dynamically adapting the proposal variance σ 2 1 and selecting regularly-spaced shape points as control points proved efficient.

Tempering

The MCMC-SAEM is proved convergent toward a local maximum of θ → q(y, z|θ) dz. In practice, the dimensionality of the energetic landscape q(y, z|θ) and the presence of multiple local maxima can make the estimation procedure sensitive to initialization. Inspired by the globallyconvergent simulated annealing algorithm, [25] proposes to carry out the optimization procedure in a smoothed version of the original landscape q T (y, z|θ). The temperature parameter T controls this smoothing, and should decrease from large values to 1, for which q T = q.

We propose to introduce such a temperature parameter only for the population variables z pop . The tempered version of the complete log-likelihood is given as supplementary material. In our experiments, the chosen temperature sequence T k remains constant at first, and then geometrically decreases to unity. Implementing this "tempering" feature had a dramatic impact on the required number of iterations before convergence, and greatly improved the robustness of the whole procedure. Note that the theoretical convergence properties of the MCMC-SAEM are not degraded, since the tempered phase of the algorithm can be seen as an initializing heuristic, and may actually be improved.

Algorithm 1: Estimation of the longitudinal deformations model with the MCMC-SAEM.

Code publicly available at: www.deformetrica.org. input : Longitudinal dataset of shapes y = (y i,j ) i,j . Initial parameters θ 0 and latent variables z 0 . Geometrically decreasing sequence of step-sizes ρ k . output: Estimation of θ map . Samples (z s ) s approximately distributed following q(z | y, θ map ). Initialization: set k = 0 and S 0 = S(z 0 ). repeat Simulation: foreach block of latent variables

z b do Draw a candidate z c b ∼ p b (.|z k b ). Set z c = (z k+1 1 , ..., z k+1 b-1 , z c b , z k b+1 , ..., z k n b ). Compute the geodesic γ : t → Exp c0,t0,t (m 0 ). ∀i, compute w i = A ⊥ m0 s i . ∀i, compute w : t → P c0,m0,t0,t (w i ). ∀i, j, compute Exp γ[ψi(ti,j )](c0) w[ψ i (t i,j )] .

Compute the acceptation ratio

ω = min 1, q(z c |y,θ k ) q(z k |y,θ k ) . if u ∼ U(0, 1) < ω then z k+1 b ← z c b else z k+1 b ← z k b . end Stochastic approx.: S k+1 ← S k + ρ k S(z k+1 ) -S k . Maximization: θ k+1 ← θ (S k+1

Sufficient statistics and maximization step

Exhibiting the sufficient statistics S 1 = y 0 , S 2 = c 0 , S 3 = m 0 , S 4 = A, S 5 = i t 0 +τ i , S 6 = i (t 0 +τ i ) 2 , S 7 = i ξ 2 i and S 8 = i j y i,jη c0,m0,t0,ψi(ti,j ) (w i ) • y 0 2 , the update of the model parameters θ ← θ in the M step of the MCMC-SAEM can be derived in closed-form. The explicit expressions are given as supplementary material.

Experiments

Validation with simulated data in R 2

Convergence study. To validate the estimation procedure, we first generate synthetic data directly from the model without additional noise. Our choice of reference geodesic γ is plotted on top line of the previously introduced Figure 1: the template y 0 is the top central shape, the chosen five control points c 0 are the red crosses, and the momenta m 0 the bold blue arrow. Those parameters are completed with t 0 = 70, σ τ = 1, σ ξ = 0.1. With n s = 4 independent components, we simulate N = 100 individual trajectories and sample n i i = 5 observations from each.

The algorithm is run ten times. Figure 2 plots the evolution of the error on the parameters along the estimation procedure in log scale. Each color corresponds to a different run: the algorithm converges to the same point each time, as it is confirmed by the small variances on the residual errors indicated in Table 1. Those residual errors come from the finite number of observations of the generated dataset

∆v 0 2 |∆t 0 | |∆σ τ | |∆σ ξ | |∆σ | ∆v i 2 i |∆ξ i | i |∆τ i | i
1.43 ±5.6% 0.89 ±0.7% 0.19 ±2.7% 0.029 ±13.2% 0.017 ±7.6% 0.11 ±0.1% 2.47 ±1.7% 0.022 ±6.7% 0.19 ±0.8%

Table 1: Absolute residual errors on the estimated parameters and associated relative standard deviations across the 10 runs. Are noted v 0 = Conv(c 0 , m 0 ) and v i = Conv(c 0 , w i ). The operator . i indicates an average over the index i. Residuals are satisfyingly small, as it can be seen for |∆t 0 | for instance when compared with the time-span max|∆t ij | = 4. The low standard deviations suggest that the stochastic estimation procedure is stable and reproduces very similar results at each run.

Figure 3: Estimated mean progression (bottom line in bold), and three reconstructed individual scenarii (top lines). Input data is plotted in red in the relevant frames, demonstrating the reconstruction ability of the estimated model. Our method is able to disentangle the variability in shape, starting time of the arm movement and speed.

and the Bayesian priors, but are satisfyingly small, as qualitatively confirmed by Figure 3. The estimated mean trajectory, in bold, matches the true one, given by the top line of Figure 1. Figure 3 also illustrates the ability of our method to reconstruct continuous individual trajectories.

Personalizing the model to unseen data. Once a model has been learned i.e. the parameters θ map have been estimated, it can easily be personalized to the observations y new of a new subject by maximizing q(y new , z new | θ map ) for the low-dimensional latent variables z new . We implemented this maximization procedure with the Powell's method [35], and evaluated it by registering the simulated trajectories to the true model. Table 2 gathers the results for the previouslyintroduced dataset with n i i = 5 observations per subject, and extended ones with n i i = 7 and 9. The parameters are satisfyingly estimated in all configurations: the recon- 1, Figure 3). The acceleration factor is the most difficult parameter to estimate with small observation windows of the individual trajectories; at least two observations are needed to obtain a good estimate.

Experience |∆σ | |∆s i | i |∆ξ i | i |∆τ i | i n i i = 5 0.

Hippocampal atrophy in Alzheimer's disease

Longitudinal deformations model on MCIc subjects. We extract the T1-weighted magnetic resonance imaging measurements of N = 100 subjects from the ADNI database, with n i i = 7.6 datapoints on average. Those subjects present mild cognitive impairements, and are eventually diagnosed with Alzheimer's disease (MCI converters, noted MCIc). In a pre-processing phase, the 3D images are affinely aligned and the segmentations of the right-hemisphere hippocampi are transformed into a surface meshes. Each affine transformation is then applied to the corresponding mesh, before rigid alignement of follow-up meshes on the baseline one. The hippocampus is a subcortical brain structure which plays a central role in memory, and experiences atrophy during the development of Alzheimer's disease. We initialize the geodesic population parameters y 0 , c 0 , m 0 with a geodesic regression [15,[START_REF] Fletcher | Geodesic regression and the theory of least squares on riemannian manifolds[END_REF] performed on a single subject. The reference time t 0 is initialized to the mean of the observation times (t i,j ) i,j and σ 2 τ to the corresponding variance. We choose to estimate n s = 4 independent components and initialize the corresponding matrix A to zero, so as the individual latent variables τ i , ξ i , s i . After 10,000 iterations, the parameter estimates stabilized. Figure 4 plots the estimated mean progression, which exhibits a complex spatiotemporal atrophy pattern during disease progression: a pinching effect at the "joint" between the head and the body, combined with a specific atrophy of the medial part of the tail. Figure 5 plots an independent component, which is orthogonal to the mean progression by construction. This component seems to account for the inter-subject variability in the relative size of the hippocampus head compared to its tail.

We further examine the correlation between individual parameters and several patients characteristics. Figure 6 exhibits the strong correlation between the estimated individual time-shifts τ i and the age of diagnostic t diag i , suggesting that the hippocampal atrophy correlates well with the cognitive symptoms. The few outliers above the regression line Classifying pathological trajectories vs. normal ageing.

We processed another hundred of individuals from the ADNI database ( n i i = 7.37), choosing this time control subjects (CN). We form two balanced datasets, each containing 50 MCIc and 50 CN. We learn two distinct longitu-dinal deformations model on the training MCIc (N = 50, n i i = 8.14) and CN (N = 50, n i i = 8.08) subjects. We personalize both models to all the 200 subjects, and use the scaled and normalized differences z MCIc i z CN i as feature vectors of dimension 6, on which a list of standard classifiers are trained and tested to predict the label in {MCIc, CN}. For several number of observations per test subject n i i configurations, we compute confidence intervals by bootstraping the test set. Table 3 compares the results with a naive approach, using as unique feature the slope of individually-fitted linear regressions of the hippocampal volume with age. Classifiers performed consistently better with the features extracted from the longitudinal deformations model, even with a single observation. The classification performance increases with the number of available observations per subject. Interestingly, from n i i = 4 pooling the shape and volume features yields an improved performance, suggesting complementarity.

Conclusion

We proposed a hierarchical model on a manifold of diffeomorphisms estimating the spatiotemporal distribution of longitudinal shape data. The observed shape trajectories are represented as individual variations of a group-average, which can be seen as the mean progression of the population. Both spatial and temporal variability are estimated directly from the data, allowing the use of unaligned sequences. This feature is key for applications where no objective temporal markers are available, as it is the case for Alzheimer's disease progression for instance, whose onset age and pace of progression vary among individuals. Our model builds on the principles of a generic longitudinal modeling for manifold-valued data [START_REF] Schiratti | Learning spatiotemporal trajectories from manifoldvalued longitudinal data[END_REF]. We provided a coherent theoretical framework for its application to shape data, along with the needed algorithmic solutions for parallel transport and sampling on our specific manifold. We estimated our model with the MCMC-SAEM algorithm both with simulated and real data. The simulation experiments confirmed the ability of the proposed algorithm to retrieve the optimal parameters in realistic scenarii. The application to medical imaging data, namely segmented hippocampi brain structures of Alzheimer's diseased patients, delivered results coherent with medical knowledge, and provides more detailed insights into the complex atrophy pattern of the hippocampus and its variability across patients. In future work, the proposed method will be leveraged for automatic diagnosis and prognosis purposes. Further investigations are also needed to evaluate the algorithm convergence with respect to the number of individual samples.

Figure 1 :

 1 Figure 1: Samples from a geodesic γ (top) and an expparallelized curve η (bottom) on S c0,y0 . Parameters encoding the geodesic are the blue momenta attached to control points and plotted together with the associated velocity fields. Momenta in red are parallel transported along the geodesic and define a deformation mapping each frame of the geodesic to a frame of η. Exp-parallelization allows to transport a shape trajectory from one geometry to another.

  ). Adaptation: if remainder(k + 1, n adapt ) = 0 then update the proposal variances (σ 2 b ) b with equation (5). Increment: set k ← k + 1. until convergence;

Figure 2 :

 2 Figure 2: Error on the population parameters along the estimation procedure, with logarithmic scales. The residual on the template shape y 0 is computed with the varifold metric.

Figure 4 :Figure 5 :

 45 Figure 4: Estimated mean progression of the right hippocampi. Successive ages: 69.3y, 71.8y (i.e. the template y 0 ), 74.3y, 76.8y, 79.3y, 81.8y, and 84.3y. The color map gives the norm of the velocity field v 0 on the meshes.

Figure 6 :Figure 7 :

 67 Figure 6: Comparison of the estimated individual timeshifts τ i versus the age of diagnostic t diag i . R 2 = 0.74.

Table 2 :

 2

		110	3.34%	37.0%	5.45%
	n i i = 7	0.095	2.98%	16.2%	3.86%
	n i i = 9	0.087	2.38%	11.9%	3.28%

Residual errors metrics for the longitudinal registration procedure, for three simulated datasets. The absolute residual error on σ is given, the other errors are given in percentage of the simulation standard deviation. struction error measured by |∆σ | remains as low as in the previous experiment (see

Table

  Distribution of the third source term s i,3 according to the categories {τ i ≤ -3}, {-3 < τ i < 3}, and {3 ≤ τ i }. Hippocampal atrophy seems to occurs later in subjects presenting a lower volume ratio of the hippocampus tail over the hippocampus head (p = 0.0049).

						No APOE4 allele
	10					One APOE4 allele
						Two APOE4 alleles
	5
	0	55 (← early onset) 60	65	70 Onset age t 0 + τ i 75	80	85 (late onset →) 90
	Figure 8: Distributions of time-shifts τ i according to the
	number of APOE4 alleles. Hippocampal atrophy occurs
	earlier in carriers of 1 or 2 alleles (p = 0.017 and 0.015).
	10.0					Early onset
	7.5					Average onset Late onset
	5.0
	2.5
	-3 0.0	-2	-1 Source parameter s i,3 0 1		2	3
	Figure 9: n i i Shape features Naive feature		All features
	1	71% ±4.5 [lr]	50% ±5.0 [nb]		58% ±5.0 [lr]
	2	77% ±4.3 [lr]	58% ±4.9 [5nn] 68% ±4.7 [dt]
	4	79% ±4.1 [svm] 67% ±4.7 [5nn] 80% ±4.0 [lr]
	5	77% ±4.2 [nn]	77% ±4.2 [lr]		82% ±3.8 [nb]
	6.86	83% ±3.7 [lr]	80% ±4.0 [lr]		86% ±3.4 [lr]
	Table 3: Mean classification scores and associated standard
	deviations, computed on 10,000 bootstrap samples from the
	test dataset. Across all tested classifiers (sklearn default
	hyperparameters), only the best performing one is reported
	in each case: [lr] logistic regression, [nb] naive Bayes,
	[5nn] 5 nearest neighbors, [dt] decision tree, [nn] neural net-
	work, [svm] linear support vector machine.

Acknowledgments. This work has been partly funded by the European Research Council with grant 678304, European Union's Horizon 2020 research and innovation program with grant 666992, and the program Investissements d'avenir ANR-10-IAIHU-06.

Appendix: supplementary material

We introduce the onset age individual random variable t i = t 0 + τ i ∼ N (t 0 , σ 2 τ ) instead of the time shift τ i . The obtained hierarchical model is equivalent to the one presented in Section 3, with unchanged parameters θ = (y 0 , c 0 , m 0 , A, t 0 , σ 2 τ , σ 2 ξ , σ 2 ) and equivalent random effects z = (z pop , z 1 , ..., z N ), where z pop = (y 0 , c 0 , m 0 , A) and ∀i ∈ 1, N , z i = (t i , ξ i , s i ). The complete log-likelihood writes:

where the densities q(y i,j |z, θ), q(z pop |θ), q(z i |θ) and q(θ) are given, up to an additive constant, by:

-2 log q(y i,j |z, θ)

noting Λ the dimension of the space where the residual y i,jη c0,m0,t0,ψi(ti,j ) (w i ) • y 0 2 is computed, and |y 0 |, |c 0 |, |m 0 | and |A| the total dimension of y 0 , c 0 , m 0 and A respectively. We chose either the current [START_REF] Vaillant | Surface matching via currents[END_REF] or the varifold [7] norm for the residuals.

Noticing the identity η c0,m0,t0,ψi(ti,j ) = η c0,m0,0,ψi(ti,j )-t0 , the complete log-likelihood can be decomposed into log q(y, z, θ) = S(y, z), Φ(θ) Id -Ψ(θ) i.e. the proposed mixed-effects model belongs the curved exponential family. In this setting, the MCMC-SAEM algorithm presented in Section 4 has a proved convergence.

Exhibiting the sufficient statistics S 1 = y 0 , S 2 = c 0 , S 3 = m 0 , S 4 = A, S 5 = i t i , S 6 = i t 2 i , S 7 = i ξ 2 i and S 8 = i j y i,jη c0,m0,t0,ψi(ti,j ) (w i ) • y 0 2 (see Section 4.5), the update of the model parameters θ ← θ in the M step of the MCMC-SAEM algorithm can be derived in closed form:

The intricate update of the parameters t 0 ← t 0 and σ 2 τ ← σ 2 τ can be solved by iterative replacement. Similarly to Equation 6, the tempered complete log-likelihood writes:

with:

-2 log q T (y i,j |z, θ)

Tempering can therefore be understood as an artificial increase of the variances σ 2 , σ 2 y , σ 2 c , σ 2 m and σ 2 A when computing the associated acceptation ratios in the S-MCMC step of the algorithm. This intuition is well-explained in [25].


Ewout Martens

Georges Gielen
email: emartens|gielen@esat.kuleuven.ac.be

Time-Domain Simulation of Sampled Weakly Nonlinear Systems Using Analytical Integration and Orthogonal Polynomial Series

This paper presents a novel method for simulation of sampled systems with weakly nonlinear behavior. These systems can be characterized by adding weakly nonlinear terms to the linear state-space equations of the system resulting in an extended state-space model. Perturbation theory is used to split these equations in an ideal linear behavior and a non-ideal small perturbation. The linear equations are solved analytically which reduces simulation time compared to numerical evaluation. The solution of the perturbation equations is approximated by orthogonal polynomials.

This methodology not only reduces simulation time compared to traditional numerical simulations, but also deals naturally with clock jitter and the discontinuous behavior of sampled systems. An implementation of the methodology has been used to analyze systems including switched filters and continuous-time ∆Σ modulators.

Introduction

Modern integrated circuits frequently contain sampled subsystems. In mixed analog-digital systems, conversions between the analog and the digital domain are key operations. Further, a part of the processing algorithm can be implemented in the discrete domain, using for example switched-capacitor circuits.

Methodologies commonly used for systematic analysis and synthesis of analog circuits require timeefficient simulations to determine and optimize the performance [START_REF] Gielen | Computer-Aided Design of Analog and Mixed-Signal Integrated Circuits[END_REF]. Using these methods for sampled systems implies executing simulations which take into account major non-idealities like clock jitter, charge injection and weakly nonlinear behavior [START_REF] Laker | Design of Analog Integrated Circuits and Systems[END_REF]. Furthermore, signals encountered in switched circuits have inherently a discontinuous time behavior making it difficult to use simulators often used by analog designers like SPICE.

Simulation of sampled systems can be performed in time or frequency domain [START_REF] Yuan | Computer Methods for Switched Circuits[END_REF]. The goal of this work is to develop more efficient time-domain simulation of weakly nonlinear sampled systems which can be used as part of simulations of larger mixed-signal systems with a clock-based computational model. Applications are, for example, time-domain analysis of sensor interfaces with switched-capacitor filters and ∆Σ converters.

The approach presented in this paper is based on splitting up the nonlinear response between two sampling moments into a linear dominant part and a small nonlinear perturbation. These parts are obtained as solutions of the linear and perturbation equations respectively. Analytical integration is employed to provide an expression for the linear response. The perturbation, on the other hand, is approximated by a series of orthogonal polynomials. Combination of the two responses results in the time-efficient evaluation of weakly nonlinear behavior.

This approach has several advantages:

• Time-domain simulations simplify the modeling of systems which contain feedback via the digital domain. This is encountered, for example, in ∆Σ modulators. Note that a transient and not only a steady-state solution is calculated. Frequencydomain approaches, for example harmonic balance [START_REF] Kundert | Simulation of Nonlinear Circuits in the Frequency Domain[END_REF] or Volterra series [START_REF] Yuan | Distortion Analysis of Periodically Switched Nonlinear Circuits Using Time-Varying Volterra Series[END_REF][START_REF] Dobrovoln | Analysis and Compact Behavioral Modeling of Nonlinear Distortion in Analog Communication Circuits[END_REF], are more suited for calculating steady-state solutions.

• It is straightforward to take into account effects of clock jitter since it is not assumed that the width of the time interval is fixed. To include jitter effects, sampling points are generated by a random process. Other time-domain simulation techniques (e.g. [START_REF] Malcovati | Behavorial Modeling of Switched-Capacitor Sigma-Delta Modulators[END_REF]) which derive an analytical expression for the state values between two sampling moments, for example using the Z-transform, can only model effects of clock jitter by approximated linear behavior.

• Sampled nonlinear circuits can be simulated with general time-marching methods, for example using VHDL-AMS simulators or MATLAB/Simulink. High accuracy can be achieved at the expense of very time-consuming simulations. The proposed approach reduces the simulation time while still achieving a high accuracy.

• Thermal noise can easily be added to the timedomain model. Random noise processes can be converted to multivariate stochastic variables if they are filtered in the same way as the input signals and subsequently sampled.

This paper is organized as follows. Section 2 introduces the model for characterizing weakly nonlinear systems. It is explained how this set of equations is split up into a dominant and a perturbation part. The analytical integration method to solve the first set is explained in section 3. The approximation method for the perturbation equations is elaborated in section 4. In section 5 two examples of the simulation method are presented. Conclusions are drawn in section 6.

Characterizing weakly nonlinear systems

A n-th order linear system can be modeled by a set of state-space equations which describe the behavior in terms of n state variables [START_REF] Chen | The Circuits and Filters Handbook[END_REF]. To describe a weakly nonlinear system, these linear equations should be extended resulting in an extended state-space model.

Extended state-space model

Analog integrated circuits exhibiting weakly nonlinear behavior can be modeled using basic nonlinear elements like a nonlinear conductance, capacitance or transconductance which are described by a polynomial [START_REF] Wambacq | Distortion Analysis of Analog Integrated Circuits[END_REF]. For example, a weakly nonlinear transconductance of order N is characterized as follows:

i(t) = g m • v(t) + α 2 v 2 (t) + . . . + α N v N (t) , ( 1
)
where α 2 , . . . , α N are distortion coefficients. v(t) is the control variable of the nonlinearity. In general, the characteristic equation is a multivariate polynomial in multiple control variables. Figure 1 depicts an example of a nonlinear integrator and its macromodel. In this example, there are three control voltages (v C , v R and v 1 ) but only two

u(t) y(t) v 1 (t) C int R u(t) R v R (t) = y(t) v 1 (t) C int (v C ) R C in R out (v R ) g m (v 1 )v 1 C int (v C ) = C int 0 1 + α 2 v C + • • • + α N v N-1 C , g m (v 1 ) = g m 0 1 + β 2 v 1 + • • • + β N v N-1 1 , G out (v R ) = G out 0 1 + γ 2 v R + • • • + γ N v N-1 R .

Figure 1. Example of a weakly nonlinear model for an integrator of order N.

state variables. In general, control voltages are linear combinations of the state and input variables which can be represented by transformation matrices:

θ(t) = T q • q(t), (2a) ρ(t) = T u • u(t), ( 2 b )
with q(t) the (linear) state variables, u(t) input signals, θ(t) and ρ(t) control voltages and T q and T u transformation matrices. Note that products of control variables which are encountered in multivariate characteristic polynomials, can be converted to a sum of powers of state variables.

For the example of Figure 1, the state transformation matrix T q is defined as follows:

θ(t) =   v 1 (t) v C (t) v R (t)   =   1 0 0 1 1 -1   • v 1 (t) v C (t) = T q q(t), ( 3
)
where v 1 (t) and v C (t) are chosen as state variables.

Based on the linear state-space equations, the mathematical description of a weakly nonlinear system with output y(t) can be written as a set of extended state-space equations:

dq(t) dt + N ∑ k=2 P k dθ (k) (t) dt = A 1 q(t) + B 1 u(t) + N ∑ k=2 A k θ (k) (t) + N ∑ k=2 B k ρ (k) (t), (4a)
y(t) = C 1 q(t) + D 1 u(t) + N ∑ k=2 C k θ (k) (t) + N ∑ k=2 D k ρ (k) (t), (4b)
where the system is characterized by the parameters of the matrices P k , A k , B k , C k and D k . N is the distortion order of the weakly nonlinear system. The superscript (k) of a vector v (k) (t) indicates that each element of the vector should be raised to the power k. Using Kirchoff's laws, the extended state-space model for the example in Figure 1 can be derived:

dq(t) dt + N ∑ k=2 0 0 0 0 α k 0 dθ (k) (t) dt =   - G+g m 0 +G out 0 C in G out 0 C in g m 0 +G out 0 C int 0 - G out 0 C int 0   q(t) + -G C in 0 u(t) + N ∑ k=2   - g m 0 β k C in 0 - G out 0 γ k C in g m 0 β k C int 0 0 G out 0 γ k C int   θ (k) n (t). ( 5
)
All parameters are defined in Figure 1. The control variables θ(t) are defined by state transformation (3).

Solving the extended state-space model

To simulate a sampled system, the response of the extended state-space model is calculated over a time interval between two sampling points provided by a sampling clock, possible with jitter. This can be simplified by taking the weakly nonlinear character into account. In this case, the nonlinear response shall be a small deviation from the linear response. So, the total response of the weakly nonlinear system can be split up into the linear response and a small perturbation:

q(t) = q(t) + q(t), (6a) y(t) = ȳ(t) + ỹ(t). ( 6 b )
When the system is part of a cascade of systems, the input signal can also be split up in a linear and a perturbation signal:

u(t) = ū(t) + ũ(t), ( 6 c )
where the perturbation input ũ(t) coincides with the perturbation output of the previous system. Substituting these expressions into the extended state-space equation (4a), and simplifying, results in two sets of equations: the linear state-space equations

d q(t) dt = A 1 q(t) + B 1 ū(t), ( 7 )
and the perturbation equations

d q(t) dt = [I n + K(t)] -1 [A 1 q(t) + B 1 ũ(t) + e(t)] , (8a)
with

e(t) = N ∑ k=2 A k θ (k) (t) + B k ρ (k) (t) -K(t) d q(t) dt , (8b)
K(t) = N ∑ k=2 k P k (diag θ(t)) k-1 T q , ( 8 c )
where 'diag x' converts vector x to a diagonal matrix.

For each of these two sets of equations a different solution strategy is followed. Whereas for the linear set an analytical solution can be found, for the perturbation set an approximation method must be applied.

Solution of the linear equations

First, the linear state-space equations ( 7) are solved to determine the dominant response of the system. This can be done analytically based on the eigenstructure of the system matrix A 1 .

There are different formulations of the solution of the differential equation possible, for example using the exponential of A 1 [START_REF] Chen | The Circuits and Filters Handbook[END_REF]. The formulation used in this work, however, is chosen to make it possible both to reuse a part of the calculation throughout the simulation and to employ analytical integration.

To find the response of (7) at time t with initial conditions q(t 0 ), one has to sum all contributions of all eigenvalues. One can prove that a real eigenvalue λ k of the system matrix A 1 gives rise to the terms

Q k,m (∆t) m e λ k ∆t q(t 0 ) + R k,m t t 0 (t -τ) m m! e λ k (t-τ) ū(τ) dτ, (9a)
for m = 0, . . . , c k . The contribution of a complex eigenvalue pair α l ± j β l are the terms P r l,p cos (β l ∆t) -P i l,p sin (β l ∆t) (∆t) p e α l ∆t q(t 0 )

+ S r l,p t t 0 (t -τ) m m! e α l (t-τ) cos [β l (t -τ)] ū(τ) dτ -S i l,p t t 0 (t -τ) m m! e α l (t-τ) sin [β l (t -τ)] ū(τ) dτ, (9b)
for l = 0, . . . , d l . Here, c k + 1 and d l + 1 are the orders of the generalized eigenspace of the real and complex eigenvalues respectively. ∆t denotes the width of the time interval. The matrices Q k,m , R k,m , P r l,p , P i l,p , S r l,p and S i l,p contain coefficients depending on the system. They are calculated at the beginning of the simulation. So, although it may be tedious to accurately calculate these matrices, overall simulation time is not increased accordingly since the calculation must be done only once. Substituting the expression for the general solution into the state-space equations, and equalizing orthogonal terms, shows that the columns of these coefficient matrices are elements of the null space of

(A 1 -λ k I n ) c k -m+1 or (A 1 -(α l + j β k )I n ) d l -p+1 .
Analytical expressions are provided to calculate the integrals in the terms [START_REF] Yuan | Computer Methods for Switched Circuits[END_REF] for commonly used input signals. These analytical calculation of the integrals reduces simulation time compared to numerical evaluations. Together with the reuse of the coefficient matrices, this results in time-efficient simulations.

Furthermore, analytical integration has the advantage that input signals may be discontinuous. For example, to model initial charge of capacitors at the start of the time interval in a sampled system, Dirac impulses can be used as input signals [START_REF] Yuan | Computer Methods for Switched Circuits[END_REF].

Solution of the perturbation equations

Once the dominant behavior has been determined, the deviation of this response should be calculated by solving equations [START_REF] Yuan | Distortion Analysis of Periodically Switched Nonlinear Circuits Using Time-Varying Volterra Series[END_REF]. The total response is then calculated as the sum of the linear and perturbation parts.

The perturbation equations can be written in the standard form of differential equations:

d q(t) dt = F ( q(t), t) . ( 10
)
Furthermore, the time interval [t 0 , t 1 ] is normalized to a standard interval (in this paper [-1, 1]) using the transformation ξ = αt + β, (11) resulting in the normalized perturbation equations in

q(ξ) = q ξ-β α : d q(ξ) dξ = 1 α F q(ξ), ξ -β α . ( 12
)
Generally, an analytical solution cannot be found for these equations. Instead, the signals can be approximated over the entire time interval:

q(ξ) = M ∑ i=0 Qi p i (ξ), (13)
where the signals {p 0 (ξ), . . . , p M (ξ)} are the base signals of the approximation. In this work, orthogonal polynomials are chosen as base functions. This set has several advantages compared to other sets:

• Using polynomials simplifies function evaluations compared to other types of base functions. • One could also use Taylor series as polynomial approximation. The advantage of series of orthogonal polynomials is that they provide an approximation over the entire time interval. Taylor series, on the other hand, approximate the signal only at the beginning of the interval. The main point of interest, however, in the case of sampled systems is the end point of the interval.

In this work, Chebyshev polynomials are used as base functions.

To find the coefficients Qi of the polynomial series (13), one could substitute the approximation into the differential equations and equalize corresponding terms. This method has some disadvantages:

• Nonlinear behavior requires the calculation of powers of the approximated signal in terms of base functions of order 0 to M. Such expressions are nonlinear and complicated. Furthermore, they are approximations which introduce an error on the differential which is not correct anywhere in the interval. Especially in the end point of the interval, a large error can be the result. • All input signals, including the linear approximation, should be approximated by series of the base functions.

Another approach consists in writing down the differential equations in M collocation points. Together with the initial conditions, a nonlinear system is obtained of as many equations as with the method described above. Since this is not easier or more difficult to solve and it does not suffer of the aforementioned disadvantages, this method is our first choice to calculate the coefficients of the approximation.

There are M collocation points to be chosen within the standard interval. For Chebyshev polynomials, it is advantageous to select the points

ξ j = cos (j -1)π M -1 , j = 1, . . . , M, (14)
if M = 1. For M = 1, the end point is selected as collocation point (which is 1 in the standard interval). Now the differential equations are written down in each point ξ j . The left-hand side of (12) becomes

d q(ξ) dξ = M ∑ i=0 Qi dp i (ξ j ) dξ , (15)
with dp i (ξ j ) dξ the derivative of the base function of order i in ξ j . For Chebyshev polynomials and the choice of collocation points of (14), one can prove that

dp i (ξ j ) dξ =    i + 2i ∑ i-1 2 l=1 cos 2l(j-1)π M-1 i odd, 2i ∑ i 2 l=1 cos (2l-1)(j-1)π M-1 i even. (16)
These values should only be calculated at the beginning of the simulation, since the interval is always normalized to the same standard interval.

For the right-hand side of (12) the values of the normalized perturbation signals q(ξ) should be written as a function of the coefficients Qi . All coefficients can be lumped together within the coefficient vector Q:

Q = QT 0 QT 1 • • • QT M T . ( 17
)
The values of q(ξ) in the collocation point ξ j can be written as a matrix multiplication:

q(ξ j ) = P (ξ j ) • Q, (18a)
with

P (ξ) = diag n p 0 (ξ) • • • diag n p M (ξ) . ( 18b
)
The operator 'diag n ' is defined as the Kronecker product between the identity matrix and a scalar, vector or matrix (i.e. diag n φ = I n ⊗ φ). The function values of a Chebyshev polynomial in the points of ( 14) are

p i (ξ j ) = cos i(j -1)π M -1 . (19)
Using these expressions, the set of differential equations written down in the collocation points together with the initial conditions can be converted to a system of (nonlinear) algebraic equations:

V ( Q) =      P (ξ 1 ) . . . P (ξ M ) P (-1)      Q - 1 α        F P (ξ 1 ), ξ 1 -β α . . . F P (ξ M ), ξ M -β α 0        = 0, (20)
with P (ξ) defined similar to P (ξ) but with derivatives instead of function values of the base functions. An approximation of the perturbation coefficients can be found by solving the set of nonlinear algebraic equations using Newton-Rhaphson iterations. Since the system is weakly nonlinear, it should have an approximately linear behavior. So, a first-order approximation is with dV (0) d Q the jacobian of the vector V ( Q). Although more accuracy can be obtained by performing multiple Newton iterations, (21) is a good approximation for weakly nonlinear circuits. No significant improvement was found when two or three iterations were performed in the experiments described in section 5. Note that this means that the coefficients are the result of the solution of one set of linear equations.

Q ≈ - dV (0) d Q -1 • V (0), (21)
u(t) x(t) Ca C b Cc C d Cg Ce φ 1 φ 1 φ 2 φ 1 φ 1 φ 2 φ 2 φ 1 φ 2 φ 2
The complexity of this algorithm depends on how efficiently the jacobian in (21) can be calculated. It is straightforward to calculate dV (0) d Q once the jacobian of F( q) of ( 12) is evaluated in all collocation points. For this evaluation, one can determine an expression similar to [START_REF] Yuan | Distortion Analysis of Periodically Switched Nonlinear Circuits Using Time-Varying Volterra Series[END_REF], i.e. the calculation of the jacobians is simplified by storing terms common to the jacobian and the function values, especially the factor [I n + K( q)] -1 . Furthermore, parts of both functions should only be calculated once at the beginning of the simulation.

Experimental results

An implementation of the above algorithm has been developed to verify its accuracy and simulation time. The model has a discrete-time computational model which makes it straightforward to implement it as a SystemC module. tions 1 . One can conclude that, although high orders of the polynomial series are required for good approximation, simulations are about six times faster.

As second experiment, the model is used to calculate the SNDR degradation in a third-order singleloop ∆Σ A/D-converter with an oversampling ratio of 128. At each time step, the SystemC module calculates the new states of the loop filter. These values are the inputs of a quantizer module which produces a bit stream. Feedback is provided by a DAC by generation of pulses at sampling moments. The integrator model is shown in Figure 1. Figure 4 depicts the drop in SNDR as a function of the distortion coefficient of the integrating capacitor. Comparison is made with a time-marching simulation of a VHDL-AMS model 2 . Due to the feedback structure of a ∆Σ modulator and its tolerance to circuit non-idealities, lower-order approximations already result in a good accuracy. Also for this experiment, a speed-up factor of about six is achieved.

Conclusions

Time-efficient simulations of sampled systems are a great aid in the systematic analysis and synthesis of switched subsystems frequently used in analog and mixed-signal integrated circuits.

The simulation algorithm described in this paper is applicable for sampled weakly nonlinear circuits and achieves a reduction of the simulation time by two means. First, the dominant behavior is calculated analytically. Then, the deviation of the linear response is 1 The Runge-Kutta algorithm is included in the SystemC module to make a fair comparison. 2 The simulator used is AdvanceMS of Mentor Graphics R . calculated as a small perturbation. Since the perturbation should be small, it can be approximated by a finite polynomial series whose coefficients are calculated by solving a set of linear equations. Experimental results of a switched-capacitor filter and a continuous-time ∆Σ converter have shown that the proposed method is more time-efficient than time-marching algorithms used by, for example, MATLAB and VHDL-AMS simulators, while giving similar results.

Figure 2 .

 2 Figure 2. Fleischer-Laker switched-capacitor biquad used in the first experiment [5].

Distortion coefficient α 3 SNDRFigure 3 .

 33 Figure 3. Effect of nonlinear capacitors C b and C d in the circuit of Figure 2.

Figure 4 .

 4 Figure 4. SNDR degradation due to nonlinear distortion in ∆Σ modulator.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)

Theory & Applications, 50(8):1013-1024, Aug. 2003.