Solving multiple windowed STFT phase retrieval problems in phase and amplitude respectively

ABSTRACT

We study the Phase Retrieval (PR) problem under the phaseless short-time Fourier transform (STFT) measurements. This paper proposes a novel algorithm named PAR to solve the STFT PR problem in phase and amplitude respectively with a milder retrieval condition compared with the original methods. First, a symmetric undirected graph of signals is proposed for the computation of the relative phase. Then the retrieval conditions of STFT PR problem are discussed for a single window case and some weaker retrieval conditions are proposed compared with the LS method. We also discuss STFT PR problem in multiple windows and establish retrieval theorems without restrictions of sliding step-size L. We give some numerical results of the PAR algorithm.

KEYWORDS:

• Phase retrieval
• short-time Fourier transform
• phase propagation
• relative phase-multiple windows

MSC:

• 94A12

1. Introduction and problem equation

Phase Retrieval(PR) is a classical problem which considers recovering a signal from the magnitude of the measurements. It can be applied in X-ray crystallography [1,2], diffraction imaging system [3], phase measurement in astronomy [4,5] and optics [6,7]. When the phase retrieval measuring vectors are chosen as Short Time Fourier Transform(STFT) measurements [6], it is an STFT PR problem. This problem serves as the model for ptychography, in which one or multiple moving probes are used to sense multiple diffraction measurements [8,9] and ultra-short laser pulse measurement techniques [10,11]. Different number of probes and laser pulses corresponds single or multiple STFT PR problem [12,13].

This paper considers the STFT phase retrieval, which recovering a signal from its STFT magnitude. Set $\mathit{x}=\left(x\right(0),x(1),\dots ,x(N-1){)}^T\in {\mathbb{C}}^N.$ The STFT can be interpreted as the Fourier transform of an N-dimensional signal $\mathit{x}\in {\mathbb{C}}^N$ multiplied by series of sliding windows ${\mathit{g}}_r\in {\mathbb{R}}^N$ with N-periodic extension: (1) $X_r(m,k)=\sum _{n=0}^{N-1}x\left(n\right)g_r(Lm-n){\mathrm{e}}^{-i2\pi kn/N},$(1) with $0\le m\le \lceil N/L\rceil -1$, $0\le k\le N-1$ and $1\le r\le R$. m represents the sliding order of the window ${\mathit{g}}_r$ and L represents the sliding step-size. STFT PR considers the problem of recovering $\mathit{x}$ from $|X_r(m,k)|$. When R = 1, the problem is a single STFT phase retrieval problem, while R>1 is a multiple STFT phase retrieval problem [14,15].

To solve STFT PR, not only the algorithm but also retrieval conditions need to be considered. It means the conditions under which the STFT magnitude uniquely identifies signals (up to a global phase) [16] under the algorithm. Researchers find that a suitable restriction of measurements including windows' number, type, length and step size enforce retrieval uniqueness.

For single STFT PR problem, different algorithms have various restrictions of windows. Denote the length nonzero window $\mathit{g}$ by W. When the measurement is with adjacent windows (L = 1), SDP approach [17] has been extended to STFT PR [18] and has uniqueness guarantees with computational complexity $\mathcal{O}\left(N^3\right)$. When the sliding step-size L = 1 and $W\ge \lceil \frac{N+1}{2}\rceil$, a Least Squares (LS) approach [19] obtains the stable solution of STFT PR. In addition, a robust method proposed in [20] uses phase polarization to solve STFT PR under a sliding window with length W = N, which is too restrictive in many applications.

As for multiple PR,  [14] discusses some sufficient and necessary conditions of multi-windows STFT PR and an algorithm based on the graph was proposed to recover the signal. But this retrieval condition is quite special ($W\le \lceil \frac{N}{2}\rceil$) and can only be satisfied in strict conditions. Greedy Angular Synchronization [21–23] is proposed to compute relative phase by autocorrelation. While it assumes that all the autocorrelation. $x\left(n_1\right)\overline{x\left(n_2\right)}$ can be obtained for all $|n_1-n_2|<\delta$ where δ is a given length, it may not easy to satisfy this condition in STFT PR problem.

No matter for single or multiple windows STFT PR, those algorithms have special retrieval restrictions of measurements, which restrict the applications. A better algorithm should consider the retrieval conditions and be suitable for broader applications. In this paper, an algorithm named PAR, which compute phase and amplitude respectively (PAR), is proposed with a milder retrieval restrictions of measurements for both single and multiple STFT PR. First, a graph is established based on STFT measurements and the retrieval conditions of PAR are transformed to the graph connectivity. Based on the analysis of the graph connectivity, some milder retrieval conditions are proposed for different types of windows, which has a broader field of application. Error estimation and experiments are also given to verify the effectiveness of the PAR algorithm.

Throughout the paper, we use the following notation. Bold-face small and capital letters denote vectors and matrices respectively. $gcd(a,b,c)$ denotes the greatest common factor between a, b and c. $\mathcal{F}$ denotes 1D discrete Fourier transform (DFT) and ${\mathcal{F}}^{-1}$ denotes 1D discrete inverse Fourier transform (IDFT). $\mathcal{F}\mathit{x}\in {\mathbb{C}}^N$ is defined as (2) $\mathcal{F}\mathit{x}\left(k\right)=\sum _{n=0}^{N-1}x\left(n\right){\mathrm{e}}^{-i2\pi kn/N},0\le k\le N-1,$(2) and ${\mathcal{F}}^{-1}\mathit{x}\in {\mathbb{C}}^N$ is defined as (3) ${\mathcal{F}}^{-1}\mathit{x}\left(n\right)=\sum _{k=0}^{N-1}\frac{1}{N}x\left(k\right){\mathrm{e}}^{i2\pi kn/N},0\le n\le N-1.$(3) The paper is organized as follows. In Section 1, we introduce the STFT PR problem and related works. Section 2 discusses the PR problem in the single-window and multi-window respectively. Section 3 gives the reconstruction algorithm and error estimation. Section 4 compares the numerical results of PAR and LS method. Section 5 is the summary of this paper.

2. Analysis of phase retrieval conditions

Define $|{\mathit{X}}_r^{\left(m\right)}{|}^2=\left(|X_r\right(m,0){|}^2,\dots ,|X_r(m,N-1\left){|}^2\right)$. Instead of recovering x from $|{\mathit{X}}_r^{\left(m\right)}{|}^2$, this paper considers the acquired data ${\mathit{Y}}_r$ by taking IDFT of $|{\mathit{X}}_r^{\left(m\right)}{|}^2$ for $0\le m\le \lceil N/L\rceil -1$, which is defined as (4) $Y_r(m,l)=\frac{1}{N}\sum _{k=0}^{N-1}|X_r(m,k){|}^2{\mathrm{e}}^{i2\pi kl/N}=\sum _{n=0}^{N-1}x\left(n\right)\overline{x(n+l)}{g}_r(mL-n)\overline{g_r(mL-n-l)},$(4) with $l=0,1,\dots ,N-1$, $0\le m\le \lceil N/L\rceil -1$ [19]. Therefore, recovering $\mathit{x}$ from ${\mathit{X}}_r$ is equivalent to recover it from ${\mathit{Y}}_{\mathit{r}}$. To recover the signal $\mathit{x}$ uniquely, the injectivity property of $\mathit{x}⟶\mathit{Y}$ should be guaranteed. Actually, there are classes of signals that have the same measurements after injection. For arbitrary $\theta \in \mathbb{C}$, $\mathit{x}$ and $\mathit{x}{\mathrm{e}}^{i\theta }$ produce the same intensity measurements. Thus we consider $\mathit{x}$ and $\mathit{x}{\mathrm{e}}^{i\theta }$ as in the same class $\left[\mathit{x}\right]\in {\mathbb{C}}^M/\sim$, where ∼ is the equivalence relation of being identical up to a global phase factor. Then the distance between two vectors in ${\mathbb{C}}^N$ is defined as (5) $d(\mathit{z},\mathit{x})=d(\mathbf{\left[}\mathit{z}\mathbf{\right]},\mathbf{\left[}\mathit{x}\mathbf{\right]})=\underset{\varphi \in [0,2\pi )}{min}\Vert \mathit{z}-\mathit{x}{\mathrm{e}}^{i\varphi }{\Vert }_2.$(5) When $d(\mathit{z},\mathit{x})=0$, $\mathit{x}$ and $\mathit{z}$ are equivalent up to a global phase shift.

