Blind Sparse Estimation of Intermittent Sources over Unknown Fading Channels

TL;DR

This paper introduces a novel blind source estimation algorithm for sparse, intermittent RF sources over unknown fading channels, combining dictionary learning and stochastic filtering to improve detection and estimation performance.

Contribution

It proposes a new two-stage algorithm integrating dictionary learning and stochastic filtering for blind estimation of sparse, intermittent sources without channel state information.

Findings

Enhanced detection and estimation performance of sources.

Robustness to different sparsity levels.

Effective handling of unknown fading channels.

Abstract

Radio frequency sources are observed at a fusion center via sensor measurements made over slow flat-fading channels. The number of sources may be larger than the number of sensors, but their activity is sparse and intermittent with bursty transmission patterns. To account for this, sources are modeled as hidden Markov models with known or unknown parameters. The problem of blind source estimation in the absence of channel state information is tackled via a novel algorithm, consisting of a dictionary learning (DL) stage and a per-source stochastic filtering (PSF) stage. The two stages work in tandem, with the latter operating on the output produced by the former. Both stages are designed so as to account for the sparsity and memory of the sources. To this end, smooth LASSO is integrated with DL, while the forward-backward algorithm and Expectation Maximization (EM) algorithm are leveraged for PSF. It is shown that the proposed algorithm can enhance the detection and the estimation performance of the sources, and that it is robust to the sparsity level.

Equations

$$\mathbf{Y}=\mathbf{HX}+\mathbf{Z},$$

$$\underset{\mathbf{H},\mathbf{X}\in \mathcal{X}}{\mathrm{minimize}}||\mathbf{Y}-\mathbf{HX}|{|}^2,$$

$${\mathbf{X}}^{(k+1)}=\underset{\mathbf{X}\in \mathcal{X}}{\mathrm{argmin}}||\mathbf{Y}-{\mathbf{H}}^{(k)}\mathbf{X}|{|}^2.$$

$$\underset{\mathbf{x}(t)}{\mathrm{minimize}}||\mathbf{y}(t)-{\mathbf{H}}^{(k)}\mathbf{x}(t)|{|}_2+\lambda ||\mathbf{x}(t)|{|}_1,t=1,\dots ,T,$$

$${\mathbf{H}}^{(j,k)}=\mathbf{Y}{\mathbf{X}}^{(j,k)T}({\mathbf{X}}^{(j,k)}{\mathbf{X}}^{(j,k)T}{)}^{-1},$$

$${\mathbf{x}}^{(j+1,k)}(t)={\mathbf{D}}^{(k)}({\mathbf{H}}^{(j,k)T}{\mathbf{H}}^{(j,k)}{)}^{-1}{\mathbf{H}}^{(j,k)T}\mathbf{y}(t),$$

$$\underset{\mathbf{X}}{\mathrm{minimize}}||\mathbf{Y}-{\mathbf{H}}^{(k)}\mathbf{X}|{|}_2^2+\lambda \sum _{t=1}^T||\mathbf{x}(t)|{|}_1+\mu \sum _{t=2}^T||\mathbf{x}(t)-\mathbf{x}(t-1)|{|}^2,$$

$$\underset{\mathbf{x}(t)}{\mathrm{minimize}}||\mathbf{y}(t)-{\mathbf{H}}^{(k)}\mathbf{x}(t)|{|}_2^2+\lambda ||\mathbf{x}(t)|{|}_1+\mu ||\mathbf{x}(t)-\mathbf{x}(t-1)|{|}^2,$$

$$\underset{\mathbf{X},{\mathbf{X}}^{\prime }}{\mathrm{minimize}}$$

$${\mathbf{x}}^{\prime }(t)=\mathbf{x}(t),t=1,\dots ,T-1,$$

$$\mathcal{L}(\mathbf{x},{\mathbf{x}}^{\prime },\mathbf{\alpha })=l_1(\mathbf{x}(1))+\sum _{t=2}^T{l}_t(\mathbf{x}(t),{\mathbf{x}}^{\prime }(t-1))+\sum _{t=1}^{T-1}{\mathbf{\alpha }}^T(t)(\mathbf{x}(t)-{\mathbf{x}}^{\prime }(t))+\rho \sum _{t=2}^T||\mathbf{x}(t)-{\mathbf{x}}^{\prime }(t-1)|{|}^2$$

$$=l_1(\mathbf{x}(1))+{\mathbf{\alpha }}^T(1)\mathbf{x}(1)+\sum _{t=2}^{T-1}(l_t(\mathbf{x}(t),{\mathbf{x}}^{\prime }(t-1))+{\mathbf{\alpha }}^T(t)\mathbf{x}(t)-{\mathbf{\alpha }}^T(t-1){\mathbf{x}}^{\prime }(t-1))$$

$$+l(\mathbf{x}(T),{\mathbf{x}}^{\prime }(T-1))-{\mathbf{\alpha }}^T(T-1){\mathbf{x}}^{\prime }(T-1)+\rho \sum _{t=2}^T||\mathbf{x}(t)-{\mathbf{x}}^{\prime }(t-1)|{|}^2,$$

$${\mathbf{\alpha }}^{(i+1)}(t)\leftarrow {\mathbf{\alpha }}^{(i)}(t)+\rho ({\mathbf{x}}^{(i)}(t)-{\mathbf{x}}^{\prime (i)}(t-1)),$$

$$\text{Pr}(s_n(t)=1|{\tilde{s}}_n(1),\dots ,{\tilde{s}}_n(T))={\alpha }_t{\beta }_t,$$

