Partial Fractional Fourier Transform (PFrFT)-MIMO-OFDM for Known Underwater Acoustic Communication Channels

Communication over doubly selective channels (both time and frequency selective) suffers from significant intercarrier interference (ICI). This problem is severe in underwater acoustic communications. In this paper, a novel partial fractional (PFrFT)-MIMO-OFDM system is proposed and implemented to further mitigate ICI. A new iterative band minimum mean square error (BMMSE) weight combining based on LDLH factorization is used in a scenario of perfect knowledge of channel information. The proposed method is extended from SISO-OFDM configuration to MIMO-OFDM. Simulation results demonstrate that the proposed PFrFT-LDLH outperforms the other methods in the SISO-OFDM scenario and that its performance can be improved in MIMO-OFDM scenarios.


Introduction
Underwater acoustic communication (UWA) suffers from significant time delays reaching fractions of seconds as well as severe Doppler spread attributed to relative motion between the transmitter and receiver. This problem is exacerbated due to the low speed of sound in water (1500 m/s). Orthogonal frequency division multiplexing (OFDM), is widely used in wireless communications that exhibits high dispersion due to its superior resistance to inter-symbol interference (ISI) [1] and its low complexity. Equipped with multiple transmitters and receivers, multiple-input and multiple output (MIMO) techniques have been applied to the OFDM system in UWA. The resulting MIMO-OFDM was shown to not only improve the data rate over the bandwidth-limited channels through spatial multiplexing, but also enhance the systems quality by means of spatial diversity. Different from spatial diversity, spatial multiplexing enables the high-rate encoded data streams to be separated into several segments with low data rate. Each segment is sent independently over different transmitters at the same frequency band and processed by receivers parallelly. In [2][3][4], a MIMO-OFDM system with spatial multiplexing was proposed in UWA communication resulting in a block-by-block receiver design.
With an increase of subcarrier numbers of MIMO-OFDM, the bandwidth of each subcarrier will be smaller, and therefore the system becomes more vulnerable to a loss in orthogonality caused by highly Doppler spread, leading to inter-carrier interference (ICI) [5]. Recent research [1,2,4,6] has focused on ICI suppression for enhanced performance with acceptable complexity. These mainly dealt with post-FFT processing methods, such as low complexity equalization. Banded minimum mean square error (BMMSE) equalization was used in [7] on the band the structure of channel frequency response. It was shown that the complexity increases linearly with the block length. However, the complexity of the BMMSE is lower than the conventional MMSE, which is proportional to the cube of the number of subcarriers [7].
In [6,[8][9][10], iterative detection and decoding techniques were introduced, including the block turbo equalizer based on the exchange of soft extrinsic information between Information 2021, 12, 469 2 of 14 MMSE equalization and a maximum a posteriori probability (MAP) decoder. The banded equalization with iterative data detection has a superior robustness against errors of channel estimation, which usually occurs in the UWA scenario.
The partial fast Fourier transform (PFFT) [1,11,12] decomposes a received signal into several segments using non-rectangular windows, followed by FFT demodulation and a weight compensation process on each segment. The working principle of the PFFT resides in the fact that the coherence time increases. This is achieved by increasing the bandwidth by splitting the signal in time domain into several segments [12]. A novel partial fractional Fourier (PFrFT) based on SISO-OFDM scenario was proposed in [13], with a band minimum mean square error (BMMSE) weight combining equalizer based on least square MINRES (LSMR) iterative algorithm. Simulation results resulted in significant BER performance improvements over traditional OFDM based methods and PFFT demodulation. In [14], the conventional FFT demodulation in the OFDM was replaced by the fractional Fourier transform (FrFT), which transforms the signal into an intermediate domain between time and frequency. Doubly selective channel response will be concentrated into a narrower band, allowing the ICI in adjacent subcarriers to concentrate around main diagonal of channel frequency matrix. The result presented in [14] showed a better performance with less complexity. FRFT could also be replaced by fractional cosine transform (FRCT), as described in [15].
In this paper, a novel method called MIMO-partial fractional Fourier transform (PFrFT) is presented in conjunction with banded MMSE equalization enhanced communications over known UWA channels. Banded MMSE equalization is used. At the analysis stage, a doubly selective UWA channel scenario is considered, extracting simulated results and showing the superior performance of the PFrFT approach compared with conventional PFFT and FrFT-OFDM ones. Its performance can be enhanced in the MIMO-OFDM scenario. Moreover, the performance obtained uses a low-cost equalizer based on the LDL H approach [16].
The remainder of the paper is organized as follows. Section 2 presents the proposed MIMO-PFrFT based algorithm that includes a MIMO-OFDM block, the partial fractional Fourier (PFrFT) demodulation, and a low complexity ICI equalization based on LDL H algorithm. A computational complexity analysis of the proposed system is also presented. In Section 3, simulation results and discussions are reported. Section 4 concludes the paper.
Notation: In this paper, transpose, conjugate, and conjugate transpose are denoted as [.] T , [.] * , and [.] H , respectively. diag{.} denotes a diagonal matrix produced by a vector [.] i,j extracts the ith row and jth column element from a matrix. Finally, < . > N is the modulo-N calculation.

PFrFT Based UWA Transceiver
The proposed MIMO-PFrFT-OFDM system is illustrated in Figure 1. Figure 1a shows the Transmitter while Figure 1b shows the Receiver. Both transmitter and receiver comprise several blocks. As seen in Figure 1a, in the transmitter serial input data is changed to a parallel stream which are transmitted through M transmitters, followed by mapping onto a constellation via quadrature phase shift keying (QPSK) modulation. Then, ∈ 1: ], the data vector at the subcarrier of transmitter, is transformed into chirp waveforms in the time domain using an inverse discrete fractional Fourier transform (IDFrFT) of order α. After inserting cyclic prefix (CP), the chirp waveforms are transmitted through the known UWA channel.
As shown in Figure 1b, after removing CP, the received data stream ) ∈ 1: ] at each of N receivers corrupted by the UWA channel with additive white Gaussian noise ) ∈ 1: ], are processed by a PFFT at each receiver, leading to B partial blocks. Subsequently, a novel factorization-based weight combining is carried out on each partial block _ , ∈ 1: ] at subcarrier, and converted to estimate the symbol stream _ using the discrete fractional Fourier transform (DFrFT) demodulation of order α 1 . The remainder of this section will describe the subcomponents in greater detail.