2.1. Single adjacent window case (R = 1, L = 1)

To find the relation between success of recovery and windows' parameters, we first discuss the STFT PR problem with a single window $\mathit{g}$. Denote $\tilde{x}(n,l)=x\left(n\mathrm{mod}N\right)\overline{x\left(\right(n+l\left)\mathrm{mod}N\right)}$ and $\tilde{g}(n,l)=g\left(n\mathrm{mod}N\right)\overline{g\left(\right(n-l\left)\mathrm{mod}N\right)}$. Then Equation (4(4) $Y_r(m,l)=\frac{1}{N}\sum _{k=0}^{N-1}|X_r(m,k){|}^2{\mathrm{e}}^{i2\pi kl/N}=\sum _{n=0}^{N-1}x\left(n\right)\overline{x(n+l)}{g}_r(mL-n)\overline{g_r(mL-n-l)},$(4) ) can be represented by (6) $Y(m,l)=\sum _{n=0}^{N-1}\tilde{x}(n,l)\tilde{g}(mL-n,l).$(6) Define ${\mathit{y}}^{\left(l\right)}=\left\{Y\right(m,l){\}}_{m=0}^{\lceil N/L\rceil -1}$, ${\mathit{x}}^{\left(l\right)}=\left\{\tilde{x}\right(n,l){\}}_{n=0}^{N-1}$, and matrix ${\mathit{G}}^{\left(l\right)}\in {\mathbb{R}}^{\lceil N/L\rceil \times N}$, where the $(m,n)th$ entry of ${\mathit{G}}^{\left(l\right)}$ is equal to $\tilde{g}(mL-n,l)$. Then Equation (6(6) $Y(m,l)=\sum _{n=0}^{N-1}\tilde{x}(n,l)\tilde{g}(mL-n,l).$(6) ) can be written as (7) ${\mathit{y}}^{\left(l\right)}={\mathit{G}}^{\left(l\right)}{\mathit{x}}^{\left(l\right)}.$(7) Due to L = 1, ${\mathit{G}}^{\left(l\right)}$ is a circulant square matrix which can be diagonalized by ${\mathit{G}}^{\left(l\right)}={\mathit{F}}^{\mathbf{†}}{\mathbf{\Sigma }}^{\left(l\right)}\mathit{F}$, where ${\mathbf{\Sigma }}^{\left(l\right)}$ is a diagonal matrix whose entries are given by the DFT of the first column of ${\mathit{G}}^{\left(l\right)}$. $\mathit{F}$ is the DFT matrix and ${\mathit{F}}^{\mathbf{†}}$ represents the conjugate transpose matrix of $\mathit{F}$ [24]. If the matrices ${\mathit{G}}^{\left(l\right)}$ are invertible for several $0\le l\le N-1$, then we can compute ${\mathit{x}}^{\left(l\right)}$ by (8) ${\mathit{x}}^{\left(l\right)}=({\mathit{G}}^{\left(l\right)}{)}^{-1}{\mathit{y}}^{\left(l\right)}.$(8) To analyse the invertibility of ${\mathit{G}}^{\left(l\right)}$, the definition of ‘non-vanishing’ is given as follows:

Definition 2.1

$\mathit{x}\in {\mathbb{C}}^N$ is non-vanishing means $x\left(k\right)\ne 0$ for all $k=0,1,\dots ,N-1$.

Denote ${\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l\right)}=\left(\tilde{g}\right(0,l),\tilde{g}(1,l),\dots ,\tilde{g}(N-1,l\left)\right)$. ${\mathit{G}}^{\left(l\right)}$ is invertible if and only if $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l\right)}$ is non-vanishing since ${\mathit{G}}^{\left(l\right)}={\mathit{F}}^{\mathbf{†}}{\mathbf{\Sigma }}_{\mathit{l}}\mathit{F}$. Based on this, LS algorithm proposed in [19] can be used to solve the PR problem when $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l\right)}$ is non-vanishing for $l=0,1,\dots ,N-1$. Nevertheless, Theorem 3.3 in [19] shows that if $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l\right)}$ is non-vanishing only for l = 0, 1, signal can still be recovered. It inspires us to consider recovering the signal under milder and wider conditions.

Suppose signal $\mathit{x}$ is non-vanishing. In fact, for any $n_1$ and $n_2$, since $\mathit{x}$ and $\mathit{x}{\mathrm{e}}^{i\theta }$ are in the same class of $\left[\mathit{x}\right]\in {\mathbb{C}}^M/\sim$, we do not need to compute all the phase, but the relative phase [6]: (9) ${\rho }_{n_1{n}_2}=\frac{\left(x\right(n_1\left)/|x\right(n_1\left)|\right)}{\left(x\right(n_2\left)/|x\right(n_2\left)|\right)}=\frac{x\left(n_1\right)\overline{x\left(n_2\right)}}{|x\left(n_1\right)\overline{x\left(n_2\right)}|}.$(9) To compute the signal phase, we define an undirected graph for STFT phase retrieval: (10) $G(\mathit{x},\mathit{g},L):=\left(V\right(\mathit{x}),E(L\left)\right),$(10) with a set of nodes: $V\left(\mathit{x}\right):=\{0\le n\le N-1|x(n)\ne 0\}$ and a set of edges: $E\left(L\right):=\left\{\right(n,n^{\prime })\in V(\mathit{x})\times V(\mathit{x}\left)|{\rho }_{n{n}^{\prime }}\text{ can be computed directly by (8) and (9)}\right\}.$ Relative phase information can be transferred through the relationship between the phase of $n_1,n_2,n_3$: (11) ${\rho }_{n_1{n}_3}={\rho }_{n_1{n}_2}{\rho }_{n_2{n}_3}.$(11) If we assign an arbitrary nonzero vertex $x\left(n_0\right)$ to have nonzero phase ${\rho }_0=x\left(n_0\right)/|x\left(n_0\right)|$, then all the phase can be computed due to phase propagation equation (11(11) ${\rho }_{n_1{n}_3}={\rho }_{n_1{n}_2}{\rho }_{n_2{n}_3}.$(11) ) if graph G is connected. In this paper, we can use the same strategy as [23] to choose the initial phase ${\rho }_0$. Therefore, the graph connectivity decides the retrieval uniqueness of phase.