$${\hat{x}}_n(t)=\{\begin{array}{ll}{\tilde{x}}_n(t)& \text{if }{\hat{s}}_n(t)=1,\\ 0& \text{if }{\hat{s}}_n(t)=0.\end{array}$$

$$p_n^{(\nu +1)}\leftarrow \frac{{\sum }_{t=2}^T\text{Pr}(s_n(t)=1,s_n(t-1)=0|{\tilde{\mathbf{s}}}_n)}{{\sum }_{t=2}^T\text{Pr}(s_n(t-1)=0|{\tilde{\mathbf{s}}}_n)},$$

$$q_n^{(\nu +1)}\leftarrow \frac{{\sum }_{t=2}^T\text{Pr}(s_n(t)=0,s_n(t-1)=1|{\tilde{\mathbf{s}}}_n)}{{\sum }_{t=2}^T\text{Pr}(s_n(t-1)=1|{\tilde{\mathbf{s}}}_n)},$$

$$p_n^{\prime (\nu +1)}\leftarrow \frac{{\sum }_{t=1}^T\text{mathbbm}1({\tilde{s}}_n(t)=1)\text{Pr}(s_n(t)=0|{\tilde{\mathbf{s}}}_n)}{{\sum }_{t=1}^T\text{Pr}(s_n(t)=0|{\tilde{\mathbf{s}}}_n)},$$

$$q_n^{\prime (\nu +1)}\leftarrow \frac{{\sum }_{t=1}^T\text{mathbbm}1({\tilde{s}}_n(t)=0)\text{Pr}(s_n(t)=1|{\tilde{\mathbf{s}}}_n)}{{\sum }_{t=1}^T\text{Pr}(s_n(t)=1|{\tilde{\mathbf{s}}}_n)},$$

$$\text{Pr}(s_n(t)=i,s_n(t-1)=j|{\tilde{s}}_n(1),\dots ,{\tilde{s}}_n(T))={\alpha }_{t-1,j+1}{\Phi }_{j+1,i+1}^{(\nu )}\text{Pr}({\tilde{s}}_n(t)|s_n(t)=i){\beta }_{t,i+1},i,j=0,1,$$

$${\mathbf{\Phi }}_n^{(\nu )}=\left[\begin{array}{cc}1-p_n^{(\nu )}& p_n^{(\nu )}\\ q_n^{(\nu )}& 1-q_n^{(\nu )}\end{array}\right].$$

$$\text{EVM}(\%)=\sqrt{\frac{{\sum }_{n=1}^N{\sum }_{t\in {\mathcal{D}}_n}(x_n(t)-{\tilde{x}}_n(t){)}^2}{{\sum }_{n=1}^N{\sum }_{t\in {\mathcal{D}}_n}{x}_n^2(t)}}\times 100\%,$$

Full text

I Introduction

Blind source separation (BSS) refers to the separation of a set of source signals from a set of mixed signals, without resorting to any a priori information about the source signals or the mixing process [1]. BSS exploits only the information carried by the received signals themselves, hence the term blind. BSS has numerous applications in speech recognition [2, 3], image extraction [4, 5], and surveillance [6, 7]. Different metrics are used to evaluate the performance of BSS methods depending on the applications. For example, signal-to-interference ratio is used in [2] for speech separation, and a performance index is introduced in [4] for image feature extraction. Based on these metrics, many approaches have been proposed to solve BSS problems, such as independent component analysis (ICA) [8], principal component analysis (PCA) [9], and singular value decomposition (SVD) [10].

This paper addresses BSS in wireless networks. We are specifically interested in the set-up illustrated in Fig. 1, in which a fusion center observes a number of radio sources via noisy sensor measurements over unknown channels. The system may model an Internet-of-Things (IoT) system, such as LoRa, Sigfox, or Narrow Band-IoT (NB-IoT) [11, 12].

In wireless networks involving multiple terminals operating over flat fading channels, the signals received at a terminal are linear mixtures of the signals emitted by the transmitting terminals. The need for BSS arises in non-collaborative applications in which the signals and the channels through which they are received at a terminal are both unknown. ICA has been widely applied to solve BSS problems in wireless networks [13, 14, 15, 16], since it yields a useful decomposition with only scaling, and permutation ambiguities [17]. To achieve signal separation, ICA relies on the statistical independence and on the non-Gaussian distribution of the components of the mixture. Key assumptions made in the implementation of the various forms of ICA are that the underlying mixing process has the same number of inputs and outputs, and that all sources are active throughout the observation interval. These assumptions are limiting and not suitable for the applications under study in this work, as discussed next.

With the aim of capturing key aspects of IoT systems, this paper focuses on practical wireless scenarios in which the number of latent sources is generally larger than the number of sensors, but the sources are active intermittently with bursty transmission patterns. The cumulative time a source is active is a small fraction of the overall observation time, and the sources’ on-off patterns vary slowly. In order to capture these properties, the sources are modeled as hidden Markov processes with known or unknown parameters. Source memory will be seen to be instrumental in enabling source separation. The wireless network operates over slow flat-fading channels; sensors communicate with a fusion center over ideal channels; and all nodes are time-synchronized to the same clock by the fusion center.

The general BSS problem with more sources than sensors can be formulated as an underdetermined linear system $\mathbf{A}\mathbf{x}=\mathbf{y}$ in the absence of noise. The columns of the matrix $\mathbf{A}$ serve as a basis for expressing the observations $\mathbf{y}$. The set of basis signals that form the matrix $\mathbf{A}$ is called a dictionary. Underdetermined linear systems of equations of the form $\mathbf{A}\mathbf{x}=\mathbf{y}$ have infinitely many solutions when the matrix $\mathbf{A}$ is full rank. Regularization may introduce additional conditions on the solutions, for example favoring smaller values of $\mathbf{x}$, leading to unique solutions of the underdetermined linear system. Sparse representations for which the solution $\mathbf{x}$ is unique have been the subject of intensive research resulting in a large body of literature [18, 19, 20]. Pre-defined dictionaries, such as based on Fourier transforms, are convenient and computationally fast, but in the BSS problem, both the dictionary $\mathbf{A}$ and the signals $\mathbf{x}$ are unknown.

Sparse representation problems for which the dictionary is unknown require dictionary learning (DL) in addition to signal recovery [21, 22]. The advantage of DL is to enable a system to learn a dictionary adaptively from a set of observations rather than assume a prescribed rigid dictionary. Among other problems, DL methods have been applied to joint direction of arrival (DOA) estimation and array calibration [23], linear transceiver design [24], cloud K-singular value decomposition (K-SVD) for big, distributed data [25], and channel representation for frequency-division duplexing (FDD) massive MIMO system [26]. In these examples, the DL algorithms solve BSS problems in which the number of sources is larger than the number of sensors, but the methods are agnostic to time variability properties of sources with memory. In [27] and [28], the authors set up a hidden Markov model (HMM) to solve a BSS problem. However, in [27] the memory of the sources is not accounted for, while in [28] a simplified model is assumed whereby only one source can appear or disappear at any given time.

In this context, we propose a two-stage algorithm for solving the BSS problem for sources with memory modeled by an HMM and observed over slow flat-fading channels. The proposed algorithm comprises a DL stage and a Per-source Stochastic Filtering (PSF) stage. The DL stage of the proposed algorithm exploits knowledge about the source sparsity and memory to aid with the source separation. The effect of source memory is introduced by a penalty term that discourages solutions with short-duration transmissions by means of a smooth LASSO algorithm. The input to the DL stage are observations from the sensors. The output from the DL stage are channel estimates, source signal estimates and source states (active or inactive). The source state estimates produced by the DL algorithm may be viewed as the output of a binary asymmetric channel in which some of the states are ”flipped” with respect to the true states. The error probabilities associated with the flipped states are referred to as flipping probabilities. The PSF stage consists of a forward-backward step along with an Expectation Maximization (EM) step that estimate the unknown HMM transition probabilities and flipping probabilities.

The main contributions of this paper are summarized as follows:

• •

A two-stage architecture is introduced for solving the problem of blind source estimation of HMM sources over slow flat-fading channel. The sources feature intermittent activity, and the number of latent source may be larger than the number of sensors;

• •

A smooth DL algorithm is proposed that is capable of exploiting source memory to support channel estimation, signal estimation and source detection. Two simplified versions of the smooth DL algorithms are introduced to reduce the computational complexity;

• •

An PSF algorithm is introduced that is capable of operating in the absence of a priori information about the HMM parameters and the state estimation flipping probabilities. This algorithm integrates an forward-backward step with an EM step.

The rest of the paper is organized as follows. The system model and a hidden Markov source model are presented in Sec. II. Background on DL algorithms is provided in Sec. III. In Sec. IV, we propose a two-stage DL-based algorithm to solve the BSS problem for sources with memory in wireless networks. The PSF is described in detail in Sec. V for the two cases with known and unknown source parameters. Finally in Sec. VI, simulation numerical results are shown to support that our proposed algorithm can separate sources and recover source signals with higher accuracy than existing DL algorithms.

II System Model

Consider a system that includes $M$ receiving antennas, or radio sensors, and $N$ sources, as illustrated in Fig. 1. The number of sources $N$ is generally larger than the number of receive antennas $M$. All sensors are connected to a fusion center, which may be implemented in a cloud processor, via backhaul links. Equivalently, the fusion center has access to $N$ receive antennas. Models with a fusion center reflect the architecture of IoT networks, such as LoRa, Sigfox, and NB-IoT [11, 12]. In our model, each source transmits intermittently, and hence is generally active only for a subset of the $T$ symbol periods over which the sensors collect data.

Assuming that all nodes are time-synchronous, the discrete-time signal received by the $M$ sensors over $T$ symbol periods is given in matrix form as

[TABLE]

where $\mathbf{Y}=[\mathbf{y}(1),\dots,\mathbf{y}(T)]$ is an $M\times T$ matrix collecting as columns the $M\times 1$ received signals $\mathbf{y}(t)$ across all $T$ symbols $t=1,\dots,T$; $\mathbf{X}=[\mathbf{x}(1),\dots,\mathbf{x}(T)]$ is an $N\times T$ matrix that gathers the $N\times 1$ signals $\mathbf{x}(t)$ transmitted from all $N$ sources over time; $\mathbf{Z}=[\mathbf{z}(1),\dots,\mathbf{z}(T)]$ contains independent zero-mean complex Gaussian noise entries with variance $\sigma^{2}$; and $\mathbf{H}$ is the $M\times N$ complex fading channel matrix. The channel matrix $\mathbf{H}$ is assumed to be constant for $T$ symbol periods.

Given the intermittent nature of the traffic pattern of the sources, the $t$-th column vector $\mathbf{x}(t)$ that collects the $M$ symbols transmitted by the sources at time $t$, is generally sparse. In other words, only a subset of the entries of $\mathbf{x}(t)$ is non-zero. The signals $\mathbf{y}(t)$, for $t=1,\dots,T,$ are collected at the fusion center. Based on the received signals $\mathbf{Y}$, the goal of the fusion center is to detect the sources’ activity and to recover the signals $\mathbf{x}(t),$ for $t=1,\dots,T$, or equivalently the matrix $\mathbf{X}$, in the absence of information about the channel matrix $\mathbf{H}$.

For each source $n$, we define the activation pattern as a binary sequence $s_{n}(t)$, where the binary state $s_{n}(t)$ indicates whether a source is active or not. Specifically, when the state of source $n$ is $s_{n}(t)=1$, then the source is active, while it is inactive when $s_{n}(t)=0$. The binary state is assumed to be described by the following models.

Intermittent and smooth deterministic model: Each source is active for a small fraction of time, and its on-off patterns tend to have few switches between on and off states. As a result, the sequence $s_{n}(t)$ is “smooth”, i.e., it has a small number of transitions between on and off states; 2. 2.

Probabilistic hidden Markov model: As illustrated in Fig. 2, the activity $s_{n}(t)$ of each source $n$ follows a two-state Markov chain. The transition probabilities of the two-state Markov chain are defined as $p_{n}=\text{Pr}(s_{n}(t)=1|s_{n-1}(t)=0)$ and $q_{n}=\text{Pr}(s_{n}(t)=0|s_{n-1}(t)=1)$. We will consider both the cases in which the probabilities $p_{n}$ and $q_{n}$ are known and unknown.

We emphasize that the two models are not mutually exclusive and can be assumed to hold simultaneously. In particular, we will propose to leverage the deterministic model in order to improve the performance of source separation and the stochastic model to refine the estimates obtained by sources separation.

When a source $n$ is active, i.e., when $s_{n}(t)=1$, it transmits an independent sample $x_{n}(t)\sim f_{n}(t)$ with a given distribution $f_{n}(t)$, e.g., Gaussian or binary. Instead, when the source $n$ is inactive, i.e., $s_{n}(t)=0$, it does not transmit and we have $x_{n}(t)=0$.

A final remark is in order concerning synchronization requirements. The described model (1) applies even in the absence of time synchronization among the sources, since the vector $\mathbf{x}(t)$ can model a generic sample at discrete time $t$ of the transmitted signals. However, if time synchronization is assumed, and the sources transmit digitally modulated signals, the vector $\mathbf{x}(t)$ can be assumed to contain the constellation points transmitted at the $t$-th symbol period. We note that time and frequency synchronization in sensor networks is a topic of great interest, and it has been investigated in blind scenarios with no pilot symbols as well [29, 30]. However, the literature about BSS often assumes that sensors and sources are synchronized in order to focus on the technical challenges of BSS [23, 24, 25, 26].

III Preliminaries: Dictionary Learning

The DL method proposed in this work leverages prior information about the memory of the sources. We start by reviewing DL methods that do not exploit such information. These methods use only the fact that the signal $\mathbf{x}(t)$ is sparse at any time $t$. Prior information about the memory of each source $x_{n}(t)$, to be considered in the next section, includes smoothness properties or statistical models.

Assuming only information about the sparseness of $\mathbf{x}(t)$ at each time $t$, a standard approach is to utilize the channel matrix $\mathbf{H}$ as a dictionary to be learned to recover $\mathbf{X}$. DL techniques approximate the solution of the maximum likelihood (ML) problem if the fusion center acts as a master clock synchronizing all other nodes.

[TABLE]

where $\mathcal{X}$ is the set of matrices with sparse columns, that is, with columns containing a limited number of non-zero entries. This problem is not convex with respect to the pair ($\mathbf{H}$, $\mathbf{X}$). DL methods use an iterative procedure, whereby the signal $\mathbf{X}\in\mathcal{X}$ and the channel $\mathbf{H}$ are optimized alternately [31]. In the following, we first discuss solutions for the optimization over the signal $\mathbf{X}$ for a given channel matrix $\mathbf{H}$, and then over the channel $\mathbf{H}$ for a given signal matrix $\mathbf{X}$.

DL methods are subject to inherent permutation and sign ambiguities [32]. The scaling ambiguity can be instead resolved if one imposes that the channel matrix columns are normalized [33, 34, 35, 36].

III-A Signal Estimation

For any fixed iterate $\mathbf{H}^{(k)}$ at the $k$-th iteration, from (2), the problem of estimating the signal $\mathbf{X}$ reduces to

[TABLE]

Standard sparse optimization estimators, such as orthogonal matching pursuit (OMP) [18, §3.1.2] can be used to address problem (3). Alternatively, one can use the LASSO algorithm [37] to solve the convex problem

[TABLE]

separately for each $t$, where the weight $\lambda$ is a parameter to be determined as a function of the sparsity of vector $\mathbf{x}(t)$.

III-B Channel Estimation

At the $k$-th iteration, for a fixed iterate $\mathbf{X}^{(k+1)}$, the channel estimation step can obtain the next channel iterate $\mathbf{H}^{(k+1)}$ by using different algorithms, such as the Method of Optimal Directions (MOD) [33], the Multiple Dictionary Update (MDU) [34], or their enhanced versions proposed in [36]. Here we summarize the enhanced MDU approach, which was shown in [36] to provide the best performance via simulation results.

The MDU approach estimates the channel matrix $\mathbf{H}$ for a given $\mathbf{X}^{(k+1)}$ by following an iterative approach. To elaborate, denote $\mathcal{S}(\mathbf{x})$ the set of indices of the non-zero elements in vector $\mathbf{x}$. Also, index $(j,k)$ the $j$-th iteration of the MDU algorithm within the $k$-th step of the DL alternate optimization scheme. At the $(j,k)$ iteration, for $j=1,2,\dots$, MDU computes

[TABLE]

and

[TABLE]

where for all $t$, $\mathbf{D}^{(k)}$ is a diagonal matrix with elements having indices in $\mathcal{S}(\mathbf{x}^{(k+1)}(t))$ equal to $1$ and zero otherwise. The iteration is initialized with $\mathbf{X}^{(1,k)}=\mathbf{X}^{(k+1)}$. Therefore, for a fixed sparsity pattern $\mathcal{S}(\mathbf{x}^{(k+1)}(t))$, MDU alternately estimates channel and signals.

The enhancement proposed in [36] substitutes at iteration $k$ the received signal $\mathbf{Y}$ in (5) and (6) with $\mathbf{Y}^{(k)}=\mathbf{Y}+\mathbf{H}^{(k)}\mathbf{X}^{(k)}-\mathbf{H}^{(k)}\mathbf{X}^{(k+1)}$.

A final remark is that the reviewed DL techniques are designated for flat-fading channels. Multipath scenarios would require convolutional dictionary learning [38], and would be the topic of a future study.

IV Smooth Dictionary Learning

In this section, we propose an improved DL-based source separation algorithm that exploits knowledge about the memory of the sources $n=1,\dots,N$. In particular, we introduce a modified DL scheme that accounts for the intermittent and smooth deterministic model discussed in Sec. II.

To account for the assumption of smoothness, the signal estimation step described in Sec. III-A is modified by substituting the LASSO algorithm (4) with a smooth LASSO algorithm [39]. Smooth LASSO adds a penalty term for on/off switches within the transmitted signals. This penalty reflects prior knowledge that the transmitted signals do not switch on/off an excessive number of times. Accordingly, at each iteration $k$, given the channel iterate $\mathbf{H}^{(k)}$, instead of using (4), the proposed method obtains the updated estimate $\mathbf{X}^{(k)}$ of the channel matrix by solving the problem

[TABLE]

where $\mu$ is a weight parameter that is set depending the level of smoothness expected in the transmitted signals. Larger values of $\mu$ indicate fewer expected changes in the transmitted signal. Solving (7) directly has a complexity in the order of $O(MN^{2}T^{2})$ [40]. To reduce the computational complexity, we introduce two approximations of (7), namely Sequential Smooth LASSO (SL-SEQ) and Alternating Direction Method of Multipliers-based Smooth LASSO (SL-ADMM). SL-SEQ solves (7) column-wise sequentially, and can be found in [41] as an approximation for (7); while SL-ADMM is an ADMM-based [42, 43] iterative approach we propose to tackle (7) with lower computational complexity.

IV-A Sequential Smooth LASSO (SL-SEQ)

SL-SEQ solves the problem

[TABLE]

sequentially for $t=1,\dots,T$, where $\mathbf{x}(t-1)$ is the solution obtained at the previous step [41]. With this scheme, the quadratic penalty term imposes on the current solution a proximity constraint with respect to the previous one. SL-SEQ has a complexity that in the order of $O(MN^{2}T)$ as opposed to $O(MN^{2}T^{2})$ since it solves the $T$ problems in (8) sequentially.

IV-B Alternating Direction Method of Multipliers-Based Smooth LASSO (SL-ADMM)

We now propose an approximation to Smooth LASSO with ADMM approach [42, 43]. The novelty of SL-ADMM is to solve BSS problems with lower computational complexity than smooth LASSO for sources that have memory. To this end, we introduce a copy $\mathbf{x}^{\prime}(t)$ of $\mathbf{x}(t)$ and rewrite problem (7) as

[TABLE]

where we have defined the matrix $\mathbf{X}^{\prime}=[\mathbf{x}^{\prime}(1),\dots,\mathbf{x}^{\prime}(T-1)].$ Define the functions $l_{1}(\mathbf{x}(1))=||\mathbf{y}(1)-\mathbf{H}^{(k)}\mathbf{x}(1)||_{2}+\lambda||\mathbf{x}(1)||_{1}$ and $l_{t}(\mathbf{x}(t),\mathbf{x}^{\prime}(t-1))=||\mathbf{y}(t)-\mathbf{H}^{(k)}\mathbf{x}(t)||_{2}+\lambda||\mathbf{x}(t)||_{1}+\mu||\mathbf{x}(t)-\mathbf{x}^{\prime}(t-1)||^{2}$, for $t\geq 2$. Then, the augmented Lagrangian [44] for problem (9) can be written as:

[TABLE]

where $\boldsymbol{\alpha}(t)$ is the $N\times 1$ Lagrange multipliers for the constraints in (9) and $\rho\geq 0$ is a parameter. The proposed method tackles problem (9) via a primal-dual subgradient method that carries out the following steps at each iteration $i$ of ADMM

• •

For the current iterates $\boldsymbol{\alpha}^{(i)}(t)$, solve in parallel the $T$ problems:

[TABLE]

• $\cdot$

$\underset{\mathbf{x}(1)}{\text{minimize}}\text{ }l_{1}(\mathbf{x}(1))+\boldsymbol{\alpha}^{(i)T}(1)\mathbf{x}(1)$, (11a)

• $\cdot$

$\underset{\mathbf{x}(t),\mathbf{x}^{\prime}(t-1)}{\text{minimize}}\text{ }l_{t}(\mathbf{x}(t),\mathbf{x}^{\prime}(t-1))+\boldsymbol{\alpha}^{(i)T}(t)\mathbf{x}(t)-\boldsymbol{\alpha}^{(i)T}(t-1)\mathbf{x}^{\prime}(t-1)+\rho||\mathbf{x}(t)-\mathbf{x}^{\prime}(t-1)||^{2},\\ \text{ for }t=2,\dots,T-1$, (11b)

• $\cdot$

$\underset{\mathbf{x}(T),\mathbf{x}^{\prime}(T-1)}{\text{minimize}}\text{ }l_{T}(\mathbf{x}(T),\mathbf{x}^{\prime}(T-1))-\boldsymbol{\alpha}^{(i)T}(T-1)\mathbf{x}^{\prime}(T-1)+\rho||\mathbf{x}(T)-\mathbf{x}^{\prime}(T-1)||^{2}$, (11c)

• obtaining the new iterates $\mathbf{x}^{(i)}(t)$ for $t=1,\dots,T$ and $\mathbf{x}^{\prime(i)}(t)$ for $t=2,\dots,T$;
• •

Update the Lagrange multipliers as

[TABLE]

where $(\cdot)^{T}$ represents the transpose of a vector or a matrix.

We emphasize that SL-SEQ has a computational complexity in the order of $O(MN^{2}T)$, and SL-ADMM has a complexity in the order of $O(MN^{2}TI)$, where $I$ is the number of iterations in the ADMM approach. In general, $I$ is chosen to be a small value compared to $T$. It is seen that the computational complexity is greatly reduced with SL-SEQ and SL-ADMM compared to solving (7) directly.

V Per-Source Stochastic Filtering

In this section, we propose a further improvement to the DL-based source separation schemes discussed thus far which as shown in Fig. 3, consists of a per-source post-processing step. Specifically, once we have obtained an estimate $\tilde{x}_{n}(t),t=1,\dots,T$, of the signal of any source $n$ based on one of the DL-based methods discussed above, a perform per-source stochastic filtering (PSF) is applied that leverages the HMM of the sources to remove some of the ”noise” in the estimates $\tilde{x}_{n}(t)$. The ”noise” flips the state of the signal from active to inactive and vice-versa. Filtering of the noise leads to a more accurate activity map of the signals.

PSF is based on the hidden Markov model of Fig. 2. We first consider the case in which the transition probabilities $(p_{n},q_{n})$ of the hidden Markov model are known, and then study the scenario in which this information is not available.

V-A Known Model Parameters

Adopting a Bayesian formulation, a first approach would be to interpret the signal $\tilde{\mathbf{x}}_{n}(t)$ obtained from each source $n$ via DL-based source separation as the output of a memoryless channel $p(\tilde{\mathbf{x}}_{n}|\mathbf{x}_{n})$ whose input is the current transmitted signal $\mathbf{x}_{n}$. Composing this channel with the hidden Markov model of Fig. 2, one could apply the forward-backward message passing algorithm on the resulting Bayesian network in order to compute the posterior probabilities $p(x_{n}(t)|\tilde{x}_{n}(1),\dots,\tilde{x}_{n}(T))$ for all $t=1,\dots,T$ [45]. Here we do not follow this approach since it is generally unclear how to define the channel $p(\tilde{\mathbf{x}}_{n}|\mathbf{x}_{n})$, and the resulting message passing algorithm typically requires computationally expensive integrations. In contrast, we take a more pragmatic low-complexity approach, which is described next.

First, we perform a binary quantization of the available estimates $\tilde{x}_{n}(t),t=1,\dots,T$, by thresholding the absolute value of $\tilde{x}_{n}(t)$. This provides an estimate $\tilde{s}_{n}(t)$ of the state $s_{n}(t)$, since small absolute values of $\tilde{\mathbf{x}}_{n}(t)$ suggest the absence of signal for source $n$ at time $t$. As a result, we obtain the binary sequence $\tilde{s}_{n}(t),t=1,\dots,T$, where $\tilde{s}_{n}(t)=1$ if $|\tilde{x}_{n}(t)|>\gamma$ for a threshold $\gamma$, and $\tilde{s}_{n}(t)=0$ otherwise. The selection of the threshold $\gamma$ will allow us to obtain different points on the trade-off between the probability of false alarm and missed detection, as discussed in Sec. VI.

To formulate the filtering problem, we model the observations $\tilde{s}_{n}(t),t=1,\dots,T$, as being received at the output of a binary asymmetric channel (BAC) with flip probabilities $p_{n}^{\prime}$ and $q_{n}^{\prime}$, as illustrated in Fig. 4. The BAC models the errors made by the DL-based source separation algorithm that produced $\tilde{x}_{n}(t),t=1,\dots,T$, in detecting the state variables $s_{n}(t)$. Parameters $p_{n}^{\prime}$ and $q_{n}^{\prime}$ can be used as additional degrees of freedom in exploring the trade-off points between the probability of false alarm and missed detection. As an alternative, we will discuss in the next subsection how to estimate them from the data.

Given the observations $\tilde{s}_{n}(t),t=1,\dots,T$, and the parameters $p_{n}$, $q_{n}$, $p_{n}^{\prime}$ and $q_{n}^{\prime}$, we compute the posterior distribution $p(s_{n}(t)|\tilde{s}_{n}(1),\dots,\tilde{s}_{n}(T))$ by using the forward-backward algorithm [46]. Accordingly, the posterior probability can be written as

[TABLE]

where the probability $\alpha_{t}=\text{Pr}(s_{n}(t)=1|\tilde{s}_{n}(1),\dots,\tilde{s}_{n}(t))$ is obtained as a result of the forward pass, while the probability $\beta_{t}=\text{Pr}(\tilde{s}_{n}(t+1),\dots,\tilde{s}_{n}(T)|s_{n}(t)=1)$ is obtained from the backward pass as explained in [46, §17.4.3]. We then estimate $s_{n}(t)$ using the maximum a posteriori (MAP) approach, i.e., we set $\tilde{s}_{n}(t)=1$ if the inequality $p(s_{n}(t)|\tilde{s}_{n}(1),\dots,\tilde{s}_{n}(T))>0.5$ holds, and $\tilde{s}_{n}(t)=0$ otherwise. To recover an estimate $\hat{x}_{n}(t)$ of the transmitted signal $x_{n}(t),t=1,\dots,T$, we finally null all entries of $\tilde{x}_{n}(t)$ corresponding to states $s_{n}(t)$ that are estimated to be zero, i.e., $\hat{s}_{n}(t)=0$, while leaving unaltered the values of $\tilde{x}_{n}(t)$ for times $t$ at which $s_{n}(t)$ is estimated to be $\hat{s}_{n}(t)=1$. This can be summarized as

[TABLE]

V-B Unknown Model Parameters

We now study the case in which the parameters $p_{n}$ and $q_{n}$ of the HMM are unknown, and need to be estimated. We also jointly estimate the parameters $p_{n}^{\prime}$ and $q_{n}^{\prime}$ in the BAC of Fig. 4 that is used to derive the forward-backward filtering algorithm. To this end, we apply the Expectation-Maximization (EM) algorithm [46], as described next.

Given the estimated states $\tilde{s}_{n}(t)$ of each source $n$ obtained from the source separation and quantization steps, we would like to estimate the state sequence $s_{n}(t),\text{ }t=1,\dots,T$, the probabilities $p_{n}$ and $q_{n}$ in the transition matrix of the Markov chain, as well as the probabilities $p_{n}^{\prime}$ and $q_{n}^{\prime}$ in the BAC. The EM algorithm can be detailed as follows.

• •

Initialization: Initialize $p^{(0)}_{n}$, $q^{(0)}_{n}$, $p^{\prime(0)}_{n}$ and $q^{\prime(0)}_{n}$.

For each iteration $\nu=1,2,\dots$:

• •

E Step: Given the probabilities $p_{n}=p^{(\nu-1)}_{n}$, $q_{n}=q^{(\nu-1)}_{n}$, $p^{\prime}_{n}=p^{\prime(\nu-1)}_{n}$ and $q^{\prime}_{n}=q^{\prime(\nu-1)}_{n}$, apply the forward-backward algorithm to calculate the probabilities $\alpha_{t}$ and $\beta_{t}$ and hence the posterior probabilities (13), for $n=1,\dots,N$.

• •

M Step: Update the probability parameters by averaging the sufficient statistics with respect to the posterior distributions identified during the E step [47]. This leads to the updates

[TABLE]

[TABLE]

where $\tilde{\mathbf{s}}_{n}=[\tilde{s}_{n}(1),\dots,\tilde{s}_{n}(T))]$, $\text{Pr}(s_{n}(t-1)=0|\tilde{\mathbf{s}}_{n})=1-\text{Pr}(s_{n}(t-1)=1|\tilde{\mathbf{s}}_{n})$ is obtained from (13), and the posterior joint probabilities in (15) and (16) are computed as detailed below; and

