Polynomial GSVD Beamforming for Two-User Frequency-Selective MIMO Channels

In this paper, we propose a generalized singular value decomposition (GSVD) for polynomial matrices, or polynomial GSVD (PGSVD). We then consider PGSVD-based beamforming for two-user, frequency-selective, multiple-input multiple-output (MIMO) multicasting. The PGSVD can jointly factorize two frequency-selective MIMO channels, producing a set of virtual channels (VCs), split into: private channels (PCs) and common channels (CCs). An important advantage of the proposed PGSVD-based beamformer, over the application of GSVD independently to each frequency bin of the orthogonal frequency division multiplexing (OFDM) scheme, is that it can facilitate different modulation and/or access schemes to various users. Using computer simulations, we characterize the bit error rate performance of our two-user MIMO multicasting system for different PCs/CCs configurations. Here, we also propose an OFDM-GSVD benchmark system, and show that our PGSVD-based beamformer compares favorably to this benchmark under erroneous and uncertain MIMO channel conditions, in addition to its advantage of facilitating heterogeneous modulation and access for various users.


Polynomial GSVD Beamforming for Two-User
Frequency-Selective MIMO Channels

I. INTRODUCTION
T HE transmission of data from one point to multiple points has become important in recent times particularly because of the increasing demand by end-users for wireless multimedia content, or multicasting. This is achieved by realizing physical-layer multicasting via adaptive beamforming. The beamformer is used to generate point-to-multipoint logical multiple-input multiple-output (MIMO) channels over the physical layer [1]- [3]. In [3], Senaratne and Tellambura propose a Manuscript received August 12, 2020; revised December 10, 2020 and January 6, 2021; accepted January 10, 2021. Date of publication January 18, 2021; date of current version February 5, 2021. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Lin Bai beamformer based on the generalized singular value decomposition (GSVD) [4] for data transmission from a single source to two users through two MIMO channels. The GSVD-based beamformer meets the needs of the end-user for multicasting by realizing virtual channels (VCs) over the physical wireless channel. The beamformer not only implements classic point-to-point VCs, called private channels (PCs), serving individual users, but also yields point-to-two point VCs, namely common channels (CCs), that address the needs of both users simultaneously. The GSVD is an extension of the singular value decomposition (SVD) [5] to two matrices. It performs joint factorization of two rectangular matrices with the constraint that the number of columns of the two matrices are equal. A number of algorithms for GSVD computation have been proposed, particularly Van Loan's approach in [6] which takes advantage of the link between GSVD and cosine-sine decomposition (CSD) [5], first established in [7]. Since each elementary decomposition is numerically stable, with mostly orthogonal transformations, this CSD-based approach is robust. A more general definition of the GSVD can be found in [8].
The GSVD-based beamformer mitigates co-channel interference (CCI), or crosstalk, between the two users and/or between the different VCs. The CCI is due to non-zero, off-diagonal elements of the two MIMO-channel matrices. The beamformer works by jointly minimizing these off-diagonal entries, and hence decoupling the MIMO channels. However, emerging wireless systems are designed to provide high data rates, which gives rise to intersymbol interference due to frequency-selective multipath fading [9], [10]. The GSVD beamformer would not be expected to combat multipath since it does not have the necessary degrees of freedom (DoFs) to realize the required temporal filters. To date, the problem of point-to-two point multicasting for frequency-selective MIMO channels has not been addressed in the literature.
One possible way forward in this regard is by way of a frequency-domain scheme based on orthogonal frequencydivision multiplexing (OFDM) [9], [11], [12]. Such an approach would utilize the discrete Fourier transform (DFT) over the entire data to convert the multipath-channel problem into a number of individual narrowband tones [13]. As such, the frequencyselective channel is converted into individual frequency-flat MIMO channels. The GSVD can then be applied to factorize the two (scalar) channel matrices at each bin. This GSVD-based MIMO OFDM (OFDM-GSVD) approach has not been investigated so far in the literature.
There are two possible drawbacks with a per tone GSVD method. Firstly, this will force all private channel users to employ an OFDM-based access scheme and applications, and OFDM frame-level synchronization. The second drawback of an OFDM-GSVD approach is that it would not be efficient in situations where the scalar MIMO-channel matrices are rank deficient in one or more frequency bins; so, for those bins, we have the underdetermined MIMO-system case [14]- [16]. Channelmatrix rank deficiency models spatial correlation and/or interbin interference. Under such conditions, the problem bin(s) would need to be identified and related-information relayed to the receiver. However, conveying this information itself may consume considerable bandwidth.
An alternative approach is one that generalizes the GSVDbased beamformer in the time-domain by utilizing iterative polynomial matrix factorization techniques. A frequency-selective MIMO channel can be represented by a matrix of filters, or a matrix with polynomial entries, where each entry now describes both the CCI and multipath in the MIMO system. Akin to the factorization of the scalar MIMO-channel matrix, a decomposition of polynomial matrices is required for successful decoupling of the two MIMO channels. Pioneering work by McWhirter has lead to an extension of the eigenvalue decomposition (EVD) to the polynomials, namely polynomial EVD (PEVD), in [17]- [19]. The PEVD decomposes a space-time covariance (parahermitian matrix) matrix into a paraunitary matrix and a diagonal parahermitian matrix.
Over the past decade, a number of time-domain, iterative algorithms for PEVD approximation have been proposed [17], [18]. These algorithms broadly fall under one of two categories: the second-order sequential best rotation (SBR2) [17] and the sequential matrix diagonalisation (SMD) algorithms [18]. Recently, their utility has been demonstrated through application to various problems, including blind source separation [20]- [22] and channel coding [23]- [26].
Analogous to the scalar GSVD, a GSVD for polynomial matrices (or PGSVD) would be a natural extension of the polynomial SVD (PSVD) [27], [28] for two polynomial matrices. Despite this innate connection and the many potential applications, work on PGSVD has not been reported on in the literature as yet. Therefore, in this paper, we first introduce the concept of a PGSVD and provide an efficient algorithm for its implementation. We then propose a broadband beamformer based on PGSVD to enable multicasting over frequency-selective MIMO channels. Here, the term broadband (or wideband) beamforming refers to the case where the (space-time) beamforming is performed for frequency-selective channels; so the beamformer characteristics are different for different frequencies, due to the frequency selectiveness of the channels.
The work considered in this paper, as well as in [3], requires the number of antennas to be of the same order as that of the intended spatial channels. However, when the number of antennas is significantly greater than the intended number of spatial channels, it may be beneficial to consider a hybrid analogue and digital beamforming, as in [29], [30], to reduce the complexity associated with the massive number of RF chains.