Discrete Fractional Fourier Transform (DFrFT)
In an OFDM scenario, DFrFT [17] is used to map the in-phase and quadrature (IQ) symbols of QPSK modulation data of different subcarriers to chirp waveforms. Figure 2 displays the chirp waveforms, spectral energy distribution and Wigner distribution of signals from the 1st and 20th subcarrier. As seen in Figure 1a, in the transmitter serial input data is changed to a parallel stream which are transmitted through M transmitters, followed by mapping onto a constellation via quadrature phase shift keying (QPSK) modulation. Then, d m k m ∈ [1 : M], the data vector at the k th subcarrier of m th transmitter, is transformed into chirp waveforms in the time domain using an inverse discrete fractional Fourier transform (IDFrFT) of order α. After inserting cyclic prefix (CP), the chirp waveforms are transmitted through the known UWA channel.
As shown in Figure 1b, after removing CP, the received data stream r n α (t) n ∈ [1 : N] at each of N receivers corrupted by the UWA channel with additive white Gaussian noise z n (t) n ∈ [1 : N], are processed by a PFFT at each receiver, leading to B partial blocks. Subsequently, a novel LDL H factorization-based weight combining is carried out on each partial block r n i_α (b), b ∈ [1 : B] at i th subcarrier, and converted to estimate the symbol streamd n i_α using the discrete fractional Fourier transform (DFrFT) demodulation of order (α − 1). The remainder of this section will describe the subcomponents in greater detail.