[TABLE]

[TABLE]

where $\mathbbm{1}(a)$ is the indicator function, i.e. $\mathbbm{1}(a)=1$ if $a$ is true and 0 otherwise.

• •

Stopping Criterion: The iteration stops when the parameters converge, i.e. $|p^{(\nu)}_{n}-p^{(\nu-1)}_{n}|<\epsilon$, where $\epsilon$ is a small value.

The posterior joint probabilities in (15) and (16) are computed as [48]

[TABLE]

where $\alpha_{t,j}$ is the $j$-th element of vector $\boldsymbol{\alpha}_{t}=[1-\alpha_{t},\alpha_{t}]^{T}$; $\beta_{t,i}$ is the $i$-th element of vector $\boldsymbol{\beta}_{t}=[1-\beta_{t},\beta_{t}]^{T}$; $\Phi_{j,i}^{(\nu)}$ is the $(j,i)$-th element of the transition matrix $\boldsymbol{\Phi}^{(\nu)}_{n}$ in the $\nu$-th iteration; and matrix $\boldsymbol{\Phi}^{(\nu)}_{n}$ is defined as

[TABLE]

VI Numerical Results

In this section, we present numerical results to obtain insights into the performance of different DL-based source separation schemes and on the advantage of PSF. We consider the following source separation schemes, all implemented with and without the post-processing PSF step: OMP, LASSO, SL-SEQ and SL-ADMM. As performance criteria, we adopt the probability of false alarm $P_{fa}$ and probability of detection $P_{d}$. The detection probability $P_{d}$ is the ratio between the number of correctly detected active sources and the total number of active sources over the $T$ symbols, and the false alarm probability $P_{fa}$ is the ratio between the number of incorrectly detected active sources and the total number of inactive sources over the $T$ time samples. We also consider the performance of source estimation in terms of error vector magnitude (EVM) to synchronized digital transmission (see Sec. II).

