INTERSPEECH 2019
September 15–19, 2019, Graz, Austria

A Modified Algorithm for Multiple Input Spectrogram Inversion
Dongxiao Wang1, Hirokazu Kameoka2, Koichi Shinoda1
1

2

Tokyo Institute of Technology, JapanNTT Communication Science Laboratories, Japan

dongxiao@ks.cs.titech.ac.jp, kameoka.hirokazu@lab.ntt.co.jp, shinoda@c.titech.ac.jp

Abstract

duced in [4] for informed source separation (ISS), where magnitudes and phases are jointly updated only when the T-F bin is“active”. Another joint estimation approach named consistentwiener ﬁltering [5] combines the consistency constraint with aloss between the magnitude estimation and classical wiener ﬁlter output, resulting in a promising performance when the given
magnitude spectra are estimated from spectral subtraction. In
recent research [6], MISI procedure is unfolded and integrated
into the chimera++ network, which improves the magnitudeonly methods. We still have a strong motivation of ﬁnding arobust algorithm for the phase estimation, since it can be used
in a source separation setting, as a post-processing procedure
after applying magnitude oriented approaches, such as [7–13].
In this paper, we propose an algorithm which is robust
against the inaccurate magnitude spectrum, by relaxing thestrict constraint as an additional term in the objective function.

We propose a new algorithm to estimate the phase of speech
signal in the mixture of audio sources under the assumptionthat the magnitude spectrum of each source is given. The previous method, multiple input spectrogram inversion algorithm(MISI), often performs poorly when the magnitude spectrograms estimated are not accurate. This may be because it imposes a strict constraint that the summation of source waveforms should be exactly the same as the mixture waveform.
Our proposing algorithm employs a new objective function inwhich this constraint is relaxed. In this objective function, thedifference between the summation of source waveforms and the
mixture waveform is the target to be minimized. The performance of our method, modiﬁed MISI is evaluated on two different experimental settings. In both settings it improves theaudio source separation performance compared to MISI.
Index Terms: Phase estimation, MISI, source separation, synthesis

2. Grifﬁn & Lim Algorithm
Let x = [x(0), . . . , x(T − 1)] ∈ RT be a time domain signal, f = 1, . . . , F and n = 1, . . . , N be frequency and frame
indices, respectively, where F is the number of frequency bins
in each frame, T is the length of the signal, and N is the number of frames. Also let wf,n = [wf,n (0), . . . , wf,n (T − 1)]be a complex sinusoid of frequency ωf modulated by the window function of the n-th frame. Note that wf,n is padded withzeros over the range outside the frame window. Then, the complex Fourier coefﬁcient cf,n of the f -th bin at the n-th frame isH
the inner product between x and wf,n , namely cf,n = wf,nx,
HFN
where · denotes Hermitian transpose. Let c ∈ C
be a vector obtained by stacking cf,n over all the frequency bins F of
all the N frames. Then the relationship between c and x can be
written as
c = Wx.
(1)

1. Introduction
Single channel source enhancement (SCSE) has been widely
applied in hearing-aids, teleconferencing and speech recognition systems. In these applications, a mixture of target signal
and background noise is collected by a single microphone, andnoise reduction methods are applied to remove the backgroundnoise.
SCSE has been intensively studied (for a detailed survey,see [1]). Most approaches use short time Fourier transform(STFT) to transfer time domain mixture signals to spectra, thenestimate the magnitude of target source by several different
ways. Next, they recover the waveform of target source by the
inverse short time Fourier transform (iSTFT) where the phaseof mixture is used for the phase of each individual sources. Of
course the phase of each source is different from that of the
mixture, and so, this process cannot be justiﬁed.
The STFT of a signal is its redundant representation when
we use an analysis window which largely overlaps with nearby
windows; The spectral coefﬁcients of successive frames are correlated. Grifﬁn and Lim [2] utilizes this fact to recover phasefrom magnitude spectrum. For multi-source situation, however,separated sources are not independent with each other, i.e. allof the individual sources should exactly add up to the mixture
signal. Gunawan and Sen [3] proposed multiple input spectrogram inversion (MISI) to utilize this information. While thisapproach works well when the true magnitude of each source isgiven, the performance signiﬁcantly degrades when inaccuratemagnitudes are given.
To address this problem, different approaches have been
proposed for different settings. A method named ISSIR is intro-

Copyright © 2019 ISCA

H
,c
Here W is a F N × T matrix in which each row is wf,n
is a complex spectrogram and W is the the matrix of complex
exponentials functions used in the short time Fourier transform(STFT) modulated by the window function.
We typically use overlapping windows, F N is larger thanT , and accordingly, c is a redundant representation of the signal
x. The inverse STFT of c can be written using the pseudoinverse matrix W+ of W:

W+ c = argmin c − Wx22x= (WH W)−1 WH c.

(2)c is consistent iff c = WW+ c. Grifﬁn and Lim [2] utilize thisconsistency constraint to estimate the phase of each frequency
component.

4569

http://dx.doi.org/10.21437/Interspeech.2019-3242

4. Modiﬁed MISI

Let a and φ be the magnitude and phase spectra of c, respectively. Namely c = a  φ, where φf,n ≡ eiθf,n and 
denotes an element-wise product. Here we estimate φ assuming a is known. Then we can obtain the phase φ by minimizing
the following objective function J (φ):
2J (φ) = a  φ − WW+ (a  φ)2 .

To mitigate the two problems above, we propose a new algorithm, Modiﬁed MISI (M-MISI) which minimizes the following
objective function:aj  φj − WW+ (aj  φj )2
I(φ) =
2

j

(3)


2





+λ y −aj  φ j  .
Instead of solving Eq.(3) directly, we introduce another variablec̃ and minimize the following function:J+(φ, c̃) ≡ a  φ − Wc̃22 .

j

We iteratively update φ and c̃ using the principle of majorization and minimization [14]:(5)

φ ← argmin a  φ − Wc̃22 = ∠Wc̃,

(6)

c̃∈RT
φ

where ∠· denotes an operation that divides each element of avector by its absolute value.

3. Multiple Input Spectrogram Inversionj

For reconstructing individual source signals that compose a
mixture signal, Gunawan and Sen [3] proposed a closed loop
method called multiple input spectrogram inversion (MISI)
where the magnitude spectra of each source is assumed to begiven.
Let ỹ be the mixture of each source, aj and φj be the magnitude and phase spectra of source j, respectively. We introduce
another variable c̃j for source j, then MISI iteratively updatesc̃j and φj :c̃j ← W+ (aj  φj ) +To solve this optimization problem, we also employ
majorization-minimization as in Section (2). We ﬁrst ﬁnd an
auxiliary function of the objective and solve it iteratively.
First let yf,n beyf,n =

j

j

(7)
(8)

=

Although it is not explicitly stated in [3], this algorithm
solves the following constrained optimization problem:=

minimizesubject toJ (φ, c̃) =jc̃j = ỹ.

(13)



2

yf,n −(9)≤  βj,f,n
j


f,n

(10)

=

j



f,n,j

The derivation can be checked in
In MISI, the summation of individual sources should exactly be the mixture signal. This constraint is often too hard
to be satisﬁed and may lead to incorrect updates when only erroneous magnitude spectra are available. Another observationfrom (7) is that the error
 between mixture signal ỹ and the sumof estimated signal j  W+ (aj   φj  ) is equally distributed
over all sources. This also may not be practical since this error
can be different from source to source. While we assume thatthe true magnitude spectrum for each source is given, the magnitude spectrum of one source may be more reliable than that ofother sources in real situations.



2aj,f,n φj,f,n

j

f,n

f,n

aj  φj − Wc̃j 22 ,

xj,f,n ,then from Jensen’s inequality, we can show that
2






aj  φ j 
y −



 +
1⎝
W (aj   φj  )⎠ ,
ỹ −
J


φj ← ∠Wc̃j ,


j

⎞

⎛

2

Here y is a vector of the complex spectrogram of the mixture
signal. The ﬁrst term in I(φ) measures how exactly aj  φjsatisﬁes the constraint from the redundancy, which correspondsto(9) and the second term represents the error between y andj aj  φj , which plays a similar role as (10), but transferred
into time-frequency domain. The motivation is that in the TF domain, we can handle each T-F bin separately. While (11)has the same consistency constraint with [5] where it is denoted
as F, the second term of (11) is different with [5] where theycalculate the quadratic loss between the current estimation andthe classical wiener ﬁlter output.
Instead of equally distributing the error into each source,
here we introduce weight β for each T-F bin:
βj,f,n = 1, 0 < βj,f,n < 1.
(12)

(4)

c̃ ← argmin a  φ − Wc̃22 = W+ (a  φ),

(11)

j

βj,f,nβj,f,n2(14)(xj,f,n − aj,f,n φj,f,n )1βj,f,n(xj,f,n − aj,f,n φj,f,n )2

1|xj,f,n − aj,f,n φj,f,n |2 ,
βj,f,nHere, the equality in (14) holds when⎛
xj,f,n = aj,f,n φj,f,n + βj,f,n ⎝yf,n −⎞aj  ,f,n φj  ,f,n ⎠ .

