Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014

FAST SIGNAL RECONSTRUCTION FROM MAGNITUDE SPECTROGRAM OF
CONTINUOUS WAVELET TRANSFORM BASED ON SPECTROGRAM CONSISTENCYTomohiko Nakamura† and Hirokazu Kameoka†‡†

Graduate School of Information Science and Technology, The University of Tokyo,
7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan.
‡
NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation,
3-1, Morinosato Wakamiya, Atsugi, Kanagawa, 243-0198, Japan
nakamura@hil.t.u-tokyo.ac.jp, kameoka@hil.t.u-tokyo.ac.jp
ABSTRACTThe continuous wavelet transform (CWT) can be seen as a filterbank having logarithmic frequency subbands spacing similar to thehuman auditory system. Thus, to make computers imitate the significant functions of the human auditory system, one promising
approach would be to model, analyze and process magnitude spectrograms given by the CWT. To realize this approach, we must beable to convert a processed or modified magnitude CWT spectrogram, which contains no information about the phase, into a time
domain signal specifically for those applications in which the aim
is to generate audio signals. To this end, this paper proposes a fastalgorithm for estimating the phase from a given magnitude CWTspectrogram to reconstruct an audio signal. The experimental results revealed that the proposed algorithm was around 100 times
faster than a conventional algorithm, while the reconstructed signals obtained with the proposed algorithm had almost the sameaudio quality as those obtained with the previous study.
1. INTRODUCTION
The continuous wavelet transform (CWT), also known as theconstant-Q transform, is used as a method for time-frequency analysis, which provides a time-frequency representation of a signal
with an equal resolution on a log-frequency scale (Fig. 1). The
human auditory filterbank is known to have an equal resolution ona log-frequency scale as with the CWT particularly in a high frequency band [1, 2]. Thus, to let computers imitate the significant
functions of the human auditory system, one promising approach
would be to model, analyze and process spectrograms obtained by
the CWT (CWT spectrogram). In fact, recent studies (see [3–6])
have shown that multiple fundamental frequency estimation performs very well in the magnitude CWT spectrogram domain. Motivated by this fact, we believe that source separation and soundmanipulation can also work well in the magnitude CWT spectrogram domain. However, in order to achieve source separation orsound manipulation, in which the goal is to produce sound, there
is a need to reconstruct an appropriate time-domain signal after
processing and modifying a magnitude CWT spectrogram. To thisend, this paper proposes a method for estimating the phase from agiven magnitude CWT spectrogram to reconstruct an audio signal.
The phase estimation algorithm from a magnitude CWT spectrogram has already been proposed by Irino et al. [7]. Irino’s algorithm consists in iteratively performing the inverse CWT and the
CWT followed by replacing the modified magnitude CWT spectrogram with a given magnitude CWT spectrogram. Since the(a) CWT spectrogram.

(b) STFT spectrogram.

Figure 1: Examples of the continuous wavelet transform (CWT)
and short-time Fourier transform (STFT) spectrograms. While
STFT spectrograms have an equal resolution on a linear frequency
scale, CWT spectrograms have an equal resolution on a logfrequency scale.

computational speed of the CWT is much slower than the shorttime Fourier transform (STFT), this algorithm needs a very longtime for computation. In practical situations, the reduction of thecomputational complexity can be extremely important.
The authors and colleagues have thus far proposed a fast method
for estimating the phase from a magnitude STFT spectrogram [8].
When the hop-size is shorter than the frame length, the waveforms
in the overlapping segment of consecutive frames must be consistent. This implies the fact that an STFT spectrogram is a redundant representation. Thus, an STFT spectrogram must satisfy
a certain condition to ensure that it is associated with a time domain signal. We have referred to this condition as the consistency
condition. In [8], we have shown that the problem of estimating
the phase from a magnitude STFT spectrogram can be formulated
as the problem of optimizing the consistency criterion describing
how far an arbitrary complex array deviates from this condition.
It became clear that the deviced algorithm is equivalent to thewell-known algorithm proposed by Griffin et al., [9]. The formulation derived from the concept of the spectrogram consistency has
provided a new insight into the Griffin’s algorithm, allowing us to
introduce a fast approximate algorithm and give a very intuitive
proof of the convergence of the algorithm. Since a CWT spectrogram is also a redundant representation of a signal [10], we may beable to make the best use of the spectrogram consistency concept
to develop a fast approximate method for phase estimation from amagnitude CWT spectrogram.
Following the idea proposed in [8], this paper derives an algo-

