Raphaële Milan

Loïc Lagadec

Lilian Bossuet

Secured-by-design systems-on-chip : a MBSE approach

Keywords: Security, CAD Tool, Model-Based Engineering, Hardware, Systems-on-Chip

Security by Design (SbD) has gained increasing interest over the past decade. While iterative processes and legacy preservation aim to reduce costs and mitigate risks through continuity, SbD encourages a break in the way we do things with a simple idea: dealing with new threats, leading to new risks, requires a complete rethink of our design processes.

In embedded systems, security has been more or less left aside for a long time, with performance being the main objective. When security concerns emerged, the response was to adapt existing solutions with security patches. This is neither sustainable (to change from simple embedded systems to complex systems-onchip) nor simply effective. It is necessary to change the mindset, which will lead to new practices. But the central question is: "How can we put security at the heart of the design process?" The aim of this paper is to contribute to this reflection by providing a rapid prototyping environment (modeling and simulation-based systems engineering) for the hardware mechanisms responsible for the deployment of rights management services.

I. INTRODUCTION

Security matters. But sometimes time-to-market constraints outweigh security, and reuse is the way to reduce design time. When writing software, reusing pieces of code opens security gaps. There are two main reasons for this. The first is the misuse or misunderstanding of the existing code. A zeroday attack may have been reported and corrected, without the solution being incorporated into the reused code. Contextrelated problems (e.g., overflows) are also likely to occur. The second problem is that reuse is sometimes based on unsourced code. In this case, incorporating code from a library is easy, but security control becomes impossible. Threats can obviously come from malicious suppliers, but also simply from poorquality code (e.g. the C system() function calls binary files, which can be silently replaced by malicious files). A corrupted process that has access to shared storage or communication devices can compromise the entire system.

One solution is to use hypervisors. These hypervisors provide an interface layer between the software and the execution This work has been funded by French AID. platform. Hypervisors check requests against legitimacy rules to grant or deny access to devices. Incorporating hypervisors, at first, promotes a secure by design approach.

While security as a whole has become a major issue, Systems-on-Chip (SoC) security in particular is a hot topic, because such systems are everywhere: computers, cell phones, IoT devices, and Supervisory Control and Data Acquisition (SCADA) systems, making them critical infrastructures.

When moving from software to system-on-a-chip, the analysis of security gaps remains. But with added complexity. On the positive side, there are some useful mechanisms. Vendors offer hardware solutions to ensure minimal rights management at a low level. For example, ARM's TrustZone [START_REF] Arm | Arm security technology : Building a secure system using trustzone® technology[END_REF] divides the system into two worlds based on the NR bit value. The rich world, dedicated to common execution, is separated from the secure world, which has additional rights. Processes in the rich world have no access to resources in the secure world. These architectures enforce strict security measures right from the powering on of the system to prevent misuse and compromise. However, on the other hand, such security measures have been found to be vulnerable where security design practices are not considered or are poorly implemented, particularly at software and hardware stack boundaries. Additionally, there is a wide range of devices, resulting in heterogeneity and therefore increased complexity. Peripherals are generally reused from previous designs or come from unreliable third parties. Particularly when considering reconfigurable implementation, for which the synthesis and optimization phase is costly, it is tempting to use off-the-shelf versions from hardware libraries.

This paper proposes an approach for Systems-on-Chip prototyping ensuring key qualities like trustworthiness and safety while maintaining performances. We show that it is possible to deploy in SoCs a hypervisor-based solution that would mediate between software and hardware but also between hardware devices, including synthesized tasks. This additional layer provides a secure by-design infrastructure. It conforms to a set of requirements expressed as rules. In such architectures, each architecture element (processor, FPGA, memory, etc.) has a wrapper that controls its security, based on a distributed set of rights (authorization) dictionaries that specify who can perform what operation on what address range in what trust mode (secure, rich, etc.) on that element. Permissions can be specified statically at synthesis time, or dynamically, i.e., they can be modified between two runs of the platform. On the contrary, by exploiting these requirements through model-to-model (M2M) transformations, test cases can be generated, which challenge the infrastructure, in search of counterexamples for validation purposes. The M2M approach also leads to additional benefits, which are summarized in the perspective section.

This paper also proposes a fast prototyping solution which conforms to modeling and simulation-based systems engineering (MBSE). A conceptual model describes the architecture (architecture elements, element wrappers, tasks, task allocation to elements, etc.). This model supports instance creation, size change, and reshaping. It offers a useful exploration feature for those who want to design an SoC, as it comes with a discrete element-based simulation framework to score metrics. This framework can be used to simulate various architecture and rights management scenarios at the model level before they are implemented at the lower level.

The remainder of this paper is organized as follows. Section II introduces the systems on the chip domain, and then the related security issues. Section III presents our two contributions: security through rights management and SoC modeling, respectively. Section IV discusses SoC prototyping using FPGA board, with an emphasis on M2M transformations to alleviate manual, repetitive, and error-prone process.

II. BACKGROUND & RELATED WORK

SoCs are complex electronic devices. They integrate heterogeneous resources (processors, memories, devices, hardware accelerators, etc.) interconnected by a shared bus (e.g. AXI for Advanced eXtensible Interface).

Designing and exploiting such devices require specific tooling, referred to as Computer-Aided Design (CAD) tools, that aid in the creation, modification, analysis, or optimization of a design. CAD processes and software have a long history and are of the highest quality. They form a solid basis for SoC implementation. However, they suffer from rigidity. This problem has been identified since the early 2000s and efforts have been made to remedy the situation. MBSE (Model-Based System Engineering) offers a range of possibilities: Madeo [START_REF] Lagadec | Object-oriented meta tools for reconfigurable architectures[END_REF] provides a model along with a DSL to support exploration of the SoC domain space. The MARTE standard [START_REF]About the uml profile for marte specification, version 1.2[END_REF] (Modeling and Analysis of Real-Time and Embedded Systems) offers a methodology for SoC co-design [START_REF] Quadri | Targeting reconfigurable fpga based socs using the marte uml profile: from high abstraction levels to code generation[END_REF], with GASPARD [START_REF] Atitallah | Gaspard2: from marte to systemc simulation[END_REF] as the SoC co-design environment to move from high-level MARTE specifications to an executable platform. Teodorov [START_REF] Teodorov | Model-driven physical-design automation for fpgas: fast prototyping and legacy reuse[END_REF] has demonstrated regression-free inheritance refactoring while offering additional debugging facilities, These works bring new functionalities (exploration, debugging, automatic generation), while preserving legacy compatibility. This approach is suitable when targeting SoC security-by-design for two reasons. First, existing tools can be used to generate parts of the SoC that are secure (known IPs such as network-on-chip or custom IPs). Second, even if the complete SoC incorporates untrusted third-party IPs, new tools can address specific needs, including the addition of security mechanisms.

Although security requires a holistic approach, targeted contributions are based on certain assumptions. For example, Behani [START_REF] Benhani | The security of arm trustzone in a fpga-based soc[END_REF] demonstrated malicious tampering with AXI signals, and Gross [START_REF] Gross | Breaking trustzone memory isolation and secure boot through malicious hardware on a modern fpgasoc[END_REF] broke the TustZone mechanism, which had a devastating impact on the SoC as a whole. However, when it comes to focusing on incorporating untrusted third-party IPs into the SoC, not all of these IP cores have passed extensive integrity tests [START_REF] Banga | Trusted rtl: Trojan detection methodology in pre-silicon designs[END_REF], so no security assumption can be made. Still, both mechanisms (Bus, TrustZone) should be considered correctly implemented.

Mal-Sarkar [START_REF] Mal-Sarkar | Hardware trojan attacks in fpga devices: threat analysis and effective counter measures[END_REF] and Jacob [START_REF] Jacob | Compromising fpga socs using malicious hardware blocks[END_REF] demonstrated that thirdparty IPs can allow attacks such as hardware Trojans and malware that can be launched within any device using the compromised IP. A use case is to compromise a SoC running a secure update. Siddiqui [START_REF] Siddiqui | Pro-active policing and policy enforcement architecture for securing mpsocs[END_REF] proposes an additional layer of hardware security when a vulnerability is found and exploited, which focuses on mitigating hardware trojans.

III. CONTRIBUTIONS

