Characterization of Analytic Wavelet Transforms and a New Phaseless Reconstruction Algorithm

TL;DR

This paper characterizes all wavelets that produce analytic wavelet transforms and introduces a new phase reconstruction algorithm that improves performance and flexibility over previous methods.

Contribution

It provides a comprehensive characterization of wavelets leading to analytic WT and develops a novel phaseless reconstruction algorithm with enhanced performance.

Findings

New phase-magnitude relationships similar to STFT

The reconstruction method outperforms previous techniques

Flexible trade-off between accuracy and complexity

Abstract

We obtain a characterization of all wavelets leading to analytic wavelet transforms (WT). The characterization is obtained as a by-product of the theoretical foundations of a new method for wavelet phase reconstruction from magnitude-only coefficients. The cornerstone of our analysis is an expression of the partial derivatives of the continuous WT, which results in phase-magnitude relationships similar to the short-time Fourier transform (STFT) setting and valid for the generalized family of Cauchy wavelets. We show that the existence of such relations is equivalent to analyticity of the WT up to a multiplicative weight and a scaling of the mother wavelet. The implementation of the new phaseless reconstruction method is considered in detail and compared to previous methods. It is shown that the proposed method provides significant performance gains and a great flexibility regarding accuracy versus complexity. Additionally, we discuss the relation between scalogram reassignment operators and the wavelet transform phase gradient and present an observation on the phase around zeros of the WT.

Tables

Table 1. Table I : Means and standard deviation of spectral convergence for the considered methods and parameter sets.

$\alpha$	method	WPGHI	FBPGHI	R-FGLIM	W-FGLIM
30	mean	$-33.9512$	$-18.2916$	$-28.9164$	$-40.2932$
	std	$7.7632$	$2.0358$	$3.2641$	$7.0526$
300	mean	$-36.6393$	$-26.8642$	$-29.3628$	$-42.7550$
	std	$7.9517$	$3.4095$	$4.0747$	$9.0264$
3000	mean	$-38.5527$	$-34.8240$	$-32.0108$	$-44.9260$
	std	$6.7961$	$5.5137$	$5.4821$	$10.6349$

Table 2. Table II : Means and standard deviation of spectral convergence for WPGHI and WPGHI-FGLIM for redundancies K / a = 30 , 10 , 5 , 3 𝐾 𝑎 30 10 5 3 K/a=30,10,5,3 , with α = 1000 𝛼 1000 \alpha=1000 .

$K/a_d$	method	WPGHI	W-FGLIM
30	mean	$-38.2012$	$-44.1790$
	std	$6.6897$	$9.8989$
10	mean	$-32.5835$	$-40.7522$
	std	$6.5999$	$6.9075$
5	mean	$-26.8400$	$-36.4934$
	std	$5.9723$	$5.6917$
3	mean	$-20.7873$	$-29.7372$
	std	$5.1843$	$6.6363$

Equations

$${\xi }^{\frac{\alpha -1}{2}}{e}^{-2\pi \gamma \xi }{e}^{i\beta \log \xi }$$

$$W_{\psi }s(x,y)=⟨s,{\mathbf{T}}_x{\mathbf{D}}_y\psi ⟩=\frac{1}{\sqrt{y}}{\int }_{\mathbb{R}}s(t)\overline{\psi \left(\frac{t-x}{y}\right)}dt,$$

$$\nabla \log (M_{\psi }^s)(x,y)=\left(\begin{array}{c}0\\ \frac{1}{2y}\end{array}\right)-\mathrm{Re}\left(\begin{array}{c}\frac{W_{{\psi }^{\prime }}s(x,y)}{y{W}_{\psi }s(x,y)}\\ \frac{W_{(\mathbf{T}\psi {)}^{\prime }}s(x,y)}{y{W}_{\psi }s(x,y)}\end{array}\right)$$

$$\nabla {\varphi }_{\psi }^s(x,y)=-\mathrm{Im}\left(\begin{array}{c}\frac{W_{{\psi }^{\prime }}s(x,y)}{y{W}_{\psi }s(x,y)}\\ \frac{W_{(\mathbf{T}\psi {)}^{\prime }}s(x,y)}{y{W}_{\psi }s(x,y)}\end{array}\right).$$

$$\begin{array}{rl}h:\{z\in \mathbb{C}:\mathrm{Im}(z)>0\}& \to \mathbb{C}\\ x+iy& \mapsto f(x,y)W_{\psi }s(x-aby,by)\end{array}$$

$$\hat{\psi }(\xi )=c{\xi }^{\frac{\alpha -1}{2}}{e}^{-2\pi \gamma \xi }{e}^{i\beta \log \xi }$$

$$x+iy\mapsto y^{-\frac{\alpha}{2}}e^{i\beta\log y}W_{\psi}s\bigg{(}x-\frac{\operatorname{Im}(\gamma)}{\operatorname{Re}(\gamma)}y,\frac{y}{\operatorname{Re}(\gamma)}\bigg{)}.$$

$$\hat{\psi }(\xi )=c_{\mathrm{n}}(-\xi {)}^{\frac{{\alpha }_{\mathrm{n}}-1}{2}}{e}^{2\pi {\gamma }_{\mathrm{n}}\xi }{e}^{i{\beta }_{\mathrm{n}}\log (-\xi )}$$

$$\frac{\partial}{\partial x}\phi_{\psi}^{s}=\frac{\alpha}{2y\operatorname{Re}(\gamma)}-\frac{\frac{\partial}{\partial y}\log\big{(}M_{\psi}^{s}\big{)}}{\operatorname{Re}(\gamma)}+\frac{\operatorname{Im}(\gamma)\frac{\partial}{\partial x}\log\big{(}M_{\psi}^{s}\big{)}}{\operatorname{Re}(\gamma)}$$

$$\frac{\partial }{\partial y}{\varphi }_{\psi }^s$$

$$\displaystyle\quad-\frac{\operatorname{Im}(\gamma)\frac{\partial}{\partial y}\log\big{(}M_{\psi}^{s}\big{)}}{\operatorname{Re}(\gamma)}.$$

$$\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x,y)=-\frac{\partial }{\partial y}\log (M_{\psi }^s)(x,y)+\frac{\alpha }{2y}$$

$$\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x,y)=\frac{\partial }{\partial x}\log (M_{\psi }^s)(x,y)-\frac{\beta }{y}.$$

$$\frac{{\partial }^2}{\partial x^2}\log (M_{\psi }^s)(x,y)+\frac{{\partial }^2}{\partial y^2}\log (M_{\psi }^s)(x,y)=-\frac{\alpha }{2y^2}$$

$$\frac{{\partial }^2}{\partial x^2}{\varphi }_{\psi }^s(x,y)+\frac{{\partial }^2}{\partial y^2}{\varphi }_{\psi }^s(x,y)=\frac{\beta }{y^2}.$$

$${\tilde{W}}_{\psi }s(x,\xi )=⟨s,{\mathbf{T}}_x{\tilde{\mathbf{D}}}_{{\xi }_{\mathrm{b}}/\xi }\psi ⟩=\sqrt{\frac{\xi }{{\xi }_{\mathrm{b}}}}{W}_{\psi }s(x,{\xi }_{\mathrm{b}}/\xi ).$$

$$\frac{\partial }{\partial x}{\tilde{\varphi }}_{\psi }^s(x,\xi )$$

$$=\frac{4\pi {\xi }^2}{\alpha -1}\frac{\partial }{\partial \xi }\log ({\tilde{M}}_{\psi }^s)(x,\xi )+2\pi \xi$$

$$\frac{\partial }{\partial \xi }{\tilde{\varphi }}_{\psi }^s(x,\xi )$$

$$=-\frac{\alpha -1}{4\pi {\xi }^2}\frac{\partial }{\partial x}\log ({\tilde{M}}_{\psi }^s)(x,\xi )+\frac{\beta }{\xi },$$

$$\frac{\partial }{\partial x}{\tilde{\varphi }}_{\psi }^s(x,\xi )$$

$$\frac{\partial }{\partial \xi }{\tilde{\varphi }}_{\psi }^s(x,\xi )$$

