{'fraction_non_alphanumeric': 0.02947968638631504, 'fraction_numerical': 0.01377049180327869, 'mean_word_length': 4.830452127659575, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Using instantaneous frequency and aperiodicity detection to estimate F0 for high-quality speech synthesis

This paper introduces a general and flexible framework for F0 and aperiodicity (additive non periodic component) analysis, specifically intended for high-quality speech synthesis and modification applications. The proposed framework consists of three subsystems: instantaneous frequency estimator and initial aperiodicity detector, F0 trajectory tracker, and F0 refinement and aperiodicity extractor. A preliminary implementation of the proposed framework substantially outperformed (by a factor of 10 in terms of RMS F0 estimation error) existing F0 extractors in tracking ability of temporally varying F0 trajectories. The front end aperiodicity detector consists of a complex-valued wavelet analysis filter with a highly selective temporal and spectral envelope. This front end aperiodicity detector uses a new measure that quantifies the deviation from periodicity. The measure is less sensitive to slow FM and AM and closely correlates with the signal to noise ratio.


Introduction
This paper describes a new F0 tracker for rapidly changing F0 trajectories with aperiodicity, which represents additive nonperiodic components.
In high-quality speech synthesis and modification applications [1][2][3], surpassing 4.2 on the 5 point MOS score, glitches in aperiodicity handling and the failure to follow rapidly changing fundamental frequencies (F0) are harmful to processed speech quality. Introducing a generative model of F0 trajectory (for example [4]) to F0 estimation provides well behaved and parametric representation. However, the estimated F0 trajectories are still not good enough for highquality speech synthesis. The actual excitation signal of speech, glottal flow, contains several sources of fluctuations [5] and consequently, the observed F0 trajectories are different from the trajectories produced by those models. To attain highly natural synthetic speech it is important to retain these fine temporal variation in F0 trajectories [6,7]. Although many F0 extractors have been proposed [8][9][10][11][12], in practice, parameter tuning and/or manual error correction is often necessary. In addition, their performance when extracting such fine temporal variations has not been investigated explicitly. That is the goal of this paper.
This paper is organized as follows. Section 2 discusses the motivation and target for designing a new F0 observer, based on a review on existing issues. It also defines aperiodicity, which is relevant for speech analysis and synthesis. Section 2.2 presents objective measures used in this paper. Based on these, section 3 introduces a general scalable architecture for F0 observer. It consists of three subsystems: front end aperiodicity detectors, the best trajectory finder, and F0 initial estimate and refinement subsystem with aperiodicity extractor. Sub-sections 3.1 and 3.3 introduce the front end and the refinement subsystems, respectively. In section 4, these subsystems are evaluated using artificial test signals. Section 5 discusses remaining issues. Example analysis results using actual speech samples and mathematical details are given in appendices.

Background
Speech synthesis requires dependable F0 values whenever producing voiced sounds. However, even for copy-synthesizing from actual speech samples, where the targets are known, this is not always easy, since voiced sounds are not purely periodic and defining F0 values for such signals is not a trivial issue. Figure 1 shows a beginning of a sentence from our speech corpus. From 0.2 s to 0.52 s, the speech signal is voiced. However, due to irregularities in glottal vibrations, defining the F0 is difficult. The lower plot shows the F0 tracks by YIN [10] and SWIPE ′ [12] to illustrate the issues. It is difficult to evaluate the relevance of these tracking results. Yet these two state of the art systems do not produce consistent results. The fact that voicing without vocal fold contact is not rare [13,14] prevents using EGG (electroglottograph) for the source of ground truth. Using the extracted trajectory and comparing the synthesized speech and the original speech is a reasonable test but it is very demanding on human resource and time to obtain reliable results.
An alternative approach for evaluating F0 extractors is to use an objectively defined artificial test signal. The ideal candidate is a speech signal, where the ground truth is available and provides wide divergence and variability. Instead, this article uses the excitation source signal defined by the L-F (Liljencrants-Fant) model [15]. The L-F model represents the time derivative of the glottal flow using a set of equations with four parameters. However, directly digitizing the L-F model, which is defined in the continuous time domain, introduces spurious components due to aliasing. To alleviate this aliasing problem this paper uses a closed-form representation of the antialiased L-F model defined in the continuous time domain [16].
Since the model is defined in the continuous time domain, it is easy to generate a signal using a given F0 trajectory that will be the ground truth used in this paper. 1 [19][20][21] Figure 2 shows an example of F0 tracking using a sinusoidally frequency modulated F0 trajectory as the test signal. This test signal has a vibrato of 16 Hz, which is large compared to the normal human voice, but demonstrates the problems due to random, cycle-by-cycle variations in the F0. The tested F0 extractors are YIN [10], SWIPE ′ [12], NDF [11], DIO [22] and the proposed method, which is described in Section 3. The trajectories obtained by YIN and SWIPE ′ are strongly distorted and attenuated, perhaps because the F0 is changing faster than these models allow. When these distorted trajectories are used to generate the excitation source for copy-synthesis, the output is perceived differently. This is because the distortion adds fast-changing modulation components that are not in the original signal. The effects of these spurious components are made worse because humans are far more sensitive to fast frequency modulations than amplitude modulations [23,24].
Voiced sounds are usually considered as periodic, and to first approximation the glottal pulses do occur at regular intervals. But due to prosodic needs the F0 of a voice is constantly 1 In an open-source implementation [17] of the anti-aliased L-F model [16], the model parameters can be controlled each glottal cycle independently to simulate the details of vocal fold behaviour [18]. It can be combined with a time varying lattice filter to simulate the dynamic speech production process, which modulates observed F0 through interaction between harmonic component and the group delay associated with resonances (formant trajectories). But these detailed simulations are for further study. changing, sometimes a simple glide as in the rise of F0 in a question, and sometimes in a regular fashion, as with vibrato. And, sometimes F0 varies in a more complex patterns, such as in tonal languages, where the F0 trajectory conveys linguistic information. On top of these intended changes in F0, there are modulations due to physiological aspects of voice production. The stochastic nature of neural pulses which drive the muscles of the vocal organ is a strong noise source and the critical conditions that produce vocal fold oscillation introduce bi-stable or chaotic vocal fold vibration, especially during voice onset and offset. Age related change and physical body status also affects the stability of vibration [5]. All these deviations from pure periodicity play important roles in speech communication and make speech a much richer media than text [25]. It is important to properly analyse and replicate these deviations from periodicity in high-quality speech synthesis and modification applications. Accurately estimating aperiodicity is still a very challenging problem. Tracking errors introduces spurious components [26,27] and they add to the original random component. These are the reasons why F0 tracking distortions as shown in Fig. 2 are harmful for high-quality speech synthesis. Two issues have to be properly solved : accurate estimate and tracking of changing F0 trajectory and accurate estimate of random components based on the accurate estimate of F0 trajectory.
These issues motivate us to develop a framework that provides a calibrated procedure to describe the amount of aperiodicity and to track F0. The primary analysis target is high quality speech corpus recorded in a quiet and acoustically controlled environment using high-fidelity microphones. The aim here is to provide accurate, certified metadata, in this case, F0 value and an index that represents the accuracy of the estimated F0 as well as a measure that represents the amount of aperiodicity. Processing speed is not the first priority of the framework described here. Note that these metadata depend only on the data in the analysis frame, because there is no reliable model yet for the dynamic behaviour of F0 and aperiodic compo-nent. Using models of dynamic F0 behaviour such as Fujisaki's model [28], or F0 continuity constraint, may introduce biases due to model mismatch. Frame-based F0 with aperiodicity information, which the proposed system produces, will help to establish certifiably accurate models of the statistical/dynamic behaviour.

What is aperiodicity?
For speech synthesis applications, amplitude and F0 are controllable parameters of the excitation source. However, only replicating amplitude and F0 precisely to the original speech yields poor quality synthetic sounds. An important attribute of excitation is missing. This missing attribute is aperiodicity. 2 In this paper, attributes that can be represented by amplitude and F0 modulation are not included in the definition of aperiodicity. What is left after removing periodic component defines "aperiodicity" in this paper. It turns out that our system's F0 estimation error is well correlated with the system's estimate of aperiodicity, described below.

Measures for objective evaluation
F0 extractors have been evaluated based on error-rate related measures; such as Gross Pitch Error (GPE), Voicing Detection Error (VDE) [29] and Pitch Tracking Error (PTE) [30]. Attaining high performance in these measures is a prerequisite for good F0 extractors. In this paper, we focus on F0 tracking fidelity, because the proposed method does not make voiced/unvoiced decision. Instead, this F0 tracker outputs a measure of aperiodicity, which closely correlates with the standard deviation of the relative F0 estimation error from the true value. This aperiodicity detector also is an informative source of the type of excitation. The voiced/unvoiced decision is left to the application, which can use the output of the proposed method to make this decision.

Architecture and subsystems
The proposed framework, YANGSAF (Yet ANother Glottal Source Analysis Framework), computes the instantaneous F0 using three steps: estimate, track, and refine. The estimation step calculates three features of the input signal over a number of bandpass channels. The maximum from the estimate stage is then tracked to produce a local estimate of the F0. Finally, an optional refinement stage combines temporal and harmonic information to produce a more accurate estimate of F0.

Estimation
The first stage of the YANGSAF algorithm analyses the signal with a number of bandpass channels, and then estimates three values for each channel as a function of time. These values are 1) the local instantaneous frequency, 2) a measure called aperiodicity that represents the amount of variability in the channel's frequency estimate, and 3) a probabilistic estimate that the channel contains a good representation of the F0. These signals are described in the subsections that follow and are used in the tracking stage described by Section 4.2. Figure 3 shows a diagram of the estimating detector in each channel. The front end breaks the input into a number of spectral channels using a bank of bandpass filters, each centered at fc. 3 The center frequencies cover the possible F0 range, with a fixed separation on the logarithmic frequency axis. The current implementation covers 400 Hz to 1000 Hz using 12 channels and detectors in each octave.
The instantaneous frequency estimate needs both the complex-valued signal and its derivative. These values are calculated starting with bandpass filter h(τ, fc) and its derivative h d (τ, fc) shown in Fig. 3 and described in Appendix A. Each bandpass filter has linear phase, is a zero-delay FIR filter, has a complex-valued response, and passes only the positive frequency components. Figure 4 shows an example of these three estimated signals for a sequence of vowels.

