{'fraction_non_alphanumeric': 0.03972386356988241, 'fraction_numerical': 0.010881647458023752, 'mean_word_length': 4.643446682073291, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Super-Resolution of Mutually Interfering Signals

We consider simultaneously identifying the membership and locations of point sources that are convolved with different low-pass point spread functions, from the observation of their superpositions. This problem arises in three-dimensional super-resolution single-molecule imaging, neural spike sorting, multi-user channel identification, among others. We propose a novel algorithm, based on convex programming, and establish its near-optimal performance guarantee for exact recovery by exploiting the sparsity of the point source model as well as incoherence between the point spread functions. Numerical examples are provided to demonstrate the effectiveness of the proposed approach.


I. INTRODUCTION
In many emerging applications in applied science and engineering, the acquired signal at the sensor can be regarded as a superposition of returns from multiple channels (or users), where the return from each channel is governed by the underlying physical field that produced it, e.g. the Green's function, the point spread function, the signature waveform, etc. The goal is to invert for the field parameters of each channel that produced the acquired signal which reflects the ensemble behavior of all channels.
Mathematically, consider the acquired signal, y(t), given as where * denotes convolution, x i (t) = Ki k=1 a ik δ(t − τ ik ) is the point source signal observed through the ith channel, g i (t) is the point spread function of the ith channel, respectively, and I denotes the total number of channels. For the ith channel, let τ ik ∈ [0, 1) and a ik ∈ C be the location and the amplitude of the kth point source, 1 ≤ k ≤ K i , respectively. In typical applications, we are interested in resolving the point sources of each channel at a resolution much higher than that of the acquired signal y(t), determined by the Rayleigh limit, or in other words, the bandwidth of the point spread functions.
The proposed model (1) occurs in a wide range of practical problems, ranging from three-dimensional super-resolution single-molecule imaging, to spike sorting in neural recording and DNA sequencing, to multi-user multi-path channel identification in communication systems, and many others.
Three-dimensional super-resolution imaging: By employing photoswitchable fluorescent molecules, the imaging process of stochastic optical reconstruction microscopy (STORM) [1] is divided into many frames, where in each frame, a sparse number of fluorophores (point sources) are randomly activated, localized at a resolution below the diffraction limit, and deactivated. The final image is thus obtained by superimposing the localization outcomes of all the frames. This principle can be extended to reconstruct a 3-D object from 2-D image frames [2], by introducing a cylindrical lens to modulate the ellipticity of the point spread function based on the depth of the fluorescent object. Therefore, the acquired image in each frame can be regarded as a superposition of returns from multiple depth layers, where the return from each layer corresponds to the convolution outcome of the point sources in that depth layer with the depth-dependent PSF, as modeled in (1). The goal is thus to recover the locations and depth membership of each point source given the image frame that records the returns from all depth layers.
Spike sorting for neural recording: Neurons in the brain communicate by firing action potentials, i.e. spikes, and it is possible to hear their communications through a single or multiple microelectrodes, which record simultaneously activities of multiple neurons within a local neighborhood. Spike sorting [3], thus, refers to the grouping of spikes according to each neuron, from the recording of the microelectrodes. Interestingly, it is possible to model the spike fired by each neuron with a characteristic shape [4]. The neural recording can thus be modeled as a superposition of returns from multiple neurons, as in (1), where the return of each neuron corresponds to the convolution of its characteristic spike shape with the sequence of its firing times. A similar formulation also arises in DNA sequencing [5].
Multi-path identification in random-access channels: In multi-user multiple access model [6], the base station receives a superposition of returns from active users, as in (1), where the received signal component for each active user corresponds to the convolution of its signature waveform with the unknown sparse multi-path channel from the user to the base station. The goal is to identify the set of active users, as well as their channel states, from the received signal at the base station.

A. Related Work and Contributions
There's extensive research literature on inverting (1) when there is only a single channel with I = 1, where conventional methods such as matched filtering, MUSIC [7], matrix pencil [8], to more recent approaches based on total variation minimization [9], can be applied. However, these approaches can not be applied directly when multiple channels exist in the observed signal, due to the mutual interference between channels. To the best of the authors' knowledge, methods for inverting (1) with multiple channels have been extremely limited. In [6], [10], [11], sparse recovery algorithms have been proposed to invert (1) by assuming the source locations lie on a fine-grain grid, which is mismatched from the actual source locations, resulting in possibly performance degradation [12]. Even when all the point sources indeed lie on the grid, existing work [11] shows that the sample complexity needs to grow logarithmically with the size of the discretized grid, which is undesirable. The continuous basis pursuit algorithm [13] can be applied to retrieve source locations off the grid, but lacks performance guarantees.
In this paper, we study the problem of super-resolving (1) when there're two channels, i.e. I = 2. We start by recognizing that in the frequency domain, the observed signal can be regarded as a linear combination of two spectrally-sparse signals, each composed of a small number of complex sinusoids. We then separate and recover the two signals by motivating their spectral structures using atomic norm minimization, which has been recently shown as an efficient convex optimization framework to motivate parsimonious structures [14], [9], [15], as well as satisfying the observation constraints.
The separation and identification of the two channels, using the proposed algorithm, denoted by convex demixing, is made possible with two additional conditions. The first condition is that the point source models satisfy a mild separation condition, such that the locations of the point sources are separated by at least four times the Rayleigh limit; this is in line with the separation condition required by Candès and Fernandez-Granda [9] even with a single channel. The second condition is that the point spread functions of different channels has to be sufficiently incoherent for separation, which is supplied in our theoretical analysis by assuming they're randomly generated from a uniform distribution on the complex circle. We demonstrate that, as soon as the number of measurements, or alternatively, the bandwidth of the point spread functions, is on the order max(K 1 , K 2 ) log(K 1 + K 2 ) up to logarithmic factors, the proposed algorithm, denoted by convex demixing, recovers the locations of the point sources for each channel exactly, with high probability. Since at least an order of K 1 + K 2 measurements is necessary, our sample complexity is near-optimal up to logarithmic factors. Moreover, the point sources can be recovered from the dual solution of the proposed algorithm, without estimating or knowing the model order a priori. Our proof is based on constructing two related polynomials that certify the optimality of the proposed algorithm. The effectiveness of the proposed algorithm is demonstrated in numerical experiments.

B. Organization
The rest of this paper is organized as follows. The proposed algorithm based on convex programming and its performance guarantee are provided in Section II. The proof of the performance guarantee is sketched in Section III. Numerical experiments are shown in Section IV, and we conclude the paper in Section V. Throughout the paper, (·) T and (·) * denote the transpose and Hermitian transpose respectively,(·) denotes the element-wise conjugate, Tr (·) denotes the trace of a matrix, and toep (u) denotes a Toeplitz Hermitian matrix with u as its first column.
which satisfies g i, Typically, M is determined by the physics, such as the aperture of the imaging device or the steering array. Taking the DTFT of (1) with I = 2, we obtain with n ∈ Ω M . The measurements y n 's in (3) can be considered as a linear combination of two spectrally-sparse signals, with g in 's determining the combination coefficients. Multiplying both sides of (3) with g −1 1n , and with slight abuse of notation, we rewrite (3) into a vector form: where y = [y −2M , . . . , y 0 , . . . , y 2M ] and ⊙ denotes the Hadamard element-wise product operator.
Define the atomic norm [14], [15] of x ∈ C 4M+1 with respect to the atoms c(τ ) as which can be regarded as the tightest convex relaxation of counting the smallest number of atoms that is needed to represent a signal x. Therefore, we seek to recover the signals x 1 and x 2 by motivating their spectral sparsity via minimizing the atomic norm, with respect to the observation constraint: The above algorithm is referred to as convex demixing. Interestingly, (6) can be equivalently rewritten with the following semidefinite programming characterization [15], which can be solved efficiently using off-the-shelf solvers, as min x1,x2,u1,u2,t1,t2 Tr (toep (u 1 )) + Tr (toep (u 2 )) 4M + 1 Define the separation condition of point sources of each channel as which is the wrapped-around distance on [0, 1], and the minimum separation of all channels as ∆ = min i ∆ i . Our main theorem is stated below.
is sufficient to guarantee that x ⋆ 1 and x ⋆ 2 are the unique solutions of (6) with probability at least 1 − η.
Theorem 1 indicates that as soon as the number of measurements, or alternatively, the bandwidth of the point spread functions, is on the order M = O(max(K 1 , K 2 ) log(K 1 + K 2 ) log M ), the proposed convex demixing algorithm recovers the locations of the point sources for each channel exactly, with high probability. This suggests that the performance of the convex demixing algorithm is near optimal in terms of the sample complexity. Remark 1. The point source separation condition ∆ ≥ 1/M is used as a sufficient condition in Theorem 1 to guarantee accurate signal demixing. It is implied in [9] that a reasonable separation is also necessary to guarantee stable superresolution. Interestingly, no separation between point sources is necessary across channels, as long as their point spread functions are sufficiently incoherent.
Remark 2. Theorem 1 assumes g n 's are generated with uniformly random phase, which may be reasonable when g n 's can be designed, such as the spreading sequences in multi-user communications. This assumption may be further relaxed. The proof procedure reveals that same results can be obtained as long as g n 's satisfy E [ḡ n ] = E ḡ −1 n = 0 and C 1 ≤ |g n | ≤ C 2 .
Remark 3. Both sign (a 1k ) and sign (a 2k ) are required to be randomly generated, which we believe are technical requirements of the proof, and may be removed with finer proof techniques.
Define the inner product of two vectors as p, x = x * p and the real-valued inner product as p, x R = Re (x * p). Then the dual norm of · A can be represented as We have the dual problem of (6) written aŝ which comes from standard Lagrangian calculations. Construct two dual polynomials using the dual solutionp: then the point source locations Υ 1 and Υ 2 can be recovered without model order estimation, by identifying the parameter τ such at |P (τ )| = 1, and |Q(τ )| = 1, respectively.

III. PROOF SKETCH OF MAIN RESULT
In this section, we proceed to sketch the proof of Theorem 1. We first provide the optimality conditions using dual polynomials to certify the optimality of the solution of (6). Illuminated by [9], [15], where the dual polynomial is constructed by the squared Fejér's kernel, we propose a construction of dual polynomials which are composed of a deterministic term and a random perturbation termed introduced by interference between channels. Finally, the remainder is to show that the constructed dual polynomials satisfy the optimality conditions with high probability when the number of measurements M is large enough.

A. Optimality Conditions using Dual Polynomials
We can certify the optimality of the primal problem (6) using the following proposition. Proposition 1. The solution of (6)x 1 = x ⋆ 1 andx 2 = x ⋆ 2 is the the unique optimizer if there exists a vector p such that the dual polynomials P (τ ) and Q(τ ) constructed from it where the sign should be understood as the complex sign.
Proposition 1 suggests that P (τ ) and Q (τ ) are dual certificates to show {x ⋆ 1 , x ⋆ 2 } is the unique primal optimizer. Therefore if we can find a vector p to construct two dual polynomials P (τ ) and Q(τ ) in (9) satisfying (10), the proposed algorithm is guaranteed to recover the ground truth.

B. Dual Certificate Construction
Our construction of dual polynomials is inspired by [9], [15], based on use of the squared Fejér's kernel. However, since the two dual polynomials are coupled together, the construction is more involved. Define the squared Fejér's kernel [9] as The value of K (τ ) is nonnegative, attaining the peak at τ = 0 and decaying to zero rapidly with the increase of the absolute value of τ .
We define K g (τ ) and Kḡ (τ ) respectively as We then construct the dual polynomials P (τ ) and Q (τ ) as and where τ 1k ∈ Υ 1 and τ 2k ∈ Υ 2 . It is straightforward to validate that there exists some corresponding vector p such that (13) and (14) can be equivalently written as (9).
The above setting, if exists, immediately satisfies the first and third conditions in (10). The rest of the proof is then to, under the condition of Theorem 1, guarantee that a solution of (15) exists with high probability, and that when existing, they satisfy the second and forth conditions in (10) with high probability, therefore completing the proof.

C. Validation of Constructed Dual Polynomials
First we want to show that the solution of (15) exists with high probability. Since the deterministic terms in constructed dual polynomials are well-conditioned if the point source separation ∆ satisfies appropriate condition [15], by showing the random perturbations are small, we can have the following proposition to guarantee the invertibility of (15) with high probability when the number of measurements is large enough.
While the random perturbation in constructed P (l) (τ ) introduced by interference effect can be verified small enough with high probability, we have the following proposition to uniformly bound the distance between P (l) (τ ) andP (l) (τ ).
Then applying Bernstein's polynomial inequality [16] and similar techniques in [15, Lemma 4.13 and 4.14], we have the following proposition.
The proof of Theorem 1 is then complete since we have established that P (τ ) and Q(τ ) constructed in (13) and (14) are indeed valid dual certificates under the condition of Theorem 1.

IV. NUMERICAL EXPERIMENTS
We carry out numerical experiments to validate the performance of the proposed convex demixing algorithm. For a fixed M , we vary the spectral sparsity level of the two channels as K 1 and K 2 . For each pair of (K 1 , K 2 ), we first randomly generate a pair of point sources Υ 1 and Υ 2 that satisfy a separation condition ∆ ≥ 1/ (2M ), which is in fact a little smaller than the theoretical constraint, with the coefficients of the point sources i.i.d. drawn from the complex standard Gaussian distribution. For each Monte Carlo trial, we then randomly generate the point spread functions g n 's in the frequency domain with i.i.d. entries drawn uniformly from the complex unit circle, and perform the algorithm by solving (7) using CVX [17]. The algorithm is considered successful when the Normalized Mean Squared Error (NMSE) satisfies Fig. 1 shows the success rate of the proposed algorithm over 20 Monte Carlo trials for each cell, when M = 8 in (a) and M = 16 in (b), respectively.

V. CONCLUSIONS
We propose a convex optimization method based on atomic norm minimization to super-resolve two point source models from the measurements of their superposition, with each convolved with a different low-pass point spread function. It is demonstrated, with high probability, that the point source locations of each channel can be simultaneously determined perfectly, from an order-wise near-optimal number of measurements, under mild conditions. The proposed framework can be extended to handle more than two channels, whose details will be discussed elsewhere.