$$\mathrm{sgn}(J_{x_0,y_0})\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x_0,y_0+\epsilon )$$

$$\mathrm{sgn}(J_{x_0,y_0})\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x_0,y_0-\epsilon )$$

$$\mathrm{sgn}(J_{x_0,y_0})\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x_0+\epsilon ,y_0)$$

$$\mathrm{sgn}(J_{x_0,y_0})\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x_0-\epsilon ,y_0)$$

$$\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x_0+\epsilon ,y_0)=\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x_0-\epsilon ,y_0)$$

$$\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x_0,y_0+\epsilon )=\underset{\epsilon \to 0}{\lim }\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x_0,y_0-\epsilon )$$

$$(x,y)\mapsto \left(x+\frac{y^2\frac{\partial }{\partial y}{\varphi }_{\psi }^s(x,y)}{{\xi }_{\mathrm{b}}},\frac{{\xi }_{\mathrm{b}}}{\frac{\partial }{\partial x}{\varphi }_{\psi }^s(x,y)}\right).$$

$$(x,y)\mapsto \left(x-\frac{y\mathrm{Im}\left(\frac{W_{(\mathbf{T}\psi {)}^{\prime }}s(x,y)}{W_{\psi }s(x,y)}\right)}{{\xi }_{\mathrm{b}}},-\frac{y{\xi }_{\mathrm{b}}}{\mathrm{Im}\left(\frac{W_{{\psi }^{\prime }}s(x,y)}{W_{\psi }s(x,y)}\right)}\right).$$

$$(x,y)\mapsto \left(x+y\mathrm{Re}\left(\frac{W_{\mathbf{T}\psi }s(x,y)}{W_{\psi }s(x,y)}\right),-\frac{y{\xi }_{\mathrm{b}}}{\mathrm{Im}\left(\frac{W_{{\psi }^{\prime }}s(x,y)}{W_{\psi }s(x,y)}\right)}\right).$$

$$\begin{array}{rl}\mathrm{Re}\left(\frac{W_{\mathbf{T}{\psi }_{\mathrm{G}}}s(x,y)}{W_{{\psi }_{\mathrm{G}}}s(x,y)}\right)& =-\mathrm{Re}\left(\frac{W_{({\psi }_{\mathrm{G}}{)}^{\prime }}s(x,y)}{W_{{\psi }_{\mathrm{G}}}s(x,y)}\right)\\ & =y\frac{\partial }{\partial x}\log (M_{{\psi }_{\mathrm{G}}}^s)(x,y)\end{array}$$

$$(x,y)\mapsto \left(x+y^2\frac{\partial }{\partial x}\log (M_{{\psi }_{\mathrm{G}}}^s)(x,y),\frac{{\xi }_{\mathrm{b}}}{\frac{\partial }{\partial x}{\varphi }_{{\psi }_{\mathrm{G}}}^s(x,y)}\right).$$

Full text

I Introduction

Time-frequency and time-scale representations are fundamental tools in many areas of signal analysis and signal processing, ranging from medical data [1, 2], damage or fault detection in materials [3] and machines [4], to image [5, 6] and audio processing [7, 8, 9]. Such representations are usually complex-valued and admit a natural decomposition into magnitude and phase components, with the phase containing crucial information about the analyzed signal. However, phase information is often discarded in favor of the magnitude information, from which, supposedly, the desired information is more readily obtained. On the other hand, whenever synthesis from the representation coefficients is desired, the phase is crucial for the quality of the synthesized signal.

In recent years, the problems of phase retrieval and phaseless reconstruction for time-frequency and time-scale dictionaries have attracted considerable attention, leading to theoretical results for the feasibility of phase retrieval in different contexts [10, 11, 12, 13, 14, 15] and various algorithms [16, 17, 18, 19, 20, 21] (see also the survey [22]). Many of these algorithms, e.g., [20], attempt to construct, iteratively or directly, an appropriate phase to match the magnitude-only representation coefficients, before finally performing a regular synthesis step.

In many applications, in particular in audio signal processing, phaseless reconstruction in time-frequency and time-scale dictionaries as proposed in, e.g., [20, 23], is successfully applied. These applications include, but are not restricted to audio and speech synthesis [24, 25, 26], source separation [27, 28, 29], as well as pitch and time-scale modifications [30, 31, 32]. In these scenarios, either reconstruction from a phaseless representation is desired or the given phase has been invalidated in the course of processing and must be replaced.

Despite a relevant research activity regarding the phase behavior of time-frequency and time-scale representations, with particular incidence in the short-time Fourier transform (STFT) [33, 34, 35], applications that consider and/or modify the phase are relatively scarce. Prime examples of analysis tools that successfully use phase information are the so-called reassignment and synchrosqueezing methods [36, 37, 38, 39, 40] that deform the signal representation using a vector field obtained from the phase gradient of the representation.

For the STFT with Gaussian generator, the notion that phase and magnitude carry equally important information has been made precise by Portnoff [34] and later by Auger and Flandrin [33]. They show that in this specific case, the phase gradient is completely characterized by the gradient of the (logarithmically scaled) magnitude.

*Contributions: *Our two main contributions are a characterization of all analytic wavelet transforms (WT) and a new method for wavelet phase reconstruction from magnitude-only coefficients. Here, the notion of analytic WTs is used in the sense that the WT is an analytic function of the upper half-plane (up to a signal-independent factor). It is common to use the same terminology for wavelet transforms generated by analytic wavelets [41], i.e., wavelets whose Fourier transform vanishes at negative frequencies. We discuss the connection between these two notions in detail.

The result characterizing all analytic wavelet transforms shows that the most general wavelet leading to an analytic transform has a Fourier transform given by

[TABLE]

for positive frequencies and we assume that the Fourier transform vanishes for negative frequencies. The appearance of the $\beta$-dependent hyperbolic chirp in (1) may be surprising since, to our knowledge, only the case $\beta=0$ has been associated with analytic functions in the literature (see [42]). It is commonly referred to as Cauchy wavelet and historically associated with affine coherent states. However, the more general wavelets given by (1) have been shown to minimize a time-scale counterpart of Heisenberg uncertainty and are also known as “Klauder wavelets” [43]. The problem of characterizing all analytic WTs has been considered and partially solved in [44], where it was also shown that the Gaussian is essentially the only window leading to an analytic STFT.

We first discuss several aspects of the WT phase, loosely following the structure of [33]. Specifically, we express the phase and log-magnitude derivatives in terms of pointwise quotients of $3$ different WTs. The mother wavelets used in these transforms can be derived directly from the mother wavelet of the original WT. These expressions can be used, e.g., to estimate the local group delay and instantaneous scale (or frequency). Furthermore, they impose conditions on the wavelets leading to analytic wavelet transforms, and this is then used to show that the class of wavelets satisfying this condition are the generalized Cauchy wavelets (1). The corresponding Cauchy Riemann (CR) equations provide a relation between the phase gradient and the (log-)magnitude gradient. Following [35], we also discuss the singular behavior of the phase close to zeros of the WT and the implications of our results for wavelet reassignment and ridge analysis. Finally, we discuss the relation between analytic wavelets, i.e., wavelets that can be extended to analytic functions on the upper half plane, and wavelets resulting in an analytic WT.

In the second part of the contribution, we implement a method for reconstruction from magnitude-only wavelet coefficients. More specifically, we use a discrete approximation of the derived phase-magnitude relations to modify the phase gradient heap integration algorithm [45, 46, 47]. The method is evaluated using Cauchy wavelets of multiple orders, obtaining favorable results. In the evaluation, we further examine (uniform) decimation and the number of scales at which the WT is sampled, controlling the transform redundancy.

Structure of the Paper: In Section II, we present our expressions for the derivatives of the the log-magnitude and the phase of the WT. These expressions are used in Section III to obtain our result characterizing the wavelets that lead to analytic WTs. Furthermore, we prove direct relations between the log-magnitude and phase derivatives. In Section IV, we provide a time-frequency interpretation of the corresponding WT. Further applications to pole behavior and scalogram reassignment of the expression provided in Section II are presented in Sections V-A and V-B, respectively. In Section V-C, we discuss the relation between analytic wavelets and our analytic WT. Finally, we apply the phase-magnitude relations to the problem of phaseless reconstruction. In Section VI, we formally introduce the discretization of our results and describe the phaseless reconstruction algorithm; in Section VII, we conduct experiments on real data and compare our method to previous approaches toward phaseless reconstruction.

