2 36

IEEE TRANSACTIONS O N ACOUSTICS, SPEECH, A N D SIGNAL PROCESSING, VOL. ASSP-32, NO. 2 , APRIL 1984

Signal Estimation from Modified Short-Time
Fourier Transform
DANIEL W. GRIFFIN A N D JAE S. LIM, SENIOR MEMBER, IEEE

Abstract-In this paper, we present an algorithm to estimate a signal
from its modified short-time Fourier transform (STFT). This algorithm
is computationally simple and is obtained by minimizing the meansquared error between the STFT of the estimated signal and the modified STFT. Using this algorithm, we also develop an iterative algorithm
to estimate a signal from its modified STFT magnitude. The iterativealgorithm is shown to decrease, in each iteration, the mean squared
error betweentheSTFT
magnitude of the estimated signal and themodified STFT magnitude. Themajor computation involved in theiterative algorithm is the discrete Fourier transform (DFT)computation, and the algorithm appears to be real-time implementable withcurrenthardware technology. The algorithm developed in this paper
has been applied to the time-scale modification of speech. The resulting system generates very high-quality speech, and appears to be betterin performancc than any existing method.

mated signal and the MSTFT. The resulting algorithm is quite
simple computationally. In Section 111, the algorithminSection I1 is used t o develop an iterative algorithm that estimatesa signal from the MSTFTM. The iterative algorithm is shownt o decrease, in each iteration, the mean squared error between
the STFTM of the estimated signal and the MSTFTM. In Section IV, we present an example of the successful application
ofourtheoreticalresults.Specifically,
we develop a timescale speechmodificationsystembymodifyingtheSTFTM
first and then estimating a signal from the MSTFTM using thealgorithmdeveloped in Section 111. The resultingsystemhas
beendemonstrated to gene]-ateveryhighquality,time-scale
modified speech.
11. SIGNALESTlblATION

I. INTRODUCTIONFROM MODIFIEDSHORT-TIMEFOUI<I~.:KTRANSFOKbl

I

N a number of practical applications [1]-[5], it is desirable
Let x ( n ) and X,(nzS, w ) denotea real sequenceandits
to modify the short-time Fourier transform (STFT) or the
short-time Fourier transform magnitude (STFTM) and then es- STFT. The variable S is a positiveinteger,whichrepresents
X,(n, w ) inthevariable n. Letthe
timate the processed signal from the modified STFT (MSTFT) thesamplingperiodof,
or the modified STFTM (MSTFTM). For example, in speech analysis window used in the STFT be denoted by ~ ( n )andenhancement by spectralsubtraction [ 2 ] , [3] theSTFT is withlittle loss ofgenerality, w ( n ) is assumed t o bereal, I,
modified by combining the STFT phase of the degraded speechpoints long, and nonzero for 0 < n < I, - 1. From the definition of the STFT
with a MSTFTM, and then a signal is reconstructed from the
MSTFT. As anotherexample,inthetime-scalemodificationmof speech, one approach is to modify the STFTM and then tox,(~?zs,o)= F~[x,(urz~,I ) ] =x,(rn~, I ) d m '(1)[=reconstruct a signal from the MSTFTM. In most applications,
including the two cited above, the MSTFT or MSTFTM is not
wherevalid in the sense that no signal has the MSTFT or MSTFTM,
x,(rnS, I ) = w(mS - I ) x(/)(2)
and therefore it is important to develop algorithms to estimate
a signal whose STFT or STFTM is close in some sense to theand Fl [xw(mS,I ) ] represents the Fourier transform ofx,(mS,to thisproblem
MSTFTor MSTFTM.Previousapproaches
I ) with respect to the variable 1.
have been mostly heuristic [6] -[8], and have been limited to
Let Y,(mS, w ) denote the given MSTFT and let y,(mS, I )estimatinga signal fromtheMSTFT [ 6 ] , 171. Inthispaper,
be given byto
we developnewalgorithmsbasedontheoreticalgrounds
estimate a signal from the MSTFT or the MSTFTM.
In addiY,(mS, w ) c j w l dw.
(3)tion, the new algorithm is applied to the problem of time-scale
Y w ( m S ,I ) = 271 w = - r
modificationofspeech.Theresultingsystem
is considerably
simpler conceptually and appears t o have better performance An arbitrary Y , ( m S . o),in general, is not a valid STFT in thethan the system described by Portnoff [ I ] .
sense that there is no sequence whose STFT isgiven by Y,(mS,The paper is organized as follows. In Section 11, we develop 0 ) . In this section, we develop a new algorithm to estimate aanalgorithm to estimate a signal fromtheMSTFT by mini- sequence X(/?) whose STFT X , ( m S , w ) is closest to Y,(rn.S,of the esti- o)in the squared error sense.
mizing the mean squared error between the STFT
~

