Kinematic model for simulating mucosal wave phenomena on vocal folds

Biomedical Signal Processing and Control 49 (2019) 328–337

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Kinematic model for simulating mucosal wave phenomena on vocal folds ∗ ˇ S. Pravin Kumar, Jan G. Svec Voice Research Lab, Department of Biophysics, Faculty of Science, Palack´ y University, Olomouc, Czech Republic

a r t i c l e

i n f o

Article history: Received 27 June 2018 Received in revised form 12 October 2018 Accepted 6 December 2018 Keywords: Mucosal wave Vocal fold vibrations Kinematic model Synthetic kymogram Videokymography

a b s t r a c t Mucosal waves have been found to be important for evaluating vocal fold vibrations in laryngological practice. While they are routinely evaluated visually, the knowledge on the physical phenomena related to mucosal wave propagation is limited. Kymographic imaging, in particular, reveals various mucosal wave features that deserve more understanding in order to advance functional diagnostics of voice disorders. Here, a kinematic model is presented which simulates mucosal waves on human vocal folds. The vocal fold geometry is based on a parametrically adjustable M5 model. A kinematic rule is used for simulating the propagation of the mucosal wave from the bottom of the vocal folds upwards and laterally over the upper vocal fold surface. The model maps the changes of the coronal shape of the vocal folds through vibration cycles. The vibration characteristics including the mucosal wave movements are then visualized using a synthetic kymogram graphically obtained through a local illumination method. The model can serve as an educational and research tool for studying the mucosal wave features and their appearance in laryngeal kymographic images. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction “Mucosal wave” (MW) is rather a general term used to describe the movement of the vocal fold mucosa during laryngeal voice production. MWs originate inferiorly and propagate upwards and around the vocal folds during phonation, creating a wave-like motion on the vocal fold surfaces [1–13]. Besides human vocal folds, mucosal waves were documented on the vocal folds of many mammals and recently also in birds [14–19]. The appearance of mucosal waves was shown to be a crucial component of the myoelasticaerodynamic theory of phonation, explaining the mechanisms of the self-sustained vocal fold oscillations [18,20–22]. In laryngology, MW is a clinically important feature revealing on the pliability of the mucosa as well as on mucosal pathology, which influence the primary voice quality [23–26]. MW analysis is useful in the diagnosis of different vocal fold lesions and has also been used, e.g., for evaluating the success of treatment procedures, hydration effects, or mucosal healing after phonosurgery [27–29].

∗ Corresponding author at: Voice Research Lab, Department of Biophysics, Faculty of Science, Palacky´ University Olomouc, 17. listopadu 12, 771 46, Olomouc, Czech Republic. ˇ E-mail address: [email protected] (J.G. Svec). https://doi.org/10.1016/j.bspc.2018.12.002 1746-8094/© 2018 Elsevier Ltd. All rights reserved.

Diverging definitions are found in the literature [23,24,30], which suggests that the term “mucosal wave” may be used in different contexts by different professional groups. Scientists usually understand the MW as a phase lag of the movement of the upper margin versus the lower margin of the vocal folds [6,21,31–33]. These vertical phase differences are actively driven by the airflow and are critical for energy transfer from the airflow to the vocal folds [21,33,34]. Clinicians, on the other hand, typically refer to MWs as an extent of displacement of the upper margin of the vocal folds (Fig. 1), mostly because routine laryngoscopic examinations using stroboscopy or high-speed video reveal only a portion of the MW propagating on the superior surface of the vocal folds [35–39]. Videostroboscopy is currently the most commonly used method for examining the vibratory patterns of the vocal folds. However, due to sampling rate limitations, it may not accurately capture the true MW movements [40–43]. Alternatively, high-speed videoendoscopy (HSV) with frame rates above 4000 images per second may be used to accurately register these movements [44–47]. However, it is a time consuming task to search every frame for the visual clues of the MW movements from the HSV recordings [48,49]. Videokymography (VKG), on the other hand, is a single-linescanning high-speed imaging technique which conveniently documents MW within a single image (e.g., a kymogram) in real time rather than using a video segment [25,26,50–52]. Alternatively, software-based approaches, such as digital kymography

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

329

Fig. 2. VKG images showing different MW features: (a) mucosal waves extending all over the upper vocal fold surface, (b) MW present but the extent of amplitude is reduced on both vocal folds and (c) MW absent on both vocal folds. Symbols used: rmw = right mucosal wave, lmw = left mucosal wave, lf = left vocal fold, and rf = right vocal fold. Total time displayed in the VKG images: (a) 22 ms, (b) 14 ms, (c) 18 ms. Fig. 1. Laryngostroboscopic image (Kay Elemetrics demo example of stroboscopy) showing upper margin (um), lower margin (lm) of the vocal folds.

(DKG) and strobovideokymography (SVKG), are available for extracting kymographic images from the previously recorded high-speed and stroboscopic videos, respectively [53]. Mucosal wave kymography (MKG) is another high-speed visualization technique based on image processing of high-speed laryngeal videos and digital kymograms, which can be used in documenting MW movements [24]. MKG pixel intensity encodes the MW velocity, and the color codes distinguish the phases (opening and closing) of motion. There have been some efforts in advancing this technique further by exploring the optical flow methodology in the analysis of the high-speech images [54]. Kymograms are useful in distinguishing the MW as being present, reduced or absent (Fig. 2). The presence of MWs indicates health and pliability of the vocal folds [26,33,55,56]. In contrast, reduced or absent MWs are typically a sign of reduced pliability of the vocal fold mucosa and thus may indicate the presence of vocal fold lesions and scarred tissue [24,57]. The kymographic techniques (e.g., VKG, DKG, SVKG, and MKG) are capable of documenting the two basic features of the MW, i.e., (a) vertical phase differences in the vocal fold movement and (b) lateral propagation or passive continuation of the flow-induced mucosal waves. Vertical phase differences are revealed in the kymograms as double edges during the closing phase of the glottal cycle (Fig. 3a) [24]. However, in many kymographic images, the upper and lower margins are not clearly distinguishable (Fig. 3b). In these cases, nevertheless, the phase differences between the upper and the lower vocal fold margins are being reflected in the sharpness of the lateral peaks shown in the kymogram [26,58]. Large or reduced vertical phase differences are seen as sharp or rounded lateral peaks in the VKG images, respectively (Fig. 3b and 3c). A rounded lateral peak is a possible indicator of a stiffened mucosa, inflammation, or scarring. The shape of lateral peaks has been reported to be a useful visual feature for assessing voice disorders using videokymography [59]. It can be seen as an easier, alternative feature for visually evaluating the vertical phase differences than the contours of the lower or upper margins. Although MWs are effectively assessed through visual inspection of the kymographic images, the physics of mucosal wave