Lemma 2.1

If there exists non-vanishing $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$, then ${\rho }_{n,n+l_0}$ can be computed by the phase of ${\mathit{x}}^{\left(l\right)}$ for any $n\in V\left(\mathit{x}\right)$.

Proof.

According to Equations (8(8) ${\mathit{x}}^{\left(l\right)}=({\mathit{G}}^{\left(l\right)}{)}^{-1}{\mathit{y}}^{\left(l\right)}.$(8) ) and (9(9) ${\rho }_{n_1{n}_2}=\frac{\left(x\right(n_1\left)/|x\right(n_1\left)|\right)}{\left(x\right(n_2\left)/|x\right(n_2\left)|\right)}=\frac{x\left(n_1\right)\overline{x\left(n_2\right)}}{|x\left(n_1\right)\overline{x\left(n_2\right)}|}.$(9) ), ${\rho }_{n,n+l_0}$ can be computed when ${\mathit{G}}^{\left(l_0\right)}$ is invertible, which is equivalent to $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$ is non-vanishing.

Since the STFT of signal is periodic, the connectivity of graph G is related to l and N. We have the following results:

Lemma 2.2

Given a graph G with N vertices, suppose there exists k such that any two vertices $n_1$ and $n_2$ are connected if and only if $(n_1-n_2)\mathrm{mod}N\equiv k$, then graph G is connected if and only if $gcd(k,N)=1$.

Proof.

Graph connectivity can be judged by the adjacency matrix $A\in \{0,1{\}}^{N\times N}$. If all the elements in $S=I+A+A^2+\cdots +A^N$ are non-zero, graph is connected. The adjacency matrix of G is in this form (12) $\begin{array}{r}A=12\cdots \cdots k+1\cdots \text{N}00\cdots \cdots 10\cdots 000\ddots \cdots 01\ddots 0⋮\ddots \cdots \ddots \ddots \ddots \ddots ⋮00\ddots \cdots \ddots \ddots \ddots 110\ddots \cdots \ddots \ddots \ddots 001\ddots \cdots \ddots \ddots \ddots 0⋮\ddots \cdots \ddots \ddots \ddots \ddots 00\cdots \cdots 10\cdots \cdots 0,\end{array}$(12) (13) $\begin{array}{r}{A}^{l+1}=12\cdots \cdots k+1+l\cdot gcd(k,N)\cdots \text{N}00\cdots \cdots 10\cdots 000\ddots \cdots 01\ddots 0⋮\ddots \cdots \ddots \ddots \ddots \ddots ⋮00\ddots \cdots \ddots \ddots \ddots 110\ddots \cdots \ddots \ddots \ddots 001\ddots \cdots \ddots \ddots \ddots 0⋮\ddots \cdots \ddots \ddots \ddots \ddots 00\cdots \cdots 10\cdots \cdots 0.\end{array}$(13) It shows that $diag(A,k)=\mathbf{1}\in {\mathbb{R}}^N$ and $diag(A^l,k+l\cdot gcd(k,N\left)\right)=\mathbf{1}\in {\mathbb{R}}^N$ for $0\le l\le N$. Since $S=I+A+A^2+\cdots +A^N$, then all elements in S are non-zero is equivalent to that $diag(A^l,k+l\cdot gcd(k,N\left)\right)$ for $0\le l\le N$ covers all the position in matrix A, which means that $gcd(k,N)=1.$

A generalized result of Lemma 2.2 is given by Lemma 2.3.

Lemma 2.3

Given a graph G with N vertices, suppose any two vertices $n_1$ and $n_2$ are connected if $(n_1-n_2)\mathrm{mod}N=k_i$ for $i=1,2,\dots ,m$. G is connected if and only if $gcd(k_1,k_2,\dots ,k_m,N)=1$.

Proof.

The lemma can be inferred by S defined in Lemma 2.2. Denote the adjacency matrix is A, then $diag(A,k_i)=\mathbf{1}\in {\mathbb{R}}^N$ and $diag(A^l,k_i+l\cdot gcd(k_i,N\left)\right)=\mathbf{1}\in {\mathbb{R}}^N$ for $0\le l\le N$ and $1\le i\le m$. Only when $gcd(k_1,k_2,\dots ,k_m,N)=1$, all elements in S are non-zero.

Lemmas 2.2 and 2.3 give the relation between k and N. Combining with Lemma 2.1, we get two theorems as follows:

Theorem 2.4

Assume there exists $l_0$, such that $gcd(l_0,N)=1$ and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$ is non-vanishing, then all phases can be recovered.

Theorem 2.5

Assume there exists $l_1,l_2,\dots ,l_m$ such that $gcd(l_1,l_2,\dots ,l_m,N)=1$, and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_i\right)}$ is non-vanishing for $i=1,2,\dots ,m$, then the phase can be recovered.

Proof.

Since $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$ is non-vanishing, ${\mathit{x}}^{\left(l_0\right)}=({\mathit{G}}^{\left(l_0\right)}{)}^{-1}{\mathit{y}}^{\left(l_0\right)}$ can be computed, which means the relative phase between $x\left(n\right)$ and $x\left(\right(n-l_0\left)\mathrm{mod}N\right)$ can be obtained. According to Lemma 2.2 and 2.3, all phases can be computed.

Theorem 2.4 and 2.5 give two phase recovery conditions. Then we turn to consider how to compute the amplitude. As is shown in Theorem 3.1 in [19] and Algorithm 1 in [23], the amplitude can be computed when $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(0\right)}$ is non-vanishing. However, if $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(0\right)}$ is vanishing, could signal $\mathit{x}$ still be recovered? Here we make an extension of proposition III.4 in [19].

Theorem 2.6