DAFX-1

Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014
where H is used to denote the Hermitian transpose. This implicitly means that the inverse CWT is defined as the solution to thefollowing minimization problem:

Inverse CWTInconsistentW + s = argmin ∥s − W f̃ ∥22 ,

CWT

(4)

f̃ ∈F

CWT

Consistent
Set of discrete- Inverse CWTtime signals

Set of CWT spectrograms

Figure 2: Illustration of the spectrogram consistency concept forthe continuous wavelet transform (CWT).

rithm for estimating the phase from a magnitude CWT spectrogram. Sec. 2 formulates the phase estimation problem as an optimization problem based on a consistency condition. Sec. 3 derives
an iterative algorithm for phase estimation based on an auxiliary
function approach, which turns out to be equivalent to the algorithm proposed by Irino [7]. Our formulation gives a very clear
proof of the convergence of the algorithm, though it should be
noted that the proof of the convergence has already been givenin [10]. Sec. 4 describes a fast approximate method for computing
each iterative step of the proposed algorithm.
2. CWT SPECTROGRAM CONSISTENCY
2.1. Consistency condition
The scale parameter of the CWT corresponds to the period (thereciprocal of the center frequency) of the wavelet basis function.
Here we consider discretizing the scale parameter such that thecenter frequencies of the wavelet basis functions are uniformly
spaced on a log-frequency scale. Let the indices of the scale parameter and time shift parameter be denoted by l ∈ [0, L−1] andt ∈ [0, T −1], respectively, and let the component of a CWT spectrogram associated with the l-th scale parameter al > 0 (hereafter,
the l-th component) be denoted by sl = [sl,0 , sl,1 , . . . , sl,T −1 ]T ∈ CT .
Given a discrete-time
signal f = [ f0 , f1 , . . . , fT −1 ]T ∈ F where F :=′T ∑ ′{ f ∈ C ; t ft = 0}, its CWT spectrogram s = [sT0 , sT1 , . . . , sTL−1 ]T ∈CLT is defined as
s = W f,(1)

where ∥ · ∥2 denotes the ℓ2 norm of a vector.
While the CWT spectrogram of an audio signal (i.e., a complex
vector that belongs to the subspace spanned by the column vectors of W) will be mapped to itself by applying the inverse CWTfollowed by the CWT, a complex vector that does not belong to
the subspace will not come back to the same point but will be projected onto the nearest point in the subspace. Thus, we can definea condition for a complex vector to be “consistent” (in the sensethat it corresponds to a CWT spectrogram of a signal) as follows:
0LT = s − WW + s,

where 0LT denotes an LT -dimensional zero vector. It is important to note that when W is replaced with a matrix in which eachrow is a basis function of the STFT, (5) becomes equivalent to theconsistency condition for an STFT spectrogram proposed in [8].
2.2. Phase estimation using spectrogram consistency
When given a magnitude CWT spectrogram a ∈ [0, ∞)LT , we canconstruct a signal by assigning phase ϕ ∈ [−π, π)LT to it to obtain a complex spectrogram s, and applying the inverse CWT, i.e.,
W + s. Here, if we assign “inconsistent” phase to the given magnitude spectrogram, the complex spectrogram s will not belong
to the signal subspace and so the spectrogram of the constructedsignal, WW + s, will be different from s. As we want to keep themagnitude spectrogram of the constructed signal consistent withthe given magnitude spectrogram, we must find “consistent” phasesuch that s satisfies the consistency condition.
2.3. Filter bank interpretation
To give a deeper insight into the consistency condition, we focus
on the filter bank interpretation of the CWT. The CWT of a signal
can be thought of as the output of a filter bank consisting of subband filters whose impulse responses are given by the the scaled
mother wavelets. Now, by applying the T -point discrete Fouriertransform (DFT) to each block of (5), (5) can be written equivalently as
0 = ŝ − Ŵ Ŵ + ŝ,

