User:NTUL/sandbox/STILIFT

Semidefinite relaxation-based algorithm for short time Fourier transform
The phase retrieval is an ill-posed problem. To uniquely identify the underlying signal, in addition to the methods that adds additional prior information like Gerchberg–Saxton algorithm, the other way is to add magnitude-only measurements like short time Fourier transform (STFT).

The method introduced below mainly based on the work of Jaganathan et al.

Short time Fourier transform
Given a discrete signal $$\mathbf{x}=(f[0],f[1],...,f[N-1])^T$$ which is sampled from $$f(x)$$. We use a window of length W: $$\mathbf{w}=(w[0],w[1],...,w[W-1])^T$$ to compute the STFT of $$\mathrm{f}$$, denoted by $$\mathbf{Y}$$:

$$Y[m,r]=\sum_{n=0}^{N-1}{f[n]w[rL-n]e^{-i2\pi \frac{mn}{N}}}$$

for $$0\leq m \leq N-1$$ and $$0\leq r \leq R-1$$, where the parameter $$L$$ denotes the separation in time between adjacent short-time sections and the parameter $$R = \left \lceil \frac{N+W-1}{L} \right \rceil$$ denotes the number of short-time sections considered.

The other interpretation (called sliding window interpretation) of STFT can be used with the help of discrete Fourier transform (DFT). Let $$w_r[n]=w[rL-n]$$ denotes the window element obtained from shifted and flipped window $$\mathbf{w}$$. Then we have

$$\mathbf{Y}=[\mathbf{Y}_0,\mathbf{Y}_1,...,\mathbf{Y}_{R-1}]$$, where $$\mathbf{Y}_r = \mathbf{x}\circ\mathbf{w}_r$$.

Problem definition
Let $${Z}_{w}[m,r]=|Y[m,r]|^2$$ be the $$N \times R$$ measurements corresponding to the magnitude-square of the STFT of $$\mathbf{x}$$, $$\mathbf{W}_{r}$$ be the $$N \times N$$ diagonal matrix with diagonal elements $$\left(w_{r}[0], w_{r}[1], \ldots, w_{r}[N-1]\right) .$$ STFT phase retrieval can be stated as:

Find $$\mathbf{x}$$ such that $$ Z_{w}[m, r]=\left|\left\langle\mathbf{f}_{m}, \mathbf{W}_{r} \mathbf{x}\right\rangle\right|^{2} $$ for $$0 \leq m \leq N-1$$ and $$0 \leq r \leq R-1$$, where $$\mathbf{f}_{m}$$ is the $$m$$-th column of the $$N$$-point inverse DFT matrix.

Intuitively, the computational complexity growing with $$N$$ makes the method impractical. In fact, however, for the most cases in practical we only need to consider the measurements corresponding to $$0 \leq m \leq M$$, for any parameter $$M$$ satisfying $$2W \leq M \leq N$$.

To be more specifically, if both the signal and the window are not vanishing, that is, $$x[n] \neq 0$$ for all $$0 \leq n \leq N-1$$ and $$w[n] \neq 0$$ for all $$0 \leq$$ $$n \leq W-1$$, signal $$\mathbf{x}$$ can be uniquely identified from its STFT magnitude if the following requirements are satisfied:


 * 1) $$L < W \leq N/2$$
 * 2) $$2W \leq M \leq N$$

The proof can be found in Jaganathan' s work.

The classic projection algorithm is proposed by Griffin and Lim, which reformulates STFT phase retrieval as the following least-squares problem:

$$\min _{\mathbf{x}} \sum_{r=0}^{R-1} \sum_{m=0}^{N-1}\left(Z_{w}[m, r]-\left|\left\langle\mathbf{f}_{m}, \mathbf{W}_{r} \mathbf{x}\right\rangle\right|^{2}\right)^{2}$$

The algorithm, although without  theoretical recovery guarantees, empirically able to converge to the global minimum when there is substantial overlap between adjacent short-time sections.

Semidefinite relaxation-based algorithm
To establish recovery guarantees, one way is to formulate the problems as a semidefinite program (SDP), by embedding the problem in a higher dimensional space using the transformation $$\mathbf{X}=\mathbf{x}\mathbf{x}^\ast$$ and relax the rank-one constraint to obtain a convex program. The problem reformulated is stated below:

Obtain $$\mathbf{\hat{X}}$$ by solving:$$\begin{align} &\mathrm{minimize}\mathrm{trace}(\mathbf{X})\\[6pt] &\mathrm{subject~to}Z[m,r]=\mathrm{trace}(\mathbf{W}_r^{\ast}\mathbf{f}_m\mathbf{f}_m^\ast\mathbf{W}_r\mathbf{X})\\[0pt] &~\mathbf{X}\succeq0 \end{align}$$for $$1 \leq m \leq M$$ and $$0 \leq r \leq R-1$$

Once $$\mathbf{\hat{X}}$$ is found, we can recover signal $$\mathbf{x}$$ by best rank-one approximation.