Assessment of Breathing Parameters Using an Inertial Measurement Unit (IMU)-Based System

Breathing frequency (fB) is an important vital sign that—if appropriately monitored—may help to predict clinical adverse events. Inertial sensors open the door to the development of low-cost, wearable, and easy-to-use breathing-monitoring systems. The present paper proposes a new posture-independent processing algorithm for breath-by-breath extraction of breathing temporal parameters from chest-wall inclination change signals measured using inertial measurement units. An important step of the processing algorithm is dimension reduction (DR) that allows the extraction of a single respiratory signal starting from 4-component quaternion data. Three different DR methods are proposed and compared in terms of accuracy of breathing temporal parameter estimation, in a group of healthy subjects, considering different breathing patterns and different postures; optoelectronic plethysmography was used as reference system. In this study, we found that the method based on PCA-fusion of the four quaternion components provided the best fB estimation performance in terms of mean absolute errors (<2 breaths/min), correlation (r > 0.963) and Bland–Altman Analysis, outperforming the other two methods, based on the selection of a single quaternion component, identified on the basis of spectral analysis; particularly, in supine position, results provided by PCA-based method were even better than those obtained with the ideal quaternion component, determined a posteriori as the one providing the minimum estimation error. The proposed algorithm and system were able to successfully reconstruct the respiration-induced movement, and to accurately determine the respiratory rate in an automatic, position-independent manner.


Introduction
Continuous monitoring of respiratory parameters such as breathing frequency (f B ), inspiratory time (T I ) and expiratory time (T E ) could foster early diagnosis of a wide range of respiratory disorders and help to track a patient's condition, discriminating between stable and at-risk patients [1,2]. Conditions of interest could be sleep breathing disorders, sudden infant death syndrome, chronic obstructive pulmonary disease (COPD) and neuromuscular disorders. The current gold standard for measuring f B is to count the number of breaths in one minute, through auscultation or observation [3,4]. Other methods for breathing function assessment currently used in clinical practice are spirometry or pneumotachograph based on airflow measurement by using mouthpiece or facemask. In overnight polysomnography, breathing activity is assessed both by measuring respiratory flow, through pressure transducer or thermistors near the nostrils, and respiratory efforts (breathing-derived chest-wall movements), by strain-gauge belts. Also, exhaled carbon dioxide sensors, transthoracic inductance and impedance plethysmography and ECG-or PPG-derived f B have been used to measure breathing signal. Despite their accuracy, these methods are uncomfortable and intrusive, and are not suitable for continuous monitoring in the clinical environment and at home. An emerging area of interest is to use motion sensors to detect the small breathing-derived movements/orientation changes of the chest wall. This method is particularly suitable for long-term breathing monitoring because it is unobtrusive, tolerable, and low-cost. The principle was first presented with a single-axis accelerometer in animal model (dog) using a pressure transducer in the trachea as reference [5]. Starting from this point, a variety of studies demonstrated the feasibility of using one accelerometer placed on the chest wall to derive respiratory signal and/or breathing frequency in different positions [6][7][8][9][10][11][12][13][14]. Morillo et al. [8] combined a piezoelectric single-axis accelerometer and a polarized capacitive microphone placed on the suprasternal notch to collect information of the cardiac, respiratory, and snoring activities for the screening of patients affected by Sleep Apnea-Hypopnea Syndrome. Measurements were limited to the supine position, that was selected to increase the sensitivity of the single-axis accelerometer, limiting the generality of the findings. The analysis method was based on the estimation of breathing frequency through the identification of the peak of the spectrum or autocorrelation; the main limitation of this approach is that, when the breathing is irregular, a main peak may not exist, and individual breaths must be identified and counted. Hung et al. [7] moved from single-axis to biaxial accelerometers. The aim of their study was to evaluate the reliability of the device in terms of detection of the onsets of expiration and expiration, and to assess the feasibility of differentiating between different breathing patterns (normal breathing, apnea, deep breathing). The signals from both axes (anteroposterior and longitudinal) of the accelerometer were summed, limiting the analysis to the sagittal plane, in sitting and lying positions. An adaptive band-pass filter was applied with a variable passband centered at the detected dominant breathing frequency.
As emerged by these studies, single or dual-axis accelerometers can be used to derive breathing signal when appropriately aligned with the major axis of rotation, which changes when the subject move from a posture to another. Contrarily, the use of a tri-axial accelerometer allows measuring inclination changes due to breathing regardless of orientation. In this case, the problem lies in the identification of the accelerometer axis to consider when posture changes. Bates et al. [13] proposed a method to track the major axis of rotation as it changes, to continuously monitor angular motion due to breathing also when subject change position/orientation. An alternative possibility to the best axis selection is fusing the axes. Jin et al. [12] proposed a posture-independent signal processing method based on three possible algorithms for accelerometer axes fusion. They demonstrated that methods based on Principal Component Analysis (PCA) obtained the highest performance in terms of Signal-to-Noise Ratio (SNR), but no results were provided about breathing rate estimation or validation against a reference method.
With the entry of tri-axial accelerometers new opportunities opened for the monitoring of breathing frequency using inertial sensors, but their use was still confined to static conditions since, when the subject is moving, the degree of the movement-related signal would exceed that due to breathing. One possible approach is to identify non-breathing motion, as proposed by Bates et al [13]. In a successive study, Mann et al. [11] furtherly developed the method proposed by Bates et al. [13], by adding activity tracking, and allowing identification of asymmetric breaths, that was not possible in the original method. An attempt to remove motion artifacts by using signal processing was made by Liu et al. [14]. They proposed an elegant method based on PCA-fusion of the three axes of an accelerometer and on filtering of the first principal component by using an adaptive filter that varied according to the energy expenditure derived by the same accelerometer. To overcome the problems of using a single accelerometer in dynamic conditions, a possibility is to fuse data from accelerometers and from other sensors, such as gyroscopes. Yoon et al. [15] investigated the feasibility of measuring breathing-related motions also during dynamic activities of the subject, by fusing data from a tri-axial accelerometer and a gyroscope and applying Kalman filter. They found that, during dynamic exercises, fusion of accelerometer and gyroscope data provided benefits in terms of reduction of estimation error. Gollee et al. used a more complex system, an inertial measurement unit (IMU) fusing accelerometer, gyroscope and magnetometer but considered only static conditions [16]. Another approach to overcome the problems related to motion artefacts is modularity. Lapi et al. [17] tried to overcome limitations deriving from the use of a single accelerometer by proposing a system based on a couple of 3-axis accelerometers placed bilaterally on the chest. Using two accelerometers permitted to detect respiration-related chest-wall movements regardless of sensor positioning with respect to the gravity vector; secondly, the breathing frequency can be obtained even when one of the two sensors is silenced by postural constraints. Recently, Gaidhani at al. [18] proposed a method that uses two IMUs composed by a 3-axis accelerometer, a 3-axis gyroscope and a 3-axis magnetometer, placed on the anterior and posterior side of the chest to decompose the motions experienced by the two IMUs into trunk movements and breathing actions. This paper presents an automatic processing algorithm to derive breathing frequency and other breathing temporal parameters from quaternion-based orientation signals recorded simultaneously at thoracic and abdominal level by using a modular, wireless, IMU-based device [19]. An important step of the processing algorithm is dimension reduction (DR) that allows the extraction of a single respiratory signal starting from 4-component quaternion data. Three different methods of DR are proposed and compared; two of them are based on the selection of one quaternion component, the third one is based on PCA-fusion of the 4 quaternion components. Results obtained using the IMU-based device, with the three different methods, are validated against optoelectronic plethysmography, an already established method to evaluate ventilation through an external measurement of the chest-wall surface motion [20][21][22][23][24][25].

