JAIST Repository
https://dspace.jaist.ac.jp/

Title

Estimation of fundamental frequency ofreverberant speech by utilizing complex cepstrum
analysis

Author(s)Unoki, Masashi; Hosorogiya, Toshihiro

Citation

Journal of Signal Processing, 12(1): 31‑44

Issue Date

2008‑01

Type

Journal Article

Text version

author

URL

http://hdl.handle.net/10119/7755

RightsCopyright (C) 2008 信号処理学会. Masashi Unoki
and Toshihiro Hosorogiya, Journal of Signal
Processing, 12(1), 2008, 31‑44.

Description

Japan Advanced Institute of Science and Technology

PAPER
Estimation of fundamental frequency of reverberant speech byutilizing complex cepstrum analysis
Masashi Unoki and Toshihiro Hosorogiya
School of Information Science, Japan Advanced Institute of Science and Technology
1-1 Asahidai, Nomi, Ishikawa 923-1292, Japan
E-mail: {unoki, t-hosoro}@jaist.ac.jp

Abstract

This paper reports comparative evaluations of twelve typical methods of estimating fundamental frequency (F0 ) over huge speech-sound datasets in artiﬁcial reverberant environments. They involve several classic algorithms such as Cepstrum, AMDF, LPC, and modiﬁed autocorrelation algorithms.
Other methods involve a few modern instantaneous amplitude- and/or frequency-based algorithms, such as
STRAIGHT-TEMPO, IFHC, and PHIA. The comparative results revealed that the percentage of correct
rates and SNRs of the estimated F0s were reduced drastically as reverberation time increased. They also
demonstrated that homomorphic (complex cepstrum) analysis and the concept of the source-ﬁlter model
were relatively eﬀective for estimating F0 from reverberant speech. This paper thus proposes a new methodof robustly and accurately estimating F0 s in reverberant environments, by utilizing the modulation transfer
function (MTF) concept and the source-ﬁlter model in complex cepstrum analysis. The MTF concept isused in this method to eliminate dominant reverberant characteristics from observed reverberant speech.
The source-ﬁlter model (liftering) is used to extract source information from the processed cepstrum. Finally, F0 s are estimated from them by using the comb-ﬁltering method. Additive-comparative evaluation
was carried out on the new approach with other typical methods. The results demonstrated that it was
better than the previously reported techniques in terms of robustness and providing accurate F0 estimatesin reverberant environments.

Keywords: F0 estimation, reverberant speech, complex cepstrum analysis, MTF concept, source-filter model

1.

Introduction

The fundamental frequency (F0 ) as well as the fundamental period (T0 ) of speech can be utilized as signiﬁcant features to represent the source information(glottal waveform or vocal-fold vibrations) of speechsounds in various speech-signal processes. These are,for example, in speech analysis/synthesis systems, automatic speech recognition (ASR) systems, and speechemphasis methods. In particular, robust and accurate
F0 can generally be used as a powerful cue to reduce
the noise component in noisy speech and/or to remove
the reverberation eﬀect in reverberant speech. Therefore, robustly and accurately estimating the F0 of target speech in real environments, which is the same as
extracting the F0 of noiseless speech, is a particularly
important issue in these applications.
Many studies on extracting or estimating the F0of target speech have been done in the literature on

speech-signal processing, and numerous methods havebeen proposed [1, 2, 3] over the last half-century. The
traditional extraction/estimation methods can be divided into processing in the time and frequency domains, or both domains. Most of these have madeuse of the periodic features of speech in the time domain (zero-cross [4, 5], periodgram [6], peak-picking[4, 7], autocorrelation [4, 8], the amplitude magnitude diﬀerence function (AMDF) [9], and maximum
likelihood [10]), or harmonic features in the frequency
domain (comb ﬁltering [11, 12], autocorrelation [13],sub-harmonic summation (SHS) [14], and cepstrum[15]).
The aim of all these methods has been to extract
the periodicity or harmonicity of source informationfrom observed speech. However, this still seems to beincompletely resolved because three main issues remain, i.e., (1) observability: the observed speech isan emitted sound passing through the mouth/nose so2that it is impossible to directly observe glottal vibrations from it without eliminating the eﬀects of the
vocal tract, (2) flexibility and irregularity: glottal
vibrations are not complete periodic signals and therange of variations in the periods is relatively wide,and (3) robustness: the observed speech signals areaﬀected by noise and reverberation so that the significant features for estimating F0 are also smeared.
Most studies have focused on the ﬁrst two issues
so that they have implicitly assumed all speech signals
are observed in clean environments or all observations
are only noiseless speech sounds. Various methodsof estimating F0 have been proposed under this assumption to solve the ﬁrst issue by suppressing the effects of ﬁlter characteristics (vocal tract), based on thesource-ﬁlter model, from the observed speech sounds.
For example, typical approaches based on this idea
have been homomorphic methods of analysis [15, 16]
and linear prediction (LP)-based methods [17, 18, 19].
A few examples of inverse ﬁltering methods are lagwindowing [20] and simpliﬁed inverse ﬁlter tracking(SIFT) [21]. Center-clipping, band-limitation [22, 23],and multi-windowing [24] techniques have also beenused in approaches based on the autocorrelation function.
A few approaches to precisely estimate the F0 oftarget noiseless speech have been established (e.g.,
STRAIGHT-TEMPO [25] and YIN [26]) by comparing electro-glottal-graph (EGG) information. The stability of the instantaneous frequency of speech has also
been used in the STRAIGHT-TEMPO method (referred to as “TEMPO” after this) to accurately estimate F0 s as signiﬁcant features to resolve the ﬁrst two
issues. This method plays an important role in controlling “pitch” related features in STRAIGHT analysis/synthesis tools [27]. YIN has also been proposed
which combines autocorrelation functions and AMDF
to resolve these. It has been reported that both methods can be used to estimate the F0 of target noiselessspeech extremely precisely so that the ﬁrst two issues
seem to be have been resolved. However, it has not
yet been clariﬁed whether these methods can precisely
estimate F0 in real environments. Hence, we need toinvestigate the last issue for realistic applications.
It is generally known that the method of estimating
F0 using periodic and/or harmonic features (e.g., autocorrelation functions and comb ﬁltering) is relatively
robust against background noise, but the estimatedF0 is not relatively accurate [2, 28, 29]. It has also
been reported that the comb-ﬁltering-based method
is more robust against background noise than the
autocorrelation-based approach [29]. The cepstrumbased method is not as robust against backgroundnoise as either of these because it is composed of homomorphic analysis so that noise components are notclearly separated in the quefrency domain [29].
The time-frequency representation of speech ob-

