Sara Kassan
email: sara.kassan@uha.fr

Jaafar Gaber
email: jaafar.gaber@utbm.fr

Pascal Lorenz
email: pascal.lorenz@uha.fr

Autonomous Energy Management system achieving piezoelectric energy harvesting in Wireless Sensors

Keywords: WSN, Network lifetime, Energy harvesting, Piezoelectric, Power consumption, Proportional Integral Derivative controller PID, energy consumption, Microprocessor Controller Unit MCU, Mica2 motes sensors

Wireless Sensor Networks (WSNs) are extensively used in monitoring applications such as humidity and temperature sensing in smart buildings, industrial automation, and predicting crop health. Sensor nodes are deployed in remote places to sense the data information from the environment and to transmit the sensing data to the Base Station (BS). When a sensor is drained of energy, it can no longer achieve its role without a substituted source of energy. However, limited energy in a sensor's battery prevents the long-term process in such applications. In addition, replacing the sensors' batteries and redeploying the sensors is very expensive in terms of time and budget. To overcome the energy limitation without changing the size of sensors, researchers have proposed the use of energy harvesting to reload the rechargeable battery by power. Therefore, efficient power management is required to increase the benefits of having additional environmental energy. This paper presents a new selfmanagement of energy based on Proportional Integral Derivative controller (PID) to tune the energy harvesting and Microprocessor Controller Unit (MCU) to control the sensor modes.

I. INTRODUCTION

The possibility to avert replacing drained batteries is extremely important in Wireless Sensor Networks to avoid the high cost and time to replace batteries and redeploy them. Therefore, energy harvesting has attracted researchers and developers to recharge low power devices as Vivo-Nano-Robots, MEMS, Claytronics and WSN [START_REF] Dolev | In-vivo energy harvesting nano robots[END_REF] [2] [START_REF] Sherazi | A comprehensive review on energy harvesting MAC protocols in WSNs: Challenges and tradeoffs[END_REF]. Energy harvesting is the method of extracting energy from the environment through different sources of energy. It is considered an emerging and reasonably mature technology to overcome the limited lifetime of battery-operated wearable devices and allows continuous recharging of the energy stored during use. The environmental energy for scavenging is mostly provided by ambient light (artificial lighting and solar lighting), radio frequency, thermal sources, kinetic and vibration sources [START_REF] Shaikh | Energy harvesting in wireless sensor networks: A comprehensive review[END_REF] [START_REF] Adu-Manu | Energy-Harvesting Wireless Sensor Networks (EH-WSNs): A Review[END_REF]. To harvest energy from vibration, there are different techniques such as electrostatic, electromagnetic, magnetostrictive and piezoelectric [START_REF] Narita | A Review on Piezoelectric, Magnetostrictive, and Magnetoelectric Materials and Device Technologies for Energy Harvesting Applications[END_REF] [7] [START_REF] Siang | Review of vibration-based energy harvesting technology: Mechanism and architectural approach[END_REF]. Table 1 compares the features of different type of vibration energy harvesting techniques. Among them in this paper, we focus on the control of the low-level energy harvesting by the most prevalent technique piezoelectric vibration system for WSN. The diagram of a piezoelectric energy harvesting system is presented in Figure Fig. 1; it can be shortened into three essential components: piezoelectric devices, converters, and electrical energy storage. In addition to these components, we introduce two controllers one to tune the energy harvesting and the other to manage the energy in the sensor according to its activities needs. This new self-organization power system will be able to insure the prolongation of the WSN lifetime.

TABLE Ⅰ: comparison of the different vibrational types of harvesting mechanisms

Figure 1: Block diagram of energy harvesting model

The paper is organized as follows: Section Ⅱ provides a review with the different mechanisms to save energy in WSN.

In section Ⅲ, we present the global proposed model to manage the harvested energy in WSN. In section Ⅳ, we describe the piezoelectric harvesting energy model. In particular, we show how to adapt the PID controller to tune power in the vibrationbased model in WSN and how to tune the residual energy in the rechargeable batteries according to the sensor's needs using MCU. In section Ⅴ, we present the simulation results. Finally, we provide a conclusion of the paper in section Ⅵ.

II. RELATED WORK

Prolonging WSN lifetime is one of the most critical challenges in WSNs. In the last few years, researchers have proposed different methods to increase the WSN lifetime via saving approaches. The main goal of these methods is saving energy, during the communication between sensors, without compensation for the energy dissipated via scavenging ambient energy. Specifically, we can classify energy methods in two types: energy saving mechanisms to decrease energy consumption and energy harvesting to recompense the energy consumed during sensors' activities.

Most of the proposed approaches in literature are saving energy-based mechanisms such as sensors deployment strategies, data-oriented techniques, topology control methods, and energy routing protocols. For example, in [START_REF] Kulkarni | A Novel Sensor Node Deployment using Low Discrepancy Sequences for WSN[END_REF], the authors propose a statistical strategy for the nodes' deployment. The statistical node deployment strategy uses the Quasi-random method of low-discrepancy sequences to increase the lifetime and the coverage of the network.

The attentiveness to data-oriented techniques has increased recently. These techniques reduce the size of information via compression. Accordingly, data-oriented techniques are efficient to save energy during transmission and reception of this data information. In addition, the reduction of data is not only efficient in saving energy but also in saving more memory in the sensors [START_REF] Wang | Energy efficient distributed compressed data gathering for sensor networks[END_REF] [START_REF] Leinonen | Distributed Distortion-Rate Optimized Compressed Sensing in Wireless Sensor Networks[END_REF]. For example, in [START_REF] Leinonen | Distributed Distortion-Rate Optimized Compressed Sensing in Wireless Sensor Networks[END_REF], the authors propose a Distributed Distortion-Rate Optimized Compressed Sensing (DQCS) method to compress data in WSN under a complexityconstrained encoding, which minimizes a weighted sum between the mean square error (MSE) signal reconstruction distortion and the average encoding rate.

A local lossless neighborhood indexing sequence (NIS) compression algorithm for data compression in WSN is proposed in [START_REF] Uthayakumar | A new lossless neighborhood indexing sequence (NIS) algorithm for data compression in wireless sensor networks[END_REF]. The NIS approach dynamically assigns shorter length code information for each character in the input sequence, by exploiting the occurrence of neighboring bits, and every data packet is decompressed independently from others. Therefore, NIS algorithm helps to decrease the network load, which results in low data packet loss.

A topology control protocol based on learning automaton is proposed in [START_REF] Javadi | Learning automaton based topology control protocol for extending wireless sensor networks lifetime[END_REF]. The mechanism chooses the proper smallest transmission range of the node using the reinforcement signal produced by the learning automaton of neighbor sensor nodes. Consequently, the choice of transmission range affects the energy consumed by nodes and the overall network lifetime.

A hierarchical clustering protocol based on sensors location and energy consumption using non-cooperative game theory (GT) approach to extend the WSN lifetime is proposed in [START_REF] Kassan | Game theory based distributed clustering approach to maximize wireless sensors network lifetime[END_REF]. The GT permits a sensor to decide between two actions: to enter a game and transmit a message or to stay out of the game and harvest energy to charge its battery to reach the Nash Equilibrium (NE) solution for mixed strategies.

In [START_REF] Wang | Genetic algorithm for energy-efficient clustering and routing in wireless sensor networks[END_REF], a Genetic algorithm-based Energy-efficient Clustering and Routing Approach (GECRA) is presented. This algorithm aims to calculate the total energy consumed by all sensor nodes where the algorithm encodes the clustering scheme and routing scheme together in the same chromosome.

It is worth noting that these research efforts aim mainly to reduce the energy consumption in the stage of wireless sensors deployment and during sensors activities, such as communication and sensing data, via energy saving mechanisms. Therefore, the WSN lifetime will be extended, but for a limited time. Instead of applying the energy saving methods, energy harvesting can be used to recharge the sensor's battery. The objective is to increase the power stored in the sensor's battery thereby extending the overall network lifetime depending on the availability of ambient energy harvesting resources. This article presents a new self-management of energy based on a Proportional Integral Derivative (PID) controller to tune the energy harvesting and a Microprocessor Controller Unit (MCU) to control the sensor modes.

III. PROPOSED SELF POWER MANAGER SYSTEM MODEL

The target is an autonomous system to harvest energy spontaneously upon reaching a threshold of energy consumption. The system proposed is composed of two main blocks as designed in figure Fig. 2: The Battery Management Subsystem block that manages a rechargeable battery or a supercapacitor of the sensor and monitors its states via two (Management Controller Unit) MCU modes. In on mode, the sensor is active and can sense information from its environment, write/read data from memory, and communicate messages with other wireless sensors. During achieving activities, the sensor consumes energy until a threshold and the MCU turns to off mode. In the off mode, the block piezoelectric energy harvesting structure is triggered. This block is controlled via a control loop feedback using a Proportional Integral Derivate controller PID that scavenges ambient energy to recharge the sensor's battery.