If there exists $l_0,l_1$ such that $l_1-l_0=M$ where $M=gcd(l_1,l_0)$, and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$, $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_1\right)}$ are both non-vanishing, then the amplitude of the signal $\mathit{x}$ can be recovered by (14) $|x\left(l_0\right){|}^2=\frac{\prod _{k=0}^{l_0/M}|x\left(Mk\right)\overline{x(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x\left(Mk\right)\overline{x(Mk+l_1)}|}.$(14)

Proof.

Without loss of generality, since $l_1=l_0+M$, then (15) $\begin{array}{rl}\frac{\prod _{k=0}^{l_0/M}|x\left(Mk\right)\overline{x(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x\left(Mk\right)\overline{x(Mk+l_1)}|}& =\frac{|x\left(0\right)\overline{x\left(l_0\right)}x\left(M\right)\overline{x(M+l_0)}\dots x\left(l_0\right)\overline{x\left(2l_0\right)}|}{x\left(0\right)\overline{x\left(l_1\right)}x\left(M\right)\overline{x(M+l_1)}\dots x(l_0-M)\overline{x(l_0-M+l_1)}|}\\ & =\frac{|x\left(0\right)\overline{x\left(l_0\right)}x\left(M\right)\overline{x(M+l_0)}\dots x\left(l_0\right)\overline{x\left(2l_0\right)}|}{|x\left(0\right)\overline{x(M+l_0)}x\left(M\right)\overline{x(M+l_1)}\dots x(l_0-M)\overline{x\left(2l_0\right)}|}\\ & =|x\left(l_0\right)\overline{x\left(l_0\right)}|=|x\left(l_0\right){|}^2.\end{array}$(15) We could compute the amplitude of $x\left(l_0\right)$ and other vertices amplitude can be computed by the same equation.

Remark 2.1

Since Equation (2.6) involves division, the signal should be non-vanishing and not too small. Here is an example. If $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(1\right)}$ and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(0\right)}$ are both non-vanishing, then (16) $|x\left(1\right){|}^2=\frac{|x\left(0\right)\overline{x\left(1\right)}||x\left(1\right)\overline{x\left(2\right)}|}{|x\left(0\right)\overline{x\left(2\right)}|}.$(16) Other vertices can also be computed by this equation.

Corollary 2.7

If there exists $l_0$ such that $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{\left(l_0\right)}$ and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}^{(l_0+1)}$ are non-vanishing, then the signal can be recovered.

Proof.

Since $gcd(l_0,l_0+1)=1$ for any integer $l_0$, then all phases can be computed according to Theorem 2.4. And it also ensures the amplitude can be computed by Theorem 2.6. Therefore, the signal can be recovered.

Remark 2.2

Corollary 2.7 gives a sufficient condition for the recovery of $\mathit{x}$. Theorem III.4 proposed in  [19] is a special case of Corollary 2.7.

Now we can prove that when $W<N/2$, our method can also recover the phase, which expands the scope of W compared with LS model.

Corollary 2.8

If $gcd(W,N)=1$ and $W<N/2$, then the signal phase can be recovered.

Proof.

Consider the matrix ${\mathit{G}}^W\in {\mathbb{R}}^N\times N$ where the $(m,n)$th entry of ${\mathit{G}}^W$ is given by $\tilde{g}(m-n,W)$. Since $W<N/2$ and $gcd(N,W)=1$, there must exist unique $n_0$ such that $\tilde{g}(n_0,W)\ne 0$. Therefore, we have (17) $Y(m,W)=\sum _{n=0}^{N-1}\tilde{x}(n,W)\tilde{g}(m-n,W)=\tilde{x}(m-n_0,W)\tilde{g}(n_0,W),m=0,1,\dots ,N-1.$(17) Since $\tilde{g}(n_0,W)\ne 0$, (18) $\tilde{x}(m-n_0,W)=\left\{\tilde{g}\right(n_0,W\left){\}}^{-1}Y\right(m,W).$(18) Then ${\rho }_{m-n_0,m-n_0+W}$ can be computed as $\begin{array}{rl}{\rho }_{m-n_0,m-n_0+W}& =\frac{\tilde{x}(m-n_0,W)}{|\tilde{x}(m-n_0,W)|}={\left\{\frac{\tilde{g}(n_0,W)}{|\tilde{g}(n_0,W)|}\right\}}^{-1}\times \frac{Y(m,W)}{|Y(m,W)|},\\ m& =0,1,\dots ,N-1.\end{array}$ Since m is arbitrary, for any $n_1$ and $n_2$, ${\rho }_{n_1{n}_2}$ can be computed if $n1-n2=W\mathrm{mod}N$. Therefore, there exists k = W so that ${\rho }_{n,n+k}$ can be computed for $0\le n\le N-1$. Together with $gcd(W,N)=1$, the signal phase can be recovered according to Lemma 2.2 where k = W.

By now, the discussion is based on L = 1. We now extend the results to L>1. In fact, since ${\mathit{G}}^{\left(l\right)}$ defined in (7(7) ${\mathit{y}}^{\left(l\right)}={\mathit{G}}^{\left(l\right)}{\mathit{x}}^{\left(l\right)}.$(7) ) is not invertible, ${\rho }_{n,n+l}$ cannot be computed by (8(8) ${\mathit{x}}^{\left(l\right)}=({\mathit{G}}^{\left(l\right)}{)}^{-1}{\mathit{y}}^{\left(l\right)}.$(8) ) for $0\le n\le N-1$. The problem is in that condition the graph G defined in (10(10) $G(\mathit{x},\mathit{g},L):=\left(V\right(\mathit{x}),E(L\left)\right),$(10) ) cannot be connected if there is only one l satisfy conditions of Theorem 2.4 when L>1. Therefore, we discuss multiple windows STFT to solve the conditions that $L\ge 1$.

2.2. Multiple discrete windows case ($R>1,L\ge 1$)

Multiple windows STFT PR provides more measurements compared with single window STFT PR. Therefore, the graph connectivity of G may be easier to satisfy.

Now we start to figure out which condition should be satisfied to ensure graph G is connected in STFT PR problem. Choose a window ${\mathit{g}}_r$ such that $W_r<N/2$. When $l=W_r$, there exists unique $n_r$ such that ${\tilde{g}}_r(n_r,W_r)\ne 0$, then $Y_r(m,W_r)=\sum _{n=0}^{N-1}\tilde{x}(n,W_r){\tilde{g}}_r(mL-n,W_r)=\tilde{x}(mL-n_0,W){\tilde{g}}_r(n_r,W_r).$ Therefore, (19) $\begin{array}{rl}{\rho }_{mL-n_r,mL-n_r+W_r}& =\frac{\tilde{x}(mL-n_r,W_r)}{|\tilde{x}(mL-n_r,W_r)|}={\left\{\frac{{\tilde{g}}_r(n_r,W_r)}{|{\tilde{g}}_r(n_r,W_r)|}\right\}}^{-1}\times \frac{Y_r(m,W_r)}{|Y_r(m,W_r)|},\\ m& =0,1,2,\dots ,N/L-1.\end{array}$(19) Note that only when $n=mL-n_r$ with $0\le m\le \lceil N/L\rceil -1$, ${\rho }_{n,n+W_r}$ can be computed by Equation (19(19) $\begin{array}{rl}{\rho }_{mL-n_r,mL-n_r+W_r}& =\frac{\tilde{x}(mL-n_r,W_r)}{|\tilde{x}(mL-n_r,W_r)|}={\left\{\frac{{\tilde{g}}_r(n_r,W_r)}{|{\tilde{g}}_r(n_r,W_r)|}\right\}}^{-1}\times \frac{Y_r(m,W_r)}{|Y_r(m,W_r)|},\\ m& =0,1,2,\dots ,N/L-1.\end{array}$(19) ). Therefore, single window STFT measurements cannot satisfy the connection of graph $G\left(V\right(\mathit{x}),E(L\left)\right)$ when L>1. If we want to compute all the relative phase in the $V\left(\mathit{x}\right)$, multiple windows $\{{\mathit{g}}_r{\}}_{r=1}^R$ are needed to compute the relative phase. As an instance shown in Figure 1, when N = 6, L = 2, $W_1=1$, $W_2=2$, the graph G is connected.

Figure 1. An example of connected graph $G\left(V\right(\mathit{x}),E(L\left)\right)$ with two sliding windows: $W_1=1$, $W_2=2$.

In this example, all the relative phase can be computed through (19(19) $\begin{array}{rl}{\rho }_{mL-n_r,mL-n_r+W_r}& =\frac{\tilde{x}(mL-n_r,W_r)}{|\tilde{x}(mL-n_r,W_r)|}={\left\{\frac{{\tilde{g}}_r(n_r,W_r)}{|{\tilde{g}}_r(n_r,W_r)|}\right\}}^{-1}\times \frac{Y_r(m,W_r)}{|Y_r(m,W_r)|},\\ m& =0,1,2,\dots ,N/L-1.\end{array}$(19) ). It shows that when L>1, the relative phase still can be obtained. Given two vertices $n_1$ and $n_2$, we want to figure out in which condition ${\rho }_{n_1,n_2}$ can be computed. Based on Theorem 2.4, we get some conclusions in multiple windows. Based on Equation (7(7) ${\mathit{y}}^{\left(l\right)}={\mathit{G}}^{\left(l\right)}{\mathit{x}}^{\left(l\right)}.$(7) ) and the solution identification theorem in linear equations, we get

Lemma 2.9

Denote the ${\mathit{G}}^{\left(l\right)}$ of ${\mathit{g}}_r$ in Equation (7(7) ${\mathit{y}}^{\left(l\right)}={\mathit{G}}^{\left(l\right)}{\mathit{x}}^{\left(l\right)}.$(7) ) as ${\mathit{G}}_r^{\left(l\right)}$, $r=1,2,\dots R$. Let ${\mathcal{G}}^{\left(l\right)}=[{{\mathit{G}}_1^{\left(l\right)}}^T{{\mathit{G}}_2^{\left(l\right)}}^T\cdots {{\mathit{G}}_R^{\left(l\right)}}^T{]}^T$. Then ${\mathit{x}}^{\left(l_0\right)}$ can be computed by ${\mathit{y}}^{\left(l_0\right)}={\mathcal{G}}_{l_0}{\mathit{x}}^{\left(l_0\right)}$ if $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathcal{G}}^{\left(l_0\right)}\right)=N$.