Fig. 3. Vertical phase differences (VPD) between lower margin (lm) and upper margin (um) are reflected in the VKG images as (a) double contours during the closing phase or (b) sharp lateral peaks when the VPD is large and (c) rounded lateral peaks when the VPD is small. Total time displayed in the VKG images: (a) 18 ms, (b) 26 ms, (c) 18 ms.

Fig. 4. Water particles undergo elliptical trajectory when the wave travels through the water surface.

propagation is not yet completely understood. McGowan [34] compared the MWs on the vocal folds to waves traveling on the surface of the water. The water wave is the combination of both longitudinal and transverse motion causing the water particles to undergo elliptical trajectories (Fig. 4). The elliptical trajectories and the phase differences between the particles in motion give an impression of a traveling wave on the water surface. The vocal fold epithelium and the Reinke’s space beneath can be considered to be analogous to the water surface layer and water, respectively. Berry, Döllinger, Bösendecker, et al. [10,11,60,61]

330

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

experimentally observed that the flesh points on the medial surface of excised canine vocal folds also exhibited approximately elliptical trajectories when they were tracked using high-speed imaging techniques. Mucosal vibration amplitudes and MW velocity are important parameters influencing the MW propagation on the vocal fold surface. The mucosal vibration amplitudes were reported to be the largest around the upper margin and decrease caudally [10,11,60,61]. The velocity of the MW traveling upwards over the medial vocal fold surface is directly related to the vertical phase differences between the motion of the upper and lower vocal fold margins. MW velocities were reported to be between 0.3 and 10 m/s and vary with the fundamental frequency of voice [11,33,62–68]. The phonation threshold pressure and the energy transferred from the aerodynamic driving force to the vocal fold tissues are influenced by the MWs [33]. Numerical models of the vocal folds and simulated synthetic kymograms allow visualizing the effect of manipulating a wide range of known physical parameters. The knowledge gained from such experiments is helpful in identifying and understanding the driving parameters of complex vocal fold dynamics. Further, it facilitates the effective interpretation of real kymograms [69]. Synthetic kymograms were used previously to document the resulting vibrations of a phase-delayed overlapping sinusoidal model [69], an aeroelastic model [70], and a finite element model [71] of the vocal folds. The present study aimed at creating a special model for simulating the MW kinematics and displaying the resulting motion of the vocal folds in the form of a digital kymogram. The parametric M5 model [72,73] of the human vocal folds was used to simulate the vocal fold geometry. The MW propagation along the model surface follows the kinematic rule described by Titze [21]. This rule was extended here to include a) vertical motion and b) lateral propagation of the waves over the upper surface of the vocal folds. A local illumination method based on the Lambertian shading model [74,75] was used in the generation of the synthetic kymogram capturing the motion of the model.

2. Method 2.1. The kinematic mucosal wave model Custom scripts written in MATLAB R2016a (MathWorks, 2016) were used for constructing the vocal fold model (VFM), simulating the MW kinematics and generating the synthetic kymogram.

The vocal fold geometry was constructed based on the specifications of the M5 model [72,73]. It was defined in two dimensions in a coronal plane by a combination of lines and arcs and was represented by surface points sampled uniformly with an equidistance of 0.01 cm. The anterior-posterior dimension was neglected here, for simplicity. The VFM points were defined in Cartesian coordinates for a single vocal fold (Fig. 5). These coordinates were simply replicated in y -axis and reversed in x -axis using another set of points in order to define the vocal fold of the opposite side. Thus, the two vocal folds were symmetric to each other. By varying the glottal angle, it was possible to gradually alter the shape of the vocal folds between divergent, uniform, and convergent. Vertical MW velocity, c, is known to be related to the vertical phase delay  of the upper vocal fold margin behind the motion of the lower vocal fold margin, expressed in degrees [33]: c = 360 fo T/

(1)

In our case, T is the vocal fold thickness as defined in the M5 model [72,73], i.e. the vertical distance between the upper (U2 ) and lower (L) margins on the vocal fold surface (Fig. 5), and fo is the fundamental frequency of vocal fold vibration. The trajectory of the motion of the surface points of the vocal folds in x (horizontal, lateral-medial) and y (vertical, caudal-cranial) direction can be approximated as



xn = x0n + An sin (2f o ) t −



d 

yn = y0n + aAn cos (2f o ) t −

n ,

c

d  c

(2)

n , n = (1, 2, . . .N)

