Assessment of Wide-Sense Stationarity of an Underwater Acoustic Channel Based on a Pseudo-Random Binary Sequence Probe Signal

Featured Application: The method described can be used to assess whether the signal received in the underwater acoustic communication system can be modelled as a stochastic wide-sense stationary process. The fulﬁllment of this condition is necessary both to calculate the basic transmission parameters of the communication channel, such as coherence bandwidth and coherence time, and to correctly detect information transmitted in a UAC systems using modern modulation and coding schemes, such as orthogonal frequency-division multiplexing or spread spectrum techniques. A fast and energy efﬁcient method of assessing whether the received signal meets the WSS assumption can be implemented in a UAC adaptive system, enabling the change of modulation parameters, in particular, the duration of the modulation symbol for which the WSS assumption is fulﬁlled, thereby increasing the UAC system’s reliability. Abstract: The performances of Underwater Acoustic Communication (UAC) systems are strongly related to the speciﬁc propagation conditions of the underwater channel. Designing the physical layer of a reliable data transmission system requires a knowledge of channel characteristics in terms of the speciﬁc parameters of the stochastic model. The Wide-Sense Stationary Uncorrelated Scattering (WSSUS) assumption simpliﬁes the stochastic description of the channel, and thus the estimation of its transmission parameters. However, shallow underwater channels may not meet the WSSUS assumption. This paper proposes a method for testing the Wide-Sense Stationary (WSS) part of the WSSUS feature of a UAC channel on the basis of the complex envelope of a received probe Pseudo-Random Binary Sequence (PRBS) signal. Two correlation coefﬁcients are calculated that can be interpreted, together, as a measure that determines whether the channel is WSS or not. A similar wide-sense stationarity assessment can be performed on the basis of the Time-Varying Impulse Response (TVIR) of a UAC channel. However, the method proposed in this paper requires fewer computational operations in the receiver of a UAC system. PRBS signal transmission tests were conducted in the UAC channel simulator and in real conditions during an inland water experiment. The correlation coefﬁcient values obtained using the method based on the envelope of a probe signal and the method of analysing the TVIR estimates are compared. The results are similar, and thus, it is possible to assess if the UAC channel can be modelled as a WSS stochastic process without the need for TVIR estimation.


Introduction
Shallow underwater acoustic communication channels are characterised by disadvantageous transmission properties. Reflections from the sea-bottom and the water's surface cause multipath propagation, which goes hand-in-hand with the refraction phenomenon, caused by a significant change in sound velocity as a function of depth [1]. The movement of the UAC system's transmitter and receiver causes a Doppler effect, resulting in the time-domain scaling of a broadband communication signal. Moreover, UAC channel propagation conditions can change over time [2]. Depending on the phenomena under consideration, the variability of the transmission properties can be in the order of several months (i.e. seasons), several days and hours (i.e. tides, times of day), minutes (i.e. internal waves), a few seconds (i.e. surface waves), or milliseconds.
The transmission properties of a shallow UAC channel can be fully described by its time-varying impulse response, measured commonly by a correlation method, and modelled as a tapped delay line h(t, τ) [3]. A single realisation of h(t, τ) is defined in a window of observation time t and delay τ. The impulse response h(t, τ) can be treated as a two-dimensional stationary random process that fulfils the WSSUS assumption. Then, it is possible to calculate the 2-dimensional scattering function S(τ, ν) = ∆t ∆ f R h (∆t, ∆ f )e −j2π(ν∆t−τ∆ f )d∆td∆ f , where R h (∆τ, ∆ f ) is the space-time-frequency correlation function of the channel, obtained on the basis of h(t, τ). It is the basis for calculation of the transition parameters; namely, multipath delay spread τ M , Doppler spread ν M , coherence time T C , and coherence bandwidth B C [4].
As it is shown in [5], UAC channels are hardly ever WSSUS, and thus the transmission parameters can be determined only in a restricted period of time, and in a limited frequency range. Determining the time interval wherein the second-order statistics of impulse response are independent of time is important for the performance of communication systems working in nonstationary conditions. Only within this time interval is it possible to calculate transmission parameters {τ M , T C , ν M , B C }. Moreover, within this time interval, the signal received by the UAC system receiver represents a WSS process, and its power spectrum can be defined via a Fourier transform, according to the Wiener-Khinchin Theorem. This is especially important for signal detection performed in the frequency domain, as is the case with communication systems using the Orthogonal Frequency-Division Multiplexing (OFDM), and the Frequency Hopping Spread Spectrum (FHSS) techniques [4,6,7].
Studies on the validity of the WSS assumption are relatively rare. Several authors have proposed methods designed to test the similarity of the spectral properties of a time-series [8,9]. The relationship between the radiocommunication signal bandwidth and the wireless channel stationarity is described in [10]. The stationarity of a UAC channel was studied in [11]. The authors described the influence of the wind-driven surface wave variability on fluctuations of the channel Power Spectral Density (PSD). A trend-stationary channel model was proposed in [12] for underwater communication system. In [13], the stationarity of PSD is analysed to verify if the UAC channel can be considered to be a WSS process. Most of these methods are based on spectral analysis.
In [14], a novel time-domain method for verifying the WSS feature of UAC channels is described. Statistical properties of the real and imaginary parts of the impulse response, measured using the correlation method, are exploited to test the WSS assumption fulfilment. In [6] a similar method is proposed to assess whether the channel is WSS/non-WSS on the basis of the complex envelope of a received probe signal. However, the method was tested only in simulation conditions.
In this paper, both methods are compared using the results of simulations and the inland water experiment. The values of correlation coefficients, which are the basis of the WSS indicator proposed in [6], are presented. These coefficients express the similarity of the autocorrelation and cross-correlation functions of the in-phase and quadrature components of the received signal and the impulse response estimated using this signal. However, the calculation of the TVIR requires performing matched filtration in the receiver, which in the case of long impulse responses, measured using dozens of repetitions of a single probe sequence, is time and energy consuming. The simulation and experimental tests have confirmed that a decision can be made if the UAC channel fulfils the WSS assumption on the basis of the analysis of the complex envelope of the received probe signal, instead of the analysis of TVIR of the channel. Skipping the TVIR calculation stage can significantly reduce the processing time and power consumption of the UAC receiver.