A. Contributions and Main Features
The main contribution of this paper is threefold: r An iterative, time-domain algorithm is proposed for computing an approximate GSVD of two rectangular polynomial matrices. The PGSVD is a novel concept which is proposed here for the first time in the literature as a polynomial-matrix equivalent of the scalar cosine-sine decomposition (CSD)-based approach in [6]; and has been achieved through the introduction of newly developed numerical algorithms.
r A broadband beamformer that enables multicasting over two frequency-selective multi-input multi-output (MIMO) channels is proposed, as illustrated in Fig. 1, which exploits PGSVD. Our approach is a non-trivial generalization of the narrowband scheme in [3].
r The proposed approach can facilitate heterogeneous modulation and/or multiple access schemes in different virtual channels (VCs), depending on user needs, which cannot be achieved with the benchmark system. The last salient feature above permits, for example, the use of code/frequency division multiple access (CDMA/FDMA) on the private channels (PCs), for secure point-to-point communications, whilst simultaneously maintaining high throughput links on the common channels (CCs) using, e.g., high-order modulation schemes, such as 256-QAM (quadrature amplitude modulation).

B. Notation and Overview
We use capital and small boldface variables, e.g. A and v, to denote matrices and vectors, respectively. The notation a(z) represents a function a : C → C; e.g., A(z) : C → C M ×N represents an M × N polynomial matrix in the indeterminate variable z −1 , and is given by  Table I. Below, in Section II, we review the scalar CSD and GSVD; in Section III, we define the PGSVD and propose an algorithm for its implementation based on an extension of CSD to polynomial matrices. Its utilization in realizing suitable beamformers for point-to-two point, frequency-selective, MIMO channels is discussed in Section IV. Our PGSVD algorithm and broadband MIMO beamformer are then evaluated in Section V. Finally, conclusions are provided in Section VI.

II. SCALAR MATRIX FACTORIZATIONS
Motivated by [6], here we describe an efficient and stable algorithm for computing the scalar GSVD by way of the CSD.

A. Generalized Singular-Value Decomposition
Consider the two matrices H 1 ∈ C M ×N and H 2 ∈ C P ×N . Let H = (H T 1 , H T 2 ) T and K = rank{H}. Define the GSVD: where U ∈ C M ×M and V ∈ C P ×P are unitary matrices and X ∈ C N ×N is a non-singular matrix; C ∈ C M ×N and S ∈ C P ×N are block diagonal matrices satisfying Depending on the relationship between M, P , and N , the matrices C and S will comprise diagonal sub-matrices Σ 1 and Σ 2 and zero matrices. These diagonal sub-matrices contain non-zero singular values, c i and s i , respectively, such that Σ 2 1 + Σ 2 2 = I (M +P −N ) . Hence, the element pairs satisfy c 2 i + s 2 i = 1, i = 1, . . . , K, and |c i | = 1, i = (K + 1), . . . , N; and are ordered: The diagonal elements of C and S define the generalized singular-value pairs (c i , s i ), or the generalized singular-values: it follows that the generalized eigenvalues are:

B. Cosine-Sine Decomposition
The CSD relates the SVD of Q 1 to that of Q 2 via the decompositions: where U ∈ C M ×M , V ∈ C P ×P and Z ∈ C N ×N are unitary matrices. The diagonal matrices C ∈ C M ×N and S ∈ C P ×N contain the non-negative sine and cosine pairs and satisfy the same properties as in GSVD. Indeed, GSVD is identical to CSD in the case where, given H in Section II-A, H H H = I N . Hence, CSD is a special case of GSVD.

C. Computing the GSVD
Highly based on the algorithm in [6], a common approach to computing the GSVD is one that uses the QRD and the CSD. The first stage of this method calculates the QRD of where Q is unitary and R is upper-triangular and non-singular. This approach proceeds by computing the CSD of Q, as in Section II-B, producing the GSVD unitary matrices U and V; and the diagonal matrices C and S. Given Z from (4), the GSVD non-singular matrix X is: For a more stable algorithm, the generalized singular-values σ i are categorized as small, A proof of this can be found in [6]. Once the diagonals of C are ordered as in (3), this categorization can be implemented by finding the index k such that There are a number of methods for computing the CSD. An efficient way involves calculating the SVD of only Q 1 from (4), yielding matrices U, C and Z; then the matrices V and S are obtained via the QRD [6], [31], [32]. In the next section, we generalize this approach to polynomial matrices.

III. POLYNOMIAL GENERALIZED SINGULAR-VALUE DECOMPOSITION
In this section, we extend the concepts in Section II and derive a stable algorithm for the computation of a PGSVD of two rectangular polynomial matrices. We first generalize the CSD algorithm to the realm of polynomials by way of a new PSVD algorithm. The proposed PSVD algorithm is based on SMD in [18] and utilizes the cost function in [23]. We then show how our new PSVD can be used to calculate the PGSVD.

A. Polynomial GSVD
Here we extend the definition of the scalar GSVD in (1) to the realm of polynomial matrices. Let H 1 (z) : C → C M ×N and H 2 (z) : C → C P ×N , and H(z) = (H T 1 (z), H T 2 (z)) T : C → C (M +P )×N . Then, assuming K = rank{H(z)} ≤ N , we define the PGSVD as: where U (z) : C → C M ×M and V (z) : C → C P ×P are paraunitary (PU) [19]; the polynomial matrix X(z) is non-singular. As in the scalar case (Section II), the structure of diagonal polynomial matrices C(z) and S(z) in (7) depends on the the relationship between M, P , and N . In this paper, we consider the more practical condition where M, P ≤ N ; in particular, if M, P = N , C(z) and S(z) take the form which is more pertinent to the MIMO systems analyzed later in Section IV.
The entries of C(z) are given by and c 2 . Hence, C(z) and S(z) satisfy

C(z)C(z) + S(z)S(z)
Their diagonal entries together can be viewed as the generalized polynomial singular-value, or eigenvalue, pairs (c i (z), s i (z)). Assuming s i (z) are stable, finite impulse response (FIR) filters, and hence invertible, then the ith generalized polynomial eigenvalue is given by: The generalized eigenvalues corresponding to zero-order coefficient matrix is: . Here, as in the scalar-matrices case, the condition s i [0] > 0 must be satisfied. In general, , for a particular frequency Ω. Then it is clear that s i (e jΩ ) must not have any spectral nulls, i.e. s i (e jΩ ) > 0, for all Ω.
Note that for a specific value of z, say z = z 0 , (7) is just the scalar GSVD as in (1). Hence, following Van-Loan in [4] we assume, without loss of generality, that X(z 0 ) is non-singular, for all z 0 , in other words, X(z) is invertible.

B. Polynomial CSD
In this subsection, we generalize the CSD to the case polynomial matrices (or PCSD). Consider the composite matrix Then we define the PCSD as: where U (z) : and S(z) : C → C P ×N satisfy the same properties as in PGSVD. The PU matrices U (z) and Z(z) can be calculated using the PSVD, described in the next subsection; whereas the PU matrix V (z) can be computed with the PQRD of S (z) = Q 2 (z)Z(z), i.e., Analogous to the scalar case, the PGSVD is identical to PCSD when