Unless stated otherwise, the numerical results were obtained for $N=30$ sources, $M=20$ sensors, signal to noise ratio (SNR) per source per sample 30 dB, $T=1000$ time samples. We assume an HMM for the state $s_{n}(t)$, which is defined by transition probabilities $p_{n}=0.0022$, $q_{n}=0.02$, for all $N$ sources, so that an average $Np_{n}/(p_{n}+q_{n})=3$ sources are active at each time sample $t$, and the average duration of transmission is $1/q_{n}=50$ time samples. Also, unless stated otherwise, the algorithms with PSF are implemented with fixed values $p^{\prime}_{n}=0.02$ and $q^{\prime}_{n}=0.27$. We optimize numerically over the multipliers $\lambda$ and $\mu$ in order to satisfy given constants on the probability of detection $P_{d}$ or the probability of false alarm $P_{fa}$. For PSF with unknown model parameters, we initialize the transition probabilities of the Markov model as $p_{n}=0.5$, $q_{n}=0.5$, and the flip probabilities of the BAC as $p^{\prime}_{n}=0.1$, $q^{\prime}_{n}=0.2$. Throughout this section, the number of iterations for SL-ADMM was fixed at $K=30$; the Lagrange multipliers $\boldsymbol{\alpha}$ in (11) are initialized as the all-one vector; $\rho$ is set as 0.1 and $\gamma$ is fixed at 0.5.