--

sn

Consider the following distance measure between
Manuscript received Dcccnlber 27, 1982; revised May 1 2 , 1983, andSeptember 26, 1983. This {vork was supported in part by the Advanced
Research Projects Agency monitored b y ONR under Contract SO00 1481-K-0742NR-049-509 and the National Scicncc I:ound;~tionunderG r a n t I<CS80-07102.
The authors arc with thc Iicscarch Laboratory of Iblcctronics. Dcpartmcnt of I~.lcctricalIknpinccring a n d ComputcrScicncc, Massachusetts
Institute ol'Tcchnology, Cambridge, M A 021 39.

a given MSTFT Y,(mS,ea

D [ X ( f Z ) ,Y,(f72S, a)] =OOc)6-35I8/84/04OO-O236$01 .OO 0 1984 I E E E

x(n) and

0):

I
2n

--

m = -m-Y,(mS,

In/X,(flZS, w )

w =-n

o)l2

dw.

(4)

GNAL
ON

LIM:
GRIFFIN AND

FROM MODIFIED STFT

237

Thedistancemeasurein(4) is thesquarederrorbetween
X,(mS, o)and Y,(mS, w ) integrated over all w and summed
over all m. I t has been written as a function of x(n) and Y ,(mS, a)to emphasize that X,(mS, o)is a valid STFT while
Y,(mS, o)is not necessarilyavalid
STFT. By Parseval’s
theorem, (4) can be written as= - _ wz(mS the window in (6) can be normalized so thatn) is unity for all n. Any nonzero window can be normalized
inthismannerformaximumwindowoverlap( S = 1). For