Journal of Signal Processing, Vol. , No. ,tained by time-frequency analysis can also adequately represent the periodic/harmonic componentsof speech. The instantaneous amplitude (IA) of speech
signals has ﬁne harmonic features that are robust
against background noise so that the comb-ﬁlteringof instantaneous amplitude has been proposed [30] toconstruct a sound-segregation model. The instantaneous frequency (IF) of speech has also been used toaccurately estimate F0 s but their stability as used in
TEMPO is sensitive to noise. More robust methods
using instantaneous frequency have been proposed byusing bandwidth equations related to instantaneous
amplitude and frequency with harmonicity [28, 31].
Other robust techniques using instantaneous amplitude and frequency-related approaches have been proposed by using periodicity and harmonicity [29]. It
has been reported that these are more robust than
TEMPO and can precisely estimate the F0 in noisy
environments.
All these methods have focused on noiseless tonoisy conditions to estimate suﬃciently accurate F0 s
of target speech. Thus, methods using instantaneous
amplitude and frequency or those with robust features
against noise such as periodicity and harmonicity havebeen regarded as accurately being able to estimate F0 s
from noisy speech. The last issue seems to be havebeen solved at this time; however, there have been nostudies on robustness in reverberant environments.
It can easily be predicted that no typical methods
will work as well and their percentage correct rates forF0 s will drastically be reduced as reverberation timeincreases. If our prediction is correct, the last issue
has not yet been completely solved and needs to beconsidered in reverberant environments and in noisy
reverberant environments. We evaluated traditional
methods of estimating F0 in terms of robustness andaccuracy in reverberant environments to investigate
this issue and discuss them in this paper. We thenpropose a method of estimating F0 from reverberant
speech without measuring impulse response in room
acoustics (i.e., blind method of estimating F0 ) by taking the characteristics of reverberation into consideration.
This paper is organized as follows. Section 2 describes the mathematical setup and then deﬁnes the
problem of estimating F0 from reverberant speech. We
discuss our evaluations of most typical methods of estimating F0 in reverberant environments in Section 3
and investigations into what the best model is. Section 4 introduces complex cepstrum analysis and investigates what the signiﬁcant features for robust estimates are. We then introduce the model concept(complex cepstrum analysis, the modulation transfer
function (MTF) concept, and source-ﬁlter model). We
ﬁnally propose a method of estimating F0 in reverberant environments. We discuss our evaluations of the
proposed method in Section 5 by comparing it with

3

other methods using the same simulations. Section 6
gives our conclusions and perspective regarding future
work.

Amplitude cepstrumCA(q,τ)Lifteringl(q)Cflt(q,τ)2.

Mathematical setup

2.1 Signal representation and STFT

Csrc(q,τ)Csrc(q,τ)A time-varying harmonic signal, x(t), can be represented as the analytic signal:
x(t) =ak (t) exp(jωk (t)t + θk (t))(1)

Cepstrum componentof filter characteristics(vocal tract)k∈K

where ak (t) is the instantaneous amplitude and θk (t)is the phase. Here, k denotes the harmonic index and
K is the number of harmonics so that ωk (t) can be
expressed as 2πkF0 (t). Fundamental frequency F0 (t)is an instantaneous frequency so that this should be
extracted from x(t) using instantaneous cues.
The short-term Fourier transform (STFT) is usually used [32] to analyze x(t) in any given short-term
segment (windowing processing):X(ω, τ ) =x(t)w(t − τ ) exp(−jωt)dt (2)A(ω, τ )==φ(ω, τ )=A(ω, τ ) exp(j arg φ(ω, τ ))|X(ω, τ )|
[X(ω, τ )]arctan[X(ω, τ )]

(3)
(4)
(5)

where w(t) is a window function and a short-term signal, x(t, τ ), is deﬁned as w(t−τ )x(t) for mathematicalconvenience. A(ω, τ ) is the amplitude spectrum andφ(ω, τ ) is the phase spectrum of X(ω, τ ).
The task of extracting/estimating fundamental frequency F0 (t) in this formulation is, therefore, to estimate the F0 in each short-term segment using theharmonicity of X(ω, τ ) or to estimate segmental T0 =1/F0 by using the periodicity of x(t, τ ). Thus, traditional methods based on waveform processing (e.g.,
zero-cross [4, 5], periodgrams [6], peak-picking [4, 7],
autocorrelation [4, 8], AMDF [9], maximum likelihood [10], STFT-based processes, and sub-harmonicsummation (SHS) [14]) estimate F0 (t) from x(t, τ ) orX(ω, τ ) by using periodicity or harmonicity.
2.2 Source-ﬁlter model
The source-ﬁlter model is a well-known concept toseparately represent the glottal (source information)and vocal-tract (ﬁlter information) characteristics forspeech production. Based on this concept, the observed clean speech signal, x(t), can be representedas(6)x(t) = e(t) ∗ vτ (t)where e(t) is the source signal related to glottal information and vτ (t) is the impulse response of the ﬁlter

Quefrency
Cepstrum component ofsource (glottal vibration)

Fig. 1 Separated representations of source and ﬁlter
characteristics in quefrency domain
related to the vocal tract at time τ . The asterisk “∗”denotes convolution. Note that the emission eﬀect has
been omitted from this formulation. Thus, Eq. (2) can
also be represented asX(ω, τ ) = S(ω, τ ) · V (ω, τ )(7)where S(ω, τ ) is the STFT of s(t, τ ) = w(t − τ )e(t)and V (ω, τ ) is that of v(t, τ ) = vτ (t). V (ω, τ ) represents ﬁlter characteristics so that the separationeﬀect of V (ω, τ ) is usually used to estimate F0 (t)from X(ω, τ ). Some traditional methods of estimation are inverse ﬁltering V −1 (ω, τ ) [21], whitening ofX(ω, τ ) using |V (ω, τ )| (or lag windowing) [20], andsubtraction on logarithmic processing log X(ω, τ ) =log S(ω, τ ) + log V (ω, τ ) [22, 23].
The linear prediction (LP) method is also one of
the most powerful techniques of analyzing speech signals. LP coeﬃcients have ﬁlter characteristics (allpole type) and the LP residue has source information. The LP coeﬃcients of x(t, τ ) can thus be used
as inverse ﬁltering V −1 (ω, τ ) in the source-ﬁlter model[17, 21]. The LP residue can also be used as a shortterm signal s(t, τ ) [19]. Waveform processing and
AMDF have also been incorporated [18].
2.3 Cepstrum representation
Cepstrum is also a well-known method of homomorphic analysis. The complex cepstrum of X(ω, τ )in Eq. (2) can be represented as
C(q, τ )= F −1 [log {|X(ω, τ )| exp(jφ(ω, τ ))}]

= F −1 log A(ω, τ ) + F −1 jφ(ω, τ )= CA (q, τ ) + Cφ (q, τ )(8)where F −1 [·] is the Fourier inverse transform, CA (q, τ )is the amplitude cepstrum, Cφ (q, τ ) is the phase cepstrum of C(q, τ ), and q denotes the quefrency (timedomain). The complex cepstrum of X(ω, τ ) in Eq. (7)4Journal of Signal Processing, Vol. , No. ,can also be represented as
C(q, τ )3.

= F −1 [log S(ω, τ )] + F −1 [log V (ω, τ )]= Csrc (q, τ ) + Cflt (q, τ )2.4 Problem with estimating F0
The task of estimating F0 in reverberant environments is to extract F0 (t) from reverberant speech signal y(t) or respective STFT Y (ω, τ ):y(t)= x(t) ∗ h(t) = e(t) ∗ vτ (t) ∗ h(t) (10)= S(ω, τ )V (ω, τ )H(ω, τ )3.1 Typical methods of estimating F0(9)where Csrc (q, τ ) is the complex cepstrum of sourceS(ω, τ ) and Cflt (q, τ ) is that of ﬁlter V (ω, τ ).
The amplitude cepstrum, CA (q, τ ), is generally
used in the traditional method so that CA,src (q, τ ) andCA,flt (q, τ ) are separately used for estimating F0 (t)from CA (q, τ ). Figure 1 outlines the concept underlying the source-ﬁlter model in the quefrency domain.
CA,flt (q, τ ) represents the dominant spectrum envelope of X(ω, τ ) (lower Fourier component in quefrency
domain) so that they are compactly located in thelower quefrency. In contrast, CA,src (q, τ ) represents
the dominant ﬁne structure of X(ω, τ ) so that they
are compactly located in the higher quefrency domain.
Therefore, the task of estimating F0 with this concept
is to ﬁnd the dominant quefrency from CA,src (q, τ ) or
to detect periodicity or harmonicity from CA,src (q, τ )by eliminating CA,flt (q, τ ) from CA (q, τ ). The lastprocessing is referred to as “liftering”. Typical approaches are Noll’s original method [15] and his clipstrum method [16].

Y (ω, τ )Evaluation of typical methods(11)

