## 4 Classification
### Feature Extraction  
Features were extracted at the utterance level. The utterances were isolated from the raw data using a Vocal Activity Detection (VAD) algorithm, which was based on the short-time energy and zero crossing rate methods [34]. A Hamming window was used to obtain NSA spectra. The NSA spectral harmonic quotients were calculated using,  
\[H_{i}=\frac{A_{i}}{A_{1}+A_{2}+A_{3}+A_{4}+A_{5}}\quad(i=1,2,3,4), \tag{1}\]where \(A_{i}(i=1,2,3,4,5)\) are the amplitudes of the first five harmonics in the NSA spectra. The sum of \(H_{i}(i=1,2,3,4,5)\) is unity according to Equation (1). The first four components, \(H_{i}(i=1,2,3,4)\), were selected as features, with \(H_{5}\) excluded because \(H_{5}\) was linearly related to \(H_{i}(i=1,2,3,4)\) and thus was redundant. The harmonics were used to approximate the envelopes of NSA spectra in a limited frequency bandwidth from 50 Hz-1500 Hz. The SE described the complexity of a signal and was defined as,  
\[\mathrm{SE}=-\sum\limits_{i=1}^{n}p_{i}\log p_{i}, \tag{2}\]  
where \(p_{i}\) is the normalized spectral density point, which is from 0-3000 Hz. For a simple signal, e.g., an ideal sinusoidal function, most of the spectral energy was concentrated within a narrow bandwidth, resulting in a small SE. On the contrary, a white noise has a broadly distributed spectrum, and the resulting SE was large. To include voice stability information in the feature set, the jitter and the shimmer were extracted from the time-domain waveforms. Since the lengths of the utterance samples were not identical, jitter (\(J_{r}\)) and shimmer (\(S_{r}\)) were calculated as relative percentages using,  
\[J_{r}=\frac{\frac{\frac{1}{N-1}\sum\limits_{i=1}^{N-1}|T_{i}-T_{i+1}|}{\frac{1 }{N}\sum\limits_{i=1}^{N}T_{i}}}{\frac{1}{N}\sum\limits_{i=1}^{N}T_{i}}, \tag{3}\]  
and:  
\[S_{r}=\frac{\frac{1}{N-1}\sum\limits_{i=1}^{N-1}|A_{i}-A_{i+1}|}{\frac{1}{N} \sum\limits_{i=1}^{N}A_{i}}. \tag{4}\]  
where \(T_{i}\) and \(A_{i}(i=1,2,...N)\) were the period of each vocal cycle and the peak magnitude in each vocal cycle, respectively. \(N\) was the number of sampling points for each voice sample. The feature vector had the form of \([H_{1},H_{2},H_{3},H_{4},\mathrm{SE},J_{r},S_{r}]\). One representative NSA waveform and its spectrum are shown in Figure 4. The feature statistics for 31 participants are shown in Figure 5a,b. Figure 5a shows the shapes of the spectral envelopes for each voice type. The breathy voice features a prominent \(H_{1}\) component, which accounts for over 60% of the total energy. The pressed voice features a prominent \(H_{2}\). The modal voice harmonics decreased monotonically with frequency from \(H_{1}\)-\(H_{4}\), but less rapidly than those of the breathy voice. The SE of pressed voice was notably different from that of the other voice types, as shown in Figure 5b. The average SE value of pressed voice was greater, which indicated that the pressed voice spectrum had a greater bandwidth than the other voice types. A comparison between the shimmer values across voice types showed that the modal voice had the smallest average value and standard deviation, which indicated that the modal voice has the most stable loudness. Figure 5a,b shows the notable feature range overlap across voice types, which caused difficulties in voice type classification. The classification results obtained using the Linear Discriminant (LD) method based on single features are shown in Table 1. The overall accuracies for every single feature classification did not exceed 65%. \(H_{1}\), SE and shimmer contributed the most to the classification accuracies for pressed (from TP to PP: 68%), breathy (from TB to PB: 93%) and modal (from TM to PM: 92%) types, according to the confusion matrix in Table 1. The modal and breathy voice had similar SE distributions. The LD method could not recognize modal voice and misclassified all modal voice samples into breathy or pressed types.