Decoupling of Broadband Optical MIMO Systems Using the Multiple Shift SBR 2 Algorithm

Polynomial singular value decomposition (PSVD) plays a very important role in broadband multiple-input multiple-output (MIMO) systems. One of its applications lies in the decoupling of MIMO convolutive mixing channel matrix in order to recover the transmitted signals corrupted by the channel interference (CI) at the receiver. In this paper, a novel algorithm, known as multiple shift second order sequential best rotation (MS-SBR2), is proposed to compute the approximate PSVD of the broadband MIMO channel matrix. Experimental examples, including a measured (2 × 2) optical MIMO channel impulse response using the multi-mode fiber (MMF) testbed, are presented to examine the proposed algorithm. Bit error rate (BER) performances are evaluated among different transmission schemes. In addition, power allocation (PA) schemes are investigated to further optimize the BER performance. Keywords—Multiple shift SBR2, polynomial SVD, precoding, equalization, multi-mode fiber, optical MIMO systems, power allocation.


I. INTRODUCTION
An explosive development of MIMO technology has been witnessed in wireless communication systems over the last decade.Compared to single-input single-output (SISO) systems, MIMO systems are capable of achieving higher data rates and transmission reliabilities by using the techniques of spatial multiplexing and transmit diversity.Aiming to increase the fiber capacity, the concept of MIMO in optical transmission systems has also attracted intensive research interests [1], [2].
Due to the multipath effect in broadband MIMO systems, the channel is characterized by frequency-selective fading, so apart from the channel interference caused by the MIMO components, there also exists inter-symbol interference (ISI) between the transmit symbols.Provided the approximate channel length is known at the transmitter, a standard approach of combating the ISI is to use multi-carrier modulation techniques, such as orthogonal frequency division multiplexing (OFDM) with cyclic prefix which can divide the spectrum into a number of narrowband channels.In other words, the frequency-selective or broadband MIMO channel is turned into a set of parallel frequency-flat or narrowband MIMO channels where the ISI does not exist anymore, and each narrowband channel can be independently addressed using Manuscript received October 30, 2016; revised March 30, 2017.Z. Wang and J. G. McWhirter are with the School of Engineering, Cardiff University, Cardiff CF24 3AA, Wales, UK (e-mail: wangz49@cardiff.ac.uk, mcwhirterjg@cardiff.ac.uk).
existing narrowband optimal techniques.This type of MIMO-OFDM system was well implemented by a spatio-temporal vector coding (STVC) [3], [4] communication structure combined with the singular value decomposition (SVD) based equalization technique [5].Another widely applied approach is based on optimal filter bank transceiver techniques [6] which involve block processing and guard intervals to eliminate interblock interference (IBI).These two traditional techniques are in common in the sense of eliminating the polynomial nature of the channel which corresponds to the IBI due to the timedispersive nature of the channel, and essentially they are all designed to convert the broadband problem into narrowband ones.