Discrete Fractional Fourier Transform (DFrFT)
In an OFDM scenario, DFrFT [17] is used to map the in-phase and quadrature (IQ) symbols of QPSK modulation data of different subcarriers to chirp waveforms. Figure 2 displays the chirp waveforms, spectral energy distribution and Wigner distribution of signals from the 1st and 20th subcarrier. A considerable amount of research has proposed different kinds to discrete FrFT (DFrFT), but many of them do not meet the requirements of all the required properties of continuous FrFT. The choice of DFrFT depends on the real application. In the work proposed in this paper, eigen-decomposition DFrFT [18] is applied due to its property of "orthogonality", "reversible" and "additivity", suitable for OFDM communication system. The DFrFT of order can be written in matrix vector multiplication as [19]  A considerable amount of research has proposed different kinds to discrete FrFT (DFrFT), but many of them do not meet the requirements of all the required properties of continuous FrFT. The choice of DFrFT depends on the real application. In the work proposed in this paper, eigen-decomposition DFrFT [18] is applied due to its property of Information 2021, 12, 469 5 of 14 "orthogonality", "reversible" and "additivity", suitable for OFDM communication system. The DFrFT of order α can be written in matrix vector multiplication as [19] r α = F α s (1) where the signal in DFrFT domain is is the N point DFrFT transformation matrix in which The DFrFT kernel F α (t, u) [19] consists of three parts, similar to the eigen-decomposition of DFT kernel. Both ϑ n [t] and ϑ n [u] are discrete Hermite-Gaussian functions, while e −j π 2 nα is the eigenvalue of DFrFT kernel.
The IDFrFT is written as According to [19], the eigen-decomposition type of DFrFT possesses the property of "additivity", so This property is used in the Partial FrFT element of the proposed system.

MIMO PFrFT-OFDM System
A DFrFT modulated MIMO-OFDM scenario with M transmitters, and N receivers is considered together with a Partial FFT demodulation at the receiver. The total block duration T = T + T g consists of OFDM symbol duration T and guard interval T g . For each transmitter m, QPSK data d m k modulated in each chirp subcarrier in the Fractional Fourier Domain of order α, where K is the total number of subcarriers IDFrFT used to convert d m k into LFM passband waveforms s m α (t) in the time domain where F −α (t, k) is the α order IDFrFT. After passing through the known time-varying channel, the received signal at the n th receiver, r n α (t), is expressed as [7] r n where z n (t) is the additive Gaussian noise and h n,m p (t) is the time varying amplitude for the p th channel of the (m, n) th transmitter-receiver pair. As shown in Figure 1b the duration of r n α (t), T, is separated into B non-overlapping intervals via non-overlapping rectangular window, followed by each segment of r n α (t) demodulated by FFT, contributing to output of b th interval of i th subcarrier written as is time varying channel frequency response in Equation (8), and z n i_α (b) denotes additive noise of b th interval of i th subcarrier.
Assuming that the UWA channel impulse response changes slowly over each segment, then the time variant channel frequency response over each partial FFT interval [ , bT B ] can be represented by its midpoint value, r n i_α (b) which can be written as where Equation (10) can be written in matrix -vector form, as (12) where The output of all K symbol subcarriers for at b th interval of n th receiver can be written as    T contains symbol of all the K subcarriers at m th receiver; and H n,m (b) is the channel frequency coupling matrix for all K symbol subcarriers at b th interval of (m, n) th transmitter-receiver pair which can be written as The majority of symbol energy for all K subcarriers is concentrated along diagonal of H n,m (b), while the ICI caused by time variant UWA channel which are distributed across neighboring subcarriers, remain to be mitigated using a successive weight combining What we need to do next is to compute w n i_α . As ICI at one subcarrier is mainly caused by subcarriers within its near neighborhood, a parameter D the "ICI span", is defined which means that each symbol only affects it D direct neighbours on each side, as illustrated in Figure 3.
contains main diagonals and two D off-diagonals (super and sub) extracted from , , while is the interference, including the ICI beyond the "ICI span" and additive noise. Then, the output vectors related to subcarrier r _ are extracted from equation (16), leading to _ , where _ = _ … _ … _ is a (2 1) size vector including not only output of subcarrier, but also its neighboring subcarriers for successive weight computation.
in equation (18)   The weight combining for the subcarrier at interval is given by where , … .
As becomes larger, the matrix inversion of , , leads to ill-conditioned matrix and a significant increase in computational complexity over a single tap weight combining. In order to solve the problem of matrix inversion in MMSE criterion an factorization is implemented Re-writing Equation (13) as follows: where contains main diagonals and two D off-diagonals (super and sub) extracted from H n,m (b), while H n r is the interference, including the ICI beyond the "ICI span" and additive noise. Then, the output vectors related to i th subcarrier r n i_α (b) are extracted from equation (16), leading to  square matrix with size of (2D + 1) × (4D + 1) extracted from H n,m d (b) related to r n i_α (b), shown below: (18) can be written as where h n,m The weight combining for the i th subcarrier at b th interval is given by Equation (23) is solved by forward substitution to obtain β n 2, i_α via the lower left triangular matrix L and a rescaling by the diagonal matrix D −1 to calculate β n 1, i_α . Finally, back substation with the upper right triangular L H yields β n i_α , which can be inserted into (20) in order to determined n i_α . After the computation ofd n i_α , the data detected is converted from frequency domain back to fractional domain via (α − 1)-order DFrFT, as followŝ