Instantaneous Frequency
The instantaneous frequency of the signal contained within each channel is calculated using Flanagsn's approach, which is based on the logarithm of a complex signal x(t) and its derivative. An AM/FM modulated signal is represented in polar form x(t) = r(t)e jθ(t) . The instantaneous (angular) frequency ωi(t) is defined as the derivative of the phase component θ(t), namely ωi(t) = dθ(t) dt . The instantaneous frequency can be derived by starting with the logarithm of the component phase and using a bit of algebra: where ℜ[x] and ℑ[x] represents the real and the imaginary part of x, respectively. The derivation of this expression is contained in Appendix A.

Aperiodicity
We also wish to calculate a measure of the aperiodicity of the signal in each channel, which will be used as a measure of the reliability of the instantaneous frequency measurement. constant sinusoid, the aperiodicity is zero, and the aperiodicity grows as the signal varies (wiggles) more within the bandpass channel. The basic idea of the periodicity detector is to calculate the amount of energy in the band-passed signal that is not the primary sinusoid. The primary sinusoidal component will have the largest energy, and when the complex signal is normalized to have unit magnitude, refiltered, and then renormalized, the primary sinusoid will still have unit magnitude. The other components will be filtered with a non-unit gain, since the filter is not an ideal brick-wall filter, and their amplitude will change. Subtracting the original and the twice-filtered and normalized response gives an estimate of the aperiodicity. Note this estimate is done without explicitly identifying the primary sinusoid and its frequency. When a signal x whose fundamental frequency is equal to fc is filtered, only the fundamental component, a complexvalued, slowly time-varying signal, is passed (appears in y1) and is normalized to become y ′ 1 . Then, by using the same filter, filtering signal y ′ 1 again, and normalizing the overall amplitude using the absolute value of the complex valued-signal, the twice filtered (and amplitude normalized) signal y ′ 2 is obtained. Subtracting this twice filtered and amplitude normalized signal y ′ 2 from the amplitude normalized first filter output y ′ 1 , yields a residual signal r. Since the signal y ′ 1 is normalized, the power of the residual represents the relative level of the other component(s).
The difference between y ′ 1 and y ′ 2 corresponds to spectral components in the channel that are not the primary sinusoid. Calculating the energy in this signal (a k ), and smoothing it gives a ks which is this system's measure of harmonic aperiodicity. Appendix B describes the relation between the SNR of the original signal and the residual aperiodicity power using equations and examples.
Placing bandpass filters having the same shape on the logarithmic frequency axis yields the detector to output higher aperiodicity value, when fc is located at harmonic frequencies other than the fundamental. This is similar to the concept "fundamentalness," which is explained in Fig. 11 of reference [31]. Appendix shows relation between filter shape examples and harmonic components.
The instantaneous frequency calculation and the aperiodicity calculation yield values at the audio sampling rate. These audio sampling rate time series are down-sampled for later processing. In this work the down-sampling is accomplished by extracting the nearest time samples from each time series, providing two sequences of instantaneous frequency and aperiodicity measure values at the frame rate (i.e. 200 Hz).

