Recovery of Sparse 1-D Signals from the Magnitudes of their Fourier Transform

Introduction

Signal extraction from the autocorrelation, or equivalently, from the magnitude of the Fourier Transform is known as phase retrieval. This problem fundamentally arises in many practical systems such as X-ray crystallography , astronomical imaging , channel estimation, speech recognition etc, and has attracted a considerable amount of attention from researchers over the last few decades . Various algorithms have been proposed to retrieve phase information and a comprehensive survey of them can be found in .

For one-dimensional signals, since the mapping from signals to autocorrelation is not one-to-one, unique recovery is not possible in general. For any given Fourier transform magnitude, every possible phase corresponds to a different signal. Hence, additional prior information on the signal is required to limit the number of valid phase combinations. uses multiple structured illuminations, in which several patterns using different masks are collected to guarantee uniqueness.

We assume that the signal is sparse, i.e., the number of non-zero entries in the signal is much less compared to the length of the signal. This constraint greatly limits the number of possible phase combinations, and research has been done recently to exploit this feature . In many applications of phase retrieval, the signals encountered are naturally sparse. For example, astronomical imaging deals with the locations of the stars in the sky, electron microscopy deals with the density of electrons, and so on.

In this work, we prove that signals can be recovered from their autocorrelation with arbitrarily high probability under certain conditions. We prove this using dimension counting, based on the ideas used in for multidimensional signals. We also propose two non-iterative recovery algorithms to extract sparse signals from their autocorrelation. Note that the phase recovery problem is inherently non-convex, and relaxations similar to the ones used in are used to develop a convex-optimization based framework.

The paper is organized as follows. In Section 2, we discuss some properties of autocorrelation and spectral factorization which we use for signal extraction. In Section 3, we prove that signals can be recovered from their autocorrelation with very high probability under certain conditions. Non-iterative recovery algorithms are proposed for extraction of the signal from their autocorrelation in Section 4. Section 5 presents the simulation results and concludes the paper.

Theory

Let $\mathbf{x}=(x_{0},x_{1},....x_{n-1})$ be a real-valued signal of length $n$. Its autocorrelation, denoted by $\mathbf{a}=(a_{0},a_{1},....a_{n-1})$, is defined as

where $A(z)$ and $X(z)$ are the $z$-transforms of $\mathbf{a}$ and $\mathbf{x}$ respectively. Since $\mathbf{x}$ is real valued, $X(z)$ is a polynomial in $z$ with real coefficients and hence its zeros occur in conjugate pairs. Also, since $A(z)=A(z^{-1})$, if $z_{0}$ is a zero of $A(z)$, then $z_{0}^{-1}$ is also a zero. Hence, the zeros of $A(z)$ appear in quadruples of the form $(z_{0},z_{0}^{\star},z_{0}^{-1},z_{0}^{-\star})$.

The extraction of $\mathbf{x}$ from $\mathbf{a}$, or equivalently $X(z)$ from $A(z)$, is known as spectral factorization and deals with the distribution of these quadruples between $X(z)$ and $X(z^{-1})$. For every quadruple $(z_{0},z_{0}^{\star},z_{0}^{-1},z_{0}^{-\star})$, we can either assign $(z_{0},z_{0}^{\star})$ to $X(z)$ and $(z_{0}^{-1}$,$z_{0}^{-\star})$ to $X(z^{-1})$, or assign $(z_{0}^{-1}$,$z_{0}^{-\star})$ to $X(z)$ and $(z_{0},z_{0}^{\star})$ to $X(z^{-1})$. The total number of different valid factorizations hence is exponential in the number of such quadruples.

If two distinct finite-length real-valued signals $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ have the same autocorrelation, then there exists finite-length real-valued signals $\mathbf{g}$ and $\mathbf{h}$ such that

Let $F_{1}(z)$, $F_{2}(z)$, $G(z)$ and $H(z)$ be the $z$-transforms of the signals $\mathbf{f}_{1}$, $\mathbf{f}_{2}$, $\mathbf{g}$ and $\mathbf{h}$ respectively. Since $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$ have the same autocorrelation, (2) gives us