where h(t) is the impulse response and H(ω, τ ) is the
STFT of h(t) in room acoustics (reverberation). Note
that, H(ω, τ ) is actually required to present all characteristics (H(ω) = H(ω, τ )) by using a long-term
Fourier transform (LTFT) so that the length of analysis (at each τ ) should be more than the reverberationtime.
The task of estimating F0 in reverberant environments is thus to select periodicity and harmonicity
from the convolved source signal, e(t), while that innoisy environments is to select them from the noisy(additive) source signal, e(t). If h(t) is simpliﬁed echo
or a minimum phase impulse response, the cepstrumbased method can be used to adequately estimate F0
from the reverberant speech signal, y(t), because homomorphic analysis is a powerful tool for dealing withsimpliﬁed echos. Realistic impulse responses in room
acoustics generally have non-minimum phase characteristics and we therefore predicted that estimating F0
robustly and accurately would be more diﬃcult thanin noisy environments.

Many methods of estimating F0 have been proposed in the literature on speech signal processing,as described in Section 1. The most comprehensive
review remains that by Hess (1983) [1] and more recent reviews are those by Hess (1992) and Cheveigné
and Kawahara (2001) [2, 3]. A few examples ofrecent approaches are instantaneous-amplitude [30],
instantaneous-frequency [28, 31], and fundamental
wave-ﬁltering [33]. There are also comparative evaluations in Atake et al.’s (2000), Ishimoto et al.’s (2001),
and Nakatani and Irino (2004) [2, 3, 28, 29, 31].
We evaluated twelve typical methods to investigate how robust estimates of F0 were in reverberant
environments: ACMWL (AutoCorrelation throughMultiple Window-Length) [24], AMDF [9], STFTACorrLog (AutoCorrelation of Log-amplitude spectrum on STFT) [22, 23, 13], STFT-ACorrLag (Lagwindowing on STFT) [20], STFT-Comb (Comb ﬁltering on STFT) [11, 12], SHS [14], Cepstrum [15], LPCresidue [19], VFWFF (Voice Fundamental Wave Filtering (Feed-forward type)) [33], TEMPO [25], IFHC(Instantaneous Frequency of Harmonic Components)[28], and PHIA (Periodicity/Harmonicity using Instantaneous Amplitude) [29]. Although other methods have been proposed, we choose these twelve because they are commonly used in comparative evaluations and the others are just modiﬁcations or heavy
revisions of them.
The characteristics, parameters, and detailed conditions on these algorithms that we used in this comparative evaluation are listed in Table 1. TEMPO
was used as a complete original version. As the othermethods implemented by the researchers were basedon their original research, they were the ﬁrst or original versions. Several parameters were set to obtainappropriate results based on preliminary simulations.
For more detailed information on their technical implementations and skill levels, please refer to their
original papers.
3.2 Sound dataset and evaluation measures
The sound dataset we used in this evaluation was
the speech database on simultaneous recordings ofspeech and EGG by Atake et al. [28]. This dataset
consisted of 30 short Japanese sentences uttered by 14
males and 14 females with voiced-unvoiced labels (total of 840 utterances, sampling frequency of 16 kHz,and quantization of 16-bits).
The reverberant speech sentences we used were created by convolving the original signals, x(t)s, with thefollowing reverberant impulse responses, h(t)s, as a

5

Table 1 Characteristics of typical methods of estimating F0 and their parameter settings
Algorithm
ACMWL [24]
AMDF [9]STFT
auto-corr. [22, 23, 13]
Lag windowing [20]Comb ﬁltering [11, 12]
SHS [14]domain
time
time
freq.
freq.
freq.
freq.
freq.

Featuresx(t, τ )x(t, τ )log |X(ω, τ )||S(ω, τ )||X(ω, τ )|log |X(ω, τ )|Cepstrum [15]quef.

CA (q, τ )LPC Residue [19]F0 ﬁltering [33]

time
time

s(t, τ )s(t, τ )TEMPO [25]freq.

IFHC [28]freq.

Fixed point analysis
on Instant. freq. (IF)Harmonicity of IFs

PHIA [29]

time/
freq.

Instant. Amp. (IA)& Dempster’s law

Parameter setting
Window candidates of nine (15 to 60-ms in 5-ms steps)
40-ms window, 2-ms shift
40-ms window, 2-ms shift
1.5-kHz band-limitation, averaged-clipping
3-ms lag windowing, 1.5-kHz band-limitation, averaged clipping
10th order of the comb ﬁlter
40-ms window length, 5-ms shift, N = 15 (harmonic order),BL = 1.25 kHz (band limitation), weigh function: hn = 0.84n−1
40-ms window, 2-ms shift, 3-ms liftering,
Noll’s method (hw = 0.54 + 0.46 cos(2πf /600))
40-ms window, 2-ms shift, LP order of 12, 1/4 down-sampling
80-ms window, 2-ms shift, 1/4 down-sampling, averaged-clipping,
1.5 kHz band-limitation, “Non-selfsupport mode”,
1-ms F0 shift length, nvo = 24, Heuristic factor: “on”
Nm = 3 (harmonic order), 1/8 down-sampling,
adaptive-windowing based on BW equation
400-channel CBFB, 64-channel CQFB,
1/4 down-sampling, averaged-clipping(1) Hanning windowing function was used in STFT-based processing.
(2) Lower limits of estimated F

0 were set at 60 Hz and higher were set at 800 Hz.

function of the reverberation time.


−6.9t
n(t)
h(t) = a exp
TR
 
1/2


T
−13.8t
a =
1
exp
dt
TR
0

(12)
(13)

where the “a” is a gain factor as the normalized
power of h(t), TR is reverberation time, and n(t)is white noise. This is the well-known stochasticapproximated impulse response in room acoustics reported by Schroeder [34, 35], in which the response hasan envelope of exponential decay and is a white noise
carrier. This formulation for the impulse response has
been used in the study of speech intelligibility in room
acoustics [36] as general artiﬁcial reverberation andthus has non-minimum phase components [34, 35]. Six
reverberation conditions (TR = 0.0, 0.1, 0.3, 0.5, 1.0,
and 2.0 s) were used in this study. There were a totalof 5, 040 stimuli.
Fine F0 error and gross F0 error within the voiced
section have been used as measures for some comparative evaluations in noisy environments [2, 28, 29, 31].
These have concentrated on error analysis. Since weconcentrated on evaluating robustness and the accuracy of F0 estimates, we used two similar measures forevaluation but not the same measures. The ﬁrst was
the percent correct rate (%) and the second was SNR(in dB).
Correct rateE =SNR = 20 log10

NF0,Est (E)× 100
NF0,Ref(F0,Ref (t) − F0,Est (t))2 dt
F0,Ref (t)2 dt

(14)

(15)

where F0,Ref (t) and F0,Est (t) are reference F0 and estimated F0 , and this integral is done in the voiced
section (t). NF0,Est (E) is the size of the correct region
that satisﬁes|F0,Ref (t) − F0,Est (t))|≤E
F0,Ref (t)(%)

within voiced section (t) where E is the error margin (%). NF0,Ref (E) is the size of region F0,Ref (t) inthe voiced section at the error margin E. In this paper, the F0 estimated by TEMPO from the EGG signal is used as the correct F0 (reference F0 , F0,Ref (t)).
F0,Est (t) was used to estimate F0 with the twelve
methods from reverberant speech signals. Two values for E (error margins of 5% and 10%) were used in
the percent correct rate.
Since gross F0 error is the ratio of the number offrames giving “incorrect” F0 values to the total number of frames, the percent correct rate approximately
indicates gross F0 error. Since ﬁne F0 error is thenormalized room mean square error between F0,Ref (t)and F0,Est (t), SNR indicates a similar measure in dB.
3.3 Results
Figure 2 plots the results of comparative evaluations for the twelve typical methods of estimating F0
from reverberant speech as a function of the reverberation time, TR . The left panels (a), (c), and (e) plot
the results for the ﬁrst six methods and the right panels (b), (d), and (f) plot them for the last six. The
top panel plots the percent correct rates (expressedas percentages) for F0 estimates within an error margin of 5 % and the middle panel plots these within