C. Polynomial SVD
Essentially, the PCSD relates the PSVD of Q 1 (z) to the PSVD of Q 2 (z) through (11). In this paper, the PSVD is implemented with the application of two PEVDs [27], [28], an approach that generalizes the well-known method for computing the conventional scalar SVD through application of two scalar EVDs [5]. Hence, the PU matrices U (z) and Z(z) in (11) can be obtained from the respective PEVDs of the two parahermitian matrices

D. PEVD Calculation
In this subsection, we propose an efficient PEVD algorithm for finding the PEVDs in (13). Sequential matrix diagonalisation (SMD) is a class of fast, iterative, time-domain algorithms that has been shown to provide a very accurate PEVD [18]. The maximum-element SMD (ME-SMD) algorithm uses a norm (correlation) based cost function to minimize the largest offdiagonal element at every iteration. However, this cost function makes the algorithm sensitive to larger correlation terms related to dominant channels. In situations where there are strong noise related correlations, there would be a tendency for ME-SMD to eliminate true cross-correlation terms rather than noise-related terms between weaker channels. This has the effect of limiting the achievable level of diagonalization [23]; thus restricting the degree to which the MIMO channels can be decoupled.
One way of alleviating this problem is to use a cost function that is, in proportion, equally responsive to correlations in all channels. The coding gain measure in [33] has this feature. Following the strategy in [23], a coding-gain optimized version of the ME-SMD, namely SMDC, can be obtained by optimizing the coding gain measure. In the following, only the cost-function aspect of the algorithm is discussed. For a full treatment of SMD, the reader if referred to [18].
Let Φ(z) : C → C M ×M be a parahermitian matrix, as in (13), with on-diagonals: Then the proposed SMDC algorithm maximizes the following sample coding-gain function at the lth iteration: where φ Note that one or more of the terms φ (l) i [0] may tend to zero, which would cause (14) to tend to infinity. This can be alleviated by applying a simple regularization technique [23]. The original proof for ME-SMD in [18] remains valid for the new SMD variant. Use of (14) is expected to enhance the spectral majorization and strong decorrelation performances of SMD. Hence, as well as being suited for MIMO-channel coding applications, our SMDC algorithm will be more applicable to subspace decomposition and data compression problems than its correlation-based counterparts.

E. Computing the Polynomial GSVD
Here we extend the approach in Section II-C to the case of polynomial matrices. The first step here is to compute the PQRD of H(z) = (H T 1 (z), H T 2 (z)) T , that is, where Q 1 (z) : C → C M ×N and Q 2 (z) : C → C P ×N , Q(z) is PU and R(z) : C → C N ×N is upper-triangular and nonsingular.
Analogous to the scalar case, this approach then calculates the PCSD of Q(z), via a PSVD and a PQRD, which yields the PGSVD canonical matrices U (z), X(z), V (z), C(z) and S(z). The matrices U (z), C(z) and Z(z) of the PCSD can be found directly by computing the PSVD of Q 1 (z). Then the PU V (z) and diagonal S(z) can be calculated from the PQRD of a submatrix defined by the first k columns of S (z) = Q 2 (z)Z(z), which yields V (z) from (7). The diagonal S(z) is simply given by: S(z) = V (z)S (z). Lastly, the PGSVD non-singular polynomial matrix X(z) is obtained thus X(z) = R(z)Z(z), where Z(z) is got from (11).
The pseudocode for the proposed PCSD algorithm is provided in Algorithm 1. Our approach makes use of the efficient, SMDbased PQRD (SM-PQRD) algorithm in [26].

IV. PGSVD APPLIED TO TWO-USER FREQUENCY-SELECTIVE MIMO CHANNELS
In this section, we show how the proposed PGSVD algorithm can be used as a broadband beamformer for point-to-two-point MIMO multicasting in a frequency-selective environment. We considered two user case only for simplicity of describing our algorithm, however, extension to point-to-multipoint MIMO multicasting would be straightforward.
upper-triangular R t (z). f: Update: (a) columns r in U (z) with those of U (z)V t (z); (b) rows r and columns j in C(z) with R t (z) diagonals.