Probability
The fundamental component in the original signal is dominant in a number of output channels because there is little else for filters centered at frequencies lower than the second harmonic can respond. Thus a number of channels will respond in the same way to the fundamental component, as seen by the blueish blob around 100Hz in the second panel of Figure 4. All channels inside this blob have information about the fundamental component, but with different reliabilities.
Given a number of (distinct) estimates of the true F0, all from different channels, a probability map indicates which channel will have the best estimate. To create this probability map, all the instantaneous frequency and aperiodicity estimates are converted into Gaussian probability masses centered at various instantaneous frequency estimates. The output of the channel's aperiodicity estimate (a ks , a measure of smoothed en-ergy) is converted into a variance σ 2 k by scaling. The scaling coefficient was empirically determined by a set of simulations. On a log-frequency scale ν, this gives a number of (independent) estimates of the instantaneous frequency, each modelled as a Gaussian mass centered at log(f k ), and with a variance of σ 2 k . Summing all these yields a probability density function pG(ν) represented as a Gaussian mixture. For each channel, integrating this distribution provides an observation probability Pr[k] that channel k should see the fundamental component in its nominal pass band [fL(k), fH (k)] is where K represents the number of filters per octave. This integrates the instantaneous frequency probability distributions between the frequency limits of filter k to arrive at an estimate of how reasonable it is for channel k to provide an estimate of the F0. An example of this result is shown in the bottom of Figure4.