UAC Impulse Response Measurement
The UAC channel impulse response measurements are usually performed by the correlation method with the use of probe signal having an impulse-like autocorrelation function [15]. One of these signals is a pseudo-random binary sequence (PRBS) being a repetition of a bipolar m-sequence c m ∈ {−1; 1}, modulated onto a binary phase-shift keyed waveform s(t) = Re[x(t)e j2π f 0 t ], and x(t) can be denoted as: where T = 1/B, B is the signal bandwidth, f c is the carrier frequency, and the sequence length M = 2 L − 1 is determined by its order L. The rect(t) denotes the rectangular function. The probe signal is constructed of K subsequent sequences s(t), each of duration T s , without any pause in between.
To produce signal s(t), m-sequence c m is upsampled by a factor R = f s /B. Next, the signal is passed through a binary phase shift keying (BPSK) block and it modulates the carrier frequency f c . At the receiver the analog signal from the ultrasonic transducer is digitised by an analog-to-digital converter. After demodulation, the signal is passed through a matched filter with coefficients corresponding to the upsampled m-sequence used in the transmitter. The matched filtration algorithm is usually implemented in a frequency domain, using a discrete Fourier transform (DFT) [? ]. As a result, at each observation time t, a single realisation of the time-varying impulse response h(t, τ) is obtained.