In this paper, we assume that the vibration resource exists continuously. Moreover, the piezoelectric harvesting energy system is modeled by an equivalent Mass-Spring-Damper (MSD) model (discussed in the section below Ⅳ) with a PID controller that adapts the quantity of energy harvested by the sensor node to still alive and continue its processes. The harvesting energy will be tuned to reduce the error 𝑒 between the setpoint or reference energy 𝑅𝐸.

When the error 𝑒 is bigger than a threshold 𝑇, the sensor node dissipates the majority of its residual energy and the microcontroller MCU switches to off mode and lets the sensors node harvest energy. Otherwise, the MCU is still in on mode and continues normally its activities. The energy consumption for one sensor depends on the energy consumed for its activities. The system block diagram is shown in Figure Fig. 2.

Figure 2: Architecture of the global system of wireless sensor composed of two main blocks one to manage energy consumed and residual energy in the battery and another one to harvest energy from the ambient source

In what follows, the model is formulated mathematically in terms of energy consumption by sensor nodes. The energy error 𝑒 can be written as follows:

𝑒(𝑡) = 𝑅𝐸 -𝐸 𝑟 (𝑡) ( 1
)
where the reference energy 𝑅𝐸 is a fixed energy harvesting to charge the sensor's battery and 𝐸 𝑟 is the current residual energy in the sensor's battery.

𝑒(𝑡) = 1 2 𝐶(𝑉 𝑟𝑒𝑓 2 -𝑉 𝑟 (𝑡) 2 ) ( 2
)
where the reference voltage 𝑉 𝑟𝑒𝑓 is the desired fixed voltage and 𝑉 𝑟 (𝑡) is the current voltage for the sensor's battery in time 𝑡.

𝐸 𝑟 (𝑡) = 𝐸 ℎ (𝑡) + 𝐸 𝑟 (𝑡 -1) -𝐸 𝑐 (𝑡)

It is required to find the total energy consumption of a node in the treatment of one data packet information. The total energy consumed by a wireless sensor node includes sensing, processing, communicating data information, switching radio model states and switching MCU modes.

 Sensing energy consumption

Sensing energy costs depends on the type of sensors. For example, the temperature sensors consumed less important energy than gas sensors. The sensor node can contain diverse sensors, and each one has its individual energy consumption attributes. Generally, the sensing energy consumption for a wireless sensor can be expressed as follows:

𝐸 𝑆 = 𝛼 × 𝑁 × 𝑉 𝑑𝑐 × 𝐼 𝑠𝑒𝑛𝑠 × 𝑇 𝑠𝑒𝑛𝑠 ( 4
)
where 𝐼 𝑠𝑒𝑛𝑠 is the required amount of current, and 𝑇 𝑠𝑒𝑛𝑠 is the duration to detect and collect 𝑁 bits data of information and 𝛼 is the percentage to compress 𝑁 bits sensing data information.

 Processing energy consumption

The sensor consumes energy to read the data message and to write it in its memory. The processing energy consumption could be calculated by:

𝐸 𝑃 = 𝛼 × 𝑁 × 𝑉 𝑑𝑐 8 × (𝐼 𝑊𝑟𝑖𝑡𝑒 × 𝑇 𝑊𝑟𝑖𝑡𝑒 + 𝐼 𝑅𝑒𝑎𝑑 × 𝑇 𝑅𝑒𝑎𝑑 ) (5)
where 𝐼 𝑊𝑟𝑖𝑡𝑒 and 𝐼 𝑅𝑒𝑎𝑑 are the necessary amount current to write and read one byte data. 𝑇 𝑊𝑟𝑖𝑡𝑒 and 𝑇 𝑅𝑒𝑎𝑑 are the necessary duration to treat the 𝐿(𝑆 𝑖 ) data information.

 Communicating energy consumption

The energy consumed to transmit and receive data messages is computed following the first-order wireless communication model for the radio hardware illustrated in Figure Fig. 3 [START_REF] Heinzelman | Energy-Efficient Communication Protocol for Wireless Microsensor Networks[END_REF].

Figure 3: First order radio transceiver model

Transmitter expends energy to turn on the radio electronics and the power amplifier. The energy required to transmit 𝐿 bits data message is: 4 𝑤ℎ𝑒𝑛 𝑑 > 𝑑 0 [START_REF] Narita | A Review on Piezoelectric, Magnetostrictive, and Magnetoelectric Materials and Device Technologies for Energy Harvesting Applications[END_REF] where 𝐸 𝑒𝑙𝑒𝑐 represents the energy consumed to transmit or receive 1 bit message, the constants 𝐸 𝑓𝑠 and 𝐸 𝑚𝑝 depend on the transmitter amplifier model. 𝐸 𝑓𝑠 is for the free space model, 𝐸 𝑚𝑝 is for multipath model, 𝑑 is the distance between the transmitter and the receiver and 𝑑 0 is a threshold distance calculated as follows:

𝐸 𝑇 (𝑑) 𝛼 × 𝑁 × 𝐸 𝑒𝑙𝑒𝑐 𝐸 𝑅 (𝑑) d 𝛼 × 𝑁 × 𝐸 𝑒𝑙𝑒𝑐 𝐸 𝑇 = { 𝛼 × 𝑁 × 𝐸 𝑒𝑙𝑒𝑐 + 𝛼 × 𝑁 × 𝐸 𝑓𝑠 × 𝑑 2 𝑤ℎ𝑒𝑛 𝑑 < 𝑑 0 𝛼 × 𝑁 × 𝐸 𝑒𝑙𝑒𝑐 + 𝛼 × 𝑁 × 𝐸 𝑚𝑝 × 𝑑
𝑑 0 = √ 𝐸 𝑓𝑠 𝐸 𝑚𝑝 ⁄ ( 7
)
and the energy consumed by the radio to receive 𝐿(𝑆 𝑖 ) bits data information is defined by:

𝐸 𝑅 = 𝛼 × 𝑁 × 𝐸 𝑒𝑙𝑒𝑐 (8)

 Switching the microcontroller (MCU) mode energy consumption

The sensor wastes energy by switching between the MCU modes. In this paper, we just take into consideration the on mode and the off mode. The energy cost for computational MCU mode can be expressed as:

𝐸 𝑆𝑤𝑖𝑡𝑐ℎ-𝑀𝐶𝑈 = 𝑉 𝑑𝑐 × (𝐼 𝐴𝑐𝑡𝑖𝑣𝑒 × 𝑇 𝐴𝑐𝑡𝑖𝑣𝑒 + 𝐼 𝑆𝑙𝑒𝑒𝑝 × 𝑇 𝑆𝑙𝑒𝑒𝑝 ) (9)

IV. PIEZOELECTRIC ENERGY HARVESTING MODEL

In this section, the piezoelectric energy harvesting model is presented and then a PID controller is adapted to the model to tune the energy harvesting and the stability of the system. The mechanical characteristics of a piezoelectric heap can be represented by an equivalent Mass-Spring-Damper (MSD) model with one degree of freedom based on the proposed model in [START_REF] Williams | Analysis of a micro-electric generator for microsystems[END_REF]. It consists of an equivalent mass 𝑀 𝑝𝑧𝑡 linked with a spring that has a constant coefficient 𝐾 𝑝𝑧𝑡 and a damper with a damping coefficient 𝐷 𝑝𝑧𝑡 and input displacement of the frame 𝑦(𝑡). The MSD model of the piezoelectric heap is presented in Figure Fig. 5.

Figure 5: Equivalent MSD model

The displacement 𝑧(𝑡) of the mass relative to the frame is characterized by the differential equation ( 3), with the mass 𝑀 𝑝𝑧𝑡 , the spring constant; the equivalent short circuit stiffness; 𝐾 𝑝𝑧𝑡 , the damping coefficient 𝐷 𝑝𝑧𝑡 and the input displacement of the frame y(t). 𝑉 is the voltage across the load resistance, 𝐶 is the capacitance, 𝑅 is the resistance and 𝑖 is the electric current. The governing equation can be found by applying Newton's second law.

𝑀 𝑝𝑧𝑡 𝑧(𝑡) + 𝐷 𝑝𝑧𝑡 𝑧(𝑡) + 𝐾 𝑝𝑧𝑡 𝑧(𝑡) = 𝑓(𝑡) [START_REF] Wang | Energy efficient distributed compressed data gathering for sensor networks[END_REF] where 𝑓(𝑡) = -𝑀 𝑝𝑧𝑡 𝑦̈ [START_REF] Leinonen | Distributed Distortion-Rate Optimized Compressed Sensing in Wireless Sensor Networks[END_REF] We can present the equivalence mechanical equation of ( 13) by the electrical equation by applying Kirchhoff's law for an RLC circuit:

𝐿𝑞̈(𝑡) + 𝑅𝑞̇(𝑡) + 1 𝐶 𝑞(𝑡) = 𝑉(𝑡) ( 12
)
Where 𝐿 is the inductance, 𝑅 is the resistance, 𝐶 is the battery's capacitance and 𝑞(𝑡) is the quantity of electricity generated by the electromotive force voltage 𝑉(𝑡).

The electro-mechanical analogy of our system is represented as follows in Table 2: Applying the Laplace transform on equation ( 13) with zero initial conditions:

𝑍(𝑠) = ℒ{𝑧(𝑡)} 𝑎𝑛𝑑 𝐹(𝑠) = ℒ{𝑓(𝑡)} (13)
𝑀 𝑝𝑧𝑡 𝑠 2 𝑍(𝑠) + 𝐷 𝑝𝑧𝑡 𝑠𝑍(𝑠) + 𝐾 𝑝𝑧𝑡 𝑍(𝑠) = 𝐹(𝑠) (

The open-loop plant transfer function between the input vibration force; 𝐹(𝑠) ant the output displacement 𝑍(𝑠) is given by:

𝑃(𝑠) = 𝑍(𝑠) 𝐹(𝑠) (15)
𝑃(𝑠) = 1

(𝑀 𝑝𝑧𝑡 𝑠 2 + 𝐷 𝑝𝑧𝑡 𝑠 + 𝐾 𝑝𝑧𝑡 ) (16)
Figure Fig. 6 shows the displacement of the mass 𝑀 𝑝𝑧𝑡 produced by an external vibration effort 𝐹. The DC gain (i.e, the amplitude ratio between the steady state response and the step input) is 1/𝐾 𝑝𝑧𝑡 , so 0.8116 is the final value of the output for a unit step input. The rise time is about 0.4 𝑠, the settling time is too long, about 6.562 𝑠 and the overshoot is high, about 50.6432.

Figure 6: Transfer function between the input vibration force and the output displacement for an open loop without the PID controller

For the closed-loop, there are three transfer functions between the setpoint and the output under the PID controller and they are given by:  The sensitivity function:

𝐸(𝑠) = 1 1 + 𝑃(𝑠)𝑈(𝑠) (17)
 The complementary sensitivity function:

𝑇(𝑠) = 𝑃(𝑠)𝑈(𝑠) 1 + 𝑃(𝑠)𝑈(𝑠) (18)
 The input sensitivity function:

𝐼(𝑠) = 𝑈(𝑠) 1 + 𝑃(𝑠)𝑈(𝑠) ( 20
)
Proportional Integral Derivative (PID) controllers are extensively used to provide robustness and optimal performance for stable, unstable, and nonlinear processes [START_REF] Ang | PID Control System Analysis, Design,and Technology[END_REF]. It can be simply implementable in analog or digital form. Further, it supports tuning and online retuning based on the performance requirement of the process to be controlled. The output of a PID controller is calculated in the time domain from the response error and it can be expressed mathematically by the equation ( 21) as follows:

𝑢(𝑡) = 𝑘 𝑑 𝑑𝑒(𝑡) 𝑑𝑡 + 𝑘 𝑝 𝑒(𝑡) + 𝑘 𝑖 ∫ 𝑒(𝑡) 𝑇 0 ( 21
)
where 𝑢(𝑡) and 𝑒(𝑡) are the control and the error signals respectively, 𝑘 𝑑 , 𝑘 𝑝 and 𝑘 𝑖 are the parameters to be tuned.

Taking the Laplace transform of equation ( 21), the transfer function for a PID controller can be expressed as follows:

𝑈(𝑠) = 𝑘 𝑑 𝑠 2 + 𝑘 𝑝 𝑠 + 𝑘 𝑖 𝑠 ( 22
)
In this section, tuning the energy harvested by the piezoelectric system is indispensable to control the quantity of energy harvested, the stability of the system. Therefore, a Proportional Integral Derivative controller (PID) is adapted to the model represented above. The PID controller uses a control loop feedback mechanism to control the piezoelectric energy harvesting system. The parameters of the PID controller are the proportional gain 𝑘 𝑝 , the integral parameter 𝑘 𝑖 and the derivative parameter 𝑘 𝑑 that affect the control of the system. The overall effects of controller parameters 𝑘 𝑝 , 𝑘 𝑖 and 𝑘 𝑑 on a closed-loop system are summarized in Table 3.

Improving proportional gain 𝑘 𝑝 has the effect of equivalently increasing the control signal for the same level of error. The fact that the controller will drive harder for a given level of error lets the closed-loop system react more speedily, but also to overshoot more. In addition, increasing 𝑘 𝑝 helps to decrease the steady-state error.

The addition of the integral parameter 𝑘 𝑖 will continually increase over time to drive the steady-state error to attain zero. If there is a steady error, the integral response will slowly increase the control signal to the error down. However, the integral term can make the system more sluggish since when the error signal changes sign. It can take a while for the integral windup phenomenon when integral action saturates a controller without the controller driving the error signal toward zero.

The derivative controller parameter 𝑘 𝑑 appends the ability of the control signal to become large, if the error begins sloping upward, even while the amplitude of the error is still quite small. This anticipation helps to add damping to the system, thereby decreasing overshoot. However, the addition of this parameter will not affect the steady-state error. There are different methods to find the PID controller parameters that satisfy these goals. The fast rise time, the minimal overshoot, the minimal steady-state error and the minimal settling time. The gains of a PID controller can be found by trial and error method. Different methods have been proposed for setting the PID controller parameters. Some of these methods are based on characterizing the dynamic response of the dynamic system to be controlled with a first-order model or second-order model with a time delay. General methods for control design can be applied to PID control. A number of special methods that are tailor made for PID control have also been developed, these methods are often called tuning methods. The most well-known tuning methods are those that are specified by Ziegler and Nichols. The Ziegler-Nichols Oscillation Method, Ziegler-Nichol Process Reaction Method and Frequency Response method, is basic Self-Tuning methods. Oscillation method is based on system gain, in other words, system gain is redounded until the system makes oscillation, then PID parameters can be found from system response graphic. In our work, we apply Ziegler-Nichols heuristic tuning method to find the PID controller parameters [19] [20]. It consists to find the ultimate gain 𝐾 𝑢 and the period 𝑇 𝑢 by setting the Integral and Derivative gains to zero. The Proportional gain 𝐺 increases until it reaches the ultimate gain, at which the output displacement of the control loop has consistent oscillations. For a classic PID controller, we can find 𝑘 𝑝 , 𝑘 𝑖 and 𝑘 𝑑 by the relations as in Table 4: The Flow Chart as shown in Figure All the values of column 1 of the table are positive. There is no sign changed. Therefore, the system controlled by pour PID controller is stable for the closed loop. The response of the MSD to a unit step input without a PID controller is shown in Figure Fig. 6, where it is observed that the MSD system without a PID controller is oscillating in nature. It is clear that the settling time is about 6.562 𝑠 and the peak overshoot is about 50.6432 𝑝. 𝑢. The rise time is about 0.4 𝑠 and the settling time is too long, about 6.562 𝑠. In addition, the terminal voltage reaches a steady value of 0.8116 𝑝. 𝑢. at about 9,2 𝑠, which is not acceptable. We can enhance the MSD system by using a controller. The using of a proportional controller can reduce the rise time of a response, but the steady-state error cannot be removed. An integral controller can improve the steady-state performance, but it can also badly affect the transient response. A derivative controller increases the performance of the transient response by reducing overshoot, thereby improving the stability margin for the system. The behavior of the MSD system can be improved by including a PID controller as shown in Figure Fig. 8.

V. SIMULATION RESULTS AND PERFORMANCE

The values of the hardware parameters used in our simulations are those of Mica2 Motes. All parameters used in our energy model are listed in Table 5. We also indicate the references where these values originated. Our results are based on existing device parameters. Moreover, that is why they can reflect a real energy model for a type of wireless sensors. Figure Fig. 9 shows the residual energy in the sensor's battery. When the battery is out of energy, the MCU switches to off mode. During the off mode, the battery is charged with a vibration harvesting energy without the adaptive PID controller. The battery is recharged to 0.1245 𝐽𝑜𝑢𝑙𝑒𝑠 and cannot reach more than this value. Therefore, there is a need to include an adaptive PID controller to the system that can enhance the stability and the performance of the system. The results illustrate that the system is still cyclical stable even if the energy consumption depends on different sensor's activities and it is not the same for each time period. Therefore, the self-power manager system must be tuned by an adaptive PID controller, to recharge the battery spontaneously (according to a critical threshold value), and stabilize the system whatever are the sensor's activities.

VI. CONCLUSION

An adaptive energy manager is essential to harvest energy in WSNs. Its core goal is to control the energy harvesting to maximize the lifetime of the sensor nodes and give the sensors the power to achieve their activities. In this paper, autonomous power management with a PID controller provides practical adaptation to the harvested and consumed energy for a sensor node. Moreover, an adaptive controller for the MCU modes, to extend a sensor lifetime depending on the existing vibration in the environment, is presented. Future work will cover a combination between simulation results and real experimental validation of the proposed autonomous energy management system. In addition, we will focus our work on the integration of some existing management energy efficient methods in our self-organized power manager system to save more energy during communications and to reduce the rapid consumption of residual energy in the sensor's battery and from the environmental energy resources. Applying energy efficient protocols for sensors' communications with the autonomous energy management system will increase the overall WSN lifetime.

For

  the piezoelectric energy harvesting model, the piezoelectric layers connected on a cantilever beam can be considered as a simple energy-harvesting device. The dimensions of piezoelectric, schematic and coordinate directions are shown in Figure Fig. 4. When the cantilever vibrates, at the first mode, the force which proceeds on the piezoelectric layers can be simplified to a 1-D model and regarded as a force 𝑓 acting on the lateral surface.

Figure 4 :

 4 Figure 4: Schematic of piezoelectric harvesting structure

Fig. 7

 7 explains how to find the PID controller parameters based on Ziegler-Nichols heuristic tuning method.

Figure 7 :

 7 Figure 7: Flowchart to find the transfer function for the PID controller based on Ziegler-Nichols heuristic tuning method

Figure Fig. 8

 8 Figure Fig. 8 shows the displacement of the mass 𝑀 𝑝𝑧𝑡 produced by an external vibration effort 𝐹 controlled by our PID controller. This controller is designed based on Ziegler-Nichols tuning method. The new system has a fast rise time 0.0311 𝑠, it reduces also the settling time 0.0509 𝑠, and minimize the overshoot 1.1887.

Figure 8 :

 8 Figure 8: Transfer function between the input vibration force and the output displacement for a closed-loop with an adaptive PID controller based on Ziegler-Nichols tuning method.

Figure 9 :

 9 Figure 9: Self-power manager behavior without the adaptive PID controller Figure Fig.10shows the residual energy in the sensor's battery. When the residual energy is less than a given threshold value (e.g.20%,) of the initial energy 0.5 𝐽𝑜𝑢𝑙𝑒𝑠 with 𝛼 = 1, the MCU switches to off mode. During the off mode, the battery is charged with a vibration harvesting energy controlled by the adaptive PID controller. The battery is recharged to a fixed threshold 0.3249 𝐽𝑜𝑢𝑙𝑒𝑠. The results illustrate that the system is still cyclical stable even if the energy consumption depends on different sensor's activities and it is not the same for each time period.

Figure 10 :

 10 Figure 10: The self-power manager controller organizes its wake-up mode to do its activities and its off mode to harvest energy controlled by the adaptive PID controller depending on its residual energy The self-power manager system including the PID controller can recharge 65% of the battery as shown in Figure Fig. 9 while the system without the adaptive PID controller cannot charge more than 25% at the same period of time as shown in Figure fig. 9.Therefore, the self-power manager system must be tuned by an adaptive PID controller, to recharge the battery spontaneously (according to a critical threshold value), and stabilize the system whatever are the sensor's activities.

TABLE Ⅱ :

 Ⅱ Electro-mechanical analogy of piezoelectric energy harvesting system

	Mechanical system	Electrical system
	Position 𝑧(𝑡)	Electrical charge 𝑞(𝑡)
	Mass 𝑀 𝑝𝑧𝑡	Inductance 𝐿
	Damping coefficient 𝐷 𝑝𝑧𝑡	Resistance 𝑅
	Spring constant 𝐾 𝑝𝑧𝑡	Inverse of capacitance 𝐶
	External force 𝑓(𝑡)	Voltage 𝑉(𝑡)

TABLE Ⅲ :

 Ⅲ The effects of PID controller parameters on a closed-loop system

	PID	Rise	Oversho	Settling	St Steady-State
	parameters	Time	ot	time		error
	𝑘 𝑝	Decrease	Increase	Small changes		Decrease
	𝑘 𝑖	Decrease	Increase	Increase		Decrease
	𝑘 𝑑	Small changes	Decrease	Decrease	No change

TABLE Ⅳ :

 Ⅳ The PID controller parameters according to Ziegler-Nichols tuning method on a closed-loop system

	𝒌 𝒑	𝒌 𝒊	𝒌 𝒅
	0.6𝐾 𝑢	1.2	𝑇 𝑢 𝐾 𝑢	0.075𝐾 𝑢 𝑇 𝑢

TABLE Ⅴ :

 Ⅴ Simulation parameters

	Parameter The supply voltage to the sensor 𝑽 𝒅𝒄 (V) Reference Energy RE (J) Threshold energy Energy dissipation: electronics 𝐸 𝑒𝑙𝑒𝑐 (nJ/bit) parameters of amplifier energy consumption 𝐸𝑚𝑝 (pJ/bit/m4) and 𝐸𝑓𝑠 (pJ/bit/m2) Data aggregation energy (J) Size of the data packet (bits) 𝑁 Compression percentage (%) 𝛼 Mass 𝑀 𝑝𝑧𝑡 (g) The damping ratio of the structure 𝐷 Parameter value The supply voltage to the sensor 𝑽 𝒅𝒄 (V) 2.7 Spring stiffness constant 𝐾 𝑝𝑧𝑡 (N. m -1 ) 0.12320 Proportional controller parameter 𝐾 𝑝 20.2366 Integral controller parameter 𝐾 𝑖 10.3729 Derivative controller parameter 𝐾 𝑑 9.8699	value 2.7 0.2 0.1 50 0.0013 and 10 5×10 -12 4000 20 182

Sensing time for one bit of data 𝑇 𝑆𝑒𝑛𝑠 (ms) 0.5

Sensing current 𝐼 𝑆𝑒𝑛𝑠 (mA) 25

Writing current 𝐼 𝑊𝑟𝑖𝑡𝑒 (mA) 18.4

Writing time for one bit 𝑇 𝑊𝑟𝑖𝑡𝑒 (ms) 12.9

Reading current 𝐼 𝑅𝑒𝑎𝑑 (mA) 6.2

Reading time for one bit 𝑇 𝑅𝑒𝑎𝑑 (µs) 565

On mode for the MicroController Unit (MCU) current 𝐼 𝐴𝑐𝑡𝑖𝑣𝑒 (mA) 8

On mode for the MicroController Unit (MCU) time 𝑇 𝐴𝑐𝑡𝑖𝑣𝑒 (ms) 1

Off mode for the MicroController Unit (MCU) current 𝐼 𝑆𝑙𝑒𝑒𝑝 (µA) 1

Off mode for the MicroController Unit (MCU) time 𝑇 𝑆𝑙𝑒𝑒𝑝 (ms) 299

Initial energy (J) 𝐸 0 0.5 𝛼 𝑝𝑧𝑡 (N. m -1 . s -1 ) 0.2

ACKNOWLEDGMENT

The authors would like to thank the editors and reviewers for their valuable comments on earlier versions of this paper.


Pierre Mahé
email: pierre.mahe@orange.com

Stéphane Ragot
email: stephane.ragot@orange.com

Sylvain Marchand
email: sylvain.marchand@univ-lr.fr

FIRST-ORDER AMBISONIC CODING WITH PCA MATRIXING AND QUATERNION-BASED INTERPOLATION

We present a spatial audio coding method which can extend existing speech/audio codecs, such as EVS or Opus, to represent first-order ambisonic (FOA) signals at low bit rates. The proposed method is based on principal component analysis (PCA) to decorrelate ambisonic components prior to multi-mono coding. The PCA rotation matrices are quantized in the generalized Euler angle domain; they are interpolated in quaternion domain to avoid discontinuities between successive signal blocks. We also describe an adaptive bit allocation algorithm for an optimized multi-mono coding of principal components. A subjective evaluation using the MUSHRA methodology is presented to compare the performance of the proposed method with naive multi-mono coding using a fixed bit allocation. Results show significant quality improvements at bit rates in the range of 52.8 kbit/s (4 × 13.2) to 97.6 kbit/s (4 × 24.4) using the EVS codec.

INTRODUCTION

Conversational applications such as telephony are typically limited to mono, with no spatial representation of the sound scene. With the emergence of new applications such as virtual reality (VR) and extended reality (XR) and the availability of devices supporting spatial audio capture and playback, there is a need to extend traditional speech/audio codecs to enable immersive communication. There are currently different spatial audio codecs developed for non-conversational applications (streaming, broadcast, etc.), including discrete coding of individual audio channels by MPEG AAC or HE-AAC, MPEG Surround [1], Dolby AC-3 or E-AC-3 [START_REF] Standard | Digital Audio Compression (AC-3, E-AC-3[END_REF], and more recently MPEG-H 3D Audio [START_REF] Herre | MPEG-H audio-the new standard for universal spatial/3D audio coding[END_REF], Dolby AC-4 [4] or DTS-UHD [START_REF][END_REF]. The most recent spatial audio codecs can handle various input formats, such as multichannel audio, object-based audio, ambisonics (also called scene-based audio [START_REF] Herre | MPEG-H audio-the new standard for universal spatial/3D audio coding[END_REF]), and multiple playback formats (e.g. mono, stereo, binaural audio, various multichannel loudspeaker setups). In this work, we investigate how spatial audio can be provided in conversational applications by extending codecs currently used in telephony or voice over IP (VoIP).

We focus in particular on reusing the Enhanced Voice Services (EVS) codec [START_REF] Bruhn | Standardization of the new 3GPP EVS codec[END_REF] which represents the state-of-art audio quality for mobile telephony applications. The EVS codec supports only mono input and output signals; a naive approach to extend this codec to spatial audio is to code each input channel by a sepa-rate instance of the mono codec -this approach is later referred to as multi-mono coding. Note that the Opus codec [START_REF] Valin | Definition of the Opus Audio Codec[END_REF] supports mono, stereo, multichannel, and recently it has been extended to code ambisonic signals, using multi-mono coding or fixed channel matrixing followed by multi-stereo coding [START_REF] Skoglund | Ambisonics in an Ogg Opus Container[END_REF].

Several methods have been proposed to code ambisonic signals. One approach is to use a parametric model of the sound field, with an assumption on the number of audio sources in the scene. For First-Order Ambisonic (FOA) coding, DirAC [START_REF] Pulkki | First-order directional audio coding (DirAC)[END_REF]Chap.5] assumes that there is a single predominant audio source in each timefrequency tile. At the encoder side, a mono downmix signal representing the predominant source is extracted together with directional parameters -in terms of direction of arrival (DoA) and diffuseness -in each time/frequency tile. The mono signal and the parameters are coded and transmitted. At the decoder side, the sound field is reconstructed by panning the predominant source (based on the DoA, typically using VBAP); this is signal is combined with an ambiance using decorrelation and the diffuseness parameters. The DirAC method has been extended to High-Order Ambisonics (HOA) in the so-called HO-DirAC [START_REF] Pulkki | First-order directional audio coding (DirAC)[END_REF]Chap.6] where the sound field is divided into angular sectors. For each angular sector, one source is extracted. More recently, Compass [START_REF] Politis | Compass: Coding and multidirectional parameterization of ambisonic sound scenes[END_REF] was proposed as a method inspired by DirAC. This method overcomes the limitation of the number of sources. The number of sources is estimated and sources are extracted by Principal Component Analysis (PCA) and the residual signal is downmixed to obtain an ambiance signal.

The Compass method is similar to the MPEG-H 3D Audio codec [START_REF] Herre | MPEG-H audio-the new standard for universal spatial/3D audio coding[END_REF] which supports ambisonic signals as one type of input format. In MPEG-H 3D Audio, the input ambisonic signal is decomposed into a number of component signals; these signals represent a number of predominant sources with an ambiance signal, they are coded using a core codec derived from MPEG USAC [START_REF] Neuendorf | The ISO/MPEG unified speech and audio coding standard-consistent high quality for all content types and at all bit rates[END_REF]. Predominant sources may be extracted using plane wave decomposition or Principal Component Analysis (PCA). When PCA is used -noting that PCA may be equivalently implemented using a Singular Value Decomposition (SVD) -components may change dramatically between consecutive frames causing channel permutations and signal discontinuities [START_REF] Zamani | Frequency domain singular value decomposition for efficient spatial audio coding[END_REF] for complex sound scenes (with many audio sources, sudden changes, etc.). The MPEG-H 3D audio codec employs channel re-alignment, overlap-add and linear interpolation to mitigate these problems. Improvements to the MPEG-H 3D Audio codec were proposed in [START_REF] Zamani | Frequency domain singular value decomposition for efficient spatial audio coding[END_REF][START_REF] Zamani | Spatial Audio Coding with Backward-Adaptive Singular Value Decomposition[END_REF]; the SVD decomposition is done in the MDCT domain to ensure smooth transitions across frames.

In this work, we wanted to avoid any assumption on the audio scene (e.g. presence of predominant sources, number of sources) based on existing mono and stereo codecs. The most basic approach is multi-mono coding; at low bit rates, due to the corre-lation of ambisonic components, the signal structure is typically degraded, which creates several spatial artifacts: spatial reduction, blur, phantom source. Multi-stereo coding with fixed matrixing was proposed in [START_REF] Skoglund | Ambisonics in an Ogg Opus Container[END_REF][START_REF]Encoding First-Order Ambisonics with HE-AAC[END_REF]. The principle consists in combining the input B-format channels by matrixing, which may be interpreted as beamforming in configurable directions. This approach has the advantage to transform the signal into a format which is more suited to multi-stereo coding, however the fixed matrixing is in general not optimal. An alternative approach proposed in [START_REF] Mcgrath | Immersive Audio Coding for Virtual Reality Using a Metadata-assisted Extension of the 3GPP EVS Codec[END_REF] is based on predictive coding of 1st-order FOA components based on the 0thorder component, which is typically efficient when the scene has few sources. In this paper, we propose an extension of multi-mono coding approach to represent FOA signals. We make no assumption on the scene content. The proposed method relies on PCA to obtain an adaptive matrixing of FOA components in time domain. To guarantee signal continuity between components across frames, two mechanisms are implemented: principal component matching (similar to MPEG-H 3D Audio and [START_REF] Zamani | Frequency domain singular value decomposition for efficient spatial audio coding[END_REF][START_REF] Zamani | Spatial Audio Coding with Backward-Adaptive Singular Value Decomposition[END_REF]) and rotation matrix interpolation in quaternion domain. We also describe an adaptive bit allocation algorithm for an optimized multi-mono coding of principal components. This paper is organized as follows. Section 2 gives a brief overview of ambisonics. Section 3 provides the relevant background on quaternions, the representation of 4D rotation matrices by double quaternions, and the interpolation of rotation matrices in quaternion domain. Section 4 describes in details the proposed coding method. Section 5 presents the results of a subjective quality evaluation comparing the proposed method and naive multimono coding, before concluding in Section 6.

AMBISONICS

Ambisonics is based on a decomposition of the sound field into an orthogonal basis of spherical harmonics. Initially limited to 1st order by Gerzon [START_REF] Gerzon | Periphony: With-height sound reproduction[END_REF], the formalism was extended to high orders by Daniel [START_REF] Daniel | Représentation de champs acoustiques, application à la transmission et à la reproduction de scènes sonores complexes dans un contexte multimédia[END_REF]. The sound field may be expressed by the equation [START_REF] Rafaely | Fundamentals of spherical array processing[END_REF]:

p(t, f, r, θ, φ) = ∞ n=0 n m=-n pnm(t, f )Ynm(θ, φ) (1)
where p(t, f, r, θ, φ) is the sound pressure at time t, frequency f , distance r, azimuth θ and elevation φ, Ynm(., .) is the spherical harmonic function, and pnm(t, f ) is the ambisonic coefficients of order n and degree m. In practice, the sound field representation is truncated to a finite order N . The higher N , the better the sound field representation accuracy. The so-called B-format corresponds to the ambisonic components pnm(t, f ). For a given order N the number of ambisonic components is (N +1) 2 in 3D and 2N +1 in 2D (in horizontal-only ambisonics). From the B-format, rendering is required to reproduce the sound field to the listener. For loudspeaker configurations the rendering computes the signal played by each loudspeaker. For headphone listening binaural rendering may be performed by decoding over virtual loudspeakers, convolving and combining the resulting feed signals by Head-Related Impulse Responses (HRIRs) [START_REF] Møller | Head-related transfer functions of human subjects[END_REF] which are filters measured for each ear.

QUATERNIONS AND 4D ROTATION MATRICES

Quaternions (also called hypercomplex numbers or hamiltonians) were introduced in 1843 by Hamilton [START_REF] Hamilton | On a new species of imaginary quantities connected with a theory of quaternions[END_REF] to define a vectorial system generalizing complex numbers. A quaternion q is defined as q = a + bi + cj + dk, where (a, b, c, d) ∈ R 4 , with the following rules for i, j, k: i 2 = j 2 = k 2 = ijk = -1. We can also write q = a + q, where a is the scalar (or real) part of q and q = bi + cj + dk is the vector (or imaginary) part of q. Note that q may be interpreted as a 3D vector by identifying i, j and k to three orthogonal Cartesian unit vectors. A quaternion q having zero scalar part is called a pure quaternion. We do not review basic operations on quaternions (e.g. addition q1 + q2, multiplication q1q2, dot product q1 • q2, conjugate q and inverse q -1 ) -see for instance [START_REF] Hamilton | On a new species of imaginary quantities connected with a theory of quaternions[END_REF][START_REF] Hanson | Visualizing Quaternions[END_REF][START_REF] Casteljau | Les quaternions[END_REF] for more details. We recall that the dot product between two quaternions q1 = a1 +b1i+c1j+d1k and q2 = a2 + b2i+c2j+d2k is q1 •q2 = a1a2 +b1b2 +c1c2 +d1d2 which corresponds to the usual dot product of 4-dimensional vectors, and the norm |q| of a quaternion is:

|q| = √ q • q = √ a 2 + b 2 + c 2 + d 2 . When |q| = 1, q is said to be a unit-norm quaternion.
Quaternions are often used as a representations of 3D rotations. We recall that the set of 3D rotations (called special orthogonal group in dimension 3 or SO(3)) can be mapped to the unit sphere in R 4 under a one-to-two mapping [START_REF] Hanson | Visualizing Quaternions[END_REF][START_REF] Casteljau | Les quaternions[END_REF][START_REF] Shoemake | Animating rotation with quaternion curves[END_REF], namely each 3D rotation matrix maps to two antipodal unit-norm quaternions: q and -q, therefore this mapping is not unique. A 3D rotation of angle θ and normalized axis n (with |n| = 1) is represented by the unit quaternion q = cos(θ/2) + sin(θ/2)n (or its opposite -q); the rotation of a 3D vector u to v can be equivalently obtained by quaternion multiplication: v = quq -1 = quq, if u and u are interpreted as pure quaternions.

Unlike other representations Euler angles or axis-angle, quaternions are interesting to interpolate rotations. One common method is called spherical linear interpolation (slerp) and consists in the following principle [START_REF] Hanson | Visualizing Quaternions[END_REF]:

slerp(q1, q2, γ) = q1(q -1 1 q2) γ (2)
where q1 and q2 are respectively the starting and ending quaternions and 0 ≤ γ ≤ 1 is the interpolation factor. The slerp interpolation may also be formulated as [START_REF] Hanson | Visualizing Quaternions[END_REF]:

slerp(q1, q2, γ) = sin((1 -γ)Ω) sin(Ω) q1 + sin(γΩ) sin(Ω) q2 (3)
where Ω = arccos(q1 • q2) is the angle between the two quaternions q1 and q2. This interpolation boils down to interpolating along the grand circle (or geodesics) on the unit 4D sphere with a constant angular speed as a function of γ. Due to the non-unique representation of 3D rotations by antipodal unit-norm quaternions, one has to ensure that the quaternion trajectory follows the shortest path [START_REF] Shoemake | Animating rotation with quaternion curves[END_REF] on the unit sphere in R 4 by selecting among the two possible choices ±q2 for the end point. Alternative interpolation methods, such as normalized linear interpolation (nlerp) or splines [START_REF] Hanson | Visualizing Quaternions[END_REF], are not reviewed here.

In this work, we used the representation of 4D rotation matrices as double quaternions. The multiplication of two quaternion matrices, an anti-quaternion matrix Q and a quaternion matrix P, associated to two quaternions q = a + bi + cj + dk and p = w + xi + yj + zk, is defined as [START_REF] Casteljau | Les quaternions[END_REF]:

R = Q .P = P.Q (4)
where

Q =      a b c d -b a -d c -c d a -b -d -c b a      (5)
and

P =      w -x -y -z x w -z y y z w -x z -y x w      (6)
It is possible [START_REF] Casteljau | Les quaternions[END_REF] to show that the product R has the properties of a rotation matrix (R.R T = R T .R = I and det(R) = +1); conversely, given a 4D rotation matrix R, it is possible [START_REF] Casteljau | Les quaternions[END_REF] to factorize this matrix as in Eq. 4. This factorization (sometimes called Cayley's factorization) may be obtained using the method described in [START_REF] Perez-Gracia | On Cayley's factorization of 4D rotations and applications[END_REF] (noting that a factor 1/4 is missing in Eqs. 64 and 65 in [START_REF] Perez-Gracia | On Cayley's factorization of 4D rotations and applications[END_REF]); this method is also described in [START_REF] Casteljau | Les quaternions[END_REF] (together with alternative factorization approaches). The factorization relies on an intermediate matrix (called associate matrix in [START_REF] Perez-Gracia | On Cayley's factorization of 4D rotations and applications[END_REF] or tetragonal transform in [START_REF] Casteljau | Les quaternions[END_REF]) to determine the two quaternions q and p up to sign (i.e. -q and -p also form a valid solution). Similar to 3D rotations, two 4D rotation matrices may be interpolated by interpolating separately the associated pairs of quaternions (for instance using slerp). However, it is important to keep the sign consistent between double quaternions when constraining the shortest path.

Note that there are other representations for 4D rotation matrices. In this work, we also used the generalized Euler angles defined in [START_REF] Hoffman | Generalization of Euler Angles to N-Dimensional Orthogonal Matrices[END_REF]. In general, a n-dimensional rotation matrix is characterized by n(n -1)/2 generalized Euler angles and for n = 4 we have 6 angles. We refer to [START_REF] Hoffman | Generalization of Euler Angles to N-Dimensional Orthogonal Matrices[END_REF] for details on the conversion between an n-dimensional rotation matrix and generalized Euler angles.

PROPOSED CODING METHOD

We describe in this section the proposed coding method. The input signal is assumed to be a first-order ambisonic signal with the ACN ordering convention (W, Y, Z, X). The n = 4 ambisonic components will be labeled respectively with an index i = 1, • • • , n. Each ambisonic component is sampled at 32 kHz. The coding method operates on successive frames of 20 ms. The different coding and decoding steps are described in the high-level diagram shown in Figure 1. We describe below the implementation for each functional blocks.

High-pass pre-processing (HPF)

The n = 4 channels of the input FOA signal are separately preprocessed by a 20 Hz high-pass IIR filter from the EVS codec defined by the transfer function

Hpre(z) = b0 + b1z -1 + b2z -2 1 + a1z -1 + a2z -2 (7)
where filter coefficients (bi, ai) are taken from [START_REF]Codec for Enhanced Voice Services (EVS); Detailed algorithmic description[END_REF]. This preprocessing is used here to avoid any bias in the subsequent estimation of the covariance matrix.

Principal Component Analysis (PCA)

In each frame, the sample covariance matrix C XX is estimated based on the 4-dimensional input X = [x1, ..., xn] :

C XX = 1 n -1 n i=1 xix T i ( 8
)
where xi is the i-th column of X and corresponds to the current frame of i-th ambisonic channel. The covariance matrix C XX is factorized by eigenvalue decomposition as:

C XX = VΛV T ( 9
)
where V is the eigenvector matrix (with eigenvectors as columns) and Λ = diag(λ1, • • • , λn) is the diagonal matrix whose coefficients are the eigenvalues.

In general the matrix V is orthogonal and it may be either a rotation matrix (det(V) = +1) or a reflection matrix (det(V) = -1). In this work, we ensured that the eigenvector matrix defines a rotation matrix by inverting the sign of vn if det(V) = -1. In the following, the resulting rotation matrix in the current frame of index t will be denoted Vt.

Re-alignment of eigenvectors

From frame to frame eigenvectors might change significantly. These modifications might create discontinuities in the signal, which can degrade audio quality. To improve signal continuity between successive frames, a signed permutation is applied to the eigenvector matrix Vt in the current frame of index t based on the eigenvector matrix Vt-1 in the previous frame. The signed permutation is obtained in two steps:

1. A permutation is found by matching eigenvectors in frames t and t -1 according to the axes (not directions) of each basis vector. This problem can be seen as an assignment problem where the goal is to find the closest eigenvector in frame t for each eigenvector in the previous frame t -1. To solve this assignment problem, we used the Hungarian algorithm to find the optimal solution in a similar way to [START_REF] Zamani | Frequency domain singular value decomposition for efficient spatial audio coding[END_REF]. After this step, eigenvectors are permuted. This allows to maximize similarity between the two basis. The similarity being defined as:

Jt = tr(|Vt.V T t-1 |) (10)
where tr(|.|) is the trace of the matrix |Vt.V T t-1 | whose coefficients are the absolute values. After applying the optimal permutation to Vt we obtain a new matrix of eigenvectors Ṽt.

The direction of each eigenvector is determined based on

the permuted eigenvector matrix in the current frame Ṽt for frame t and the rotation matrix Vt-1 in the previous frame t -1

Γt = Ṽt.V T t-1 (11)
A negative diagonal value in Γt indicates a direction inversion between two frames. The sign of the respective columns of Ṽt is inverted to compensate for this change of direction.

Note that the resulting matrix Ṽt in the current frame will be saved for the next frame processing to become Vt-1. Quant.

-1

Interp.

FOA

Figure 1: Overview of proposed coding method.

Quantization of the 4D PCA rotation matrix

In general a n-dimensional rotation matrix has n(n-1)/2 degrees of freedom. In [START_REF] Briand | Études d'algorithmes d'extraction des informations de spatialisation sonore : application aux formats multicanaux[END_REF] 2D and 3D rotation matrices were quantized with an angle representation: one angle in 2D, three Euler angles in 3D. In this work we used a similar idea by quantizing 6 generalized Euler angles. The 6 angles were obtained from the 4D rotation matrix [START_REF] Hoffman | Generalization of Euler Angles to N-Dimensional Orthogonal Matrices[END_REF] and coded by scalar quantization with a budget of respectively 8 and 9 bits for angles defined over a support of length π and 2π, with an overall budget of 51 bits per frame for 6 angles.

Interpolation of quantized rotation matrices by subframe

The rotation matrices are interpolated by subframes to smooth variations of PCA rotation matrices and avoid signal discontinuities. The rotation matrix is not suitable for interpolation. Instead, we converted 4D rotation matrices to double quaternions as explained in Section 3. The current frame of length L = 640 samples is divided into K sub-frames. We used K = 128 which gives a subframe length of L/K = 10 samples (0.3125 ms). For each subframe of index 1 ≤ k ≤ K in the current frame, the left (qt-1, qt) and right quaternions (pt-1, pt) are interpolated by the slerp algorithm with an interpolation factor given by γ = k/K as defined in Eq. 3. The interpolated double quaternions are converted back to a 4D matrix using Eqs. 4, 5, 6.

PCA matrixing

The pre-processed FOA signal is transformed into 4 principal components by applyinng the interpolated rotation matrix in each subframe.

Adaptive bit rate allocation to multi-mono EVS coding

After PCA matrixing the n = 4 channels could have been coded using a fixed bit allocation, as in naive multi-mono coding. Depending on the input signals, the signals after PCA matrixing may vary significantly in importance and it was found experimentally that an adaptive bit allocation is necessary to optimize quality. We use a greedy bit allocation algorithm which aims at maximizing the following score:

S ( b1, • • • , bn) = n i=1 Q(bi).E β i ( 12
)
where bi and Ei are respectively the bit allocation and the energy of the i th channel in the current frame and Q(bi) is a quality score reflecting the estimated quality of the codec for a bit rate corresponding to bi bits per 20 ms frame. This optimization is subject to the constraint b1 +• • •+bn ≤ B where B is the budget allocated for multi-mono coding. The term Q(bi) is defined here in Table 1 to reflect the fact that the EVS codec quality does not increase [START_REF]Codec for Enhanced Voice Services (EVS); Performance Characterization[END_REF] linearly with increasing bit rate and it was found experimentally the theoretical rate-distortion function that would apply for a usual source model (e.g. Gaussian source) would not be suitable. The scores Q(bi) correspond here to average MOS (Mean Opinion Score) values found during the performance characterization of the EVS codec by 3GPP [START_REF]Codec for Enhanced Voice Services (EVS); Performance Characterization[END_REF]. Note that if another core codec than EVS was used, the values Q(bi) could be adjusted accordingly; for instance, a quality evaluation of Opus can be found in [START_REF] Rämö | Voice quality characterization of IETF Opus codec[END_REF]. To guarantee maximum compatibility with the EVS codec, the bit allocation to individual audio channels is restricted to the set of EVS bit rates: 9.6, 13.2, 16.4, 24.4, 32, 48, 64, 96, 128 kbit/s. Note that a minimum bit allocation of 9.6 kbit/s was defined to ensure a super-wideband coded bandwidth. The energy term E β i is raised to the power β, where 0 ≤ β ≤ 1 is defined as an extra tuning parameter: when β is close to 1, channels with more energy will dominate the bit allocation and when β = 0 channels would receive an equal allocation. The value of β = 0.5 was selected experimentally. The bit allocation b1, • • • , bn selected in the current frame is coded and transmitted to the decoder. rate (corresponding to b1 + b2 + b3 + b4) does not always use the available budget and in this case padding bits have to be added to the bitstream.

Multi-mono EVS coding

The bit rate to code the meta-data (rotation matrices, bit allocation) is subtracted to the target bit rate and the remaining bit rate is used for the adaptive bit allocation to multi-mono coding. The transformed FOA channels are coded by separate instances of the EVS codec and the associated bitstreams are transmitted to the decoding part. In this work we used the fixed-point implementation of EVS (V14.2.0) with discontinuous transmission disabled.

A bitstream structure example is shown in Figure 3 for the bit rate of 4 × 24.4 = 97.6 kbit/s; the bitstream is divided in several sections: bit allocation (6 bits), quantized generalized Euler angles (51 bits), four coded channels (b1 + b2 + b3 + b4 bits) and padding (when needed).

Decoding part

The decoding part is similar to the encoding part. The bit allocation is demultiplexed and the 4 coded channels are decoded by separate instances of the EVS decoder. In addition, the generalized Euler angles are decoded and converted into a 4D rotation matrix. The interpolation in quaternion domain done in each subframe is performed. The signal in the current frame divided in K subframes is then transformed by inverse PCA matrixing to retrieve the reconstructed ambisonic components.

EVALUATION

Test methodology

We conducted subjective tests according to the MUSHRA methodology [START_REF]Method for the subjective assessment of intermediate quality level of coding systems[END_REF] to compare the perceptual performance of naive multimono coding and the proposed coding method. For each item, subjects were asked to evaluate the quality of conditions with an integer grading scale ranging of 0 to 100. This interval is divided in 5 sections of 20 points: bad (0-20) to excellent (80-100). The test conditions included three specific items: the hidden reference (FOA) and two anchors. Traditional MUSHRA tests for mono signals typically use a low anchor (reference processed by a 3.5kHz low-pass filter) and a medium anchor (reference processed by a 7kHz low-pass filter). For MUSHRA tests with stereo, it is suggested to use "reduced stereo image" as degradations in anchors [START_REF]Method for the subjective assessment of intermediate quality level of coding systems[END_REF]. There is no clear recommendation on spatial alterations for MUSHRA tests with ambisonic signals. In this work we used anchors with some spatial deformation. We use the following spatial reduction defined by:

F OA =      W αX αY αZ      , α ∈ [0, 1] (13)
with α = 0.65 and α = 0.8 for the low and medium anchors, respectively.

Experimental setup

The subjects were placed in an audio listening room. The room complied with ITU-R recommendations [START_REF]Methods for the subjective assessment of small impairments in audio systems including Multichannel Sound Systems[END_REF] in terms of acoustics properties, reverberation time and background noise. All subjects conducted the listening test with the same hardware, which included an external sound-card (Focusrite Scarlett 6i6) and highquality headphones (Sennheiser HD 650) and a test interface running on a MacBook pro computer.

The test items consisted of 10 challenging ambisonic items: 4 voice items, 4 music items and 2 ambient scenes. Six items were real recordings made by ambisonic microphones (EigenMike or SoundField SPS200), others were synthetic items. A description of these items can be found in Appendix 7. Each item was about 10 s long. All FOA items were binauralized with Resonance Audio renderer [START_REF]Resonance audio : Rich, immersive, audio[END_REF] using generic HRTFs corresponding to the KU100 manikin. Original FOA signals were normalized in loudness using the following procedure: the rerefence FOA signals were binauralized using the Resonance audio renderer; the resulting binaural signal was normalized to -23 dBLUFS according to ITU-R BS.1770 [START_REF]Algorithms to measure audio programme loudness and true-peak audio level[END_REF], and after verifying that there was no saturation in binaural signals, the obtained scale factor was re-applied to the respective coded FOA signal. In total 11 listeners participated in the test; all of them are expert or experienced listeners without hearing impairments. Each item was codec at three bit rates for multi-mono coding: 52.8, 65.6, 97.6 kbit/s which corresponds to a fixed allocation of 13.2, 16.4 and 24.4 kbit/s per channel. For the proposed coding method, as explained in Section 4, the bit rate was dynamically distributed between channels; however the target (maximum) bitrate was set to the same bit rate as multi-mono coding for a fair comparison. All test conditions are summarized in Table 2.

Test results

The subjective evaluation results, including the mean and 95% confidence intervals, are presented in Figure 4. They show that the proposed coding method improves quality when compared with multi-mono coding at the same bit rate. It is particularly noticeable that multi-mono coding at 65.6 kbit/s is equivalent to the proposed coding method at 52.8 kbit/s. These quality improvements are largely due to the suppression of spatial artifacts. These artifacts are present at every bit rate in multi-mono coding and they can be classified into three categories: diffuse blur, spatial centering, phantom source. These three types of artifacts in multi-mono coding result from the degraded structure between ambisonic components. With the proposed coding method, these artifacts are mostly removed because the structure is less important after PCA matrixing. This explanation was supported by the feedback from some subjects, after they conducted the subjective test.

It is possible to analyze MUSHRA scores for different item categories. Figure 5 shows the scores for recorded and synthetic scenes. As can be seen, the proposed coding method brings significant improvements for synthetic items. The result for each item is the same whatever the position of sources or the number of sources. Two assumptions may explain this result. The first one is that the sources do not interact together, and the PCA can decorrelate each source, consequently the proposed method avoided spatial artifacts. The second one is that spatialization is more pronounced in synthetic items (wider displacement, more localized source). This spatialization puts the emphasis for the listener on spatial artifacts.

H R E F L O W _ A N C H O R M E D _ A N C H O R M U L T I 5 2 M U L T I 6 5 M U L T I 9 7 P C A 5

CONCLUSION

This paper presented a spatial extension of a mono speech/audio codec. The proposing coding method operates in time domain to avoid extra delay and allow maximum compatibility with existing codecs (e.g. EVS or Opus) which are used as a black box. The ambisonic components are transformed by adaptive matrixing depending on the audio scene.

For each frame, a PCA allows to find a new basis where the components are decorrelated. To guarantee smooth transition between consecutive frames, rotation matrices are interpolated in quaternion domain. Subjective test results show that the PCA matrixing together with the adaptive bit allocation give significant improvements over naive multi-mono coding for bit rates from 4 × 13.2 to 4 × 24.4 kbit/s. For future work, it will be interesting to characterize the spatial artifacts for various coding methods and to apply the proposed method to Opus.

H R E F L O W _ A N C H O R M E D _ A N C H O R M U L T I

APPENDIX: DESCRIPTION OF MUSHRA ITEMS

Synthetic items were generated in Orange Labs, recorded items were captured and mixed by Orange Labs or done jointly with partners (see details below). HOA items were truncated to FOA for testing.

• Drums (synthetic) Drum tracks were generated with a MIDI sequencer. HOA mixing, spatialization and near-field simulation were done with an internal Matlab library. This item is a shorterned version of the H_02_Drums1 test item [START_REF]Submission and Evaluation Procedures for 3D Audio[END_REF] contributed by Orange Labs to MPEG-H standardization.

• Modern Music (synthetic) Single source (jazz music) crossing left to right with near-field effect. The mono source is an except of "The Present" in Laurent de Wilde's album "Time for Change". Spatial rendering was done with an internal Matlab library. This item is a shorterned version of the H_03_Modern test item [START_REF]Submission and Evaluation Procedures for 3D Audio[END_REF] contributed by Orange Labs to MPEG-H standardization.

• Opera (recorded) Italian Opera with female singer, harpsichord and strings in concert hall from the Ambisonia database [START_REF]Ambisonia[END_REF]. Players were in front of a SoundField SPS200 microphone.

• Orchestra (recorded) Orchestra in concert hall. An EigenMike microphone was placed in middle of the orchestra and spot microphones were also used for ambisonic mixing.

• Theater (recorded) Theater play with 3 moving talkers (in French), one in near-field and two in far-field. Recorded with an EigenMike microphone during the Bili Project [START_REF] Roncière | Compte-rendu de captation (Rapport projet BiLi) -Fiction CNSMDP : 2 femmes pour un fantôme de René de Obaldia[END_REF].

• Kids Playground (synthetic) Two talkers (female and male in English) mixed with an ambisonic recording of kids playground. The mix was done with Reaper, the mono speech signals were spatialized with the Ambix framework [START_REF] Kronlachner | ambix v0.2.8 -ambisonic plug-in suite[END_REF], the ambisonic kids playground signal came from the Ambisonia database [START_REF]Ambisonia[END_REF].

• Little Prince (recorded) Man reading excerpt of "Le Petit Prince" (in French) at a fixed (left, top) position in a living room. Recorded with an EigenMike microphone.

• Talks (recorded) People at different positions talking (in French) in a large room. Recorded with an EigenMike microphone during the Bili Project [START_REF] Roncière | Compte-rendu de captation (Rapport projet BiLi) -Fiction CNSMDP : 2 femmes pour un fantôme de René de Obaldia[END_REF].

• Nature (synthetic) Artificial scene with bumblebee sound moving in space, mixed with an ambisonic brook noise and a bird singing at fixed position. The mix was done with Reaper, the mono signals of bumblebee and bird were spatialized with the Ambix framework [START_REF] Kronlachner | ambix v0.2.8 -ambisonic plug-in suite[END_REF], the bird noise came from the Ambisonia database [START_REF]Ambisonia[END_REF].

• Applause (recorded) Applause at the end of a concert. A Sound-Field SPS200 microphone was placed in the middle of the crowd. The recording came from the Ambisonia database [START_REF]Ambisonia[END_REF].

Figure 2 :

 2 Figure 2: Bit rate allocation example at 4 × 24.4 kbit/s.

Figure 2

 2 Figure2shows an example of bit rate allocation to components (denoted Comp. 1, 2, 3, 4) in successive 20 ms frames for one test item. The maximum allocated bit rate (94.75 kbit/s) corresponds to the difference between 4 × 24.4 kbit/s and the bit rate used for side information (2.85 kbit/s). The overall allocated bit

Figure 3 :

 3 Figure 3: Bitstream structure example at 4 × 24.4kbit/s.

Figure 4 :

 4 Figure 4: MUSHRA mean scores with 95% confidence intervals.

Figure 5 :

 5 Figure 5: MUSHRA mean scores with 95% confidence intervals for different types of scenes.

Table 1 :

 1 Bit allocation parameters.

	bi	192 264 328 488 640 960 1280 1920 2560
	rate kbit/s 9.6 13.2 16.4 24.4 32 48 64	98 128
	Q(bi) 3.62 3.79 4.25 4.60 4.53 4.82 4.83 4.85 4.87

4 Coded Signals (1895 bits)

	Side information (57 bits)
	Bit allocation (6 bits)	6 Euler Angles (51 bits)	Comp. 1	Comp. 2	Comp. 3	Comp. 4 Padding

Table 2 :

 2 List of MUSHRA conditions.

	Short name	Description
	HREF	FOA hidden reference
	LOW_ANCHOR 3.5 kHz LP-filtered and spatially-reduced FOA (α = 0.65)
	MED_ANCHOR 7 kHz LP-filtered and spatially-reduced FOA (α = 0.8)
	MULTI52	FOA coded by multimono EVS at 4 × 13.2 kbit/s
	MULTI65	FOA coded by multimono EVS at 4 × 16.4 kbit/s
	MULTI97	FOA coded by multimono EVS at 4 × 24.4 kbit/s
	PCA52	FOA coded by proposed method at 52.8 kbit/s
	PCA65	FOA coded by proposed method at 65.6 kbit/s
	PCA97	FOA coded by proposed method at 97.6 kbit/s

Proceedings of the 22 nd International Conference on Digital Audio Effects (DAFx-19), Birmingham, UK, September 2-6, 2019

ACKNOWLEDGMENTS

The authors would like to thank all participants in the audio test. They also thank Jérôme Daniel for helpful discussions on test items and spatial reduction for MUSHRA anchors.