where t is the time instant, n is a counter for the sample points defining the vocal fold surface, (x0n , y0n ) are the initial coordinates of the vocal fold surface points prior to MW propagation, d is the distance between the samples on the vocal fold surface (0.01 cm here), An is the vibration amplitude in the horizontal direction and a is a factor relating the vertical vibration amplitude to the horizontal amplitude. Time is incremented in steps of 1/sampling rate till it reaches the desired upper limit. Here, the sampling rate was set to 7200 Hz (identical with the VKG frame rate). The default upper time limit was set to be the same as the duration of one videokymographic frame, i.e. 40 ms, corresponding to a rate of 25 frames per second (CCIR standard). The vocal fold model surface was divided into 4 regions. The first region was below the so-called mucosal upheaval point (the

Fig. 5. Kinematic mucosal wave model of the human vocal folds. (a) prephonatory position; (b) model shape during vibration. The circles in panel (b) show trajectories of selected points of the model surface. The vibration amplitudes increase from the upheaval point upwards to the lower vocal fold margin. The prephonatory geometry of the model in (a) is based on default M5 parameters [73] with the vertical thickness of 0.2 cm, glottal convergence angle of 0◦ , and prephonatory glottal width of 0.07 cm. The horizontal size of the model is 1 cm. Indices used: Z = mucosal upheaval point at which the subglottal mucosa starts vibrating, L = vocal fold lower margin, U = vocal fold upper margin (U1 and U2 refer to the indices where the curvature of upper margin starts and ends, respectively), S = most lateral point of the vocal fold upper surface, W = width of the upper vocal fold surface, T = vertical thickness of the vocal fold, and r = radius of curvature at the upper margin.

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

331

standard deviation of the Gaussian function, which is related here to the extent E of the laterally travelling MW through the formula: =

Fig. 6. Light conditions considered for producing the kymographic image of the vocal fold model surface. Diffuse reflections from the slopes mab , mbc and mcd between the exemplary VFM points a − d are considered here. The slopes determine the light intensity which is virtually detected by a simulated camera positioned above the vocal folds. (Note that only a few VFM points are selectively shown here in order to emphasize the slopes between them – the VFM points in the model were set to be much denser, as shown in Fig. 5).

E 2



W +r

 2



−1

(4)

where W denotes the width of the upper vocal fold surface (horizontal distance between points U1 and S, see Fig. 5) and r is the radius of curvature of the upper margin (exit radius) defined in the M5 model. For most of the demonstrations in this paper (except of Figs.12 and 13) very long mucosal wave extent was used by setting E = 100, which theoretically allowed the MW (if generated) to travel over a surface 100x larger than the width of the upper vocal fold surface before its amplitude decreased to ca. 14% (e−2 ) of its initial value. This resulted in the MW propagating over the entire upper vocal fold surface till the lateral end without noticeable amplitude declination, as shown, e.g., in Fig. 8e or Fig.11. The Eq. (2) is in the form of a parametric wave equation with a phase delay term. The motion along the surface exhibits both vertical and horizontal components and approximates circular (for a = 1) (Fig.  5b) or  elliptic   (for a =/ 1) trajectories [10,60]. The phase delay term t − dc n defines the wave travelling across the vocal fold surface over the distance d*n. The mucosal wave speed c is directly related to the phase differences  between the upper and lower vocal fold margins through a Eq. (1). The glottal midline was placed exactly at x = 0 cm and the VFM points were constrained from crossing this boundary during the vocal fold collision. 2.2. The simulated kymogram

boundary between the vibrating and non-vibrating subglottal vocal fold surface, see Fig. 5). In this region, the vibratory amplitudes were set to zero [7,32,76]. The second region was from the mucosal upheaval point up to the lower margin of the vocal folds where the amplitude was considered to be linearly increasing. The third region included the medial surface between the lower and upper vocal fold margin, where the amplitude of vibration was determined by linear interpolation between the amplitudes of the lower and upper vocal fold margins. The fourth region extended from the upper margin laterally over the upper vocal fold surface. Here, the MW amplitude was allowed to progressively decrease thus simulating the chosen mucosal wave extent. Mathematically, the vibration amplitudes in these regions were defined according to the equations An = 0 ∀n ∈ (K, Z), An =

AL .(n − 1) ∀n ∈ (Z, L), L−Z −1

An = AL +

AU1 − AL . (n − 1) ∀n ∈ (L, U1 ) U1 − L − 1

1 An = AU1 e 2 −

(3)

 nd 2 

∀n ∈ (U1 , S)

where n is the counter which increases from the bottom of the vocal fold surface upwards and sidewards, K is the index of the lowermost point of the model surface, Z is the index of the mucosal upheaval point at which the subglottal mucosa starts vibrating, L is the index of the vocal fold lower margin at which y = 0, as defined by Li et al. [72], U1 is the index of the start of the vocal fold upper margin curvature, and S is the index of the most lateral point of the vocal fold upper surface. The surface points belonging to the indexes used in Eq. (3) are shown in Fig. 5a. AU and AL denote the amplitudes at the upper and lower vocal fold margins respectively (for most of the demonstrations, except of Figs. 12d and 13 a, they were set to be equal), d is the distance between the adjacent samples on the vocal fold surface defined earlier. The parameter  denotes the

Videokymographic images were simulated as if the vocal fold vibration was viewed from above with a 2nd generation VKG camera [50,77]. In the VKG examinations, the light from a laryngoscope illuminates the vocal fold surface and a high-speed VKG camera attached to it registers the light reflected from the vocal fold surface along a single scan-line and stacks consecutive scan-line image strips together to form the VKG image. The default VKG rate is 7200 image lines/s and there are 288 image lines per kymographic frame [77]. Each of these image lines is duplicated here in order to imitate the behavior of the 2nd generation videokymograpic camera which stores each scan line twice, in order to comply with the CCIR/PAL standard requirements [77]. To obtain the simulated VKG image, a local illumination model, based on diffuse reflection, was adopted to simulate illumination of the VFM and to calculate the amount of light reflected from its surface [74,75]. The VFM surface was presumed to be Lambertian (i.e., a perfect diffusely reflecting surface which scatters incident light equally in all directions [78]). Therefore, its brightness does not change with viewing direction [79] and depends only on the angle of incident light and the distance of that surface from the light source. The model surface was two-dimensionally represented by the slopes between the VFM points. A virtual laryngoscopic camera-view was established, orthogonal to the superior surface of the VFM. For brevity, a few selected VFM points are shown in the illustration on Fig. 6. Here, L is the unit vector pointing to the direction of the light source and towards the virtual camera, and N is the unit vector normal to the surface. The intensity of diffuse reflection from the VFM according to Lambertian law depends on the surface normal to the slope between the VFM points and the respective intensity of the point light source Ip . This can be modeled as I=

kd Ip dy2

(L.N) =

K dy2

cos

(5)

where K = kd Ip , kd is the diffuse coefficient of reflection, Ip is the intensity of the light source, dy is twice the vertical distance of light

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

332

Fig. 7. Mapping the VFM physical coordinates to the virtual videokymographic camera’s pixels. The vertical lines show the boundaries of the individual pixels. A low pixel resolution of Npx = 50 pixels per 0.5 cm vocal fold surface is adopted here for illustration purposes. Representative synthetic kymograms were generated with a pixel resolution of Npx = 100 pixels per 0.5 cm vocal fold surface. The VFM-mapped pixel values (from the point of view of a virtual line camera on the glottal axis orthogonally above the VFM) were algorithmically computed at each time instance and stored as single line images (shown on top of the image). The surface illumination was used for determining pixel values. Here, an exemplary physical range bounded by bi and bi + dw is magnified to show the slopes which were considered in calculating the light intensity detected by a single pixel.