Testing Wide-Sense Stationary Assumption Fulfilment
A random process is wide-sense stationary if its mean value and its autocorrelation function are time invariant. The relation between the WSS bandpass process r(t), and its complex envelope y(t), has been described in [17] as: The complex envelope y(t) must be a zero-mean process in order that E[r(t)] be independent of t. Moreover, the autocorrelation function R rr (t + ∆t, t) = E[r(t + ∆t)r(t)] of bandpass signal r(t) is related to its complex envelope y(t) in the following manner: where If r(t) is WSS, then the t-dependent term in Equation (3) vanishes; i.e., R yy * (∆t) = 0. The process y(t) = I(t) + jQ(t) is a complex value; thus, R yy * (∆t) can be rewritten as: where are the autocorrelation functions of I(t) and Q(t), and R IQ (∆t) = E[I(t + ∆t)Q(t)] and R QI (∆t) = E[Q(t + ∆t)I(t)] are their cross-correlation functions. Thus, the condition for wide-sense stationarity of bandpass signal r(t) requires that: Moreover, the complex envelope y(t) should be a zero-mean process. According to [6], the similarity of autocorrelation functions R I I (∆t, τ) and R QQ (∆t, τ) can be expressed as the cross-correlation coefficient: and similarity of cross-correlation functions R IQ (∆t, τ) and R QI (∆t, τ), can be described as coefficient: In [14] the method of testing, if the impulse response measured with the probe signal s(t) satisfies the assumption of WSS, it is presented. It is based on testing the similarity of the autocorrelation functions: , defined similarly as in Equations (6) and (7) for in-phase h I (t, τ) and quadrature h Q (t, τ) components of the impulse response h(t, τ). The cross-correlation coefficients c h IQ and c h IQQI are calculated as functions of delay τ: For each TVIR, mean values of these coefficients can be obtained as: The calculation of these coefficients requires the calculation of the TVIR estimate first. As mentioned in Section 2, the DFT transform is usually used to implement the matched filtering algorithm. As it is shown in Figure 1, to calculate K realisations of h(t, τ), each of which has N samples, the N-point DFT algorithm should be run 2K times.To be more precise, one DFT is used for frequency domain transition, and one inverse DFT is needed to return to the time domain. Next, the autocorrelation functions (ACFs): R h I I and R h QQ , and cross-correlation functions (CCFs): R h IQ and R h QI are also calculated using DFT. There are 2N calls of K-point DFT in the case of calculation of ACFs, and 3N calls of DFT of the same size in the case of calculation of CCFs. The difference is due to the fact that in the case of ACF, the DFT is calculated for one input vector, and in the case of CCF, it is calculated for two different input vectors. In order to determine the correlation coefficients, 3N DFT calls are needed, each of the size 2K − 1, resulting from the size of input ACFs and CCFs. In the case of calculating the c IQ and c IQQI coefficients directly from the complex envelope y(t) of the received signal, the matched filtering step is omitted and there is 2K of N-point DFTs less to calculate.

Results
Simulation and experimental tests have been performed in order to compare the values of {c IQ , c IQQI } coefficients calculated on the basis of probe PRBS signal and {c h IQ ,c h IQQI } coefficients calculated on the basis of measured impulse response.

Simulation Tests
The simulation tests were performed using the Watermark simulator. It is a freely available benchmark for physical-layer schemes for underwater acoustic communications [18]. Its core is a replay channel simulator driven by at-sea measurements of the time-varying impulse response. Three of the communication channels available in Watermark were selected, represented by impulse responses measured in: Norway-Oslofjord (NOF1), Norway-Continental Shelf (NCS1), and Brest Commercial Harbor (BCH1). The NOF1 channel was measured at a range of 750 m and the water depth was 10 m. The probe signal frequency band was 10-18 kHz. The NCS1 channel was measured at a range of 540 m and the water depth was 80 m. The probe signal frequency band was the same as in the case of NOF1 channel measurements. Finally, the BCH1 channel was measured at a range of 800 m, where the water depth was 20 m. The frequency band of the probe signal was 32.5-37.5 kHz. The transmitter was bottom-mounted in the case of NOF1 and NCS1, and suspended in the water column in BCH1. The receiver was a bottom-mounted hydrophone (NOF1, NCS1) or a vertically suspended hydrophone array (BCH1). A more detailed description of Watermark channels can be found in [18]. Figure 2 shows the absolute values of exemplary impulse responses for each channel.
During the simulation tests performed, a PRBS signal was transmitted in the NOF1, NCS1, and BCH1 channels. For the construction of the PRBS signal, m-sequences of order 8 and 10 were used. The spectrum of the signal was in the 10-18 kHz band for the NOF1 and NCS1 channels, and 32.5-37.5 kHz for the BCH1 channel. The sampling frequency was equal to 200 kHz. Sixty transmission tests were carried out for each configuration of the signal and channel parameters. During a single transmission test, the probe sequence was repeated 32 times. On the basis of the received signals, after performing a matched filtration, TVIR estimates and the corresponding coefficientsc h IQ andc h IQQI were determined. Next, the complex envelope y(t) of each of the received signal was divided into segments of a duration corresponding to the duration of a single PRBS sequence. The coefficients c IQ and c IQQI were calculated as mean values over all segments of the probe signal. The minimum value of the coefficient is 0, and the maximum is 1. It is assumed that a value of the cross-coefficient above 0.5 means a strong similarity, and below 0.5 a weak similarity. This is consistent with the assumption of a stochastic model of the wireless communication channel, that two received signals are similar if their correlation coefficient is equal or greater than 0.5. This rule applies, inter alia, when calculating coherence time and the coherence bandwidth of the channel as a time interval and a frequency shift, respectively, for which the correlation coefficient of such delayed or frequency-shifted signals is equal or greater than 0.5 [19]. It is reasonable to assume that if values of both cross-coefficients: c IQ and c IQQI are greater than 0.5, the channel can be described as a WSS process.
In order to compare the coefficients obtained for the TVIR and for the complex envelope of the received signal, a euclidean distance measure d was calculated as: It has a meaning of error of estimation of c IQ and c IQQI parameters, assuming thatc h IQ andc h IQQI are the reference values. The distance d was averaged over 60 transmission tests for each channel. Table 1 presents the values of the coefficients calculated on the basis of TVIR, the coefficients obtained on the basis of the envelope of the PRBS signal, and values of distance measure d. It can be seen that the smallest differences between the correlation coefficientsc h IQ ,c h IQQI , and {c IQ , c IQQI } were observed for the NCS1 channel, and the largest differences were recorded for the BCH1 channel.