II

Log-Magnitude and Phase Derivatives of the Wavelet Transform

Fix a function $\psi\in\mathbf{L}^{2}(\mathbb{R})$ such that its Fourier transform $\widehat{\psi}$ vanishes almost everywhere on $\mathbb{R}^{-}$. The continuous WT (CWT) of a function (or signal) $s\in\mathbf{L}^{2}(\mathbb{R})$ with respect to the mother wavelet $\psi$ is defined as

[TABLE]

for all111Although we restrict here to positive scales, there is no technical obstruction to allowing $y\in\mathbb{R}\setminus\{0\}$. $x\in\mathbb{R}$, $y\in\mathbb{R}^{+}$. Here, $\mathbf{T}_{x}$ and $\mathbf{D}_{y}$ denote the translation and dilation operators, respectively, given by $(\mathbf{T}_{x}s)(t)=s(t-x)$, and $(\mathbf{D}_{y}s)(t)=y^{-1/2}s(t/y)$ for all $t\in\mathbb{R}$.

The CWT can be represented in terms of its magnitude $M_{\psi}^{s}:=|W_{\psi}s|\geq 0$ and phase $\phi^{s}_{\psi}:=\arg(W_{\psi}s)\in\mathbb{R}$ as usual. With this convention, $\log(W_{\psi}s)=\log(M_{\psi}^{s})+i\phi^{s}_{\psi}$. This straightforward relation is the basis of the following expressions for the partial derivatives of the log-magnitude and phase components, derived in Appendix A.

Theorem 1**.**

Let $\psi\in\mathbf{L}^{2}(\mathbb{R})$ with $\widehat{\psi}(\xi)=0$ for $\xi<0$ and assume that $\psi$ is continuously differentiable with $\psi^{\prime},\mathbf{T}\psi^{\prime}\in\mathbf{L}^{2}(\mathbb{R})$, where $\mathbf{T}$, without subscript, denotes the time-weighting operator $(\mathbf{T}s)(t)=ts(t)$. Then, for all $x\in\mathbb{R}$ and $y\in\mathbb{R}^{+}$ satisfying $M_{\psi}^{s}(x,y)>0$,

[TABLE]

and

[TABLE]

The partial phase derivatives of the WT are often related to the local instantaneous scale, as well as the local group delay, of the analyzed signal, although at least the latter notion is not entirely clear for WTs, see Section V-B. The second order derivatives, which can be obtained similar to the first order derivatives in Appendix A, describe the variation of these quantities across phase space and thus proved useful as well, see [48, 49, 50].

Formulas (3)–(4) can be used for the computation of the partial derivatives by using efficient implementations of the WT. Furthermore, the possible accuracy of direct numerical differentiation of the wavelet transform is limited if the WT can be computed only at certain positions and not everywhere, e.g., in the presence of decimation in either coordinate. In contrast, even if a closed form expression for the derivative of the (time-weighted) wavelet is not available, numerical differentiation of the wavelet is not limited by this constraint.

III The Phase-Magnitude Relationship

In general, the observations in Section II do not yield a direct connection between the partial derivatives of the log-magnitude and phase components. However, if the WT is analytic, we can characterize the phase gradient by the log-magnitude gradient. Based on the Cauchy-Riemann (CR) equations, we can construct conditions on the mother wavelet $\psi$ such that the WT of any signal is an analytic function. Similar to the analysis of analytic STFTs studied in [44] and the partial study of analytic wavelet transforms in the same contribution, we allow for an $(x,y)$-dependent factor $f(x,y)$ that is independent of the signal $s$. This factor can easily be accounted for when applying results from complex analysis to the transformed signal and leads to a significantly less restrictive class of analyticity-inducing wavelets. Furthermore, we want the class of analyticity-inducing wavelets to be invariant under the natural transforms associated with the WT, namely dilation and translation. Thus, we also allow for a constant dilation by $b\in\mathbb{R}^{+}$ and a translation specified by $a\in\mathbb{R}$ in our analysis.

Theorem 2**.**

Let $\psi\in\mathbf{L}^{2}(\mathbb{R})$ with $\widehat{\psi}(\xi)=0$ for $\xi<0$. There exist constants $a\in\mathbb{R}$, $b\in\mathbb{R}^{+}$ and a $\mathcal{C}^{\infty}$ function $f\colon\mathbb{R}\times\mathbb{R}^{+}\to\mathbb{C}$ with $f(x,y)\neq 0$ such that

[TABLE]

is analytic for all $s\in\mathbf{L}^{2}(\mathbb{R})$, if and only if

[TABLE]

for all $\xi>0$ and some constants $c\in\mathbb{C}$, $\alpha>-1$, $\beta\in\mathbb{R}$, and $\gamma\in\mathbb{C}$ with $\operatorname{Re}(\gamma)>0$.

A proof of the theorem is provided in Appendix B. The wavelets $\psi$ specified by (6) are known as Klauder wavelets [43] and are a minor generalization of Cauchy wavelets $\psi^{(\alpha)}$ [42], which are recovered for the choice $\beta=0$ and $\gamma=1$. Among other effects, modifying $\beta$ results in a proportional shift of the temporal concentration of the mother wavelet away from time zero. Furthermore, a change in $\gamma$ results only in a scale change, dependent on $\operatorname{Re}(\gamma)$, and a time shift, dependent on $\operatorname{Im}(\gamma)$. Disregarding the constant factor $c\in\mathbb{C}$ in (6), we denote the generalized Cauchy wavelets in Theorem 2 by $\psi^{(\alpha,\beta,\gamma)}$ or $\psi^{(\alpha,\beta)}$, if $\gamma=1$.

If $\psi$ is given by (6), a corresponding choice for $h$ being analytic is $f(x,y)=y^{-\frac{\alpha}{2}}e^{i\beta\log y}$, $a=\operatorname{Im}(\gamma)$, and $b=1/\operatorname{Re}(\gamma)$, i.e.,

[TABLE]

We note that the wavelets are admissible only for $\alpha>1$.

Theorem 2 can easily be modified to allow for wavelets $\psi$ where $\widehat{\psi}$ does not vanish for negative frequencies. In this case, $\widehat{\psi}$ must satisfy (6) for $\xi>0$ and

[TABLE]

for $\xi<0$ and some constants $c_{\mathrm{n}}$, $\alpha_{\mathrm{n}}$, $\gamma_{\mathrm{n}}$, and $\beta_{\mathrm{n}}$ satisfying the same constraints as $c$, $\alpha$, $\gamma$, and $\beta$, respectively.

Based on the CR equations, we obtain a phase-magnitude relation for the WTs using $\psi^{(\alpha,\beta,\gamma)}$.

Theorem 3**.**

*Let $\psi$ be given by (6) with $c=1$. Then *

[TABLE]

and

[TABLE]

For $\gamma=1$, these relations simplify to

[TABLE]

and

[TABLE]

A proof of the theorem is provided in Appendix C. As a simple consequence of Theorem 3, we obtain for the second order derivatives in the case $\gamma=1$

[TABLE]

and

[TABLE]

For Cauchy wavelets, it is known that the magnitude uniquely determines the phase up to a constant phase factor. Moreover, this statements even holds after decimation in the scale component and in certain discretized settings [10].

IV The Cauchy Wavelet Transform as Time-Frequency Representation

If the mother wavelet $\psi$ is frequency-localized around frequency $\xi_{\mathrm{b}}$, then we can interpret $W_{\psi}s(x,y)$ as a time-frequency measurement at frequency $\xi=\xi_{\mathrm{b}}/y$. In the case of the wavelets $\psi^{(\alpha,\beta)}$, this leads to a particularly convenient form of the phase-magnitude relationship. Here, we consider the unique peak of $\lvert\widehat{\psi^{(\alpha,\beta)}}\rvert$ as center frequency, i.e., $\xi_{\mathrm{b}}=\frac{\alpha-1}{4\pi}$. Additionally, instead of the unitary dilation $\mathbf{D}_{y}$, we consider the dilation $\tilde{\mathbf{D}}_{y}s(t)=y^{-1}s(t/y)$ to define