where W ∈ CLT ×T denotes the CWT matrix, defined as W0  W1 W :=  .  , .. 


WL−1


 ψl,0ψl,T −1
Wl :=  .
 ..

ψl,1

ψl,1
ψl,0
ψl,2

···
···
..
.
···


ψl,T −1 ψl,T −2 ..  .
. 
ψl,0

W + := (W H W)−1 W H ,

(6)

where


 Ŵ0 
 Ŵ 
 1 
Ŵ =  .  , .. 


ŴL−1

(2)

Here, ψl,t := ψ(t∆/al )/al is a scaled mother wavelet with thescale of al and the time shift of t∆, where ∆ denotes the sampling period of the time shift parameter and ψ(·) ∈ C denotes themother wavelet satisfying the admissibility condition. Each row ofWl ∈ CT ×T contains the wavelet basis function of scale al with adifferent time shift parameter. The inverse CWT can be defined by
the pseudo-inverse of W:
f = W + s,

(5)

(3)

Ŵl = FT Wl FTH ,Ŵ + = (Ŵ H Ŵ)−1 Ŵ H ,(7)

FT ∈ CT ×T is the DFT matrix and ˆ· denotes the DFT of a variable.
Since Wl is a circulant matrix, Wl is diagonalized by FT and FTH .
The diagonal elements of Ŵl represent the frequency response of
the l-th subband filter associated with the scale parameter al . The
k-th diagonal element of (6) is explicitly written as
1 ∑0 = ŝl,k −ψ̂l,k ψ̂∗l′ ,k ŝl′ ,k ,(8)
C k l′DAFX-2

Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014
where k ∈ [0, T − 1] denotes the angular frequency index, Ck is anormalization constant, and ∗ is used to denote the complex conjugate.
If the subbands of the filter bank overlap each other (moreprecisely, if there exists a pair of channels such that the productof their frequency responses is non-zero at every frequency), i.e.
∀k, ∃l , l′ , ψ̂l,k ψ̂l′ ,k , 0, (5) becomes a nontrivial condition for acomplex vector s ∈ CLT to correspond to a consistent CWT spectrogram. Otherwise, all the elements of CLT trivially satisfy (5),implying that the consistency condition cannot be used as a criterion for phase estimation. Therefore, care must be taken in choosing the quantization intervals of the scale parameter and the typeof the mother wavelet function. The Morlet [11], the log-normal
wavelet [4] and the wavelets used in the auditory wavelet transform [7] satisfy the above requirement when the quantization intervals of the scale parameter are appropriately chosen. We hereafter
assume to use a filter bank that satisfies ∀k, ∃l , l′ , ψ̂l,k ψ̂l′ ,k , 0.
The requirement for the subbands of the CWT to overlap eachother is analogous to the requirement for the short time frames of
the STFT to overlap. The consistency condition of STFT spectrograms can be understood as implying that the waveforms within
the overlapping segment of consecutive frames must be consistent [8]. The consistency condition of CWT spectrograms, on theother hand, can be interpreted as implying that the outputs of adjacent channels within the overlapping subbands must be consistent.

Set of complex vectors withthe given magnitude CWT spectrogram

Inverse CWT

CWT
Update phases

Figure 3: Illustration of the iterative phase estimation algorithm.
The red and blue arrows correspond to (14) and (15).
shall be omitted owing to space limitations. Thus, I+ (ϕ, s̃) should
be designed as a function that can be minimized analytically withrespect to ϕ and s̃. Such a function can be constructed as follows.
Recall that the operator WW + is an orthogonal projection onto
the subspace spanned by the column vectors of W and so WW + s
indicates the closest point in the subspace from s. Thus, we canshow that
I(ϕ) = min ∥s(a, ϕ) − W f̃ ∥22

3. PHASE ESTIMATION BASED ON CWTSPECTROGRAM CONSISTENCYf̃ ∈F= min ∥s(a, ϕ) − s̃∥22 ,s̃∈W