Tracking
Given the three instantaneous maps (as a function of frame time and spectral channel) computed in Section 3.1, an initial estimate of the single best F0 at each frame is calculated by finding the channel with the highest probability. This is done in four steps: estimate the pitch range for this utterance, smooth the probability map, find the highest probability F0, and then refine the F0 estimate. The result is a smooth estimate of the true F0 based on the instantaneous frequency calculated in each channel.
First, the F0 search range is estimated by a weighted average of the instantaneous frequencies seen in the utterance. The temporal weighting is calculated from the energy in the original signal, after filtering it between 40-1000Hz, which is the prospective pitch range. Then each frame of the instantaneous frequency map is weighted and combined to form an overall instantaneous frequency histogram. By weighting by the signal's amplitude at each point in time, the high-energy portion of the utterance (vowels) are treated with more importance.
The median of this instantaneous frequency distribution (marginal distribution) defines the center point of the F0 search range. The tracker looks for peaks in the probability distribution within 1.2 octaves above this center point, and 1.3 octaves below, a total of a 2.5 octave range.
Second, in order to better estimate the F0 at the start and end of voicing the probability map computed in Section 3.1.3 is smoothed in time using a 45ms Hanning window with amplitude weighting. Smoothing is done before tracking so that we extend the F0 estimates at the start and end of voiced segments. For example, at the onset of voicing, the probability at F0 is not high, because the signal level is low and the SNR is low. Smoothing using amplitude weighting increases the probability at F0, because at frames after the onset the level grows and consequently the SNR become higher. In other words, the probability distribution of the onset frames become more like the probability distribution of later frames. This way smoothing reduces tracking error at the beginning of voicing. The same thing happens at the voice offset.
Thirdly, given the F0 range and the smoothed probability map the best channel across time can be tracked. For a range of channels that are within the 2.5 octave range defined for the entire utterance, and 0.7 octaves of the last frames best channel, the channel with the highest smoothed probability is chosen.
Finally, this channel selection is further refined by returning to the original probability map computed in Section 3.1.3 and choosing the channel with the highest probability closest to that bin chosen from the smoothed estimate. The following provides the initial F0 estimate fOI .
where the best weights bm are calculated from σ 2 m in V[k].

