Structure-Based Subspace Method for Multi-Channel Blind System Identification

TL;DR

This paper introduces a novel subspace-based approach for blind identification of multichannel FIR systems that leverages the channel's Toeplitz structure, offering advantages over standard methods in challenging scenarios.

Contribution

It presents a new structure-based subspace (SSS) method that directly exploits the Toeplitz structure for improved blind system identification.

Findings

Outperforms standard subspace method in adverse conditions

Effectively estimates multichannel FIR systems

Can be extended to other linear structures like Hankel

Abstract

In this work, a novel subspace-based method for blind identification of multichannel finite impulse response (FIR) systems is presented. Here, we exploit directly the impeded Toeplitz channel structure in the signal linear model to build a quadratic form whose minimization leads to the desired channel estimation up to a scalar factor. This method can be extended to estimate any predefined linear structure, e.g. Hankel, that is usually encountered in linear systems. Simulation findings are provided to highlight the appealing advantages of the new structure-based subspace (SSS) method over the standard subspace (SS) method in certain adverse identification scenarii.

Figures

Tables

Table 1. TABLE I: Duality Table

Method	Toeplitz Structure	Orthogonality
SS	forced	minimized
SSS	minimized	forced

Equations

$$y(t)=\sum _{k}h(t-k)s(k)+e(t),t=0,\dots ,N-1$$

$$\mathbf{y}(t)=\sum _{i=0}^{L-1}\mathbf{h}(i)s(t-i)+\mathbf{e}(t)$$

$${\mathbf{y}}_M(t)={\mathcal{H}}_M(\mathbf{h}){\mathbf{s}}_{M+L-1}(t)+{\mathbf{e}}_M(t)$$

$$\displaystyle{\boldsymbol{\cal H}}_{M}({\mathbf{h}})=\left[{\begin{array}[]{*{20}{c}}{{\mathbf{h}}(0)}&\cdots&{{\mathbf{h}}(L-1)}&\cdots&0\\ \vdots&\ddots&\ddots&\ddots&\vdots\\ 0&\cdots&{{\mathbf{h}}(0)}&\cdots&{{\mathbf{h}}(L-1)}\end{array}}\right]$$

$${\mathbf{v}}^H{\mathcal{H}}_M(\mathbf{h})=0$$

$$\displaystyle{{\mathbf{h}}^{H}}\left[{\begin{array}[]{*{20}{c}}{{{\mathbf{v}}_{1}}}&\cdots&{{{\mathbf{v}}_{M}}}&{0}&{0}\\ {0}&\ddots&{}\hfil&\ddots&{0}\\ {0}&{0}&{{{\mathbf{v}}_{1}}}&\cdots&{{{\mathbf{v}}_{M}}}\end{array}}\right]={{\mathbf{h}}^{H}}{\boldsymbol{\cal V}}=0$$

$${\mathcal{V}}_i^H{\mathbf{h}}^{\prime }=0,\mathrm{for}i=1,\cdots ,d⟺{\mathbf{H}}^{\prime }(z)=\alpha \mathbf{H}(z)$$

$${\mathbf{R}}_y={\mathcal{H}}_M(\mathbf{h}){\mathbf{R}}_s{\mathcal{H}}_M^H(\mathbf{h})+{\sigma }_e^2\mathbf{I}$$

$${\mathbf{R}}_y={\mathbf{V}}_s\mathrm{diag}\left({\lambda }_1^2,\cdots ,{\lambda }_{\mathrm{M}+\mathrm{L}-1}^2\right){\mathbf{V}}_{\mathrm{s}}^{\mathrm{H}}+{\sigma }_{\mathrm{e}}^2{\mathbf{V}}_{\mathrm{e}}{\mathbf{V}}_{\mathrm{e}}^{\mathrm{H}}$$

$$\hat{\mathbf{h}}=\arg \min \Vert {\mathbf{V}}_e^H{\mathcal{H}}_M(\mathbf{h}){\Vert }^2=\arg \min {\mathbf{h}}^H\left[\sum _i{\mathcal{V}}_i{\mathcal{V}}_i^H\right]\mathbf{h}$$

$$\displaystyle\begin{array}[]{l}{\cal J}={{\cal J}_{1}}+{{\cal J}_{2}}+{{\cal J}_{3}}\\ =\;{\left|{\sum\limits_{j=1}^{K-1}{\sum\limits_{i=1}^{m(M-1)}{\hat{w}(i,j)-\hat{w}(i+m,j+1)}}}\right|^{2}}\\ +{\left|{\sum\limits_{j=L+1}^{K}{\hat{w}(1:m,j)}}\right|^{2}}+{\left|{\sum\limits_{i=m+1}^{mM}{\hat{w}(i,1)}}\right|^{2}}\end{array}$$

$${\mathcal{J}}_1={\Vert {\mathbf{I}}_L\hat{\mathbf{W}}{\mathbf{I}}_R-{\mathbf{J}}_L\hat{\mathbf{W}}{\mathbf{J}}_R\Vert }^2$$

$$\displaystyle\begin{array}[]{l}{{\cal J}_{1}}={\|{\left({{{\mathbf{I}}_{R}}\otimes{{\mathbf{I}}_{L}}-{\mathbf{J}}_{R}^{T}\otimes{{\mathbf{J}}_{L}}}\right){\rm{vec}}(\hat{\mathbf{W}})}\|^{2}}\\ \quad={\|{\left({{{\mathbf{I}}_{R}}\otimes{{\mathbf{I}}_{L}}-{\mathbf{J}}_{R}^{T}\otimes{{\mathbf{J}}_{L}}}\right)\left({{\mathbf{I}}\otimes{\mathbf{V}}_{s}}\right){\mathbf{q}}}\|^{2}}={\|{{{\mathbf{K}}_{1}}{\mathbf{q}}}\|^{2}}\end{array}$$

$$\displaystyle\begin{array}[]{l}{{\cal J}_{2}}={\|{\hat{\mathbf{W}}\left({1:{m},L+1:end}\right)}\|^{2}}=\|{{\mathbf{V}}_{s,row}}{\mathbf{Q}}{{\mathbf{I}}_{row}}\|^{2}\\ \quad\quad={\|{\left({{{\mathbf{I}}_{row}}\otimes{{\mathbf{V}}_{s,row}}}\right){\mathbf{q}}}\|^{2}}={\|{{{\mathbf{K}}_{2}}{\mathbf{q}}}\|^{2}}\end{array}$$

$$\displaystyle\begin{array}[]{l}{{\cal J}_{3}}={\|{\hat{\mathbf{W}}\left({m+1:mM,1}\right)}\|^{2}}=\|{{\mathbf{V}}_{s,col}}{\mathbf{Q}}{{\mathbf{I}}_{col}}\|^{2}\\ \quad\quad={\|{\left({{{\mathbf{I}}_{col}}\otimes{{\mathbf{V}}_{s,col}}}\right){\mathbf{q}}}\|^{2}}={\|{{{\mathbf{K}}_{3}}{\mathbf{q}}}\|^{2}}\end{array}$$

$$\underset{\mathbf{q}}{\min }{\mathbf{q}}^H{\mathbf{K}}^H\mathbf{K}\mathbf{q}$$

$$\displaystyle\begin{array}[]{l}{{\bf{h}}_{1}}={\left[{\begin{array}[]{*{20}{c}}1&{-2\cos(\theta)}&1\end{array}}\right]^{T}},\\ {{\bf{h}}_{2}}={\left[{\begin{array}[]{*{20}{c}}1&{-2\cos(\theta+\delta)}&1\end{array}}\right]^{T}}\end{array}$$

$$SNR(dB)=10{\log }_{10}E\frac{{\Vert {\mathcal{H}}_N{\mathbf{s}}_{N+L-1}\Vert }^2}{mN{\sigma }_v^2}$$

$$MSE(dB)=20lo{g}_{10}\left(\frac{1}{||\mathbf{h}||}\sqrt{\frac{1}{N_{mc}}\sum _{i=1}^{N_{mc}}||{\hat{\mathbf{h}}}_i-\mathbf{h}|{|}^2}\right)$$

Full text

I Introduction

Blind system identification (BSI) is one of the fundamental signal processing problems that was initiated more than three decades ago. BSI refers to the process of retrieving the channel’s impulse response based on the output sequence only. As it has so different applications, such as mobile communication, seismic exploration, image restoration and other medical applications, it has drawn researchers’ attention and resulted in a plethora of methods. Since then, a class of subspace-based methods dedicated to BSI has been developed, including the standard subspace method (SS) [1, 2], the cross-relation (CR) method [3] and the two-step maximum likelihood (TSML) method [4]. According to the comparative studies which have been done early in [5] and [6], the SS method is claimed to be the most powerful one.

In this paper, we introduce another subspace-based method based on the channel’s Toeplitz structure which is employed directly to formulate our cost function. The Toeplitz structure is an inherent nature that exists in most of the linear systems due to their convolutive nature.

The paper presentation focuses at first on the development of the proposed structure-subspace (SSS) method. Then, we highlight the improvement that is obtained by the SSS method over SS in the case of channels with closely spaced roots. The SSS method sounds to be a promising technique, yet it has a higher computational complexity that needs to be addressed in a future work.

Notation: The invertible column vector-matrix mappings are denoted by ${\rm{vec}}\{.\}:\mathbb{C}{{}^{a\times b}}\to\mathbb{C}{{}^{ab\times 1}}$ and ${\rm{mat}}_{a,b}\{.\}:\mathbb{C}{{}^{ab\times 1}}\to\mathbb{C}{{}^{a\times b}}$. $\left({{{\mathbf{A}}}\otimes{\mathbf{B}}}\right)$ is the Kronecker product. ${\mathbf{A}}^{T}$ and ${\mathbf{A}}^{H}$ denote the transpose and Hermitian transpose, respectively.

II Problem Formulation

II-A Multi-channel model

Multichannel framework is considered in this work. It is obtained either by oversampling the received signal or using an array of antennas or a combination of both [7]. To further develop the multi channel system model, consider the observed signal $y(t)$ from a linear modulation over a linear channel with additive noise given by

[TABLE]

where $h(t)$ is the FIR channel impulse response, $s(k)$ are the transmitted symbols and $e(t)$ is the additive noise. If the received signal is oversampled or recorded with $m$ sensors, the signal model in (1) becomes $m$-variate and expressed as

[TABLE]

where ${\mathbf{y}}(t)=[y_{1}(t),\cdots,y_{m}(t)]^{T}$, ${\mathbf{h}}(i)={[h_{1}(i),\cdots,h_{m}(i)]^{T}}$, ${\mathbf{e}}(t)={[e_{1}(t),\cdots,e_{m}(t)]^{T}}$. Define the system transfer function ${\mathbf{H}}(z)=\sum\nolimits_{k=0}^{L-1}{{\mathbf{h}}(k){z^{-k}}}$ with $(L-1)=deg({\mathbf{H}}(z))$. Consider the noise to be additive independent white circular noise with $E[{\mathbf{e}}(k){{\mathbf{e}}^{H}}(i)]=\delta_{k,i}\sigma_{e}^{2}{\mathbf{I}_{m}}$. Assume reception of a window of $M$ samples, by stacking the data into a vector/matrix representation, we get:

[TABLE]

where ${\mathbf{y}}_{M}(t)={[{\mathbf{y}}^{H}(t),\cdots,{\mathbf{y}}^{H}(t-M+1)]^{H}}$, ${\mathbf{s}}_{M+L-1}(t)={[{s}(t),\cdots,{s}(t-M-L+2)]^{T}}$, ${\mathbf{e}}_{M}(t)$ is stacked in a similar way to as ${\mathbf{y}}_{M}(t)$, and ${{\boldsymbol{\cal H}}_{M}({\mathbf{h}})}$ is an $mM\times(M+L-1)$ block Toeplitz matrix defined as

[TABLE]

${\mathbf{h}}$ is the desired parameter vector containing all channels taps, i.e. ${\mathbf{h}}=[{\mathbf{h}}(0)^{T},\cdots,{\mathbf{h}}(L-1)^{T}]^{T}$. Using the observation data in (3), our objective is to estimate the different channels’ impulse responses, i.e, recover ${\mathbf{h}}$ up to a possible scalar ambiguity. In the following subsection, we describe the subspace method, briefly.

II-B Subspace method revisited

For consistency and reader’s convenience, the SS method [1] which is also referred to as noise subspace method, shall be reviewed hereafter. The SS method implicitly exploits the Toeplitz structure of the filtering matrix ${\boldsymbol{\cal H}}_{M}({\mathbf{h}})$. Let ${\mathbf{v}}=[{\mathbf{v}}_{1}^{H},\cdots,{\mathbf{v}}_{M}^{H}]^{H}$, where ${\mathbf{v}}_{i}=[v_{(i-1)m+1},\cdots,v_{im}]^{T},\;i=1,\ldots,M$, be in the orthogonal complement space of the range space of ${\boldsymbol{\cal H}}_{M}({\mathbf{h}})$ such that

[TABLE]

Using the block Toeplitz structure of ${{\boldsymbol{\cal H}}}_{M}({\mathbf{h}})$, the above linear equation can be written in terms of the channel parameter ${\mathbf{h}}$ as

[TABLE]

The former equation can be used to estimate the channel vector ${\mathbf{h}}$ provided that (12) has a unique solution. Moulines et al. [1] proposed the SS method which is based on the following theorem:

Theorem 1**.**

Assume that the components of ${\mathbf{H}}(z)$ have no common zeros, and $M\geq L$. Let $\left\{{{\mathbf{v}}_{i}}\right\}_{i=1}^{d}$ be a basis of the orthogonal complement of the column space of ${\boldsymbol{\cal H}}_{M}({\mathbf{h}})$, then for any ${\mathbf{H}}^{\prime}(z)$ with $deg({\mathbf{H}}^{\prime}(z))=L-1$ we have

[TABLE]

where $\alpha$ is some scalar factor.

One of the encountered ways to estimate the orthogonal complement of ${\boldsymbol{\cal H}}_{M}({\mathbf{h}})$, i.e. noise subspace, is the signal-noise subspace decomposition. From the multi-channel model and noise properties, the received signal covariance matrix ${{\mathbf{R}}}_{y}=E[{\mathbf{y}}_{M}(t){\mathbf{y}}_{M}^{H}(t)]$ is given as

[TABLE]

The singular value decomposition of ${{\mathbf{R}}}_{y}$ has the form

[TABLE]

where $\lambda_{i}^{2}$, $i=1,\cdots,M+L-1$ are the principal eigenvalues of the covariance matrix ${\bf R}_{y}$. Also, the columns of ${\mathbf{V}}_{s}$ and ${\mathbf{V}}_{e}$ span the so-called signal and noise subspaces (orthogonal complement), respectively. After having the basis of the noise subspace, the channel identification can be performed based on the following quadratic optimization criterion:

[TABLE]

In brief, the SS method achieves the channel estimation by exploiting the subspace information (i.e. ideally, (${\rm Range}{{\boldsymbol{\cal H}}}_{M}({\mathbf{h}}))={\rm Range}({\mathbf{V}}_{s})~{}\bot~{}{\rm Range}({\mathbf{V}}_{e})$) as well as the block Sylvester (block-Toeplitz) structure of the channel matrix. More precisely, it enforces the latter matrix structure through the use of relations (8) and (12) and minimizes the subspace orthogonality error in (16). In the sequel, we propose a dual approach which enforces the subspace information (i.e. ${\rm Range}({{\boldsymbol{\cal H}}}_{M}({\mathbf{h}}))={\rm Range}({\mathbf{V}}_{s})$ where ${\mathbf{V}}_{s}$ refers to the principal subspace of the sample covariance matrix) while minimizing a cost function representing the deviation of ${{\boldsymbol{\cal H}}}_{M}({\mathbf{h}})$ from the Sylvester structure as indicated in Table I.

III Structure-Based SS method (SSS)

In the proposed subspace method, one searches for the system matrix ${{\boldsymbol{\cal H}}}_{M}$ in the form $\hat{{\boldsymbol{\cal H}}}_{M}={\mathbf{V}}_{s}{\mathbf{Q}}$ so that the orthogonality criterion in (16) is set equal to zero, i.e. $\|{\mathbf{V}}_{e}^{H}\hat{{\boldsymbol{\cal H}}}_{M}\|^{2}=0$ while ${\mathbf{Q}}$ is chosen in such a way the resulting matrix is close to the desired block Toepliz structure. This is done by minimizing w.r.t. ${\mathbf{Q}}$ the following structure-based cost function (informal Matlab notions are used):

[TABLE]

where $K=M+L-1$ and $\hat{{\mathbf{W}}}$ refers to $\hat{{\boldsymbol{\cal H}}}_{M}$. The cost function in (20) is inspired and matched to the Toeplitz structure introduced in (7). It is a composite of three parts; ${\cal J}_{1}$ seeks to force Toeplitz structure on the possibly non-zero entries, while ${\cal J}_{2}$ and ${\cal J}_{3}$ account for the zero entries in the first $m$ rows and first column, respectively.

Starting with ${\cal J}_{1}$, one can express it in a more compact way as follows:

[TABLE]

where:

${\mathbf{I}}_{L}$ is the $(mM)\times(mM)$ left identity square matrix with setting the last $m$ diagonal entries to zeros.

${\mathbf{I}}_{R}$ is the $K\times K$ right identity square matrix with setting the last diagonal entry to zero.

${\mathbf{J}}_{L}$ is a $(mM)\times(mM)$ square translation matrix with ones on the sub-diagonal and zeros elsewhere, i.e., ${\left[{{J_{L}}}\right]_{i,j}}={\delta_{i+{m},j}}$.

${\mathbf{J}}_{R}$ is a $K\times K$ square translation matrix with ones on the super-diagonal and zeros elsewhere, i.e., ${\left[{{J_{R}}}\right]_{i,j}}={\delta_{i,j+1}}$.

Now, using the Kronecker product property ${\rm{vec}}({\mathbf{A}}{\mathbf{G}}{\mathbf{B}})=\left({{{\mathbf{B}}^{T}}\otimes{\mathbf{A}}}\right){\rm{vec}}({\mathbf{G}})=\left({{{\mathbf{B}}^{T}}\otimes{\mathbf{A}}}\right){\mathbf{g}}$, one can write ${\cal J}_{1}$ as follows:

[TABLE]

where ${\mathbf{q}}={\rm{vec}}({\mathbf{Q}})$. In a similar way, ${\cal J}_{2}$ can be expressed as

[TABLE]

where ${\mathbf{V}}_{s,row}$ is the sub-matrix of ${\mathbf{V}}_{s}$ given by its first $m$ rows, and ${\mathbf{I}}_{row}$ is the $K\times K$ square identity matrix with setting the first $L$ diagonal entries to zero. Finally, ${\cal J}_{3}$ can also be set up as

[TABLE]

where ${\mathbf{V}}_{s,col}$ is the sub-matrix of ${\mathbf{V}}_{s}$ given by its last $m(M-1)$ rows, and ${\mathbf{I}}_{col}$ is the $K\times K$ square diagonal matrix with one at the first diagonal entry and zeros elsewhere.

As a result of (24), (27) and (30) the optimization problem in (20) is reduced to the minimization of the following quadratic equation

[TABLE]

where {\mathbf{K}}={\left[{\begin{array}[]{c|c|c}{{\mathbf{K}}_{1}^{T}}&{{\mathbf{K}}_{2}^{T}}&{{\mathbf{K}}_{3}^{T}}\end{array}}\right]^{T}}.

The optimal solution ${\mathbf{q}}$ of (31), under unit norm constraint of ${\mathbf{q}}$, is the least eigenvector that corresponds to the smallest eigenvalue of ${\mathbf{K}}^{H}{\mathbf{K}}$. The square matrix ${\mathbf{Q}}$ can be constructed by reshaping the obtained solution ${\mathbf{q}}$ from a vector into the matrix format, such that ${\mathbf{Q}}={\rm{mat}}_{K,K}\{{\mathbf{q}}\}$. Once matrix ${\mathbf{Q}}$ is obtained, the channel taps are estimated by averaging over the non-zero diagonal blocks of matrix ${\mathbf{V}}_{s}{\mathbf{Q}}$.

IV Discussion

In this section, we provide some insightful comments in order to highlight the advantages and drawbacks of the proposed subspace method.

• •

As explained earlier the proposed approach consists of neglecting the subspace error (i.e. considering Range(${\mathbf{V}}_{s}$) as perfect in the sense one searches for the desired solution within that subspace) while minimizing the system matrix (Toeplitz) structure error. The motivation behind this choice resides in the fact that the subspace error at the first order is null as shown in [8] and hence it can be neglected at the first order in favor of more flexibility for searching the appropriate channel matrix. This explains the observed gain of the SSS over SS method in certain difficult scenarii including the case of closely spaced channels roots.

• •