Device Architecture and Hardware Description
The system used in this study is composed by three IMU-sensor units that communicate via Bluetooth with a smartphone; here data are pre-processed and saved. Two of the three sensor units (peripheral units) are dedicated to the recording of chest-wall respiratory-related movements and are placed on the thorax and on the abdomen to record respiratory information about both the compartments; the third sensor unit (reference central unit) is placed on a body area that is integral with the chest wall, but not involved in respiratory movements (e.g., coccyx or anterior superior iliac crest). The measurement of chest-wall movements, related to both abdominal and thoracic compartments, allows the consideration of the two-degree-of-freedom (DoF) model of chest-wall breathing movements [26], that considers abdomen and rib cage (thorax) as acting independently. Moreover, the compartmental contribution to total chest-wall volume changes according to posture and the breathing strategy adopted by each subject. Thus, the recording of chest-wall movements at different levels provides on the one hand, a more accurate estimation of the breathing signal, and on the other hand allows investigation of asynchronies between compartments, typical of different pathological conditions. The reference central unit, in addition to performing a central role within the Bluetooth piconet, can be used to discriminate between static and dynamic conditions and to map the activity state of the subject. Moreover, it could be used to reduce movement information not linked to breathing by means of frequency domain analysis or by referring orientation change experienced by the peripheral units to the coordinate frame of the reference unit. Each unit is composed by a printed circuit board equipped with a low-power microcontroller, a Bluetooth Low Energy (BLE) module, a 9-DoF IMU (3-axis accelerometer, a 3-axis gyroscope, and a 3-axis magnetometer) and lithium polymer rechargeable battery. A voltage regulator circuit, and Li-Po battery recharge circuit with mini USB port are also included in the design. Differently from the peripheral units, the reference central unit is equipped with a different BLE module, able to support simultaneous central/peripheral role and also brings a Micro Secure Digital (SD) Memory Card Connector for data logging. The dimensions of each peripheral unit, comprehensive of the 3D-printed housing, are 41 mm × 33 mm× 19 mm (LWH), Sensors 2019, 19, 88 4 of 24 and the weight is 25 g, including the battery, while the reference central unit measures 45 mm × 45 mm × 15 mm (LWH), and weighs 35 g. A prototypal version of this device has been described in [19].