j

(15)
In the same way as in Section 2, we have
aj  φj − WW+ (aj  φj )22
2≤ aj  φj − Wc̃j 2 =|aj,f,n φj,f,n − cj,f,n |2 ,
f,n

(16)

4570

+

c̃j = W (aj  φj ).

relative SDR

Hwhere cj,f,n = wf,nc̃j . Here the equality in (16) holds when

(17)

From (14) and (16), we can use
J (Φ, x, c̃) =|aj,f,n φj,f,n − cj,f,n |2

j,f,n

j,f,n

(18)

λ
|xj,f,n − aj,f,n φj,f,n |2
βj,f,nλ|x − aφ|2 .
β−10 dBMISI G&L 10−2 10−1

+ d,|xj,f,t − aj,f,t φj,f,t |j  xj ,f,t − aj ,f,t φj ,f,t

−10 dBMISI G&L 10−2 10−1⎛

xj ← aj  φj + βj  ⎝y −

aj   φ j  ⎠

0 dB

103

104

10 dB

105

relative SIR

0.5
0.0
−0.5

−10 dB

1

0 dB

103

104

10 dB

105

where t and f are the time and frequency indices respectively,
the local SNR for each T-F bin is deﬁned by
SNR(t, f ) = 10 log10(22)|S(t, f )|2,|N (t, f )|2(28)

where S and N are the spectrogram of the clean speech andnoise, respectively. θ is a local criterion(LC) in dB, which is
chosen to be 0 dB. To add errors to the IBM values, we randomly ﬂip the value with different rate between 0 and 20%,with 5% step. In 2), we only assume knowing the time averageof power spectrum of the noise signal, and use power spectrum
subtraction to recover the value of speech signal.
For both conditions, the mixtures are generated by mixing
a speech signal taken from the ATR speech database [15], and astationary noise from the noise database published by JapaneseStandards Association (JSA) [16]. All of the conversational
sentences in ATR (115 audio clips in total) are used for generating mixtures. The noise dataset contains 20 types of in-door
noises including air conditioner, washing machine, and dryer.
All utterances are re-sampled to 16kHz. The overlap ratio of
each frame is set to 50%. The window size is set to 20ms.
For each condition, the Signal-to-Distortion, Interference
and Artifacts ratio (SDR, SIR and SAR) are used as evaluation
metrics. The performance of the proposed method is comparedwith MISI and Grifﬁn & Lim algorithm.



⎞

1

We evaluate the signal reconstruction performance of M-MISIin two different conditions, 1) Binary Mask and 2) Spectrum
Subtraction. In both conditions, we assume that only imperfectmagnitude spectra is known beforehand. In 1), The magnitude
spectrum of each source is estimated by applying a binary mask
to the mixture magnitude spectra. Here the ideal binary mask iscomputed with

1 if SNR(t, f ) ≥ θIBM(t, f ) =(27)
0 otherwise(20)

Thus we can iteratively update the parameters:c̃j ← W+ (aj  φj )

102

−0.5

5.1. Setup

(21)

.

101

10 dB

105

5. Experimental Evaluation

where d is other terms that does not related to the optimal solution of φ. (21) achieves its minimum value when both terms inside have the same phase which is φ = ∠ βc+λx. For β, by minβ+λimizing (14) with constraint (12) when ﬁxing xj,f,t and φj,f,n ,we have
βj,f,t = 

102

104

Figure 1: Output SDR, SAR and SIR of different λ values inSpectral Subtraction case, averaged over all mixtures for each
input SNR, subtracting the mean value for each input SNR.
Numbers on the x-axis stand for different λ values in M-MISI.

Here ¯· denotes the complex conjugate. Then by completing thesquare, (20) becomes
2

101

0 dB

103

(c)

.

(β + λ) βc + λx
− aφ
β
β+λ

102

0.5

By expanding it to a second order polynomial of aφ, (19) becomes

β |c|2 + λ |x|2+
(β + λ)101

0.0

MISI G&L 10−2 10−1(19)(β + λ) 2 2 (βc̄ + λx̄)aφ + (βc + λx)aφ̄
a |φ| −
β(β + λ)1(b)

as an auxiliary function for optimizing I(φ). where d is otherterms that are not related to the solution.
We have discussed how to update c̃j and xj respectively.
Next we explore how to update the parameters φj , and βj .
First, for φj , since this term appears in both of the termsin (18), we merge those two terms together. Since the optimal value of φj,f,n can be found independently, we omit thesubscript in order to keep the notation uncluttered. Then theobjective function for each (j, f, n) is given by|aφ − c|2 +