source from the VFM point (in our model, the virtual light source was chosen to be 7 cm away from the uppermost surface of vocal fold before it was set to vibration), and  is the angle of the surface with respect to the position of the light source and the camera. Here, the light source was assumed to be placed vertically above the vocal folds and the camera orientation was considered to be the same as the direction of the light source. The angle  can be calculated from the slope between the VFM , where y and x are the vertical and points i.e.,  = tan−1 y x horizontal differences between the coordinates of the adjacent points on the VFM surface. The exemplary VFM points and their slopes in Fig. 6 show the relationship between the L vs. N angle and the corresponding intensity of the reflection expressed on a grayscale (i.e. scale between black and white). The value of K = 0.97 max(dy2 ) was empirically determined to provide good-looking kymograms. In order to generate the synthetic kymogram, the physical coordinates of the VFM points in the x -axis were mapped onto the image coordinates of the virtual VKG camera with the total of Npx pixels in a row (Fig. 7). The virtual camera was set to cover the vocal fold region horizontally (along the medio-lateral axis) from -0.5 cm to 0.5 cm, which included most of the upper vocal fold surface. At a given time instance, pixel array elements (columns) from 1 to Npx with an incremental step size of 1 were mapped onto the x -axis of the VFM from -0.5 cm to 0.5 cm with an incremental step size of dw cm, corresponding to an image size of 1 pixel. pi → [bi < ∀x ≤ bi + dw] , b = [−0.5, −0.5 + dw, . . .0.5] , i



= 1, 2. . ..2Npx



(6)

where pi is a single pixel mapped onto the range bi to bi + dw, in the x –-axis. At each time instance, single line kymographic images were computationally generated by dynamically mapping the momen-

tary surface slopes to pixels. Because the horizontal distance between the neighboring surface points was generally different from the pixel image size, there was a need to determine representative surface slope values for each pixel. The pixel intensity values in a mapped physical range were calculated from the intensities of diffuse reflection from the surfaces between the VFM points and the normalized differences of their positions in the x -axis. Calculation of the intensity value from an exemplary range of bi to bi + dw, corresponding to the size of one pixel, is illustrated in Fig. 7. The pixel value corresponding to the physical range bi and bi + dw was obtained as Pi =

 1  d I + dfg Ifg + dgl Igh dw uf ef

(7)

where dj is the horizontal distance of a VFM point in x-axis from the neighboring point or from the physical bounds, andIj (j ∈





ef, fg, gh ) is the normalized reflected light intensity calculated





from the Eq. (5) using the term j = tan−1 mj , where mj is the surface slope. A visibility constraint was applied so that the VFM points hidden from the laryngoscopic view were not included in the computation. Depending on the phases of the glottal cycle, there were some instances when no VFM points existed in the lateral and medial ends of the mapped physical coordinates. If such a case occurred at the lateral ends (usually during the closing phase), the respective pixel arrays were assigned with the immediate non-empty adjacent pixel values. If there were no medial VFM points (i.e., when glottal gap was present between the vocal folds) then the respective pixel arrays were assigned with zeros representing black color. Since the VFM motion was symmetric, it was sufficient to calculate the pixel array values (1 to Npx ) for one side of the VFM alone. These values were then simply reversed and duplicated in the pixel arrays (Npx + 1 to 2N px ) mapped to other side of the VFM.

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

333

Fig. 8. Exemplary model kinematics selected at the time instances 8.7 ms, 10.6 ms, 12.1 ms, and 15.0 ms to illustrate the (a) opening, (b) maximally open, (c) closing, and (d) complete closure phases of the glottal cycle and their respective computationally generated single line images A, B, C and D. The kymographic image (e) was obtained by concatenating the single line images from all the successive time steps. The dotted white lines in panel (e) show the location and time instances of the exemplary line images A, B, C and D. To generate the VFM motion, the following VFM parameters were used: fo = 100 Hz, Au = 0.07 cm, a = 1, and  = 90◦ . The resulting mucosal wave velocity was 0.8 m/s. The M5 parameters were the same as indicated in Fig. 5.

Fig. 9. Montage of selected VFM frames illustrating MW dynamics during the closing and closure phases of the glottal cycle - the mucosal wave peaks highlighted by dashed circles progressively move towards the lateral ends in subsequent frames. This phenomenon was reflected in the lateral shifting of white pixels in the consecutive line images. The VFM settings were identical to those in Fig. 8.

At each given time instance, pixel values for all the mapped physical ranges were stored in a single row of the synthetic kymographic image. Subsequent rows were then concatenated together to form the final image.

3. Results Fig. 8 demonstrates the shapes of the VFM at four time instances corresponding to the opening, maximally open, closing and

334

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

Fig. 10. Model kinematics with highly reduced MW during the (a) opening (b) maximally open (c) closing and (d) complete closure phases of the glottal cycle, and (e) the resulting kymogram. The model shows rather uniform movement of the vocal fold surface with no visible signs of mucosal wave traveling laterally across the upper vocal fold surface. In the kymogram, notice also the absence of visibility of the lower margin, as well as the rounded shape of the lateral peaks. The model settings are identical to those in Figs. 8 and 9, except for the mucosal wave velocity which was increased about 9 times (from 0.8 to 7.2 m/s) by changing the vertical phase difference  from 90◦ to 10◦ .

completely closed phases of a glottal cycle, as well as their respective algorithmically computed line images A, B, C and D (Fig. 8a–d). The synthetic kymogram, obtained by sequentially concatenating all the single line images (generated at all time instances), is shown in Fig. 8e. The start and maximum separation of the VFM edges from each other were seen as opening (Fig. 8a) and maximally open phases (Fig. 8b) of the vibration cycle, respectively. During the closing phase, the VFM edges moved medially towards each other, with the motion of the lower margin preceding that of the upper margin, as expected (Fig. 8c), thus exhibiting vertical phase differences. During the closing phase, the mucosal wave propagated from the upper vocal margin laterally on the superficial layer of the VFM (Fig. 9). This was reflected in the kymographic line image as the spatial shifting of the white pixels from the glottal edge laterally. The mucosal wave velocity and the phase difference between the lower and upper margin can be varied by changing the values of . For instance, in Fig. 8, the chosen vertical phase difference of 90◦ between the lower and upper margins determined the mucosal wave speed c = 0.8 m/s and the resulting synthetic kymogram exhibited sharp lateral peaks and laterally traveling MWs. In contrast to Fig. 8, Fig. 10 shows the MW kinematics with a greatly reduced phase difference (10◦ ) between the motions of the lower and upper margins, making the mucosal wave velocity 9 times higher, while the remaining parameters were the same. Due to the very small vertical phase differences, the vocal fold margin moved rather uniformly. In the resulting synthetic kymogram (Fig. 10e), the lateral peaks appeared rounded and the lateral propagation of MW was not evident. Fig. 11 demonstrates different degrees of sharpness of the lateral peaks (sharp towards rounded) in the kymograms due to the increase in the MW speed by reducing the vertical phase differences  (recall the relationship in expression (1)). Besides of the rounding of the lateral peaks, the increasing speed of the MW is visible in the kymogram through the varying slope of the white