D [ x ( n ) , Y,(mS,all =2 2 [x,(mS, 0

m =--- y,(mS,

012

f

(’1

partial window overlap, however, the window ismore restricted.
Several windows which have this property for partial window
overlap are discussed below.
When thewindowshift ( S ) divides thewindowlength (I,)
evenly, the rectangular
windowdefinedby

O<n<L
w,(n) =

Since (5) is in the quadratic form of~ ( n )minimization,ofotherwiseD[x(n), Y,(mS, w ) ] canbeaccomplishedbysettingthe.
solving for x(n) hastheproperty
gradientwithrespectt o ~ ( n to) zeroand
which leads to the following result:(LIa-1 S
w,”(mS-n)=- = 1.
(1 0 )
w(mS - n)y,(mS, n )??I = -mm = O I,m =-m
( 6 ) We canfurthershowwithsomealgebrathat
if thewindow
x(n) =
length ( L ) is a multiple of four times the window shift( S ) thenw2(mS- n )m =-mbydefinedwindowsinusoidal the

2

5

5

This solution is similar in form to the standard overlap-add
procedure [6], [ 7 ] , ortheweightedoverlap-addprocedure[ 9 ] , [ I O ] . The overlap-add procedure can be expressed as

has the property given by (10). In addition, we require that
this class of sinusoidal windows be symmetric so that w(n) =w(L - 1 - n). Thisrequirementcanbesatisfiedbychoosing
qb = v/L. By choosing values for a and b , windows similar totheHammingwindowandtheHanningwindowcanbeobtained. Thus, the modified Hamming window used for timeIV will be defined as
scale modification of speech in Section
The weighted overlap-add procedure can be expressed as(1 1) for a = 0.54, b = - 0.46, and qb = v/L. The major differencebetweenthisdefinitionandthestandarddefinitionof
the Hamming window is that the period of the sine wave is I
x(n) =f ( m S - n)v,(mS, n )m=in the modified Hamming window as opposed to I - 1 for the
standardHammingwindow.Similarly,amodifiedHanningforsome“synthesis”filterf ( n ) . Themajordifferencebetween (6) and (7) is that (6) specifies that y w ( m S ,n ) should window can be defined as (1 1) for a = 0.5, b = -0.5, and qb =n/L. Use of these modified windows eliminates the need forbewindowedwiththeanalysiswindowbeforebeingover=-mw 2 (mS - n ) in ( 6 ) : which
reduces
lap added and w(mS - n ) should be squared before summation normalizingby
computation and/or storage requirements for partial window
over the variable m for normalization. The difference between(6) and (8) is that (6) explicitly specifies what f ( n ) is and has overlap.
Estimating x(n) based on (6) minimizesthesquarederror
the normalization constant. In addition, the major difference
between ( 6 ) , and (7) and ( 8 ) , is that (6) was theoretically de- between X,(mS, o)and Y,(rnS, a),andthereforecanbe
As will
rived explicitly for the purpose of estimating a signal from the used directly to estimate a sequence from a MSTFT.
MSTFT based on the least squares error criterion of (4). Equa- bediscussed in the next section, (6) canalsobe used to designal fromthetions (7) and (8), however, were derived to reconstruct a signal velopaniterativealgorithmthatestimatesa
MSTFTM.
to estimate a signal from the MSTFT
from its exact STFT or
for a very restricted class of modifications, and were sometimes
111. S I G N A L
ESTIMATIONm O M M O D I F I I ~used as ad hoc methods to estimate a signal from the MSTFT.
STFTMAGNITUDE
From the computational point ofview, the differences cited
In this section, we consider the problem of estimating x ( n )above are minor in terms of both the number of arithmetic opfrom the modified STFT magnitude
Y,(nzS, o ) l . The algoerations and the amount of on-line-storage required. For exrithm we develop is an iterative procedure based on the LSEEample, (6) can be implemented with little on-line storage andis similar in styleto several otherdelay, in the same manner [ I O ] as the standard-overlap proce- MSTFTalgorithmwhichof (8). iterative algorithms [ 1 I 1 , [ 121. In this algorithm, the squared
dureof (7) ortheweightedoverlap-addprocedure
and Y,(mS, o) is decreased inSince the algorithm represented by (6) minimizes the distance error between IX,(mS, o)/measure of (4), it will be referred to as LSEE-MSTFT, meaning each iteration. Let x’(n) denote the estimated x(n) after thei f 1st estimate x‘+’(n) is obtained by takith iteration. The
least squares error estimation from the MSTFT.
In the standard overlap-add method, the window
is usually ing the STFT of x‘(n), replacing the magnitude of Xk(nzS, w )mnormalized so that X,,=-, w ( m S - n ) is unityfor all n in with the given magnitude 1 Y,(nzS, w)l and then finding theorder to reduce computation. As in the overlap-add method, signal with STFT closest to this modified STFT using (6). The

2

-.X

~

I

1

238

IEEE
TRANSACTIONSGiven 1 Y,(mS,O)

O N ACOUSTICS, SPEECH, A N D SIGNAL PROCESSING, VOL. A S S - 3 2 , NO. 2, APRIL 1984

?L(rnS?

1

wherew ) is given by (13). Thisalgorithm will becalledOA(overlap-add)-MSTFTM to distinguish it fromthe
LSEE-MSTFTMalgorithm.Although
OA-MSTFTM requires
less windowing
fewermultiplicationsperiterationsinceone
step is required, it is not guaranteed to converge to the critical paints of D M . As willbe shown in Section IV, however,OA-MSTFTM does appear to reduce DM enough t o produce areasonable signal estimate for the purposes of time-scale modification of speech.

Initial Estimate of x(n)

>

J

x'(n)

One method of decomposing a speech signaly(n) is to represent it as the convolution of an excitation function with the
vocal tract impulse response. Consequently, the STFT
magniY,(mS, w ) can be written as thetude of this speech signal
IP,product of a component due to the excitation function(nzS, w)l andacomponentduetothevocaltractimpulse
Thisdecomposition is valid if the
response IH,(mS, o)/.
1,'ig. 1. LSI1E-MSTFTM algorithm.
analysis window is long enough to include several vocal tract
impulse responses and short enoughso that the speech signal
iterative algorithm, which is illustrated in Fig. 1. results in the is approximatelystationaryoverthewindowlength.Under
following update equation:/P,(mS, w)l will correspondtheseconditions,thefunction
to the rapidly varying portion of IY,(mS, w ) / with w ?taking
on an harmonic structure for voiced speech or noise for unIH,(mS, w)l will correspond tovoiced speech. The function
the slowly varying portion ofY,(mS, w)I with w , and will
LV2 ( m s - n )includetheformant
informationofthespeechsignal. Sincem = -mthespeech signal is assumed to beapproximatelystationary
over the window length, IP,(mS, w)l and IH,(mS, w)l willchange slowly with the time index mS as the pitch pcriod andwhere
vocal tract impulse response change.
The goal of time-scale modification is to modify the rate at
which /P,(mS, w)l and /H,(mS, w ) ' vary withtime,and*.
hence the rate a t which Y,(mS, w)I varies with time, without
In (13): if IX:v(mS, w)1 = 0, then X k ( m S , w ) is set t o Y,
affecting the spectral characteristics. This can be accomplished(nzS, a)/.
It can be shown (see Appendix) that the algorithm
by estimating a signal with STFT magnitude close to a tirnein Fig. 1 decreases in eachiteration the followingdistancescale modified version of Y,(mS, w)i. A time-scale modificameasure:tion of SI:S2 can be performed by calculating Y,(mS,, w)Ist
the window shift SI and XL>(mS2,o)at the window shiftlT
4RI[x(n), /Yw(mS,C3)II =[ : X w ( m S , w ) I 5'2 in the LSEE-MSTFTM or OA-MSTFTM algorithms.Form
2n w =-71"line up atthe
example, Y , . ( ~ s , ,w j l forthesentcnce
Y,(nzS, w ) l ] %d o .
(14) screcn door." saclpled at I O kHz is shown in Fig. 2 for a 256
pointmodifiedHammingwindowand
a windowshift S1 of
It can also be shown (see Appendix) that the algorithm always
128. Fig. 3(a) shows a
128 : 64 time-scale modified version of
converges to a set consisting of the critical points of the disY,(/nS,, o ) l produced by displaying these samples of ' YW(n,tance measure D,v as a function of x ( M ) . This algorithn~ willw ) ; withaspacing of 64 samplesinstead of 128 samples.A
be referred t o as LSEE-MSTFTM.
signal withSTFTMclose
to this MSTFT" was estimated byX(.)
It is possible t o develop ad hoc methods to estimate
starting with an initial white Gaussian noise sequence and thenfromtheMSTFTM
by modifyingtheiterativealgorithminiterating with LSEE-MSTFTM until the distance measure DAM
Fig. 1. For example, suppose we use in one step of the iterawas decreased to the desired level. The Fourier transforms intive procedurethestandardoverlap-addmethodratherthan
thcalgorithm were implemcntedwith 5 I ?-point FFT's. Fig.
i obtainingthenextestimate
theLSEE-MSTFTmethod nw ) / for S 2 = 64 after 100 iterations.
3(b)shows IX{v(n~S2,
x i + ' ( / ? ) fromtheMSTFT
XL,(n7S. w). Thisresults i n the3(c) shows IX{,,(mS,, w ) / after 100 iterations
Similarly,
Fig.
following update equation:of the OA-MSTFTM algorithm using thesame initial estimate.
Comparisons of Fig. 3(b) and 3(c) with Fig. 3(a) indicatc that
the STFTM of the signal estimate is very close t o the desired
and
MSTFTM andthattheperformanceofLSEE-MSTFTM
0.4-MSTFTM is similar. In Fig. 4. the distance measure Div isshown as ;I function o f thcnumber o f itcrations f o r LSEE-

1

1

I

I

1

I

2

J

1

1

I

I

GRIFFIN AND LIM: SIGNAL ESTIMATION FROM MODIFIED STFT

Fig. 2. STFTMthe
ofat“line
up

2 39

I\

screen door.”

LSEE-MSTFTM

0.2 x IO4
0

Fig. 4. DM

I

20

I

I
60

I

40
BO
Number of Iterations

versus iterationnumberofLSEE-MSTFTMand
MSTFTM.

I
IO0

C

OA-

of reverberati0n.k detectable in the signal estimate due to thenonstationarity of the 2: 1 time-scale compressed speech over
the window length.
In addition to the above three examples, other speech material includingnoisyspeechhasbeenprocessedbythetwo
methods at various compression and expansion ratios. Informallisteningappears
toindicatethattheperformanceof(C)
Fig. 3. (a) 128 :64 time-scale compressedSTFTM of originalspeech.
these methods is superior to that of the system by Portnoff(b) STFTM of LSEE-MSTFTM estimate. (c) STFTM of OA-MSTFTM[ l ] . It should be noted that this approach to time-scale modiestimate.
fication of speech differs’considerably from that of Portnoff.
phaseof Y,(rnS, o)is explicitly
MSTFTM and OA-MSTFTM. Although OA-MSTFTM performs In Portnoffs method, theobtainedbyphase
unwrapping
whichisundesirableduetosomewhat better during the initial iterations, LSEE-MSTFTM
variousconsiderations
includingthecomputational
aspect.
eventually surpasses it. This same performance difference wasIn
theLSEE-MSTFTMorOA-MSTFTM
algorithms,
the phase
noted inallof
the examples where these two methods wereofYw(mS,o)is
implicitly
estimatedintheprocessof estimatcompared.Inaddition,LSEE-MSTFTM
wasobserved t o alY,(rnS,
o
)
l
and
no phaseing
asignal
with
STFTM
close
to
ways decrease DM whereas OA-MSTFTM usually stopped deunwrappingisperformed.
creasing DM after about 100 iterations and in some cases inEven though we used a large number of iterations (100) for
creased DM as more iterations were performed.
theexamples illustrated in this paper,
we have observed thatTo show that these methods perform as well for noninteger
essentiallythesameresults
interms
ofspeechquality canbecompression or expansion factors, the second example shows
25
to
100
iterations.
In
addition,
wehave
obobtained
after
a 35 :64 expansion. Fig. 5(a) shows a 35 : 6 4 time-scale modiserved that speech quality improves rapidly initially and then
IY,(rnSI, o ) l ’ calculatedfromtheoriginal
fiedversionof
more
slowly as the number of iterations increases. This is evispeech signal. As in the first example, the initial estimate wasdenced,to some extent, in Fig. 4, where DM decreases rapidly
a white Gaussian noise sequence. Fig. 5(b) and 5(c) show theinitiallybutmore slowly as the number of iterations increases.
100 iterations using a 256
STFTM of the signal estimate afterWithabetterchoice of the initial estimate ofx(n) than a GauspointmodifiedHammingwindowforLSEE-MSTFTMand
sian
noise
sequence,
it may be possible t o reduce the number
OA-MSTFTM,respectively.Inbothoftheseexamples,theofiterations
requiredto achieve a certain performance.
resultantsignalestimatewascleanhighqualityspeechand
Despitethelargenumberof iterations’ required, real time’the estimates produced by LSEE-MSTFTM and OA-MSTFTM
were indistinguishable in listening tests.
‘Due t o iterations, the total number of computations is considerably
The final exampleconsists of a 1 : 2 time-scaleexpansion
larger than Portnoff‘s method [ 1 1. In a multiprocessor environment,of the 2: 1 time-scale compressed speech generated in the first however, the computational requirement of each processor is comparable or perhaps less than that of Portnoff‘s method.
example.TheSTFTMofthesignalestimatesproducedare’The definition of “real time” for time-scale modification depends onthencomparedwiththeSTFTMoftheoriginalspeechsig- the application. In applications where the input to the algorithm is from
nal.Fig.6(a)and6(b)showtheSTFTMofthesignalestisomestorage device andtheoutput
is convertedtoan
analog signal
mates after100 iterations of LSEE-MSTFTM and OA-MSTFTM,which the user listens to, the algorithm must produceoneoutput sample
in an average time less than T I w’here T1 is the sampling period associrespectively.ComparisonsofFig.6(a)and6(b)withFig.
ated with the digital to analog converter used to gcnerateoutput speech.
2 show that both BEE-MSTFTM and OA-MSTFTM produce a In applications where the input
to the algorithm is digitized directly
signal estimate with STFTM close to the STFTM of the original from the user’s speech and the output is placed on some storage device,
the algorithm must process an input data sample in an average time lessspeech signal. Theprimarydifferencebetweenthesesignal
the analog tothan Tz where T2 is the sampling period associated withestimates and the original speech signal is that a small amount digital converter used to digitize the input speech.

1

240

IEEE
TRANSACTIONSON ACOUSTICS, SPEECH,
AND
SIGNALPROCESSING,
VOL.
ASSP-32,

NO. 2, APRIL 1984

Fig. 5 . (a) 35 :64 time-scale expanded STFTM of original speech. (b) STFTM of LSEE-MSTFTM estimate. (c) STFTM of
OA-MSTFTM estimate.

tion. The inherent delay associated with each iteration is only
L data points. This is due
the length of the analysis window,
to the fact that the computational aspect of each iteration of
the algorithm is essentially the same as the weighted overlapadd method [lo], in which the delay between the input andoutpat data is L points assuming the required computation foreachwindoweddatasegmentcanbeperformedduringthe
time corresponding to the window shift, S data points. As an
example that illustrates the computational requirements anddelay involved, suppose SI = Sz = 64, L = 256, the size of the
DFT used is 512, the number of iterations required and thenumber of processors available is 50, and speech is sampled atin thea 10 kHz rate. Since the major computations involved
algorithm are due to the DFT and IDFT,if each processor can
DFT’sonceevery6.4ms,thenthe
Fig. 6 . I: 2 expansion of 2 : 1 compressed speech for (a) LSEE-MSTFTM compute two 512-point
iterativealgorithmcanbeimplementedinrealtimewitha
and (b) OA-MSTFTM.
1.3 s. Currenthardwaretechnology
is moredelayofabout
implementation appears possible if enough processors are used thancapableofhandlingsuchcomputationalrequirements,
and a delay of a few seconds is not a serious problem in mostin series. Specifically, as input data are received, the ith proi + 1st processor applications of time-scale modification of speech.
cessor can perform the ith iteration and the
Even though LSEE-MSTFTM and OA-MSTFTM had similari + 1st iterawhich follows the ith processor can perform the

GRIFFIN AND LIM: SIGNAL ESTIMATION FROM MODIFIED STFT

performanceinthecontextoftime-scalemodificationof
speech,itshould
be pointedoutthatLSEE-MSTFTMdecreases the distance measure DM of (14) in each iteration until
it converges to a critical point, while OA-MSTFTM sometimes
increases D M . In all cases we considered so far, LSEE-MSTFTMDM thanOA-MSTFTMaftera
alwaysproducedasmaller
sufficiently large number of iterations. This difference may besignificant in other applications.
In this paper, we considered the application of the theoreti-

_X_ =

24 1

a) ifx'(n) 4r, then D [ x i + ' ( n ) ] <D[x'(n)]b) ifx'(n)Er, thenD[x"'(n)]<D[x'(n)]and iii) the mapping A is closed at points outside r. Then thelimit of any convergent subsequence of {xi(n)}is in the solution set.
The first requirement of the global convergence theorem
isthat all estimates are contained in the compact set X. Define

2 w2(mS n )
m

-

m =-m
~~~

where O,(mS, w ) E [-T,T] for all m , w .
cal resultsinthispaperonlytotheproblemoftime-scaleof theseresultsto
modificationofspeech.Theapplication
other problems such as enhancement of speech degrade'd byhelium is currently under study and these results will be
re- We willshow that X is compact since it is both closed andbounded. In ordertoensurethatx(n)r-isafiniteiengthseported in a later paper.
quence, the given MSTFTM Y w ( m S ,h)]will be assumed tov. S U M M A R Y
be zero outside of agiven range of m,j In ( A l ) , X has been
Three new algorithms have been presented in this paper. The expressed as a continuous function of the &sed set consistingof the phase angles Q,(rnS, w ) which indicates that X is closed.
firstalgorithm,LSEE-MSTFT,estimatesasignalwithSTFT
closest to a MSTFT and is similar to the overlap-add method. We can further show thatX is bounded as follows:

1

Thesecondalgorithm,LSEE-MSTFTM,
is aniterativealgo-Equation(A2)leads
torithmbased on LSEE-MSTFTwhichwasshownt o converge
criticalsolution
consisting
pointsthesetaof tomagniof a1=
w(mS - n) 2n-.
= li Yw(mS7
tude-only distance 'measure. A third algorithm, OA-MSTFTM,
m = -mIx(n)lis heuristicallydeveloped
based
on
the
overlap-add
method.
. .
LSEE-MSTFTMOA-MSTFTM
time-scale
andapplied
towerew2(mS-- n)modification
speechofwithresults
appearthat
tosuperiorbem =- m
to
[I] .

I

E

1

w)l

2

API,L:NI)tX
In thisAppendix,
we show[hattheiterativealgorithm(LSEE-MSTFTM) in Fig. 1 decreases in each,iterationthe
toacriticaldistancemeasure of (14) andalwaysconverges
point where the gradient of the distance measure of (14) withrespect to x ( n ) is zero. I t should be noted that convergence
to a critical point does not necessarily imply convergence tot h e global m i n i m u m .
To show the above, we use the following global convergence
theorem [ 131 .
Let A be an algorithm on R N , and suppose that, given x o ( n )the sequence { ~ ' ( n ) } r is
= generated~
satisfying'

x i + ( n ) = A [x'(n)].

Lct a solution set P C RN be given, and suppose i) all signal
X C R N , ii)
estimates x'(n) arecontained in acompactset
there is a continuous distance measure D on R N such that

(A31

where

x ( n ) E X.

Jz 1

Therefore, since (1/2n)= -,Y,(mS, a)/dw is bounded
and the sum over m reduces to a finite sum for any single valueof n, then x ( n ) is bounded and so is the set X .
The second requirement is the existence of a distance measure D for a solution set r and the algorithm A that, satisfiesii)oftheglobal_convergencetheorem.
Using the distan_ce
measure of (4), Xh(mS, w ) o,f (1 3) minimizes D [x'(n), X ,( m S , w ) ] for x'(n) fixed and X,(mS, w ) constrained to havemagnitude lY,(rnS, o ) l . Thus, we must have

D [x'(n), '?k(mS, a)]< D [x'(n), 2;- ' (mS, w ) ]

(A41

and

D [ x ' + ' ( n ) ,i $ ' ( m S , w ) ] G D [ x ' + ' ( n ) ,&(rnS, w ) ] .
(A51

IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-32, NO. 2 , APRIL 1984

242

Sincethegradientisthedifferencebetweensuccessive
estimates multiplied by a constant, the solution
set corresponds
to the zeros of the gradient of DM. So, if x’(n) is not an element of r, then x”’(n) f xi(n) and DM[x’+’(n), IYw(mS,w ) < DM [x’(n>, Y,(rn~, w ) ] . If xi(n) is an element ofr, then x”’(n) = x’(n) and DM [xi+’(n), Y,(mS, o)l]<
DM [x’(n>, Yw(mS, 11 .
The final requirement for convergenceis that the mapping
A be closed. Since A is a continuous function of x(n), it mustFig. 7. Successive iterations of LSEE-MSTFTM.
be a closed mapping which
satisfies iii) of the global convergence theorem. Thus, LSEE-MSTFTM converges to a solution
Now, the modified STFT i k ( m S , o)is available which allows set consisting of thecriticalpointsoftheSTFTmagnitude
estimation of thenext signal x’+l(n) usingLSEE-MSTFT.
distance measure DM.
D[x(n), z h ( m S , w ) ] for
Sincethisprocedureminimizesz k ( m S , w ) fured, we have
REFERENCES

I

I]

1

D[x““(n), l?f(mS, o)]< D[x‘(n),i k ( m S , o)]and equality holds if and only(A5) and (A6), we obtain

D [~“‘(n),2:(A6)

if xi+’( n ) = xiin). Combining’(mS, o)]< D [x‘ln),i k ( m S , o)] (A7)and equality holds if and only if xi+’( n ) = xi(.).
Fig. 7 shows
D [x’(n>,l ? t < m ~ 0, 1 1 as segment a, D [x’ (a), i k ( m ~011,
assegmentb,andD[x’+l(n),?k’(mS, a)]assegment c.
Since Xk(mS, 0)and ?k(mS, o)have the same phase, thedistance represented by segment a is equivalent to the distance
between Xk(mS, a)/and Yw(mS,o) This can be shown,u ) ] explicitly:
by writing ~ [ x ’ ( n )2L(ms,
+

I

I

I.

which reduces to-

1 Yw(mS,
w)l]’

(A9)

dw.

This leads to the definition of the distance function based on
STFTmagnitudeA given by (14). Since DM [xi(n), Y,(mS,o ) l ] = D[x’(n), XL(mS,a ) ](A7), can be written as

1

Dn/r[x‘++’(n),I Y d m S , 4 1 1 < DM[Xi(4, iYw(mS,o)ll(A1 0)and equality holds if and only if xi+’ ( n ) = x’(n). Taking thegradient of DM with respect to x(n) yields

I

I

M. R. Portnoff, “Time-scale modificationofspeechbasedon
short-time Fourier analysis,” ZEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 374-390, June S981.
J. S. Lim,“Enhancementandbandwidthcomprcssionofnoisy
speech,”Proc. IEEE, vol. 67, pp. 1586-1604, Dec. 1979.
J. S. Lim, Ed., Speech Enhancement. Englewood Cliffs, NJ:
Prentice-Hall, 1983.
J. B. Allen, D. A. Berkley,and J. Blauert,“Multimicrophonesignal-processing techniques to remove room reverberation fromspeech signals,” J. Acoust. Soc. Amer., vol. 62,pp.912-915,
Oct. 1977.
M. A. Richards,“Heliumspeechenhancement
using theshorttimeFouriertransform,”
ZEEE Trans. Acoust.,Speech, Signal
Processing, vol. ASSP-30, pp. 841-853, Dec. 1982.
L. R. Rabiner and R. W. Schafer, Digital Processing of Speech
Signals. Englewood Cliffs, NJ:Prentice-Hall,1978.
J. B. Allen, “Short term spectral analysis, synthesis, and modifiIEEEDans. Acoust.,
cationbydiscreteFouriertransform,”
Speech, SignalProcessing, vol. ASSP-25, pp. 235-238, Junk 1977.
S. H . Nawab, T. F. Quatieri, and J. S. Lim, “Signal reconstruction
ZEEE Trans.
fromshort-timeFouriertransformmagnitude,”
Acoust., Speech, Signal Processing, vol. ASSP-31, pp. 986-998,
Aux. 1983.
[ 91 M.R. Portnoff, “Time-irequency representation of digital signals
ZEEE Trans.
and systems based on short-time Fourier analysis,”
Acoust.,Speech, Signal Processing, vol. ASSP-28, pp.55-69,
Feb. 1980.
[ 101 R. E. Cyochicre, “A weighted overlap-add method of short-time
Fourieranalysis/synthesis,” IEEE Trans. Acoust.,Speech, Signal Processing, vol. ASSP-28, pp. 99-102, Fcb. 1980.
[ 11 ] R. W. Gerchberg and W. 0. Saxton, “A practicalalgorithmfor
thedeterminationofphasefromimageanddiffractionplane
pictures,” Optik, vol. 35, pp. 237-246, 1972.
[12] A. V. Oppenheim, M. H. Hayes, and J. S. Lim, “Iterative procedures for signal reconstruction from phase,”in Proc. 1980 Znt.
Opt. Comp. Conb, Apr. 1980, vol. SPIE-23 I , pp. 12 1-1 29.
[ 131 D. G. Luenbergcr, Introduction to Linear and Nonlinear Programming. Reading, MA: Addison-Wesley, 1973.
Daniel W. Griffin was born in Detroit,MI, onDccember 18,1960. He received thc B.S. degree in copputcr engineering from the University of Michigan, Ann Arbor, in April 1981.
Hc is currently a Rcsearch Assistant in the
Department of Iclectrical Engineering and Computer Science at the Massachusetts institutc ofTechnology,Cambridge. His rcsearch intcrests
include digital signal processing and its application to specch processing.

IEEE TRANSACTIONS ON ACOUSTICS,
SPEECH,AND SIGNAL PROCESSING, VOL. ASSP-32, NO. 2, APRIL 1984

Jae S. Lim (S’76-M’78-SM’83) wasbornon
December 2, 1950. He received the S.B., S.M.,
E.E., and Sc.D. degrees in electrical engineering
andcomputer science fromthe Massachusetts
Institute of Technology,Cambridge,in1974,
1975, 1978, and 1978,respectively.
He joinedthe
MIT faculty in 1978as an
AssistantProfessor, and is currently Associate
Professor in the Department of Electrical Engineering andComputer Science. His research

243interests include digital signal processing and its applications to image
and speech processing. He has contributed more than sixty articles tojournals and conference proceedings, and is the editor of a reprint book,Speech Enhancement (Prentice-Hall, 1982).
Dr. Lim is the winner of two prize paper awards, one from the Boston
chapter of the Acoustical Society of America in December 1976, andone from the IEEE ASSP Society in April 1979. He is a member of Eta
Kappa Nu, Sigma Xi, and is the chairman of the IEEE ASSP Society’s
Technical Committee onDigital Signal Processing.

Maximum Likelihood Spectral Estimation and Its
Application to Narrow-Band Speech Coding
ROBERT J. McAULAY, MEMBER,IEEE

Abstract-Itakura and Saito[ 11 used the maximum likelihood (ML)method to derive a spectral matching criterion for autoregressive (i.e.,
all-pole) random processes. In this paper, their results are
generalized
to periodic processes having arbitrary model spectra. For the all-polemodel, Kay’s [2] covariance domain solutionto the recursive ML (RML)
problem is cast into the spectral domain and used to obtain the RML
solution for periodic processes. When applied to speech, this leads to amethod for solving the joint pitch and spectrum envelope estimation
problems. It is shown that if the number of frequency power measurements greatly exceeds the model order, then the
RML algorithm reduces
to a pitch-directed, frequency domain version of linear predictive (LP)
spectral analysis. Experiments o n a real-time vocoder reveals that the
RML synthetic speech has the quality ofbeing heavily smoothed.

constraint on the model spectrum. The analysis in Section
I1
leads toaspectralmatchingcriterionidenticaltothatobtained by Itakura and Saito, which shows that the criterion isa fundamental property of the maximum likelihood method.
Furthermore, it is shown that in the periodic case, the model
spectrum is fitted to the power measurements at the harmonicfrequencies.
In Section I11 extensive use is made of results obtained bycase inwhich
Kay [ 2 ] to specialize the MLcriteriontothe
IV, lattice
thespectralenvelope is all-pole.TheninSection
methodsareusedtoderivearecursivemaximumlikelihood(RML)algorithmforestimatingtheARparameters.InSection V the application of the RML technique to speech coding
I . INTRODUCTION
is described, and the results of a perceptual evaluation using a
TAKURA and Saito [ l ] have shown that spectral envelopedis(LPC) real-time analysis/synthesis system are discussed. A brief
estimation using linear predictive coding techniques
VI on the application of thecussion is presented in Section
hasafundamentaltheoreticalbasisinmaximumlikelihoodof thepitchandvocal
ML criterion to the joint estimation(ML) estimation. Furthermore, they have used this theory to
tract spectral parameters. Finally, in Section VII, some generaldevelopaspectralmatchinginterpretationintermsoftheML
conclusions regarding the usefulness of the application of
Itakura-Saito criterion. Their basic mathematical model dealt
with speech waveforms that were sample functionsof an auto- techniques to speech analysis are discussed.
regressive (AR) random process.
While this is an appropriate
11. THEORETICAL FORMULATION
model for the class of unvoiced sounds, one wonders if perhaps
By definition a real, stationary random processs(n) is pethe criterion is valid for voicedspeechsoundsaswell,sinceriodic with period N if its autocorrelation function R ( m ) =in this case the waveforms are periodic. Thisis the problem
E[s(n)s(n+ m ) ] is periodic with period N [4] . Then R(m)addressed in this paper.
In setting up the formalism for applying the ML method for can be expanded by Fourier series asperiodic processes, it was not necessary to impose the all-pole
N- 1
R(m)=
Pk exp ( i m w k )I

k=O

Manuscript received March 23, 1982; revised October 25, 1982, April
29, 1983, and August 26, 1983. This work was sponsored by the Dcpartmcnt of the Air Force. The U S . Government assumes no responsibility for the information presented.
The author is with the Lincoln Laboratory, Massachusetts Institute ofTechnology, Lexington, MA 02173.

where w k = 2nk/N, and where0096-3518/84/0400-0243S01.OO 0 1984 IEEE