[TABLE]

Using the relations in Section III, it is easy to derive

[TABLE]

and

[TABLE]

where $\tilde{M}_{\psi}$ and $\tilde{\phi}_{\psi}$ denote the magnitude and phase of $\tilde{W}_{\psi}$, respectively. Interestingly, for standard Cauchy wavelets, i.e., $\beta=0$, these formulas coincide with the phase-magnitude relations for the STFT with a dilated Gaussian (see [45, Section III]) up to the simple change that the constant time-frequency ratio $\lambda$ is replaced by the frequency-dependent term $\frac{\alpha-1}{4\pi\xi^{2}}$. The above form (IV)–(IV) will enable us to adapt the phase reconstruction method presented in [47] to the WT more easily in Section VI below.

The additive term $\beta/\xi$ in (IV) compensates for the fact that $\psi^{(\alpha,\beta)}$ is time-localized around $\frac{2\beta}{\alpha-1}$. In other words, the frequency (or scale) bands in $\tilde{W}_{\psi^{(\alpha,\beta)}}$ (or $W_{\psi^{(\alpha,\beta)}}$) are not time-aligned. Time-alignment can be restored by choosing $\operatorname{Im}(\gamma)=\frac{2\beta}{\alpha-1}$, thus removing the additive term at the cost of introducing directional derivatives of $\log(\tilde{M}_{\psi}^{s})$ in the expression of the phase gradient:

[TABLE]

with $d_{1}(\xi)=\left(\frac{2\beta}{\alpha-1},\frac{4\pi\xi^{2}}{\alpha-1}\right)$ and

[TABLE]

with $d_{2}(\xi)=\left(\frac{\alpha-1}{4\pi\xi^{2}}+\frac{\beta^{2}}{(\alpha-1)\pi\xi^{2}},\frac{2\beta}{\alpha-1}\right)$.

V Derivatives and Analyticity of the Continuous Wavelet Transform—Further Observations

V-A The Phase Around Zeros of the Wavelet Transform

For the STFT with Gaussian window, it was remarked by Auger et al. [33], that the phase has characteristic poles where the STFT is zero and that this fact can be derived from the analyticity of the Bargman transform. In [35], the characteristic pole behavior was proven under weaker conditions, as long as the STFT is smooth enough, more specifically $\mathcal{C}^{2}$ or $\mathcal{C}^{3}$.

In fact, the techniques used therein apply to any complex-valued function of two real variables, as long as its higher-order partial derivatives are continuous. In the case of the WT, this can be ensured by selecting a sufficiently smooth and decaying mother wavelet $\psi$. In particular, if the $k$-th derivative of $\psi$ weighted by $t^{l}$ is square integrable, i.e., \mathbf{T}^{l}\big{(}\frac{\partial^{k}}{\partial t^{k}}\psi\big{)}\in\mathbf{L}^{2}(\mathbb{R}), for all $l,k\in\{0,\ldots,K\}$, then $W_{\psi}s\in\mathcal{C}^{K}(\mathbb{R}\times\mathbb{R}^{+},\mathbb{C})$, for all $s\in\mathbf{L}^{2}(\mathbb{R})$, cf. Appendix D. This implies the following result.

Theorem 4**.**

Let $\psi\in\mathbf{L}^{2}(\mathbb{R})$ and assume that \mathbf{T}^{l}\big{(}\frac{\partial^{k}}{\partial t^{k}}\psi\big{)}\in\mathbf{L}^{2}(\mathbb{R}), for all $l,k\in\{0,\ldots,2\}$. If $W_{\psi}s(x_{0},y_{0})=0$ and the Jacobian determinant $J_{x_{0},y_{0}}:=det(D(W_{\psi}s))(x_{0},y_{0})$ of $W_{\psi}s$ at $(x_{0},y_{0})$ is nonzero, then for $\varepsilon>0$ converging to [math]

[TABLE]

and

[TABLE]

If even \mathbf{T}^{l}\big{(}\frac{\partial^{k}}{\partial t^{k}}\psi\big{)}\in\mathbf{L}^{2}(\mathbb{R}), for all $l,k\in\{0,\allowbreak\ldots,\allowbreak 3\}$, then the limits

[TABLE]

and

[TABLE]

exist and are finite.

After noting that the assumptions imply that $W_{\psi}s\in\mathcal{C}^{2}(\mathbb{R}\times\mathbb{R}^{+},\mathbb{C})$ (or $\mathcal{C}^{3}(\mathbb{R}\times\mathbb{R}^{+},\mathbb{C})$), the proof of the above result is identical to the proofs of [35, Theorem 4.7–4.9], which only rely on continuous differentiability locally.

V-B Scalogram Reassignment and Ridge Points

Reassignment is a technique for sharpening time-frequency and time-scale representations [36]. The reassignment map is a vector field that is used to deform the representation of choice and depends on the representation and the input signal. It has been derived using different methods, e.g., as a constant phase deformation derived from group theoretical properties [37, 38] or by a center of gravity argument relying on the Wigner distribution [36]. In the case of the spectrogram, both notions can be shown to lead to the same reassignment map. However, this is no longer true for the scalogram, i.e., the squared modulus of the WT.

The reassignment map given in [38] relies only on the phase gradient. It is defined as

[TABLE]

Using (4), the map can be rewritten leading to

[TABLE]

Similar to the expressions provided in [36], (26) can prove useful for the efficient calculation of reassigned scalograms. In particular, a direct computation of the phase gradient from samples of the phase may be quite inaccurate. If the phase-magnitude relations (9) and (10) are satisfied, then inserting these relations into the reassignment map (25) allows scalogram reassignment from the scalogram itself.

In [36], a different reassignment map is given by222Note that the equality $\psi=h(-\bullet)\exp(i\xi_{\mathrm{b}}(\bullet))$ converts between the different WT conventions used.

[TABLE]

The two reassignment maps coincide in the second coordinate, which can be considered a notion of local instantaneous scale. However, in the first coordinate, which is often considered an estimate of the local group delay, they are quite different. In particular, there seems to be no connection between the first coordinate of (27) and partial derivatives of the WT $W_{\psi}s$ in general. However, assuming $\mathbf{T}\psi=d_{0}\psi+d_{1}\psi^{\prime}$ for some constants $d_{0},d_{1}\in\mathbb{C}$, enables an expression of the first coordinate in (27) as a linear combination of the partial derivatives of $\log(M_{\psi}^{s})$ and $\phi^{s}_{\psi}$. A mother wavelet satisfying this differential equation is, e.g., the Gabor wavelet $\psi_{\mathrm{G}}(t)=e^{-t^{2}/2+i\xi_{\mathrm{b}}t}$, with $\mathbf{T}\psi_{\mathrm{G}}=i\xi_{\mathrm{b}}\psi_{\mathrm{G}}-(\psi_{\mathrm{G}})^{\prime}$. In this case, we can rewrite

[TABLE]

and therefore (27) becomes

[TABLE]

Note that the Gabor wavelet does not vanish at negative frequencies. However, the proof of Theorem 1 does not rely on this property, such that the derivation of (29) from (27) remains valid.

The expression of the first coordinate in (29) is not too surprising, since $\mathbf{D}_{y}\psi_{\mathrm{G}}$ is simply a dilated, modulated Gaussian, and thus this expression could also be obtained using the reassignment operators and phase-magnitude relationship for the Gaussian STFT, as observed in [33].

Similar to reassignment, wavelet ridge analysis [41, 51] attempts to identify a skeleton of essential time-scale positions in the wavelet transform. The notion of magnitude ridge points is defined as the points $(x,y)$, such that $\frac{\partial}{\partial y}\log(y^{-\frac{1}{2}}M_{\psi}^{s})(x,y)=0$ and $\frac{\partial^{2}}{\partial y^{2}}\log(y^{-\frac{1}{2}}M_{\psi}^{s})(x,y)<0$. Similarly, phase ridge points are defined as the points $(x,y)$, such that $\frac{\partial}{\partial x}\phi^{s}_{\psi}(x,y)-\frac{\xi_{\mathrm{b}}}{y}=0$ and $\frac{\partial}{\partial y}\left(\frac{\partial}{\partial x}\phi^{s}_{\psi}(x,y)-\frac{\xi_{\mathrm{b}}}{y}\right)>0$. It is straightforward to verify that Theorems 2 and 3 imply that the phase and magnitude ridge points coincide if and only if $\psi=c\psi^{(\alpha,\beta,\gamma)}$, with $c\neq 0$, $\alpha>-1$, $\beta\in\mathbb{R}$ and $\gamma\in\mathbb{R}^{+}$.