The SoC must be immune to single failure points and must then contain an active response or mitigations for circumstances where a compromise occurs. The threat model consists of one or more IPs that are malicious or misprogrammed and attempt to perform transactions that violate the SoC access control policy. Our solution is based on wrappers as introduced by Siddiqui [START_REF] Siddiqui | Pro-active policing and policy enforcement architecture for securing mpsocs[END_REF]. Unlike IPs, wrappers appear as relatively simple white boxes, hence are supposed to be immune to external non-physical attacks. Wrappers grant/deny access to remote address ranges depending on the requester, its world, and the requested operation (R/W). Compared to [START_REF] Siddiqui | Pro-active policing and policy enforcement architecture for securing mpsocs[END_REF], wrappers are not bounded to filtering-out malicious IPs that alter the NS bit, but also cover broader use cases, as any device, including hard or soft IPs (implemented on FPGAs), is free to emit requests of his own. Such requests expose ID, trust mode, R/W mode, and target address; destination wrappers are responsible for reducing and ensuring system assets on the hardware, based on predefined rules.

Wrappers also have to initiate procedures if malicious activity is discovered, then in order to cope with new threats, they should be easily extendable to incorporate new functionalities. The deployment of new features (e.g. redundancy, reallocation, etc.) or policies (e.g. performing collaborative checks, preventing denial of service attacks) is beyond the scope of this article.

A. System-on-Chip executable models

Capturing SoCs as executable models allows fast prototyping of both architectures and protection mechanisms thanks to rich expressiveness and simulation. Such models are dual: on the one hand, they offer an abstraction with pure concepts, such as those expressed in UML (Figure 1), on the other hand, they provide an executable platform for efficient execution, verification, scenario elicitation, and change tracking for nonregression during refactoring. The automatic generation of code from the executable model and the deployment of this code directly as requirement conformance (within wrappers) and compatibility tests (test-driven development), make the process fast and reliable. Supporting a gradation of abstraction levels arbitrates between simulation accuracy and computation time (e.g. the different burst types as shown in Figure 2). The simplification of the bus protocol (direct addressing of nodes, MMUbased memory shifting, etc.) supports debugging of the node implementation (mixed behavioral level validation). On the contrary, serialized bursts totally reflect the binary information that transits over the bus.

Finally, the scenario captures both dimensions of the model and operates by invoking a simulator on the tasks. The scenario also supports the exploration of the domain space using a parametric architecture and a pseudorandom distribution of tasks. One given architecture is specified as an instance of the model. Listing 1 illustrates a piece of an action language (Smalltalk Pharo dialect [START_REF]Pharo: The immersive programming experience[END_REF]) code, which defines a simple architecture, made up of a processor and a memory, interconnected via a bus. During the architecture setup phase, the elements connect to the bus. This results in creating a wrapper that will be further in charge of rights management. The wrapper is the only interface between the nodes and the bus that they both reference. A wrapper has no ID but delegates all its incoming requests to its node.

A task requesting remote access delegates it to its wrapper. The wrapper generates a burst to be transmitted to the bus, then over the bus. The bus dispatches the burst to the proper destination wrapper. This wrapper is responsible for checking the rights (see Section III-B) before serving the request.

Figure 2 represents as an activity diagram the sequence of operations when receiving a burst. Three abstraction levels, ranging from high-level to serialized bursts, are considered. High-level bursts consist of a receiver, data and a relative address at the receiver, as well as a mode (R/W) and a trust mode (rich, secure, extended 1. . . n). In other words, these bursts contain the elements of the model (arbitrary complex data, real counterpart node, relative address instead of absolute address). This makes it possible to avoid numerous decoding operations, to easily visualize the content of the transmissions, and to be able to quickly verify the correct operation of the implemented mechanisms (rights, but also actual operations at the receiver). Deserializing a binary stream goes through slicing it into words based on the packet's size, and then using the first word to reconstruct a control packet. This packet carries all the structural information to rebuild a low-level burst. Then the following words are deserialized to rebuild the data fragments.

B. Security mechanism relying on rights management

When receiving the burst, the wrapper evaluates the security policy and delegates the operation to its node -be the operation legitimate. Figure 3 represents the sequence of operations when a write request is initiated: delegation of the platform on which the task is executed to the wrapper, invocation of the bus by the wrapper, which transmits the data to the appropriate destination wrapper. Then, two possibilities must be considered: either the operation is granted, in which case the target wrapper delegates to its device the execution of the operation, or it filters the request. In both cases, a return burst is emitted.

Rights are described as a dictionary whose associations link a tuple consisting of (requestor ID, mode, trust mode) to an address range. Authorization is granted if and only if the address range allows the task running in the mode (rich, secured, or an extended version of modes based on a set of security flags) to perform the operation (read/write). The absence of an entry in the dictionary means that there ∃!(task, mode, trustmode -→ @ ) =⇒ REJECT ∪{ ----------→ @ right≥requested }∩ ------→ @ requested = ------→ @ requested =⇒ OK is no authorization. When it receives a request, a wrapper compares it with the registered authorizations and either serves the request or returns a missing-authorization warning.

The algorithm 1 summarizes how the decision is made. Address ranges, denoted as -→ @ support set operators ∪, ∩, and ⊂.

Algorithm 2: Expressing the SR MEMORY2 rule: Secure and non-secure application memory pages are physically isolated from each other ∀(task, mode, trustmode), (task, mode, trustmode ′ ) ∃(task, mode, trustmode,

-→ @ A ) ∃(task, mode, trustmode ′ , -→ @ B ) =⇒ -→ @ A ∩ -→ @ B = ∅
Conforming to the security rules goes through deploying authorizations accordingly, as illustrated in Algorithm 2 for memory segregation. These rules are formulated as OCL [START_REF] Cabot | Object Constraint Language (OCL): A Definitive Guide[END_REF] code, the default language for expressing all types of (meta) model query, manipulation, and specification requirements. Listing 2 encodes the property of Algorithm 2.

Context Right Inv SR Memory2

Select (i | i key.task = self.key.task and r.key.trustMode <> self.key.trustMode and r.value = self.value ) -> isEmpty Listing 2: Encoding the SR MEMORY2 rule as OCL code

C. Simulation

Executable models can be easily simulated using a discreteevent engine, which represents the operation of a system as a list of events ordered by date [START_REF] Goldberg | Smalltalk-80: The Language and its Implementation[END_REF]. Each event occurs at a given time t and marks a state change in the system. Between two consecutive events, no change occurs in the system; therefore, the simulation time can jump directly to the timestamp of the next event and then resume.

The tasks and respective duration δ are known; a task is triggered as soon as all its predecessor tasks have been completed. In the context of this work, events are model activation: task start/stop, plus internal events (burst, dispatch, permission checks, etc.). Starting a task generates its end as a future event at t + δ, and triggers its associated action, modeled as a closure that takes the platform on which the task is running, as a parameter. Figure 3 provides an illustration in which the values δ are reflected by rectangles of respective height. The simulation loop ends when there are no more pending tasks. The task graph can be reshaped at will, on the fly, to add/remove dependencies (for example, adding "virtual" dependencies can be used to force sequential execution), modify task duration (including interactive change), define triggers to stop execution on specific events such as rights violation, and the simulation supports past snapshot restoration, to test multiple execution paths.

Rules implement rights policies in wrappers, but also serve as a basis for automatic test generation [START_REF] Verhaeghe | Usage of tests in an open-source community: A case study with pharo developers[END_REF] via scenarios, as illustrated by Listing 3. Tests are used to validate compliance with the policy. As an example, Listing 4 illustrates the compliance check of the OCL code of Listing 2, with isConsistent an atomic test ensuring non-overlapping memory addresses. Compliance-oriented scenarios are pseudogenerated and aim to find a counterexample that would identify the violation of a rule. Scenario parameters are generated randomly: topology, tasks, tasks' implementation (IP, process), events distribution (read/write), etc.Once generated, the scenario is saved to be replayed again.

Determinism-oriented scenarios are based on a new randomization step. They are derived from scenarios that do not violate any property. A random mix of abstraction levels (thanks to the M2M transformation) is used to run the same scenarioExample: param | arch proc ip1 ip2 bus mem data read pData | "param is (@W1, @W2, @W3, grants)" arch := TSArchitecture architectureExample: param last.