contour representing the MW traveling laterally over the upper surface of the vocal folds. Notice that, with the increased MW speed, the slope of the contour is moving towards the horizontal. Furthermore, it can also be observed that the increasing MW speed causes the width of the MW contour, which corresponds to the width of the mucosal wave peak, to increase and to become less distinguishable visually. This is due to the increase of the wavelength of the mucosal wave, which is directly proportional to the mucosal wave velocity. At very large speeds of the mucosal wave (i.e. for vertical phase differences smaller than ca. 30 ◦ in our model), its theoretical wavelength becomes so large that the mucosal wave on the upper surface becomes practically invisible (recall Fig.10). Fig. 12 shows the capability of the model to produce numerous vibration patterns that are analogous to those observed videokymographically in various phonation types [26,80–86], such as: a) soft voice, which can be obtained in the model by lowering the vibration amplitude, b) pressed voice, obtained by adjusting the prephonatory glottal width to be zero or negative, (c) breathy voice, obtained by enlarging prephonatory glottal width, or d) head/falsetto register, obtained by increasing fundamental frequency, slightly increasing prephonatory glottal width and decreasing vertical phase differences in comparison to the default values of the model (as indicated in Figs. 5 and 8, approximately representing modal male phonation at comfortable pitch and loudness). In contrast to the previous simulations with unlimited mucosal wave extent (Figs. 5,7–11), in Fig.12 the MW extent was reduced to half of the upper vocal fold surface by setting E = 0.5. 4. Discussion The developed kinematic model allows relating mathematically defined mucosal wave features to coronal shapes of the vocal folds during their vibratory behavior. The MW kinematics of the model approximates the trajectories measured by Berry et al [60] and Döllinger et al [10]. The model allows parametrically changing vocal

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

335

Fig. 11. Examples of synthetic kymograms with (a)  = 90◦ , (b)  = 70◦ and (c)  = 50◦ , simulating an increased velocity of the mucosal wave at 0.8, 1.03, and 1.14 m/s respectively. Notice the changes in diagonality of the mucosal wave contour, and in its breadth. Except of the , the VFM settings were identical to those in Figs. 8–10.

Fig. 12. Synthetic kymograms simulating the vocal fold vibratory patterns in various phonation types: a) soft voice, b) pressed voice, c) breathy voice, d) head/falsetto register. The model parameters were as follows: fo = 100 Hz (a–c), 400 Hz (d); Au = 0.02 cm (a), 0.07 cm (b–d); Al = Au (a–c), 0.7Au (d);  = 70◦ (a–c), 40◦ (d); E = 0.5 (a–d); prephonatory glottal width = 0.07 cm (a), -0.05 cm (b), 0.18 cm (c), 0.1 cm (d); M5 glottal convergence angle = 15◦ (a–d). The rest of the parameters were the same as in Figs. 5 and 8. Total time displayed in the kymograms, from top to bottom: 40 ms.

fold shape and vibration characteristics and to observe the influence of these parameters on the resulting kymogram. Simulation of the local illumination with a virtual laryngoscopic camera view enables capturing the vocal fold dynamics in synthetic kymograms by gray levels. The gray levels are particularly informative for distinguishing the glottis (black), the medial surface of the vocal folds (light to dark gray during the closing phase) and the upper vocal fold

Fig. 13. Synthetic kymograms with (a) sharp and (b) rounded lateral peaks obtained from the model comparable to the clinically observed VKG images showing the same respective features in (c) and (d). Total time displayed in the kymograms, from top to bottom: 40 ms.

surface (light gray to white) during vibration. This monochromatic light intensity also enables displaying the lateral propagation of the mucosal waves. The synthetic kymogram shows reasonable similarity to the laryngoscopically observed real VKG images (Fig. 13). The mucosal wave activity is biomechanically related to the vertical phase difference between the motions of the upper and lower margins. Large vertical phase differences are laryngoscopically recognizable in the VKG images as sharp lateral peaks. Recently, sharpness and roundedness of the lateral peaks in the VKG images have been found to be clinically helpful for diagnosis and treatment of voice disorders [59]. The developed model allows simulating various vertical phase differences and their influence on the sharpness of the lateral peaks in kymograms (Fig. 13a and b).

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337

336

MWs originate from the inferior surface of the vocal fold and travel medially upwards and laterally around the vocal fold [1–13]. They are clinically useful in assessing the pliability or stiffness of the vocal fold. However, the medial dynamics of the vocal fold is usually hidden from the endoscopic view and therefore only the superior surface vibrations are typically assessed and reported [35–39]. The inaccessibility of information on the medial dynamics severely restricts the scientific and clinical expertise on understanding the physical phenomena of the mucosal waves and associated features in voice diagnosis. Better comprehension and interpretation of mucosal wave dynamics from the medial surface may result in better diagnostic tools. The model allows simulating and displaying the medial dynamics of the vocal folds in the coronal section. This view offers the benefit of visualizing the mucosal waves as they propagate from the lower vocal fold margin upwards and further to the lateral end. Furthermore, it allows relating the vibratory dynamics of the vocal folds to the resulting kymograms. As such, the model can be used to better understand the different features of kymographic images, and consequently help to interpret features seen in kymographic images as obtained in clinical practice. Besides its educational value, the model can also be useful for scientific purposes, such as for testing and verifying algorithms for image analysis, glottal segmentation and parameterization of the vibratory contours of the vocal folds from videokymographic images [87–90]. Though MWs are a complex phenomenon, simplified equations and basic parameters are used in this simulation, neglecting the biomechanical tissue characteristics and their complex aerodynamic interactions. For the sake of brevity, only left-right symmetric and regular vocal fold vibrations are discussed here. A modified model, exhibiting asymmetric and abnormal vibratory patterns, is deferred to future studies. No doubt, the model can further be improved by implementing further knowledge from biomechanical experiments and detailed observations of mucosal wave dynamics under physiologic and pathologic conditions. Sophisticated experiments involving kymographic characterization of influential mucosal wave parameters and subsequently verifying their relationship through faithful reproduction of the synthetic kymograms should also be in the focus of future investigations. 5. Conclusion The developed kinematic model allows a) simulating mucosal waves traveling over the medial vocal fold surface upwards and then over the upper vocal fold surface laterally; b) parametrically changing the vocal fold shape and glottal configuration; c) displaying the dynamic changes of the vocal fold shape in the coronal plane, which are normally hidden in laryngoscopic observations; and d) creating synthetic kymograms with similar appearance to those obtained in clinical laryngologic practice, capturing the simulated vocal fold behavior in a single image. We find the model useful for deeper understanding of the appearance of mucosal waves on the vocal folds and their variability when changing the driving parameters. As such, this model can be used as an educational and scientific tool to better understand the physical phenomena of the mucosal waves and their appearance in kymographic images. This is desirable for advancing the diagnostic possibilities of kymographic imaging in laryngology. The model in the form of MATLAB scripts can be freely accessed at [91].