3.1. Formulation of phase estimation problem
Assume that we are given a magnitude CWT spectrogram, arranged as a non-negative vector a ∈ [0, ∞)LT . We would liketo estimate the phase of the given magnitude spectrogram suchthat it meets the consistency condition. To allow for any vectora ∈ [0, ∞)LT as the input, we consider finding a phase estimateϕ ∈ [−π, π)LT that minimizes the consistency criterion
I(ϕ) := ∥s(a, ϕ) − WW + s(a, ϕ)∥22 ,(9)where s(a, ϕ) denotes the estimated CWT spectrogram defined by jϕ0  e  ejϕ1 s(a, ϕ) := a ⊙  .  .
 ..  jϕ 
e LT −1

Set of CWT spectrograms

Set of discretetime signals⊙ denotes the element-wise product. I(ϕ) describes how fars(a, ϕ) deviates from the consistency condition. Namely, the moreconsistent s(a, ϕ) becomes, the smaller I(ϕ) becomes.
3.2. Iterative algorithm with auxiliary function approach
Unfortunately, the optimization problem of minimizing I(ϕ) withrespect to ϕ is difficult to solve analytically. However, we can
invoke the auxiliary function approach to derive an iterative algorithm that searches for the estimate of ϕ, as with [8]. To apply theauxiliary function approach to the current optimization problem,the first step is to construct an auxiliary function I+ (ϕ, s̃) satisfying I(ϕ) = min s̃ I+ (ϕ, s̃). We refer to s̃ as an auxiliary variable.
It can then be shown that I(ϕ) is non-increasing under the updates
ϕ ← argminϕ I+ (ϕ, s̃) and s̃ ← argmin s̃ I+ (ϕ, s̃). The proof of this(12)

where W denotes the set of consistent CWT spectrograms (thesubspace spanned by the column vectors of W). Therefore, we canconfirm that
I+ (ϕ, s̃) := ∥s(a, ϕ) − s̃∥22 , s̃ ∈ W,

(13)

satisfies I(ϕ) = min s̃∈W I+ (ϕ, s̃). (13) can thus be used as
an auxiliary function for I(ϕ). We can thus monotonically decrease I(ϕ) by iteratively performing s̃ ← argmin s̃ I+ (ϕ, s̃) andϕ ← argminϕ I+ (ϕ, s̃). Here, s̃ ← argmin s̃ I+ (ϕ, s̃) and ϕ ←argminϕ I+ (ϕ, s̃) can be written explicitly ass̃ ←WW + s(a, ϕ),ϕ ←∠s̃,

(10)

(11)

(14)
(15)

respectively, where ∠ denotes an operator that gives the argumentsof the components of a complex vector as a real vector in [−π, π)LT .
(14) means applying the inverse CWT followed by the CWT
to s(a, ϕ). Here, when s(a, ϕ) is already a complex vector corresponding to a consistent spectrogram, this update simply becomess̃ ← s(a, ϕ). (15) means replacing the phase estimate ϕ with thephase of s̃. A schematic illustration of these updates is shown in
Fig. 3. I(ϕ) = 0 indicates that s(a, ϕ) lies in the intersection of
the set of consistent CWT spectrograms and the set of complex
vectors that are equal to a up to a phase factor.
3.3. Relation to previous work
The present algorithm is equivalent to an algorithm proposed by
Irino [7]. In addition, when W is replaced with a matrix in which
each row is a basis function of the STFT, the present algorithm
becomes equivalent to the phase estimation algorithm for a magnitude STFT spectrogram proposed in [8].

DAFX-3

Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014
outside the range k ∈ [B, B + D − 1] (0 ≤ B, 0 < D ≤ T ), we canapproximate sl reasonably well by using the elements of {ψ̂l,k fˆk }konly within that range and neglecting the remaining elements. Thisimplies the possibility of computing an approximation of sl with alower computational cost.
For simplicity of notation, let us put Gl,k = ψ̂l,k fˆk . We are concerned with computing an approximation of the full-band inverse
DFT of Gl,k :Figure 4: Example of the frequency responses of different subband
filters (i.e., the scaled mother wavelets). The mother wavelet is the
log-normal wavelet [4].

sl,t =

T −1
∑

2πkt

Gl,k e j T .

(17)

k=0