−0.5

(a)
relative SAR

+



0.5
0.0

(23)
(24)

j



φj ← ∠ (βj  Wc̃j + λxj )  (βj + λ1)
(25)
⎛
⎞

abs(xj  − aj   φj  )⎠ ,βj ← ⎝abs(xj − aj  φj ) j(26)
where 1 is an all-ones vector,  denotes the element-wise division and abs(·) denotes an operation that takes the absolutevalue of each element of a vector. As for the value of λ, it is set
to be a hyperparameter to be optimized.

4571

Table 1: Comparison of the SDR, SAR and SIR for the G&L[2] algorithm, MISI[3] and proposed Modiﬁed-MISI, with λ set to 103 foreach input SNR, in 1) Binary Mask and 2) Spectral Subtraction conditions−10 dB
SAR
10.59
10.11
10.27
5.28
3.79
4.55

Input SNR
BinaryMask
Spectral
Subtractioin

G&L[2]MISI[3]
MMISI (Proposed)
G&L[2]MISI[3]
MMISI (Proposed)-10 dB

SDR
5.74
5.76
6.01
−6.46
−6.34
−5.52

0 dB

-4

10 dB

5

14

SDR

4
13
3

-8

12

2
-10

MISI

G\&L MMISI

1

MISI

G\&L

MMISI

11

MISI

SDR
13.62
13.40
13.68
3.83
3.93
4.79

0 dB
SAR
16.81
16.07
16.46
9.98
9.30
10.25

SIR
17.67
17.87
18.02
5.55
6.40
6.91

SDR
20.19
19.03
19.71
13.25
13.44
14.26

+10 dB
SAR
22.11
20.36
21.27
17.96
17.93
18.82

SIR
25.00
25.10
25.20
15.26
15.85
16.46

is conducted with different error rate, as mentioned in Section
5.1, and the upper rows show the average values of all experiments. Since G&L algorithm does not involve error distribution, it achieves the best SAR values in most of the times. For
condition 2), we show further the statistics in Figure 2 for each
input SNR. While in both conditions, the separation quality isdominated by the quality of magnitude estimations, our proposed method still achieved 0.4 dB in condition 1) and 1 dB in
condition 2) compared with MISI.

15

6

-6

SIR
9.88
10.18
10.33
−4.91
−3.88
−3.44

G\&L MMISI

Figure 2: The box plot showing the SDR performance of theSpectral Subtraction case

6. Conclusion
7.2

We have proposed a new phase recovery algorithm, M-MISI,
which employs a soft objective function. This algorithm isrobust against imperfect magnitude spectra of the estimatedsources. In our experiment, the proposed method achieved 1 dB
improvement when we only use spectral subtraction to recover
the magnitude spectrograms, and 0.4 dB on average when weuse binary mask as the magnitude estimator. In future, we needto combine our method with recent source separation methods
involving deep learning. Also, we would like to tackle the casewhen the number of recording channels (microphones) is morethan one.

10 -5

λ = 102
λ = 1037

10 -4

λ = 105
λ = 104

10

-3

10 -2

6.8

10 -1

SDR

0.0

G&L

6.6

1
10
10 2

6.4

10 3

MISI

10 4
10 5

6.2

misi

6

G&L

0

50

100

150

200

250

300

# of iterations

7. Appendix

Figure 3: Output SDR of MISI, G&L(λ = 0), and M-MISI withdifferent λ, as a function of number of iterations.

7.1. The derivation of MISI
To conﬁrm that (9) and (10) lead to the solution given in (7) and(8), we ﬁrst deﬁne the Lagrangian asHc̃j − ỹ .
L(φ, c, γ) = J (φ, c̃) + γ

5.2. Experimental Results
Figure 1 compares different λ values in terms of SDR, SAR and
SIR, relative to the average value of each input SNR. Note that
the G&L algorithm corresponds to the case where λ = 0. It isshown for each input SNR, the proposed method outperforms
MISI when λ > 10. While the proposed method achieves similar performances when λ > 102 , we ﬁnd that the convergence
becomes slower when λ becomes greater, which could be further conﬁrmed from Figure 3. In other experiments, we choose
λ = 103 for a best trade off between execution time and performance. Another observation from Figure 3 is that it is difﬁcult
to stop MISI when it achieves the best performance, since theperformance degrades drastically after it reaches the peak at thevery beginning of the execution, while M-MISI is assured to
improve the SDR by running for more iterations, at a cost of
more computational time.
Table 1 shows the comparison between M-MISI, MISI, andGrifﬁn & Lim algorithm. For each experiment, the phase spectra of each source is computed by running 200 iterations of
all of the 3 algorithms. In condition 1), the same experimentj