Acknowledgments The work has been supported by the Czech Science Foundation (GA CR) project no. GA16-01246S. S.P.K. programmed the model and made the simulations; J.G.S. provided the model concepts, verified the model functionality and supervised the work. Both authors wrote the article. The authors greatly acknowledge Prof. Jan Nauˇs, Department of Biophysics, Faculty of Science, Palacky´ University, Olomouc, Czech Republic, and Mr. Lukáˇs Vrajík, Zonky s.r.o., Brno, Czech Republic, for their helpful suggestions on the calculation of illumination of the synthetic kymogram. The authors also express their thanks to Dr. J. Vydrová for the clinical VKG images. S.P.K. also thanks Sri Sivasubramaniya Nadar (SSN) College of Engineering, Kalavakkam, India, for granting the leave of absence to the author during his sabbatical stay (2016-2018) at Palacky´ University in Olomouc, Czech Republic.

References [1] L. Réthi, Experimentelle Untersuchungen über den Schwingungstypus und den Mechanismus der Stimmbänder bei der Falsettstimme, in: Sitzungsberichte der mathematisch-naturwissenschaftlichen Classe der kaiserlichen Akademie der Wissenschaften, CV Band, Abtheilung III, Vienna, 1896, pp. 197–212. [2] D.W. Farnsworth, High-speed motion pictures of the human vocal cords, in: Bell Lab Record, 1940, pp. 203–208. [3] S. Smith, Remarks on the physiology of the vibrations of the vocal cords, Folia Phoniatr. (Basel) 6 (3) (1954) 166–178. [4] P. Moore, H. Von Leden, Dynamic variations of the vibratory pattern in the normal larynx, Folia Phoniatr. (Basel) 10 (4) (1958) 205–238. [5] I. Hiroto, Vibration of Vocal Cords: an Ultra High-Speed Cinematographic Study (film), Department of otolaryngology, Kurume University, Kurume, Japan, 1968. [6] T. Baer, Investigation of Phonation Using Excised Larynxes, Massachusetts Institute of Technology, Cambridge, Mass, 1975. [7] H. Matsushita, The vibratory mode of the vocal folds in the excised larynx, Folia Phoniatr. (Basel) 27 (1) (1975) 7–18. [8] S. Saito, et al., X-ray stroboscopy, in: K.N. Stevens, M. Hirano (Eds.), Vocal Fold Physiology, University of Tokyo Press, Tokyo, 1981, pp. 95–106. [9] M. Döllinger, N. Tayama, D.A. Berry, Empirical eigenfunctions and medial surface dynamics of a human vocal fold, Methods Inf. Med. 44 (3) (2005) 384–391. [10] M. Döllinger, D.A. Berry, S. Kniesburges, Dynamic vocal fold parameters with changing adduction in ex-vivo hemilarynx experiments, J. Acoust. Soc. Am. 139 (5) (2016) 2372–2385. [11] A. Boessenecker, et al., Mucosal wave properties of a human vocal fold, Acta Acust. United Acust. 93 (5) (2007) 815–823. [12] T. Vampola, J. Horacek, I. Klepacek, Computer simulation of mucosal waves on vibrating human vocal folds, Biocybern. Biomed. Eng. 36 (3) (2016) 451–465. [13] G.S. Berke, B.R. Gerratt, Laryngeal biomechanics: an overview of mucosal wave mechanics, J. Voice 7 (2) (1993) 123–128. [14] C.G. Garrett, J.R. Coleman, L. Reinisch, Comparative histology and vibration of the vocal folds: implications for experimental studies in microlaryngeal surgery, Laryngoscope 110 (5) (2000) 814–824. [15] M.F. Regner, M.J. Robitaille, J.J. Jiang, Interspecies comparison of mucosal wave properties using high-speed digital imaging, Laryngoscope 120 (6) (2010) 1188–1194. [16] C.T. Herbst, et al., Hemi-laryngeal setup for studying vocal fold vibration in three dimensions, J. Vis. Exp. 129 (2017) e55303. [17] C.T. Herbst, et al., Biophysics of vocal production in mammals, in: R.A. Suthers (Ed.), Vertebrate Sound Production and Acoustic Communication, Vol. 53, Springer International Publishing, Switzerland, 2016, pp. 159–189. [18] C.P.H. Elemans, et al., Universal mechanisms of sound production and control in birds and mammals, Nat. Commun. 6 (2015). [19] C.T. Herbst, et al., Complex vibratory patterns in an elephant larynx, J. Exp. Biol. 216 (21) (2013) 4054–4064. [20] I.R. Titze, Comments on the myoelastic-aerodynamic theory of Phonation, J. Speech Hear. Res. 23 (3) (1980) 495–510. [21] I.R. Titze, The physics of small-amplitude oscillation of the vocal folds, J. Acoust. Soc. Am. 83 (4) (1988) 1536–1552. [22] I.R. Titze, The Myoelastic Aerodynamic Theory of Phonation, National Center for Voice and Speech, Denver CO and Iowa City IA, 2006. [23] C.R. Krausert, et al., Mucosal wave measurement and visualization techniques, J. Voice 25 (4) (2011) 395–405. [24] H.S. Shaw, D.D. Deliyski, Mucosal wave: a normophonic study across visualization techniques, J. Voice 22 (1) (2008) 23–33. ˇ ˇ [25] J.G. Svec, F. Sram, Videokymographic examination of voice, in: E.P.M. Ma, E.M.L. Yiu (Eds.), Handbook of Voice Assessments, Plural Publishing, San Diego, CA, 2011, pp. 129–146.