VI-A Tuning Multipliers $\lambda$ and $\mu$

We start by discussing the choice of the multipliers $\lambda$ and $\mu$ in optimization problems (8) and (9) for SL-SEQ and SL-ADMM. These parameters penalize the number of non-zero elements in the solution vector and the number of state changes of each source, respectively. The choice of these two multipliers is scenario-dependent, and a discussion on their optimal selection can be found in [49]. In our extensive simulations, we found that the rule-of-thumb choices $\lambda=1/$SNR and $\mu=0.1/q_{n}$, where $q_{n}$ is the transition probability from an active to an inactive state (see Fig. 2), works well in practice. The selection $\lambda=1/$SNR reflects a decrease in the relevance of the sparsity prior as the quality of the observations increases [49]. In contrast, the multiplier $\mu\propto 1/q_{n}$ increases as the average transmission duration $1/q_{n}$ increases, implying fewer state changes.

As an example, Fig. 6 and Fig. 7 show the probability of detection $P_{d}$ for a fixed probability of false alarm $P_{fa}=0.1$ versus $\lambda$ and $\mu$, respectively. From Fig. 6, we observe that, when varying the SNR level, the optimal selection of $\lambda$ changes according to the mentioned inverse proportionality rule $1/$SNR. In a similar manner, Fig. 7 confirms the validity of our choice for multiplier $\mu$.