As mentioned above, Gl,0 , . . . , Gl,T −1 can be approximately viewed
as a band-limited spectrum. In general, the inverse DFT of a bandlimited spectrum can be computed by taking the inverse DFT overthe finite support. In the time domain, this process corresponds
to downsampling the signal given by the “full-band” inverse DFT.
The proposed method uses this idea to approximate the inverse
DFT of the full-band spectrum Gl,0 , . . . , Gl,T −1 . Now, if we choose
D such that T/D becomes an integer, we can approximate thedownsampled version of sl,t byAngularfrequency䠖

(i) Bandlimiting

Frequencyindexs̃l,d =

(ii) Circular shifting

B+D−1
∑

2πkd

Gl,k e j D =

k=B

B+D−1
∑

Gl,k e j

2πk(T/D)d
T

.

(18)

k=B

By comparing (17) and (18), we can confirm thatsl,(T/D)d ≃ s̃l,d (d ∈ [0, D − 1]),

if we assume Gl,k ≃ 0 outside the range k ∈ [B, B + D − 1]. Sinces̃l,d can be rewritten as

Frequency
index

Angularfrequency䠖s̃l,d =

D−1
∑

2πk

B

B

Gl,k+B e j( D +2π D )d = e j2π D d

k=0Figure 5: A circularly shifted version of Gl,B , . . . , Gl,B+D−1 [12].

4. FAST PHASE ESTIMATION ALGORITHM
4.1. Fast approximate continuous wavelet transform
The CWT and the inverse CWT are computationally expensive
compared to the STFT and the inverse STFT. Here we briefly describe the fast approximate method for computing CWT proposedin [12]. The proposed fast approximate CWT uses the fact that thedominant part of the frequency response of each subband filter isconcentrated around its center frequency (as shown in Fig. 4), as
is common in many types of mother wavelets including the Morlet
and log-normal wavelets [4].
According to the filter bank interpretation of the CWT, theCWT of an input signal, sl = [sl,0 , . . . , sl,T −1 ]T = Wl f , canbe computed by multiplying the DFT of the entire signal, i.e,f̂ = [ fˆ0 , . . . , fˆT −1 ]T = FT f , by the frequency response of the lth subband, i.e., Ŵl = diag(ψ̂l,0 , . . . , ψ̂l,T −1 ), and then computing
the inverse DFT of Ŵl f̂ . This can be confirmed fromsl = Wl f = FT FT Wl FT FT f = FT Ŵl f̂ .
H

H

H

(19)

(16)

Note that the second equality follows from the fact that the DFT
matrix FT is a unitary matrix, i.e., FTH FT = IT . Here, if we canassume that the elements of {ψ̂l,k }k are dominant within and near 0

D−1
∑

2πkd

Gl,k+B e j D ,

(20)

k=0

we notice that s̃l,d can be computed by multiplying the inverse DFT
Bof Gl,B , . . . , Gl,B+D−1 by e j2π D d . Note that this is equivalent to computing the inverse DFT of a circularly shifted version of Gl,k (seeFig. 5):(k = 0, . . . , B − (n − 1)D − 1)
Gl,k+nDG̃l,k = (21),Gl,k+(n−1)D (k = B − (n − 1)D, . . . , D − 1)
where n is an integer such that
n−1<

B
≤ n.
D

(22)

We consider invoking the fast Fourier transform (FFT) algorithm
for computing the inverse DFT and so we assume the size D to
be a power of 2. Since D < T , the computational cost for computing s̃l,0 , . . . , s̃l,D−1 is obviously lower than that for computing
sl,0 , . . . , sl,T −1 .
4.2. Fast phase estimation algorithm
The processes of bandlimiting and circular shifting can be represented by a matrix K:
[
][
]
0
IB0
0
ID 0D×(T −D−B)
K := B0 ×(D−B0 )
(23)
ID−B0
0(D−B0 )×B0 | D×B
{z
}
|
{z
}
(i) Bandlimiting

DAFX-4(ii) Circular shiftingProc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014
where ID and 0D×B are the D × D identity matrix and the D × B
zero matrix. The downsampled version of sl obtained with theabovementioned fast approximate CWT can be described as