S.P. Kumar, J.G. Sˇ vec / Biomedical Signal Processing and Control 49 (2019) 328–337 ˇ ˇ [26] J.G. Svec, F. Sram, H.K. Schutte, Videokymography in voice disorders: what to look for? Ann. Otol. Rhinol. Laryngol. 116 (3) (2007) 172–180. [27] M. Kaneko, et al., Optimal duration for voice rest after vocal fold surgery: randomized controlled clinical study, J. Voice 31 (1) (2017) 97–103. [28] A. Yamauchi, et al., Characterization of vocal fold vibration in sulcus vocalis using high-speed digital imaging, J. Speech Lang. Hear. Res. 60 (1) (2017) 24–37. [29] L. Li, et al., Quantitative study for the surface dehydration of vocal folds based on high-speed imaging, J. Voice 29 (4) (2015) 403–409. [30] D. Voigt, et al., Objective detection and quantification of mucosal wave propagation, J. Acoust. Soc. Am. 128 (5) (2010) EL347–EL353. [31] B.H. Story, An overview of the physiology, physics and modeling of the sound source for vowels, Acoust Sci Tech 23 (2002) 195–206. [32] E. Yumoto, H. Kurokawa, H. Okamura, Vocal fold vibration of the canine larynx: observation from an infraglottic view, J. Voice 5 (4) (1991) 299–303. [33] I.R. Titze, J.J. Jiang, T.-Y. Hsiao, Measurement of mucosal wave propagation and vertical phase difference in vocal fold vibration, Ann. Otol. Rhinol. Laryngol. 102 (1) (1993) 58–63. [34] R. McGowan, An analogy between the mucosal waves of the vocal folds and wind waves on water, Haskins Lab. Status Rep. Speech Res. 101 (102) (1990) 243–249. [35] M. Hirano, D.M. Bless, Videostroboscopic Examination of the Larynx, Singular Publishing Group, San Diego, California, 1993. [36] G. Cornut, M. Bouchayer, Assessing dysphonia: the role of video stroboscopy, in: An Interactive Video Textbook, Kay Elemetrics Corp., Lincoln Park, NJ, 2004. [37] M.M. Abaza, Abnormalities of the mucosal wave, in: R.J. Leonard (Ed.), Indirect Laryngoscopy to High-Speed Digital Imaging, Thieme, New York, 2010. [38] B.J. Poburka, R.R. Patel, D.M. Bless, Voice-vibratory assessment with laryngeal imaging (VALI) form: reliability of rating stroboscopy and high-speed videoendoscopy, J. Voice 31 (4) (2017) 513.e1–513.e14. [39] R.R. Patel, et al., Recommended minimum protocols for instrumental assessment of voice: American Speech-Language Hearing Association Committee on Instrumental Voice assessment protocols, Am. J. Speech. Pathol. 27 (3) (2018) 887–905. [40] D.D. Mehta, R.E. Hillman, Current role of stroboscopy in laryngeal imaging, Curr. Opin. Otolaryngol. Head Neck Surg. 20 (6) (2012) 429–436. [41] P. Woo, Stroboscopy, Plural Publishing, San Diego, CA 92123, 2009. [42] K.A. Kendall, R.J. Leonard, Laryngeal Evaluation: Indirect Laryngoscopy to High-speed Digital Imaging, Thieme, New York, 2010. [43] D. Bless, R. Patel, N. Connor, Laryngeal imaging: stroboscopy, high-speed digital imaging, and kymography, in: M. Fried, A. Ferlito (Eds.), The Larynx, Plural Publishing, San Diego, CA, Oxford, and Brisbane, 2009, pp. 182–191. [44] S. Hertegard, What have we learned about laryngeal physiology from high-speed digital videoendoscopy? Curr. Opin. Otolaryngol. Head Neck Surg. 13 (3) (2005) 152–156. [45] H. Hirose, High-speed digital imaging of vocal fold vibration, Acta Otolaryngol. (1988). [46] K.A. Kendall, High-speed digital imaging of the larynx: recent advances, Curr. Opin. Otolaryngol. Head Neck Surg. 20 (6) (2012) 466–471. [47] D.D. Deliyski, et al., Experimental investigation on minimum frame rate requirements of high-speed videoendoscopy for clinical voice assessment, Biomed. Signal Process. Control 17 (2015) 21–28. [48] T. Ikuma, M. Kunduk, A.J. McWhorter, Preprocessing techniques for high-speed videoendoscopy analysis, J. Voice 27 (4) (2013) 500–505. [49] J. Lohscheller, U. Eysholdt, Phonovibrogram visualization of entire vocal fold dynamics, Laryngoscope 118 (4) (2008) 753–758. [50] Q. Qiu, H.K. Schutte, A new generation videokymography for routine clinical vocal fold examination, Laryngoscope 116 (10) (2006) 1824–1828. ˇ [51] J.G. Svec, H.K. Schutte, Videokymography: high-speed line scanning of vocal fold vibration, J. Voice 10 (2) (1996) 201–205. ˇ ˇ H.K. Schutte, Videokymography, in: M. Fried, A. Ferlito [52] J.G. Svec, F. Sram, (Eds.), The Larynx, Plural Publishing, San Diego, CA, 2009, pp. 253–271. ˇ [53] J.G. Svec, H.K. Schutte, Kymographic imaging of laryngeal vibrations, Curr. Opin. Otolaryngol. Head Neck Surg. 20 (6) (2012) 458–465. [54] G. Andrade-Miranda, N.H. Bernardoni, J.I. Godino-Llorente, Synthesizing the motion of the vocal folds using optical flow based techniques, Biomed. Signal Process. Control 34 (2017) 25–35. [55] J.J. Jiang, et al., Quantitative study of mucosal wave via videokymography in canine larynges, Laryngoscope 110 (9) (2000) 1567–1573. [56] A. Yamauchi, et al., Quantitative analysis of vocal fold vibration in vocal fold paralysis with the use of high-speed digital imaging, J. Voice 30 (6) (2015) 766.e13–766.e22. [57] D.M. Bless, M. Hirano, R.J. Feder, Videostroboscopic evaluation of the larynx, Ear Nose Throat J. 66 (7) (1987) 289–296. [58] J. Sundberg, C. Högset, Voice source differences between falsetto and modal registers in counter tenors, tenors and baritones, Logoped. Phoniatr. Vocol. 26 (1) (2001) 26–36. [59] K.V. Phadke, et al., Evaluation of clinical value of videokymography for diagnosis and treatment of voice disorders, Eur. Arch. Oto-rhino-laryngology 274 (11) (2017) 3941–3949. [60] D.A. Berry, D.W. Montequin, N. Tayama, High-speed digital imaging of the medial surface of the vocal folds, J. Acoust. Soc. Am. 110 (5) (2001) 2539–2547. [61] M. Döllinger, D.A. Berry, Visualization and quantification of the medial surface dynamics of an excised human vocal fold during phonation, J. Voice 20 (3) (2006) 401–413.