In the favorable cases where the channel matrix is well conditioned, the two subspace methods lead to similar performance as illustrated next by the simulation example of Fig. 1.

• •

For the SS method to apply one needs that the noise subspace vectors generate a minimal polynomial basis of the rational subspace orthogonal to ${\rm Range}({\mathbf{H}}(z))$ (see [1] for more details) and so the condition $M\geq L$ is considered to guarantee such requirement to hold. As the SSS does not explicitly rely on the orthogonality relation in (16), the latter condition might be relaxed as illustrated by the simulation example of Fig. 5.

• •

The proposed subspace method has a higher numerical cost as compared to the SS method. However, the cost might be reduced by taken into account the Kronecker products involved in building matrix ${\mathbf{K}}$. This issue is still under investigation together with an asymptotic statistical performance analysis of SSS.

• •

In the case $M\geq L$, the solution of (31) can be shown to be unique (up to a constant) thanks to the identifiability result of Theorem 1. Indeed, if ${\mathbf{q}}^{\prime}$ is another solution zeroing criterion (20), then the FIR filter associated to matrix ${\boldsymbol{\cal H}}^{\prime}={\mathbf{V}}_{s}{\mathbf{Q}}^{\prime}$ satisfies all conditions of Theorem 1, which leads to ${\mathbf{V}}_{s}{\mathbf{Q}}^{\prime}=\alpha{\mathbf{V}}_{s}{\mathbf{Q}}$ or equivalently ${\mathbf{Q}}^{\prime}=\alpha{\mathbf{Q}}$.

V SIMULATION RESULTS

In this section, the devised SSS method will be compared to the standard SS method as a benchmark. Three different experiments will be examined to illustrate the behavior of SSS in different contexts.

Two FIR channels are considered, each has a second order impulse response given by [6]:

[TABLE]

where $\theta$ is the absolute phase value of ${\bf{h}}_{1}$’s zeros and $\delta$ indicates the angular distance between the zeros of the two channels on the unit circle. Small $\delta$ results into an ill-conditioned system. In all simulations, the excitation signal is a 4-QAM, each channel receives $N=100$ samples, and the noise is additive white Gaussian. Note that the SNR is defined as

[TABLE]

The performance measure is the mean-square-error (MSE), given as

[TABLE]

where $N_{mc}=100$ refers to the number of Monte-Carlo runs and $\hat{{\mathbf{h}}}_{i}$ is the channel vector estimate at the $i$-th run.

In the first experiment given by Fig. 1, we show that for a well-conditioned system ($\delta=\pi$), both methods have a comparable performance. In the second one, we consider ill-conditioned systems (i.e. poor channel diversity). In that case, the devised SSS method outperforms the SS method at a moderate ill-conditioned system ($\delta=\pi/10$), and its performance gain becomes more obvious at severely ill-conditioned case ($\delta=\pi/50$) as shown in Fig.’s 2 and 3, respectively. When the system is ill-conditioned, the difference becomes more pronounced at low and moderate SNR values. Also, at the severe ill-conditioned case, the SS methods becomes unresponsive to the changes in the signal’s SNR, as revealed in Fig. 3, since the effect of ill-conditioning becomes prominent at low SNR. Figure 4 depicts the consequence of varying $\delta$ on the MSE for SNR=10dB.

In the last experiment, the number of channels is $m=3$ and the number of taps in each channel is $L=5$, the transfer function of the channels are given in [9]. In this experiment, we are primarily interested to look at the impact of the processing window length on the estimation performance. As can be seen from the results reported in Fig. 5, the performance of the SS method gets worse and degrades when the processing window length $M$ becomes less than the number of the channels’ taps $L$, while our proposed SSS is weakly affected by the window length condition, i.e. $M\geq L$. This allows us to reduce the dimension of the channel matrix ${{\boldsymbol{\cal H}}}_{M}$ with smaller window size values, especially for large dimensional systems where $m\gg 1$.

VI CONCLUSION

In this paper, we proposed a dual approach to the standard subspace method, whereby the channel matrix is forced to belong to the principal subspace of the data covariance matrix estimate while its deviation from Toeplitz structure is minimized. By doing so, we show that the channel estimation is significantly improved in the difficult context of weak channels diversity (i.e. channels with closely spaced roots). Interestingly, the principle of the proposed approach can be applied for estimation problems with other matrix structures where subspace method can be used.