Computation time[s/iteration]šl = F DH K Ŵl FT f .

1000

(24)

Similarly to the inverse CWT, the fast approximate version of
the inverse CWT can be defined by the pseudo-inverse matrix of
F DH K Ŵl FT . It is important to note that convergence of the phaseestimation algorithm in which the CWT and inverse CWT steps
are replaced with the fast approximate versions is still guaranteed.

Proposed
Conventional

100
10
1
0.1

0.01
0.001

4.3. Time and space complexity0

The computational costs for the CWT and the fast approximate
CWT mainly depend on the number of the points for the inverse
DFT. Since the computational complexity of the full band inverse
DFT is O(T log2 T ), the total computational complexity of the
CWT is O(T log2 T + LT log2 T ). By contrast, the computational
complexity of the band-limited DFT is O(D∑
log2 D) and so the total
L−1computational complexity is O(T log2 T + l=0
Dl log2 Dl ).
The space complexity of the proposed algorithm is small compared to Irino’s algorithm [7]. When the signal length T is longenough, the space complexity depends primarily on the size of
the CWT spectrogram. While the size of the CWT spectrogramof∑ Irino’s algorithm is LT , that of the proposed algorithm is onlyl Dl .
5. EXPERIMENTAL EVALUATIONS
5.1. First experiment: Computation time and audio quality
5.1.1. Experimental conditions
To evaluate the computation time and the audio quality of the reconstructed signals by the phase estimation algorithms, we conducted an objective experiment, and compared the proposed algorithm with the Irino’s algorithm [7].
We used the magnitude CWT spectrograms of 16kHz-sampledacoustic signals of the 113 male and 115 female speeches in the
ATR Japanese speech database A-set [13]. The FFT performs
faster for acoustic signals with a power of 2 length than those withthe other length, and the used signals were filled by 0 till eachlength reached a power of 2. Phases were initialized randomly,
and both the algorithms were finished at 1000 iterations. As themother wavelet, we used the log-normal wavelet [4], which is defined in the Fourier-transformed domain:()
(log ω)2



exp
−
(ω > 0)
4σ2(25)ψ̂(ω) := 0(ω ≥ 0)where ω is an angular frequency and σ is a standard deviation. σ
was set at 0.02 and the analysis frequencies ranged 27.5 to 7040 Hz
with 20 cents interval (i.e. uniform interval in the log-frequency
domain). In the proposed algorithm, we computed the elements
within ±3σ around the central frequencies in the log-frequency
domain. The used computer had the Intel Xeon CPU E31245 (3.3GHz) and a 32 GB RAM.
We employed the perceptual evaluation of speech quality(PESQ) [14] as the evaluation measure for audio quality, which is
the world-standard objective evaluation measure for speech quality. It ranges −0.5 to 4.5 and speech quality is higher as the PESQ2

4

6 8 10 12 14 16 18
Signal length [s]Figure 6: The averaged computation time per iteration with standard errors with respect to various signal lengths.

becomes larger. As an evaluation measure of the computation
speed, the computation time per iteration was used.
5.1.2. Results
The averaged PESQ with standard errors were 4.20 ± 0.08 for the
Irino’s algorithm and 4.1 ± 0.1 for the proposed algorithm. The
result indicates that the speech qualities of the reconstructed signals were high enough for practical use. The difference between
the Irino’s and proposed algorithms was negligible practically1 .
Fig. 6 shows the results for the computation speed with respect
to the signal length, since the computational complexity of the algorithms primarily depends on the signal length. The proposed
algorithm was around 100 times faster than the Irino’s algorithm
in the computation time. For example, the averaged computation
time per iteration by the Irino’s algorithm was around 10 s/iteration
for the 15 s signal. In contrast, that by the proposed algorithm wasaround 0.1 s/iteration.
5.2. Second experiment: Relation between approximation accuracy and audio quality
5.2.1. Experimental conditions
The proposed algorithm includes the approximation, and we next
evaluated the relation between the approximation accuracy and theaudio quality of the reconstructed signals. We used the 5 s from30 s of 102 music audio files with 16 kHz sampling frequency in
the RWC music genre database [15]. As the mother wavelet, thelog-normal wavelet with σ = 0.02 was chosen. The approximation accuracy of the proposed algorithm corresponds to the calculated range by the downsampling step, and we used the elements
within ±Pσ (P = 1, 2, 3, 5) around the central frequencies in the
log-frequency domain. The number of iterations was set at 500 for
the proposed algorithm and at 100 for the Irino’s algorithm. The
used computer had the Intel Core i3-2120 CPU (3.30 GHz) and a
8 GB RAM. The other experimental conditions were the same as
in Sec. 5.1.1.
An evaluation measure for audio quality was the objective differential grade (ODG) by the perceptual evaluation of audio qual1 Audio samples are available at http://hil.t.u-tokyo.ac.
jp/~nakamura/demo/fastCWT.html.