Refinement of the initial estimate
The third stage further improves this F0 estimate by adding two refinements. First, and most importantly, the higher harmonics of an F0 estimate can refine the estimate. Secondly, adaptive time warping of the original signal, combined with further refinement using higher harmonics of the warped signal, reduces the amount of F0 trajectory deviation for better analyses. The first procedure uses harmonic frequencies and their variance. Each harmonic component, from first to m-th, has corresponding aperiodicity detector. Each bandpass filter of the detector has the same shape on the linear frequency axis and does not cover neighbouring harmonic components. Each detector yields instantaneous frequency f k and its aperiodicity a k , where k represents the harmonic number. These values are converted to F0 estimate f k /k and its variance σ 2 k . The weighted average m k=1 b k f k /k provides the refined F0 estimate. Variance values {σ 2 k } m k=1 are used to calculate the best mixing weights {b k } m k=1 (Appendix D). However, this refinement does not properly make use of higher harmonic information when the F0 trajectory is rapidly changing. This is because a rapid movement of higher frequencies generates strong side-band components and they smear the analysed harmonic structure [26,27,32].
Thus, the second procedure uses F0 adaptive time axis warping to alleviate this problem. Stretching the time axis, proportional to an instantaneous F0 value makes the observed F0 value constant [26,27,32] and keeps the harmonic structure intact. Then, placing aperiodicity detectors on harmonic frequencies, from first to m-th, the weighted average of F0 information yields the F0 estimate on the warped time axis. Converting this estimate value to the value on the original time axis provides the further improved F0 estimate.
These two procedures are applied serially as well as recursively. Let Hm represent the operation of harmonic based refinement using the first through m-th harmonic components and Tm represent the operation of F0 adaptive time warping-based refinement using the first through m-th harmonic components. Let PX[x; Θ] represent the function of initial estimate F0 where x represents the input signal and Θ represents a set of the associated design parameters for analysis. The following equations describes the configurations of the two trackers tested: where T • H represents the composite function of the functions T and H. Finally, by placing aperiodicity detectors on all harmonic frequencies in the warped time axis, estimated SNR around each harmonic component provides the excitation source information for speech synthesis. Because any F0 trajectories on this warped time axis are constant in time, aperiodicity values which detectors output are consistent with the aperiodicity definition of this paper.

Evaluation using test signals
This paper uses two measures of performance. Most importantly, the standard deviation of the relative error tells us the total distortion of the estimated F0 trajectory from the ground truth. The second performance measure is the frequencymodulation amplitude transfer function (FMTF), which expresses how well a F0 tracker follows fast F0 modulations. The test signal uses sinusoidal modulation on the logarithmic frequency axis, since F0 dynamics is better described on the logarithmic frequency axis [28]. Consequently, both FMTF and distortion evaluation measures use logarithmic frequency to calculate their value.
The proposed algorithms are implemented using MATLAB and tested using synthetic signals. Only representative results are described below. In the following tests, the test signals were generated using the aliasing-free L-F model [16]. 4 The sampling frequency fs was 22050 Hz and the "modal" voice quality parameters [33] for the L-F model were used in the following examples.
We test this new F0 tracker in two different ways: additive noise and FM modulation.