A. MIMO-Channel Decoupling
Consider a transmitter with N antennas and two receivers, consisting of M and P antennas respectively, as shown in Fig. 1. Suppose that the source signal is s(n) whose Z-transfrom is s(z) : C → C. We stack the signals transmitted via N antennas into a signal vector whose Z-transform is s(z) : C → C N . This signal propagates through two different frequency-selective (FS) MIMO channels whose Z-transforms are H 1 (z) : C → C M ×N and H 2 (z) : C → C P ×N with the same number of columns. So without any transmitter precoding, the signals at the two receive-antenna arrays can be represented in Z-domain as To begin with, we assume perfect knowledge of the channelstate information. Later, in Section V-B, we assess the impact of this assumption on the performance of the proposed system.
We assume that the two MIMO frequency selective channels undergo block fading, so the channels do not change within a block, but can change between blocks. We can jointly factorize H 1 (z) and H 2 (z) by way of a PGSVD: where W (z) : C → C N ×N is invertible and U (z) : C → C M ×M and V (z) : C → C P ×P are PU. Here, Σ(z) and Λ(z), satisfy the same properties as the diagonal matrices C(z) and S(z) for the PGSVD defined in Section III-A. The diagonal elements of Σ(z) and Λ(z) represent the gains of the non-interfering VCs (CC or PC) from the source. The PC is a VC from S to U 1 (or U 2 ), having unit gain. The CC, unlike the PC, is point-to-two point VCs that cater to both the users simultaneously, with channel gain factors satisfy (9). The total number of PCs and CCs for each user depends on the number of antennas at the three terminals, however, the numbers of VCs are fixed. In Table II, we show how the numbers of PCs/CCs vary for the various antenna configurations. Note that the only system configuration that provides CCs and PCs for both users simultaneously is M, P < N under the constraint of (M + P ) ≥ N . Hence, we consider this class of two-user MIMO-system configuration for the purposes of investigating our PGSVD-based beamformer.
NB: An important feature of the proposed PGSVD beamformer is that it allows the use of different modulation or access schemes on different spatial channels. For example, CDMA/FDMA can be used for PCs while higher-order modulation schemes, such as 256-QAM, can be used for CCs, facilitating the flexibility to have heterogeneous applications and services.

B. Precoding Stage
The factor W (z) in (17) can be used to enable joint precoding of the two MIMO subsystems at the transmitter. Specifically, assuming W (z) is invertible, the two coded transmit signals are obtained by: However, matrix inversion is often fraught with difficulties, numerical instability being the main problem, especially for matrices with large condition numbers. Another potential issue with this approach is that the precoding is non-energy preserving, so it may attenuate or amplify the transmit power, which is undesirable. This may be alleviated somewhat by normalizing the transmit signal such that the average total transmit power is maintained at the desired level. One way of realizing the normalization is with a diagonal invertible polynomial matrix that is chosen to provide the needed transmit power normalization. However, this requires updated information about the average channel power. An alternative approach is to work with the pseu- is the generalization of the Moore-Penrose pseudoinverse to polynomial matrices. The matrix Σ † W (z) is formed by taking the reciprocal of all the non-zero entries of Σ W (z); i.e., So the precoding is performed with the application of W † (z) in (18), instead of W −1 (z). Hence, the coded broadcast channels are given by: An approximation to the pseudoinverse in (19) can be computed via the frequency domain, where a scalar Moore-Penrose pseudoinverse is applied at each frequency bin. However, this approach would ignore any phase-coherence that may exist between the bins.
A more accurate pseudoinverse can be obtained via a timedomain PSVD, such as the algorithm in [28]. After calculating the PSVD, Σ † W (z) can then be obtained simply by taking the reciprocal of the diagonals at each frequency. This method is naturally energy preserving, so it does not attenuate or amplify the transmit power. The divide-by-zero problem can be averted by adding an arbitrarily small value to the divisor. Of course, the success of this simple approach relies on the accuracy of the decomposition, and hence the diagonality of PSVD-generated Σ W (z). However, simple checks can be incorporate into this approach that enable one to maintain sufficiently accurate PSVDs.

C. Decoding Stage
The corresponding decoder, or reconstruction, matrices, applied at the receiver, are given by U (z) and V (z) from (17). These are PU matrices so they do not enhance the noise in the system. Therefore, the k decoded, non-interfering, broadband channels are simply given by: and

D. System Complexity
The computational complexity of the proposed MIMOchannel broadband beamformer is dominated by the number of arithmetic operations required by the PGSVD algorithm over a given number of iterations. As explained in Section III, computation of PGSVD via the PCSD involves application of both the SMDC (to calculate the PSVD) and SM-PQRD algorithms. Here, we see that the most expensive component of our PGSVD-based system is the polynomial matrix multiplications, which incurs the M 3 factor. Consequently, the computational burden increases significantly for larger spatial dimensions. However, in recent years, both hardware and algorithmic solutions have been proposed that make, respectively, polynomialmatrix multiplication [34] and iterative PEVD algorithms [35] more applicable to real-world problems.
As discussed in Section I, point-to-two point MIMO multicasting can also be achieved via a frequency-domain approach. The DFT-processed channel data has K uniformly spaced frequency bins. The OFDM-GSVD broadband beamformer operates by applying the scalar GSVD to the (scalar) channel matrices at each bin, independently, through use of OFDM. With the assumption that N = M , the baseline OFDM-GSVD scheme is found to perform O(M 2 K log (K)) + O(KM 3 ) work. Note that the O(M 2 K log (K)) term is due to computing the DFT of the M × M MIMO-channel matrix; and dominates the complexity for K >> M. For large M , the bulk of the computational cost is due to the matrix-matrix multiplication operations required by the GSVD [5].
Comparing the complexities of the two methods, with M common in both expressions, we see that the proposed PGSVDbased approach is computationally more demanding for large T and/or L; whereas the latter becomes more exhaustive for large K.