DAFX-5

Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014
4.5
Perceptual evaluationof speech quality

Objective difference grade

Conventional

-0.5
-1
-1.5
-2
-2.5
-3

3.5
3
2.5

-3.5

2

-4

1.5

0

100
200
300
400
The number of Iterations

500

Conventional0

Figure 7: Evolution of the averaged objective difference gradeswith standard errors by perceptual evaluation of audio quality withrespect to the number of iterations for the proposed algorithms
with various approximations ([−Pσ, Pσ] (P = 1, 2, 3, 5)) and the
Irino’s algorithm [7].

20

40
60
80
Computation time [s]100Figure 9: Evolution of the perceptual evaluation of speech quality
and the computation time with respect to the proposed algorithms
with various approximations ([−Pσ, Pσ] (P = 1, 2, 3, 5)) and the
Irino’s algorithm [7].

the Irino’s algorithm. The similar result for the speech signal(fafsc110 in the ATR Japanese speech database A-set [13], the 7s signal) was shown in Fig. 9. Therefore, we conclude that the
proposed algorithm with around P = 3 provides the reconstructedsignals with high audio quality in a reasonable computation time.

Conventional

Objective difference grade

4

-0.5
-1
-1.5
-2

5.3. Demonstration of phase estimation

-2.5
-3
-3.5
-4
0

10

20
30
40
Computation time [s]50

Figure 8: Evolution of the objective difference grades by perceptual evaluation of audio quality with respect to the computation time for the proposed algorithms with various approximations([−Pσ, Pσ] (P = 1, 2, 3, 5)) and the Irino’s algorithm [7].

We demonstrate pitch transposition of acoustic signals to confirm
effectiveness of the proposed algorithm for sound manipulation.
When the analysis frequencies are located uniformly in the logfrequency domain and D0 = D1 = · · · = DL−1 in the proposed
algorithm, we simply shift the components of the CWT spectrograms to the lower or higher analysis frequency components, andthe blank components by the move are filled by zero. However, theshifts cause the mismatches of phases, and the use of the originaland zero phases leads to failure of the pitch transposition, hence
we need to use the phase estimation for synthesizing the pitchtransposed acoustic signals. By the proposed algorithm, we obtained the synthesized signals 3 as we expected.

ity (PEAQ) [16]. It ranges −4 to 0, and the acoustic quality ishigher as the ODG becomes larger.

6. CONCLUSION

Fig. 7 illustrates the averaged ODGs with standard errors. The
ODGs by the proposed algorithms with P = 3, 5 were larger than
−2.0 after 100 iterations, and the results for P = 3, 5 shows highaudio quality2 . The results does not significantly differ from thatby the Irino’s algorithm in audio quality. We can thus say that the
proposed algorithm with around P ≥ 3 reconstructs the acoustic
signals with almost the same audio quality as the Irino’s algorithm.
In a viewpoint of the computation speed, the computation time
becomes shorter as P is smaller. Fig. 8 shows the result for one of
the acoustic signals (RWC-MDB-G-2001 No. 1), and the ODGs
by the proposed algorithms quickly become higher than those by