Additive noise
Firstly, the quality of the F0 estimate in the face of additive white noise was tested using the configuration given by Eq. 8 (H10 • H3). The F0 extractor for the initial estimate (Section 3.2) (PX[x; Θ]) was tested to clarify the effects of refinement (Section 3.3). Four popular F0 extractors were also evaluated for reference; YIN [10], SWIPE ′ [12], NDF [11] and DIO [22,34]. They were tested using their default or recommended settings. A constant F0 trajectory was used in this test. Figure 5 shows the results for a 120 Hz F0. The vertical axis represents the relative RMS error. When the SNR is larger than 5 dB, YIN yielded the best results. But, YIN's performance is obtained at the cost of poor temporal resolution, which will be shown in the following test. DIO was designed for high-quality recordings and is not tolerant to noise. While SWIPE ′ showed good performance from 0 to 20 dB SNR, performance saturated there after. The harmonic refinement procedure reduced the error in the initial estimate by a factor of 8, even in high noise, 4 The original L-F model [15] is anti-aliased using a closed form representation. The MATLAB implementation of this function and GUIbased interactive application for speech science education are open source [17]. Spurious levels around the fundamental component of the model's output are lower than -120 dB.  Figure 6: Frequency modulation transfer function for F0 modulation. The higher tracking frequency limit of the initial F0 estimate (triangle) is expanded two times by the proposed refinement using F0 adaptive time warping (circle).
because the standard deviation of error in n-th harmonic component is 1/n as described in previous paragraph. In total, this is the second best result.

Frequency modulation of F0
Measuring the ability of a F0 tracker to follow F0 modulation is a more relevant test for speech sounds with rapid changes. The instantaneous frequency of the aliasing-free L-F model output was controlled at audio sampling rate (22050 Hz) resolution. The average F0 was 120 Hz with 100 musical cent peak-to-peak modulation depth roughly to 6% frequency modulation peak-topeak in frequency. 5 In the two tests described in this section, a bit of white noise (SNR 100 dB) was added.   Figure 6 shows the frequency modulation transfer function for the four F0 trackers that serve as a benchmark and two variations of the F0 tracker described in this paper. For very low vibrato frequency (low modulation frequency) all F0 trackers work well at high SNR. At higher modulation frequencies all F0 trackers except for T10 • T10 • H3 fail to follow the full modulation, which shows up as a reduced gain when considering the output vs input modulation deviation. For higher F0 signals, the 3 dB point increased proportionally to the F0 value, except YIN. Figure 7 shows the RMS error of the F0 trajectories as a function of the modulation frequency. The dashed line and dash dot line show the RMS error of the best approximation to the true F0 using piece-wise linear function with segment lengths 1 ms and 5 ms respectively. SWIPE ′ and YIN yielded large RMS error, corresponding to the strong distortion shown in Fig. 2. The refinement performance without time warping is comparable to NDF. DIO showed the best performance among popular methods. The refined F0 trajectory using F0 adaptive time warping reduces the RMS error by a factor of 10 or more over the range from 2 Hz to 16 Hz modulation. For higher F0 values, RMS errors of other methods decrease inversely proportionally to the F0 value.
The F0 adaptive time warping also reduced spurious component due to FM substantially. For example, for a test signal with 16 Hz frequency modulation and 100 musical cent p-p depth, the refined F0 by the analysis configuration T10 •T10 •H3 reduced spurious residual levels lower than −40 dB. This is perceptually negligible.

Discussion
The goal of this paper is to estimate F0 trajectories, which consist of rapidly changing components, accurately for high-quality speech synthesis. The proposed set of procedures provide a prospective framework. However, the following aspects of F0 estimation were not exploited here. Investigations of the following issues could be important for improving synthesis quality further.
Plosive sounds such as /k/, /t/ sometimes sound like fricative by smearing temporal sharpness due to the smoothing effect of time windowing. This is a common degradation found in STRAIGHT.
Some speakers and languages frequently use "creaky voice." Representing these sounds using periodic signal plus noise results in poor reproduction. Relevant analysis and representations have to be investigated.
Temporal variation of F0 consists of effects caused by interactions between harmonic components and group delay in vocal tract transfer function. It is desirable to compensate this effect for speech synthesis applications, because this effect can be accumulated in each analysis and synthesis cycle.

Conclusions
This paper introduced a framework for intantaneous estimates F0 and aperiodicity. It is able to improve the ability of F0 extractors to temporally follow varying F0 trajectories by a factor of 10. It may serve as an useful infrastructure for speech research and applications.