data := (Array new:2000 withAll:1).

"arbitrary addressing" ip1 write:(1 to:250) at: (param at:1). ip1 write:(1000 to:1250) at: (param at:2).

read := proc read: (TSAddressRange from: param at:1) +50 to: (param at:2) + 200).

"Simple example here" pData := (read collect:[:a| a * 2]). "arbitrary addressing" proc write: pData at: (param at:3). ˆarch

Listing 3: Parametric scenario definition scenario repeatedly. Examples of distinct abstraction levels are high-/low-level burst transactions, mock computation (XOR vs. AES), etc. In addition, data are moved into new address ranges, and all relative addresses are modified accordingly. Faulty behavior means getting different results for variants of the scenario. Such behaviors are detected through known result assertion violation, as illustrated in Listing 5. to: (param at:3)+150=(2 to: 300 by:2). self assert: (mem contents from: (param at:3) +51...

Listing 5: Scenario behavioral check IV. HARDWARE DESIGN

A. Wrappers

Wrappers serve as a firewall, filtering requests to IPs based on a permissions table, as shown in Figure 4. The experiments were conducted on a ZYBO (ZYnq BOard) platform [START_REF] Digilent | Zybo Zynq-7000 ARM/FPGA SoC Trainer Board[END_REF], an entry-level digital circuit and embedded software development platform built around Xilinx Z-7010 [START_REF] Xilinx | Zynq-7000 SoC Data Sheet: Overview[END_REF]. Vivado Design Fig. 4. Permissions table within a wrapper Suite [START_REF] Xilinx | Vivado[END_REF], from Xilinx, addresses productivity bottlenecks in system-level integration and implementation. It provides synthesis and analysis facilities of hardware description language (HDL) designs, with features for system-on-a-chip development and high-level synthesis. Considering a commercial platform for prototyping purposes brings additional complexityas some practical concerns must be addressed. As an example, Figure 5 illustrates that the AXI bus connecting the various IPs is separate from the one associated with the processor. In the ZYBO platform the processor is hardwired and connected to an AXI-3 bus. IPs are synthesized on the FPGA, as is their bus, with an AXI-4 version. Therefore, an AXI-3 to AXI-4 converter must be inserted into the design. Adding this module has no meaning from a conceptual point of view, but relates to low-level platform-specific constraints.

B. Design

A first proof-of-concept has been implemented. Wrappers have been added to the SoC nodes at design time, as shown in Figure 5. Experiments have shown an average increase of 0.34% in terms of LUTs, and 0.1% in terms of registers, without impact on frequency, over a set of five representative IPs. The design has been specified using the Vivado graphical user interface. In addition to being nonintuitive for nonspecialists (e.g. computer scientists), a significant part of our target audience, such interfaces degrade productivity, still reflect users' traditional way of doing things. Experiments demonstrated that the wrappers adequately mitigated unauthorized access.

Hardware validation goes through three steps. First, we need to make sure that the unprotected hardware is working properly and that the processor is accessing the devices. This validation is essential, as communication must pass through the AXI-3 to AXI-4 converter in both directions. This check is carried out by executing a C code that writes and then reads back values (i.e., a simplified Listing 5-like mechanism). The second step then consists of reproducing the same approach, at a lower level of abstraction, by using a TCL code, as illustrated in Listing 6 which describes a test comparing the values written and read according to a burst. Last, if securitybased filtering is added, the burst operation can be allowed or denied, as shown in Figure 5. For the sake of conciseness, the corresponding test code is not provided in this article.

The experiments showed a strict equivalence between the high-level executable model and the hardware implementation. However, adapting the hardware design demonstrated to be time-consuming, and hence preserves from design space exploration. Fortunately, model-based engineering (MBE) helps capture the expertise of hardware design while promoting a higher level of abstraction. Two key points strengthen this need: first, rapidly evolving security issues lead to short refactoring cycles, and second, the high modularity of chip systems requires point-topoint solutions combined with global composition. This work emphasizes the need to master several levels of abstraction. First, in modeling, we facilitate a detailed understanding of the interactions between the elements while preserving agility through transparent interoperability between these levels. MBE is part of the solution, and targeted development, debugging, and testing can take advantage of. When it comes to the application of rights policies, this leads to real productivity gains. Designing complex rights policies requires ease of specification, simulation facilities, and automatic property validation prior to moving to implementation.

In this work, the Vivado tool was used for demonstration purposes. It was used to design SoC nodes and wrappers.Then to incorporate parts of HDL code from libraries and interconnect them. Among these codes, we obviously find the wrapper code, which is secure by design. This hardware demonstrator has been proven to fully conform to the highlevel model and to detect and prevent unauthorized access to devices. The execution of models allowed us to expand the testing by quickly exploring hardware topologies along with the distribution of variable rights. This demonstrated that the combined approach is valuable and viable.

This work opens up many perspectives. Modern commercial EDA tools usually provide end-users with a framework for application-specific customization. This aims at easy interfacing between general-purpose programming languages and underlying circuit object models. Vivado takes advantage of TCL, which offers an M2M transformation back-end process to accelerate the synthesis of parametric SoC and allow defining scenarios. This work perspective also includes the interface between LiteX from Enjoy-digital [START_REF] Kermarrec | LiteX: an open-source SoC builder and library based on Migen Python DSL[END_REF], a SoC builder/IP library, and utilities to create SoCs and complete FPGA designs. In addition to being open source and BSD-approved, its unique feature is that all IP components are described using Migen Python's internal DSL, which simplifies its design in depth. These two perspectives claim for code generation as an enabler to either prototype or implement the model as HDL. The application on the processor requires a C code, which [START_REF]About the uml profile for marte specification, version 1.2[END_REF] will help generate. The work presented in this paper allows to quickly prototype mitigation hardware mechanisms, but also will, through M2M, support Vivado TCL, LiteX, and Marte as backends to port global design to many technologies.

The agile methodology encourages consideringsimple cases, before tackling realistic and complex cases. Accordingly, we will validate the refinement of high-level abstractions against hardware description language code, as well as the accuracy of functional analyses. Once done, relevant measurements will be available and flexible processes will be used to examine increasingly complex and realistic cases. The objective of measurement is productivity, i.e. the time difference between designing a complex system by hand and modeling the same system in our environment. Also is the ability to interrupt a simulation, or even to go back in time and modify a parameter to measure its impact. This refers both to the hardware aspect (number of cores, heterogeneity, rights, reconfiguration, cryptography facilities, etc.) and to attack scenarios; taking these into account ensures that our protection mechanisms are effective enough to prevent known threats across the platform as a whole, as well as to prevent future threats, due to their malleability and scalability. The changes we make require very little hardware development time (e.g. rewriting wrappers).

Fig. 1 .

 1 Fig. 1. Capturing the Security-oriented SoC domain as executable model

Fig. 2 .

 2 Fig. 2. Receiving high level, low level and serialized bursts

Fig. 3 .

 3 Fig. 3. Write operation sequence

Listing 4 :

 4 TSArchitecture architectureExample: param last. self assert: (arch allDevices reject:[:a| a wrapper rights isConsistent not]) isEmpty. Scenario compliance check behavioralCompliance: grantedParameters | arch | arch := TSArchitecture scenarioExample: grantedParameters. self assert: (mem contents from: (param at:3)

Fig. 5 .

 5 Fig. 5. Low level of architecture (3 IPs + Bus)

#

  Test: S00_AXI # Create a burst write transaction at # s00_axi_addr address create_hw_axi_txn w_s00_axi_addr [get_hw_axis $jtag_axi_master] -type write -address $s00_axi_addr -len 4 -data $wdata_2 -burst INCR # Create a burst read transaction at # s00_axi_addr address create_hw_axi_txn r_s00_axi_addr [get_hw_axis $jtag_axi_master] -type read -address $s00_axi_addr -len 4 -burst INCR # Initiate transactions run_hw_axi r_s00_axi_addr run_hw_axi w_s00_axi_addr run_hw_axi r_s00_axi_addr set rdata_tmp [get_property DATA [get_hw_axi_txn r_s00_axi_addr]] # Compare read data if { $rdata_tmp == $wdata_2 } { puts "Data comparison test pass for -S00_AXI" } else { puts "Data comparison test fail for -S00_AXI, expected-$wdata_2 actual-$rdata_tmp" inc ec } # Check error flag if { $ec == 0 } { puts "PTGEN_TEST: PASSED!" } else { puts "PTGEN_TEST: FAILED!" } Listing 6: Tcl Code for burst test V. CONCLUSION Effective hardware design requires specific skills and often focuses on low-level implementation considerations. This can be mitigated by providing a broader view of the whole process.


Nathalie Herve

Nathalie Hervé-Fournereau

WATER REPRESENTATIONS AND WATER SERVICES IN RURAL AREAS : PARADIGM SHIFTS AND SOCIO-ECOLOGICAL TRANSITIONS

Objectives

Ø Study the evolution of the representations of water of local actors in rural territories Ø Focus on the singularity of 'autonomous' water services for the supply of drinking water Ø Assess the way in which the objectives of the socio-ecological transition guide the action of local water 'user-managers' in rural areas and modify their representations of water and aquatics ecosystems. Ø Study the representations of water & aquatics ecosystems in different local planning documents and the "visibility" of 'autonomous' water services in such territorial planification Ø Analyse the evolutions or transformations of water & aquatics ecosystems governance and the contribution of 'autonomous' water services in the territory of agglomeration of Redon.

The near future

Why not deepen these questioning in a comparative framework with Irish researchers ? Isn't it a Great Opportunity !

  Spaces and Societies » Ø 11 researchers CNRS -Staff of 6 engineers CNRS and 15 engineers and techniciens of other institutions -Total workforce : 280 (120 PhD students)Ø 3 axes :

representations and water services in rural areas : paradigm shifts and socio-ecological transitions Research Project Spirit

  Local Stakeholders: local actors in particular in charge of 'autonomous' water services Ø Arts & Research: ideographer to illustrate the research & exposition citizen association Redon Ø Exploratory research (one year 2023) : recent legal study on rights of Nature and EU Law (European Parliament) -National Research Program One water, water as public good (first workshop on representations of water and aquatic ecosystems in Rennes 2023).

	Emmanuelle HELLIER, Professor, Space and Urban Planning
	Main fields : public utilities, intermunicipality approaches, water resources, uses
	and management; France, Italy, Romania, Sénégal
	Ex. Research projects : Mobilization for wetlands conservation (Ramsar
	procedure) with Laurence Le Du (2019-2022); Water representations & water

o Production of inequalities and ranking in societies o Practices, experiences and representations of social space o Relocation of polity and policies : a conflicting process? services with N. Hervé-Fournereau; Rainwater uses and policies for sufficiency (interdisciplinary projects) Water Ø Senior researchers & PhD students : training through research Ø Interdisciplinary : Law -Geography -Biogeochemistry Ø Case study : agglomeration in the south of Ille et Vilaine (Brittany): the Redon agglomeration Ø

-Water representations & water services with E. Hellier -Research methodologies & socio-ecological transitions with A. Pomade -EU


Johan Medrano
email: johan.medrano@ucl.ac.uk

Abderrahmane Kheddar

Sofiane Ramdani

Assessing Time Series Correlation Significance: A Parametric Approach with Application to Physiological Signals

Correlation coefficients play a pivotal role in quantifying linear relationships between random variables. Yet, their application to time series data is very challenging due to temporal dependencies. This paper introduces a novel approach to estimate the statistical significance of correlation coefficients in time series data, addressing the limitations of traditional methods based on the concept of effective degrees of freedom (or effective sample size, ESS). These effective degrees of freedom represent the independent sample size that would yield comparable test statistics under the assumption of no temporal correlation. We propose to assume a parametric Gaussian form for the autocorrelation function. We show that this assumption, motivated by a Laplace approximation, enables a simple estimator of the ESS that depends only on the temporal derivatives of the time series. Through numerical experiments, we show that the proposed approach yields accurate statistics while significantly reducing computational overhead. In addition, we evaluate the adequacy of our approach on real physiological signals, for assessing the connectivity measures in electrophysiology and detecting correlated arm movements in motion 1 capture data. Our methodology provides a simple tool for researchers working with time series data, enabling robust hypothesis testing in the presence of temporal dependencies.

Introduction

Correlation coefficients, such as Pearson's or Spearman's, are fundamental tools for assessing linear relationships between random variables. Although originally not designed for dependent samples, these coefficients are widely used with time series in diverse fields. Notable examples can be found in neuroimaging, where Pearson's correlation coefficient is used to construct functional brain networks from functional Magnetic Resonance Imaging (fMRI) data [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Van Dijk | Intrinsic functional connectivity as a tool for human connectomics: theory, properties, and optimization[END_REF], electrophysiology data [START_REF] Zhang | A lightweight solution to epileptic seizure prediction based on eeg synchronization measurement[END_REF][START_REF] Naira | Classification of people who suffer schizophrenia and healthy people by eeg signals using deep learning[END_REF][START_REF] Ji | Correlation analysis of eeg alpha rhythm is related to golf putting performance[END_REF][START_REF] Zhong | Spatiotemporal evaluation of epileptic intracranial eeg based on entropy and synchronization: A phase transition idea[END_REF] or to relate brain signals to other behavioral parameters, such as movement data [START_REF] Lu | Evaluation of classification performance in human lower limb jump phases of signal correlation information and lstm models[END_REF]. In this paper, we address the estimation of the statistical significance of correlation coefficients for time series.

Fisher's variance-stabilizing transformation is useful to obtain simple bounds and significance of a correlation coefficient r. Under the null hypothesis that two sets of n independent data points are uncorrelated, the null distribution of the Fisher-transformed variable z = arctanh(r) is approximately normally distributed with mean 0 and standard deviation 1/ √ n -3. In other words, we have

√ n -3 arctanh(r) ∼ N (0, 1) (1)
In real time-series data, the assumption of pairwise independence underlying traditional correlation coefficient tests often falters due to dependencies between consecutive observations.

When data exhibit temporal dependencies, Eq. ( 1) is overconfident and fails to account for the loss of degrees of freedom due to temporal correlations. A pivotal solution, pioneered by [START_REF] Bartlett | Some aspects of the time-correlation problem in regard to tests of significance[END_REF][START_REF] Bartlett | Some aspects of the time-correlation problem in regard to tests of significance[END_REF], adjusts for this change of degrees of freedom by introducing a number of "effective degrees of freedom", or "effective sample size" (ESS), representing the independent sample size producing comparable test statistics. Under the null hypothesis that two sets of n dependent data points, with ESS ν (ν ≤ n), are uncorrelated, we can simply rewrite Eq. ( 1) to obtain the approximate distribution as:

√ ν -3 arctanh(r) ∼ N (0, 1) (2)
This enables hypothesis testing even when actual samples show temporal dependencies.

The approach pioneered by Bartlett, that had various extensions, e.g., [START_REF] Bartlett | On the theoretical specification and sampling properties of autocorrelated time-series[END_REF][START_REF] Bayley | The" effective" number of independent observations in an autocorrelated time series[END_REF][START_REF] Quenouille | Notes on the calculation of autocorrelations of linear autoregressive schemes[END_REF][START_REF] Pyper | Comparison of methods to account for autocorrelation in correlation analyses of fish data[END_REF][START_REF] Afyouni | Effective degrees of freedom of the pearson's correlation coefficient under autocorrelation[END_REF], relies on estimating the ESS from the sum of the product of the autocorrelation functions (ACFs) of the time series. Let ρ k and γ k be the autocorrelation at lag k of two time series, resp. x and y,

of n samples each. The ESS of the correlation coefficient between x and y is [START_REF] Quenouille | Notes on the calculation of autocorrelations of linear autoregressive schemes[END_REF]:

ν = n ρ 0 γ 0 + 2 n-1 k=1 (n -k) n ρ k γ k -1 (3)
While this estimator of the ESS has been widely adopted, it requires computing the sample ACF which can be computationally demanding and can also yield inaccurate estimates of the ESS due to accumulating noise (in the sum of Eq. ( 3)).

In this work, we assume a parametric form for the ACF to simplify the computation and estimation of the ESS. In particular, we show that the sum in Eq. (3) converges to an integral which can be analytically solved under a Gaussian (Laplace) approximation. The resulting expression relates the ESS to the average second spectral moment, also called roughness, of the pair of time series to correlate. Importantly, the roughness of a series can be estimated from the variance of its temporal derivatives. This scaffolds the central result of this work: for two processes x and y of length n with temporal derivatives ẋ and ẏ, the number of effective degrees of freedoms ν is approximately:

ν = n var( ẋ) + var( ẏ) 2π (4)
We show that this formula yields accurate estimation of the ESS and effectively provides a computationally effective method to evaluate the significance of correlation coefficients.

This article is organized as follows. First, we present the key contributions of our work, which encompass the derivation of the asymptotic expression for ESS, a straightforward approximation for Gaussian autocorrelation, and the introduction of an estimator based on the variance of the temporal derivatives of the processes. Then, we validate our approach through numerical experiments, demonstrating its effectiveness and its computational efficacy. Finally, we apply our methodology to real data to (i) highlight the significance of power-based connectivity measures in electrophysiology, and (ii) distinguish correlated and uncorrelated arm captured motion data.

A parametric estimator of the effective sample size

Our work builds on estimating the ESS using a Gaussian approximation to the ACF. In this section we show that the ESS converges to an integral; then that a Laplace approximation of this integral relates the ESS to the signals roughness; and finally, that common roughness estimators can be used to build an estimator of the ESS.

This work focuses on estimating the ESS for the correlation of smooth, wide-sense stationary (WSS) stochastic processes. A smooth WSS stochastic process x is defined by convolving white noise w with a smooth (C ∞ ), square-integrable function K : R → R + , in other word:

x(t) = K(t -s)w(s)ds (5)
For clarity and brevity, we reuse the same notation throughout this paper. By default x and y are two smooth WSS stochastic processes with zero mean, variances resp. σ x and σ y , and autocorrelations functions ρ and γ. For simplicity, we restrict our result section to the case where both x and y have the same autocorrelation, i.e., ρ = γ. This restriction is moreover motivated in Section 2.2.

Asymptotic expression for the ESS

This section introduces an asymptotic form for the ESS. We consider Eq. ( 3) for the ESS in the case of infinitely large n and infinitely small sampling interval of the ACF. Under these conditions, the sum in Eq. ( 3) converges to

lim n→+∞ n ν = +∞ -∞ ρ(τ )γ(τ )dτ (6)
This yields a new asymptotic ESS expression, ν ∞ :

ν ∞ = n +∞ -∞ ρ(τ )γ(τ )dτ -1 (7)
The resulting integral formulation is pivotal in deriving analytical expressions for the ESS of stochastic processes with known ACF. As for processes with known autocorrelation function, we can directly evaluate Eq. ( 7). In the next section, we develop a generic approximation based on a Laplace approximation of the integral.

Closed-form expression under a Laplace approximation

A general approach to approximate the integral in the asymptotic ESS expression is to use Laplace's method. In this section, we derive a general approximation to asymptotic ESS (Eq. ( 7)), based on a Laplace approximation of the integral (Eq. ( 6)). In brief, Laplace's method is used to approximate the integral of a arbitrary function with a unique global maximum using a simpler, Gaussian function around its peak [START_REF] Penny | Statistical parametric mapping: the analysis of functional brain images[END_REF]. Typically in physiological signals, ACFs have a mode at lag 0 and decrease to 0 for large lags, although they might exhibit long-range correlations. Thus, the area under the curve of the product of two such ACFs is expected to be contained around 0, with an attenuation of long-range correlations. This justifies using Laplace's method to approximate the asymptotic ESS expression.

To derive the expression under Laplace approximation, note that the second-order expansion of the ACFs ρ and γ around their mode at τ = 0 is similar to processes with Gaussian ACFs

ρ(τ ) = exp(- 1 2 |ρ ′′ (0)|τ 2 ) γ(τ ) = exp(- 1 2 |γ ′′ (0)|τ 2 ) (8)
where ρ ′′ (0) and γ ′′ (0) are the second spectral moment of the processes, i.e., the second-order derivatives of the ACFs at 0. The second spectral moment is a universal measure of "roughness" in the literature of stochastic processes [START_REF] Friston | Dem: a variational treatment of dynamic systems[END_REF][START_REF] Cox | The theory of stochastic processes[END_REF]. From this Gaussian form, we can explicitly evaluate the integral in Eq. ( 7):

∞ -∞ ρ(τ )γ(τ )dτ = 2π |ρ ′′ (0)| + |γ ′′ (0)| (9)
By substituting Eq. ( 9) in Eq. ( 7), the ESS under a Laplace approximation is

ν ∞ = n |ρ ′′ (0)| + |γ ′′ (0)| 2π (10)
To summarise, we can use a Laplace (Gaussian) approximation to approximate the ESS. This is equivalent to fitting the ACF using a one-parameter family of Gaussian functions parameterised by their second spectral moment, or roughness. This allows to derive a simple expression for the ESS that only involves the second spectral moment of each process.

Remark 1. When ρ ′′ (0) = γ ′′ (0), Eq. (10) reduces to n |ρ ′′ (0)|/π. When ρ ′′ (0) ̸ = γ ′′ (0), the ESS of the pair of signals given by Eq. ( 10) is similar to that of a fictive pair of signals with equal roughness (|ρ ′′ (0)| + |γ ′′ (0)|)/2. Thus, there is no need to consider furthermore the case with heterogeneous roughness, as it can be mathematically reduced to the case with a (homogeneous) roughness given by the arithmetic mean of the original roughnesses.

Estimating the second spectral moment

A well-established result from stochastic process theory is that the roughness can be conveniently estimated from the variance of the first-order temporal derivatives of the process [START_REF] Cox | The theory of stochastic processes[END_REF][START_REF] Adler | The geometry of random fields[END_REF][START_REF] Worsley | An unbiased estimator for the roughness of a multivariate gaussian random field[END_REF]. In this section, we leverage this result to construct an estimator of the asymptotic ESS.

Let x and y be two processes with autocorrelation functions ρ and γ and temporal derivatives

ν ∞ = n var( ẋ) + var( ẏ) 2π (11)
By construction, this new parametric estimator for the ESS depends only on the variance of the temporal derivatives of the time series. As such, its consistency and unbiasedness derive directly from those of the variance estimators (from direct application of the continuous mapping theorem).

Remark 2. The second spectral moment is also related to the expected number of zero-crossings of a process by Rice's formula [START_REF] Rice | Mathematical analysis of random noise[END_REF]. For didactic purposes, we show an alternative ESS estimator based on Rice's formula in Appendix A.

Numerical validation

We conduct numerical experiments on Gaussian processes with Gaussian autocorrelation functions (GPGA) to assess the validity and limitations of our approach. The numerical results are organised in four parts. First, we evaluate the quality of fitting the squared autocorrelations from the process roughness. Second, we compare the quality of the estimates of the ESS, corresponding to the integral of the squared ACFs, with existing methods. Then, we compare the statistics yielded by our method with existing ones. Finally, we investigate the computational performance of our approach.

Assessment of parametric autocorrelation function estimates

The method we present leverages the analytical evaluation of the ESS of a pair of stochastic processes from their roughness, estimated from the variance of the temporal derivatives. Here, we verify whether this method can accurately estimate the ACF of the process, when the ACF is known to be Gaussian.

We generate 2000 sample paths from GPGA with different roughness levels. We then retrieve the squared autocorrelation using the variance of the derivatives. To sample a sample path from a GPGA with roughness r, we convolve a unit white Gaussian noise sample with a Gaussian kernel of variance 2/r as per [START_REF] Friston | Dem: a variational treatment of dynamic systems[END_REF]. We then apply Eq. ( 4) to determine its second spectral moment from the variance of the process temporal derivatives.

Figure 1 illustrates that low roughness results in significant biases in squared sample autocorrelation functions, leading to ESS underestimation. This bias becomes clear in Figure 1

(f),
where the grey curves diverge significantly from their expected near-zero values. This figure underscores the merit of predetermining the functional form of the ACF, thereby filtering out spurious long-term correlations. Furthermore, our method consistently aligns more closely with the theoretical ACF relatively to the direct sample squared ACF.

In conclusion, explicitly formulating the functional form of the ACF filters out errant longterm correlations otherwise amplified by traditional ESS estimators based on sample autocorrelations. This leads to more accurate and consistent ACF estimates across varying roughness scales. However, we note that the variance of the squared ACF appears to be inversely proportional to roughness. This issue is further investigated in subsequent sections. Figure 1: Sample paths and squared autocorrelation functions for roughness values ρ ′′ (0). Rows one and two display, respectively, samples of Gaussian processes for different roughness values and their squared autocorrelation functions. Notably, the squared ACF plays a crucial role in ESS estimation. The plots showcase various ACF estimations: theoretical (plain yellow), sample (plain grey), average of sample (dotted black), Gaussian fit based on estimated roughness (plain red), and Gaussian fit based on average roughness (dashed blue). The bias in sample ACFs, evident in the grey curves, increases as roughness decreases due to increasing random long-range correlations. On the other hand, the Gaussian fits (plain red) remain unbiased at large lags, ensuring the ACF retains its Gaussian form.

Assessment of parametric ESS estimates 1

This section deals with the influence of the roughness estimator's variance on the variance of 2 the estimated sum of squared autocorrelation and that of the ESS. We also examine how these 3 relate to the length of series.

4

We sampled 1000 sample paths from GPGA, adjusting roughness between 10 -6 and 1. Each 5 sample path's roughness is then estimated. Using this roughness estimate, we compute the 6 ESS using Eq. ( 11). Furthermore, we determine the ACF of the process using the inverse 7 Fourier transform of the sample path's power spectral density (PSD). The latter is estimated 8 either through the FFT or the Welch periodogram with a 256-point window having 128-point 9 overlap [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF]. The ESS is then obtained from the ACF using Eq. (3). We replicate this 10 procedure for series lengths of 500, 1000, and 2000 points. For comparison, the ESS is normalized From Figure 2(a), it is clear that roughness estimates share a similar trend. For minuscule 1 roughness values, bias emerges. Extended series better approximate low roughness, hinting that 2 biases result from series that are too short relative to their roughness. As expected from their 3 relationship (Eq. ( 7)), ESS factors derived from roughness mirror this trend. 11

The trends in Figure 3 are more comprehensible when compared with the ESS estimates presented in Figure 2. The transition from a roughness of 10 -2 to 10 -4 equates to a shift from an ESS factor of roughly 10 -1 to 10 -2 . Such changes mirror the variances and biases seen in our ESS estimates, which directly affect p-value computation. Given that our method consistently provides superior ESS estimations across roughness values, it is hence inferred to offer enhanced statistics.

In summary, our numerical analyses confirm the paramountcy of a parametric autocorrelation function, where the unique parameter is estimated via the process's second spectral moment.

This method circumvents biases introduced by random long-range dependencies, resulting in more accurate sample p-values for correlation coefficients. In subsequent sections, we advocate for GPGA's aptness in real-world applications, such as neuroimaging data connectivity analysis and movement trajectory correlations.

Comparison of computational performance

In this section, we evaluate the performance gain of our approach. The main difference with existing methods is that we avoid explicit evaluation of the ACF by using a parametric ACF, which only requires to estimate the average series roughness.

To evaluate the impact of series length on the computation time and speedup, we generate pairs of sample paths from a GPGA process with roughness 10 -3 having varying length, from 100 to 10 7 points with 10 fold increments. For each pairs of paths, we compute the ESS using FFT-based and Welch-based computation of the ACFs, and with our approach. The FFT-based ACF is obtained using numpy's fft and ifft functions [START_REF] Harris | Array programming with NumPy[END_REF]. The Welch-based ACF is computed using scipy's welch function from the signal package with a window length of 256 points (100 points for the 100 points series) and numpy's ifft [START_REF] Harris | Array programming with NumPy[END_REF].

Our approach uses numpy's diff and var to compute the variance of the temporal derivatives of the process. The results are obtained on a laptop equipped with an 12-core Intel Xeon W-10855M CPU and 32Gb of RAM running Ubuntu 22.04 and Python 3.8. For each approach, we computed the timings from 100 loops and present the mean and standard deviation across 7 different runs.

We report the results on Speedup of our approach, obtain as the average computation time of each method divided by the average computation time of the proposed method. For instance, a speedup of 5 means that our approach divides by 5 the average computation time.

tational advantage of the proposed approach. In average across different number of points, the 1 proposed approach is 7.4 times faster than an FFT-based approach. This speedup is negligible 2 for short time series (resp. 1.3 and 1.6 for 100 and 1000 points) but significant for long time se-3 ries (10.2 to 13.4 for 10 5 to 10 7 ). The speedup from Welch-based approach is relatively constant 4 across time series length due to the fixed window length, with an average speedup of 5.1.

5

To sum-up, our proposed approach, based on the temporal derivatives of the time series,

6
gives a significant speedup as compared to other approaches based on the computation of the 7 ACF. This computational argument motivates using a parametric Gaussian form of the ACF in 8 applications where timing is critical or where large amounts of data have to be processed; for 9 instance, when evaluating correlations of brain activity between a large number of brain regions. Step 4.

Test against the null hypothesis

Step 1.

Compute the time-frequency decomposition

Step 2.

Extract the 4th roots of power time series p-value

Step 3. uate the functional connectivity between brain regions. A standard approach, power-based connectivity analysis, measures functional connectivity by evaluating the correlation between the time-frequency power of the signal at two different locations [START_REF] Cohen | Analyzing neural time series data: theory and practice[END_REF]. Several methods exist to perform power-based connectivity analysis, but the correlation between wavelet-based time-frequency representations is particularly interesting within the scope of this work.

Compute the correlation coefficient and the ESS

The time course of power can be obtained using a continuous wavelet transform with a Morlet wavelet. A Morlet wavelet combines a Gaussian kernel with a complex sinusoid at a specified frequency. Therefore, a continuous wavelet transform with a Morlet wavelet introduces a form of temporal smoothing with a Gaussian kernel. Intuitively, if the width of the Gaussian smoothing kernel is large enough relative to the autocorrelation of the signal, the resulting time series will have roughly a Gaussian autocorrelation. In this case, one can consider testing the strength of the correlation between two power time series under the null hypothesis that the signals are two uncorrelated GPGA. The overall procedure is illustrated in Figure 5. This approach allows a more efficient analysis of connectivity significance without relying on computationally demanding methods.

This section aims at gaining some insight on whether the null hypothesis of uncorrelated GPGA processes might yield appropriate statistics. In contrast to the numerical experiments section, we cannot sample from the null distribution of correlation coefficient between bandpower time series. Thus, we can only investigate a few characteristics of the data to assess the validity of our approach. Indeed, this section echoes the first step of statistical analysis in real world application, where statisticians have no access to the data generative distribution and shall assess the suitability of their model.

Data presentation and preprocessing

We analyse a dataset containing EEG recordings of a subject listening to continuous, naturalistic speech [START_REF] Di Liberto | Low-frequency cortical entrainment to speech reflects phoneme-level processing[END_REF][START_REF] Crosse | The multivariate temporal response function (mtrf) toolbox: a matlab toolbox for relating neural signals to continuous stimuli[END_REF]. The subject listened to a classic work read by an English speaker. EEG signals were recorded from a single participant with a Biosemi system having 128 channels and a sampling rate of 512 Hz. Prior to analysis, the data has been filtered between 1 and 15 Hz and downsampled to 128Hz. We then computed the time-frequency representation of the EEG signals using a continuous wavelet transform. We used Morlet wavelet with 7 cycles to produce time-frequency representation with a sufficient temporal smoothness to fall under our Gaussian assumptions. Finally, we take the fourth root of the power to make the data marginally Gaussian [START_REF] Hawkins | A note on the transformation of chi-squared variables to normality[END_REF].

Remark 3. Note that selecting the number of cycles in the wavelet effectively selects the degree of smoothness of the signal. This shows that the number of cycles is related to the ESS. In Appendix B, we give a formula giving the ESS as a function of the number of cycles.

Adequacy of GPGA assumptions for wavelet-based EEG power

Here, we show that our approach yieds appropriate statistics for power-based connectivity analysis. We first look at whether the signals have an approximately Gaussian autocorrelation and whether our approach provides a good fit of the squared autocorrelation function. We compare results given by FFT-based and Welch-based approaches with our approach, and analyse the results in the light of our numerical results. Finally, we compare the ESS and 97.5% quantile given by all three approaches with real data and sample GPGA sample paths with matching roughness.

Figure 6 summarizes the properties of the data. In particular, we show a time-frequency plot (Fig. 6 (a) and (b)) and some samples power time series at different frequencies. Intuitively, the time course of power at lower frequencies is smoother than for higher frequencies. This is confirmed by looking at the roughness of the signals across channels, whose mean increase with frequency (Fig. 6 (d)). In addition, we display (Fig. 6 (c)) the marginal distribution of the signals (e-g) Squared ACF at 4 Hz, 10 Hz, and 14 Hz. Gray lines indicate different channels, the black line indicates the mean squared ACF across channels, the dotted blue line indicates the mean fit using the proposed approach, and the filled blue area shows the area that is integrated in the denominator of the ESS expression.

across channels to inspect its normality -a key element in mandating using Pearson's correlation 1 coefficient instead of Spearman's rank correlation. Finally, we inspect (Fig. 6 (e), (f), and(g))

2
the squared ACF at 4 Hz, 10 Hz, and 14 Hz, and compare it against the Gaussian fit obtained 3 from the signals roughness. We observe that the squared ACF is roughly Gaussian, despite a 4 slower decrease (lower kurtosis) which might come from autocorrelations inherent to the nature 5 of the signal. Overall, these results suggest that GPGA might provide a good approximation 6 for the data and that our approach gives appropriate statistics.

Assessment of ESS and quantiles for wavelet-based EEG power

To confirm these results, we evaluate the ESS factor for all possible between-channel pairs connectivity at each given frequency. In addition, for each signal we generate a random GPGA sample path with matching length and roughness, and compute the ESS factor for all pairs of random sample paths. This gives us a hint on how the ESS factor would behave if the sample paths where effectively sampled from the null hypothesis, i.e., GPGA with matching length and roughness. We complement this analysis by looking at the 97.5% quantile of both real and simulated data, under the null hypothesis. Any correlation coefficient greater than the 97.5% quantile -or lower than its opposite -would be effectively considered as significant with p < 0.05. The 97.5% quantile is computed by the inverse Fisher transform of the 97.5% quantile of the Gaussian distribution:

Q 97.5% = tanh(1.96/ √ ν ∞ ) (12)
The ESS factor and 97.5% quantile for both real and simulated data are presented on Figure 7.

We observe that the ESS factor for real EEG signal increases with frequency. This is expected as the roughness of power time series increases with the frequency, as shown in Figure 6. This increased ESS factor results in a decreased 97.5% quantile. Indeed, as the effective sample size increases and statistical power increases, we can arbitrate on the rejection of the null hypothesis at lower values of correlation coefficient. This is conform to the role played by the ESS in Eq. ( 12).

Interestingly, we see that both FFT-based and Welch-based approaches give lower values of ESS, and higher 97.5% quantile value, than our proposed approach. This effect is the strongest at 10 Hz and 12 Hz. Comparing each figure with its analogous generated under GPGA assumptions, we see that under GPGA assumptions all three methods yield approximately the same values of ESS factor and 97.5% quantile. Because these differences are stronger at 10 Hz and 12 Hz, i.e., in the well-identified alpha band, we hypothesis that task-related modulations of the power task cause variations in the autocorrelation of the power; which induce variability and shift in the ESS factor and related quantiles. Note that the size of this effect on the 97.5% quantile is relatively small (less than 0.02). We argue that this is negligible in most cases, especially given the variability of the sample estimates of the roughness or ACFs. Gaussian. Thus, we aim at evaluating our proposed approach on less 'well-behaved' biological 7 signals.

8

The motion capture data analysed in this section comes from a single subject performing self-paced unilateral elbow flexion/extension. The data was acquired during a study on the EEG correlates of arm movements. The experimental procedure conformed the Declaration of Helsinki and was approved by the local ethics committee. The subject was instructed to move continuously in a self-paced pseudo-random manner for 40 runs (one per side) of 23.5 seconds.

The position of body segments was recorded using an XSens Awinda suit at a 60 Hz rate and the elbow flexion angle was obtained as the angle between the upper arm and forearm.

We computed the marginal distribution, squared ACF, and roughness of each trial. In addition, we computed the ESS factors and the 97.5% quantile for the correlation coefficient between every possible pair of trials. Similarly to the previous section, we generated for each trial a random GPGA sample path with matching roughness, and computed the same quantities as for the real data. Results are presented in Figure 8.

Looking first at the statistics of the joint trajectories, we see that the joint angle distribution is not Gaussian, nor even unimodal. In practice, this would motivate to use Spearman's correlation coefficient instead of Pearson's. Additionally, we observe that the squared ACFs display correlations for lags around 200 points, that are not matching the Gaussian profile. On the other hand, the main peak of the squared ACF seem well approximated by the Gaussian fit.

We also observe that the signal roughness is quite low, with an average roughness of 3.7 × 10 -3 .

Interestingly, the ESS factors and 97.5% quantiles computed for both the real joint trajectories and random GPGA processes are relatively similar, with differences being much smaller than the spread of the values. This result highlight that overall, the ESS of pairs of joint trajectories from this dataset are quite similar to that of GPGA with matching roughness. Consequently, quantiles are also relatively similar. This comparative results shows that, despite being not marginally Gaussian and having long-range autocorrelations not fitting to the Gaussian profile, the significance of the correlation between joint trajectories from this dataset could be accurately computed under GPGA assumptions.

Conclusion and future directions

In this work, we considered the problem of testing for significant correlation between autocorrelated time series. Several works have proposed to use an ESS to correct the test distribution.

The ESS correspond to the size of a sample of independent observations that would produce the approach retrieves the statistics of correlation coefficients under the null hypothesis that signals are GPGA. For this particular type of process, our approach outperforms classical approaches, yielding robust estimate of the statistics under a large range of roughness. In addition, our approach shows relatively higher computational performances compared to alternative methods.

In a second part, we show that some data issued from biological signals satisfy GPGA assumptions. Our approach seem adequate to test for significant correlations between brain regions, when used with power-based connectivity analysis of electrophysiological data. Our method also applies well to random joint trajectories, as measured from a motion capture system.

Overall, obtained results suggest that our approach yields accurate and robust statistics under assumptions that can be found in real biological signal. More generally, we claim that using a parametric form of autocorrelation can yield more accurate statistics, without being overly restrictive on the form of the process. In particular, we think that more advanced methods for parameter estimation could allow better estimation of the autocorrelation function. However, as often with parametric approaches, there is no one-size-fits-all form of the autocorrelation function, and the statistician holds the responsibility to evaluate the adequacy of the model for the data under study.

In addition to correlation coefficients, our approach can be straightforwardly extended to linear regression. A pivotal assumption of linear regression is the independence of residuals.

An assumption often violated in time series data due to inherent temporal dependencies. The parametric ESS estimators derived here can be seamlessly adapted to assess the significance of regression coefficients in linear regression models applied to time series data. Thus, our approach not only fortifies the foundation of correlation analysis in time series but also extends its robustness to the broader landscape of linear regression.

A Link with Rice's formula

The second spectral moment is related to the average number of zero-crossings N 0 of the process by Rice's formula [START_REF] Rice | Mathematical analysis of random noise[END_REF] (for derivation, see (Cox and Miller, 2017, Ch. 7.4, Example 7.9)):

N 0 = n |ρ ′′ (0)|/π (13)
This expression can be used to express the ESS in terms of expected number of zero-crossings.

We observe that

ν ∞ = √ πN 0 (14)
This gives a simple way to estimate the ESS of empirical time series from the mean number of zero-crossings.

B Analytical ESS for Wavelet-based correlations

This section concerns testing correlations in band-limited power using time-frequency representations, particularly with Morlet wavelets, which align well with the ESS theory for quasi-Gaussian autocorrelation processes.

Morlet wavelets, defined as the product of a complex sinusoid with a Gaussian, provide a time-frequency representation of a time series through the wavelet transform. Specifically, a wavelet w is given by (Cohen, 2019):

w = exp(-j2πf t -t 2 /(2σ 2 )) ( 15
)
where t is time, f is frequency, and σ is tied to the number of cycles N c as σ = N c /2πf .

For signals with a rapidly decaying autocorrelation compared to the Gaussian kernel width in the wavelet, the time-frequency representation remains relatively smooth. Considering a time series derived from the log-power (or other normalisation function) of a wavelet transform at a given frequency, the autocorrelation ρ is approximately Gaussian with second-spectral moment

|ρ ′′ (0)| ≈ π 2 f 2 N 2 c ( 16
)
of cycles in Morlet wavelets can be adjusted, providing a direct assessment of the ESS for the correlations. If two series, derived from the same number of cycles but different frequencies f 1 and f 2 , are correlated, their ESS is:

ν ∞ = n N c π 2 f 2 1 + f 2 2 (17)
For f 1 = f 2 = f , the ESS is simply ν ∞ = √ πnf /N c . Importantly, the chosen number of cycles determines the ESS of power correlations.

ẋ

  and ẏ. The second spectral moments |ρ ′′ (0)| and |γ ′′ (0)| can be estimated from the variance of the temporal derivatives of the process, i.e., |ρ ′′ (0)| = var( ẋ) and |γ ′′ (0)| = var( ẏ). Pluging in Eq. (10), we obtain an estimator of the ESS:

  Squared ACF, ρ ′′ (0) = 10 -4

Figure 2

 2 Figure 2: a) Estimated roughness as a function of the process roughness, for series of 500 points (dark blue), 1000 points (light blue), and 2000 points (yellow). b-d) Estimated ESS factor as a function of the process true ESS, computed from Eq. (10), for series of 500 points (dark blue), 1000 points (light blue), and 2000 points (yellow) for the FFT-based (b),Welch-based (c), and proposed approaches (d). Shaded areas correspond to the 95% confidence intervals.