We have proposed a fast and convergence-guaranteed algorithm of
the phase estimation by using the fast approximate CWT [12]. The
phase estimation problem has been formulated based on the consistency condition, and the iterative algorithm has been derived byapplying the auxiliary function method, which is the same as the
Irino’s algorithm [7]. Furthermore, we show the requirement on
scale factors and mother wavelets for the phase estimation by using
the consistency condition. The experimental results have shownthat the proposed algorithm was about 100 times faster than thealgorithm provided in [7]. The audio quality of the reconstructedsignals for music and speech data was high enough for practicaluse, and the difference between the results by the proposed algorithm and the algorithm provided in [7] was negligible.

2 c.f.) When the used audio signals were converted into MPEG-3 files
with 160 kbps, the averaged ODGs with standard errors were −3.68 ± 0.03.

3 The synthesized signals are available at http://hil.t.
u-tokyo.ac.jp/~nakamura/demo/fastCWT.html.

5.2.2. Results

DAFX-6

Proc. of the 17th Int. Conference on Digital Audio Effects (DAFx-14), Erlangen, Germany, September 1-5, 2014

We plan to combine the phase estimation with source separations for magnitude CWT spectrograms for music acoustic signal
manipulation such as conversions of chords, keys and scales. To
increase convenience, developing the online version of the proposed algorithm is important.

[15] M. Goto, “Development of the RWC Music Database,” inProc. Int. Congress Acoust., 2004, pp. l–553–556.
[16] “ITU-T recommendation BS.1387-1, Perceptual evaluationof audio quality (PEAQ): Method for objective measurements of perceived audio quality,” Sep. 2001.

7. ACKNOWLEDGMENTS
This work was supported by JSPS Grant-in-Aid for Young Scientists B Grant Number 26730100.
8. REFERENCES[1] H. Fletcher, “Auditory patterns,” Rev. Mod. Phys., vol. 12,
pp. 47–61, 1940.
[2] R. D. Patterson, “Auditory filter shape,” J. Acoust. Soc. Am.,
vol. 55, no. 4, pp. 802–809, 2005.
[3] M. N. Schmidt and M. Mørup, “Nonnegative matrix factor
2-D deconvolution for blind single channel source separation,” in Independent Component Analysis and Blind Signal
Separation, pp. 700–707. Springer, 2006.
[4] H. Kameoka, Statistical Approach to Multipitch Analysis,
Ph.D. thesis, The University of Tokyo, Mar. 2007.
[5] M. Muller, D. P. W. Ellis, A. Klapuri, and G. Richard, “Signal processing for music analysis,” IEEE J. Sel. Topics. Signal Process., vol. 5, no. 6, pp. 1088–1110, 2011.
[6] J. P. de León, F. Beltrán, and J. R. Beltrán, “A complex
wavelet based fundamental frequency estimator in singlechannel polyphonic signals,” in Proc. Digital Audio Effects,
2013.
[7] T. Irino and H. Kawahara, “Signal reconstruction from modified auditory wavelet transform,” IEEE Trans. Signal Process., vol. 41, no. 12, pp. 3549–3554, 1993.
[8] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama, “Fastsignal reconstruction from magnitude STFT spectrogram
based on spectrogram consistency,” in Proc. Int. Conf. Digital Audio Effects, Sep. 2010, pp. 397–403.
[9] D. Griffin and J. Lim, “Signal estimation from modified
short-time fourier transform,” IEEE Trans. Acoust., Speech,
Signal Process., vol. 32, no. 2, pp. 236–243, 1984.
[10] D. M. Lopes and P. R. White, “Signal reconstruction fromthe magnitude or phase of a generalised wavelet transform,”in Proc. Eur. Signal Process. Conf., 2000, pp. 2029–2032.
[11] C. Torrence and G. P. Compo, “A practical guide to wavelet
analysis,” B. Am. Meteorol. Soc., vol. 79, no. 1, pp. 61–78,
1998.
[12] H. Kameoka, T. Tabaru, T. Nishimoto, and S. Sagayama,“(Patent) Signal processing method and unit,” Nov. 2008, inJapanese.
[13] 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 Commun., vol. 9, no. 4, pp. 357–363, 1990.
[14] “ITU-T recommendation P.862, Perceptual evaluation of
speech quality (PESQ): an objective method for end-to-end
speech quality assessment of narrow-band telephone networks and speech codecs,” Feb. 2001.

DAFX-7