VI-B Source Activity Detection

We then investigate the trade-off between the probability of detection $P_{d}$ versus the probability of false alarm $P_{fa}$ in Fig. 5. The transmitted signal $x_{n}(t)$ is assumed here to be distributed as $x_{n}(t)\sim\mathcal{N}(0,1)$ whenever $s_{n}(t)=1$. We also assume known source statistics $p_{n}$ and $q_{n}$. A first observation is that both SL-SEQ and SL-ADMM significantly outperform all other source separation schemes, with the latter obtaining some performance gain over the former. Furthermore, PSF provides a significant performance boost for all algorithms. For example, at $P_{fa}=0.07$, the probability of detection with SL-ADMM is increased from 0.87 to 0.98.

We now explore the impact of different operating regimes on the probability of detection when the probability of false alarm is constrained to be smaller than 0.1. In particular, Fig. 8 shows performance with respect to different SNR values, and Fig. 9 investigates the impact of the sparsity level, which is defined as the ratio between the average number $Np_{n}/(p_{n}+q_{n})$ of active sources and the total number $N$ of sources. The former is modified by changing $p_{n}$. Both figures confirm the main conclusions obtained above in terms of the relative performance of the considered schemes. Moreover, Fig. 8 suggests that the performance undergoes a threshold phenomenon with respect to the SNR, particularly if implemented without PSF. It is also seen that SL with PSF is able to obtain a vanishing probability of missed detection as the SNR increases, unlike the other schemes whose probability of detection reaches a ceiling lower than one. Finally, Fig. 9 indicates that SL with PSF is robust to the sparsity level, while the other schemes are extremely sensitive to an increase in the average number of active sources.