Acknowledgements
The authors appreciate insightful discussions with Prof. Roy Patterson on human auditory perception, especially on fine temporal structure and detection of interfering sounds. Malcolm Slaney provided editorial assistance. He and Dan Ellis also provided productive as well as critical comments.

A. Note on the Flanagan's equation
Flanagan uses the time derivative of the logarithm of a complex signal x(t) to estimate the instantaneous frequency. By introducing a logarithmic function, the phase component is linearly separable from amplitude.
To make derivation simpler, as far as no ambiguity is introduced, time dependency representation by (t) is omitted afterwards.
Which is the Flanagan's equation.
The complex-valued signal x in Eq. 12 is a filtered output of h(t). It is a function of the center frequency fc and time. Let explicitly represent x using X(ωc, t) and its time derivative using X d (ωc, t). Then the following holds. where Substituting these two time windows w(t) and w d (t) into Eq. 12 removes time derivatives: Note that the TKEO (Teager Kaiser Energy Operator [35]) is not relevant for estimating rapidly changing F0 trajectories, since it uses an approximation, which requires slowly changing AM and FM. Using Flanagan's equation is relevant, since it does not rely on this approximation.

B. Residual calculation in each detector
This section shows how the aperiodicity detector in Fig. 3 works. The input to this detector is x(t). Let h(t, fc) represent the complex valued impulse response of each band pass filter centered around fc. a k (t, fc) = |r(t, fc)| 2 (18) y2(t, fc) = where the integration interval (−2/fc, 2/fc) is for the Nuttall window (Eq. 25). For Hann window the interval is (−1/fc, 1/fc) and for Blackman window the interval is (−1.5/fc, 1.5/fc). Band pass filters having these impulse response lengths have first spectral zeros at 0 and 2fc.  Smoothing the relative residual level a k (t, fc) yields the aperiodicity parameter a ks (t, fc). a ks (t, fc) = 2/fc −2/fc |h(τ, fc)|a k (t − τ, fc)dτ.
B.1. Operation and implementation of the procedure Figure 8 illustrates the process use to calculate the aperiodicity component. The impulse response of the filter h(t, fc) is h(t, fc) = w(t) exp(2πjfct), where j = √ −1 and the coefficients {a k } 3 k=0 are (0.338946, 0.481973, 0.161054, 0.018027). This is the 11-th item in Table  II of Nuttall's work [36].  The detector is designed to cancel the primary periodic component in the input signal by adjusting the filter gain at the frequency of the primary component. This is done by normalizing the output by its RMS level. In a high SNR case, total RMS level of the filtered signals are approximately equal to the RMS level of the periodic component. The RMS level of the lower level components are affected by this suppression process. Since the equivalent filter gain from this suppression process is the difference of two filters, it yields the filter shape shown in the red curve of left plot of Fig. 8. The right plot of Fig. 8 shows the output aperiodicity parameter a k as a function of the location of the primary component and the level of the lower level components. Figure 9 shows detector filter shapes of front end and the third stage. In the front end, the filter width is proportional to the center frequency. In the refinement stage, the filter width is constant. The filters in the refinement stage are designed using the estimated F0.

D. Mixing F0 information
The band of estimators in the front end independently estimate instantaneous frequency and an estimate of the quality of this estimate in the form of an aperiodicity measure. We need to consolidate these estimates to get a single estimate of F0 and we do this with a weighted average. Assume a set of random variables X k , k = 1, . . . , N having zero mean (E[X k ] = 0) and variances σ 2 k (Var[X k ] = σ 2 k ). We wish to generate a new estimate from all the noisy estimates by weighting the individual estimates to arrive at an answer with the minimum estimated variance. Thus, assume the following cost function.
where b k represents the mixing coefficient. When mixing F0 estimates derived from different sources, the sum of weights has to satisfy the condition ( N k=1 b k = 1). The optimum coefficientsâ k for k = 1, . . . , N − 1 are derived by solving the following set of equations.
bn for (k = 1, . . . , N − 1). (28) The final coefficientâN is given bŷ Other source of F0 information can be used to improve this estimate further, if the variance of the estimate is available.