of Institute Phonetic Sciences, University of Amsterdam, Proceedings 17 (1993), 97-110.
ACCURATE SHORT-TERM ANALYSIS OF THE FUNDAMENTAL FREQUENCY AND THE HARMONICS-TO-NOISE RATIO OF A SAMPLED SOUND
Paul Boersma
Abstract
We present a straightforward and robust algorithm for periodicity detection, working in the lag (autocorrelation) domain. When it is tested for periodic signals and for signals with additive noise or jitter, it proves to be several orders of magnitude more accurate than the methods commonly used for speech analysis. This makes our method capable of measuring harmonics-to-noise ratios in the lag domain with an accuracy and reliability much greater than that of any of the usual frequency-domain methods.
By definition, the best candidate for the acoustic pitch period of a sound can be found from the position of the maximum of the autocorrelation function of the sound, while the degree of periodicity (the harmonics-to-noise ratio) of the sound can be found from the relative height of this maximum. However, sampling and windowing cause problems in accurately determining the position and height of the maximum. These problems have led to inaccurate timedomain and cepstral methods for pitch detection, and to the exclusive use of frequency-domain methods for the determination of the harmonics-to-noise ratio. In this paper, I will tackle these problems. Table 1 shows the specifications of the resulting algorithm for two spectrally maximally different kinds of periodic sounds: a sine wave and a periodic pulse train; other periodic sounds give results between these.
Table 1. The accuracy of the algorithm for a sampled sine wave and for a correctly sampled periodic pulse train, as a function of the number of periods that fit in the duration of a Hanning window. These results are valid for pitch frequencies up to 80% of the Nyquist frequency. These results were measured for a sampling frequency of 10 kHz and window lengths of 40 ms (for pitch) and 80 ms (for HNR), but generalize to other sampling frequencies and window lengths (see section 5).
Periods per window >3 >6 > 12 > 24
Pitch determination error F/F sine wave pulse train < 5104 < 5105 < 3105 < 5106 7 < 410 < 2107 < 2108 < 2108
Resolution of determination of harmonics-to-noise ratio sine wave pulse train > 27 dB > 12 dB > 40 dB > 29 dB > 55 dB > 44 dB > 72 dB > 58 dB
IFA Proceedings 17, 1993
97
1 Autocorrelation and periodicity
For a time signal x(t) that is stationary (i.e., its statistics are constant), the autocorrelation rx() as a function of the lag is defined as r x ( ) x (t ) x (t + ) dt (1)
This function has a global maximum for = 0. If there are also global maxima outside 0, the signal is called periodic and there exists a lag T0, called the period, so that all these maxima are placed at the lags nT0, for every integer n, with rx(nT0 ) = rx(0). The fundamental frequency F 0 of this periodic signal is defined as F0 = 1 T0 . If there are no global maxima outside 0, there can still be local maxima. If the highest of these is at a lag max , and if its height rx(max) is large enough, the signal is said to have a periodic part, and its harmonic strength R 0 is a number between 0 and 1, equal to the local maximum r x ( max ) of the normalized autocorrelation rx ( ) rx ( ) r x (0) (2)
We could make such a signal x(t) by taking a periodic signal H(t) with a period T0 and adding a noise N(t) to it. We can infer from equation (1) that if these two parts are uncorrelated, the autocorrelation of the total signal equals the sum of the autocorrelations of its parts. For zero lag, we have r x (0) = r H (0) + r N (0) , and if the noise is white (i.e., if it does not correlate with itself), we find a local maximum at a lag max = T0 with a height r x ( max ) = r H (T0 ) = r H (0). Because the autocorrelation of a signal at zero lag equals the power in the signal, the normalized autocorrelation at max represents the relative power of the periodic (or harmonic) component of the signal, and its complement represents the relative power of the noise component: r x ( max ) = r H (0) r (0) ; 1 - r x ( max ) = N r x (0) r x (0) (3)
This allows us to define the logarithmic harmonics-to-noise ratio (HNR) as HNR (in dB) = 10 10 log r x ( max ) 1 - r x ( max ) (4)
This definition follows the same idea as the frequency-domain definitions used by most other authors, but yields much more accurate results thanks to the precision with which we can estimate rx(). For perfectly periodic sounds, the HNR is infinite. For non-stationary (i.e., dynamically changing) signals, the short-term autocorrelation at a time t is estimated from a short windowed segment of the signal centred around t. This gives estimates F0(t) for the local fundamental frequency and R0(t) for the local harmonic strength. If we want these estimates to have a meaning at all, they should be as close as possible to the quantities derived from equation (1), if we perform a short-term analysis on a stationary signal. Sections 2 and 3 show how to cope with the windowing and sampling problems that arise. Section 4 presents the complete algorithm. Sections 5, 6, and 7 investigate the performance of the algorithm for three kinds of stationary signals: periodic signals without perturbations, with additive noise, and with jitter.
98
IFA Proceedings 17, 1993
2 Windowing and the lag domain
Candidates for the fundamental frequency of a continuous signal x(t) at a time tmid can be found from the local maxima of the autocorrelation of a short segment of the sound centred around tmid . In figure 1, we summarize the algorithm for the speech-like signal x(t) = (1 + 0.3 sin 2 140 t) sin 2 280 t, which has a fundamental frequency of 140 Hz and a strong `formant' at 280 Hz. The algorithm runs as follows: Step 1. We take from the signal x(t) a piece with duration T (the window length, 24 ms in figure 1), centred around tmid (12 ms in figure 1). We subtract from this piece its mean x and multiply the result by a window function w(t), so that we get the windowed signal a(t) = x(tmid - 1 T + t) - x w(t) 2
(
)
(5)
The window function w(t) is symmetric around t = 1 T and zero everywhere outside 2 the time interval [0, T]. Our choice is the sine-squared or Hanning window, given by w(t ) =
1 2
- 1 cos 2
2 t T
(6)
We will see how the Hanning window compares to several other window shapes. Step 2. The normalized autocorrelation ra() (we suppress the primes from now on) of the windowed signal is a symmetric function of the lag :
T -
ra ( ) = ra (- ) =
a(t) a(t + ) dt
0 T
(7)
2
a
0
(t) dt
In the example of figure 1, we can see that the highest of these maxima is at a lag that corresponds to the first formant (3.57 ms), whereas we would like it to be at a lag that corresponds to the F 0 (7.14 ms). For this reason, Hess (1992) deems the autocorrelation method "rather sensitive to strong formants". Moreover, the skewing of the autocorrelation function makes the estimate of the lag of the peak too low, and therefore the pitch estimate too high (e.g., for 3 periods of a sine wave in a Hanning window, the difference is 6%). One method commonly used to overcome the first problem, is to filter away all frequencies above 900 Hz (Rabiner, 1977), which should kill all formants except the first, and estimate the pitch from the second maximum. This is not a very robust method, because we often run into higher formants below 900 Hz and fundamental frequencies above 900 Hz. Other methods to lose the formant include centre clipping, spectral flattening, and so on. Such ad-hoc measures render the method speech- and speaker-dependent. All these patches to the autocorrelation method are unnecessary, for there is a simple remedy: Step 3. We compute the normalized autocorrelation rw() of the window in a way exactly analogous to equation (7). The normalized autocorrelation of a Hanning window is
2 1 2 1 2 rw ( ) = 1 - + cos sin + T 3 3 T 2 T
(8)
IFA Proceedings 17, 1993
99
x(t)
1
multiplied by
1
w(t)
gives
1
a(t)
-1
0
Time (ms) ->
24
0
0
Time (ms) ->
24
-1
0
Time (ms) ->
24
ra ()
1
divided by
1
rw()
gives
1
rx()
-1
0
7.14 Lag (ms) ->
24
0
0
Lag (ms) ->
24
-1
0
7.14 Lag (ms) ->
24
Fig. 1. How to window a sound segment, and how to estimate the autocorrelation of a sound segment from the autocorrelation of its windowed version. The estimated autocorrelation rx() is not shown for lags longer than half the window length, because it becomes less reliable there for signals with few periods per window.
To estimate the autocorrelation rx() of the original signal segment, we divide the autocorrelation ra( ) of the windowed signal by the autocorrelation rw( ) of the window: rx ( ) ra ( ) rw ( ) (9)
This estimation can easily be seen to be exact for the constant signal x(t) = 1 (without subtracting the mean, of course); for periodic signals, it brings the autocorrelation peaks very near to 1 (see figure 1). The need for this correction seems to have gone by unnoticed in the literature; e.g., Rabiner (1977) states that "no matter which window is selected, the effect of the window is to taper the autocorrelation function smoothly to 0 as the autocorrelation index increases". With equation (9), this is no longer true. The accuracy of the algorithm is determined by the reliability of the estimation (9), which depends directly on the shape of the window. For instance, for a periodic pulse train, which is defined as x (t ) =
n=-
(t - t0 - nT0 )
+
(10)
where T0 is the period and t 0 /T (with 0 t0 < T0 ) represents the phase of the pulse train in the window, our estimate for the relevant peak of the autocorrelation is r x (T 0 ) =
w(t0 + nT0 ) w(t0 + (n + 1)T0 )
n
rw (T0 ) w 2 (t0 + nT0 )
n
(11)
This depends on the phase t 0 /T. If the window is symmetric and the pulse train is symmetric around the middle of the window, the derivatives with respect to t0 of both
100
IFA Proceedings 17, 1993
1
1
0
Time ->
0
Time ->
Fig. 2. Example of a windowed signal showing the two phases of a pulse train that yield extrema in the HNR estimation of the autocorrelation peak at a lag that equals the period.
the numerator and the denominator are 0; the extrema of rx(T0 ) as a function of t0 are thus found for the two phases exemplified in figure 2 for 3.0 periods per window. If such an extremum is greater than 1, it must be reflected through 1 to give a mathematically possible value of the autocorrelation, e.g., an initial estimate of 1.01 must be converted to 1/1.01 before computing its final HNR estimate, which is 20 dB. Figure 3 shows the worst-case HNR values for a perfectly periodic pulse train, calculated with equation (11) for a Hanning window, and for the rectangular window w(t ) = 1 ; r w ( ) = 1 - and for the Welch window w(t ) = sin
T
(12)
t T
1 ; rw ( ) = 1 - cos + sin T T T
(13)
as well as for the Hamming window 2 t T 2 2 1 sin + 0.3910 1 - 0.2916 + 0.1058cos T T T 2 rw ( ) = 0.3974 w(t ) = 0.54 - 0.46 cos
(14)
As we can see from figure 3, the Hanning window performs much better than the other three window shapes. Furthermore, the Hanning window is the `narrowest' of the four window shapes, which makes it the least vulnerable of the four to rapidly changing sounds. That makes two reasons for forgetting about the other three. In our implementation, the autocorrelations of the windowed signal and the window are numerically computed by Fast Fourier Transform. This is possible thanks to the fact that the autocorrelation can be obtained by first computing the Fourier transform of the windowed signal, which gives in the frequency domain ~ a( ) = a(t ) e -i t dt
2
(15)
~ and then computing the inverse Fourier transform of the power density a( ) , which brings us to the lag domain ~ ra ( ) = a( ) ei
2
d 2
(16)
IFA Proceedings 17, 1993
101
80 Hanning
60 Smallest HNR (dB) -> Hamming (drawn), Welch (dotted) 40 rectangular
20
0 3.0
6.0
12.0 24.0 Number of periods per window ->
48.0
Fig. 3. The sensitivity of several window shapes to the phase of a pulse train. Every point on a curve represents the worse of the two HNR values belonging to the phases in fig. 2.
This procedure allows us to try two other functions in the lag domain, besides the autocorrelation. The first of these functions is what we will call the `zero-phased' windowed signal, which is the sum of all the Fourier components of a(t), reduced to cosines with a starting phase of 0. This function is obtained by computing in the ~ frequency domain the absolute value a( ) instead of the power density. This conserves the relative amplitudes of the components of the windowed signal, which is nice because it gives the formants the peaks that they deserve. The second function is known as the cepstrum (Noll, 1967) and is obtained by computing in the frequency domain the logarithm of the power: ~ log 1 + c a( )
(
2
)
(17)
for large enough c > 0. The cepstral pitch-detection tactic was very common in the days that equation (9) was unknown, because it was the only one of the three methods that could raise the second peak of ra() (see figure 1) above the first peak. However, for both the zero-phased signal and the cepstrum, the addition of noise strongly suppresses all peaks relatively to the one at zero lag to a degree that depends on the frequency distribution of the noise. This makes these two methods unsuitable both for voiced-unvoiced decisions and for determination of the harmonics-to-noise ratio in the lag domain. Also, the pitch estimates are less accurate by several orders of magnitude as compared to the autocorrelation method. With equation (9) at our disposal, the advantages of these two alternative methods have vanished.
102
IFA Proceedings 17, 1993
3 Sampling and the lag domain
Consider a continuous time signal x(t) that contains no frequencies above a certain frequency fmax. We can sample this signal at regular intervals t 1 (2 f max ) so that we know only the values xn at equally spaced times tn: xn = x (tn ) ; tn = t0 + nt (18)
We lose no data in this sampling, because we can reconstruct the original signal as x(t) =
n=-
+
xn
sin (t - tn ) / t (t - tn ) / t
(19)
The autocorrelation computed from the sampled signal is also a sampled function: rn = r (n ) There is a local maximum in the autocorrelation between (m1) and (m+1) if rm > rm-1 and rm > rm+1 (20)
A first crude estimate of the pitch period would be max m , but this is not very accurate: with a sampling frequency of 10 kHz and = t , the pitch resolution for fundamental frequencies near 300 Hz is 9 Hz (which is the case for most time-domain pitch-detection algorithms); moreover, the height of the autocorrelation peak (rm) can be as low as 2/ = 0.636 for correctly sampled pulse trains (i.e., filtered with a phasepreserving low-pass filter at the Nyquist frequency prior to sampling), which renders HNR determination impossible and introduces octave errors in the determination of the fundamental period. We can improve this by parabolic interpolation around m :
max
1 r ( -r ) m + 2 m+1 m-1 2rm - rm-1 - rm+1
; rmax rm +
8(2rm - rm-1 - rm+1 )
(rm+1 - rm-1 )2
(21)
However, though the error in the estimated period reduces to less than 0.1 sample, the height of the relevant autocorrelation peak can still be as low as 7/(3) = 0.743. Now for the solution. We should use a `sin x / x' interpolation, like the one in equation (19), in the lag domain (we do a simple upsampling in the frequency domain, so that = t/2). As we cannot do the infinite sum, we interpolate over a finite number of samples N to the left and to the right, using a Hanning window again to taper the interpolation to zero at the edges: r( ) rnr -n
n=1 N N
(22) sin ( r + n - 1) 1 1 ( r + n - 1) rnl +n ( + n - 1) 2 + 2 cos + N r r n=1 where nl largest integer ; nr nl + 1 ; l - nl ; r 1 - l
sin ( l + n - 1) 1 1 ( l + n - 1) 2 + 2 cos + N + ( l + n - 1) l
IFA Proceedings 17, 1993
103
In our implementation, N is the smaller of 500 and the largest number for which (nl + N) is smaller than half the window length. This is because the estimation of the autocorrelation is not reliable for lags greater than half the window length, if there are few periods per window (see figure 1). Note that the interpolation can involve autocorrelation values for negative lags. The places and heights of the maxima of equation (22) can be determined with great precision (they are looked for between (m1) and (m+1) ). We can show this with long windows, where the windowing effects have gone, but the sampling effects remain. E.g., with a 40-ms window, any signal with a frequency of exactly 3777 Hz, sampled at 10 kHz, will be consistently measured as having a fundamental frequency of 3777.00000 0.00001 Hz (accuracy 108 sample in the lag domain, N=394) and a first autocorrelation peak between 0.99999999 and 1. The measured HNR (80-ms window) is 94.0 0.1 dB. This looks like a real improvement.
4 Algorithm
A summary of the complete 9-parameter algorithm, as it is implemented into the speech analysis and synthesis program praat, is given here: Step 1. Preprocessing: to remove the sidelobe of the Fourier transform of the Hanning window for signal components near the Nyquist frequency, we perform a soft upsampling as follows: do an FFT on the whole signal; filter by multiplication in the frequency domain linearly to zero from 95% of the Nyquist frequency to 100% of the Nyquist frequency; do an inverse FFT of order one higher than the first FFT. Step 2. Compute the global absolute peak value of the signal (see step 3.3). Step 3. Because our method is a short-term analysis method, the analysis is performed for a number of small segments (frames) that are taken from the signal in steps given by the TimeStep parameter (default is 0.01 seconds). For every frame, we look for at most MaximumNumberOfCandidatesPerFrame (default is 4) lag-height pairs that are good candidates for the periodicity of this This frame. number includes the unvoiced candidate, which is always present. The following steps are taken for each frame: Step 3.1. Take a segment from the signal. The length of this segment (the window length) is determined by the MinimumPitch parameter, which stands for the lowest fundamental frequency that you want to detect. The window should be just long enough to contain three periods (for pitch detection) or six periods (for HNR measurements) of MinimumPitch. E.g. if MinimumPitch is 75 Hz, the window length is 40 ms for pitch detection and 80 ms for HNR measurements. Step 3.2. Subtract the local average. Step 3.3. The first candidate is the unvoiced candidate, which is always present. The strength of this candidate is computed with two soft threshold parameters. E.g., if VoicingThreshold is 0.4 and SilenceThreshold is 0.05, this frame bears a good chance of being analyzed as voiceless (in step 4) if there are no autocorrelation peaks above approximately 0.4 or if the local absolute peak value is less than approximately 0.05 times the global absolute peak value, which was computed in step 2. Step 3.4. Multiply by the window function (equation 5). Step 3.5. Append half a window length of zeroes (because we need autocorrelation values up to half a window length for interpolation). Step 3.6. Append zeroes until the number of samples is a power of two. Step 3.7. Perform a Fast Fourier Transform (discrete version of equation 15), e.g., with the algorithm realft from Press et al. (1989). Step 3.8. Square the samples in the frequency domain.
104
IFA Proceedings 17, 1993
Step 3.9. Perform a Fast Fourier Transform (discrete version of equation 16). This gives a sampled version of ra(). Step 3.10. Divide by the autocorrelation of the window, which was computed once with steps 3.5 through 3.9 (equation 9). This gives a sampled version of rx(). Step 3.11. Find the places and heights of the maxima of the continuous version of rx(), which is given by equation 22, e.g., with the algorithm brent from Press et al. (1989). The only places considered for the maxima are those that yield a pitch between MinimumPitch and MaximumPitch. The MaximumPitch parameter should be between MinimumPitch and the Nyquist frequency. The only candidates that are remembered, are the unvoiced candidate, which has a local strength equal to (local absolute peak ) ( global absolute peak ) (23) R VoicingThreshold + max 0, 2 - SilenceThreshold (1 + VoicingThreshold ) and the voiced candidates with the highest (MaximumNumberOfCandidatesPerFrame minus 1) values of the local strength R r ( max ) - OctaveCost 2 log( MinimumPitch max ) (24)
The OctaveCost parameter favours higher fundamental frequencies. One of the reasons for the existence of this parameter is that for a perfectly periodic signal all the peaks are equally high and we should choose the one with the lowest lag. Other reasons for this parameter are unwanted local downward octave jumps caused by additive noise (section 6). Finally, an important use of this parameter lies in the difference between the acoustic fundamental frequency and the perceived pitch. For instance, the harmonically amplitude-modulated signal with modulation depth dmod x (t ) = (1 + dmod sin 2 Ft ) sin 4 Ft (25)
has an acoustic fundamental frequency of F, whereas its perceived pitch is 2F for modulation depths smaller than 20 or 30 percent. Figure 1 shows such a signal, with a modulation depth of 30%. If we want the algorithm's criterion to be at 20% (in order to fit pitch perception), we should set the OctaveCost parameter to (0.2) 2 = 0.04; if we want it to be low (in order to detect vocal-fold periodicity), say 5%, we should set it to (0.05)2 = 0.0025. The default value is 0.01, corresponding to a criterion of 10%. After performing step 2 for every frame, we are left with a number of frequencystrength pairs (Fni, R ni), where the index n runs from 1 to the number of frames, and i is between 1 and the number of candidates in each frame. The locally best candidate in each frame is the one with the highest R. But as we can have several approximately equally strong candidates in any frame, we can launch on these pairs the global path finder, the aim of which is to minimize the number of incidental voiced-unvoiced decisions and large frequency jumps: Step 4. For every frame n, pn is a number between 1 and the number of candidates for that frame. The values {pn | 1 n number of frames} define a path through the candidates: {(Fnpn, Rnpn) | 1 n number of frames}. With every possible path we associate a cost cost ({ pn }) =
numberOfFrames n=2
transitionCost Fn-1, pn -1 , Fnpn -
(
)
numberOfFrames n=1
Rnp
n
(26)
IFA Proceedings 17, 1993
105
where the transitionCost function is defined by (F = 0 means unvoiced) 0 transitionCost( F1, F2 ) = VoicedUnvoicedCost F1 2 OctaveJumpCost log F2 if F1 = 0 and F2 = 0 if F1 = 0 xor F2 = 0 (27) if F1 0 and F2 0
where the VoicedUnvoicedCost and OctaveJumpCost parameters could both be 0.2. The globally best path is the path with the lowest cost. This path might contain some candidates that are locally second-choice. We can find the cheapest path with the aid of dynamic programming, e.g., using the Viterbi algorithm described for Hidden Markov Models by Van Alphen & Van Bergem (1989). For stationary signals, the global path finder can easily remove all local octave errors, even if they comprise as many as 40% of all the locally best candidates (section 6 presents an example). This is because the correct candidates will be almost as strong as the incorrectly chosen candidates. For most dynamically changing signals, the global path finder can still cope easily with 10% local octave errors. For many measurements in this article, we turn the path finder off by setting the VoicedUnvoicedCost and OctaveJumpCost parameters to zero; in this way, the algorithm selects the locally best candidate for each frame. For HNR measurements, the path finder is turned off, and the OctaveCost and VoicingThreshold parameters are zero, too; MaximumPitch equals the Nyquist frequency; only the TimeStep, MinimumPitch, and SilenceThreshold parameters are relevant for HNR measurements.
1 Spectrum (dB) -> -1 1 Spectrum (dB) -> 0 0.01 80 60 40 20 0 Frequency (Hz) -> 5000
Time (seconds) ->
-20 0 80 60 40 20 0
-1
0
Time (seconds) ->
0.01
-20 0
Frequency (Hz) ->
5000
Fig. 4. At the left: two periodic signals, sampled at 10 kHz: a sine wave and a pulse train, which was squarely low-pass filtered at 5000 Hz (acausal, phase-preserving filter). Both have a fundamental frequency of 490 Hz. At the right: their spectra.
106
IFA Proceedings 17, 1993
5 Accuracy in measuring perfectly periodic signals
The formula for a sampled perfect sine wave with frequency F is xn = sin 2 Ftn (28)
and the formula for a correctly sampled pulse train (squarely low-pass filtered at the Nyquist frequency) with period T is xn = sin Fs (tn - mT ) m=- Fs (tn - mT )
+
(29)
These two functions form spectrally maximally different periodic signals. Figure 4 shows examples of these signals, together with their spectra. The spectrum of the sine wave is maximally narrow, that of the pulse train is maximally wide. Table 1 (page 97) shows our algorithm's accuracy in determining pitch and HNR. We see from table 1 that for pitch detection there should be at least three periods in a window. The value of 27 dB appearing in table 1 for a sine wave with the worst phase (symmetric in window) and the worst period (one third of a window), means that the autocorrelation peak can be as low as 0.995, which means that the signal xn = (1 + dmod sin 2 Ftn ) sin 4 Ftn (30)
(see also figure 1), whose fundamental frequency is F, can be locally ambiguous for F near MinimumPitch, if the modulation depth dmod is less than = 7%. The 10.995 critical modulation depth, at which there are 10% local octave errors (detection of 2F as the best candidate), is 5%, for the lowest F (equal to MinimumPitch). Note that the global path finder will not have any trouble removing these octave errors. We also see from table 1 that for HNR measurements, there should be at least 6.0 periods per window. The values measured for the HNR of a pulse train are the same as those predicted by theory for continuous signals, as plotted in figure 3. This suggests that the windowing effects have the larger part of the influence on HNR measurement inaccuracy, and that the sampling effects have been effectively cancelled by equation (22). For very short windows (less than 20 samples in the time domain, MinimumPitch greater than 30% of the Nyquist frequency), the HNR values for pulse trains do not deteriorate, but those for sine waves approach the values for pulse trains; the relative pitch determination error rises to 10-4. The problems with short-term HNR measurements in the frequency domain, are the sidelobes of the harmonics and the sidelobes of the Fourier transform of the window: they occur throughout the spectrum. Pitch-synchronous algorithms try to cope with the first problem, but they require prior accurate knowledge of the period (Cox et al., 1989; Yumoto et al., 1982). Using fixed window lengths in the frequency domain requires windows to be long: the shortest window used by Klingholz (1987) spans 12 periods, De Krom (1993) needs 8.2 periods; with the shortest window, both have a HNR resolution of apx. 30 dB for synthetic vowels, as opposed to our 37 dB with only 6 periods (48 dB for 8.2 periods, 52 dB for 12 periods). In the autocorrelation domain, the only sidelobe that could stir trouble, is the one that causes aliasing for frequency components near the Nyquist frequency; that one is easily filtered out. This is the cause of the superior results with the present method.
IFA Proceedings 17, 1993
107
60 50 measured HNR (dB) -> 40 30 20 10 0 -10 -30 measured HNR (dB) ->
60 50 40 30 20 10 0 -10 -30 -20 -10
-20
-10
0 10 20 SNR (dB) ->
30
40
50
0
10 20 30 SNR (dB) ->
40
50
60
Fig. 5. Measured HNR values for a sine wave (left) and a pulse train (right) sampled at 10 kHz, both with a periodicity of 103 Hz, with additive noise. The figures show the 10%, median, and 90% curves. The window length was 80 ms.
6 Sensitivity to additive noise
The formula for a sampled sound consisting of a sine wave with frequency F and additive `white' noise (squarely low-pass filtered at the Nyquist frequency) is xn = 2 sin 2 Ftn + 10 -SNR/ 20 z n (31)
where SNR is the signal-to-noise ratio, expressed in dB, and zn is a sequence of real numbers that are independently drawn from a Gaussian distribution with zero mean and unit variance. The formula for a sampled sound consisting of a correctly sampled pulse train (squarely low-pass filtered at the Nyquist frequency) with period T and additive `white' noise is xn = Fs / F 1 - F / Fs sin Fs (tn - mT ) + 10 -SNR/ ...