6

Journal of Signal Processing, Vol. , No. ,

90

90

80

80

70

70

Correct rate (%)

100

Correct rate (%)

100

60
50
40
30

10

50
40
30

Cepstrum
ACMWL
AMDF
STFT−ACorrLog
STFT−ACorrLag
STFT−Comb

20

60

20

(a)
10

5 % error margin

100

100

90

90

80

80

70

70

Correct rate (%)

0

Correct rate (%)

0

60
50
40
30
20
10

10

30
Cepstrum
ACMWL
AMDF
STFT−ACorrLog
STFT−ACorrLag
STFT−Comb

25

SHS
VFWFF
LPC−residual
IFHC
PHIA
TEMPO

25

20

SNR (dB)

20

SNR (dB)

(d)

0

30

15

10

15

10

5

0

10 % error margin

40

20

0

SHS
VFWFF
LPC−residual
IFHC
PHIA
TEMPO

(b)

50

(c)
10 % error margin

5 % error margin

60

30

Cepstrum
ACMWL
AMDF
STFT−ACorrLog
STFT−ACorrLag
STFT−Comb

SHS
VFWFF
LPC−residual
IFHC
PHIA
TEMPO

5

(e)

0

0.5

1

Reverberation time TR (s)

1.5

2

0

(f)

0

0.5

1

Reverberation time TR (s)1.5

2

Fig. 2 Estimation results: (a)–(b) percent correct rate within error margin of 5 %, (c)–(d) percent correct rate
within error margin of 10 %, and (e)–(f) SNR (s: original, n: error between original and estimated F0 ) of F0
estimates from reverberant speech using twelve typical methods as function of reverberation time, TR

7an error margin of 10 %. The bottom panel plots the
SNRs. The correct rates and SNRs of all 12 methods are drastically reduced as the reverberation timeincreases. The correct rates within the 5 % error margin for all methods were less than 50 % and the SNRs
were less than about 15 dB, especially when reverberation time TR was 2.0 s. Moreover, the correct rates
within the 10 % error margin as an approximate evaluation were also less than 70 %. We hence concluded
that none of these methods worked as well as robust
and accurate F0 estimates and they had drawbacks inestimating F0 from reverberant speech.
However, we found a few clues in doing this evaluation for improving these methods. We can see from
Fig. 2 that the cepstrum method is the most accurate excluding the clean condition (TR = 0.0). Cepstrum analysis is homomorphic and this can deal withconvolution processing as additive (subtractive) processing. Although the impulse responses we used in
evaluations were not minimum-phase characteristics,
the cepstrum method seemed to reduce the eﬀect ofreverberation for estimating F0 since this can treat a
direct sound and a reﬂected sound as the same signal.
Therefore, the cepstrum method has the possibility of
estimating F0 from reverberant speech if it is not affected too much by reverberation. The comb-ﬁltering
method is slightly more robust against reverberation
as we can see from Figs. 2(c) and (e). Maximizationof matched harmonicity may have the eﬀect of tracking stationary ﬂuctuations in harmonics that are notoften aﬀected by reverberation.
4.

Proposed method

4.1 Complex cepstrum analysis
Let us overview the results in Subsection 3.3 byreconsidering the complex cepstrum representation ofreverberant speech y(t). From Eqs. (9)-(11), the complex cepstrum of y(t) can be represented as
CY (q, τ ) = Csrc (q, τ ) + Cflt (q, τ ) + CH (q, τ )(16)where CH (q, τ ) is the complex cepstrum of the reverberant impulse response, h(t). These cepstra can also
be represented as all amplitude and phase cepstra (denoted by subscripts “A” and “φ”).
Complex cepstrum analysis, on the other hand, isusually used to separate minimum and non-minimum
phase characteristics. The complex cepstrum, C(q, τ ),can also be separately represented as
C(q, τ )=Cmin (ω, τ ) + Call (ω, τ )=CA,min (q, τ ) + Cφ,min (q, τ )+CA,all (q, τ ) + Cφ,all (q, τ ) (17)where the subscripts “min” and “all” indicate minimum and non-minimum phase characteristics. Here,respective spectra can be represented asX(ω, τ ) ==Xmin (ω, τ ) · Xall (ω, τ )|Xmin (ω, τ )| exp(jφmin (ω, τ ))×|Xall (ω, τ )| exp(jφall (ω, τ )) (18)where |Xall (ω, τ )| = 1 and CA,all (q, τ ) = 0. Hence,CY (q, τ ) can be separately represented as
CY,A,min(q, τ ) + CY,φ,min (q, τ ) + CY,φ,all (q, τ )= Csrc,A,min(q, τ ) + Csrc,φ,min (q, τ ) + Csrc,φ,all (q, τ )+Cflt,A,min (q, τ ) + Cflt,φ,min (q, τ ) + Cflt,φ,all (q, τ )+CH,A,min(q, τ ) + CH,φ,min (q, τ ) + CH,φ,all (q, τ )(19)Note that the amplitude cepstrum of all-pass phasecharacteristics has been omitted from this equation.
According to Eq. (16), an optimal F0 estimate is
only used to extract Csrc (q, τ ) from CY (q, τ ) to deal
with the periodicity/harmonicity of source information as a ﬁlter and the reverberation characteristics
are eliminated. It is too diﬃcult only to deal with
Csrc (q, τ ) in this estimation task, without measuring
h(t) or CH (q, τ ) (i.e., blind F0 estimation). In addition, long-term CH (q, τ ), in which the length of analysis is more than the reverberation time, is needed toaccurately extract Csrc (q, τ ).
We did a preliminary investigation into which
component, CH,min (q, τ ) or CH,all (q, τ ), most aﬀected
dealing with Csrc (q, τ ) for estimating F0 , using Eq.
(19). Figure 3 shows the process of estimating one of
the reverberant speech signals (/Tokushima-To-IebaAwa-Odori-Ga-Yuumei-Desu/, female speaker, reverberation time TR of 2.0 s) we used in the evaluations.
Speech signals (x(t) and reverberant y(t)) are shownin Figs. 3(a) and (b). The reference F0 (F0,Ref (t) by
TEMPO from the EGG signal) and the F0 (F0,Est (t))
estimated by the cepstrum method from y(t) correspond to the dashed and solid lines in Fig. 3(c).
Note that F0,Ref (t) obtained by TEMPO are completely within the voiced part because TEMPO has
an Unvoiced/Voiced (U/V) decision. As can be seen,
the estimated F0 was not close to the reference in thevoiced part. This method, however, can be used toaccurately estimate F0 from y(t) by eliminating theeﬀect of h(t) from y(t) on the complex cepstrum inthe long-term Fourier transform, as plotted in Fig.
3(d). At the same time, two comparative F0 s wereobtained as plotted in Figs. 3(e) and (f) by estimating F0 from y(t) by eliminating minimum phase or the
all-pass phase component from y(t).
The all-pass phase component of the reverberant
impulse response, h(t), we used appears to have adominant eﬀect from these comparisons of robust andaccurate F0 estimates. Although the same comparisons of all the other stimuli are not presented in this8

Journal of Signal Processing, Vol. , No. ,paper, the same trends were observed. Hence, we concluded that eliminating the all-pass phase characteristics of h(t) would enable eﬀective estimates of F0
from reverberant speech y(t). In addition, the cepstrum method with the all-pass component eliminatedraised the possibility of achieving robust and accurate
estimates of F0 since we knew homomorphic analysis
could easily deal with minimum phase characteristics
such as simpliﬁed echoes.