4Figure 3 :6

 3 Figure3: Estimated probability against empirical probability for the different approaches. The first row illustrates a roughness of 10 -2 , while the second showcases 10 -4 . Each column represents different estimation methods: FFT-based (green),Welch-based (red), and our method (blue).

Figure 4 .Figure 4

 44 Figure4: a) Comparison of the computation time for varying number of points from 100 to . b) Speedup of our approach, obtain as the average computation time of each method divided by the average computation time of the proposed method. For instance, a speedup of 5 means that our approach divides by 5 the average computation time.

  ) captures the electric potentials on the scalp. Its low-cost and 15 high-temporal-resolution make EEG a particularly appealing neuroimaging technique to eval-

FOIFOIFigure 5 :

 5 Figure 5: Procedure to assess the connectivity of two EEG channels at a particular frequency of interest (FOI).

Figure 6 :

 6 Figure 6: (a) Time-frequency plot of the fourth root of power for channel A9. (b) Sample fourthroot of power time series at different frequencies. (c) Marginal distribution of fourth-root of power time series across channels. (d) Roughness across channels at different frequency. (e-g) Squared ACF at 4 Hz, 10 Hz, and 14 Hz. Gray lines indicate different channels, the black line indicates the mean squared ACF across channels, the dotted blue line indicates the mean fit using the proposed approach, and the filled blue area shows the area that is integrated in the denominator of the ESS expression.

  7