V. SIMULATION RESULTS
We provide two sets of simulation studies to demonstrate the performance of the proposed PGSVD algorithm and PGSVDbased, point-to-two point MIMO-channel equalization method described in the previous sections. The first set of results, in Section V-A, relate to evaluating the convergence behavior and efficiency of proposed PGSVD algorithm. The second set of simulations assess the efficacy of proposed broadband MIMOchannel beamformer for multicasting applications. Results are averaged across an ensemble of 1000 random scenarios of the two types of randomly generated system models, for the two sets of experiments.
Unless stated otherwise, the following PQRD/PSVDalgorithm parameter settings were used for all simulations, since they provided a good trade-off between decomposition accuracy and computational cost and/or communications efficiency. The algorithms were allowed to run for l = 200 iterations. We found that beamformer performance did not improve by a notable amount for l ≥ 200. The algorithm convergence parameter was set to = 10 −5 , which is an upper bound on the dominant offdiagonal element. The polynomial-matrix trim function in [17] was applied to the updated canonical matrices with μ = 10 −5 . This function truncates the matrix coefficients at large lags (in the tails of polynomial matrices) containing a μ-th of the total power.

A. Algorithm Convergence
We will now illustrate the validity and investigate the decomposition quality of the proposed CSD-based PGSVD algorithm by means of some simple but insightful experiments. The algorithm was applied to two order-5 polynomial matrices A : C → C 3×3 and B(z) : C → C 3×3 . The coefficients of A(z) and B(z) were drawn from a zero-mean, complex-Gaussian, wide-sense stationary process.
In the first experiment, the input polynomial matrices were reconstructed from the canonical-form matrices generated by the PGSVD algorithm; i.e., A(z) ≈ U (z)C(z) X(z). In Fig. 2(a) we show the input A(z) along with the reconstructed A(z) in  Fig. 2(b). Notice the strong likeness between the two polynomial matrices. The maximum off-diagonal energy in the matrices C(z) and S(z) was of the order of 10 −6 , which meant Fig. 3(a) and (b), respectively. These results demonstrate that the proposed PGSVD algorithm is able to provide a joint factorization of the two polynomial matrices A(z) and B(z) to a very good approximation.

C(z)C(z) + S(z)S(z) = I was approximately satisfied, as indicated by the product terms C(z)C(z) and S(z)S(z) shown in
The PQRD and PSVD algorithms operate by means of iterative procedures and are known to converge [27]. As described in Section III, the PGSVD algorithm operates by using those routines within the well established structure of the CSD algorithm. As such, the PCSD algorithm, on which our PGSVD algorithm is based, does not involve any iterations; so the concept of convergence does not apply to the PCSD and, by construction, nor does it apply to the PGSVD. Therefore, in the following, we explore the convergence of the underlying PQRD and PSVD algorithms, since these reflect the decomposition quality of the proposed PGSVD/PCSD algorithm.
The PSVD and PQRD operate by transferring, respectively, the off-diagonal and beneath diagonal coefficients onto the main diagonal. For a fair comparison, we investigate the magnitude of the largest off-diagonal/beneath-diagonal coefficient at each iteration. Fig. 4 shows the magnitude, g, of the dominant off-diagonal/beneath-diagonal term as a function of iteration number, l. It can be seen that both algorithms converge in ∼200 iterations. Also, note that the PSVD algorithm converges slightly faster than the PQRD algorithm, which could be because the former converges monotonically, whereas the latter converges non-monotonically.