Based on Lemma 2.9, Theorem 2.10 gives an extension of Theorem 3.1 in [19] and Lemma 2.1.

Theorem 2.10

Let $\{{\mathit{g}}_r{\}}_{r=1}^R$ be a family of sliding windows with the same sliding step-size L, where L be a separation parameter with $N/L\in \mathbb{Z}$. Define ${\mathit{B}}_m^{\left(l\right)}\in {\mathbb{R}}^{R\times L}$, where ${\mathit{B}}_m^{\left(l\right)}(r,j)=\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+jN/L)$. ${\rho }_{n,n+l}$ can be computed if ${\mathit{B}}_m^{\left(l\right)}$ has rank L for all $0\le m\le N/L-1$.

Proof.

Since $Y_r(m,l)=\frac{1}{N}\sum _{k=0}^{N-1}|X_r(m,k){|}^2{\mathrm{e}}^{i2\pi kl/N}$ with $0\le m\le N/L-1$ and $1\le r\le R$. Denote ${\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}=\left(\tilde{x}\right(0,l),\tilde{x}(1,l),\dots ,\tilde{x}(N-1,l\left)\right)$, according to Fourier convolution property, we have $\begin{array}{rl}{Y}_r(m,l)& =\sum _{n=0}^{N-1}\tilde{x}(n,l){\tilde{g}}_r(mL-n,l)={\mathcal{F}}^{-1}\{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\cdot \mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\}\left(mL\right)\\ & =\frac{1}{N}\sum _{k=0}^{N-1}\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\left(k\right)\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\left(k\right){\mathrm{e}}^{2\pi imkL/N}.\end{array}$ Denote ${\mathit{Y}}_{\mathit{r}}\left(l\right)=\left(Y_r\right(0,l),Y_r(1,l),\dots ,Y_r(N/L-1,l\left)\right)$. Since ${\mathit{B}}_m^{\left(l\right)}=\left(\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\right(m+jN/L){)}_{1\le r\le R,0\le j\le L-1}$ and $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\left(k\right)=\frac{1}{N}\sum _{k=0}^{N-1}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\left(k\right){\mathrm{e}}^{-i2\pi kl/N}$. Then according to Fourier transform property, for $0\le m\le N/L-1$ and $1\le r\le R$, (20) $\begin{array}{rl}{\mathcal{F}{\mathit{Y}}_{\mathit{r}}}^{\left(l\right)}\left(m\right)& =\sum _{m^{\prime }=0}^{N/L-1}{Y}_r(m^{\prime },l){\mathrm{e}}^{-2\pi im{m}^{\prime }L/N}\\ & =\sum _{m^{\prime }=0}^{N/L-1}\left(\frac{1}{N}\sum _{k=0}^{N-1}\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\right(k\left)\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\right(k\left){\mathrm{e}}^{2\pi i{m}^{\prime }kL/N}\right){\mathrm{e}}^{-2\pi im{m}^{\prime }L/N}\\ & =\frac{1}{L}\sum _{k=0}^{N-1}\sum _{m^{\prime }=0}^{N/L-1}\frac{L}{N}\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\left(k\right)\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}\left(k\right){\mathrm{e}}^{2\pi i{m}^{\prime }(k-m)L/N}\\ & =\frac{1}{L}\sum _{j^{\prime }=0}^{L-1}\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}(m+j^{\prime }N/L)\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L).\end{array}$(20) Based on the equation of Lemma III.17 in [14], we get $\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}(m+jN/L)=L\sum _{j^{\prime }=0}^{L-1}\sum _{m^{\prime }=0}^{N/L-1}{b}_m(j,j^{\prime }){\mathrm{e}}^{-2\pi im{m}^{\prime }L/N}\left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}{Y}_r\right(m^{\prime },l\left)\right)$ where $\left(b_m\right(j,j^{\prime }){)}_{0\le j,j^{\prime }\le L-1}={({\mathit{B}}_m^{\left(l\right)}}^H{{\mathit{B}}_m^{\left(l\right)})}^{-1}$ for all $0\le m\le N/L-1$ and $0\le j\le L-1$.

Thus together with $x={\mathcal{F}}^{-1}\mathcal{F}x$, (21) $\begin{array}{rl}\tilde{x}(n,l)& ={\mathcal{F}}^{-1}\left\{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\right\}\left(n\right)\\ & =\frac{L}{N}\sum _{j^{\prime },j=0}^{L-1}\sum _{m,m^{\prime }=0}^{N/L-1}{\mathrm{e}}^{-2\pi i\left(m\right(m^{\prime }L-n)/N-jn/L)}{b}_m(j,j^{\prime })\left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}{Y}_r\right(m^{\prime },l\left)\right)\end{array}$(21) Then ${\rho }_{n,n+l}$ can be computed by the phase of $\tilde{x}(n,l)$. Since $\left({{\mathit{B}}_m^{\left(l\right)}}^H{\mathit{B}}_m^{\left(l\right)}\right)$ is invertible for all m, l, which means $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{B}}_m^{\left(l\right)}\right)=L$, then we finished the proof.

Remark 2.3

Now consider some special cases of the condition given by Theorem  2.10.

1. When L = 1, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{B}}_m^{\left(l\right)}\right)=L$ is equivalent to $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\sum _{r=1}^R|\mathcal{F}{\tilde{g}}_r\right(m,l\left){|}^2\right)=1$ for all $0\le m\le N-1$, which means there exists $r_0$, so that $\mathcal{F}{\tilde{g}}_{r_0}(n,l)\ne 0$, which is consistent with Theorem 3.1 in [19]. If R = 1, which means in single window case, $\mathcal{F}\tilde{g}(n,l)$ is non-vanishing.