where $A(z)$ is the $z$-transform of the autocorrelation of $\mathbf{f}_{1}$ and $\mathbf{f}_{2}$. For every quadruple ($z_{0}$,$z_{0}^{\star}$,$z_{0}^{-1}$,$z_{0}^{-\star}$) which are zeros of $A(z)$, $(z_{0},z_{0}^{\star})$ has to be assigned to $F_{1}(z)$ or $F_{1}(z^{-1})$, and $F_{2}(z)$ or $F_{2}(z^{-1})$. Let $P_{1}(z)$, $P_{2}(z)$ and $P_{3}(z)$ be the polynomials constructed from such conjugate pairs of zeros which are assigned to $(F_{1}(z),F_{2}(z))$ and $(F_{1}(z),F_{2}(z^{-1}))$ and $(F_{1}(z^{-1}),F_{2}(z))$ respectively. Note that $P_{2}(z)=P_{3}(z^{-1})$. We have

and hence $F_{1}(z)$ and $F_{2}(z)$ can be written as

where $G(z)=P_{1}(z)$ and $H(z)=P_{2}(z)$, or equivalently

Unique Recovery

In this section, we establish the fact that within the class of signals with non-uniform support (defined later), there is a one-to-one mapping between signals and their autocorrelation almost surely.

If $f:\mathcal{A}\rightarrow\mathcal{B}$ is a map from $\mathcal{A}$ to $\mathcal{B}$, where $\mathcal{A}$ is a manifold of dimension $d_{a}$ and $\mathcal{B}$ is a manifold of dimension $d_{b}$, then the image of $f$ is measure zero in $\mathcal{B}$ if $d_{a}<d_{b}$.

Note that any signal of length $n$ can be represented as a vector in $\mathcal{R}^{n}$. Let $\mathbf{f}$ be a finite-length real-valued signal of length $n$. Let $I$ represent its support, defined as the set of locations where the $\mathbf{f}$ can have non-zero entries. We say that a signal $\mathbf{f}$ has uniform support if the indices of the elements belonging to the support are periodic, i.e., in an arithmetic progression. The size of the set $I$ denotes the sparsity of $\mathbf{f}$. Let $\mathcal{F}_{k}$ denote the set of signals with sparsity $k$. Observe that $\mathcal{F}_{k}$ is a manifold of dimension $k$.

Suppose $\mathbf{g}$ and $\mathbf{h}$ are finite-length real-valued signals with support set $I_{g}$ and $I_{h}$ of sparsity $k_{g}$ and $k_{h}$ respectively. If $\mathcal{F}_{gh}$ denotes the set of signals $\mathbf{g}*\mathbf{h}$, and $I_{gh}$ its support. Then

The set $\mathcal{F}_{gh}$ is a manifold of dimension $k_{g}+k_{h}-1$.

$I_{gh}$ has sparsity $k_{gh}\geq k_{g}+k_{h}-1$, with equality iff $\mathbf{g}$ and $\mathbf{h}$ have uniform support.

If $\mathbf{f}=\mathbf{g}*\mathbf{h}$, where $I$, the support of $\mathbf{f}$, is a subset of $I_{gh}$ with sparsity $k$. The set of such $\mathbf{f}$ is a manifold of dimension $k_{g}+k_{h}-1-\gamma$, where $\gamma=k_{gh}-k$.

We refer the readers to for the proof of $(i)$ and $(iii)$. $(ii)$ directly follows from the properties of convolution. ∎

Suppose $\mathbf{f}=\mathbf{g}*\mathbf{h}$, with $\mathbf{f}$ having non-uniform support where as $\mathbf{g}$ and $\mathbf{h}$ have uniform support, also has the additional property that $\mathbf{f}^{\prime}$ has non-uniform support. Then, the set of such signals is a manifold of dimension strictly lesser than $k_{g}+k_{h}-1-\gamma$.

The idea of the proof is similar to , based on dimension counting. We saw in Lemma 3.2 that the set of signals $\mathbf{f}$ which can be represented as $\mathbf{g}*\mathbf{h}$ with sparsity $k$ can be written as a manifold of dimension $k_{g}+k_{h}-1-\gamma$. The new set of constraints introduced by terms in $\mathbf{f}^{\prime}$ being result in a further reduction in dimension. Hence the set of such signals belong to a manifold of dimension strictly lesser than $k_{g}+k_{h}-1-\gamma$.

Signals can be uniquely recovered from their autocorrelation, or equivalently, from the magnitudes of their Fourier Transforms almost surely iff they have non-uniform support.