V-C Analytic Wavelets and the Analytic Wavelet Transform

To prevent confusion related to other works on WTs, we want to point out the connection between our analytic WT and the WT using analytic wavelets [41, 52]. A wavelet is called analytic if it vanishes almost everywhere on $\mathbb{R}^{-}$, i.e., $\widehat{\psi}(\xi)=0$ for $\xi<0$. The reason for this terminology is that these wavelets can be extended to an analytic function on the upper half-plane by the Paley-Wiener theorem. Furthermore, the WT using an analytic wavelet $\psi$ and at a fixed scale $y_{0}$ also has the property that it can be extended to an analytic function on the upper half-plane in a complex variable $w$, i.e., the function $x\mapsto W_{\psi}s(x,y_{0})$ can be extended to an analytic function

[TABLE]

Now, for $\psi=\psi^{(\alpha,\beta,\gamma)}$, the function in (7) is also analytic. The two functions in (7) and (30) coincide (up to a constant) for $y=\operatorname{Re}(\gamma)y_{0}$ and $x=w+\operatorname{Im}(\gamma)y_{0}$ for all $w\in\mathbb{R}$. More specifically,

[TABLE]

for all $w\in\mathbb{R}$ and where both sides are analytic functions in $w$ on the upper half-plane. Thus, they have to coincide everywhere and we see that $h$ describes the analytic continuation $W^{(\textrm{a})}_{\psi}s(w,y_{0})$ of an arbitrary scale up to some constants and shifts. In particular, for the case of Cauchy wavelets $\psi^{(\alpha)}$, we have

[TABLE]

and the analytic continuation at a given scale only differs by a constant multiple from the analytic continuation at any other scale. This equivalence of all analytic continuations (30) is unique to the wavelets $\psi^{(\alpha,\beta,\gamma)}$.

VI Application—Phaseless Reconstruction for the Discrete Continuous Wavelet Transform

In the following, we propose and evaluate a method for signal reconstruction from magnitude-only wavelet coefficients. More specifically, the proposed algorithm computes a phase estimate from the given magnitude-only coefficients. After combining the magnitude-only coefficients with the estimated phase, the wavelet transform must be inverted to obtain a time-domain signal. Here, any method that implements reconstruction from wavelet coefficients can be used.

Since arbitrary dilations cannot be naturally transferred to the discrete domain, discrete implementations of the WT can be quite different from each other, see [53, 54, 55, 56, 57] and references therein. For illustrative purposes and to clarify notation, we shortly sketch an implementation of the discrete WT that follows [58, 59, 60] closely, with some of the modifications introduced in [9]. In particular, we mimic the dilation operator by sampling the continuous frequency response of the mother wavelet $\psi\in\mathbf{L}^{2}(\mathbb{R})\cap\mathbf{L}^{1}(\mathbb{R})$ at the appropriate density. The following description reflects the implementation used in our experiments, see Section VII. We use the terminology discrete continuous WT (DCWT) to distinguish this type of discrete WT from methods based on a multiresolution analysis and wavelet bases [53, 57], that are commonly known as discrete wavelet transform (DWT).

VI-A Discrete Continuous Wavelet Transform

We will denote discretizations of continuous signals by brackets, e.g., the discretized signal $s_{\mathrm{d}}[l]\in\mathbb{C}$ for $l\in\{1,\dots,L\}$ and some $L\in\mathbb{N}$. In this discrete domain, the translation operator acts circularly, i.e., $s_{\mathrm{d}}[l-m]$ is interpreted as $s_{\mathrm{d}}[\mathop{\operatorname{mod}}(l-m,L)]$. Assuming the sampling rate $\xi_{\mathrm{s}}$, the wavelet at scale $y=\xi_{\mathrm{b}}/\xi$ is derived from the frequency response of the mother wavelet $\psi\in\mathbf{L}^{2}(\mathbb{R})\cap\mathbf{L}^{1}(\mathbb{R})$ as

[TABLE]

for $k\in\{-\lfloor L/2\rfloor,\dots,\lceil L/2\rceil-1\}$. Naturally, only a finite range of scales can be considered before the wavelet deteriorates either due to the sampling density being too coarse ($y$ large) or its bandwidth approaching $\xi_{\mathrm{s}}$ ($y$ small). Hence, in order to cover the entire frequency range, we introduce an additional low-pass function in the style of [9, Section 3.1.2].

The entire wavelet system is characterized by the minimum scale $y_{\mathrm{m}}\in\mathbb{R}^{+}$, the scale step333We choose the commonly used geometric spacing of center frequencies, but the proposed phase reconstruction method remains valid for any center frequency spacing. $2^{1/B}$, with $B\in\mathbb{R}^{+}$, the number of scales $K\in\mathbb{N}$, and the decimation factor $a_{\mathrm{d}}\in\mathbb{N}$, with $a_{\mathrm{d}}|L$. The corresponding scaled and shifted wavelets are given as

[TABLE]

for $k\in\{0,\dots,K-1\}$ and $n\in\{0,\dots,L/a_{\mathrm{d}}-1\}$. A plateau function $P_{\mathrm{lp}}\in\mathbb{C}^{L}$, centered at [math], specifies the low-pass function as

[TABLE]

where

[TABLE]

An analysis with the constructed system yields $LK/a_{\mathrm{d}}$ complex-valued coefficients for the wavelet scales and additional $L/a_{\mathrm{d}}$ real-valued coefficients for the low-pass function, for a total redundancy of $(2K+1)/a_{\mathrm{d}}$ when analyzing real-valued signals. With a slight abuse of terminology, we will from now on refer to the proportional quantity $K/a_{\mathrm{d}}$ as the redundancy.

If $\psi$ is smooth and $a_{\mathrm{d}},1/B$ are small enough, then the results in [58] imply that the DCWT is invertible. Inversion can be achieved by interpreting the wavelet transform as a filter bank analysis and invoking the frame theory of uniform filter banks [61, 62, 63, 64] to compute a dual filter bank synthesis. This can be done either directly using dual filters $\widetilde{\psi}_{k}$, or iteratively [65, 9] using conjugate gradient iterations. The consideration of general uniform filter banks is necessary: It is not always possible to find a dual filter bank, which is required to achieve perfect reconstruction, with wavelet structure. Nonetheless, the dual filter bank shares the number of channels $K+1$ and the decimation factor $a$ of the wavelet analysis.

VI-B Application to Phaseless Reconstruction

For the phase-magnitude relations presented in Theorem 3 to hold, we have to assume a wavelet $\psi^{(\alpha,\beta,\gamma)}$ as in (6). In particular, we will restrict to the case $\gamma=1$ for simplicity and drop the superscript $(\alpha,\beta)$ for notational convenience. For the WT and the phase-magnitude relations, we will use the convention introduced in Section IV. The generalization to the full class of wavelets described by (6) is straightforward.

Note that Theorem 3 provides only the phase-derivative and indeed reconstruction can at best be expected to be accurate up to a global phase factor. Furthermore, the reconstruction quality is expected to be worse for low-magnitude areas and thus the proposed algorithm only reconstructs the phase down to a certain magnitude-threshold. Coefficients below that threshold are expected to have little effect on the synthesis and can thus be assigned a random phase. As a consequence of the local, adaptive integration scheme, the reconstructed phase is in fact only expected to be consistent locally with changes by a constant phase factor between these local components. On audio signals, such as the chosen corpus of test data, this change is not expected to have notable perceptual effects. Nonetheless, it is visible in the phase difference (between original and reconstructed phase) in Fig. 3. At low redundancy, which is not well-suited for phase reconstruction in general, the phase distortion may become more severe (see the lower right corner in Fig. 3), sometimes leading to perceivable distortion.