2. When L = N, $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left({\mathit{B}}_m^{\left(l\right)}\right)=L$ is equivalent to $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\right\{{\tilde{g}}_r(n,l){\}}_{1\le r\le R,0\le n\le N-1})=N$, which is consistent with Corollary 2.9.

Based on Theorem 2.10, we can get theorems as follows:

Theorem 2.11

Suppose there exists $l_1,l_2,\dots ,l_t$ so that $(l_1,l_2,\dots ,l_t,N)=1$, and ${\mathit{B}}_m^{\left(l_i\right)}$ has rank L for $i=1,2,\dots ,t$, then the signal phase can be recovered.

Theorem 2.12

If there exists $l_0,l_1$ such that $l_1-l_0=(l_1,l_0)$, ${\mathit{B}}_m^{\left(l_0\right)}$ and ${\mathit{B}}_m^{\left(l_1\right)}$ have rank L, then the amplitude of the signal can be recovered.

Proof.

According to Theorem 2.10, if ${\mathit{B}}_m^{\left(l\right)}$ has rank L, ${\rho }_{n,n+l}$ can be computed by Equation (21(21) $\begin{array}{rl}\tilde{x}(n,l)& ={\mathcal{F}}^{-1}\left\{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\right\}\left(n\right)\\ & =\frac{L}{N}\sum _{j^{\prime },j=0}^{L-1}\sum _{m,m^{\prime }=0}^{N/L-1}{\mathrm{e}}^{-2\pi i\left(m\right(m^{\prime }L-n)/N-jn/L)}{b}_m(j,j^{\prime })\left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}{Y}_r\right(m^{\prime },l\left)\right)\end{array}$(21) ). Thus Theorem 2.11 and 2.12 can be inferred based on Theorems 2.5 and 2.6.

3. Reconstruction and error estimation

Theorems 2.11 and 2.12 extend retrieval conditions from a single window $\mathit{g}$ to multiple windows $\{{\mathit{g}}_r{\}}_{r=1}^R$, and the sliding size L = 1 to L>1. Based on these analyses, we propose the following reconstruction algorithm (PAR).