B. BER Performance
In the following, we consider the frequency-selective, twouser, MIMO communications system in Fig. 1, where, data is transmitted from a single source S to two users U 1 and U 2 . We simulated transmitting 1100 coded BPSK/QPSK modulated symbols through each of the two channels. We assume block fading and perfect knowledge of channel-state information. Convolutional source coding with the code rate of 1 2 was used along with appropriate equalization at the two receivers, as in [25]. The propagation environment was modeled using channels possessing an exponential power-delay profile, which are typically seen in macro-cellular communications systems [25]. A five-path fading model was used for the frequency-selective channel from the nth source antenna element to the mth/pth antenna sensor of U 1 /U 2 , where m = 1, . . . , M, n = 1, . . . , N and p = 1, . . . , P . The channel matrices have polynomial entries with exponentially decaying envelopes, which are: for τ = 1, 2, . . ., 5, Here, h mn [τ, t] is the τ th multipath coefficient between the nth transmit antenna and mth receive antenna w.r.t time interval t, η[t] is a zero-mean, complex Gaussian random variable; and the parameter α is chosen to normalize the average channel gain to unity. Here, ψ denotes the decay factor, which was set to ψ = 0.8 in generating H i (z), i ∈ {1, 2}.
The following simulation results are for different antenna configurations, chosen to explore the performance of realistic CCs and/or PC scenarios. Note that, for the same configurations, BER results for the frequency-selective environments observed here analogously mirror the narrowband results reported in [3].  Note also that the reason behind very low bit errors in our results is partly due to the fact that we use error-correcting coding followed by equalization.
In Fig. 5, we show coded BER performances corresponding to an antenna configuration given by {N, M, P } = {3, 3, 3}. Our PGSVD beamformer produces three CCs for each of the two users. It can be seen that the BERs experienced by U 2 for CC k are similar to those of U 1 for CC (4−k) , k ∈ {1, 2, 3}. A striking result is the symmetry observed here, which is due to the fact that M = P and the coefficients H i (z) are drawn from a symmetric and independent distribution. Notice that the BER shown for U 1 improves from CC 1 to CC 3 , whereas, the converse is true for U 2 . This observation is consistent with the fact that the cosine coefficients c i (z), from PGSVD, appear in descending order, as in (9).
The coded BER performances of the PC/CCs for the {N, M, P } = {4, 3, 3}-antenna case is shown in Fig. 6. We see that the: (i) CCs, for each user, have a different BER performance; (ii) BERs of U 1 for CC 1 /CC 2 are very similar, if not the same, to those of U 2 for CC 2 /CC 1 ; (iii) BER relating to the PCs of both users are approximately the same; and (iv) the PCs outperform the CCs. Results (i) to (iii) are again a consequence of  the statistical symmetry in the two MIMO channels. Result (iv) is because of the extra spatial DoF possessed by the source coupled with the fact that each PC caters to just one user. Fig. 7 shows coded BER performances for {N, M, P } = {4, 3, 2} case. As with previous configurations, the PCs show very similar BERs; moreover, the PCs outperform the CCs, which is because most of the available DoFs go towards enabling the PCs. This PC-to-CC beam distribution provided by the PGSVD-based beamformer maximizes the multiplexing gain. Notice also that the CCs of the two users experience slightly different BERs. The reason for this is that the CC 1 of U 1 is equivalent to the third sub-channel gain factor which is smaller than U 2 's CC 1 corresponds to the second sub-channel gain factor. Fig. 8, shows the coded BER performance for the {N, M, P } = {4, 2, 2} antenna configuration. This is an interesting case in that the source has as many spatial DoFs as the combined DoFs of the two users. In this case, one would expect the PGSVD-based beamformer to approximate a broadband zero-forcing beamformer -able to minimize CCI equally well across all PCs. However, the beamformer generates four PCs with differing BERs; notably, U 1 's PC performs significantly worse than U 1 's PC 2 and U 2 's PC 1 , especially at high SNR. This is mainly due to the fact that the PSVD produces large-order U (z) matrices in (11) for encoding this PC. Recall from Section III-C that our approach to implementing the PSVD involves the application of two PEVDs, which effectively doubles the order of the PSVD polynomial matrices. One possible way of alleviating this problem is to consider a PQRD-based PSVD. The PQRD of Foster et al. in [27] could be applicable, once modified to handle underdetermined systems.