Selection of Optimal Fractional Order
As indicated in Equations (13)- (18), the majority energy of H n,m (b) is spread in "ICI span", which includes main diagonals and two D off-diagonals extracted from H n,m (b), denoted as H n,m d (b). The allocation of its power changes according to fractional order α of PFrFT-OFDM. When the system reaches its optimal fractional order α opt , the power in H n,m d (b) will reach maximum value, which means power beyond "ICI span" contributes to the least interference. Consequently, the search of optimal order α opt is based on exploiting carrier-to-interference ratio (CIR), defined as where ||.|| 2 F denotes the Frobenius Norm. The optimal order can be estimated as a maximum CIR problem, as follows The coarse to precise iterative search approach [20] can be applied within the range of α from −1 to1. It is obvious that the selection of optimal order α depends on the UWA channel properties, including number of subcarriers K, Doppler distortion, number of resolvable paths, and channel power delay profile.

Complexity Consideration
The computational cost analysis of whole system of both PFrFT-LDL H and G-PFFT [1] are provided in Table 1. Table 1. Complexity analysis of whole system.

Whole System
Computational Complexity The complexity of PFrFT-LDL H consists of three parts, including (i) calculation of LDL H based weight combining algorithms O(B(8D 2 + 22D + 4)K), (ii) computation of PFFT and (α − 1)-order DFrFT O(K 2 ), and (iii) B times of FFT demodulation O(BKlogK). As seen, there are three key parameters affecting the complexity, including the subcarrier number K, ICI span D and partial interval number B. As seen in Table 1 the complexity of the whole G-PFFT system consists of (i) calculation of weight combining algorithm O(KB 3 ) + O(K 2 B 2 ) and (ii) computation of B times of FFT demodulation O(BKlogK). It can be seen that the PFrFT-LDL H outperforms G-PFFT with a decrease in complexity. According to the weight combining part of the whole system, a comparative computational cost analysis per iteration between LDL H , BMMSE and the weight combining of G-PFFT system is provided in Table 2. Table 2. Complexity analysis of different weight combining algorithms.

Weight Combining Computational Complexity
It can be seen that the implementation of LDL H requires O(Q(8D 2 + 22D + 4)K) flops per iteration, which is less than the cost of BMMSE and general partial FFT. The complexity of weight combining algorithms of G-PFFT scenario rely on B times computation of one matrix inversion and two covariance matrices, with the cost of O(KB 3 ) + O(K 2 B 2 ) altogether.

Simulation Result and Discussion
In this section, the performance of the proposed algorithms are evaluated using simulations.