Assume that the continuous-time signal $s$ is approximately band- and time-limited on $[0,\xi_{\mathrm{s}})$ and $[0,L/\xi_{\mathrm{s}})$, respectively. Then, with $s_{\mathrm{d}}[l]=s(l/\xi_{\mathrm{s}})$, for $l\in\{0,\dots,L-1\}$, $a_{\mathrm{d}}=a\xi_{\mathrm{s}}\in\mathbb{N}$, and $\xi_{k}=2^{-k/B}\xi_{\mathrm{b}}/y_{\mathrm{m}}$, we obtain the approximation

[TABLE]

Note that, after taking the logarithmic derivative of (36), the normalization by $\xi_{\mathrm{s}}$ becomes irrelevant.

Hence, we have

[TABLE]

and

[TABLE]

Here, $\Delta_{n}$ and $\Delta_{k}$ are appropriate discrete differentiation schemes. For $\Delta_{n}$, we can use centered differences, i.e.,

[TABLE]

The sampling step in the scale coordinate changes depends on $k$ and weighted centered differences can be used:

[TABLE]

Now, from $\Delta^{\tilde{\phi},x,s}_{\psi}$ and $\Delta^{\tilde{\phi},\xi,s}_{\psi}$, an estimate of the phase of $\tilde{W}_{\psi}s$ at the sampling points $\{(na,\xi_{k})\}_{n,k}$ can be obtained using a quadrature rule considering the variable sampling intervals. That even simple $1$-dimensional trapezoidal quadrature provides satisfactory results is illustrated by our experiments, see Section VII.

The integration itself can be performed by a slightly modified Phase Gradient Heap Integration (PGHI) algorithm [45, 66, 47], see Algorithm 1, using, e.g., the following integration rule on the set of neighbors $(n_{\mathrm{n}},k_{\mathrm{n}})\in\mathcal{N}_{n,k}:=\{(n\pm 1,k),(n,k\pm 1)\}$ of $(n,k)$

[TABLE]

When inserting (39) and (VI-B) into (41), the absolute scale of the center frequencies $\xi_{k}$ and sampling rate $\xi_{\mathrm{s}}$ becomes unimportant and only their ratio enters the quadrature (41). Hence, by considering relative frequencies $\xi_{k}/\xi_{\mathrm{s}}$, the algorithm is valid independent of the assumed sampling rate.

Once the phase estimate $(\tilde{\phi}^{s}_{\psi})_{\text{est}}$ has been computed, it is combined with the magnitude by $W_{s,\text{est}}:=M_{s}e^{i(\tilde{\phi}^{s}_{\psi})_{\text{est}}}$. Subsequently, a time-domain signal can be synthesized as usual, e.g., using a dual filter bank.

VII Experiments