4

2

x 10

(a)

x(t)

1
0
−1
−2 x 1042(b)

4.2 Estimates of h(t) based on MTF concept

y(t)

1

The MTF concept was proposed by Houtgast andSteeneken [36] to account for the relation between thetransfer function of frequency in an enclosure in termsof the envelopes of input and output signals (x(t) andy(t)), and characteristics of the enclosure such as reverberation. This concept was introduced as a measure in room acoustics to assess what eﬀect enclosure
had on the intelligibility of speech [36]. The complex
modulation transfer function, m(ω), is deﬁned as

0
−1
−2
500
(c)

F0(t) (Hz)

400
300
200
100

m(ω) =

500∞
20 h(t) exp(jωt)dt∞
h(t)2 dt
0

(20)

(d)

F0(t) (Hz)400

where h(t) is the impulse response of the room acoustics and ω is the radian frequency. This equationmeans the Fourier transform of the squared impulse
response is divided by its total energy.
When reverberant impulse response h(t) as deﬁned
in Eq. (12) is substituted into the equation above, theMTF, m(ω), can be obtained as

300
200
100
500
(e)

F0(t) (Hz)

400




2 −1/2TRm(ω) = |m(ω)| = 1 + ω
13.8

300
200
100
500
(f)

F0(t) (Hz)

400
300
200
100
0

0.5

1

1.5

2

2.5

3

Time (s)Fig. 3 Examples: (a) original speech x(t), (b) reverberant speech y(t) (reverberation time of 2.0 s),(c) reference F0 using TEMPO from EGG of x(t) indicated by dashed-line and estimated F0 using cepstrum method from y(t) indicated by solid line, (d)estimated F0 from dereverberated y(t) using h−1 (t),(e) F̂0 from y(t) eliminated by minimum phase characteristics, and (f) F̂0 from y(t) eliminated by all-passphase characteristics(21)

where TR is the reverberation time, i.e., the time required for the power of h(t) to decay by 60 dB [36].
Figure 4 plots the MTF, m(ω), as a function of themodulation frequency, Fm , (i.e., the dominant frequency in the temporal envelope). These theoreticalcurves were calculated by substituting ﬁve reverberation times (TR = 0.1, 0.3, 0.5, 1.0, and 2.0 s) andω = 2πFm into Eq. (21). Here, m(ω) can also be
regarded as the modulation index with respect to Fm .
These curves reveal how much the modulation index
of the envelope will be reduced from 1.0 to 0.0 depending on TR at a speciﬁc Fm . In other words, TR can bepredicted from a speciﬁc m(2πFm ) at a speciﬁc Fm .
Therefore, the temporal envelope of the reverberant
impulse response, a exp(−6.9t/TR), can also be predicted using TR and the “a” in Eq. (13).
Based on the MTF concept, we can establish howmuch reverberation aﬀects a reduction in m(ω) and we
can then predict the characteristics of room acoustics(TR ) using inverse MTF. MTF-based power-envelope
inverse-ﬁltering methods, which have aimed at restoring the reduced MTF in the temporal envelope of the91
TR = 0.1 s

0.9
0.8

MTF, m(2π Fm)

0.7
T = 0.3 s

0.6

R

0.5
0.402

T = 0.5 s

0.4

R

0.3
T =1s
R

0.2

TR = 2 s

0.1
0

0

5

10
15
Modulation frequency, Fm (Hz)20

Fig. 4 Theoretical curves representing modulationtransfer function m(2πFm ) for various conditions with
TR = 0.1, 0.3, 0.5, 1.0, and 2.0 s
signals, have been proposed by the present authors[37, 38]. A technique of predicting TR using the inverse
MTF from the temporal envelope of observed signaly(t) has also been proposed [37, 38]; this was used as aconstraint for deriving the optimal restoration of thepower envelope for the original signal based on the deconvolution relationship between the power envelopesof y(t) and h(t). (For details, see Appendix):

 T


2
T̂R = max

arg minTR,min ≤TR ≤TR,max

0min êx,TR (t) , 0

dt

(22)
2

where T is the signal duration and êx,TR (t) is the setof candidates for the restored power envelope of clean
signal x(t) via inverse MTF [37] as a function of TR .
Note that the operation of “max(arg min{·})” meansthat the maximum argument of TR needs to be determined from a timing point where the negative areaof êx,TR (t)2 approximately equals zero or a particularminimum area. Here, TR,min and TR,max are the lower
and upper limited regions of TR [37].
According to Eqs. (12), (13) and (22), h(t) canbe estimated by utilizing â exp(−6.9t/TˆR) with simulated white noise n̂(t). This is referred to as ĥ(t). In
this case, long-term CH (q, τ ) can be directly obtained
from ĥ(t). Although this does not completely equal
the original h(t) we used in the evaluation, long-term
amplitude cepstrum CH,A (q, τ ) can only be matchedto the original. This is because the MTFs of h(t) andĥ(t) are the same if T̂R is a complete value, and Eqs.
(20) and (21) can be regarded as having characteristics of CH,A (q, τ ) that can be indirectly obtained fromF −1 [log |m(ω)|] with the power factor on the LTFT.
Therefore, it can be easily predicted that CH,A (q, τ )becomes a cepstral shape of exponential decay withrespect to quefrency (dominant at lower quefrencies).

In contrast, although it is diﬃcult to obtain a complete value with regard to phase cepstrum CH,φ (q, τ ),long-term CH,φ,all (q, τ ) can be estimated from themby using Eqs. (17) and (19). As explained in Sec.
4.1, using estimated CH,φ,all (q, τ ) from ĥ(t) to eliminate the all-pass phase component from reverberant
speech y(t) on the basis of LTFT should be done toestimate F0 . Although the estimated CH,A,min (q, τ )can also be canceled out in Eq. (19) on LTFT, theelimination of minimum-phase characteristics in Eq.
(19) on LTFT is not as eﬀective for eliminating allpass phase characteristics so that this has not beenused in this paper. The short-term CH,A,min (q, τ ) andCH,φ,min (q, τ ) to be canceled out in Eq. (19) on STFT
will be considered in the next section.
4.3 Liftering on complex cepstrum
CH,φ,all (q, τ ) is canceled out in Eq. (19) on LTFT
as explained in the previous section, so that the remaining terms are Cflt (q, τ ) and CH,min (q, τ ) to extract Csrc (q, τ ). Complex cepstrum analysis and thesource-ﬁlter model are used to cancel out the remaining terms in Eq. (19) on STFT to take the best advantage of homomorphic processing.
There is a Hilbert transform relationship between
CA,min (q, τ ) and Cφ,min (q, τ ), and the latter has thesame characteristics in the positive quefrency domain
based on the minimum phase characteristics. However, short-term CH,A,min (q, τ ) and CH,φ,min (q, τ ) are
not the same as the long-term versions when STFT
analysis is shorter than the reverberation time. However, amplitude cepstrum CH,min (q, τ ) in the lower
quefrency parts is generally larger than that in thehigher parts and this exponentially attenuates as thequefrency increases. Therefore, the minimum phasecharacteristics, CH,min (q, τ ), can be assumed to concentrate on the lower quefrency parts.
The cepstrum components of the source characteristics are separately concentrated on the higher quefrency parts and those of the ﬁlter are separately concentrated on the lower based on the advantages of thesource-ﬁlter model, as shown in Fig. 1. Therefore, if
a component on the low quefrency part can only beremoved by liftering, the ﬁlter characteristics as wellas the dominant components of the minimum-phasecharacteristics of reverberation can be canceled out in
Eq. (19). Thus, the following lifter, l(q), is used inthis paper to cancel them out in Eq. (19).