II. THE STATE OF THE ART
The work which we present here is related with the broadband MIMO decoupling method proposed in [7], [8].This method is essentially distinct from the traditional techniques like OFDM or optimal filter bank transceiver, as it can be applied to decouple the broadband (polynomial) MIMO channel directly, instead of converting it into narrowband channels.It mainly consists of two steps.The first step is based on the PSVD that was used to diagonalize the broadband MIMO channel matrix in order to remove the CI, which results in the frequency-selective MIMO channel decoupled into a number of independent frequency-selective SISO channels.The second step involves removing the remaining intersymbol interference (ISI) for each SISO channel, which can be implemented by further equalization techniques, such as zero-forcing equalization or maximum likelihood sequence estimations (MLSE).
There are different ways of calculating the PSVD of a polynomial matrix, such as using polynomial matrix QR decomposition to formulate the PSVD [9], PSVD based on generalized Kogbetliantz transformations [10], and PSVD by polynomial matrix eigenvalue decomposition (PEVD) method [11], which is analogous to how the scalar matrix eigenvalue decomposition (EVD) can be used to generate the singular value decomposition (SVD) of a matrix.In terms of the PSVD by PEVD method, the second order sequential best rotation (SBR2) algorithm [12] has been used in the existing literatures.However, an improved version of the SBR2 algorithm, i.e.MS-SBR2 [13], has been recently proposed by the authors for calculating the PEVD of polynomial matrices.The improved algorithm can provide much faster convergence than the SBR2 algorithm when dealing with high dimension polynomial matrices.In other words, the diagonalization of bigger MIMO doi: 10.11601/ijates.v6i1.207channel matrices can be implemented faster than that of using the SBR2 algorithm.
The main contributions of this work are the exploitation of the proposed PSVD by MS-SBR2 method in the application of solving the broadband MIMO decoupling problem, comparisons against the PSVD by SBR2 method, the discussion of the accuracy of the PSVD approach.The results presented in this paper include a simulated broadband channel matrix example which is designed to test how the PSVD method works.More importantly, a (2 × 2) optical MIMO channel which comprises a 1.4 km MMF and optical couplers at both ends is designed to examine the BER performance of a optical MIMO system in which the channel impulse responses are measured for the operating wavelength of 1576 nm [8].In particular, transmission and power allocation schemes are employed to bring further improvement with respect to the BER performance.
The remaining parts of the paper are structured as follows.The convolutive MIMO channel model with polynomial matrix representation is described in Sec.III.In Sec.IV we introduce the idea of broadband MIMO channel decomposition, i.e.PSVD.Sec.V presents the proposed MS-SBR2 algorithm for calculating the PSVD.Simulation results and conclusions are shown in Sec.VI and Sec.VII, respectively.

III. MIMO CONVOLUTIVE MIXING MODEL
Given a frequency selective MIMO link with n T inputs and n R outputs, the convolutive mixing channel can be modelled as a polynomial matrix with an indeterminate variable z −1 given by where τ, T ∈ Z and C[τ ] ∈ C nR×nT denotes the polynomial coefficient matrix at time lag τ and c νµ (z), , is the polynomial matrix entity which represents the channel impulse response between the µ-th input and the ν-th output.It takes the form of where c νµ [τ ] denotes a non-zero element of the symbol rate sampled overall channel impulse response at the τ -th lag.In this case there are T + 1 lags in total for each SISO channel.Throughout this paper, polynomial matrices and vectors are denoted as underscored boldface letters.
In the context of optical MIMO systems, it should be noticed that the group delays in a MMF optical channel belong to a fixed set of values in contrast with that of the wireless channel which can change from one realisation to another [14].Assuming that the transmit signal is represented by s (z), the convolutively mixed received signal x (z) can be expressed as where n(z) represents the additive noise which has variance of σ 2 I nR .
IV. BROADBAND MIMO CHANNEL DECOMPOSITION VIA PSVD One potential application of PSVD is to enable communication over a broadband MIMO system in which the channel matrix is represented by a polynomial matrix as shown in (1).In this case, provided the channel matrix has firstly been estimated, the PSVD then can be used to simplify a MIMO channel equalization problem into a set of independent SISO problems.In other words, the CI can be removed by performing the PSVD to the channel matrix C(z), which can be expressed as [11] where we assume n R ≥ n T , and Γ(z) is a diagonal polynomial matrix with n = n T diagonal elements, i.e.
Here the notation { } over a polynomial matrix denotes the paraconjugate operation which is computed by performing Hermitian transpose {•} H to all the polynomial coefficient matrices U[τ ] and timereversing all entries inside, i.e.U(z) = U H (1/z).
Note that U(z) and V(z) are acting as the multichannel all-pass filters which can transform the frequency selective MIMO channel into a number of independent frequency selective SISO channels while still preserving the total signal energy [15].In this paper, the PSVD in ( 4) is implemented by calculating the PEVD of two para-Hermitian polynomial matrices C(z) C(z) and C(z)C(z), which take the form as and Further details about the PEVD algorithms will be discussed in the following section.To eliminate the CI, the source signal s(z) is filtered by the paraunitary transformation matrix V(z) at the transmitter, i.e. s (z) = V(z) s(z), and the received signal x (z) is pre-multiplied by U(z) such that x(z) = U(z) x (z) at the receiver, which results in where w(z) = U(z)n(z).Note that neither the transmit power is increased, nor the channel noise is enhanced here.Unlike the conventional SVD-based method, each diagonal element (also called layer) in Σ(z) is frequency-selective and hence ISI occurs.In order to remove the ISI, layer-specific T-spaced zero forcing equalizers [8] are adopted in this paper and therefore this equalization scheme is entitled T-PSVD.A block diagram of the proposed communication system can be depicted in Fig. 1.
Fig. 1.Block diagram of the proposed communication system using the PSVD.