Figure 7 :2

 7 Figure 7: (a) Distribution of the ESS factor evaluated with the FFT-based (green), the Welch-based (red), and our approach (blue) for pair of power time series at a given frequency. (b) Similar to (a) for random GPGA sample paths with roughness matching that of every pair of signals at a given frequency. (c) Distribution of the 97.5% quantile obtained from the ESS of each pair of real time series. (d) Similar to (c) for random GPGA sample paths with matching roughness.

same statistics. These approaches require to estimate the autocorrelation function, which can 1 be problematic when testing correlation between two unique and smooth series.

2 Here, we derive the asymptotic expression of the ESS. This asymptotic expression can be used 3 with parametric forms of ACF to analytically derive the ESS of a given process. In particular, 4 we show that the expression for the ESS takes a simple form when considering a Gaussian 5 approximation of the ACF. This simple expression depends on the second-order derivatives of 6 the autocorrelation function at its mode, which is a measure of the roughness of the process 7 and relates to several of its statistical properties. Roughness can be easily estimated from the 8 variance of the first-order temporal derivatives of the process.

9

We conducted numerical experiments to validate the proposed method. We observe that our 10

Supporting Information

The Wellcome Centre for Human Neuroimaging is supported by core funding from Wellcome [203147/Z/16/Z]. This research was partly supported by a grant from the French Ministère de l'Éducation Nationale, de l'Enseignement Supérieur et de la Recherche.