0, q ≤ qlif
(23)l(q) =
1, q > qlif
where qlif = 1.25 ms. This means the upper limit for
estimated F0 is 800 Hz.

10

Journal of Signal Processing, Vol. , No. ,Reverberant speech y(t)Long-term Fouriertransform (LTFT)

TR estimation

Y(ω)

TR

Complex cepstrum
analysis (CCA)

Estimation of MTF
based imp. response

proposed method treats a complex cepstrum, the restored short-term waveform s(t, τ ) from Csrc (q, τ ) canbe used to estimate F0 with the autocorrelation function and/or AMDF. As the aim of this paper was topropose a model concept for robustly estimating F0
in reverberant environments, these kinds of considerations with regard to modiﬁcations in processing are
beyond the scope of this paper.

CY(q)Elimination CH,φ,all(q)from CY(q)LTFT & CCA
CH(q)Inversed LTFTy(t)Sort-term Fouriertransform (STFT)

CCA

Liftering based on
source-filter modell(q)Csrc(q,t)5.

Evaluation of proposed method

5.1 Method
We evaluated the proposed method with (labeled “Proposed(Est)”) and without (labeled “Proposed(Org)”) TR estimates by using the same procedure and sound dataset described in Section 3. Withand without comparisons of the proposed method
were done to ﬁnd out how accurate the TR estimates
were. We compared them with TEMPO, the cepstrum
method, and a modiﬁed complex cepstrum method
based on the source-ﬁlter model (labeled “SrcFlt”).
The SrcFlt method was used to ﬁnd out how eﬀectively CH,φ,all (q, τ ) was eliminated on the LTFT withthe proposed method.

FFT

5.2 Results and discussion
logXsrc(ω,t)

Comb filtering

F0 estimationF0(t)Fundamental frequency

Fig. 5

Algorithm for proposed method

4.4 Proposed method of estimating F0
The algorithm for estimating F0 based on complex
cepstrum analysis, the MTF concept, and the sourceﬁlter model are explained by Fig. 5. This method iscomposed of three main processes: (1) estimating the
MTF-based reverberation impulse responses and eliminating them from reverberant speech, (2) extracting
Xsrc (ω, τ ) from the processed reverberant speech byusing liftering on the complex cepstrum based on thesource-ﬁlter model, and (3) estimating F0 from themby using a ﬁnal decision block.
Comb ﬁltering was employed in the ﬁnal two blocks
in Fig. 5. As these are commonly used in classical methods of estimation, such as in comb-ﬁlteringand autocorrelation functions, they can be replacedby the autocorrelation function. In addition, since the

Figure 6 plots the results for the comparative evaluations. The correct rates within error margins of 5 %and 10 % for the proposed and the other methods areplotted in Figs. 6(a) and (b). Their SNRs are plottedin Fig. 6(c). The results for the cepstrum method indicate the baselines in the evaluations while those for
TEMPO (dashed-line) indicate the lower limits.
Although the overall accuracy of F0 estimatestended to be reduced as reverberation time increased,
about a 10 % improvement in the correct rates andabout a 5 dB improvement in the SNR could be obtained with the new method. There are fewer diﬀerences in the results for both the proposed methods
with and without TR estimates. This means the TR
estimates can work as well. Since a correct rate of 60% within an error margin of 5 %, a correct rate of
75 % within an error margin of 10 %, and an SNR of
17 dB at TR = 2.0 s were achieved with the method
we propose, we concluded that MTF-based impulse
responses can be precisely estimated by utilizing TR
estimates. For example, the results for extracting F0
at TR = 2.0 with the proposed method, with and without TR estimates, from the same reverberant speech(Fig. 3(b)) are plotted in Figs. 6(d) and (e).
The SrcFlt method results indicate a slight improvement (about 3 % in the correct rate) to that withthe cepstrum method. In contrast, there were about
7 % and 5 dB improvements in the percent correct
rate and in the SNR by using the new method. We

11concluded that the use of complex cepstrum analysis with regard to non-minimum phase characteristics
eﬀectively enabled F0 to be estimated in reverberant
environments.
As shown in Figs. 3(b), 6(d), and 6(e), most F0 estimation methods including the proposed method,which does not have the U/V decision output of redundant information in the unvoiced and silent sections, while TEMPO does. An important issue is how
to determine the U/V decision rules for applicationsof speech-signal processing in our next stage of research. This issue should be able to be resolved byincorporating power discrimination such as that usedin TEMPO with the U/V decision rules, but this isbeyond the scope of this paper.
6.

Conclusion and Future Perspectives

We evaluated the robustness and accuracy of
twelve typical methods of estimating F0 (i.e., classic ACMWL, AMDF, STFT-based, cepstrum, LPC,
and SHS algorithms, and modern IFHC, PHIA, andTEMPO algorithms) in artiﬁcial reverberant environments using huge speech datasets. The results revealed that none of these methods could accurately estimate F0 in reverberant environments and that theiraccuracies drastically decreased as reverberation timeincreased. The results also demonstrated that the best
method was cepstrum-based and that the worst was
the instantaneous frequency-based model. We foundthat periodicity and/or harmonicity on the complex
cepstrum eﬀectively enabled F0 to be estimated in reverberant environments.
We proposed a robust and accurate method of
estimating F0 that was based on the source-ﬁlter
model concept and the MTF concept in complex cepstrum analysis. This method included (1) eliminating the dominant reverberation eﬀect from observed
speech by estimating MTF-based reverberant impulse
responses and (2) extracting source information fromthem by subtracting the remaining cepstrum related
to ﬁlter characteristics and the remaining reverberation through liftering. We demonstrated that our new
method is robust against reverberation and can accurately estimate F0 from observed reverberant speech,
using the same comparative evaluations.
Additional improvements may be possible by modifying the F0 determination block. Further evaluations using real reverberant impulse responses in room
acoustics are required for real applications. Extensional improvements may be possible by incorporating
the proposed method into the U/V decision rules andby considering the robustness of the new method inboth noisy and reverberant environments, but theseare beyond the scope of this paper.
In future work, we hope to incorporate our new
robust and accurate approach to estimate F0 into the

MTF-based process of speech dereverberation [39, 40]
and then to improve them so that they can become amore complete blind-dereverberation method. Thisis because the current method, which consists of
MTF-based power-envelope inverse ﬁltering and carrier restoration using F0 information can dereverberate reverberant speech in power envelopes as well as inresynthesized waveforms, assuming that F0 can be accurately estimated. As mentioned in the Introduction,robust and accurate F0 is a signiﬁcant feature of our
new method so that it should be able to contributeto resolving the problems we previously experienced
with speech dereverberation [39, 40].
7.

Acknowledgments

This work was supported by a Grant-in-Aid for
Scientiﬁc Research from the Ministry of Education,
Culture, Sports, Science, and Technology of Japan(No. 18680017). It was also partially supportedby the Strategic Information and COmmunications
R&D Promotion ProgrammE (SCOPE) (071705001)
of the Ministry of Internal Aﬀairs and Communications (MIC), Japan.