Let $\mathcal{F}^{\prime}_{k}$ be the set of all signals $\mathbf{f}$ with non-uniform support of sparsity $k$ which have another signal $\mathbf{f}^{\prime}$ with non-uniform support and same autocorrelation. Note that $\mathcal{F}_{k}^{\prime}$ is the set of signals of sparsity $k$ which cannot be recovered uniquely from their autocorrelation. Lemma 2.1 showed the existence of signals $\mathbf{g}$ and $\mathbf{h}$ such that

From Lemma 3.2, we note that the dimension of $\mathcal{F}^{\prime}_{k}$ is less than or equal to $k_{g}+k_{h}-1-\gamma$

Suppose $\mathbf{f}$ or $\mathbf{f}^{\prime}$ have uniform support, there will be no additional reduction in dimension. This case is equivalent to recovering a one-dimensional signal uniquely with no additional constraints, which is almost surely not possible. ∎

Recovery Algorithms

In this section, we develop two non-iterative recovery algorithms for the extraction of sparse signals from their autocorrelation.

Algorithm 1 is based on combinatorial analysis. We propose a method to recover the support of the signal from the support of the autocorrelation, and prove that recovery is possible with very high probability if the sparsity of the signal is $o(n^{1/3})$. Using this support knowledge, we show that signals can be recovered from the autocorrelation with very high probability.

Suppose $\mathbf{x}$ is a signal of length $n$ such that each element in $\mathbf{x}$ belongs to the support with a probability $\frac{s}{n}$, where $s=n^{\alpha},\alpha\leq 1$, independent of each other. Let $\mathbf{a}$ denote its autocorrelation, $k$ denote its sparsity and $D=\{d_{1},d_{2},.....d_{k}\}$ be the set of indices of the elements belonging to the support. Also, let $d_{ij}$ be defined as $|d_{i}-d_{j}|$ for $(i,j)=\{1,2,....k\}$. If $A$ is the set of indices of elements belonging to the support of the autocorrelation, then $A=\{\bigcup_{i,j}{d_{ij}}\}$. Note that $d_{i,i+1}$ is a geometric random variable with parameter $\frac{s}{n}$. Without loss of generality, let us assume $d_{k-1,k}\geq d_{12}$, otherwise we could just flip the signal and consider the flipped signal. Define $A_{1}=\{d_{ij}-d_{12}|d_{ij}\in A\}$ and $A_{2}=\{d_{ij}-d_{k-1,k}|d_{ij}\in A\}$.

The algorithm for signal recovery is described below. In what follows, we give a sequence of lemmas to justify various steps of the algorithm.

The sparsity $k$ of the signal satisfies $(1-\epsilon)s\leq k\leq(1+\epsilon)s$ with very high probability for any $\epsilon>0$, $n>n(\epsilon)$.

For three independent random variables $X_{1}$, $X_{2}$ and $X_{3}$ where $X_{1}$ and $X_{2}$ are geometric random variables with parameter $\frac{s}{n}$ , $P(X_{1}-pX_{2}=qX_{3})\leq\frac{s}{n}$ if $s=n^{\alpha},\alpha<1$ for $n>n(\epsilon)$, where $p$ and $q$ are integers.

$P(d_{k-1,k}-d_{12}\in A)\leq(1+\epsilon)\frac{s^{3}}{n}$ for any $\epsilon>0$, $n>n(\epsilon)$.

Note that the $d_{ij}$’s for $i\neq 1,j\neq k$ are independent of $d_{12}$ and $d_{k-1,k}$. Hence Lemma 4.2 can be applied and each term in the first summation can be upper bounded by $\frac{s}{n}$. Since $d_{k-1,k}<d_{ik}$ and $d_{12}>0$, all the terms in the second summation are zero. The terms in the third summation can be equivalently written as $P(d_{k-1,k}-2d_{12}=d_{2j})$, and Lemma 4.2 can be used to upper bound every term by $\frac{s}{n}$. Since $d_{1k}$ is the largest sum, we need not consider it in the summation. Hence, we get

$d_{12}$ and $d_{k-1,k}$ can be recovered from the autocorrelation with very high probability if $s=o(n^{1/3})$.