Declaration of Generative AI and AI-assisted technologies in the writing process

During the preparation of this work the author(s) used OpenAI ChatGPT 4 in order to improve the readability and phrasing. After using this tool/service, the author(s) reviewed and edited the content as needed and take(s) full responsibility for the content of the publication.


Arnaud Quérel

Denis Quélo

Yelva Roustan

Anne Mathieu

Cloud diagnosis impact on deposition modelling applied to the Fukushima accident

At the Fukushima airport, a deposit as large as 36 kBq.m -2 of Cs-137 was measured. Both dry and wet deposition were probably involved since a raining event occurred on the 15th of March when the plume was passing nearby.

Radar RAP vs cloud diagnosis

Comparison between observation and diagnosis at one station

  In case of a nuclear accident, they are two phases concerning the dispersion of radioactive materials: 1. Forecast: Anticipating the consequences of an atmospheric release of radioactive material. 2. Aftermath: Understand the soil contamination and the possible harm suffered by the populations during the event. Wet deposition modeling is important to achieve these goals. The cloud diagnosis is a key issue for wet deposition modelling since it allows distinguishing between two processes: in-cloud scavenging: the collection of radioactive particles into the cloud below-cloud scavenging: the removal of radioactive material due to the falling drops. Which cloud diagnosis to use for the atmospheric transport models ? Cloud diagnosis choice have a major impact to the volume of the atmosphere considered as -. Then, the repartition between in-cloud and below-cloud may be strongly impacted by the cloud diagnosis. Cloud water mixing ratio (Q C ) is the most interesting variable, which describes only the cloud water. Q c provides satisfactory results and is not sensitive to the threshold. We therefore recommend to use Q c to distinguish in-cloud and below-cloud scavenging in the atmospheric transport modelling. Comparison of observed and diagnosed cloud maps, the 16 th March 00:00 EGU 2017 8845