Parameter of MIMO-OFDM
We design a MIMO-OFDM system with subcarriers of K = 256, in which carrier frequency and baseband frequency is f c = 10 kHz and B w = 12 kHz respectively. The bandwidth of each subcarrier spacing is ∆f = B w K = 46.875 Hz, while the OFDM symbol duration is T = 1 ∆f = 21.33 ms. A cyclic prefix interval T g = T 8 = 2.67 ms is inserted to eliminate inter-symbol interference (ISI).

Parameter of UWA Channel
The channels of MIMO system are assumed independent. The multipath delay for each path is distributed exponentially duration of maximum delay spread approximately τ max = 2.5 ms. The average power reduces exponentially with delay, while the amplitudes are Rayleigh distributed. The normalized Doppler spread f d T = 0.1, 0.2, 0.3 and 0.4 respectively. The signal to noise ratio (SNR) per channel ranges from to 0 to 30 dB.  It is assumed that the time-variant channel information is perfect with 0.4 and number of partial segments are selected as B = 1, 2, and 4, respectively. The DFrFT optimal order is is searched in range of [0, 2] by fine-to-coarse method [20] and = 0.7, The ICI span for channel frequency matrix and LDL equalization is D = 3. It is seen that the performance of PFrFT-LDL is superior to that of the PFFT with a BER improvement of 13 dB for B = 1, 9 dB for B = 2, and 3 dB for B = 4. The reasons for positive results are attributed to increase of ICI span and better energy aggregation along the diagonal of channel frequency coupling matrix for DFrFT than that of conventional DFT.