Here we discuss the computation complexity of PAR algorithm. To finish calculating $Y_r(n,l)$ and $\mathcal{F}\left\{{\tilde{g}}_r\right(n,l\left)\right\}$, we can use fast Fourier transforms with $\mathcal{O}(N\log N)$. Computing the rank of ${\mathit{B}}_m^{\left(l_i\right)}$ needs $\mathcal{O}\left(LRN\right)$ operations, while computing the highest common factor of no more than N number needs $\mathcal{O}(N\log N)$ by division algorithm in Theorem 2.11. Assume the number of l meeting conditions in Theorem 2.11 is k, then judging the condition of Theorem 2.12 needs at most $\mathcal{O}\left(kN\right)$. According to Equations (21(21) $\begin{array}{rl}\tilde{x}(n,l)& ={\mathcal{F}}^{-1}\left\{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\right\}\left(n\right)\\ & =\frac{L}{N}\sum _{j^{\prime },j=0}^{L-1}\sum _{m,m^{\prime }=0}^{N/L-1}{\mathrm{e}}^{-2\pi i\left(m\right(m^{\prime }L-n)/N-jn/L)}{b}_m(j,j^{\prime })\left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}{Y}_r\right(m^{\prime },l\left)\right)\end{array}$(21) ) and (14(14) $|x\left(l_0\right){|}^2=\frac{\prod _{k=0}^{l_0/M}|x\left(Mk\right)\overline{x(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x\left(Mk\right)\overline{x(Mk+l_1)}|}.$(14) ), computing relative phase and amplitude needs $\mathcal{O}\left(N\right)$ respectively. Thus the total runtime complexity of PAR algorithm is $\mathcal{O}(N\log N)$ in general.

Based on our phase retrieval procedure, we now present the following guarantee of stable performance. Suppose the measurements are corrupted by random noise $ϵ(m,k)$ with level $|ϵ|$ according to Equation (4(4) $Y_r(m,l)=\frac{1}{N}\sum _{k=0}^{N-1}|X_r(m,k){|}^2{\mathrm{e}}^{i2\pi kl/N}=\sum _{n=0}^{N-1}x\left(n\right)\overline{x(n+l)}{g}_r(mL-n)\overline{g_r(mL-n-l)},$(4) ), then the problem become recovering the signal $\mathit{x}$ from $|X_r(m,k){|}^2+{ϵ}_r(m,k)$. For a window family $\{{\mathit{g}}_r{\}}_{r=1}^R$ with period N extension, we define $\mathit{g}=[{\mathit{g}}_1,{\mathit{g}}_2,\dots ,{\mathit{g}}_R]$, ${\mathit{B}}^{\left(l\right)}=[{\mathit{B}}_0^{\left(l\right)},{\mathit{B}}_1^{\left(l\right)},\dots ,{\mathit{B}}_{N/L-1}^{\left(l\right)}]$, $\begin{array}{rl}\Vert \mathit{g}{\Vert }_2& ={\left({\mathrm{\Sigma }}_{r=1}^R{\mathrm{\Sigma }}_{n=0}^{N-1}|g_r\right(n\left){|}^2\right)}^{1/2},\\ \Vert {\mathit{B}}^{\left(l\right)}{\Vert }_1& ={\mathrm{\Sigma }}_{m=0}^{N/L-1}{\mathrm{\Sigma }}_{j,j^{\prime }=0}^{L-1}|b_{m,l}(j,j^{\prime })|.\end{array}$

Theorem 3.1

Let $C=\underset{n\in V\left(x\right)}{min}|x\left(n\right){|}^2$, $A=\underset{n\in V\left(x\right)}{max}|x\left(n\right){|}^2$. Denote $\Vert \mathit{B}\Vert =\underset{l\in \mathcal{L}}{max}\Vert {\mathit{B}}^{\left(l\right)}{\Vert }_1$, $T=(2\frac{l_0}{M}+1)\{\frac{A}{C}{\}}^{(\frac{l_0}{M}+1)}$, where $\mathcal{L}$, $l_0$ and M are defined in Equation (14(14) $|x\left(l_0\right){|}^2=\frac{\prod _{k=0}^{l_0/M}|x\left(Mk\right)\overline{x(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x\left(Mk\right)\overline{x(Mk+l_1)}|}.$(14) ) and Algorithm 1. Let $x_{ϵ}=\left(x_{ϵ}\right(0),\dots ,x_{ϵ}(N-1){)}^T$ be the approximation obtained by Algorithm 1 from noisy measurements $\{|{\mathit{X}}_r{|}^2+{\mathit{ϵ}}_{\mathit{r}}{\}}_{r=1}^R$, $|ϵ|=\underset{r,m,k}{max}|{ϵ}_r(m,k)|,$ then there exists $\text{β}\in \mathbb{R}$ such that for all $n=0,1,\dots ,N-1$, (22) $|\frac{x_{ϵ}\left(n\right)}{|x_{ϵ}\left(n\right)|}-{\mathrm{e}}^{2\pi i\text{β}}\frac{x\left(n\right)}{|x\left(n\right)|}|\le \frac{2NL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C},$(22) and (23) $||x_{ϵ}\left(n\right){|}^2-|x\left(n\right){|}^2|\le TL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|.$(23)

Proof.

For any $l\in 0,1,\dots ,N-1$, $\sum _{r=1}^R|\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}|\le \sum _{r=1}^R\sum _{k=0}^{N-1}|{\tilde{g}}_r(k,l)|=\sum _{r=1}^R\sum _{k=0}^{N-1}|g_r\left(k\right)\overline{\left(g_r\right(k-l\left)\right)}|\le \Vert \mathbf{g}{\Vert }_2^2.$ The last inequality comes from Cauchy–Schwarz inequality. Denote $Y_{r,ϵ}(m^{\prime },l)$ is the IDFT of $\{|{\mathit{X}}_{\mathit{r}}{|}^2+{\mathit{ϵ}}_{\mathit{r}}{\}}_{r=1}^R$ in (4(4) $Y_r(m,l)=\frac{1}{N}\sum _{k=0}^{N-1}|X_r(m,k){|}^2{\mathrm{e}}^{i2\pi kl/N}=\sum _{n=0}^{N-1}x\left(n\right)\overline{x(n+l)}{g}_r(mL-n)\overline{g_r(mL-n-l)},$(4) ). Combine with Equation (21(21) $\begin{array}{rl}\tilde{x}(n,l)& ={\mathcal{F}}^{-1}\left\{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{x}}}^{\left(l\right)}\right\}\left(n\right)\\ & =\frac{L}{N}\sum _{j^{\prime },j=0}^{L-1}\sum _{m,m^{\prime }=0}^{N/L-1}{\mathrm{e}}^{-2\pi i\left(m\right(m^{\prime }L-n)/N-jn/L)}{b}_m(j,j^{\prime })\left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}{Y}_r\right(m^{\prime },l\left)\right)\end{array}$(21) ), then we get (24) $\begin{array}{rl}|x_{ϵ}\left(n\right)\overline{x_{ϵ}(n+l)}-\tilde{x}(n,l)|& =\frac{L}{N}\sum _{j^{\prime },j=0}^{L-1}\sum _{m,m^{\prime }=0}^{N/L-1}{\mathrm{e}}^{-2\pi i\left(m\right(m^{\prime }L-n)/N-jn/L)}{b}_m(j,j^{\prime })\\ & \ast \left(\sum _{r=1}^R\overline{\mathcal{F}{\stackrel{\mathbf{~}}{\mathit{g}}}_r^{\left(l\right)}(m+j^{\prime }N/L)}\right(Y_{r,ϵ}(m^{\prime },l)-Y_r(m^{\prime },l)\left)\right)\\ & \le \frac{L}{N}\ast L^2\ast {\left\{\frac{N}{L}\right\}}^2\Vert {\mathit{B}}^{\left(l\right)}{\Vert }_1\Vert \mathit{g}{\Vert }_2^2\frac{1}{N}|ϵ|=L\Vert {\mathit{B}}^{\left(l\right)}{\Vert }_1\Vert \mathit{g}{\Vert }_2^2|ϵ|\\ & \le L\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|.\end{array}$(24) According to the error transfer formula: If $\overline{N}=f(\overline{x_1},\overline{x_2},\overline{x_3})$, then error estimation (25) $\mathrm{\Delta }N=\left|\frac{\mathrm{\partial }f}{\mathrm{\partial }{x}_1}\right|\mathrm{\Delta }{x}_1+\left|\frac{\mathrm{\partial }f}{\mathrm{\partial }{x}_2}\right|\mathrm{\Delta }{x}_2+\left|\frac{\mathrm{\partial }f}{\mathrm{\partial }{x}_3}\right|\mathrm{\Delta }{x}_3,$(25) where $\mathrm{\Delta }N,\mathrm{\Delta }{x}_1,\mathrm{\Delta }{x}_2,\mathrm{\Delta }{x}_3$ represent the error bound of each variable. (26) $\begin{array}{rl}||x_{ϵ}\left(n\right){|}^2-|x\left(n\right){|}^2|& \le |\frac{\prod _{k=0}^{l_0/M}|x_{ϵ}\left(Mk\right)\overline{x_{ϵ}(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x_{ϵ}\left(Mk\right)\overline{x_{ϵ}(Mk+l_1)}|}-\frac{\prod _{k=0}^{l_0/M}|x\left(Mk\right)\overline{x(Mk+l_0)}|}{\prod _{k=0}^{l_0/M-1}|x\left(Mk\right)\overline{x(Mk+l_1)}|}|\\ & \le (2\frac{l_0}{M}+1){\left\{\frac{A}{C}\right\}}^{(\frac{l_0}{M}+1)}L\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|=TL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|.\end{array}$(26) Then according the triangle inequality, we get that (27) $\begin{array}{rl}|\frac{x_{ϵ}\left(n\right)\overline{x_{ϵ}(n+l)}}{|x_{ϵ}\left(n\right)\overline{x_{ϵ}(n+l)}|}-\frac{\tilde{x}(n,l)}{|\tilde{x}(n,l)|}|& \le \frac{|x_{ϵ}\left(n\right)\overline{x_{ϵ}(n+l)}-\tilde{x}(n,l)|+||x_{ϵ}\left(n\right)\overline{x_{ϵ}(n+l)}|-|\tilde{x}(n,l)||}{|\tilde{x}(n,l)|}\\ & \le \frac{2L\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C}.\end{array}$(27) Since the relative phase propagation in graph G, we have the estimation equation (22(22) $|\frac{x_{ϵ}\left(n\right)}{|x_{ϵ}\left(n\right)|}-{\mathrm{e}}^{2\pi i\text{β}}\frac{x\left(n\right)}{|x\left(n\right)|}|\le \frac{2NL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C},$(22) ).

Let $\mathit{p},{\mathit{p}}_{\mathit{ϵ}}\in {\mathbb{C}}^N$ be the vectors of phases of the entries of ${\mathit{x}}_{\mathit{ϵ}}$ and $\mathit{x}$. ° represents Hadamard product. The reconstruction error (28) $\begin{array}{rl}d(\mathit{x},{\mathit{x}}_{\mathit{ϵ}})& =\underset{\varphi }{min}{\Vert e^{i\varphi }\mathit{x}-{\mathit{x}}_{\mathit{ϵ}}\Vert }_2\\ & \le {\Vert \mathit{x}|\circ \mathit{p}-|\mathit{x}|\circ {\mathit{p}}_{\mathit{ϵ}}\Vert {}_2+\Vert |\mathit{x}|\circ {\mathit{p}}_{\mathit{ϵ}}-|{\mathit{x}}_{\mathit{ϵ}}|\circ {\mathit{p}}_{\mathit{ϵ}}\Vert }_2\\ & =\sqrt{\sum _{n=1}^N{\left|x\right(n\left)\right|}^2{|\frac{x_{ϵ}\left(n\right)}{|x_{ϵ}\left(n\right)|}-\frac{x\left(n\right)}{|x\left(n\right)|}|}^2}+\sqrt{\sum _{n=1}^N{\Vert x\left(n\right)|-|x_{ϵ}\left(n\right)\Vert }^2}\\ & \le \Vert \mathit{x}{\Vert }_2\frac{2NL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C}+\sqrt{NTL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}.\end{array}$(28) Theorem 3.1 ensures that when error ε is in low level and the minimal of $x\left(n\right)$ is not small, the error between retrieval phase and exact phase can be controlled. In addition, the norm of $\mathit{g}$ and $\mathit{B}$ influences the upper bound of error. Therefore, the number of windows does not need to be too large but should meet the retrieval conditions. It means that if we choose some suitable windows, the signal can be recovered well by PAR.

4. Numerical results

We now show the results of PAR algorithm in different types of signals and sliding windows. Here we choose three different signals, Gaussian complex signal, chirp signal and real temperature signal. A Gauss complex signal $\mathit{x}\in {\mathbb{C}}^N$ with N = 101, which is defined as (29) $x\left(n\right)\sim \mathcal{N}(\mu ,1)+i\mathcal{N}(0,1),n=0,\dots ,N-1.$(29) Chirp signal is a typical non-stationary signal in communication, sonar, radar. Figure 2 shows a chirp signal with length N = 101 and average value is 0. Earth abnormal temperature signal is a real signal describing the temperature variation in Figure 3.

Figure 2. Chirp signal.

Figure 3. Earth abnormal temperature.

In addition, we choose some different typical sliding windows in STFT including rectangle windows, triangle windows and hamming windows. These windows have different main lobes and side lobes, which are suitable for different situations. All the sliding windows satisfy the retrieval condition in PAR method. These simulations display the relation between relative error and sliding window as well as signal-noise ratio (SNR), where the relative error between estimation ${\mathit{x}}_{\mathit{ϵ}}$ and real signal $\mathit{x}$ is defined as $\mathrm{R}\mathrm{e}({\mathit{x}}_{\mathit{ϵ}},\mathit{x})=\frac{d({\mathit{x}}_{\mathit{ϵ}},\mathit{x})}{\Vert \mathit{x}{\Vert }_2}.$ We reproduce our procedure in https://github.com/zhouxianchen/PAR-algorithm with python.

4.1. Relative error and length of window

Since LS method in [19] requires $W\ge \lceil \frac{N+1}{2}\rceil$, we first discuss the relationship between relative error and sliding window with two signals. Here we choose parameter $\mu =10$ and $\mathrm{S}\mathrm{N}\mathrm{R}=60$ in (29(29) $x\left(n\right)\sim \mathcal{N}(\mu ,1)+i\mathcal{N}(0,1),n=0,\dots ,N-1.$(29) ). Figure 4 and Figure 5 show PAR method has less than 0.025 relative error without considering the length of window W in three types of windows. And two figures show that PAR method has no restriction with length of window and can be used in wider range compared with the LS method. Relative error in chirp signal and temperature signal is relatively smaller than that in Gauss complex signal, which might be because that it is no relative phase error in the reconstruction procedure with real signal (Figure 6).

Figure 4. Relative error in recovering Gaussian complex signal by different windows using PAR.

Figure 5. Relative error in recovering chirp signal by different windows using PAR.

Figure 6. Relative error in recovering temperature signal by different windows using PAR.

4.2. Relative error and SNR

Here we discuss PAR performance in different SNR. We choose the length of window W = 7 and $\mu =10$ in (29(29) $x\left(n\right)\sim \mathcal{N}(\mu ,1)+i\mathcal{N}(0,1),n=0,\dots ,N-1.$(29) ). Figure 7 shows that PAR can retrieve signal when $\mathrm{S}\mathrm{N}\mathrm{R}\ge 40$. When $\mathrm{S}\mathrm{N}\mathrm{R}\le 40$, different windows have different performance in the same level SNR. The triangle window has better performance compared with rectangle and hamming windows. We infer it is due to that $\Vert \mathit{B}\Vert$ and $\Vert \mathit{g}{\Vert }_2$ in (22(22) $|\frac{x_{ϵ}\left(n\right)}{|x_{ϵ}\left(n\right)|}-{\mathrm{e}}^{2\pi i\text{β}}\frac{x\left(n\right)}{|x\left(n\right)|}|\le \frac{2NL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C},$(22) ) are influenced by the type of windows, which lead to different error estimation (Figures 8 and 9).

Figure 7. Relative error in recovering noisy Gaussian signal using PAR (W = 7).

Figure 8. Relative error in recovering noisy chirp signal using PAR (W = 7).

Figure 9. Relative error in recovering temperature signal using PAR (W = 7).

4.3. Relative error and μ

To verify our error estimation, we discuss the relationship between relative error and μ in (29(29) $x\left(n\right)\sim \mathcal{N}(\mu ,1)+i\mathcal{N}(0,1),n=0,\dots ,N-1.$(29) ) using a Gauss complex signal. Here we choose a Hamming window STFT. Figure 10 shows that when $\mathrm{S}\mathrm{N}\mathrm{R}\ge 40$, $\mu \ge 20$ can retrieve signal correctly. And when $\underset{n\in V\left(x\right)}{min}|x\left(n\right){|}^2$ is larger, PAR method has robuster retrieval performance. It verifies our error estimation formula (22(22) $|\frac{x_{ϵ}\left(n\right)}{|x_{ϵ}\left(n\right)|}-{\mathrm{e}}^{2\pi i\text{β}}\frac{x\left(n\right)}{|x\left(n\right)|}|\le \frac{2NL\Vert \mathit{B}\Vert \Vert \mathit{g}{\Vert }_2^2|ϵ|}{C},$(22) ) and inspires us when we deal with some data in reality, we can use a signal translation to enlarge that $\underset{n\in V\left(x\right)}{min}|x\left(n\right){|}^2$ (Figure 11).

Figure 10. Relative error in recovering noisy Gaussian signal using PAR in different μ.

Figure 11. Relative error in recovering Gauss complex signal with Hamming window using PAR and LS method (SNR=100).

4.4. Relative error between LS and PAR algorithm

In this section, we display the performance of PAR algorithm and the classical LS [19] algorithm under different length of window. PAR method has good performance when the length of W is small compared with the LS method. Therefore, PAR method has less restriction of the window length than the LS method.

As is shown in the experiments, no matter for Gaussian complex, chirp or real temperature signal, our method obtains better performance than the LS method. However, note that the retrieval conditions of PAR based on the assumption that the signal is non-vanishing and graph connectivity based on windows. Some signals and windows do not satisfy these assumptions, and our method is not suitable for those situations.

5. Conclusion and further work

Previous research proposed STFT PR methods in relatively strict retrieval conditions. This paper proposes a PAR reconstruction algorithm with a milder retrieval conditions, which solved the PR by computing phase and amplitude respectively. Compared with the previous method, PAR reconstruction Algorithm 1 has a milder constraint of windows and can be used to solve the STFT phase retrieval problems both in single window and multiple windows cases, when prior knowledge is known including that the noise is low and the minimal magnitude of nonzero components of original signal is enough large. Numerical results show that when noise is in low level, PAR has better performance and a translation shift of initial signal can lead to a robuster result. In addition, since the computation cost is $\mathcal{O}(N\log N)$ when l is selected in Algorithm 1, it can also be used in initialization of other non-convex method.

Disclosure statement

There are no relevant financial or non-financial competing interests to report. And no potential competing interest was reported by the authors.

Additional information

Funding

This work was supported in part by The National Key Basic Research Program (Grant No. 2020YFA0713504) and National Natural Science foundation (China): 61977065,11971489.