Q

  Models cannot currently reproduce the Fukushima accident in a completely satisfactory manner. One of the potential improvement may concern the cloud diagnosis on this case. The period studied is around the 15 of March 2011; This is the period where the largest deposits occurred and several cloud observations were available. The meteorological model used here is a 3km/10min resolution simulation provided by JMA. To model wet deposition, the variables of interest are: *Cloud diagnoses: Relative humidity, Liquid Water Content (LWC) and cloud water mixing ratio (Q C ) are used for diagnoses. Two thresholds are used for each diagnosis. **Precipitations are only used here to verify the consistency with cloud presence. this case, the meteorological model does fairly reproduce the cloud cover for low and middle layer. Upper layer clouds are not well reproduced. But, these types of clouds are not important for deposition modelling. At this station, cloud thickness is very sensitive to the cloud diagnosis, contrary to the cloud presence. Relative humidity at 80% and LWC at 10 4 Kg.m 3 give poor performances. LWC at 10 4 Kg.m 3 is not usable. Other diagnosis are equally coherent with RAP data, satellite diagnosis include. The RH 80% is the best one, but only due to its cloud overestimation. ., Shimbori, T., Draxler, R., 2015. regional atmospheric transport model calculations for the WMO technical task team on meteorological analyses for Fukushima Daiichi Nuclear Power Plant accident. Journal of Environmental Radioactivity 139, 185 199. doi:10.1016/j.jenvrad.2014.02.007 Thanks to Dr Sekiyama and Dr Kajino for the Himawari, AMeDAS and meteorological modelling data providing. Thanks to Japan Meteorological Agency Sendai Airport Observatory to the Grant-in-Aid for Scientific Research on Innovative Areas ISET-R modeling team and Pr Tsuruta for Fukushima airport data. product (for the cloud top) (Himawari) Ceilometer (for cloud base) Radar, rain gauge, radar/rain gauge-analyzed precipitation (RAP, Saito et al., 2015) Meteorological model Simulation output Diagnosed* with help of simulation output Simulation output