V. PSVD USING THE MS-SBR2 ALGORITHM
As mentioned above, the PEVD method can be used to formulate the PSVD problem in (4), and the idea of PEVD has been generalized as [12] where R(z) ∈ C M ×M is a para-Hermitian matrix, such that R(z) = R(z).H(z) is a paraunitary matrix which aims to diagonalize R(z) by means of paraunitary similarity transformation, and D(z) is (ideally) a diagonal matrix.The approximation sign in (8) indicates that a PEVD does not necessarily exist if the paraunitary matrix H(z) only contains FIR components.However it has been proved that a very close approximation can be achieved by letting the polynomial order of H(z) grow arbitrarily large [16].This is an iterative process which transforms all the off-diagonal elements in R(z) onto the diagonal.Several algorithms [12], [17], [18] exist for calculating the PEVD, however, this paper is concerned only with the MS-SBR2 algorithm previously presented by the authors in [13].
The MS-SBR2 algorithm is an improved version of the SBR2 algorithm in terms of the algorithm convergence speed.Basically it adopts the advantages of less computational cost from SBR2 and the faster convergence from MSME-SMD [18], which seems to provide a compromise between the SBR2 and the SMD algorithm family.For the following part of this section, the SBR2 algorithm is firstly introduced before we move forward to the MS-SBR2 algorithm, and a numerical example is chosen to assess the performance of the algorithms, and then followed by a discussion of the computational accuracy.

A. Second Order Sequential Best Rotation (SBR2) Algorithm
The SBR2 algorithm was designed to iteratively eliminate the off-diagonal elements for para-Hermitian matrices by using the paraunitary similarity transformations shown in (8).At the i-th iteration, the SBR2 algorithm [12] starts by locating the maximum off-diagonal element r To find the maximum off-diagonal element, we define a matrix S (i) [τ ], which contains only the upper triangular elements in R (i−1) [τ ] with the remaining elements set to zero.Thus the location of r where j (i) , k (i) and τ (i) are the corresponding row, column and time lag index.An elementary delay matrix P (i) (z) and Jacobi rotation Q (i) are sequentially applied to bring r Step 1 Step 2 Step 3 coefficient matrix R (i−1) [0], and then rotate its energy onto the diagonal.This results in R (i) (z) given by where P (i) (z) is expressed as A 3-dimensional illustration which shows the procedure of the i-th iteration in SBR2 is described in Fig. 2. Thus the elementary paraunitary matrix E (i) (z) can be expressed as The algorithm continues its iterative process until all the offdiagonal elements in R (i) (z) are smaller than a given threshold which can be set to a very small value to achieve sufficient accuracy.Assuming that the algorithm has converged at the N -th iteration, the diagonalized para-Hermitian matrix in (8) takes the form of and the generated paraunitary polynomial matrix is given by