C. BER Performance -Channel Errors and Uncertainties
Note the absence of a bit-error floor in all the BER results presented here thus far. This is indicative of the fact that the PGSVD-based beamformer minimizes the CCI between the two users. Elimination of CCI is to be expected given exact knowledge of the two MIMO-channel polynomial matrices. Perfect channel knowledge is typically assumed, in the literature, when evaluating channel coding schemes [1], [2], [25]. However, in practice, the channel state information may have to be estimated and tracked, especially in a fast varying environment. The channel-estimation process gives rise to estimation errors which negatively impact the BER performance of the MIMOchannel communication systems. Therefore, it would be prudent to investigate the robustness of our PGSVD-based beamformer to errors and uncertainties involving the MIMO channels.
Firstly, we assess system performance under channel estimation errors. The errors in the two estimated MIMO-channel matrices can be represented by additive channel-estimation error matrices E 1 (z) : C → C M ×N and E 2 (z) : C → C P ×N with entry-coefficients drawn from a zero-mean complex Gaussian process with variance σ 2 . Then the erroneous channel estimates, given byĤ are used to compute the PGSVD-based beamforming matrices instead of H i (z). In Fig. 9, we show the effect of channel-estimation errors on the coded BER performance of U 1 's CC 1 for the {N, M, P } = {3, 3, 3}-antenna configuration. Expectedly, total bit errors increase markedly as σ 2 increases, producing BER floors at high channel-estimation errors. For example, there is over an order of magnitude reduction in BER performance, for σ 2 = 0.01, at 10-dB SNR. This performance degradation highlights the reliance of PGSVD-based beamforming on MIMO-channel estimation accuracy. The semi-blind broadband MIMO-channel estimator of Hassan et al. in [25] can be used to provide good channel estimates when considering such MIMO systems; the investigation of which is left for future work as it would detract from the main goals of this work.
With the next set of results we demonstrate the robustness of our broadband beamformer to rank-deficient MIMO links, which models MIMO-channel spatial correlations. Experiments conducted were for the channel model in Section V-B and the {N, M, P } = {3, 3, 3}-antenna configuration, for which the channel matrices have a (full) rank of 3. The rank of the channel matrices was reduced to 2 for 16 of the frequency bins, uniformly distributed across the spectrum.
As discussed in Section I, a competitor method for point-totwo point MIMO multicasting can be obtained by performing per tone GSVD beamforming via OFDM (or OFDM-GSVD). In this frequency-domain approach, the scalar GSVD is used to jointly factorize the channel matrices at each frequency bin. For the OFDM-GSVD method used here, we found that working with 512 DFT points struck a good balance between BER performance and PU filter lengths.
In Fig 10, we compared the BER performances of the proposed approach with that of the OFDM-GSVD method under full-rank and rank-deficient channel conditions. Here, we only show the performance of the dominant CC for U 1 , which was CC 3 in this case. Clearly, both systems show similar performance for the full-rank case, with the former attaining slightly lower bit errors at higher SNRs. However, a striking result is that, for rank-deficient MIMO links, the PGSVD-based system outperforms the OFDM-GSVD approach; this performance gap increases significantly at higher SNRs. The performance degradation of OFDM-GSVD is because information of each transmitted symbol is constrained to one frequency bin. This problem can be alleviated with feedback of channel condition Fig. 11. BER performance comparison of the proposed PGSVD-based and OFDM-GSVD beamformers, for Jakes MIMO channels. information to the receiver; however, at the cost of reduced usable bandwidth and higher overhead. Contrastingly, with our time-domain approach, the entire channel spectrum is used to communicate the information in each symbol, making the system more robust to this type of channel corruption.
Lastly, we provide validation of our approach for the case of non-flat fading channels using the Jakes fading model. The frequency-selective transmission channel from the nth source antenna element to the mth/pth output of a 3 × 3 MIMO channel is described by a seven-path fading model. The multipath problem is modeled by using the velocity variation of the source from 3 to 4 km/h. Fig. 11 shows the BER performances of the proposed approach applied in the case of a Jakes MIMO channel. For comparison, we include the bit errors made by the OFDM-GSVD baseline scheme. As before, the BERs of the dominant CC for U 1 are shown. Results show that the proposed PGSVD-based beamformer compares favorably to the OFDM-GSVD approach. The reason for this is that, as in the case of the rank-deficient channel problem, the OFDM-based method decodes on a binby-bin basis, and so will produce more bit errors for bins that coincide with parts of the spectrum where the channel fading is great. Conversely, the proposed system is more robust to these types of errors since the information in each symbol is spread across the entire frequency bandwidth.

VI. CONCLUSION
In this paper, we have proposed a broadband beamformer that can generate point-to-two point MIMO channels for multicasting in frequency-selective environments. To this end, we proposed a generalized singular value decomposition (GSVD) for polynomial matrices, achieved by extending the cosinesine decomposition (CSD) to the polynomials. The polynomial GSVD (PGSVD) is used to design the required broadband beams for efficient broadband MIMO multicasting. The main features of the proposed approach are: r Yields broadband virtual channels that simultaneously accommodate for both point-to-point private channels (PCs) and point-to-two point common channels (CCs).
r Robust to channel-estimation errors, rank-deficient links and non-flat fading channels; and produces favorable BER performances as compared to a baseline OFDM-based scheme, which is also proposed in this paper.
r Supports the use of various and varied modulation/access schemes on the different virtual channels (VCs), depending on end-user needs. As an example for the last feature, our MIMO system can allow CC users to benefit from the new multi-gigabit wireless personal and local area network specifications, such as the IEEE 802.11ad/ay, whilst simultaneously providing flexibility to use other conventional access schemes, such as point-to-point CDMA/FDMA, for the PCs.