The first and second highest terms in $A$ are $d_{1k}$ and $d_{2k}$ respectively since $d_{12}\leq d_{k-1,k}$. Note that $d_{1k}-d_{2k}=d_{12}$, hence $d_{12}$ can be recovered from the autocorrelation. The only terms that can be higher than $d_{1,k-1}$ in $A$ are $\{d_{3k},d_{4k},.....d_{k-1,k}\}$. Note that $d_{2k}-d_{ik}=d_{2i}$, which belongs to $A$ for all $i=\{3,.....k-1\}$. So if $d_{2k}-d_{1,k-1}$ doesn’t belong to $A$, we can recover $d_{1,k-1}$ by considering the highest term which when subtracted from $d_{2k}$ produces a value which doesn’t belong to $A$. The probability of failure can hence be written as $P(d_{k-1,k}-d_{12}\in A)$ which goes to zero if $s=o(n^{\frac{1}{3}})$, as seen in Lemma 4.3. Hence both $d_{12}$ and $d_{k-1,k}$ can be recovered with very high probability if $s=o(n^{1/3})$. ∎

With the knowledge of $d_{12}$ and $d_{k-1,k}$, we can construct the sets $A_{1}$ and $A_{2}$. Consider the intersection of $A$ and $A_{1}$. All entries of the form $d_{2i}$ for $i=\{3,4,...k\}$ will survive trivially for any signal. Similarly, all entries of the form $d_{i,k-1}$ for $i=\{1....k-2\}$ will survive the intersection of $A$ and $A_{2}$ for any signal. If we subtract the survivors of the intersection of $A$ and $A_{2}$ from $d_{2,k-1}$, we get $d_{2i}$ for $i=\{3,4,...k-1\}$. Hence the elements $d_{2i}$ for $i=\{3,4,...k-1\}$ will survive $(A\cap A_{1})\cap(d_{2,k-1}-(A\cap A_{2}))$.

No other $d_{ij}$ will survive $(A\cap A_{1})\cap(d_{2,k-1}-(A\cap A_{2}))$ and hence the support can be recovered with very high probability if $s=o(n^{1/3})$

Suppose you choose $d_{ij}$ such that $i$ and $j$ are picked at random. The probability that $d_{ij}$ is a particular value can be upper bounded by $\frac{1}{n}$. For a non-trivial $d_{ij}$ in $A$ to survive $A\bigcap A_{1}$, $d_{ij}+d_{12}$ has to be in $A$. Similarly, $d_{2k}-d_{ij}$ and $d_{2,k-1}-d_{ij}$ has to be in A for it to survive $d_{2k}-A\bigcap A_{2}$. Using union bounds, we see that the probability of survival of some other $d_{ij}$ goes to when $s=o(n^{1/3})$. Note that we have information about $d_{k-1,k}$ upto $s=o(n^{1/3})$.

If no other elements survive, from $d_{2i}$ for $i=\{3,4,...k-1\}$, we can extract $d_{i,i+1}$ for $i=\{3,4,...k-2\}$ and since we already know $d_{12}$ and $d_{k-1,k}$, we have the support of the signal.

Suppose we have the support of the signal, $D=\{d_{1},d_{2},.....d_{k}\}$ being the indices of the elements belonging to the support. Define a pair $(d_{i},d_{j})$ as a good pair if they are the only pair separated by $|d_{i}-d_{j}|$. Note that for such a pair, $a_{|d_{i}-d_{j}|}=x_{d_{i}}x_{d_{j}}$

Consider a graph $G$ with $k$ vertices, each vertex representing an element of the support. Draw a weighted edge between every good pair, the weight being the value of the corresponding autocorrelation. If the graph $G$ has an odd cycle and is connected, then the signal can be extracted from the autocorrelation upto a global sign.

Consider an odd cycle with $2r-1$ vertices ${i_{1}},{i_{2}},...i_{2r-1}$. The term $\frac{x_{i_{1}i_{2}}x_{i_{3}i_{4}}....x_{i_{2r-1}i_{1}}}{x_{i_{2}i_{3}}...x_{i_{2r-2}i_{2r-1}}}$ gives $x_{i_{1}}^{2}$, from which $x_{i_{1}}$ can be extracted upto a sign, and from it the other terms in the odd cycle can be extracted using the weight corresponding to the edges. Since the graph is connected, all the other terms can be calculated.

The graph $G$ has an odd cycle and is connected with very high probability for $s=o(n^{1/3})$.