B. Multiple Shift SBR2 (MS-SBR2) Algorithm
The MS-SBR2 algorithm uses a distinguishing search strategy of the off-diagonal elements which is akin to that of the MSME-SMD algorithm, so that it can achieve the diagonalization with less iterations than the SBR2 algorithm.For the i-th iteration, the MS-SBR2 algorithm involves multiple shifts operations P (i) (z), followed by a sequence of Jacobi rotations Q (i) .Therefore the resulting para-Hermitian matrix is computed by where l=1 Q (l,i) and L (i) denotes the total number of off-diagonal elements shifted to the zero-lag coefficient matrix at the i-th iteration (L (i) ∈ Z, 1 ≤ L (i) ≤ M/2 ).Accordingly the delay matrix at the l-th delay stage within i-th iteration is represented by and the elementary paraunitary matrix can be expressed as . Note that when L (i) = 1, the MS-SBR2 algorithm is identical to the SBR2 algorithm.Another motivation of introducing the multiple shift idea into the SBR2 algorithm is that it permits us to minimize the order growth of polynomial matrices by making all row (column) shifts in the same direction, which can potentially reduce the computational cost of the algorithm [19].The PEVD algorithms are assessed in terms of the normalized offdiagonal energy η (i) at the i-th iteration, and it is defined as where the notation • F denotes the Frobenius norm.The comparison between these two PEVD algorithms is calculated via Monte Carlo simulations over an ensemble of 100 different random 10×10 para-Hermitian matrices of order 5, which can be generated from matrices A(z) ∈ C 10×10 of order 3 with i.i.d.zero mean unit variance complex Gaussian entries, such that R(z) = A(z) A(z).Fig. 3 shows the normalized off-diagonal energy η (i) versus the iteration index i.Obviously the MS-SBR2 algorithm requires much fewer iterations than the conventional SBR2 algorithm to achieve the same level of diagonalization.However, it should be noticed that each iteration within MS-SBR2 involves more rotation steps, which means the computational costs between them are comparable.Nonetheless, the MS-SBR2 algorithm was found to converge faster than SBR2 as shown in Fig. 4. For further details of the MS-SBR2 algorithm, see [13].

C. Accuracy of the Decomposition
There are two main factors which can affect the accuracy of the decomposition.Firstly, since the decomposition is performed upon the two para-Hermitian matrices C(z) C(z) and C(z)C(z) as shown in ( 5) and ( 6), the resulting diagonal matrix Σ(z) might be less accurate than that found by the way of operating the decomposition directly upon the channel matrix C(z).Secondly, in the sense of broadband MIMO Comparison of normalized off-diagonal energy η (i) between SBR2 and MS-SBR2 algorithms, showing ensemble averages versus mean execution time (measured in MATLAB R2014a on a PC with configurations Intel Core i7-3770T CPU@2.50GHz and 16 GB RAM).
decoupling, a strictly diagonalized channel matrix is required.However, the proposed PSVD method can only generate an approximately diagonal matrix subject to the pre-specified stop condition of the algorithm, so there will be errors when assuming all the off-diagonal elements of Σ(z) are equal to zero.In addition, due to the fact that the orders of the polynomial matrices increase as the iteration goes throughout the PEVD process, the equalization becomes difficult when the order of Σ(z) is too large.Therefore polynomial order truncation operations [20], [21] are usually required in order to keep orders as small as possible and reduce the computational cost of the algorithm, which can also cause a very small proportion of the total Frobenius norm of the matrix being eliminated.To assess how well the proposed PSVD method performs, the error metric of the PSVD is defined as where Σ(z) is equal to Σ(z) with all the off-diagonal elements set to zero.

