A pneumatic Bionic Voice prosthesis—Pre-clinical trials of controlling the voice onset and offset

Despite emergent progress in many fields of bionics, a functional Bionic Voice prosthesis for laryngectomy patients (larynx amputees) has not yet been achieved, leading to a lifetime of vocal disability for these patients. This study introduces a novel framework of Pneumatic Bionic Voice Prostheses as an electronic adaptation of the Pneumatic Artificial Larynx (PAL) device. The PAL is a non-invasive mechanical voice source, driven exclusively by respiration with an exceptionally high voice quality, comparable to the existing gold standard of Tracheoesophageal (TE) voice prosthesis. Following PAL design closely as the reference, Pneumatic Bionic Voice Prostheses seem to have a strong potential to substitute the existing gold standard by generating a similar voice quality while remaining non-invasive and non-surgical. This paper designs the first Pneumatic Bionic Voice prosthesis and evaluates its onset and offset control against the PAL device through pre-clinical trials on one laryngectomy patient. The evaluation on a database of more than five hours of continuous/isolated speech recordings shows a close match between the onset/offset control of the Pneumatic Bionic Voice and the PAL with an accuracy of 98.45 ±0.54%. When implemented in real-time, the Pneumatic Bionic Voice prosthesis controller has an average onset/offset delay of 10 milliseconds compared to the PAL. Hence it addresses a major disadvantage of previous electronic voice prostheses, including myoelectric Bionic Voice, in meeting the short time-frames of controlling the onset/offset of the voice in continuous speech.