References[1] W. J. Hess: Pitch Determination of Speech Signals,
Springer-Verlag, New York, 1983.
[2] W. J. Hess: Pitch and Voicing Determination in Advancesin speech signal processing, Eds. S. Furui and M. M.
Sondhi, pp. 3–48, Marcel Dekker. Inc. New York, 1992.
[3] A. de Cheveigné and H. Kawahara: Comparative evaluation of F0 estimation algorithms, Proc. Eurospeech2001,
pp. 2451–2454, Sep. 2001.
[4] B. Gold and L. Rabiner: Parallel processing techniques
for estimating pitch periods of speech in the time domain,
J. Acoust. Soc. Am., Vol. 46, No. 2, pp. 442–448, Aug.
1969.
[5] N. C. Geckinli and D. Yavuz: Algorithm for pitch extraction using zero-crossing interval sequence, IEEE Trans.
Acoustics, Speech, and Signal processing, Vol. ASSP-25,
No. 6, pp. 559–564, Dec. 1977.
[6] M. R. Schroeder: Period histogram and product spectrum:
new methods for fundamental frequency measurement, J.
Acoust. Soc. Am., Vol. 43, No. 4, pp. 829–834, Apr. 1968.
[7] D. M. Howard: Peak-picking fundamental period estimation for hearing prostheses, J. Acoust. Soc. Am., Vol. 86,
No. 3, pp. 902–910, Sep. 1989.
[8] M. M. Sondhi: New methods of pitch extraction, IEEE
Trans. Audio and Electroacoustics, Vol. AU-16, No. 2, pp.
262–266, Jun. 1968.
[9] M. J. Ross, H. L. Shaﬀer, A. Cohen, R. Freudberg and H.
J. Manley: Average magnitude diﬀerence function pitch
extractor, IEEE Trans. Acoust., Speech, Signal Processing, Vol. ASSP-22, No. 5, pp. 353–361, Oct. 1974.

12

Journal of Signal Processing, Vol. , No. ,

90

90

80

80

70

70

Correct rate (%)

100

Correct rate (%)

100

60
50
40
30

50
40
30

TEMPO
Cepstrum
SrcFltProposed(Org)
Proposed(Est)

20
10
0

60

0

0.5

TEMPO
Cepstrum
SrcFltProposed(Org)
Proposed(Est)

20

(a)
10

5 % error margin
1

T (s)

1.5

0

2

0

0.5

(b)
10 % error margin
1

T (s)

R

1.5

2

R

30
500

TEMPO
Cepstrum
SrcFltProposed(Org)
Proposed(Est)

(d)
400

F0(t) (Hz)25

SNR (dB)

20

300
200
100

15
500

(e)
400

5

0

F0(t) (Hz)

10

(c)

300
200
100

0

0.5

1

Reverberation time TR (s)

1.5

2

0

0.5

1

1.5
2
Time (s)2.5

3Fig. 6 Evaluation results: (a) percent correct rate within error margin of 5%, (b) percent correct rate withinerror margin of 10%, (c) SNR of F0 estimates from reverberant speech using proposed method, and examples ofextracted F0 using proposed model (d) without TR estimation and (e) with TR estimation[10] J. D. Wise, J. R. Caprio and T. W. Parks: Maximum likelihood pitch estimation, IEEE Trans. Acoustics, Speech,
Signal Processing, Vol. ASSP-24, No. 5, pp. 418–423, Oct.
1976.
[11] K. Nishi and S. Ando: An optimal comb ﬁlter for timevarying harmonics extraction, IEICE Trans. Fundamentals, Vol. E81-A, No. 8, pp. 1622–1627, Aug. 1998.
[12] T. Miwa, Y. Tadokoro and T. Saito: The pitch estimationof diﬀerent musical instruments sounds using comb ﬁlters
for transcription, IEICE, Trans. D-II, vol. J81-D-II, no. 9,
pp. 1965–1974, Sep. 1998.
[13] T. Shimamura and H. Kobayashi: Weighted autocorrelation pitch extraction of noisy speech, IEEE Trans. Speech
and Audio Processing, Vol. 9, No. 7, pp. 727–730, Oct.
2001.
[14] D. J. Hermes: Measurement of pitch by subharmonic summation, J. Acoust. Soc. Am., Vol. 83, No. 1, pp. 257–264,
Jan. 1988.

[15] A. M. Noll: Ceptrum pitch determination, J. Acoust. Soc.
Am., Vol. 41, No. 2, pp. 293–309, Aug. 1966.
[16] A. M. Noll: Clipstrum pitch determination, J. Acoust.
Soc. Am., Vol. 44, No. 6, pp. 1585–1591, Jul. 1968.
[17] L. R. Rabiner, M. J. Cheng, A. E. Rosenberg and A. McGonegal: A comparative study of several pitch detection
algorithms, IEEE Trans. Acoustics, Speech, Signal Processing, Vol. ASSP-24, pp. 399–413, 1976.
[18] C. K. Un and S. C. Yang: A pitch extraction algorithm
based on LPC inverse ﬁltering and AMDF, IEEE Trans.
Acoust., Speech, Signal Process. Vol. ASSP-25, No. 6, pp.
565–572, Dec. 1977.
[19] T. V. Ananthapadmanabha and B. Yegnanarayana:
Epoch extraction from linear prediction residual for indentiﬁcation of closed glottis interval, IEEE Trans. Acoustics,
Speech, Signal Processing, Vol. ASSP-27, No. 4, pp. 309–
319, Aug. 1979.
[20] H. Singer and S. Sagayama: Pitch dependent phone mod-

13

(a)

(b)
4
2
x(t)

x

e (t)2

1
0.5

0
−2

0
−3

(c)

x 10

(d)
0.2
TR =0.5 (s)

h(t)

h

e (t)2

4
2

0
−0.2

0
(e)

(f)
4
2
y(t)

y

e (t)2

1
0.5

0
−2

0

0

(g)
2

x

^ e (t)2

eling for HMM-based speech recognition, J. Acoust. Soc.
Jpn. (E), Vol 15, No. 2, pp. 77–86, Mar. 1994.
[21] J. D. Markel: The SIFT algorithm for fundamental frequency estimation, IEEE Trans. Audio, Vol. AU-20, No.
5, pp. 367–377, Dec. 1972.
[22] N. Kunieda, T. Shimamura and J. Suzuki: Pitch extraction by using autocorrelation function on the log spectrum, IEICE Trans. A, Vol. J80-A, No. 3, pp. 435–443,
Mar. 1997.
[23] H. Kobayashi and T. Shimamura: An extraction methodof fundamental frequency using clipping and band limitation on log spectrum, IEICE Trans. A, Vol. J82-A, No. 7,
pp. 1115–1122, Jul. 1999.
[24] T. Takagi, N. Seiyama and E. Miyasaka: A Method forpitch extraction of speech signal using autocorrelation
functions through multiple window-length, IEICE Trans.
A, Vol. J80-A, No. 9, pp. 1341–1350, Sep. 1997.
[25] H. Kawahara, H. Katayose, A. de Cheveigné and R. D.
Patterson: Fixed Point analysis of frequency to instantaneous frequency mapping for accurate estimation of F0
and periodicity, Proc.Eurospeech99, No. 6, pp. 2781–2784,
Sep. 1999.
[26] A. de Cheveigné and H. Kawahara: Yin, a fundamental
frequency estimator for speech and music, J. Acoust. Soc.
Am., Vol. 111, No. 4, pp. 1917–1930, Apr. 2002.
[27] H. Kawahara, I. Masuda-Katsuse and A. de Cheveigné:Restructuring speech representations using a pitchadaptive time-frequency smoothing and an instantaneousfrequency-based F0 extraction: Possible role of a repetitive structure in sounds, Speech Communication, Vol. 27,
pp. 187–207, Apr. 1999.
[28] Y. Atake, T. Irino, H. Kawahara, J. Lu, S. Nakamura
and K. Shikano: Robust fundamental frequency estimation using instantaneous frequencies of harmonic components, Proc of ICSLP2000, Vol. 2, pp. 907–910, Oct. 2000.
[29] Y. Ishimoto, M. Unoki and M. Akagi: A Fundamental
Frequency Estimation Method for Noisy Speech Basedon Instantaneous Amplitude and Frequency, Proc. EuroSpeech2001, pp. 2439–2442, Sep. 2001.
[30] M. Unoki and M. Akagi: A Method of Signal Extraction from Noisy Signal based on Auditory Scene Analysis,
Speech Communication, Vol. 27, No. 3, pp. 261–279, Apr.
1999.
[31] T. Nakatani and T. Irino: Robust and accurate fundamental frequency estimation based on dominant harmoniccomponents, J. Acoust. Soc. Am. Vol. 116, No. 6, pp.
3690–3700, Dec. 2004.
[32] P. P. Vaidyanathan: Multirate systems and Filter Banks,
Prentice-Hall, New Jersey, 1993.
[33] H. Ohmura and K. Tanaka: Fine pitch contour extraction
by voice fundamental wave ﬁltering method, J. Acoust.
Soc. Jpn, Vol. 51, No. 7, pp. 509–518, Jul. 1995.
[34] M. R. Schroeder: Modulation transfer functions: deﬁnition and measurement, Acustica, Vol. 49, pp. 179–182,
1981.
[35] H. Kuttruﬀ: Room Acoustics, Taylor & Francis, fourth
edition, London, 2000.
[36] T. Houtgast and H. J. M. Steeneken: The modulationtransfer function in room acoustics as a predictor of speech
intelligibility, Acustica, Vol. 28, pp. 66–73, (1973).