A. Example 1
To demonstrate the proposed PSVD method, a polynomial or broadband channel matrix C 1 (z) ∈ C 4×3 was generated to describe the propagation of three source signals onto four sensors.Each of the polynomial entries of the matrix was chosen to be order-5 FIR filter, where both the real and imaginary parts were drawn randomly from a uniform distribution in the range [−1, 1].A graphical representation of this channel matrix is plotted in Fig. 5.With the truncation and stopping parameters set as µ = 10 −4 and = 10 −3 , the diagonalized channel matrix Σ 1 (z) from the PSVD by MS-SBR2 method is shown in Fig. 6, and the performance comparison against the PSVD by SBR2 method is summarized in Tab.I.The results presented in the table demonstrate that with the same level of decomposition achieved, i.e. g = 9.91 × 10 −4 , the PSVD by MS-SBR2 method outperforms the PSVD by SBR2 method in terms of the number of iterations, relative error and computational time.
In addition, akin to an ordered SVD with singular values in descending order, the power spectral densities (PSDs) σ mm (e jΩ ) = σ mm (z)| z=e jΩ , m = 1, 2, • • • , M of the diagonalized matrix Σ 1 (z) has the spectral majorisation property [22] such that for all normalised angular frequency values Ω Compared to the original PSDs of C 1 (z) shown in Fig. 7, the PSDs of the diagonalized channel matrix Σ 1 (z) are spectrally majorised, which intuitively indicates that the MIMO channel has been decoupled into a set of independent SISO channels.

B. Example 2
This example shows the practical application of the MS-SBR2 algorithm to a measured (2 × 2) optical MIMO channel.In fiber-optic systems one method to realize a MIMO transmission is to carry the data streams on different optical modes through a few-mode or multi-mode fiber (MMF) [2],    9 demonstrate that spatially diverse channels are generated.However, an ideal separation of the two channels is hard to achieve since mode mixing usually occurs in the mode multiplexing and demultiplexing process which is implemented by fusion couplers, and during the transmission through the fiber.An overview of the testbed used for measuring optical MIMO impulse responses is shown in Fig. 10.Here the impulse responses of the (2 × 2) optical MIMO channel, consisting of a 1.4 km MMF, fusion couplers and differently aligned SMFs, are measured at an operating wavelength of 1576 nm with the aid of signal deconvolution [8].The measured impulse responses are then sampled at the symbol rate of 620 MHz and used to constitute the channel matrix C 2 (z) as plotted in Fig. 11.
Applying PSVD to this frequency-selective MIMO channel F , is given by 1.26 × 10 −6 .The value of ε is negligibly small compared with the input energy, which means that the CI has been significantly eliminated.
The BER quality is studied by using fixed transmission modes with a spectral efficiency of 8 bit/s/Hz, and the analyzed quadrature amplitude modulation (QAM) constellation arrangements are depicted in Tab.II.In addition to bit-loading, the allocation of the transmit power to the activated layers is needed to optimize the BER performance.In the optimum case, the PA aims at equalizing the BERs over all layers.However, this approach is computationally complex and hence it is simplified by just equalizing the SNRs as a sub-optimum solution.After the application of the T-PSVD, the decoupled MIMO layers exhibit decreasing SNRs at higher layers.This conforms with the spectral majorisation property shown in example 1.The SNR conditions subsequent to the PSVD equalization and the conditions after power allocation are illustrated in Fig. 13.The BER performance results for T-PSVD equalization, obtained by applying the MS-SBR2 algorithm for calculating the PSVD, are depicted in Fig. 14 for a range of SNRs.Here E s is the transmit signal energy and N 0 denotes the constant noise power spectral density of the additive white Gaussian noise.As seen from the graph, the (256, 0) QAM transmission scheme shows the best performance results.It should be noted that no PA is needed for the (256, 0) QAM transmission mode.In addition, when activating multiple numbers of layers the benefit of using the equal SNR power allocation method is clearly visible.As the dimension of the polynomial matrices C 2 (z) C 2 (z) and C 2 (z)C 2 (z) is 2 × 2, this means that only one off-diagonal element can be eliminated for each iteration when computing the PEVD via the MS-SBR2 algorithm, i.e.L (i) = 1, ∀ i.Therefore, the MS-SBR2 algorithm operates the same as the SBR2 and there is no difference between them in terms of the BER performance for this example.However, in real world implementations involving more sources and sensors, the benefits of the MS-SBR2 algorithm could come into play.
It is expected that PSVD based MIMO systems can offer the same BER performance compared to systems based on STVC with SVD equalization, as it is suggested by the achievable spectral efficiencies presented in [24].The major advantage of MIMO systems based on PSVD is that they do not require a block-wise transmission.