Underwater Experiment
The inland water channel measurements were conducted on the 4th and 5th of May 2017, in the Wdzydze Lake on the northern edge of the Bory Tucholskie forest complex (53 • 58 31 N 17 • 54 19 E). The bottom of the lake falls steeply into the depths of the water. In the shore zone, it is lined with gravel-stony material, and in the deeper parts, covered with a layer of mud. Wdzydze Lake is in many respects very similar to the Baltic Sea. A significant part of the lake area is more than 40 m deep, and the deepest area is more than 70 m deep. The shape of the bottom of the Wdzydze lake is more diverse than the bottom of the Baltic Sea; thus, the propagation conditions are even more difficult than those in the sea.
The weather was windy and it was raining on the 4th May, when the measurements were performed at a distance of 550 m. The next day, during the measurements at distances of 340 m and 1035 m, the weather was windless; it was not raining and the water surface was calm.
Each probe signal was generated using a computer with MATLAB software. Conversion from digital to analogue signal was performed using a digital-to-analogue converter from an NI USB-6363 device. Next, the signal was amplified and transmitted to water by an underwater telephone HTL-10 [20]. The transmission stand was placed on board the boat. The same equipment was used in the receiving stand, but was configured differently. The boat was not anchored, but drifted very slowly. Therefore, the Doppler effect had an effect on the transmitted signal, but to a much lesser extent than strong multipath propagation. The receiving equipment was placed in a measuring container whose position was fixed. The transmission transducer was sunk to a depth of 10 metres, regardless of the water depth of about 20-30 metres. The receiving transducer was sunk to a depth of 4 metres at a constant water depth of 7 metres; thus, even a slight rippling of the water had an effect on the stationarity of the received signal. A more detailed description of the experimental set-up can be found in [21]. Figure 3 shows the absolute values of exemplary impulse responses for each of measured channels. During the experiment, the PRBS probe signal was transmitted, based on m-sequences of order 8 and 10, repeated 500 times, and occupied frequency bands of 5 kHz, 8 kHz, and 10 kHz. The carrier frequency and the sampling frequency were equal to 30 and 200 kHz, respectively. Due to different propagation conditions and different transmission powers of the measuring system, the received signal to noise ratio (SNR) was different for each of the three measured channels. At a distance of 340 m, the SNR was in the range of 10 to 13 dB; at a distance of 550 m it was in the range of 13 to 18 dB; and at a distance of 1035 m the SNR was in the range of 32 to 35 dB.
On the basis of the received signals, similar to the simulation experiment, after performing a matched filtration, the impulse response estimates were determined. Next, the coefficientsc h IQ ,c h IQQI , were calculated. Simultaneously, the complex envelope y(t) of each of received signals was divided into segments of duration corresponding to the duration of a single PRBS sequence, and the coefficients c IQ and c IQQI were determined. Next, the values of the distance measure d according to Equation (10) were calculated.
The results of analysis of the transmitted PRBS signals and estimated TVIRs are shown in Table 2. As seen in Figure 3 there is a small Doppler shift observed in the impulse responses measured at distances of 340 and 550 m. The impulse responses were resampled to remove the influence of the Doppler effect. The cross correlation coefficients were calculated for both cases: before and after resampling. The values of the coefficients and the distance measure calculated after resampling the TVIR are written in brackets. It is clearly seen that a small Doppler shift due to the drifting of the boat with the transmitting stand does not significantly affect the obtained results. The values of cross-correlation coefficients {c h IQ ,c h IQQI }, and {c IQ , c IQQI } are also presented in Figure 4. As the measurement system range increases, the difference between the cross-correlation coefficients {c h IQ ,c h IQQI } and {c IQ , c IQQI } decreases. At a range of 340 m, the distance measure d is up to 0.28, while for a distance of 1035 m it does not exceed the value of 0.04. At the same time the highest correlation coefficient results were obtained for the channel at a distance of 1035 m. This channel turned out to be the strongest wide-sense stationary. The lowest coefficient results were obtained for the 550 m channel, in which transmission was carried out in bad weather conditions. The c IQQI values below 0.5 do not allow it to be classified as WSS. In order to estimate the impact of SNR (which was different during measurements at different distances) on the results obtained, additional simulation tests were performed using the signals measured at a distance of 1035 m (SNR was the largest in this case). Additive Gaussian noise of a different level was added to the signal, and then the correlation coefficients were calculated. The results are shown in Figure 5. It is clearly seen that if the SNR is equal or greater than 10 dB, the effect of additive noise on the calculation results is negligible.