337

[62] T. Baer, in: K.N. Stevens, M. Hirano (Eds.), Observation of Vocal Fold Vibration: Measurements of Excised Larynges, University of Tokyo Press, Tokyo, 1981. [63] M. Hirano, et al., Vibratory behavior of the vocal folds in a case with a unilateral polyp, Folia Phoniatr. (Basel) 33 (5) (1981) 275–284. [64] S.H. Sloan, et al., Determination of vocal fold mucosal wave velocity in an in vivo canine model, Laryngoscope 103 (9) (1993) 947–953. [65] S. Nasri, J.A. Sercarz, G.S. Berke, Noninvasive measurement of traveling-wave velocity in the canine larynx, Ann. Otol. Rhinol. Laryngol. 103 (10) (1994) 758–766. [66] D.G. Hanson, et al., Clinical measurement of mucosal wave velocity using simultaneous photoglottography and laryngostroboscopy, Ann. Otol. Rhinol. Laryngol. 104 (5) (1995) 340–349. [67] Y.W. Shau, et al., Noninvasive assessment of vocal fold mucosal wave velocity using color Doppler imaging, Ultrasound Med. Biol. 27 (11) (2001) 1451–1460. [68] B. Jing, et al., Visualizing the mechanical wave of vocal fold tissue during phonation using electroglottogram-triggered ultrasonography, J. Acoust. Soc. Am. 143 (5) (2018) El425–El429. [69] J. Schoentgen, P. Aichinger, Synthetic kymograms and glottal area waveforms in simulated non-neutral phonation, in: C. Manfredi (Ed.), MAVEBA 2015: Models and Analysis of Vocal Emissions for Biomedical Applications, 9th International Workshop (2015) 75–78. ˇ ˇ [70] J. Horáˇcek, J.G. Svec, P. Sidlof, Numerical simulation of videokymographic images of self-oscillating vocal folds, in: J.I. Godino-Llorente, P. Gomez Vilda, R. Fraile (Eds.), AVFA 09, 3rd Advanced Voice Function Assessment International Workshop (2009) 13–16. ˇ [71] P. Svancara, et al., Numerical simulation of videokymographic images from the results of the finite element model, 20th International Conference Engineering Mechanics (2014) 640–643. [72] S. Li, et al., The effect of entrance radii on intraglottal pressure distributions in the divergent glottis, J. Acoust. Soc. Am. 131 (2) (2012) 1371–1377. [73] R.C. Scherer, et al., Intraglottal pressure profiles for a symmetric and oblique glottis with a divergence angle of 10 degrees, J. Acoust. Soc. Am. 109 (4) (2001) 1616–1630. [74] J.D. Foley, et al., Computer Graphics: Principles and Practice, second ed., Addison-Wesley Publishing Company, 1996. [75] S. Marschner, P. Shirley, Fundamentals of Computer Graphics, fourth ed., CRC Press, 2015. [76] E. Yumoto, Y. Kadota, T. Mori, Vocal fold vibration viewed from the tracheal side in living human beings, Otolaryngol. Head. Neck Surg. 115 (4) (1996) 329–334. [77] Q. Qiu, H.K. Schutte, Real-time kymographic imaging for visualizing human vocal-fold vibratory function, Rev. Sci. Instrum. 78 (2) (2007) 1–6. [78] M. Pharr, G. Humphreys, Physically Based Rendering, second ed., Morgan Kaufmann, 2010. [79] M. Oren, S.K. Nayar, Generalization of Lambert’s reflectance model, SIGGRAPH ’94 Proceedings of the 21st Annual Conference on Computer Graphics and Interactive Techniques (1994) 239–246. ˇ [80] J.G. Svec, On Vibration Properties of Human Vocal Folds: Voice Registers, Bifurcations, Resonance Characteristics, Development and Application of Videokymography, University of Groningen, Groningen, the Netherlands, 2000. ˇ ˇ [81] M. Friˇc, F. Sram, J.G. Svec, Voice registers, vocal folds vibration patterns and their presentation in videokymography, Proceedings ACOUSTICS High Tatras 06. 33rd International Acoustical Conference - EAA Symposium, Strbske Pleso, Slovakia (2006). ˇ [82] J.G. Svec, J. Sundberg, S. Hertegard, Three registers in an untrained female singer analyzed by videokymography, strobolaryngoscopy and sound spectrography, J. Acoust. Soc. Am. 123 (1) (2008) 347–353. ˇ [83] C.T. Herbst, S. Ternström, J.G. Svec, Investigation of four distinct glottal configurations in classical singing-a pilot study, J. Acoust. Soc. Am. 125 (3) (2009) EL104–EL109. [84] M. Echternach, et al., High-speed imaging and electroglottography measurements of the open quotient in untrained male voices’ register transitions, J. Voice 24 (6) (2010) 644–650. [85] P.H. Dejonckere, et al., Automated tracking of quantitative parameters from single line scanning of vocal folds: a case study of the’ messa di voce’ exercise, Logoped. Phoniatr. Vocol. 40 (1) (2015) 44–54. [86] H. Yokonishi, et al., Relationship of various open quotients with acoustic property, phonation types, fundamental frequency, and intensity, J. Voice 30 (2) (2016) 145–157. [87] C. Manfredi, et al., Videokymographic image processing: objective parameters and user-friendly interface, Biomed. Signal Process. Control 7 (2) (2012) 192–201. [88] C. Piazza, et al., Quantitative analysis of videokymography in normal and pathological vocal folds: a preliminary study, Eur. Arch. Oto-rhino-laryngology 269 (1) (2012) 207–212. [89] Q. Qiu, et al., An automatic method to quantify the vibration properties of human vocal folds via videokymography, Folia Phoniatr. Logop. 55 (3) (2003) 128–136. [90] S.P. Kumar, et al., Visual and automatic evaluation of vocal fold mucosal waves through sharpness of lateral peaks in high-speed videokymographic images, J. Voice (2018), http://dx.doi.org/10.1016/j.jvoice.2018.08.022 (Early Online). [91] S.P. Kumar, J.G. Svec, Simulation of Vocal Fold Mucosal Waves & Synthetic Kymograms, MATLAB (2018), Central File Exchange Retrieved December 18, from https://in.mathworks.com/matlabcentral/fileexchange/69753simulation-of-vocal-fold-mucosal-waves-synthetic-kymograms.