By setting the partial derivative of the Lagrangian with respect
to c̃∗j at zero, we obtain
WH Wc̃j = WH (aj  φj ) − γ.
By summing this expression over j, we get H
W (aj  φj ) − WH Wỹ.
Jγ =jHence, when {φj }j is ﬁxed, (9) is minimized with respect to
c̃j subject to (10) when


 +
1
+
W (aj  φj ) .
c̃j = W (aj  φj ) +ỹ −JjIt is straightforward to show that when {c̃j } is ﬁxed, (9) is minimized with respect to φj whenφj = ∠Wc̃j .

4572

8. References[1] R. C. Hendriks, T. Gerkmann, and J. Jensen, “Dft-domain based
single-microphone noise reduction for speech enhancement: A
survey of the state of the art,” Synthesis Lectures on Speech andAudio Processing, vol. 9, no. 1, pp. 1–80, 2013.
[2] D. Grifﬁn and J. Lim, “Signal estimation from modiﬁed shorttime fourier transform,” IEEE Transactions on Acoustics, Speech,
and Signal Processing, vol. 32, no. 2, pp. 236–243, April 1984.
[3] D. Gunawan and D. Sen, “Iterative phase estimation for the synthesis of separated sources from single-channel mixtures,” IEEE
Signal Processing Letters, vol. 17, no. 5, pp. 421–424, May 2010.
[4] N. Sturmel and L. Daudet, “Informed Source Separation using
Iterative Reconstruction,” arXiv:1202.2075 [cs], Feb. 2012,arXiv: 1202.2075. [Online]. Available: http://arxiv.org/abs/1202.
2075[5] J. L. Roux and E. Vincent, “Consistent Wiener Filtering for Audio Source Separation,” IEEE Signal Processing Letters, vol. 20,
no. 3, pp. 217–220, Mar. 2013.
[6] Z.-Q. Wang, J. L. Roux, D. Wang, and J. R. Hershey, “End-to-end
speech separation with unfolded iterative phase reconstruction,”
arXiv preprint arXiv:1804.10204, 2018.
[7] A. Hyvärinen and E. Oja, “Independent component analysis: algorithms and applications,” Neural networks, vol. 13, no. 4-5, pp.
411–430, 2000.
[8] M. Davies and C. James, “Source separation using single
channel ica,” Signal Processing, vol. 87, no. 8, pp. 1819 –
1832, 2007, independent Component Analysis and Blind Source
Separation. [Online]. Available: http://www.sciencedirect.com/
science/article/pii/S0165168407000151[9] M. N. Schmidt and R. K. Olsson, “Single-channel speech
separation using sparse non-negative matrix factorization,” inInterspeech, sep 2006. [Online]. Available: http://www2.imm.
dtu.dk/pubdb/p.php?4511[10] T. Virtanen, “Monaural sound source separation by nonnegative
matrix factorization with temporal continuity and sparseness criteria,” IEEE transactions on audio, speech, and language processing, vol. 15, no. 3, pp. 1066–1074, 2007.
[11] J. R. Hershey, Z. Chen, J. L. Roux, and S. Watanabe, “Deepclustering: Discriminative embeddings for segmentation and separation,” in 2016 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP), March 2016, pp. 31–35.
[12] Y. Isik, J. L. Roux, Z. Chen, S. Watanabe, and J. R. Hershey,“Single-channel multi-speaker separation using deep clustering,”
arXiv preprint arXiv:1607.02173, 2016.
[13] M. Kolbæk, D. Yu, Z.-H. Tan, J. Jensen, M. Kolbaek, D. Yu, Z.-H.
Tan, and J. Jensen, “Multitalker speech separation with utterancelevel permutation invariant training of deep recurrent neural networks,” IEEE/ACM Transactions on Audio, Speech and Language
Processing (TASLP), vol. 25, no. 10, pp. 1901–1913, 2017.
[14] D. R. Hunter and K. Lange, “A tutorial on mm algorithms,” The
American Statistician, vol. 58, no. 1, pp. 30–37, 2004.
[15] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara,
and K. Shikano, “Atr japanese speech database as a tool of speech
recognition and synthesis,” Speech communication, vol. 9, no. 4,
pp. 357–363, 1990.
[16] J. S. Association, “Design guidelines for alarm sound of consumer
products-environmental sound database,” JIS/TR S 0001, 2002.
[Online]. Available: https://ci.nii.ac.jp/naid/10011785283/

4573