Discussion
In the simulation tests carried out, the differences between the pairs of cross-correlation coefficients, {c h IQ ,c h IQQI } and {c IQ , c IQQI }, were different depending on the channel tested. The highest distance measurement values (up to 0.22) were obtained for the BCH1 channel, where the transmitting and receiving transducers were suspended in the water. In case of the two other tested channels, namely, NOF1 and NCS1, the transmitting and receiving transducers were bottom-mounted. The distance measure d for NOF1 and NCS1 is less than 0.1. Another difference between the BCH1 and NCS1/NOF1 channels is the fact that in the case of BCH1 channel the ratio of the bandwidth (5 kHz) to the carrier frequency (35 kHz) is 1/7, while in the case of the NCS1 and NOF1 channels, this ratio is equal to 8 kHz/14 kHz = 4/7. Thus, the BCH1 channel is much more narrowband than the two other simulated UAC channels. Assuming that values of cross-correlation coefficients c IQ and c IQQI greater than 0.5 indicate that the channel is WSS, using both pairs of coefficients leads to the same Watermark channels classification results as being WSS.
During the inland water experiment, the largest differences between the pairs of cross-correlation coefficients were obtained for the channel with the smallest range (340 m). In this case the distance measure d was up to 0.28. This channel is characterised by the strongest multipath propagation. In case of the longest range, namely, 1035 m, the distance measure between the coefficients was the lowest (up to 0.039). Both the coefficientsc h IQQI and c IQQI for the channel measured in bad atmospheric conditions over a distance of 550 m take values of less than 0.5. Thus, the channel was classified as non-WSS.
Assessment of channel wide-sense stationarity based on the coefficients c IQ and c IQQI , calculated directly from the complex envelope y(t) of the received signal, makes it possible to obtain the same classification result as using the coefficientsc h IQ andc h IQQI . At the same time, calculating c IQ and c IQQI requires fewer computational operations, due to the fact that the TVIR channel estimation step is omitted. Table 3 shows the sample calculation times and power consumption related to the estimation of the UAC channel TVIR by a fixed-point digital signal processor TMS320VC5505 by Texas Instruments [22]. The size of the impulse response corresponds to the results of measurements obtained during the inland-water experiment at a distance of 1 km. For each impulse response, the size of the complex FFT and the corresponding calculation time and the amount of energy consumed by the processor were determined. The calculation time for each estimate is about 1.5-2 s. In the case of very short range systems, it is longer than the time at which the acoustic wave travels the path between the transmitter and receiver. It may be an undesirable, significant delay in the case of an adaptive UAC system. The signal processor needs 63.21 to 83.10 mJ of energy to estimate the TVIR. If the calculation takes about 2 s, these values correspond to about 30-40 mW of power consumption. In case of the off-the-shelf acoustic modems, the power consumption on the receiving side usually does not exceed 1 W. Thus, every few dozen Watts of less power consumed might be a significant improvement of UAC system power efficiency.

Conclusions
An underwater acoustic communication channel can be classified as WSS or non-WSS on the basis of analysis of a baseband, complex value PRBS signal received by the UAC system receiver. If the received signal can be modelled as a WSS stochastic process, then its in-phase and quadrature components have similar autocorrelation functions, and its cross-correlation function is symmetric. The correlation coefficients of these components are being proposed as indicators for assessment of wide-sense stationarity of the UAC channel. Two algorithms of WSS assessment have been compared in this paper. One of them is used to calculate the correlation coefficients on the basis of the TVIR of the channel. The second one skips the TVIR estimation stage, and thus allows the channel to be evaluated with fewer computational operations. Analysis of the results of PRBS signal transmission in the UAC channel under simulation conditions and during an inland water experiment showed that although the two algorithms lead to slightly different values of correlation coefficients, the results of classification of the measured channel as WSS/non-WSS are the same in both cases.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: