Blind Identiﬁcation of Convolutive MIMO Systems with 3 Sources and 2 Sensors

We address the problem of blind identiﬁcation of a convolutive Multiple-Input Multiple-Output (MIMO) system with more inputs than outputs, and in particular, the 3-input 2-output case. We assume that the inputs are temporally white, non-Gaussian distributed, and spatially independent. Solutions for the scalar MIMO case, within scaling and permutation ambiguities, have been proposed in the past, based on the canonical decomposition of tensors constructed from higher-order cross-cumulants of the system output. In this paper, we look at the problem in the frequency domain, where, for each frequency we construct a number of tensors based on cross-polyspectra of the output. These tensors lead to the system frequency response within frequency dependent scaling and permutation ambiguities. We propose ways to resolve these ambiguities, and show that it is possible to obtain the system response within a scalar and a linear phase.


INTRODUCTION
The goal of blind r-input n-output (n × r) system identification is to identify an unknown system H(z), driven by r unobservable inputs, based on the n system outputs. Blind identification of a Multiple-Input Multiple-Output (MIMO) system is of great importance in many applications, such as speech enhancement in the presence of competing speakers, digital multiuser/multiaccess communications systems, biomedical engineering [1,2,3,4,5].
Most of the literature on n × r MIMO problems refers to the case of n ≥ r. In that case, system identification can lead to recovery of the inputs via deconvolution. Here we consider the case of more inputs than outputs, that is, n < r. In such a scenario, recovery of the input is generally not possible, except in cases where some a priori information about the inputs is available, such as the finite alphabet property [6,7]. Very few results exist for the convolutive MIMO problem with more inputs than outputs. In [8], the identifiability of a Moving Average (MA) system with possibly more inputs than outputs has been studied. A special case of a blind 2 × 3 convolutive system, where the cross-channels are simple delay elements, has been studied in [9]. The delays were estimated via a polyspectra based method.
The scalar 2 × 3 MIMO case has been approached in [6,10], based on the canonical decomposition of tensors, which were constructed from higher-order cross-cumulants of the system output. That approach yields the system within scaling and permutation ambiguities.
In this paper, we address the blind identification of 2 × 3 convolutive systems. We look at the problem in the frequency domain, where, for each frequency we construct a number of tensors based on cross-polyspectra of the system output. Based on these tensors, the problem at each frequency can be formulated as that of scalar MIMO, thus the method of [10] can be applied to yield the system frequency response within scalar and column permutation ambiguities. The latter ambiguities, however, now vary between different frequencies. We exploit the information that is available in the polyspectra domain to resolve these ambiguities, and recover the MIMO system frequency response within a scalar and linear phase ambiguities.

EXISTING RESULTS ON THE SCALAR MIMO CASE
In [6,10] a solution for the 3-input 2-output scalar MIMO problem has been proposed within scaling and permutation ambiguities. The following system was considered: T is 2 × 1 vector representing the output signals; H is the instantaneous mixing matrix of dimensions 2×3; s = s 1 (t) s 2 (t) s 3 (t) T is 3×1 vector representing the input signals and t denotes discerte time. The signals in (1) satisfy the following assumptions: (A1) s i (t), i = 1, 2, 3, are zero-mean, stationary, white, non-Gaussian distributed and independent of each other, with finite cumulants up to fourth-order; (A2) the s i (·)'s have unit variance, for i = 1, 2, 3. Define where CUM{·} represents the cross cumulant of signals, that is, Let Ꮿ 40 , Ꮿ 31 , and Ꮿ 30 be tensors with elements Let H p denote the pth column of the mixing matrix H. It holds that where and "•" denotes the tensor outer product [10]. Equations (4), (5), and (6) are referred to as "canonical decomposition" of Ꮿ 40 , Ꮿ 31 , and Ꮿ 30 . They are decompositions in a minimal sum of rank-1 terms [10]. Thus, the problem of estimating the matrix H is that of decomposing Ꮿ 40 , Ꮿ 31 , and Ꮿ 30 in a minimal number of rank-1 terms, in which each rank-1 term is the contribution of one source signal. For higher-order tensors, as opposed to matrices (second order arrays), the number of rank-1 terms is not bounded by the dimension of the column nor the row space. This enables the identification of systems with more inputs than outputs.
In the sequel, we outline the process of obtaining the decomposition of the above tensors [11]. Consider the cumulant based matrix equation where G is a 4×1 vector. This system of linear equations has to be solved in a least-square sense for a nontrivial G * . Assuming that G * has unit norm, the solution can be found as the right singular vector of the coefficient matrix in (8), corresponding to the smallest singular value. We view the elements of G as the coefficients of a polynomial, and let r 1 , r 2 , r 3 , be the three roots of this polynomial. Also, we defineH asH r 1 r 2 r 3 1 1 1 .
It was shown thatH is an estimate of the mixing matrix H up to an unknown column scaling and permutation, that is, where Λ is a diagonal matrix representing the real and positive column scaling, e jΦ is also a diagonal matrix representing the phase ambiguity of each column, while P is a permutation matrix representing the column permutation. These ambiguities are acceptable for the instantaneous mixture case.

THE PROPOSED APPROACH FOR THE 2 × 3 CONVOLUTIVE MIMO CASE
Consider the case of convolutive mixtures. The MIMO system output is given by where s(t) is the 3 × 1 input vector containing mutually independent entries with unit variance and non-Gaussian distribution, h(l) is the 2 × 3 impulse response matrix with elements {h i j (l)}, i = 1, 2, j = 1, 2, 3, x(t) is the 2 × 1 vector of observations, and t denotes discrete time.
In addition to assumptions (A1) and (A2) introduced for the scalar MIMO case, we further assume the following: (A3) there exist a nonempty subset of ω's, denoted by ω * , and a nonempty subset of the indices 1, . . . , n, denoted by l * , so that for l ∈ l * and ω ∈ ω * , the lth row of the matrix H(ω) has elements with magnitudes that are mutually different.
At first look, an extension of the scalar MIMO case to the convolutive case would appear feasible by observing that, in the frequency domain, it holds that where x(ω), s(ω), and H(ω) are the Discrete-Time Fourier transform of x(t), s(t), and h(l), respectively. Thus, (12) is similar to (1) for a fixed ω.
To extend the idea of [10] to this case, one would need to estimate the cross cumulants of x i (ω 1 ), x j (ω 2 ), . . . . Although cumulants of Fourier transforms has been considered in [12], there is a problem right there. According to Brillinger [13, Theorem 4.3.2 on page 93], where K is an integer. Thus, C 40 and C 30 would be identically zero. When considering discrete Fourier transforms the cumulants will not be zero, but they will be very small and thus sensitive to estimation errors.
In the following, we propose an approach that does not involve cumulants of x(ω). We define a set of tensors based on cross-polyspectra of the system output, which lead to expressions similar to those of (4) and (5).
First, we define three types of the fourth-order cross cumulants of the received signals [14]: where are the three types of the fourth-order cumulants of s p (·), respectively. The corresponding fourth-order crosspolyspectra, defined as the Fourier transform of c 40 , equals [14]: Taking (16) and (17), we get, respectively, These two equations enable us to construct two tensors -1 31 , with elements C 40 jkli (ω, ω, ω), and -2 31 , with elements C 31 jkli (ω, ω, ω), where j, k, l, i = 1, 2. We can show that these two tensors correspond to the tensor C 31 in the instantaneous case, that is, where H p (ω) denotes the pth column of H(ω). For system with real impulse response h(n), the two tensors are identical . On the other hand, for complex systems, the two tensors are nonidentical. This observation implies that we need to treat the real and complex cases separately.
The complex system case From the two tensors -1 31 (ω) and -2 31 (ω) we can derive four independent equations, to be used in the estimation of the polynomial G as defined in (8): Note for the convolutive MIMO system, the G is a function of frequency ω. The above equation provides enough equations for the unique estimation of the polynomial G(ω).

The real system case
This leaves one complex degree of freedom in the solution.
There is no easy way to get another equation like in -1 31 (ω), nor to get an equivalent of tensor Ꮿ 40 as in the instantaneous case. However, here we propose a tensor -22 , which can provide enough information to reach the solution. Taking ω 1 = ω 2 = −ω 3 = ω in (18), we get (23) For the real system case, we have This enables us to construct a tensor -22 based on the crosspolyspectra C 22 i jkl (ω, ω, −ω). Then, for the tensor -22 it holds that The tensor -22 can be utilized following the approaches proposed in [6,15]. The idea in [6] is to start an exhaustive search, such that the structure of -22 is also taken into account. For each possible solution of (22), the corresponding values of {H p (ω)} (1≤p≤3) and of the rank-1 tensors are computed. At that point, the "goodness" of the approximation by a sum of rank-1 tensors in (25) is assessed, and the global optimum is sought. Each step amounts to computing the roots of a polynomial of degree 3 (computation of {H p (ω)} (1≤p≤3) ) and verifying how close a vector in a real 35-dimensional space is to a subspace spanned by three other vectors (checking (25)). In [15], an alternating least squares (ALS) method was proposed to simultaneously solve the equations based on the two tensors; here one alternates between the computation of the roots of a polynomial of degree 3 and 2, respectively. Note here that the proposed method fails to estimate the system transfer function H(ω) at frequency ω = 0 and ω = π, because the tensors -1 31 and -22 are identical at these two frequencies. Since discrete frequencies will be used in the implementation, one possible remedy is to obtain the system transfer function H(ω) at these two frequencies by interpolation using the estimate in surrounding frequencies. Simulations examples show that the interpolation method works well with Finite Input Response (FIR) systems.
Based on the above methodology, we can get the estimate of H(ω), that is,H(ω), up to some ambiguities as stated in (10). These ambiguities are acceptable for the instantaneous mixture case, but not for the convolutive mixture case. As in the latter case the ambiguities are frequency dependent, one cannot combine the estimates of H(ω) at different frequencies to get a final estimation. We next propose some steps to solve this problem.

DEALING WITH FREQUENCY DEPENDENT AMBIGUITIES
The solutionH(ω) is related to H(ω) as

Estimation of Λ(ω)
Based on assumption (A2), it holds that where P X (ω) is the cross power spectrum matrix of the received signal x(t), thus can be estimated;H p (ω) is the pth column ofH(ω); and λ p (ω) is the pth diagonal element of Λ(ω).

Resolving frequency dependent permutation ambiguity P(ω)
Resolving frequency dependent column permutation ambiguity, P(ω), amounts to reducing it to a constant permutation matrix P. This can be achieved as follows.
Based on the definition of cross-polyspectra in (18), we construct a polyspectra matrix C 22 l1l2 (ω 1 , ω 2 , ω 3 ) whose (i, j)th element equals C 22 l1i jl2 (ω 1 , ω 2 , ω 3 ), then we can rewrite (18) in matrix form as follows: Based on assumption (A3), it holds that where P(p) represents the unknown column permutation. By solving (32), we get an estimate of the γ 22 sp |H lp (ω 3 )| 2 in some permuted order, where each γ 22 sp |H lp (ω 3 )| 2 , p = 1, 2, 3, is associated with a column vectorH P(p) (ω). By sorting the columns ofH p (ω) according to a predefined order of the estimated γ 22 sp |H lp (ω 3 )| 2 , which is the same for all ω's, we can achieve a constant permutation. Up to this point, we have the estimateH(ω) of the system transfer function with only phase ambiguity and constant permutation, that is,

H(ω)e jΦ(ω) P = H(ω).
(33) Resolving phase ambiguity e jΦ(ω) Equation (33) would suffice for the identification of minimum phase systems only. For nonminimum phase systems the phase ambiguity can be resolved as follows. The recovery of Φ(ω) is also based on the special structure of (29). Combining (29) and (33), we get Define which is a diagonal matrix. Since C 22 l1l2 (ω, −ω − α, ω 3 ) can be estimated, andH(ω) andH H (ω + α) are known at this stage, e jΨ(ω) can be estimated based on the following equation: in the same way as we did for solving the frequency dependent scaling and permutation. Based on (35), we can get a recursive equation of the unknown phase ambiguity Φ(ω) as follows: where Θ is Note that for fixed l 1 , l 2 , α, and ω 3 , Λ 2 (ω, −ω − α, ω 3 ) is independent of ω. Equation (37) can be solved as in [16,17], to obtain Φ(ω) up to a linear phase and constant phase ambiguity, which, respectively, correspond to a time delay and a complex scaling of the columns of the system impulse response matrix. Both of the latter ambiguities are acceptable for blind system identification. Finally, we obtain a solutionĤ(ω) up to a constant permutation, complex scaling and linear phase ambiguity, that is,Ĥ where M is a 3 × 3 diagonal matrix with integer elements that represents the linear phase ambiguity, Φ(0) represents the remaining constant phase ambiguity or the complex scaling. As a summary, the computation ofĤ(ω) is carried out in the following steps: (S1) estimate the cross-power spectrum matrix P X (ω) of the received signal vector x(t); (S2) estimate the fourth-order cross-polyspectra slices C 40 i jkl (ω, ω, ω), C 31 i jkl (ω, ω, ω), and C 22 i jkl (ω, ω, −ω) of the received signal vector x(t) using the indirect class method [14]. Then construct the matrix equation (21) for complex system. For real system, construct matrix equation (22) and use the method proposed in [6] or [15] to estimate the roots of the polynomials constructed by the columns of H(ω) at each discrete frequency ω. At this stage we get the estimateH(ω); (S3) solve for the frequency dependent scaling based on the estimated cross-power spectrum matrix P X (ω) by solving (27). This step yieldsH(ω); (S4) solve for the frequency dependent permutation ambiguity using the cross-polyspectra matrix C 22 ll (ω, −ω, ω 3 ) in (32), to obtainH(ω); (S5) compute the phase ambiguity using (37) (see [16,17] for more details). At this stage, we have the estimatê H(ω) of the system transfer function up to a constant permutation, complex scaling, and linear phase ambiguity; (S6) estimate the time domain impulse response H(l) using inverse Fourier transform ofĤ(ω).