VI-C Signal Estimation

As seen above, the SL schemes with PSF have the best performance in terms of source activity detection. Here, we study the performance in terms of the quality of signal estimation. To this end, we adopt the criterion of the EVM. EVM is defined as

[TABLE]

where $\mathcal{D}_{n}$ is the set of time samples in which user $n$ is correctly detected as active. Here we assume that, when $s_{n}(t)=1$, the transmitted signal is binary, that is $x_{n}(t)=$1 or $-1$ with equal probability.

Assuming again known sources’ transition probabilities, Fig. 10 shows the EVM for all the considered algorithms with respect to SNR where the probability of missed detection $P_{md}$ is constrained to be smaller than 0.2. The behavior of the EVM is in line with the discussion above regarding source activity detection. In particular, we observe a 40% decrease at 30 dB SNR that are achievable with SL and PSF, as well as the threshold behavior as a function of the SNR.

VI-D Parameter Estimation

Here, we evaluate the performance loss incurred when the parameters $p_{n}$ and $q_{n}$ in the Markov model defining the sources’ activities, are not known. Using the framework of Sec. V-B, We jointly estimate the parameters $p_{n}$, $q_{n}$, $p^{\prime}_{n}$ and $q^{\prime}_{n}$ with the EM algorithm. Fig. 11 investigates the performance of the estimation of the hidden parameters $p^{\prime}_{n}$ and $q^{\prime}_{n}$ in the BAC. From the figure, the estimate is seen to be close to the real value. The parameters $p_{n}$ and $q_{n}$ in the HMM have similar estimation accuracy as $p^{\prime}_{n}$ and $q^{\prime}_{n}$.