Pick any three vertices randomly. Choose any path of length $k-3$ from one of those vertices to cover all the remaining vertices randomly. If all the edges exists between the three vertices and the chosen $k-3$ length path exists, we are through. If any of the $k$ edges doesn’t exist, it implies that the distance between that pair of vertices occurs more than once. Since there are less than $k^{2}$ pairs, the probability of a pair of vertices not having an edge can be union bounded by $\frac{k^{2}}{n}$. Since there are $k$ edges to be considered, the probability of failure can be upper bounded by $\frac{k^{3}}{n}$. Hence if $s=o(n^{1/3})$, any chosen triangle and path exists with very high probability.

2 Algorithm 2

Algorithm 2 is developed using a convex optimization based framework. Semidefinite relaxation is used to convert the non-convex constraints into a set of convex constraints. We break the problem into two stages. First, the support of the signal is recovered from the autocorrelation and then we solve for the signal in the support.

We have to extract $\mathbf{u}$ from the autocorrelation of the signal. We will assume that the support of the signal is a subset of the support of the autocorrelation. This is the same as assuming there is no cancellation of support in the autocorrelation, which is a very weak requirement and holds with probability one if the coefficients of the signal are chosen randomly from a non-degenerate distribution. With this assumption, $a_{i}=0$ implies that no two elements in the support are separated by a distance $i$, and if $a_{i}$ is non-zero, there is atleast one pair of elements in the support separated by a distance $i$, i.e.,

where $\mathbf{u}$ is the binary support vector. This is clearly non-convex as the constraints are non-convex and $\mathbf{u}$ is binary. Define $\mathbf{S}=\mathbf{u}\mathbf{u}^{T}$, which is allowed to be positive semidefinite, as it is the smallest convex set containing all rank one matrices. The entries of $\mathbf{S}$ are allowed to be in $$, which is the best convex relaxation for binary variables. The trace of$\mathbf{S}$is given by$\sum_{i}{u_{i}^{2}}=\sum_{i}u_{i}=k$, the sparsity of the signal. Also, note that$\sum_{i}S_{ij}=\sum_{i}u_{i}u_{j}=u_{j}\sum_{i}u_{i}=ku_{j}=ku_{j}^{2}=kS_{jj}$and similarly$\sum_{j}S_{ij}=kS_{ii}$. Since flipped version of the support also satisfies all the constraints, a random matrix$\mathbf{V}$ is used to bias the cost. The support estimation problem becomes

Note that we assume apriori knowledge of the sparsity of the signal, i.e., the number of non-zero locations of the signal is known.

2.2 Signal Recovery

Note that the autocorrelation constraints are non-convex. As we did in the support extraction, we use the semidefinite relaxation $\mathbf{X}=\mathbf{x}\mathbf{x}^{T}$. We append $n$ zeros to the signal so that cyclic indexing scheme can be applied, hence a $m=2n$ order DFT matrix is required. Suppose $\mathbf{M_{k}}$ is the $m\times m$ matrix defined by $\mathbf{M}_{k}=\mathbf{f}_{k}\mathbf{f}_{k}^{T}$, where $\mathbf{f}_{k}$ is the $k^{th}$ column of the DFT matrix for $k=\{0,1,....m-1\}$. The autocorrelation constraints can be written in the Fourier domain as

where $\mathbf{Y}=\{|y_{0}|^{2},|y_{1}|^{2},......|y_{m-1}|^{2}\}$ is the vector containing the squared magnitude of the Fourier transform of $\mathbf{x}$. We can solve for the signal using L1-minimization .

Simulation Results

Figure 1 shows the success rate of signal recovery using Algorithm 1 as a function of the sparsity of the signal. We see that signals with $s=o(n^{1/3})$ are recovered successfully with very high probability. While the algorithm is computationally very cheap, it is not robust to noise due to error propagation.

Figure 2 demonstrates the performance of Algorithm 2 as a function of the sparsity of the signal. Numerical simulations strongly suggest that signals with sparsity upto $s=o(n^{1/2})$ can be recovered using this algorithm. It is also very robust to noise and hence more practical. We hope to provide theoretical guarantees in a future publication.

Appendix

For a pair of geometric random variables $X_{1}$ and $X_{2}$ with parameter $\frac{s}{n}$ each, $P(X_{1}-pX_{2}=c)\leq\frac{s}{n}$ if $s=n^{\alpha},\alpha<1$ for $n>n(\epsilon)$, where $p$ and $c$ are integers.

2 Proof of Corollary 4.2

References