0.2

0.4

0.6

T =0.3 (s)R
TR=0.5 (s)TR=1.0 (s)

1
0
0

0.2

0.4

0.6

time (s)Fig. 7 Examples of relationship between power envelope of system based on MTF concept: (a) powerenvelope e2x (t) of (b) original signal x(t), (c) powerenvelope e2h (t) of (d) impulse response h(t), (e) powerenvelope e2y (t) derived from e2x (t) ∗ e2h (t), (f) reverberant signal y(t) derived from x(t)∗h(t), and (g) restored
power envelope ê2x (t)[37] M. Unoki, M. Furukawa, K. Sakata and M. Akagi: An
improved method based on the MTF concept for restoring
the power envelope from a reverberant signal, Acoustical
Science and Technology, Vol. 25, No. 4, pp. 232–242, Apr.
2004.
[38] M. Unoki, K. Sakata, M. Furukawa, and M. Akagi: A
speech dereverberation method based on the MTF concept in power envelope, Acoustical Science and Technology, Vol. 25, No. 4, pp. 243–254, Apr. 2004.
[39] M. Unoki, K. Sakata and M. Akagi: A speech dereverberation method based on the MTF concept, Proc. EuroSpeech2003, pp. 1417–1420, Sep. 2003.
[40] M. Unoki, M. Toi, and M. Akagi: Development of the
MTF-based speech dereverberation method using adaptive time-frequency division, Proc. Forum Acusticum
2005, pp. 51–56, Aug. 2005.

Appendix MTF-based power envelope inverse
filtering
In the model of inverse-ﬁltering of the power envelope [37], the observed reverberant signal, the original
signal, and the stochastic-idealized impulse response
in room acoustics are denoted as bf y(t), x(t), andh(t), respectively, and are modeled as:
y(t)
x(t)

= x(t) ∗ h(t)= ex (t)n1 (t)h(t)= eh (t)n2 (t) = a exp(−6.9t/TR)n2 (t)(26)

(24)
(25)

where the asterisks “*” denote the convolution operation, ex (t) and eh (t) are the envelopes of x(t) and14h(t), and n1 (t) and n2 (t) are the mutually independent respective white noise functions, i.e.,< nk (t), nk (t − τ ) >= δ(τ ). The parameters of theimpulse response, a and TR , correspond to a constant
amplitude term and the reverberation time.
In this model, the power envelope of the reverberant signal, ey (t)2 , can be determined as
 2 
(27)
y (t) = e2x (t) ∗ e2h (t) = e2y (t)where < · > is the ensemble average operation. Thisequation shows that there is a signiﬁcant relationshipbetween the envelopes; i.e., the MTF concept. Basedon this result, e2x (t) can be recovered by deconvoluting
e2y (t) with e2h (t). Here, the transmission functions of
power envelopes Ex (z), Eh (z), and Ey (z) are assumed
to correspond to the z-transforms of e2x (t), e2h (t), ande2y (t). Thus, the transmission function of the powerenvelope of the original signal, Ex (z), can be determined fromEy (z)
13.8−1
Ex (z) =
1
−
exp
−
z
(28)
a2
T R · fs
where fs is the sampling frequency. Finally, the powerenvelope, e2x (t), can be obtained from the inverse ztransform of Ex (z).
Figure 7 shows an example of how power-envelope
inverse ﬁltering is related to the MTF concept. Figure 7(a) shows a sinusoidal power envelope as original power envelope e2x (t) (= 0.5(1 + sin(2πFm t)); themodulation frequency, Fm , was 10 Hz and the modulation index, m, was 1.0). Figure 7(b) shows theoriginal signal, x(t), calculated from e2x (t) and a white
noise carrier, n1 (t), using Eq. (25). Figure 7(c) shows
power envelope e2h (t) calculated using Eq. (26) with
TR = 0.5 s. Figure 7(d) shows the impulse response,h(t), of Eq. (12), calculated from e2h (t) and a white
noise carrier, n2 (t). Figures 7(e) and (f) show thepower envelope, e2y (t), obtained from a convolution ofe2x (t) ∗ e2h (t) and the observed reverberant signal, y(t),obtained from a convolution of x(t) with h(t), respectively. The left panels ((a), (c), and (e)) show thepower envelopes of the signals and the right panels((b), (d), and (f)) show the corresponding signals. In
this ﬁgure, the modulation index decreased from 1.0(in Fig. 7(a)) to 0.404 (maximum deviation of envelope between the dotted lines in Fig. 7(e) relative tothat in Fig. 7(a)). Since the MTF concept shows themodulation index as a function of Fm and TR , it canalso be shown that the decreased modulation index is
derived from m(2πFm ) = 0.402 using Eq. (21) by substituting TR = 0.5 s and Fm = 10 Hz into Eq. (21).
The solid line in Fig. 7(g) indicates the restored powerenvelope ê2x (t) obtained from reverberant power envelope e2y (t) (Fig. 7(e)) using Eq. (28) with TR = 0.5 s.
We can see that power-envelope inverse ﬁltering canaccurately restore the power envelope from a reverberant signal in terms of shape and magnitude.

Journal of Signal Processing, Vol. , No. ,Masashi Unoki
was born in
Akita Pref., Japan, in 1969. Hereceived his M.S. and Ph.D. (Information Science) from the Japan
Advanced Institute of Science andTechnology (JAIST) in 1996 and1999. His main research interests
are auditory-motivated signal processing and the modeling of auditory
systems. He was a JSPS research
fellow from 1998 to 2001. He wasassociated with the ATR Human Information Processing Laboratories as a visiting researcher during 1999–2000, and from2000 to 2001 he was a visiting research associate at CNBH in
the Dept. of Physiology at the University of Cambridge. He
has been on the faculty of the School of Information Science at
JAIST since 2001 and he is now an associate professor. He is amember of the Research Institute of Signal Processing (RISP),
the Institute of Electrical and Electronic Engineering (IEEE),
the Institute of Electronics, Information and Communication
Engineers (IEICE) of Japan, the Acoustical Society of America(ASA), the Acoustical Society of Japan (ASJ), and the International Speech Communication Association (ISCA). Dr. Unoki
received the Sato Prize for an Outstanding Paper from the ASJ
in 1999 and the Yamashita Taro Prize for Young Researcher
from the Yamashita Taro Research Foundation in 2005.

Toshihiro Hosorogiya
wasborn in Ishikawa Pref., Japan in
1980. He received his B.E. from
Nagoya University in 2005, and his
M.S. from the Japan Advanced Institute of Science and Technology in2007. He is a member of the Research Institute of Signal Processing(RISP) and the Acoustical Societyof Japan (ASJ).
(Received June 17, 2007; revised September 11, 2007)