Fig. 12 shows the probability of detection $P_{d}$ versus the probability of false alarm $P_{fa}$ under the same conditions as Fig. 5. A first observation is that, even with unknown source parameters, PSF can provide a significant performance boost. Nevertheless, the gain is somewhat reduced as compared to the case with perfect knowledge. For instance, for SL-ADMM, when $P_{fa}=0.07$, the probability of detection $P_{d}$ is reduced from 0.98 to 0.89. The loss is generally increases when $T$ is smaller.

VII Conclusions

In this paper, we have introduced a two-stage Dictionary Learning (DL)-based algorithm for solving the BSS problem as the presence of radio sources with memory that are observed over slow flat-fading wireless channels. The DL stage of the proposed algorithm exploits source memory information to aid with the source separation. The effect of source memory is accounted for by a penalty term that discourages short-duration transmissions. The PSF stage utilizes source model information to learn about unknown source model parameters to further enhance source estimation. Numerical results show that the proposed algorithm outperforms existing DL algorithms in terms of source activity detection as well as source signal estimation. As a representative numerical example, even with unknown source statistics, it is shown that with 30 potential sources, 20 sensors and average 3 sources active in each time sample, it is possible to increase the probability of detection from 0.9 without PSF to 0.96 with PSF at 0.1 probability of false alarm and 30 dB SNR. The proposed algorithm can also decrease the EVM from by 40% at 0.8 probability of detection at 30 dB SNR. In addition, the proposed algorithm is robust to the source sparsity level, while existing DL algorithms are extremely sensitive to an increase in the number of active sources.