To test and evaluate the proposed method, we performed two experiments, described and discussed below. Both experiments were run on the first $5$ seconds of all $70$ test signals from the Sound Quality Assessment Material recordings for subjective tests provided by the European Broadcasting Union (SQAM database) [68]. For wavelet analysis and synthesis, we used the filter bank methods in the open source Large Time-Frequency Analysis Toolbox (LTFAT [69], http://ltfat.github.io/), where our implementation of Wavelet Phase Gradient Heap Integration (WPGHI) is available by using the ’wavelet’ flag in filterbankconstphase. A function to generate the wavelet filters and scripts for generating the individual experiments and figures are provided on the manuscript website http://ltfat.github.io/notes/053/, where the resulting audio files for all experiment conditions can be found as well. Experimental conditions were restricted to classical Cauchy wavelets, i.e., $\beta=0$ and $\gamma=1$.

Thus, the WT parameters used in the experiments are $(\alpha,a_{\mathrm{d}},K)$, where $\alpha$ is the order of the Cauchy wavelet, $a_{\mathrm{d}}$ is the decimation step and $K$ is the number of frequency channels (or scales) used before adding the lowpass filter. As quantitative error measure, we employ (wavelet) spectral convergence [70], i.e., the relative mean squared error (in dB) between the wavelet coefficient magnitude of the target signal $s_{\mathrm{t}}$ and the proposed solution $s_{\mathrm{p}}$:

[TABLE]

It should be noted that the wavelet coefficient magnitude in the above formula was computed using the same parameter set $(\alpha,a_{\mathrm{d}},K)$ for which phaseless reconstruction was attempted.444Although spectral convergence is in some cases sensitive to parameter changes, preliminary tests showed that, for fixed $\alpha$, comparable spectral convergence is achieved with respect to representations with varying $a_{\mathrm{d}}$ and $K$. On the other hand, when $\alpha$ is changed as well, then the value of $SC(s_{p},s_{t})$ may change dramatically, such that comparing the results across different choices of $\alpha$ may be misleading.

VII-A Experiment I—Comparison to Previous Methods

To study the performance of the proposed algorithm in comparison with previous methods for phaseless recovery from wavelet coefficients, we selected three settings of the WT parameters $(\alpha,a_{\mathrm{d}},K)$. For all settings, the channel center frequencies where geometrically spaced in $\frac{\xi_{s}}{20}\cdot[2^{-6},2^{3.3}]$. We considered the following tuples of parameters: $(30,5,100)$, $(300,12,240)$, and $(3000,20,400)$. Here, the ratio $K/a_{\mathrm{d}}=20$ was fixed in all cases, but $a_{\mathrm{d}}$ and $K$ were adjusted to accommodate for bandwidth variations with changing $\alpha$.

The dimensionality of the considered audio data renders a systematic comparison to existing implementations of some established methods, e.g., [19, 16] unfeasible, such that we resort to fast Griffin-Lim [20, 23] as baseline method. We compare four different methods: wavelet PGHI (WPGHI, proposed), filter bank PGHI (FBPGHI, [47]), fast Griffin-Lim with random initialization (R-FGLIM, [23]) and fast Griffin-Lim initialized with the result of WPGHI (W-FGLIM, proposed). Fast Griffin-Lim was restricted to at most $150$ iterations.555Although we are mainly interested in the reconstruction quality and not in computational performance, it is worth mentioning that the solutions of plain WPGHI and FBPGHI are computed in a small fraction of the time required for executing either R-FGLIM or W-FGLIM, even if, at the cost of reconstruction quality, the maximum number of iterations was significantly reduced. Spectral convergence of the four methods on all test signals is shown in Figure 1 for the different parameter sets $(\alpha,a_{\mathrm{d}},K)$. The means and standard deviation across all signals, for every method and parameters set are shown in Table I.

It can be seen that on average, plain WPGHI (proposed) outperforms both R-FGLIM and FBPGHI on all parameter sets, although FBPGHI approaches the other methods for larger values of $\alpha$. Moreover, W-FGLIM (proposed) shows significant improvements over either WPGHI or R-FGLIM. Looking at the individual signals more closely, we see in Figure 1 that there are only very few cases in which R-FGLIM yields a better result than W-FGLIM. For $\alpha=3000$, all methods show comparable performance, with the exception of W-FGLIM, which still provides a clear advantage.

The figures and the computed standard deviations both suggest that methods that perform well on average are prone to larger performance fluctuation between individual signals. However, there is a small set of signals on which all methods perform badly, indicating that the fault is with the wavelet representation rather than the method applied for phaseless reconstruction. Notably, signal $65$ (corresponding to the rightmost signal in Figures 1–2) from the SQAM database, which, in the considered range, only contains a sustained, extremely low-pitched note is badly resolved by the employed WT and yields the worst spectral convergence of all signals, for any of the employed methods.

Informal listening showed that, for $\alpha=3000$ and an untrained listener, the obtained reconstructions are mostly indistinguishable from the original signal. For $\alpha=300$ and, more prominently, $\alpha=30$, FBPGHI often introduces a characteristic pitch-shift, likely due to a wrongly estimated time-direction phase derivative, while R-FGLIM suffers from undesired frequency modulation artifacts; both types of distortion are most easily audible in simple signals, such as signal $1$ (sine wave) and $4$ (electronic gong) of the SQAM database and not present in the reconstructions provided by WPGHI. For $\alpha=30$ and some select cases, e.g., signals $16$ (clarinet) and $32$ (triangle), audible distortions were present in the solutions by W-FGLIM, but not those by plain WPGHI, indicating that the observed improvement in terms of spectral convergence does not necessarily provide a perceptual improvement. To confirm our observations and to form their own opinion, the reader is invited to visit the manuscript webpage http://ltfat.github.io/notes/053/.

As a side note, the audio examples666Avalabile at http://ltfat.github.io/notes/040/ provided with [45] have clearly audible artifacts for signal $54$ (male German speech), for STFT-based PGHI and some competing algorithms. These artifacts are not present in any of the reconstructions we obtained using WPGHI, R-FGLIM, or W-FGLIM, for any considered parameter set, indicating that in some cases, usage of the WT may provide a genuine advantage over the STFT.

VII-B Experiment II—Changing the Redundancy

In a second set of experiments, we investigate the influence of the redundancy $K/a_{\mathrm{d}}$ on the performance of the proposed methods WPGHI and W-FGLIM. For this purpose, we fixed an intermediate value for the order parameter, setting $\alpha=1000$ and consider redundancies $K/a_{\mathrm{d}}\in\{3,5,10,30\}$. Here, however, the reconstruction quality does not only depend on the accuracy of WPGHI on the given magnitude coefficients, but also on the robustness of the synthesis by the dual system. This robustness can be quantified by the so-called frame bound ratio of the respective wavelet system (for details see [71, 72, 53]).

In the ranges considered and for fixed $K/a_{\mathrm{d}}$, the number $K$ of frequency channels has a larger influence on WPGHI performance than the decimation step $a_{\mathrm{d}}$, which was generally small. On the other hand, the wavelet frame bound ratio deteriorates very quickly777Large frame bound ratios also decrease numerical stability, such that audio file generation from the obtained reconstructions is prone to clipping artifacts. for too large decimation steps $a_{\mathrm{d}}$. Hence, the choice of wavelet parameters was a trade-off between the two factors with no clear optimal solution. After some preliminary testing, we fixed the following parameter sets $(\alpha,a_{\mathrm{d}},K)$: Low redundancy $(1000,30,90)$ (low), Medium redundancy $(1000,25,125)$ (medium), Medium high redundancy $(1000,18,180)$ (medhigh), High redundancy $(1000,10,300)$ (high).

Similar to Experiment I, mean value and standard deviation over all signals are presented in Table II, for all parameter sets, with detailed results for all test signals shown in Figure 2. Additionally, Figure 3 shows an example of the difference between the target phase and the WPGHI-proposed phase estimate at different redundancies.

As expected, performance of both proposed methods increases with increasing redundancy. These improvements are apparent in the average performance over the whole signal set, but in most cases also on the level of individual signals, as can be seen in Figure 2. Similar to Experiment I, the average improvement in terms of spectral convergence of W-FGLIM over WPGHI is significant, ranging between $10$ (low) and $6$ dB (high).

The performance difference between the different redundancies is also apparent in Figure 3. The phase reconstruction quality can be visually estimated from the characteristics of the difference between the target phase and the proposed estimate. Large areas of flat color indicate good quality; as the quality decreases, the phase difference becomes more patchy, with stronger fluctuation within patches. The figure shows that this patchiness is closely linked to redundancy of the underlying wavelet representation, or more generally, the employed sampling scheme. This behavior is characteristic for WPGHI and was previously observed for STFT-based PGHI [45] as well.

For all redundancies, informal listening has shown good to excellent reconstruction quality. Usually, the reconstructions using medium, medium high, or high redundancy and either of the proposed methods were indistinguishable from one another. In the case where distortions are audible at all, they were never found to be irritating. At low redundancy, the results are generally good as well, but a larger number of examples has clearly audible artifacts, e.g., signals $14$ (oboe), $15$ (cor anglais) and $16$ (clarinet).

VIII Conclusion

We have highlighted a number of interesting properties of the partial derivatives of the continuous WT. In particular, we showed how the log-magnitude and phase partial derivatives are related to wavelet transforms with modified mother wavelets. We characterize the class of wavelets that generate analytic WTs and obtain as a result a generalization of the Cauchy wavelets. Based on the analyticity of these WTs, we obtain explicit phase-magnitude relations. Similar to the Gaussian STFT, the phase-magnitude relationship can be used as a basis for implementing magnitude-only reassignment and phaseless reconstruction. We explored the second application, providing very good results when applied to complex audio data, in most cases free of any perceptible distortion. We demonstrated that reconstruction quality is often notably improved over the established Griffin-Lim algorithm and our own previous implementation, relying on approximate phase-magnitude relations for Gaussian filter banks.

Future work will be concerned with real-time implementation of the proposed algorithm in the style of RTPGHI [66] and its generalization to nonuniform decimation schemes. The new phase-magnitude relations for WTs are crucial for the derivation of more appropriate phase approximation schemes for more general time-frequency filters, improving and extending the previously proposed filter bank PGHI algorithm [47]. We also plan to investigate phase-magnitude relations for the polyanalytic generalizations of the Cauchy WT [73]. Finally, the proposed scheme can serve as starting point for a wavelet-based phase vocoder for time-stretching and pitch-shifting of audio in the spirit of [74].

Acknowledgments

We wish to thank Andrés Marafioti for performing a small listening test confirming the reported observations on the perceptual quality of the reconstructed audio samples. Furthermore, we would like to thank Patrick Flandrin for pointing us to the literature on Klauder wavelets which coincide with the analyticity inducing wavelets given by (6). Finally, we thank the reviewers, whose comments helped us improve our results and their presentation.

Appendix A Proof of Theorem 1

Proof.

Under the assumption that $\psi\in\mathbf{L}^{2}(\mathbb{R})$ is continuously differentiable with $\psi^{\prime},\mathbf{T}\psi^{\prime}\in\mathbf{L}^{2}(\mathbb{R})$, we can exchange differentiation and integration in (2), see Appendix D. Thus, the partial derivatives of $W_{\psi}s$ can be expressed as

[TABLE]

and

[TABLE]

Here, we used that $(\mathbf{T}\psi)^{\prime}=\psi+\mathbf{T}\psi^{\prime}$ and that that the WT is conjugate linear with respect to the chosen wavelet. Using that $\frac{\partial}{\partial x}\log(W_{\psi}s)=\frac{\frac{\partial}{\partial x}W_{\psi}s}{W_{\psi}s}$ and similarly for the partial derivative with respect to $y$, we obtain for all $(x,y)$ with $W_{\psi}s(x,y)\neq 0$,

[TABLE]

Taking real and imaginary parts in (44) results in (3) and (4), respectively. ∎

Appendix B Proof of Theorem 2

We first argue that analyticity of $h$ and differentiability of $f$ already imply that $\psi$ satisfies the assumptions in Theorem 1, i.e., $\psi$ is continuously differentiable with $\psi,\psi^{\prime},\mathbf{T}\psi^{\prime}\in\mathbf{L}^{2}(\mathbb{R})$. To this end, we note that the assumptions imply that $W_{\psi}s(x-aby,by)$ must be $\mathcal{C}^{\infty}$ in $x$ and $y$ for an arbitrary $s\in\mathbf{L}^{2}(\mathbb{R})$. The same holds true for $W_{\psi}s(x,y)$.

At $x=0$ and $y=1$, the derivative $\frac{\partial}{\partial x}W_{\psi}s$ can be written as

[TABLE]

which converges for every fixed $s\in\mathbf{L}^{2}(\mathbb{R})$. Thus, by a variant of the Banach-Steinhaus theorem [75, Ch. II.1, Corollary 2], the limit $\lim_{x_{0}\rightarrow 0}\frac{\psi-\mathbf{T}_{-x_{0}}\psi}{x_{0}}$ represents a continuous linear functional and hence, by Riesz representation theorem, an element $\psi_{x}$ in $\mathbf{L}^{2}(\mathbb{R})$. Rewriting the derivative for compactly supported, smooth $s$ alternatively as

[TABLE]

we see that $\psi_{x}$ is the weak derivative of $\psi$, i.e., the weak derivative of $\psi$ exists and belongs to $\mathbf{L}^{2}(\mathbb{R})$. Repeating the argument for higher derivatives guarantees that weak derivatives of arbitrary order exist and by standard Sobolev embeddings so do continuous derivatives.

Similarly, the derivative $\frac{\partial}{\partial y}W_{\psi}s$ at $x=0$ and $y=1$ can be written as

[TABLE]

Again we have the convergence $\lim_{y_{0}\rightarrow 0}\frac{\psi-\mathbf{D}_{-y_{0}}\psi}{y_{0}}=\psi_{y}$. Rewriting the derivative for $s\in\mathcal{C}^{\infty}_{00}$ alternatively as

[TABLE]

Now the operator $s\mapsto\mathbf{T}s^{\prime}$ is well defined for compactly supported, smooth $s$ with adjoint $s\mapsto-(\mathbf{T}s)^{\prime}$. Furthermore, $\psi$ is in the domain of this operator because $\langle-(\mathbf{T}s)^{\prime},\psi\rangle=\langle s,\psi_{y}\rangle$ for all $s$ in a dense subset. Thus, $\mathbf{T}\psi^{\prime}=\psi_{y}\in\mathbf{L}^{2}(\mathbb{R})$. Hence, we established all assumptions of Theorem 1.

Analyticity of the function $h$ is equivalent to it satisfying the CR equations that can be compactly expressed as $\frac{\partial}{\partial x}h=-i\frac{\partial}{\partial y}h$. To rewrite the CR equations for the function $h$ in (5), we use (42), to obtain

[TABLE]

Similarly, by (A), we have

[TABLE]

Inserting these expressions into the CR equations results in

[TABLE]

We note that this condition depends on $f$ only via the function $g(x,y)=\frac{y\frac{\partial}{\partial x}f(x,y)+iy\frac{\partial}{\partial y}f(x,y)}{f(x,y)}+\frac{i}{2}$. Moreover, using the definition of the WT in (51) implies that

[TABLE]

for all $s\in\mathbf{L}^{2}(\mathbb{R})$ and thus

[TABLE]

as a function of $t$ in $\mathbf{L}^{2}(\mathbb{R})$. In particular, this implies that $g(x,y)$ must be a constant $w\in\mathbb{C}$ and further

[TABLE]

To solve this differential equation, it is more convenient to consider the Fourier transformed equivalent of (54). Using the standard properties of the Fourier transform $\widehat{\psi^{\prime}}=2\pi i\mathbf{T}\widehat{\psi}$ and $\widehat{\mathbf{T}\psi}=-(2\pi i)^{-1}(\widehat{\psi})^{\prime}$, this is easily seen to be given by

[TABLE]

We first note that our assumption $\widehat{\psi}(\xi)=0$ for $\xi<0$ satisfies this differential equation on this domain. For $\xi>0$ and $\widehat{\psi}\neq 0$, we can easily solve the differential equation by rewriting

[TABLE]

which gives

[TABLE]

for an arbitrary constant $c\in\mathbb{C}$. To obtain the parameters used in the theorem, we substitute $-\operatorname{Im}(w)=\frac{\alpha-1}{2}$, $\frac{1}{b}=\operatorname{Re}(\gamma)$, $a=\operatorname{Im}(\gamma)$, and $\operatorname{Re}(w)=-\beta$. Based on our assumptions, we obtain the constraints $\operatorname{Re}(\gamma)>0$ and $\alpha>-1$ to guarantee $\widehat{\psi}\in\mathbf{L}^{2}(\mathbb{R})$.

Appendix C Proof of Theorem 3

We can use the CR equations to obtain relationships between the derivatives of real and imaginary parts of the WT, which we will show to result in (11) and (12). For an arbitrary analytic function $h=u+iv$ the CR equations hold and are given by $\frac{\partial}{\partial x}u=\frac{\partial}{\partial y}v$ and $\frac{\partial}{\partial y}u=-\frac{\partial}{\partial x}v$. Writing $h=Me^{i\phi}$, the CR equations imply that

[TABLE]

and

[TABLE]

Using (58) and (59) for the function given by (7), yields

[TABLE]

and

[TABLE]

These are equivalent to

[TABLE]

and

[TABLE]

Inserting (62) into (63) finally results in

[TABLE]

which is equivalent to

[TABLE]

This concludes the proof.

Appendix D On Differentiability of the Wavelet Transform

Denote by $\Phi_{z}$ and $\Psi_{z}$ the difference operators

[TABLE]

and note that

[TABLE]

as well as

[TABLE]

provided that the right-hand sides converge.

We first show that, for $z\rightarrow 0$

[TABLE]

as functions in $\mathbf{L}^{2}(\mathbb{R})$, provided that $\psi$, $\psi^{\prime}$, and $\mathbf{T}\psi^{\prime}$ are elements of $\mathbf{L}^{2}(\mathbb{R})$.

Since $\psi^{\prime}$ is assumed to be in $\mathbf{L}^{2}(\mathbb{R})$, there is for every $\tilde{\epsilon}>0$ an $r_{\tilde{\epsilon}}>0$, such that \big{\lVert}\psi^{\prime}\chi_{\mathbb{R}\setminus\overline{B_{r_{\tilde{\epsilon}}}(0)}}\big{\rVert}<\tilde{\epsilon}. Moreover, by the fundamental theorem of calculus,

[TABLE]

Using (67), we obtain

[TABLE]

where we used Jensen’s inequality and Fubini’s theorem. Furthermore, we have $|\psi^{\prime}(t-s)-\psi^{\prime}(t)|<\epsilon_{r}(|s|)$ for all $|t|<r$ and $|s|<1$ by uniform continuity of $\psi^{\prime}$ on the compact set $B_{r+1}(0)$, where $\epsilon_{r}(\delta)\searrow 0$ for $\delta\rightarrow 0$. Thus, similar to (68), we obtain

[TABLE]

for $|z|<1$. Finally, for every $\tilde{\epsilon}$, there is a $z_{\tilde{\epsilon}}\in(0,1)$, such that $\epsilon_{r_{\tilde{\epsilon}}+1}(z_{\tilde{\epsilon}})<\tilde{\epsilon}/\sqrt{2(r_{\tilde{\epsilon}}+1)}$, which implies $\|\Phi_{z}\psi+\psi^{\prime}\|<\sqrt{7}\tilde{\epsilon}$ for all $|z|\leq z_{\tilde{\epsilon}}$.

A similar, but slightly more complicated argument shows that $\Psi_{1+z}\psi\rightarrow-(\psi/2+\mathbf{T}\psi^{\prime})$, provided $\psi,\mathbf{T}\psi^{\prime}\in\mathbf{L}^{2}(\mathbb{R})$. Here, we start with an $r_{\tilde{\epsilon}}>0$, such that $\|\mathbf{T}\psi^{\prime}\chi_{\mathbb{R}\setminus\overline{B_{r_{\tilde{\epsilon}}}(0)}}\|<\tilde{\epsilon}$ and $\|\psi\chi_{\mathbb{R}\setminus\overline{B_{r_{\tilde{\epsilon}}}(0)}}\|<\tilde{\epsilon}$. As above, we use the fundamental theorem of calculus to obtain

[TABLE]

Similar to (68), this results in

[TABLE]

Furthermore, we can bound for $t\in\overline{B_{2r}(0)}$

[TABLE]

for any $\lvert s\rvert\leq 1/2$, where $\epsilon_{r}(\delta)\to 0$ monotonically for $\delta\to 0$. Analogously, we obtain

[TABLE]

Equations (70)–(72) imply

[TABLE]

and thus

[TABLE]

for $\lvert z\rvert<1/(4r_{\tilde{\epsilon}})$ and where we assumed for simplicity $2r_{\tilde{\epsilon}}>1$. Again, choosing $\lvert z\rvert$ sufficiently small implies the proposed convergence.

Thus, we finished the proof of (66), which implies

[TABLE]

and

[TABLE]

Because convergence is in $\mathbf{L}^{2}(\mathbb{R})$, we can exchange limit and integral by continuity of the inner product.

Clearly, this argument can be repeated to obtain higher order derivatives, provided that $\psi$ has sufficient regularity and decay.