Simulation Analysis
The channel frequency matrices for both PFFT and PFrFT are shown in Figure 5a,b respectively. In Figure 5 the channel frequency matrix H is scalar into a grey level image, the value of which is in the range of [0, 255]. The value "0" means black and "255" means white. The bands, which are diagonal lines in both images with white color represent the energy of H . The larger the energy, the closer the pixel value will be to "255". It is can be appreciated that the distribution of the band of PFFT distributes over more It is assumed that the time-variant channel information is perfect with f d T = 0.4 and number of partial segments are selected as B = 1, 2, and 4, respectively. The DFrFT optimal order is α is searched in range of [0, 2] by fine-to-coarse method [20] and α = 0.7, The ICI span for channel frequency matrix and LDL H equalization is D = 3. It is seen that the performance of PFrFT-LDL H is superior to that of the PFFT with a BER improvement of 13 dB for B = 1, 9 dB for B = 2, and 3 dB for B = 4. The reasons for positive results are attributed to increase of ICI span and better energy aggregation along the diagonal of channel frequency coupling matrix for DFrFT than that of conventional DFT.
The channel frequency matrices for both PFFT and PFrFT are shown in Figure 5a,b respectively. In Figure 5 the channel frequency matrix H df is scalar into a grey level image, the value of which is in the range of [0, 255]. The value "0" means black and "255" means white. The bands, which are diagonal lines in both images with white color represent the energy of H df . The larger the energy, the closer the pixel value will be to "255". It is can be appreciated that the distribution of the band of PFFT distributes over more frequencies is than that of PFrFT, which means the energy of H df based on PFrFT concentrate closer to the diagonal, leading to less ICI and better BER performance. span for channel frequency matrix and LDL equalization is D = 3. It is seen that the performance of PFrFT-LDL is superior to that of the PFFT with a BER improvement of 13 dB for B = 1, 9 dB for B = 2, and 3 dB for B = 4. The reasons for positive results are attributed to increase of ICI span and better energy aggregation along the diagonal of channel frequency coupling matrix for DFrFT than that of conventional DFT.
The channel frequency matrices for both PFFT and PFrFT are shown in Figure 5a,b respectively. In Figure 5 the channel frequency matrix H is scalar into a grey level image, the value of which is in the range of [0, 255]. The value "0" means black and "255" means white. The bands, which are diagonal lines in both images with white color represent the energy of H . The larger the energy, the closer the pixel value will be to "255". It is can be appreciated that the distribution of the band of PFFT distributes over more frequencies is than that of PFrFT, which means the energy of H based on PFrFT concentrate closer to the diagonal, leading to less ICI and better BER performance.   Figure 6 displays how the signal band to interference ratio (SBIR) changes with increase of ICI span D for partial segments B = 2. As D increases, the SBIR improves, but the rate of increase decreases as SNR increases. The reason is that symbol energy congregates among desired signal subcarriers, while the leakage energy of ICI spread into subcarriers alongside. As the index of ICI relative subcarriers increases, the normalized ICI power decrease. As ICI span D increases, more ICI in its neighboring subcarriers are considered for further LDL H equalization, contributing to decrease of residual energy out of ICI span, and finally increases BER performance.
Information 2021, 12, x FOR PEER REVIEW 11 of 13 Figure 6 displays how the signal band to interference ratio (SBIR) changes with increase of ICI span D for partial segments B = 2. As D increases, the SBIR improves, but the rate of increase decreases as SNR increases. The reason is that symbol energy congregates among desired signal subcarriers, while the leakage energy of ICI spread into subcarriers alongside. As the index of ICI relative subcarriers increases, the normalized ICI power decrease. As ICI span D increases, more ICI in its neighboring subcarriers are considered for further LDL equalization, contributing to decrease of residual energy out of ICI span, and finally increases BER performance. Comparing the performance of PFrFT-LDL in Figure 4 with B of 1,2, and 4, also confirms that an improved BER is achievable using larger values of B. Figure 7 displays the SIR of PFrFT in PFrFT with B from 0 to 16. The system suffers from the Doppler distortion from 0.1 to 0.4 and SNR is set at 30 dB. SIR grows with an increase of B, as the growth rate becomes slower. It is obvious that partial segments play a significant role in mitigating residual ICI out of span. Comparing the performance of PFrFT-LDL H in Figure 4 with B of 1,2, and 4, also confirms that an improved BER is achievable using larger values of B. Figure 7 displays the SIR of PFrFT in PFrFT with B from 0 to 16. The system suffers from the Doppler distortion from F d T = 0.1 to 0.4 and SNR is set at 30 dB. SIR grows with an increase of B, as the growth rate becomes slower. It is obvious that partial segments play a significant role in mitigating residual ICI out of span.
Comparing the performance of PFrFT-LDL in Figure 4 with B of 1,2, and 4, also confirms that an improved BER is achievable using larger values of B. Figure 7 displays the SIR of PFrFT in PFrFT with B from 0 to 16. The system suffers from the Doppler distortion from 0.1 to 0.4 and SNR is set at 30 dB. SIR grows with an increase of B, as the growth rate becomes slower. It is obvious that partial segments play a significant role in mitigating residual ICI out of span.    It can be seen that an improved BER is achieved in MIMO-OFDM compared to that in the SISO-OFDM due to diversity gains with cost of increase of computational complexity.  It can be seen that an improved BER is achieved in MIMO-OFDM compared to that in the SISO-OFDM due to diversity gains with cost of increase of computational complexity. DFrFT based serial Band MMSE equalization (D = 3) and conventional Partial FFT (B = 4) achieves a moderate performance and their performance are close to each other. The proposed PFrFT-LDL outperforms the other methods in SISO OFDM scenario and its performance can be improved in MIMO OFDM scenario. At SNR = 15 dB the performance of the B = 4 for MIMO system with N M 2 2 is 25 dB better than B = 4 and 35 dB better than PFFT B = 4 and DFRFT SBMMSE D = 3.

Conclusions
In this paper, a novel concept of PFrFT-MIMO-OFDM based on the hybrid use of the discrete fractional Fourier transform (DFrFT), partial Fourier transform (PFFT) demodulation, and low complexity banded MMSE equalization has been proposed. The LDL H based equalizer was applied not only to improve the performance, but also to reduce the computational complexity. The ICI is significantly mitigated under doubly selective UWA channel compared to the DFrFT-OFDM, PFFT-OFDM at moderate complexity providing an improvement in the overall Bit Error Rate.