Introduction
The Bionic Voice source is an electronic prosthesis which substitutes for the voice generation function of the missing vocal folds of a laryngectomy patient [1]. The theory of human speech generation is based on a source-filter model in which a source, (generally a combination of aspiration noise and glottal pulse train generated by vocal folds), is filtered through a resonance model of the vocal tract. In that sense, laryngectomy patients lose the ability to generate and control this voice source signal, while retaining a functional vocal tract filter. The primary PLOS ONE | https://doi.org/10.1371/journal.pone.0192257 February 21, 2018 1 / 20 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 A) The hybrid electronic voice source to generate voiced and airflow components. The voiced source is generated by a miniature loudspeaker. A simulated PAL (or vocal folds) model provides the voiced component signal to substitute the voice generation of the mechanical reed element of the PAL. The airflow component is generated by a micro-blower and is proportional to the airflow of the stoma. The resulting voiced and airflow components are combined and applied to the vocal tract.
B) A set of pressure sensors to monitor the intra-oral and stoma pressures (the pressure above and below the missing larynx respectively). As a mechanical reed, the PAL is driven by the variations of these two pressures at its input and output sides [25].
C) A control unit which uses the monitored pressure values to predict the onset and offset instances and the pitch of the voiced component together with the magnitude of the added airflow.
The main reported disadvantage of the traditional PAL device is its cumbersome design [9] and the need for a tube to transfer the sound to the vocal tract resulting in an unhygienic direct coupling of the stoma to the mouth [24]. This is addressed in the design of Fig 2 by providing a wireless link between the pressure sensing unit at the stoma and the intra-oral source. The PAL mechanism of voice generation in speech When used in speech as a voice source, the PAL complies with two main hypotheses that underlie a source-filter separation of speech: A) The length of PAL's resonant tube is small compared to the quarter wavelength of the resonance of the reed. The tube has an acoustic output impedance, much larger compared to the vocal tract [26]. Hence, the PAL reed oscillates close to its natural frequency [27][28][29]. B) The natural frequency of the reed (100 ± 20Hz) is significantly lower than the vocal tract first formant so, the acoustic coupling of the PAL and the vocal tract can be assumed to be insignificant [30,31].
With the source-filter separation assumption in place, similar to the vocal folds, the PAL can be modelled in two different modes of operation in speech: when the source undergoes small-amplitude oscillations in voiced speech and 2) when it reacts to large amplitudes of the driving pressure in the transition between voiced/unvoiced speech. The latter is directly relevant in understanding PAL's mechanism of controlling the voice onset and offset. In both cases, the reed oscillations are driven by the pressure difference between the two sides of the reed P g = P s −P m [32,33] (with pressure inside the mouth P m and the stoma P s ).
The small amplitude, self-oscillation of the PAL reed in voiced speech can be described by a simple lumped element model of a free reed or a single mass-spring model of vocal folds [25,34,35]. Increasing the driving pressure (P g ), the PAL exits the small amplitude oscillation regime and moves toward a switched behaviour turning the voice on or off. A hysteresis effect is experimentally observed in the voice onset/offset transition of the PAL. Such hysteresis phenomenon is well-known to vocal folds vibrations with a higher value of pressure threshold for the onset of the voice compared to the offset [36,37]. The voice onset/offset hysteresis is also experimentally observed in the oscillations of the vocal folds in excised larynx [36,38], the mechanical replica of vocal folds (with or without vocal tract coupling) [30,[39][40][41] and pneumatically driven mechanical free reeds [42].
Lucero et al. have proven both theoretically [43] and empirically [44] that a mechanical replica of the vocal folds still demonstrates the voice onset/offset hysteresis in the absence of vocal tract coupling. With the source filter separation in place, this framework [43] equally describes the case of an ideal PAL where a mechanical replica of the vocal folds (placed external to the body) substitutes the PAL's reed [44]. The significance of this configuration is that the use of a vocal fold replica in an ideal PAL is expected to generate a more natural sound compared to a PAL with a simple reed. Hence, implementing the hysteresis behaviour will remain a general requirement for designing Pneumatic Bionic Voice prostheses to approach the performance of an ideal PAL in controlling the voice onset and offset.

A numerical hysteresis model for the PAL onset/offset control
The PAL device used in this study is the DSP-8 (Fig 3). It has a rectangular reed element vibrating inside a cylindrical cavity and demonstrates a voice onset/offset hysteresis effect similar to the vocal folds [37] when used in speech. To provide a precise model of the onset/offset control of this PAL which can be also extended to other PAL designs with different shapes of the reed (such as the Tokyo artificial larynx [21]), a numerical optimization approach is proposed. Such optimization is specifically advantageous as the mechanical attributes and the elasticity of the PAL reed seems to change over time and the parameters of a hysteresis model to describe the onset/offset control need to be updated.
A bi-state hysteresis model H(ψ,;,P g ) is proposed to be adapted to the time-varying attributes of the onset/offset control of this PAL source. In this model, P g (t), the pressure inside the source that drives the vibrations at time t, is the rectified difference between stoma P s (t) and intra-oral P m (t) pressure values. The rectified difference (half brackets) is used in (1) since the source will not generate any voice for negative values of the driving pressure. The hysteresis function H(ψ,P g ) (2) with ψ = {θ L ,θ H } has an onset threshold: θ H and offset threshold: θ L (θ H > θ L ).
HðtÞ ¼ Hðc; P g ðtÞÞ ¼ 1 where P g > y H 0 where P g < y L Hðc; P g ðt À 1ÞÞ else The hysteresis function output, H(ψ,P g (t)), is smoothed using three time delays, for attack, release, and hold, ; = {T A , T R , T H } (in milliseconds) to from H s (;,t) in (3). When an onset (offset) instance is detected by (2) an attack (release) counter, C A (C R ) is initiated to keep track of the onset (offset) time. The smoothed Hysteresis output is then defined as H s (t) = H s (;,t) [45]: The smoothed model (3) becomes more robust to noise when P g approaches the onset and offset thresholds by delaying the decision by increasing the hold time T H . The main task of the modeling is then to optimize the five parameters of the onset and offset thresholds (θ H , θ L ) and time delays (T A , T R , T H ) when the PAL is used in speech.

Pre-clinical trials of controlling a pneumatic Bionic Voice source onset and offset
The parameters of the model (2, 3) were optimized and its predictions were verified against recordings of the mechanical PAL device through pre-clinical trials. The trials were supported by an ethics approval from the Sydney Local Health District ethics committee, Royal Prince Alfred Hospital zone (protocol number X14-0276, HREC/14/RPAH/362). One laryngectomy patient was enrolled following his expression of interest to the public advertisement of the research. The recruitment continued from January 2015 to September 2018 and written consent was obtained. The participant was a proficient user of the DSP-8 PAL device (Fig 3). The use of only one patient is consistent with the aim of the trials to achieve a precise numerical model of the mechanical PAL device rather than study patients' interaction with the device.

The dataset
The laryngectomy patient sat in a quiet room and used the PAL source as his voice prosthesis to speak, while his supra-glottal and subglottal pressures were recorded simultaneously. Recordings of respiration (the patient's stoma and his intra-oral pressure) and speech were performed at 1 kHz and 48 kHz respectively. The respiration recordings were manually labelled to identify onset/offset instances of the source during speech. For verification, the performance of the model (2,3) in estimating the onset/offset was evaluated against hand-labelled PAL recordings. The criterion for evaluation has been the correlation coefficient between the estimated and target voiced/unvoiced labels. To elicit words and sentences, standard speech tokens used in the evaluation of intelligibility of dysfunctional voice were used [46][47][48]. The spoken phrases were the rainbow passage [49] (first paragraph) at slow and normal conversational speech rates, the Diagnostic Rhyme Test (DRT) [50], Modified Rhyme Test (MRT) and the Phonetically Balanced (PB) word list [51]. The length of each recording was 45 seconds. Overall, more than five hours (20,000 seconds) of data were collected. An anonymized dataset underlying the results of this study is available as a Supporting Information File (S1 File).

Pre-processing
When generating a voice for the patient, the PAL source creates air pressure vibrations which are picked up by the pressure sensors. This contaminates the patient's recorded stoma and intra-oral pressures. To estimate the voice onset/offset from the respiratory data to control the Pneumatic Bionic Voice source in the absence of the PAL, a pre-processing step is required to remove the PAL source vibrations. Considering P s raw ðtÞ, P m raw ðtÞ as the raw recorded stoma and the intra-oral pressures respectively, a low-pass Butterworth filter (with a cut-off frequency of 70 Hz) is applied to these signals to calculate P s (t) and P m (t) in (1).

Threshold optimisation on offline data
Two scenarios (an offline and a concurrent situation) were considered for finding the parameters of the hysteresis model (2,3). Fig 4 shows the threshold optimization for the offline scenario. The preprocessing (low-pass filtering) of raw stoma and mouth pressure recordings is depicted in Fig 4a and 4b. The calculated difference between the low-pass filtered pressure of stoma and mouth, P s (t), P m (t) in (1) was the input to the hysteresis model (Fig 4c). The parameters of the model were optimized for each recording of 45 seconds. The attack, release and hold times (T A , T R , T H ) were initially set to 1ms for all recordings to minimize the response delay of the controller at 1kHz. The two thresholds (θ H , θ L ) were then optimized to minimize the Mean Square Error (MSE) between the manual-labels of each recording of 45 seconds with the model predictions. The search range for (θ H , θ L ) was modified for each recording: using the manual labels, the set of values of P g in (1) at which a voice transitions from voiced to unvoiced or vice versa were selected, and their median and Mean Absolute Deviation (MAD) were calculated. For each parameter (θ H , θ L ), the range of median ± 3 MAD was then divided into 25 segments, and the MSE search was then performed on the resulting 625 (25×25) points. The calculated thresholds together with the pre-set values of (T A , T R , T H ) provided a sufficiently low MSE error so further fine-tuning of (T A , T R , T H ) was not needed. Fig 4d shows the offline MSE threshold fitting.
When fitted individually for each recording, the average correlation coefficient of the predicted and manual voiced/unvoiced labels is 98.45 ±0.54% on the recorded dataset. Offline MSE results are promising, but the approach is computationally expensive and needs the full length (45 seconds) of the recording to maintain its accuracy. In addition, the elasticity of the source changes, resulting in variations of θ H and θ L between different recordings of the PAL source. However, it is safe to consider θ H , θ L constant for the full length of each recording (45s). The attack, release and hold times (; = {T A , T R , T H }) remained consistent along the entire database.

Threshold optimisation on concurrent data
To overcome the shortcomings of the offline MSE approach, a concurrent scenario was proposed to minimize the length of data required for estimating the thresholds of each recording.
The concurrent scenario was designed for online evaluation of the model against the PAL when the patient is using the PAL to speak. The concurrent approach (which may not be necessarily in real-time) has to adapt to the changes in PAL hysteresis thresholds as the new recording data is accumulated. Thus, it can no longer use a manual labelling of the data as opposed to the offline mode.
The block diagram of the concurrent method is reflected in Fig 5. The system uses a lowpass filtering of the pressure values P s (t), P m (t) (as described in the Pre-processing section) to provide input to the hysteresis model. In a parallel pathway, an auto label detection module is designed to provide the ground truth of the onset/offset labels at each time step. The auto-labels are detected using the unfiltered PAL recording of the stoma P s raw ðtÞ (which is contaminated by the vibrations of the source). With the ground truth established, an optimization algorithm adjusts the model thresholds by minimizing a cost function as a function of the thresholds.
Automatic voiced/unvoiced label detection. Since the vocal tract acts as a resonant cavity, the raw intra-oral pressure of the mouth P m raw ðtÞ will maintain resonances even after the PAL source relinquishes vibrations. Hence detecting the onset/offset labels automatically is more precise from the raw stoma recording P s raw ðtÞ ( Fig 5). P s raw ðtÞ shows pseudo harmonic oscillations added to a slow-varying envelope (Fig 6). The auto-label detection removes this slow-varying envelope from the signal. Applying a moving average filter with a 20-sample window to P s raw ðtÞ results in the slow-varying envelope ma Ps raw ðtÞ. Next, the difference signal d t is calculated (4) which shows a pseudo harmonic oscillation pattern around zero when the PAL is generating voice.

d t ¼ P s raw ðtÞ À ma Ps raw ðtÞ ð4Þ
To estimate the onset and offset, the moving standard deviation of the rectified differencẽ d t is defined where w t is the moving average window centered at time t. The difference signal d t oscillates in voiced speech, so the moving standard deviation v t shows a bi-modal histogram divided by a threshold. In the unvoiced speech, v t falls in the lower amplitude range of the histogram, moving towards the larger amplitudes for voiced speech. The adaptive histogram thresholding by Otsu et al. [52] is performed to determine the point to divide the v t histogram into voiced/unvoiced regions. The auto-thresholding method clusters amplitudes of v t in the histogram into two voiced/unvoiced classes. The class with values higher than the Otsu's optimal threshold v t is labelled voiced, and the other unvoiced.
The advantage of the auto-labelling method is that the histogram is dynamic, changing over time as the PAL pressure data is accumulated and the Otsu thresholding [52] is updated over time. This makes the auto-labeler similarly strong in estimating voices/unvoiced instances at continuous speech (Fig 6). When applied to the training dataset, the method provides a high matching of 98.67 ± 0.58% with manual labels for isolated words and 98.2 ±0.9% for the recordings of the rainbow passage at a conversational speech rate.
Concurrent optimization of the thresholds. Starting from the beginning of the recording and without prior knowledge about the recorded data, at time steps of t = T, the optimization of a cost function f t (θ H , θ L ) (6) is performed and the thresholds of the hysteresis model (2, 3) are updated. The output of the model has to match PAL ground-truth labels L t 2 [0,1], which are calculated through automatic label detection (Fig 5). The cost function of this optimization, f t (θ H , θ L ), is defined as the variance of the weighted error between estimated onset/offsets of the hysteresis model H(ψ, P g (t)) and the ground truth (L t ). f t ðy H ; y L Þ ¼ Var t¼1:T ½w t ðL t À Hðc; P g ðtÞÞÞ ð6Þ The error is weighted by w t , i.e., the rectified driving pressure of the PAL, giving more significance to errors observed in larger values of the driving pressure P g (t).
The optimization algorithm searches for the threshold values θ H , θ L to minimize f t (θ H , θ L ) at each time step. To facilitate the convergence of the algorithm and reduce the execution time, the history of ground truth labels L t is used to limit the range of values searched for (θ H , θ L ) at each time step. Despite being bistable, if the hysteresis system (2, 3) detects onset at P g (t) = θ H0 , it will remain in a stable "on" state for increased values of the driving pressure. These values (P g (t) > θ H0 and rP g (t) > 0) should be ruled out from the search domain of the upper threshold θ H for the next time-step. A median of the elimination candidates is chosen to limit the upper bound of the search domain. For the lower threshold θ L (which has a smaller dynamic range), the search range is proven sufficient to be limited to {0, median (θ Lt )} where θ Lt are the calculated lower thresholds of previous time steps.
The cost function f t (θ H , θ L ) (6) is non-smooth, discrete and non-differentiable so derivative based approaches may not necessarily converge during optimization. Fig 7 explains the case for convergence of an optimization algorithm in the search domain to minimize f t (θ H , θ L ). As observed in this figure, the variations of cost function f t (θ H , θ L ) in the recorded dataset with changing thresholds show a global minimum. Yet, the presence of local optima supports the choice of non-gradient based or heuristic optimization techniques. Accordingly, two direct search methods are employed for threshold optimization which do not require any derivative information (explicit or implicit) from f t (θ H , θ L ). These are the pattern search [53] and the Nelder-Mead method [54]. The Nelder-Mead method is a widely used multidimensional minimization and is inherently unconstrained. The search domain of the pattern search, however, was bounded as described before.
In both approaches the optimization starts from initial values of the thresholds θ H0 and θ L0 and recursively updates these, using the thresholds values derived at a previous step. At   Fig 7. The logarithmic error, cost function f t (θ H , θ L ), calculated over the search domain for θ H , θ L for the dataset.
https://doi.org/10.1371/journal.pone.0192257.g007 each iteration step (t), the pattern search algorithm places a stencil (pattern) centred around the values of the thresholds derived in iteration (t−1) in the R = {θ H , θ L } domain. This pattern includes a set of search directions to cover the points adjacent to the center. If the cost function value f t (θ H , θ L ) is decreased in any of these directions, the pattern is moved to the new center point, otherwise, the stencil size is reduced. The pattern search is terminated when the iteration error is less than the desired accuracy, or the available number of iterations has expired. Fig 8 demonstrates the convergence behavior of the methods when applied to 24 isolated and continuous speech recordings of 45 seconds sampled at 1kHz with a total of 1,080,000 decision points. The accuracy is measured by the percentage of correctly estimated onset samples over the total number of samples. The thresholds have been updated at intervals of T = 1s. To avoid any confusion, the reported accuracy is against original manual labels in this figure.
As observed in Fig 8, the overall accuracy of the optimization method for pattern search and bounded pattern search is better than the unconstrained Nelder-Mead method with bounded pattern search slightly higher than pattern search. The algorithm needs at least 5 seconds of recording of the patient at 1kHz (using the PAL source to speak) to reach 98 ± 2% of accuracy. However, to reduce the standard deviation of the error for all recordings, 10 seconds of recorded data is needed to limit the variation of the accuracy to remain close to 98.3±1% approaching to 98.7 ± 0.4 at 15 seconds.  Curves are constructed by first finding the global optimum for both thresholds and then varying each while keeping the other one fixed at its optimal value. A reverse of this strategy has been adopted in the pattern search algorithm to derive the thresholds by optimizing one of the thresholds first and next using its value to derive the optimal value for the other.

Real-time performance of the voice onset/offset control
To assess the real-time performance, the control algorithm of the Pneumatic Bionic Voice prosthesis (with fixed thresholds) was implemented in parallel with a hardware emulation of the PAL. The PAL hardware emulation was implemented using an Arduino Uno™ which streamed out the previously recorded pressure values of the PAL (raw pressures of the mouth and stoma) at 1kHz. The manual labels of the respiration data were loaded to the Arduino™ to generate a beep sound in the voice onset states. The Pneumatic Bionic Voice control was implemented using C++ and compiled on a Mac computer with 64 GB RAM and a 3.5 GHz 6 core Intel™ Xeon E5 processor and similarly generated a sine wave sound whenever a voice onset was detected in the real-time streamed respiration data. The two audio outputs (from the PAL emulator and the Bionic Voice control algorithm) were isolated and simultaneously recorded in two separate channels of an audio stereo recorder. The measured delay between Bionic Voice and the PAL emulator varied in the range of 5-17ms with an average of 10ms which was consistent between the onset and offset.

Discussion
This study defines a new framework for using PAL as an excellent reference for designing a Pneumatic Bionic Voice prosthesis. As the first step, a precise model of the onset and offset control of a Pneumatic Bionic Voice prosthesis is provided with an accuracy of 98.45 ± 0.54% when compared to the PAL performance in a mixed (continuous and isolated) speech database. The onset and offset thresholds, mainly depend on the physical attributes of the PAL reed. However, these may also be speaker-dependent as each speaker uses the PAL in a different range of respiration pressure values. As the attributes of the reed change with time, threshold determination will be session-dependent as well. The fact that the optimal thresholds θ H , θ L of the offline scenario remain constant in the duration of each session (of 45 seconds) enables source implementation with fixed values of these optimal thresholds. Having uniform thresholds allows the system to perform faster in real-time. Nevertheless, a concurrent threshold optimization is also implemented to enable online evaluations of the source against changes of the PAL attributes over time.
In terms of an intuitive pneumatic control of a voice prosthesis, this study complements Takahashi et al. [55] who implemented a real-time voicing control of an Electrolarynx which terminated the voice after observing a peak in the intra-oral pressure. The concept [55] is similarly valid for the PAL as an increase in the intraoral pressure decreases the driving pressure of a PAL and leads to a voice offset. The real-time voicing control approach of Takahashi et al. [55], has a delay of less than 20ms. The system was tested on one laryngectomy patient and reported to correctly detect 90% of unvoiced consonants, improving the misidentification of patient's voiced/unvoiced stop consonants by 50% [56].
The other significant research direction in intuitive control of an artificial voice source is the myoelectric Bionic Voice prostheses. These are essentially motor control prostheses similar to a Bionic Arm, in which the intended motor function is voice generation. Goldstein et al. [57,58] established the concept of using surface electromyography (sEMG) of neck strap muscles as a reliable indicator of the voice onset and macro-variations of the pitch of the voice. Fuchs et al. [59] improved these results using a database of sEMG/speech recordings of 19 minutes of phonetically balanced sentences of healthy male and female subjects. With an approach similar to the benchmark study of Goldstein et al. [58], they reached an average onset/offset error range of 9.5 ±5.6% and 6.6 ±4.7% for the female and male subjects respectively. However, the original algorithm [58] provided an accuracy significantly lower than this for laryngectomy patients [46,58].
Using myoelectric Bionic Voice prostheses to control the voice onset/offset faces two major challenges ahead. The first is accessing neck strap muscles to control the prosthesis since these muscles are normally excised at the time of laryngectomy to minimize the risk of cancer spread [60]. Heaton et al. [61] have demonstrated the feasibility of preserving these muscles although this will require a modification of the standard laryngectomy surgical procedure. Stepp et al. [62] demonstrated the utilization of other residual face/neck muscles for non-invasive myoelectric control of a voice prosthesis, obviating the need for modification of the standard laryngectomy surgical procedure. However, these alternative muscles have not yet provided a reliable substitute for the neck straps in terms of controlling voice onset/ offset and the pitch [62].
The second challenge is the significantly larger time scales of controlling muscle contraction/release compared to the real-time control of voice. A myoelectric controller for Bionic Arm can have a reaction delay of up to 300ms before being perceived as sluggish by the user [63]. A window size of 150-250ms is recommended as optimal when segmenting the sEMG data to control a Bionic Arm in a simple two-class task [64,65], resulting in an optimal (minimum) controller delay of 100-125ms [64,66]. These are much larger than the voice onset and offset time frames of continuous speech where the delay between the start of a voiceless consonant and the start of voicing for the next vowel is in the range of 10-70ms [55,[67][68][69]. The inability to match this temporal resolution will translate to the voicing of unvoiced phonemes, i.e., the largest source of reduced intelligibility of EL speech [70][71][72]. Meeting these time limits has been particularly challenging in the myoelectric Bionic Voice when a subject is asked to actively control voice termination by relaxing their neck muscles at the end of a phrase [46,73] (leading to large voice termination delays of 400-700ms for healthy subjects [74] and 1120-1970ms for the laryngectomy participant [46]). Myoelectric Bionic Voice sources, however, seem to be reliable controllers of the pitch (macro-prosody) of the voice which works on a larger time scale at sentence level [74].
These two challenges make the Pneumatic Bionic Voice prosthesis a convenient substitution for the myoelectric Bionic Voice prostheses in terms of onset/offset control. 1) It is a noninvasive, prosthesis which does not require modifications in laryngectomy surgery, 2) In terms of the prosthesis controller delay, the PAL already provides voice onset/offset control in slow and fast speech rates with high intelligibility scores [6-13, 21, 22]. The authors have also achieved a 10ms average delay in a real-time implementation of Pneumatic Bionic Voice onset/offset control when compared to the PAL. A direct comparison is not possible due to variations between recorded data sets. However, the precision and speed of the Pneumatic Bionic Voice on the large recorded corpus seem to be uniquely ahead of previous intuitive voice prostheses which enabled a laryngectomy patient to directly control the onset/offset of the voice in continuous speech in real-time [46,55,75].
It is worth mentioning that, despite its high quality and the ability to provide pitch variations at syllable level [8,9], the PAL device (Fig 3) used in this study as the ground truth is a simplified version of an ideal PAL (where a mechanical replica of the vocal folds substitutes the PAL's simple reed). Hence following an ideal PAL as the reference may improve the results of this study in terms of estimating a speaker's intended voice onset/offset and the pitch. In addition, even in an ideal configuration, the PAL will lack the potential to modulate the macro-prosody (pitch variation at sentence level). The respiratory control (of the subglottal pressure) may need to be combined with laryngeal control (of the tension of vocal folds) to enable modulating macro-prosody [76]. Meltzner et al. [4] have reported that natural modulation of the pitch is the most influential attribute to improve the quality of speech after laryngectomy [4,77]. This limitation seems to be the strength of the Myoelectric Bionic Voice solutions [78,79] and potentially ties the two research directions of Pneumatic and Bionic Voice prostheses together to provide an intuitive control over the onset/offset and the pitch of the voice.

Conclusion
The PAL can be considered as a simple model of the human larynx with a fixed pair of vocal folds driven exclusively by the variations of the intraoral and subglottal pressure values and without any neural/neuro-muscular input from the missing larynx. The quality of PAL speech is comparable to the existing gold standard of TE voice prostheses and far better than the Electrolarynx [6][7][8][9][10][11][12]. The traditional PAL also holds a significant advantage over the existing gold standard as being non-invasive.
These advantages advocate defining a new pathway in designing Pneumatic Bionic Voice prosthesis as electronic adaptations of the PAL. This study aims to be the first in this direction and provides a precise model that describes the PAL voice onset/offset control with a low computational cost suitable for real-time implementations. The next step for the authors is to combine this solution with a PAL pitch modulation model in real-time and evaluate the quality of the resulting speech against the PAL and the existing gold standard.
Supporting information S1 File. This is the dataset of recorded intra-oral and stoma pressure values, manually labelled for voiced/unvoiced detection. The zip archive contains anonymized recordings of time-aligned intra-oral and stoma pressure (P m raw ðtÞ, P s raw ðtÞ respectively), with their voice/ unvoiced labels (L t ) as the ground truth. (ZIP) Validation: Farzaneh Ahmadi, Farzad Noorian, André van Schaik.