SIMULATION RESULTS
In this section, we provide two simulation examples to demonstrate the feasibility of the proposed algorithm. Matlab code for the simulation example can be found at http://www.ece.drexel.edu/CSPL/.
To demonstrate the feasibility of the approach, we first provide an example based on true cumulants, rather than cumulant estimates and true power spectrum matrix. Example 1. The impulse response matrix H(l) is taken to be a 2 × 3 nonminimum phase system with transfer function The length of the DFT used in the computation of the cross-power spectrum and fourth-order cross-polyspectra was taken to be F = 128. Figure 1 shows the absolute value of the estimated roots of the polynomials constructed by the columns of H(ω) at each discrete frequency. As it can be seen in Figure 1, the estimation is good at most frequencies, except at frequencies 0 and π, where errors occur in the estimation of roots. These errors occur due to the failure of the proposed method as mentioned in Section 4. Figure 2 shows the computed frequency domain magnitude estimation after the frequency dependent scaling and permutation ambiguity have been recovered by the proposed approach. We can see that there are errors at certain frequencies due to errors in the estimation of the roots. Simulations show that after the interpolation of the system estimate at the frequencies 0 and π, the magnitude estimation will be almost perfect. Figure 3 shows the phase estimation result using the proposed approach. Since the phase is estimated by a recursive equation (see (37)), the errors at frequencies 0 and π tend to accumulate. However, simulations showed that the phase estimate can be improved greatly if the interpolated version of the system estimate is used. The time-domain impulse response h(l) was estimated by taking the inverse Fourier transform of the estimated channel frequency domain responseĤ(ω). Since we do not know the channel order L in advance, we can use an overestimated channel order L e to truncate the estimated impulse response using the inverse Fourier transform. In this experiment, the extended channel order was taken to be L e = 11. For comparison purpose, proper alignment and scaling with the true impulse response was performed. Figure 4 shows the estimated channel response estimation.
Numerical simulations also showed that the proposed methods work well with complex MIMO systems with 3input 2-output given true cumulants. When it comes to applying the same approach to cumulant estimates, the result is rather sensitive to estimation errors. Errors are mainly caused by the rooting step in obtaining an initial estimateH(ω). The reason is that the algorithm in [6] which we used to get the initial estimate, is rather sensitive to cumulant estimation errors and has local minima. Unfortunately, no better method exists at this time.
In the following, an example using cumulant estimates is given.
The inputs {s j (k)}, j = 1, 2, 3, were mutually independent, zero-mean i.i.d. signals, single-side exponentially distributed, with length T = 8192. The cross power spectrum matrixP x (ω) was estimated using the Blackman-Tukey method [18]. The polyspectra slices used in the algorithm were estimated via the indirect class method [14], and the sample cross-cumulant sequence was windowed by the Kaiser window with parameter 6 [14]. The DFT length in the computation of the cross-power spectrum and fourth-order cross-polyspectra was taken to be F = 64.
The extended channel order was taken to be L e = 5. Proper alignment and scaling with the true impulse response was also performed for comparison purpose. Figure 5 shows the estimated channel response estimation.

CONCLUSION
We proposed a polyspectra based frequency domain method to show the feasibility of MIMO system identification with more inputs than outputs. The method proposed in [6] and [10] for the instantaneous case was extended to the convolutive case, and the frequency dependent ambiguities related with frequency domain method were resolved using power spectrum and polyspectra matrices. The method was shown to work well when true cumulant were provided, while it was in general sensitive when cumulant estimates were used. No comparisons were provided because no other methods exist for the same problem.