Quaternion-Based Orientation Estimation and Fusion Algorithm
The final goal is to derive breathing signal by measuring orientation changes during the respiratory movements, both at thoracic and abdominal level. The IMUs provide 3D-acceleration, 3D-magnetic field, and 3D-angular rate. These measures are combined to provide accurate 3D orientation data aboard each unit. The orientation is represented with quaternions [27,28], that even though may suffer from problems of interpretation in terms of meaningfully physical angles, are interesting mathematical entities (four-dimensional complex number (q = [q 0 q 1 q 2 q 3 ]), since they require less computing time and avoid the singularity problems (i.e., "gimbal lock") typical of other orientation descriptors, e.g., Euler angles. The fusion of the data collected from the sensors is done by using the sensor fusion algorithm proposed by Madgwick et al. [29], based on an analytically derived and optimized gradient descent algorithm enabling levels of accuracy exceeding that of the Kalman-based algorithm, with low computational (277 scalar arithmetic operations each filter update) load and low sampling rates (e.g., 10 Hz); this orientation filter also provides an online magnetic distortion compensation algorithm and gyroscope bias drift compensation. The sensors data were collected at 40 Hz and the fusion algorithm was updated with the same rate, but due to limited buffer of the BLE module and to the stricter timings used for the Bluetooth communication, just one quaternion out of 4 computed is considered (10 Hz); nevertheless, the final sampling rate was considered appropriate given the relative low frequency of the respiratory signal [0.1 ÷ 1 Hz]. Thus, the microprocessor of each unit, receives data from accelerometer, gyroscope and magnetometer that are on board and implements Madgwick fusion filter [29] to compute a quaternion representing the change of orientation of each unit relative to the earth frame ( Th Earthq , Ab Earthq ,

Re f
Earthq ), or more correctly the change of orientation of the earth relative to each unit frame [29]. In fact, in quaternion form, an arbitrary orientation of a coordinate frame B relative to coordinate frame A, achieved through a rotation of angle θ around an axis A r (r x , r y , r z ) defined in frame A, is univocally represented through the normalized quaternion A Bq defined by Equation (1):

Quaternion-Derived Breathing Frequency
All the elaborations and computations needed to extract breathing parameters from data collected by the device were performed offline using MATLAB, the processing took on average 1.027 ± 0.129 seconds for the analysis of signals of 1071 ± 270 samples. A signal processing procedure was designed to extract the breathing frequency starting from quaternions representing the change of orientation of each unit relative to the earth frame ( Th Earthq , Ab Earthq , Figure 1. Block Diagram of the Analysis algorithm that allows derivation of breathing temporal parameters (fB. TI, TE) from quaternion-based orientation change signals recorded on Thorax, Abdomen and Reference point. (2) and (3), that are composed by 4 components each, and provides as output 2 single-component signals (1 for the abdomen and 1 for the thorax) representing chest-wall respiratory-related orientation change signals. These two signals represent the input of the power spectrum block and of the processing block. To reduce dimension from 4 components to 1, two possibilities were investigated as shown in Figure 2:

Dimension-reduction block takes the quaternions obtained from Equations
(i). Best quaternion component selection (ii). PCA-based fusion of the quaternion components To select the best component among the 4 components representing the orientation quaternion, two different methods were proposed, both based on spectrum analysis. The idea was to choose the component with the highest breathing information, computing the power spectral density estimate (PSD) between 0.5-2 Hz for each component and selecting the component with: (1) maximum PSD peak ("Peak" method) or (2) maximum area under the PSD ("Area" method). To assess the goodness of these two methods in predicting the best quaternion component, the ideal component ("Ideal") was determined a posteriori, case by case, based on minimum breathing frequency estimation error (see Section 2.5).
Since more than one quaternion component is supposed to convey breathing information, the possibility to maximize this information fusing the 4 components of the quaternion by means of PCA was investigated. PCA is a mathematical procedure that transforms an original set of correlated variables into a (smaller) number of uncorrelated variables by determining a set of orthogonal vectors called principal components, which are defined by a linear combination of the original variables [30,31]. To do this, the directions in the data with the most variation, i.e., the eigenvectors corresponding to the largest eigenvalues of the covariance matrix, are computed and the data are projected onto these directions. To compute the eigenvectors, data were arranged into a twodimensional matrix X(m × n), where m was the number of observations of the time series and n the number of variables (quaternion components). Then, the univariate means were subtracted from the n columns, to center the data. Singular Value Decomposition (SVD) was used to compute the eigenvectors (V = [v1, v2, v3, v4]) and corresponding eigenvalues (λ1, λ2, λ3, λ4). Original data were  (2) and (3), that are composed by 4 components each, and provides as output 2 single-component signals (1 for the abdomen and 1 for the thorax) representing chest-wall respiratory-related orientation change signals. These two signals represent the input of the power spectrum block and of the processing block. To reduce dimension from 4 components to 1, two possibilities were investigated as shown in Figure 2: To select the best component among the 4 components representing the orientation quaternion, two different methods were proposed, both based on spectrum analysis. The idea was to choose the component with the highest breathing information, computing the power spectral density estimate (PSD) between 0.5-2 Hz for each component and selecting the component with: (1) maximum PSD peak ("Peak" method) or (2) maximum area under the PSD ("Area" method). To assess the goodness of these two methods in predicting the best quaternion component, the ideal component ("Ideal") was determined a posteriori, case by case, based on minimum breathing frequency estimation error (see Section 2.5).

Dimension-reduction block takes the quaternions obtained from Equations
Since more than one quaternion component is supposed to convey breathing information, the possibility to maximize this information fusing the 4 components of the quaternion by means of PCA was investigated. PCA is a mathematical procedure that transforms an original set of correlated variables into a (smaller) number of uncorrelated variables by determining a set of orthogonal vectors called principal components, which are defined by a linear combination of the original variables [30,31]. To do this, the directions in the data with the most variation, i.e., the eigenvectors corresponding to the largest eigenvalues of the covariance matrix, are computed and the data are projected onto these directions. To compute the eigenvectors, data were arranged into a two-dimensional matrix X(m × n), where m was the number of observations of the time series and n the number of variables (quaternion components). Then, the univariate means were subtracted from the n columns, to center the data. Singular Value Decomposition (SVD) was used to compute the eigenvectors (V = [v 1 , v 2 , v 3 , v 4 ]) and corresponding eigenvalues (λ 1 , λ 2 , λ 3 , λ 4 ). Original data were finally projected in the new coordinate system (Y = XV) and the first principal component, accounting for the largest possible variance, was selected and passed to other blocks [30,31]. finally projected in the new coordinate system (Y = XV) and the first principal component, accounting for the largest possible variance, was selected and passed to other blocks [30,31]. Spectrum Analysis block include a set of steps needed to optimize the subsequent processing phase. The two signals representing chest-wall (abdominal and thoracic) respiratory-related orientation obtained downstream of the dimension-reduction block underwent the following steps ( Figure 1): A low-frequency threshold (fLOW) was determined based on a first estimate of the breathing frequency (fB). The rough estimate of fB was done by identifying maxima points of the signal and computing the fB, breath by breath, as reciprocal of the temporal distance between consecutive maxima points. Then, the mean (fB_Rough) and the standard deviation (fB_Rough_SD) of the fB over the entire trial were computed. To facilitate maxima points identification, signals were at first bandbass filtered using a first-order infinite impulse response (IIR) Butterworth filter [0.05H-2 Hz] and smoothed with a third-order Savitzky-Golay [32] finite impulse response (FIR) filter (fixed window length = 31 samples). Low thresholds fLOWAb and fLOWTh were determined for the abdominal and thoracic signals respectively as difference fB_Rough − fB_Rough_SD. Then the minimum value between fLOWAb and fLOWTh was chosen as final low-frequency threshold, named fLOW, and it was used in the next step. (ii).
PSD estimate (Welch's method, Hamming window size: 300 samples, overlapping: 50 samples) was computed and the spectrum frequency corresponding to the breathing rate was identified, both for the thorax (fpeak_T) and the abdomen (fpeak_A), by looking for the local peak of the PSD within the window [fLOW ÷ 2 Hz]. The use of a low threshold, based on a rough estimate of the breathing frequency, supports the selection of the PSD peak linked to breathing rate and avoid selecting wrong peaks, often related to low-frequency oscillation artifacts. Spectrum Analysis block include a set of steps needed to optimize the subsequent processing phase. The two signals representing chest-wall (abdominal and thoracic) respiratory-related orientation obtained downstream of the dimension-reduction block underwent the following steps ( Figure 1): (i). A low-frequency threshold (f LOW ) was determined based on a first estimate of the breathing frequency (f B ). The rough estimate of f B was done by identifying maxima points of the signal and computing the f B , breath by breath, as reciprocal of the temporal distance between consecutive maxima points. Then, the mean (f B_Rough ) and the standard deviation (f B_Rough_SD ) of the f B over the entire trial were computed. To facilitate maxima points identification, signals were at first band-bass filtered using a first-order infinite impulse response (IIR) Butterworth filter [0.05 Hz-2 Hz] and smoothed with a third-order Savitzky-Golay [32] finite impulse response (FIR) filter (fixed window length = 31 samples). Low thresholds f LOW Ab and f LOW Th were determined for the abdominal and thoracic signals respectively as difference f B_Rough − f B_Rough_SD . Then the minimum value between f LOW Ab and f LOW Th was chosen as final low-frequency threshold, named f LOW , and it was used in the next step. (ii). PSD estimate (Welch's method, Hamming window size: 300 samples, overlapping: 50 samples) was computed and the spectrum frequency corresponding to the breathing rate was identified, both for the thorax (f peak_T ) and the abdomen (f peak_A ), by looking for the local peak of the PSD within the window [f LOW ÷ 2 Hz]. The use of a low threshold, based on a rough estimate of the breathing frequency, supports the selection of the PSD peak linked to breathing rate and avoid selecting wrong peaks, often related to low-frequency oscillation artifacts. (iii). The breathing frequency derived by the spectrum was used to set an adaptive band-pass filter, as proposed in a previous study [7], centered on f peak frequency. For the abdomen, upper (f U ) and lower (f L ) cut-off frequency points for the band-pass filter were defined, by applying Equations (4) and (5) respectively [7]: For the thorax, Equations (6) and (7) were applied: Moreover, based on f peak , a set of parameters was selected to optimize subsequent smoothing and minima/maxima detection phases of the processing block.
Processing block includes all the steps needed to extract breathing frequency and temporal parameters from the signals obtained downstream of the dimension-reduction block. Chest-wall respiratory-related orientation change signals (abdominal and thoracic) underwent the following steps: (i). Adaptive band-pass filter. The signals were band-pass filtered (first-order IIR Butterworth filter), with f U and f L cut-off frequency points determined within the spectrum analysis block. (ii). Smoothing. Filtered signals were furtherly smoothed (third-order Savitzky-Golay FIR filter) to simplify subsequent identification of maxima and minima points. The level of smoothing (window length) was automatically selected based on f peak , i.e., increasing window length for decreasing f peak . Relation between optimal window length values and f peak values has been determined empirically. (iii). Minima and maxima points detection. A set of optimized parameters (i.e., minimum peak distance (MPD) and minimum prominence threshold (MPT)) was automatically selected based on f peak to optimize recognition of minima and maxima points of the smoothed signals. Optimal MPD and MPT values depending on f peak were experimentally determined.
(iv). Breathing frequency extraction. Breath by breath, inspiratory time (T I ) was computed as the temporal distance between a minimum point (m i ) and the consecutive maximum point (M i ); Expiratory time (T E ) was computed as the temporal distance between the maximum point (M i ) and the consecutive minimum point (m i + 1); total time (T TOT ) was computed as T TOT = T I + T E [s], duty cycle (DC) was computed as T I T TOT × 100 [%] and breathing frequency was computed as 60 (T TOT ) [breaths/minute]. A mean value for each of the above-mentioned parameter was computed for each trial (~3 min).

Experimental Setup
To evaluate the capability of the device and of the proposed methods to correctly estimate breathing frequency (and temporal parameters) 8 healthy volunteers (4 males, 4 females) were enrolled. All subjects gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol (Project identification code n • 534) was approved by the Ethics Committee of Scientific Institute IRCCS Medea (date of approval: 25 January 2018). Chest-wall movements during breathing, in seated and supine position, were measured using the proposed device and optoelectronic plethysmography (OEP) simultaneously. OEP [21] is a technique based on a similar functioning principle of the proposed device; in fact, it allows assessment of ventilatory and breathing pattern by measuring chest-wall movements related to breathing, by using motion capture principles. The system is composed of eight infrared video cameras working at a sampling rate of 60 Hz. It can compute the 3D coordinates of retro-reflective markers positioned on the chest wall in specific anatomic points. From the three-dimensional coordinates of Sensors 2019, 19, 88 8 of 24 the markers, it is possible to obtain the volume enclosed by the chest-wall surface, by applying the Gauss's theorem. The chest wall is modelled by a bicompartmental model, composed of rib cage and abdomen, and thus it is possible to investigate the contribution of both the compartments. This is an advantage for the validation of the proposed device, in fact, using OEP as reference method it is possible to compare the data recorded with the thoracic and abdominal units of the device with those obtained by using OEP for the thoracic and abdominal compartments, respectively. OEP has been widely validated against spirometer, in healthy subjects, in different conditions and positions, also during submaximal and maximal exercise on cycle ergometer, obtaining discrepancies in tidal volume measurements always <5% [22][23][24]33,34].
The subjects were prepared, placing the reflective markers according to the 89-marker protocol (previously described in [24,35]) used for seated position and the 52-marker protocol (previously described in [36,37]) used for supine position or, more generally, when a back support is present. Then peripheral IMU-units were placed on the thorax and on the abdomen, while reference IMU-unit was placed on the coccyx in seated position, and on the bed in supine position ( Figure 3). dimensional coordinates of the markers, it is possible to obtain the volume enclosed by the chest-wall surface, by applying the Gauss's theorem. The chest wall is modelled by a bicompartmental model, composed of rib cage and abdomen, and thus it is possible to investigate the contribution of both the compartments. This is an advantage for the validation of the proposed device, in fact, using OEP as reference method it is possible to compare the data recorded with the thoracic and abdominal units of the device with those obtained by using OEP for the thoracic and abdominal compartments, respectively. OEP has been widely validated against spirometer, in healthy subjects, in different conditions and positions, also during submaximal and maximal exercise on cycle ergometer, obtaining discrepancies in tidal volume measurements always <5% [22][23][24]33,34]. The subjects were prepared, placing the reflective markers according to the 89-marker protocol (previously described in [24,35]) used for seated position and the 52-marker protocol (previously described in [36,37]) used for supine position or, more generally, when a back support is present. Then peripheral IMU-units were placed on the thorax and on the abdomen, while reference IMUunit was placed on the coccyx in seated position, and on the bed in supine position ( Figure 3). Subjects were then asked to seat or lie on a bed and were invited to perform a slow vital capacity maneuver (SVC) and then to start breathing with the following patterns: (I) quiet breathing (QB), (II) increasing fB but same tidal volume of QB (↑fB, VT=), (III) increasing fB and reducing tidal volume (↑fB, VT↓), (IV) decreasing fB with the same tidal volume of QB (↓fB, VT=), (V) decreasing fB increasing tidal volume of QB (↓fB, VT↑). QB trial was repeated two times, thus, each subject performed 6 trials of the duration of 3 minutes each. The SVC maneuver was used to align OEP signal and device signals during data analysis, since it is generally recognizable with respect to QB. In fact, SVC requires a maximal inspiration followed by a complete expiration without forced or rapid effort.
The subjects were asked to maintain the same breathing pattern (namely, QB, ↑fB, VT=, ↑fB, VT↓, ↓fB, VT=, ↓fB, VT↑) until the end of the trial; in case of fatigue they were asked to perform a second SVC before returning to QB. This procedure was repeated in seated position and in supine position. Subjects were then asked to seat or lie on a bed and were invited to perform a slow vital capacity maneuver (SVC) and then to start breathing with the following patterns: (I) quiet breathing (QB), (II) increasing f B but same tidal volume of QB (↑f B , V T =), (III) increasing f B and reducing tidal volume (↑f B , V T ↓), (IV) decreasing f B with the same tidal volume of QB (↓f B , V T =), (V) decreasing f B increasing tidal volume of QB (↓f B , V T ↑). QB trial was repeated two times, thus, each subject performed 6 trials of the duration of 3 min each. The SVC maneuver was used to align OEP signal and device signals during data analysis, since it is generally recognizable with respect to QB. In fact, SVC requires a maximal inspiration followed by a complete expiration without forced or rapid effort. The subjects were asked to maintain the same breathing pattern (namely, QB, ↑f B , V T =, ↑f B , V T ↓, ↓f B , V T =, ↓f B , V T ↑) until the end of the trial; in case of fatigue they were asked to perform a second SVC before returning to QB. This procedure was repeated in seated position and in supine position.

Statistical Analysis
For each trial, mean values of f B , T I , T E and DC were extracted from the best quaternion components identified online by using "Area" and "Peak" methods, and from the signal obtained with the PCA-based fusion method, both for the thoracic and abdominal tracings. Moreover, to evaluate the performance of the selection methods ("Area" and "Peak") and their ability to select the best component, the same parameters were obtained for all the quaternion components (q 0 , q 1 , q 2 , q 3 ) and compared with those obtained by OEP, on the abdominal and thoracic compartment, respectively. The "Ideal" quaternion component was identified a posteriori, trial by trial, as the one providing the minimum estimation error of the breathing frequency. Obviously, the "Ideal" component cannot be identified during online analysis, or when a reference method is not present. Thus, for each trial, 5 sets of parameters were available: Among the entire set of trials, those with f B_OEP < 6 breaths/minute or f B_OEP > 60 breaths/minute were discarded. Then, the absolute (Equation (8)) and relative (Equation (9)) errors of estimation in static conditions (supine and seated position) were computed for each parameter: Relative For all the dimension-reduction methods ("Area", "Peak", "PCA"), mean and standard deviation (SD) were computed for E and E% considering all the subjects and all the trials, for the supine and seated position and compared with those obtained considering the "Ideal" component. The error obtained with the "Ideal" component, identified a posteriori, is thus the minimum error obtainable using a single quaternion component, and represents the performance that the other methods ("Area", "Peak", and PCA-fusion) should achieve or beat. For E% obtained in f B and DC estimation, non-parametric alternative to the one-way Analysis of variance (ANOVA) with repeated measures (Friedman test) was performed to assess if significant differences between methods ("Area", "Peak", "PCA") and "Ideal" component occurred, "Ideal"); post-hoc analysis was done performing Wilcoxon signed-rank tests on the different combinations of related methods, applying the correction for multiple comparisons using false discovery rate (FDR) method [38,39].
For f B , T I , T E , linear regression analysis and correlation analysis (Pearson's product-moment correlation r P , or Spearman's rank-order correlation r S , if data were not normally distributed) were performed between measurements obtained with the device and measurements obtained with the OEP, for the supine and seated position, respectively.
To assess the agreement between measurements obtained with the device and with the OEP, Bland-Altman analysis was performed plotting the difference of the two paired measurements (device-OEP) against the mean of the two measurements [40][41][42]. Mean of the differences (d) and limits of agreement (LOA: from d − (1.9 × SD) to d + (1.9 × SD)) were calculated. The presence of heteroscedasticity was always examined to assess the presence of proportional biases and/or the correlation between differences and mean values. As proposed by Brehm et al. [43], to determine if data were heteroscedastic a visual inspection of Bland-Altman plots was performed at first. If the errors (y-axes: absolute differences) increased with increasing measured values (x-axes: mean), the data were suspected of being heteroscedastic. Then Kendall's tau (τ) correlation between the absolute differences and the corresponding means was computed to assess the degree of heteroscedasticity. Data were denoted heteroscedastic when a positive, significant correlation (τ > 0.1 and p-value < 0.05) was found, for other cases data were considered homoscedastic [43].
When heteroscedasticity was present the "classical" 95% confidence and tolerance limits cannot be constructed, thus the approach based on the construction of V-shaped limits was applied: the regression line (ordinary least squares (OLS) best fit) was constructed for differences on mean values and the V-shaped confidence limits (upper confidence limit: UCL, lower confidence limit: LCL) were constructed modelling the variability in the SD of the differences directly as a function of the level of the measurement, using a method based on absolute residuals from a fitted regression line [44,45]. Table 1 presents the mean and SD of breathing rate for each breathing pattern (QB 1 e QB 2 , ↑f B , V T =, ↑f B , V T ↓, ↓f B , VT=, ↓f B , V T ↑) estimated with OEP and device, using "PCA", "Area", and "Peak" methods and the "Ideal" component, for all subjects, in supine and seated position. Sample size (n) of each condition is reported in Table 1 for the breathing pattern ↑f B , V T ↓, just one thoracic tracing was available for seated position (n = 1). It can be noticed that each subject demonstrated a different breathing frequency for each breathing pattern and SD in the forced breathing patterns is higher than those obtained for QB, meaning that subjects interpreted the required speed differently.     Across subject mean ± SD of the breathing frequency (f B , [breaths/minute]) measurements with OEP and the device, using best component-selection methods ("Area" and "Peak"), PCA-fusion method and "Ideal" component, for the requested patterns. Data are reported for the supine and seated positions, subdivided in abdominal and thoracic contributions.

Accuracy Errors
Relative errors of estimation in supine and seated position computed for best component-selection methods ("Peak", "Area"), for PCA-fusion method and for the "Ideal" quaternion component for f B and DC are presented in Figure 4. For what concerns f B estimation in supine position, relative errors obtained using PCA were similar or even better than those provided by the "Ideal" component; on the contrary, both component-selection methods, namely "Area" and "Peak", provided errors higher than 10%, both for the abdominal and the thoracic compartments. Errors obtained with PCA resulted significantly lower than those obtained with the "Area" method, both for the abdominal (Wilcoxon post-hoc test FDR-adjusted, p = 0.038) and thoracic compartment (Wilcoxon post-hoc test FDR-adjusted, p = 0.015); also, PCA was significantly better than "Peak" method considering abdominal compartment

Accuracy Errors
Relative errors of estimation in supine and seated position computed for best componentselection methods ("Peak", "Area"), for PCA-fusion method and for the "Ideal" quaternion component for fB and DC are presented in Figure 4. For what concerns fB estimation in supine position, relative errors obtained using PCA were similar or even better than those provided by the "Ideal" component; on the contrary, both component-selection methods, namely "Area" and "Peak", provided errors higher than 10%, both for the abdominal and the thoracic compartments. Errors obtained with PCA resulted significantly lower than those obtained with the "Area" method, both for the abdominal (Wilcoxon post-hoc test FDR-adjusted, p = 0.038) and thoracic compartment (Wilcoxon post-hoc test FDR-adjusted, p = 0.015); also, PCA was significantly better than "Peak" method considering abdominal compartment (Wilcoxon post-hoc test FDR-adjusted, p = 0.038).  In seated position, fB estimation errors obtained with component-selection methods were lower on average than those obtained in supine position, while PCA performances declined. This led to a sort of equalization effect, confirmed also by the statistical analysis: significant differences remained only for comparisons "Ideal" vs. "Area" method (Wilcoxon post-hoc test FDR-adjusted, AB: p = 0.102, In seated position, f B estimation errors obtained with component-selection methods were lower on average than those obtained in supine position, while PCA performances declined. This led to a sort of equalization effect, confirmed also by the statistical analysis: significant differences remained only for comparisons "Ideal" vs. Absolute estimation errors of f B , T I and T E obtained with the device using different methods (Area, Peak, Ideal, PCA) relative to OEP are reported in Table 2. Table 2. Absolute errors of breathing frequency (E_f B ), Inspiratory time (E_T I ) and expiratory time (E_T E ) obtained for the device with respect to OEP, using best component-selection methods ("Area" and "Peak"), PCA-fusion method and "Ideal" component. Data are reported as mean ± SD, in supine and seated position for thoracic (TH) and abdominal (AB) compartments.

Linear Regression and Correlation Analysis
Scatter plots showing the relationship between measurements obtained with the OEP and with the device, using "Area", "Peak", "Ideal" components, and PCA-fusion respectively are presented for f B (Figure 5), T I ( Figure 6) and T E (Figure 7). For each scatter plot the regression line is computed, both for thorax and abdomen, and the relative equations are reported.
Correlation coefficients for the comparisons Device vs. OEP are reported in Table 3. Regarding the main parameter, f B , results obtained with correlation analysis confirmed what emerged from estimation error analysis: in supine position, PCA exhibited the best performances in terms of correlation with OEP measurements both in terms of regression line and correlation coefficient. In seated position, "Ideal" component was the one with the highest correlation with OEP measurements, followed by PCA.       Regarding TE measurements obtained with the IMU-device, three dimension-reduction methods were considered. Area, Peak and PCA-fusion. The performance obtained by using these three methods is benchmarked against that obtained with the Ideal quaternion component determined a posteriori based on the minimum estimation error. The regression line between measurements done by OEP and the proposed device is plotted, and the relative equation presented, both for the thorax and the abdomen. Table 3 Correlation outcomes across subjects and breathing patterns. Coefficient of correlation (r) between measurements obtained using Device vs. OEP are reported for fB. TI. and TE using best component-selection methods ("Area" and "Peak"), PCA-fusion method and "Ideal" component, in supine (Thorax: n = 37. Abdomen: n = 37) and seated (Thorax: n = 35. Abdomen: n = 39) position. With reference to TI estimation in supine position, "Ideal" component provided the best performances, followed by PCA method; "Peak" and "Area" methods provided comparable, poor performances. In seated position, measurements of TI provided by component-selection methods Regarding T E measurements obtained with the IMU-device, three dimension-reduction methods were considered. Area, Peak and PCA-fusion. The performance obtained by using these three methods is benchmarked against that obtained with the Ideal quaternion component determined a posteriori based on the minimum estimation error. The regression line between measurements done by OEP and the proposed device is plotted, and the relative equation presented, both for the thorax and the abdomen. Table 3. Correlation outcomes across subjects and breathing patterns. Coefficient of correlation (r) between measurements obtained using Device vs. OEP are reported for f B . T I . and T E using best component-selection methods ("Area" and "Peak"), PCA-fusion method and "Ideal" component, in supine (Thorax: n = 37. Abdomen: n = 37) and seated (Thorax: n = 35. Abdomen: n = 39) position.

Supine Seated
Thorax With reference to T I estimation in supine position, "Ideal" component provided the best performances, followed by PCA method; "Peak" and "Area" methods provided comparable, poor performances. In seated position, measurements of T I provided by component-selection methods were on average more correlated with OEP measurements than measurements obtained using PCA-fusion method. The "Ideal" component presented the best results, followed by "Area" and "Peak" methods; PCA provided the worst performance considering the abdominal compartment, while correlation between measurements obtained with the thoracic unit and OEP measurements was good.
Estimation of T E was on average more problematic. In terms of regression lines in fact, slope values were far from the unity for all the considered methods, highlighting a proportional error leading to an overestimation for low values of expiratory time and an underestimation at high expiratory times, as shown in Figure 7. For what concerns supine position, correlation coefficients were good both for "Ideal" component and PCA-fusion method; on the contrary, correlation coefficients were low both for "Area" and "Peak" methods. Also, in seated position correlation coefficients provided by "Ideal" and PCA-fusion method were higher than those provided by "Area" and "Peak" methods, especially with respect to the thoracic compartment.

Bland-Altman Analysis
Bland-Altman plots representing the agreement between measurements obtained with the OEP and with the device, using "Area", "Peak", "Ideal" components, and PCA-fusion respectively are presented for f B (Figure 8), T I ( Figure 9) and T E (Figure 10). In Bland-Altman plots, the difference of the two paired measurements (device-OEP) is plotted against the mean of the two measurements (device+OEP)⁄2. Results of agreement analysis, including evaluation of heteroscedasticity (Kendall's τ correlation and relative p-value) are reported in Table 4. As shown there, for homoscedastic data, the mean of the differences representing the fixed bias, and LOAs were computed. On the other hand, for heteroscedastic data, OLS line of best fit representing the proportional bias and upper and lower 95% V-shape confidence limits (UCL and LCL) are reported. were on average more correlated with OEP measurements than measurements obtained using PCAfusion method. The "Ideal" component presented the best results, followed by "Area" and "Peak" methods; PCA provided the worst performance considering the abdominal compartment, while correlation between measurements obtained with the thoracic unit and OEP measurements was good. Estimation of TE was on average more problematic. In terms of regression lines in fact, slope values were far from the unity for all the considered methods, highlighting a proportional error leading to an overestimation for low values of expiratory time and an underestimation at high expiratory times, as shown in Figure 7. For what concerns supine position, correlation coefficients were good both for "Ideal" component and PCA-fusion method; on the contrary, correlation coefficients were low both for "Area" and "Peak" methods. Also, in seated position correlation coefficients provided by "Ideal" and PCA-fusion method were higher than those provided by "Area" and "Peak" methods, especially with respect to the thoracic compartment.

Bland-Altman Analysis.
Bland-Altman plots representing the agreement between measurements obtained with the OEP and with the device, using "Area", "Peak", "Ideal" components, and PCA-fusion respectively are presented for fB (Figure 8), TI ( Figure 9) and TE ( Figure 10). In Bland-Altman plots, the difference of the two paired measurements (device-OEP) is plotted against the mean of the two measurements (device+OEP)⁄2. Results of agreement analysis, including evaluation of heteroscedasticity (Kendall's τ correlation and relative p-value) are reported in Table 4. As shown there, for homoscedastic data, the mean of the differences representing the fixed bias, and LOAs were computed. On the other hand, for heteroscedastic data, OLS line of best fit representing the proportional bias and upper and lower 95% V-shape confidence limits (UCL and LCL) are reported. Figure 8. Agreement analysis between OEP and the IMU-based device for breathing frequency (fB, expressed as breaths/minuteute) measurements, in supine (top panels) and seated (bottom panels) position. In each Bland-Altman plot the differences between measurements of fB obtained by using the IMU-based device and by using OEP are plotted against the mean of the two measurements. For homoscedastic data, the mean of the differences (bias: -) and limits of agreement (black dotted line) from mean − 1.96 s to mean + 1.96 s are represented by lines parallel to the X axis. For heteroscedastic data, the proportional bias (-) is represented by the ordinary least squares (OLS) line of best fit for the difference on mean values; V-shaped upper and lower 95% confidence limits (---) are calculated according to Bland [44]. . Agreement analysis between OEP and the IMU-based device for breathing frequency (f B , expressed as breaths/minuteute) measurements, in supine (top panels) and seated (bottom panels) position. In each Bland-Altman plot the differences between measurements of f B obtained by using the IMU-based device and by using OEP are plotted against the mean of the two measurements. For homoscedastic data, the mean of the differences (bias: -) and limits of agreement (black dotted line) from mean − 1.96 s to mean + 1.96 s are represented by lines parallel to the X axis. For heteroscedastic data, the proportional bias (-) is represented by the ordinary least squares (OLS) line of best fit for the difference on mean values; V-shaped upper and lower 95% confidence limits (---) are calculated according to Bland [44].      With respect to the main parameter (f B ), agreement between OEP and the device is very strong when the "Ideal" component or the PCA-fusion are used, both in supine and seated position. In relation to time estimation, the agreement decreases for all the methods considered. In particular, for what concerns inspiratory times, a significant relationship between errors and mean value emerged, with a general increase of the difference (device-OEP) at higher time values (overestimation of the device), both in supine and seated position. Also, for expiratory times absolute errors increased with increasing time values, but in this case the device underestimated at high time values (negative slope of the OLS).

Quaternion Component Selection
Considering the quaternion components selected by the "Area" and "Peak" methods as best component or identified as "Ideal" component, a clear rule did not emerge. In fact, there was not a quaternion component that was selected as best component with a considerable frequency. Relative frequencies of quaternion component selection with the different methods are presented in Figure 11. It is interesting to notice that quaternion component q 0 was never selected by "Area" and "Peak" methods, while the "Ideal" component was q 0 in 14.86 % of cases (n = 74) in supine position and 6.76% of cases (n = 74) in seated position. In seated position, the component q 1 was selected more frequently as best component both by using "Area" (44.59) and "Peak" (39.19%) methods. In contrast, in supine position, the components q 2 ("Area" 51.35%, "Peak": 41.89) and q 3 ("Area": 39.19% and "Peak": 40.54%) were selected more frequently.
In regard to the ability of the two component-selection methods ("Area" and "Peak") to identify the "Ideal" component, i.e., the component providing the minimum f B estimation error, in supine position, the "Area" method was able to identify the "Ideal" component in 45.94% of the cases (relative frequency for the event "the component selected by "Area" method and the "Ideal" component corresponded"), the "Peak" method identified the "Ideal" component in 52.70% of the cases, while for the 45.94% of the cases neither the "Area" method nor the "Peak" method were able to identify the "Ideal" component. In 44.59% of cases, "Area" method and "Peak" method selected the same quaternion component, that was also identified as "Ideal" component.
what concerns inspiratory times, a significant relationship between errors and mean value emerged, with a general increase of the difference (device-OEP) at higher time values (overestimation of the device), both in supine and seated position. Also, for expiratory times absolute errors increased with increasing time values, but in this case the device underestimated at high time values (negative slope of the OLS).

Quaternion Component Selection
Considering the quaternion components selected by the "Area" and "Peak" methods as best component or identified as "Ideal" component, a clear rule did not emerge. In fact, there was not a quaternion component that was selected as best component with a considerable frequency. Relative frequencies of quaternion component selection with the different methods are presented in Figure  11. It is interesting to notice that quaternion component q0 was never selected by "Area" and "Peak" methods, while the "Ideal" component was q0 in 14.86 % of cases (n = 74) in supine position and 6.76% of cases (n = 74) in seated position. In seated position, the component q1 was selected more frequently as best component both by using "Area" (44.59) and "Peak" (39.19%) methods. In contrast, in supine position, the components q2 ("Area" 51.35%, "Peak": 41.89) and q3 ("Area": 39.19% and "Peak": 40.54%) were selected more frequently.
In regard to the ability of the two component-selection methods ("Area" and "Peak") to identify the "Ideal" component, i.e., the component providing the minimum fB estimation error, in supine position, the "Area" method was able to identify the "Ideal" component in 45.94% of the cases (relative frequency for the event "the component selected by "Area" method and the "Ideal" component corresponded"), the "Peak" method identified the "Ideal" component in 52.70% of the cases, while for the 45.94% of the cases neither the "Area" method nor the "Peak" method were able to identify the "Ideal" component. In 44.59% of cases, "Area" method and "Peak" method selected the same quaternion component, that was also identified as "Ideal" component.

Discussion
In this study, we presented an automatic, position-independent processing algorithm to derive breathing signal, and subsequently breathing temporal parameters, from chest-wall orientation changes acquired using an IMU-based device previously developed by our group [19], composed of three sensor units. Even if the modular configuration of the device was designed to reduce non-breathing movements, the aim of this work was neither to demonstrate the effectiveness of this approach nor to support the presence of the reference unit. On the contrary the focus is on the analysis algorithm, which uses quaternion form to represent orientation, avoiding singularity problem that affects Euler representation; thus, thoracic and abdominal orientation change signals are 4-dimensional entities. The proposed algorithm includes therefore a dimension-reduction block to obtain a 1-dimension signal representing chest-wall orientation changes due to breathing activity.
Another aim of this study was, therefore, to compare three different dimension-reduction methods. The first two methods ("Peak" and "Area") are based on the selection of the quaternion component with the highest breathing information. The third method is based on the fusion of the 4 components of the quaternion using PCA.
The PCA-fusion method performed better than the best component-selection methods ("Peak" and "Area") as regards breathing frequency estimation, both in supine and seated position. In supine position, it provided better results than "Ideal" component, while in seated position it provided closer performances to those obtained by using the "Ideal" component with respect to "Area" and "Peak" methods.
About estimation of other temporal parameters (T I , T E and duty cycle), "Ideal" component provided the best results, while PCA-fusion method gave results comparable to the best component-selection methods. Thus, a quaternion component providing the best performance exists, the problem lies in its a priori identification. In fact, both "Area" and "Peak" methods failed to identify it on the basis of spectral analysis (in 45.94% of the cases neither the "Area" method nor the "Peak" method were able to identify the "Ideal" component), and no quaternion component emerged as the most selected as "Ideal" component (supine q 0 :14.86%, q 1 : 22.97%, q 2 : 32.43%, q 3 : 29.73%; seated q 0 : 6.76%, q 1 : 33.78, q 2 : 22.97%, q 3 : 36.49%).
Geometrical or morphological considerations to determine which quaternion component is more involved in breathing movement are problematic when considering quaternions, and are position-and IMU-placement dependent, thus not suitable in dynamic conditions. On the contrary, PCA-fusion method represents an interesting solution to this problem because it fuses the information of the four quaternion components regardless the position of the subject or the IMUs placement, avoiding the necessity to select a best component/axes, as reported in previous studies [7,11,13].
In this study, we found that PCA-fusion method provided the best f B estimation performance in terms of mean absolute errors (<2 breaths/minute), correlation (r > 0.963) and agreement (see Table 4) with the reference method. Comparing our results in terms of accuracy errors with those obtained by previous studies is difficult because in most cases only relative errors were reported, but these errors depend on the breathing frequency adopted. Liu et al. [14] reported a mean absolute error of 15.45 breaths/minute (thus about 7 times higher than the error obtained in this study) during quiet sitting. Bates et al. [13] obtained an RMS error of 0.38 breaths/minute and a peak error of 3 breaths/minute in a postoperative patient during sleep. Considering comparable conditions (abdominal compartment in supine position) we obtained an RMS error of 1.51 breaths/minute, using PCA-fusion method, but in our study different, forced, breathing patterns were included, leading to higher mean error.
Regarding correlation between f B measurements obtained with the proposed method and OEP, our results are comparable to those obtained by Bates et al. [13] that reported a correlation coefficient equal to 0.985 between measurements of f B obtained with the accelerometer and nasal cannula. Mann et al. [11] obtained a correlation coefficient of 0.97 between measurements of f B obtained with a tri-axial accelerometer and with a system based on oxygen consumption measurement (Oxycon Mobile). In both cases [11,13], breath-by-breath analysis was not possible. Liu et al. [14] reported correlation coefficients lower than 0.6 between f B computed with a 3-axis accelerometer and with the reference (Airflow CO 2 analysis).
There are few studies in the literature that performs Bland-Altman analysis to assess the agreement between breathing frequency measurements obtained by using inertial sensors and other validated methods. In the study from Morillo et al. [8] agreement analysis using Bland Altman plots was done against PSG thermistor. The authors reported a mean difference (fixed bias) of 0.02 (SD = 1.09) breaths/minute and LOAs from −3.05 to +2.11 breaths/minute in the range~12 ÷ 35 breaths/minute. In that case, the use of a single-axis accelerometer, prevent the use of that method during postural changes. Dehkrodi et al. [9] used a tri-axial accelerometer placed on the suprasternal notch extending the validation presented in [8] to different sleeping positions and breathing conditions (Deep: 13.5 ± 4.3, Normal: 16.5 ± 5.2 and Shallow: 39.7 ± 30.3 breath/minute). They reported a mean difference (fixed bias) of 0.042 breaths/minute and LOAs from -0.65 to 0.74 breaths/minute. Lapi et al. [17] performed agreement analysis with Bland Altman plots for measurements of breathing frequency (range 12 ÷ 26 breaths/minute) obtained with the accelerometer and with the standard method (counting breaths by visual inspection), in supine position. They reported a mean difference (fixed bias) of 0.33 breaths/minute and LOA from −1.92 to 2.60 with 3.2% of data outside that range. For all the above-mentioned studies heteroscedasticity of data was not considered or reported making it difficult compare them directly with our results. In fact, taking into account heteroscedasticity of data, for f B in supine position, we built Bland-Altman plot with proportional bias (OLS: y = 0.008x + 0.130, thus going from 0.18 (at x = 6 breaths/minute) to 0.61 (at x = 60 breaths/minute) breaths/minute) and V-shaped limits (LCL: y = −0.038x − 2.039 thus the lower limit goes from −2.26 (at x = 6 breaths/minute) to −4.32 (at x = 60 breaths/minute) breath/minute; UCL: y = 0.054x + 2.299; thus the upper limit goes form 2.62 (at x = 6 breaths/minute) to 5.539 (at x = 60 breaths/minute) breaths/minute) ). Thus, considering comparable breathing frequency ranges our results are closer to those obtained by Morillo et al. [8] and Lapi et al. [17]. On the contrary, Dehkrodi et al. [9] obtained better results; unfortunately, the steps to obtain the acceleration derived respiratory (ADR) signal are not described in detail.
For the best of our knowledge, this is the first study that provides a detailed analysis of respiratory timing measurements obtained by using inertial sensor systems, validating them against an established method.

Conclusions
PCA-fusion method provided overall best performances with respect to selecting the best quaternion component identified based on spectrum analysis. In supine position results obtained fusing the 4 quaternion components were even better than those obtained with the "Ideal" component, identified a posteriori considering the minimum breathing frequency estimation error. Performance in seated position were worse than those obtained in supine position, probably because subjects were seated without the back support and some oscillations of the trunk were more likely to occur. This could particularly affect PCA-based method where the first principal component selected for further analysis is the one with the largest variance, and thus more subject to larger body motions. This must be taken into account in dynamic conditions; signal baseline removal prior to PCA-fusion should be considered in this case.
The analysis algorithm proposed in this work, applying PCA-fusion as dimension-reduction method, can be used to analyze further data. In fact, an extended validation of the proposed device and method is needed also in dynamic conditions, during daily activities, considering not only healthy subjects but also patients that could particularly take advantage of this system (e.g., COPD, neuromuscular patients, sleep apnea, etc.). This would also allow study of asynchronies of thoraco-abdominal compartments taking advantage of the modular configuration of the device.
Another key step will be the adaptation of our method, currently implemented as an offline analysis, to online monitoring, moving the computation process aboard the smartphone. This enhancement could allow immediate computation of an average breathing frequency over a certain period (e.g., 60 s) directly aboard the smartphone, fostering the use of the device in other applications such as sport and fitness, exercise testing, breathing training to use different respiratory muscles, rehabilitation protocols and treatment evaluation where respiratory assessment could be of great interest.

Patents
The present work is partially described in the International Patent application n • PCT/IB2018/ 054956, priority date 11 July 2017, title "A wearable device for the continuous monitoring of the respiratory rate". Inventors: Ambra Cesareo, Andrea Aliverti, Assignee: Politecnico di Milano.
Author Contributions: A.C. and A.A. contributed IMU circuit design and fabrication; A.C. conceived and designed the experiments; A.C. and Y.P. performed the experiments and analyzed the data; A.C. and E.B. prepared the original draft; E.B. and A.A. reviewed the paper; E.B. and A.A. supervised the project and contributed to funding acquisition.
Funding: This research received no external funding.