VII. CONCLUSION
We have investigated how the proposed MS-SBR2 algorithm can be used in the application of decomposing the channel matrix of a measured (2 × 2) broadband optical MIMO system.Furthermore, different transmission schemes have been employed to illustrate the BER simulations.In particular, the power allocation scheme has been utilized to further optimize the BER performance.Simulation results have shown that the activation of all transmission layers does not necessarily lead to the best BER performance.On the contrary, the (256, 0) QAM with the T-PSVD equalization scheme seems to achieve the best performance in the studied example.

Fig. 2 .
Fig. 2. A 3-dimensional illustration of a 5 × 5 polynomial matrix example, showing the i-th iteration process using SBR2; Assuming the maximum offdiagonal element r (i) jk [τ ] found is at the location of {1, 5, 2} represented in green color, step 1 shows the location information; Step 2 describes the corresponding row and column shift operations; Step 3 is to transfer the pairwise maximum elements r (i) jk [τ ] and r (i) kj [−τ ] onto diagonal (only zerolag coefficient matrix is shown here for visibility purpose)[12],[18].
Fig. 4.Comparison of normalized off-diagonal energy η (i) between SBR2 and MS-SBR2 algorithms, showing ensemble averages versus mean execution time (measured in MATLAB R2014a on a PC with configurations Intel Core i7-3770T CPU@2.50GHz and 16 GB RAM).

Fig. 5 .Fig. 6 .
Fig. 5.The stem plot of the 4 × 3 broadband MIMO channel matrix C 1 (z) showing the magnitude of channel impulse response at different time lag.

[ 23 ]
. For the excitation of different modes, certain singlemode fiber (SMF) to MMF alignments with varying radial offsets δ are used in this work.The spatial diversity of the optical MIMO channel is visualized by showing the measured spatial intensity distributions when exciting the two optical MIMO inputs of the (2 × 2) system separately.The measured patterns depicted in Fig.

Fig. 9 .
Fig. 9. Intensity distribution patterns at the end face of a MMF fiber when launching centric δ = 0 µm (left) and launching with an eccentricity of δ = 15 µm (right) measured at a wavelength of 850 nm; the dashed line represents the 50 µm core diameter.
results in layers having a time-dispersive characteristic and hence ISI occurs on each layer as shown by the diagonalized matrix Σ 2 (z) in Fig.12.The remaining ISI is removed by applying the T-spaced zero forcing equalization as mentioned before.The equalizers modify the noise power on each layer differently, which is expressed by the weighting factors θ , with denoting the layer index.These factors determine the layer specific SNRs and hence also the total BER performance[8].In this example, the noise weighting factors for each layer are computed as θ 1 = 37.22 and θ 2 = 4243.46.In addition, the remaining off-diagonal energy ε, defined as ε

Fig. 13 .
Fig. 13.Illustration of the remaining SNRs in T-PSVD systems without applying PA (left) and with PA (right).The color black refers to high and white to low SNR values.

Fig. 14 .
Fig. 14.BER with PA (dotted line) and without PA (solid line) by applying the T-PSVD equalization scheme, showing the comparisons among different transmission modes when transmitting over the 2 × 2 optical MIMO channel.

TABLE I RESULTS
OF APPLYING THE PSVD BY PEVD ALGORITHMS TO C 1 (z) IN EXAMPLE 1 Fig. 7.The PSDs c mm (e jΩ ), m = 1, 2, 3 of the channel matrix C 1 (z).