A polynomial phase model for
estimation of underwater acoustic
channels using superimposed pilots
Felix Trulsson
Engineering Physics and Electrical Engineering, master's level 2019
Luleå University of Technology
A
BSTRACT
In underwater acoustic communications the time variation in the channel is a huge chal-lenge. The estimation of the impulse response at the receiver is crucial for the decoding of the signal to become accurate. One way is to transmit a superimposed pilot sequence along the unknown message, and by the knowledge of the sequence have the possibility to continuously track the variation in the channel over time.
This thesis investigates if it is possible by the aid of superimposed pilot sequences to separate the taps in the channel impulse response and using a parametric method to describe the taps as polynomial phase signals.
The method used for separation of the taps was a moving least squares estimator. Thereafter each tap was optimised to a polynomial phase signal (PPS) using a weighted non-linear least squares estimator. The non-linear parameters of the model was then determined with the Levenberg-Marquardt method. The performance of the method was evaluated both for simulated data as well as for data from field tests. The performance was determined by calculating the mean squared error (MSE) of the model over different frame lengths, signal to noise ratio (SNR), weights for the superimposed pilots, rapidness of time variation and impulse response lengths.
The method was not sensitive to the properties of the channel. Even though the model had high performance, the complexity of the computations generated long compilation times. Hence, the method needs further work before a real time implementation could be possible.
P
REFACE
I am grateful for the opportunity to perform the work of this thesis at FOI. I would personally like take the opportunity to give a special thanks to Dr. Magnus Lundberg Nordenvad for introducing me to the thesis topic as well as the overall interesting field of underwater communication and to Professor Jaap van de Beek for the interesting discussion the last few months. Both of which have been my supervisors during the work.
I would also like to thank Bernt Nilsson for helping me extract field test data and for the discussions of the content and Professor Johan Carlson for providing the LaTeX template.
The fantastic illustrations in Fig. (1.3) and on the cover deserves some extra attention, which all goes to Frida Georgsson.
Last but absolutely not least, a great thanks to my friends and family for all the support during my time at the university.
Felix Trulsson
C
ONTENTS
Chapter 1 – Introduction 1
1.1 Purpose and outline . . . 1
1.2 Introduction to digital communication . . . 2
1.3 Underwater acoustic communication channels . . . 4
1.3.1 Attenuation . . . 4
1.3.2 Noise . . . 4
1.3.3 Multipath . . . 5
1.3.4 Time variability . . . 5
1.3.5 The Doppler effect . . . 6
1.4 A review of underwater acoustic communications . . . 6
1.5 Contributions of the thesis . . . 9
Chapter 2 – Theory 11 2.1 Signal representation . . . 11
2.1.1 Superimposed pilots . . . 13
2.2 Channel representation . . . 15
2.2.1 Channel impulse response . . . 15
2.2.2 Polynomial phase signals . . . 15
2.3 Moving linear least squares estimator . . . 17
2.4 Weighted non-linear least squares estimator . . . 20
2.5 Exhaustive search method . . . 22
2.6 Levenberg-Marquardt method . . . 23 Chapter 3 – Method 27 3.1 Simulation setup . . . 27 3.1.1 Transmission . . . 27 3.1.2 Channel . . . 28 3.2 Experimental setup . . . 28
3.3 Separation of taps in time variant impulse response . . . 29
3.4 Parameter estimation . . . 29
Chapter 4 – Results 31 4.1 Simulation results . . . 31
Chapter 5 – Discussion 39 5.1 Future work . . . 41 5.2 Ethics . . . 42
Chapter 6 – Conclusion 43
C
HAPTER
1
Introduction
1.1
Purpose and outline
Underwater wireless communication has played an important role in commercial systems as well as in the military. It is necessary for national security and defence, the off-shore oil industry, pollution monitoring, collection of scientific data, disaster detection and early warnings, and also for the discovery of new resources. As the technology develops rapidly these fields grows explosively. Hence, the need for stable, high rate underwater communication is and will be a field that have to evolve [1]. This thesis serves as an attempt to accurately estimate all taps in the impulse response from the channel of an underwater acoustic communication channel.
The thesis is composed as follows. Chapter 1 first, apart from the purpose and outline, gives a review in the field of underwater communication, its applications and the proper-ties of the communication channel. It then introduces the concept of channel estimation and shortly describes different different method which have been used in present time. The last part of this chapter presents the scope of this thesis. It describes in which way the specific work presented fills a former gap in the field of underwater acoustic communication.
In Fig. (1.1) a visualisation of a communication system is presented which also con-cludes the work performed during the thesis. The theory describing this is presented in chapter 2. The chapter therefore begins by explaining how a signal s can be constructed before transmitted from the transmitter T x. Thereafter the convolution between the sig-nal s and the channel H is thoroughly described in section 2.2. Also how a time variant channel can be described using polynomial phase signals is presented in section 2.2.2. The channel H and the disturbance or noise e will distort the signal. Thus, when the receiver Rx samples the signal r as a distorted version of s an estimation of the channel H is needed in order to reverse the distortion. The last part of chapter 2 will therefore describe the theory of how an estimation of H can be performed where this estimation
2 Introduction
can be represented as ˆH.
Tx s H r Rx Hˆ
e
Figure 1.1: Visualisation of a communication system describing the thesis content.
The testing of the algorithms were in a large extent done by simulation. Therefore, the simulation of the systems transmitter and the propagation through its channel is described it chapter 3. Both of which refers back to the theory presented in the be-ginning of chapter 2. The second half of chapter 3 describes how the time variant taps of the channel impulse response were separated and how the full impulse response was estimated. This section connects the theory to the specific problem.
The results generated from performing the methods described in chapter 3 will be presented in chapter 4. Results both from simulated data as well as from data collected during field tests are presented.
Chapter 5 evaluates how well the method performs, which are its strengths and which improvements can be made. The chapter does also include a section describing future work suited as the next step regarding the method presented.
Chapter 6 will present the conclusion of the methods performance and in which state the method is when presented.
1.2
Introduction to digital communication
Telecommunication has been evolving since the discovery of electricity in the 1800’s. An example is the line telegraphy which was perfected by Morse in 1844 [2]. The mathe-matical theory of communication and its possibilities was published by Shannon in 1948 [3], where after the amount of applications started to arise [2]. The development of com-munication systems and the capacity restrictions in analogue systems made the need of digital communications during the later half of the 1900’s.
A digital communication system consists of a transmitting part, a channel and a receiv-ing part, roughly described. Both the transmitter and the receiver consists of multiple blocks, making communication through the system possible. A wireless communication system divided into the essential blocks is shown in Fig. (1.2).
1.2. Introduction to digital communication 3
Source Source Encoder Channel Encoder Modulator
channel Demodulator Channel Decoder Source Decoder Destination Transmitter Receiver
Figure 1.2: Block diagram of a wireless communication system [2].
techniques when transmitting analogue waveforms can be categorised into three main modulation methods [2]. These are phase-shift keying (PSK), amplitude-shift keying (ASK) and frequency-shift keying (FSK). Among these the modulation technique most often used in communication is PSK. An other method called QAM, which is a hybrid between PSK and ASK is also widely used in communications [2]. There are multiple variations within the modulation techniques giving a variety in robustness and complexity of the system.
From modulation, the symbols are still discrete, the transmitter therefore uses a mod-ulator translating the discrete symbols into continuous time analogue waveforms to be transmitted over the physical channel [4].
The physical characteristics of a channel varies widely as the channel may be a wire, a band of radio or audio frequencies, or a beam of light [3]. The characteristics of the channel distorts the transmitted wave. Communication over wireline can be well modulated as a linear time-invariant system [4]. The transfer function can then, by a feedback loop from the receiver, be assumed to be known at the transmitter. For wireless communications this is not the case. Due to mobility between the transmitter and receiver the channel may vary, leading to that an accurate channel feedback will be unavailable [4].
4 Introduction
1.3
Underwater acoustic communication channels
Due to the electromagnetic- and optical waves poor propagation properties underwater [5], acoustical waves are most often used for underwater communication. This even though the method suffers from limited bandwidth, the mediums refractive properties, Doppler-shift, rapid time variation, extended multipath and severe fading [6]. Thus, water is referred to as one of the most difficult medium working with, in use today [7].
One of the significant properties is the low propagation speed of approximately 1500 m/s. Since the wave also has multipath propagation, these properties lead to a delay spread in the order of tens or even hundreds of milliseconds. For an acoustical wave to propagate in water over long distance a low carrier frequency must be used. Therefore, even if the bandwidth is in the kHz region, no narrowband approximation can be done.
For a broader understanding the channel properties attenuation, noise, multipath, time variability and Doppler effect are described below in separate sections.
1.3.1
Attenuation
A signal transmitted in water will experience losses of different forms, one being the loss due to absorption. The absorption of the signal will affect the propagation range and is frequency dependent. Apart from the loss due to absorption the signal will suffer from spreading loss. A loss increasing by the propagation distance. The total loss can be modelled as a function of the signal frequency f and the distance of the transmission as,
A(l, f ) = l lr
k
a(f )l−lr, (1.1)
where lr represents some reference distance and models the loss due to spread and a(f )
the frequency dependent absorption coefficient[7].
1.3.2
Noise
1.3. Underwater acoustic communication channels 5
Figure 1.3: The multiple paths of transmitted signals towards the receiver.
1.3.3
Multipath
Two effects that will induce multipath formation of a transmitted signal is reflection and sound refraction in water. Reflection of the signal is due to reflection at the surface, at the bottom and also at different kind of objects such as large rocks, islands etc [7]. Refraction in the water is a result of the sound speeds spatial variability. The variability is a consequence of that the sound speed depends on the temperature, salinity and pressure. In shallow water these factors depends on environmental circumstances such as season, depth and location. Breaking waves among others will in many cases also induces multipath propagation due to scattering [8].
As is shown in Fig. (1.3) these different effects will affect a transmitted signal beam of rays, each ray might take a slightly different path to the receiver. The receiver then observes multiple signal arrivals. In accordance to Snell’s law a ray will always bend towards a region of lower propagation speed. Therefore, even though a ray has a longer propagation distance it might reach the receiver before one of shorter distance [7].
1.3.4
Time variability
6 Introduction
a very long time scale such as monthly changes in temperature, which does not affect the signal instantaneously, as well as waves breaking which will affect the signal due to scattering of the signal and Doppler spread. The Doppler spread is induced by the changing in path length [7].
1.3.5
The Doppler effect
The motion of transmitter or receiver does not only induce time variability, it also con-tributes to the change in channel response through the Doppler effect. It causes both frequency shifting and frequency spreading with a magnitude proportional to the ratio between the sound speed and the velocity the transmitter and receiver have in relation to each other. This effect is in radio- and electromagnetic communication negligible since the velocity of the signal is in a much greater order then any transmitter-receiver motion. In acoustic communication this is not the case. A small motion such as if a transmit-ter/receiver moves unintentionally, the relation to the very low speed of sound is not even then negligible. An underwater autonomous vehicle can move in the order of a few m/s which creates a severe distortion[7].
1.4
A review of underwater acoustic communications
Historical perspective
The upcoming of submarines started the development of underwater acoustic communi-cation as there was a need for communicommuni-cation. An underwater telephone using analogue modulation was developed in 1945 [9]. The phone was developed by the U.S Navy and was called ”Gertrude” [6]. Similar systems have been used even in present time [6].
The technique uses single side band modulation as carrier with analogue filters for pulse spectral shaping around the human voice band at the transmitter. The filtering, reproduction and demodulation performance is often of low accuracy. However, as the human mind is able in some extent to process distorted speech the system does work [6]. In the 1960’s the awareness of signalling and modulation in imperfect channels in-creased as the development of digital communication did as well. The throughput in these channels were limited leading to extreme low data rates. From then researchers have tried to increase the throughput for underwater propagation [6]. Hence, increasing the performance in digital underwater acoustic communication.
1.4. A review of underwater acoustic communications 7
Outline of underwater acoustic communication
The earlier sections stated the multiple challenges regarding communication under water. Due to these complications techniques developed for wired or wireless communications for terrestrial use does have to be significantly modified for underwater use [6].
Before transmission in water a robust system is crucial. In usage of coherent signalling PSK and QAM, mentioned in section 1.2, does for the medium provide a robust modu-lation [11].
To estimate the channel in communication, sequences known to both the transmitter and the receiver are used to estimate the change in the channel. These sequences are called pilot sequences. Two competitive symbol based schemes are conventional pilot sequences and superimposed pilot sequences [12]. The difference between them will be described thoroughly in the theory. A brief description, however, will be made here. The approach of conventional pilots is often used, meaning a sequence of know pilot symbols will be added both before and after a transmitted signal, so called preamble and post amble. This method restricts the length of a signal possible to decode in the receiver depending on how rapid the channel is shifting. Another concern is the waste of bandwidth during the pilot sequence as it bears no information [12]. When using superimposed pilots the preamble and post amble are still added. However, parallel with the unknown signal a known pilot sequence is transmitted. This makes it possible to continuously track change in the signal along the unknown symbol sequence [13]. Therefore, the throughput of the channel can be increased in comparison with other conventional pilot distributions as the length of the signal between preamble and post amble can be increased [12].
To estimate the channel for the signal to propagate in, one choice to be made is if a parametric or non-parametric estimation is to be used. The parametric estimation is defined as, to compute the estimate of a sample in a received symbol the gain and delay of different paths are estimated and used to describe the sample. A non-parametric estimator on the other hand estimates the sample directly from the received symbol and ignores the propagation path and the structure underlying the symbol [14].
Methods used for estimation
The section highlight major challenges in underwater acoustic communication and how researchers in present time have developed methods to avoid them.
In channels induced by high Doppler spread PSK modulation together with decision feedback equalisers (DFE) and spatial diversity have been a computationally complex but effective combination for communication [15].
Due to discrete arrivals of multipath spread have lead to lower complexity, enhanced performance and faster channel tracking. In these cases a sparse structure of the equaliser has been used [16], [17], [18].
equali-8 Introduction
sation. These do not require training sequences for the equaliser to converge. However, the method converges slower than conventional equalisers resulting in limitations in long or continuous data streams [19].
As the DFE due to the inaccurate decisions in the feedback loop the method suffers from error propagation. Therefore, to ensure low bit error rate (BER) forward error correction codes have been used by researchers while developing turbo equalisation techniques. The techniques have iteratively interacted between decoder and equaliser resulting in joint estimation, equalisation and decoding [20].
Maximum a posteriori probability (MAP) equalisers are equalisers of high computa-tional complexity often used in turbo code. Previous work have developed a soft input DFE for each receiver with a linear equaliser both lowering the complexity and still achieving high performance [21].
1.5. Contributions of the thesis 9
1.5
Contributions of the thesis
In this thesis a parametric method was used to estimate a time variant channel impulse response. The method is based on that the transmitted signal did include a superim-posed pilot sequence [12]. The prosuperim-posed model was to represent each tap in the impulse response by a polynomial phase signal (PPS) [24]. The PPS-model have previously been used in other areas [24],[25], [26].
The scope was to separate the taps in the impulse response by a moving least squares method [27]. This step resolves the multipath of the channel, see Fig. (1.4). When separated, each tap was estimated as a PPS with constant amplitude using the weighted non-linear least squares method [28]. The estimation describes the time variance of the channel, see Fig. (1.4). The constant amplitude was used due to the properties of the channel and the fact that a PSK modulation only depends on the phase and not the amplitude.
Hence, the following were investigated in the thesis:
• Is it possible to describe an underwater acoustical communication channel using polynomial phase signals?
• How well does the proposed method perform with different channel properties? • In what extent is the model affected when extending the message block length?
Multipath
Time variation
r
H
ˆ
1H
ˆ
2C
HAPTER
2
Theory
The theory chapter is divided into several sections starting with the basics of signal repre-sentation. It will then describe the different parts required to perform the work presented in the thesis. The order the sections are structured makes it possible to understand one part at a time and at the end have the full picture.
2.1
Signal representation
In communications, passband channels are employed, which implies that the transmitter and receiver must be able to handle passband signals. However, all information carried in a real passband signal is possible to represent by a complex baseband signal. The signal representation in a complex baseband is of profound practical significance since a markedly lower sampling rate can be used and still receive an accurate discrete time signal representation [4].
QPSK modulation
Wireless digital communication needs to modulate data bits, taking the values of 0 or 1, into an analogue waveform as mentioned in section 1.2. Modulation is used to create an alphabet which can represent digital information in the analogue wave [2]. With the alphabet a bit-to-symbol map can be created, were a symbol either just represents a 0 or a 1 alternative a sequence of bits. As stated in section 1.2 one well used form of modulation is PSK, a method which holds the modulus of symbols constant and only lets the phase vary, see Fig. (2.1). The figure describes a constellation called Quadratic PSK (QPSK), each symbol can then carry the information of two bits [4].
12 Theory Re Im 01 00 10 11 ◦ ◦ ◦ ◦ 1 √ 2 1 √ 2 −1 √ 2 −1 √ 2 − − − −
Figure 2.1: The QPSK-constellation in the complex plane.
Root-raised cosine
Transmitting and receiving a digital signal, intersymbol interference (ISI) is a factor that has to be minimised. One way to achieve this is by applying a matched filter in the transmitter and/or in the receiver. In communications one of the most frequently used is the raised cosine filter [4]. In practice the filter should be evenly divided between the transmitter and receiver to minimise the ISI [29]. The result is applying a root-raised cosine both in the transmitter and the same in the receiver [30].
Regarding the Nyquist criterion it applies to the transmitter, channel and receiver. Therefore the filter in the transmitter GT X(f ) and in the receiver GRX(f ) is constructed
such that the product GT X(f ) GRX(f ) is Nyquist in frequency domain [4]. As the raised
cosine filter GRC given as,
GRC(f ) = T |f| ≤ 1−a2T T 21 − sin((|f|) πT a ) 1−a 2T ≤ f ≤ 1+a 2T 0 otherwise, (2.1)
where a is the fractional excess bandwidth is Nyquist [4]. The root raised cosine is then given as,
GRRC(f ) =
p
GRC(f ), (2.2)
2.1. Signal representation 13
Upconversion and downconversion
After the signal is modulated and filtered in the complex baseband it must be upconverted to passband representation. The passband signal which is transmitted over the channel is described as a real sinusoid with a carrier frequency fc representing of the complex
baseband signal.
At the receiver the passbandsignal is downconverted from the real passband to the com-plex baseband. As the digital signal processing preferably is performed on the baseband representation of the signal, downconversion is the first step at the receiver minimising amount of analogue processing. Extraction of the complex signal components are done separately. For further explanation the reader is recomended to read [4].
2.1.1
Superimposed pilots
In radio communication, adding a preamble and a post amble is commonly used, where the preamble and post amble contains data symbols known to both the transmitter and the receiver, these are called pilot sequences. By this knowledge it is possible for the receiving side to decode the unknown symbols in between, see Fig. (2.2) for setup. The transmitted signal s[n] will then be formed as,
s[n] = ( p[n] 1≤ n ≤ N p d[n] Np+ 1≤ n ≤ N − Np p[n] N− Np + 1≤ n ≤ N, (2.3)
where N is the total number of samples, p[n] is a known pilot sequence with length Np
and d[n] is the unknown data of length N − 2Np.
As for the case of underwater acoustic communication the coherence time of the wave is too short for this method to be efficient. Meaning that, the channel the wave propagates through change too rapidly during the time it takes for the unknown symbols to reach the receiver for the preamble and post-amble to decode the received data. A method to handle this problem is then to use a small fraction of the energy used to transmit the unknown data to add a sequence of known pilots, see Fig. (2.3). This sequence is called, as the title reviles, superimposed pilots. Mathematically the transmitted signal s[n] is described as, s[n] = ( p[n] 1≤ n ≤ Np− Nk (1− α)d[n] + αp[n] Np− Nk+ 1≤ n ≤ N − Np+ Nk p[n] N − Np+ Nk+ 1≤ n ≤ N, (2.4)
14 Theory
If the energy is to be kept constant in the signal Nk is the number of samples which
the preamble and post-amble is shortened with. Nk can therefore be described as,
Nk = α2(N − 2N p) 2(1− α2) [13]. (2.5)
If there is no restrictions for the signal energy, Nk is set to zero.
Data Preamble Post-amble time N Np Np N - 2Np
Figure 2.2: Sequence of transmitted data using preamble and post-amble to decode at the receiver.
Superimposed pilots Preamble Post-amble Data time N Np- 2Nk Np - 2Nk N - 2Np+ 2Nk
2.2. Channel representation 15
2.2
Channel representation
2.2.1
Channel impulse response
A signal s[n] propagating through a channel will be distorted by the channel impulse response due to the physical properties described in section 1.3. If a channel with Additive White Gaussian Noise (AWGN) is considered the signal also becomes affected by the noise in the channel. Therefore, the signal from a such a channel can be defined as an distorted AWGN-signal and the received signal can be expressed as,
r[n] =
K−1
X
k=0
hk[n− k]s[n] + e[n]. (2.6) Where K is the length of the channel response, the notation k denotes a response from the channel at a specific delay and e[n] is the AWGN.
As Eq. (2.6) describes definition of a convolution between the signal and the channel it can be written in matrix form as,
r = Hs + e, (2.7)
where H is the N × N convolution matrix for the channel specified as,
H = h0[1] 0 0 . . . 0 h1[2] h0[2] 0 . . . 0 h2[3] h1[3] h0[3] . . . 0 .. . ... ... . . ... hK[K + 1] hK−1[K + 1] . . . h0[K + 1] . . . 0 .. . ... ... ... . . . ... 0 . . . hK[N ] . . . h0[N ]. (2.8)
The diagonals in the matrix H is then the time variant taps of the impulse response of the channel. Each tap in an impulse response represents a path taken by the transmitted signal.
2.2.2
Polynomial phase signals
16 Theory Re Im 01 00 10 11 ◦ ◦ ◦ ◦ 1 √ 2 1 √ 2 −1 √ 2 −1 √ 2 − − − − ◦ ◦ ◦
Figure 2.4: Showed in the figure is the effect on a symbol due to a time variant channel.
x[n] = γe(jθ[n]). (2.9) θ[n] is a polynomial of any given order and can be written accordingly,
θ[n] =
M −1
X
m=0
θm(n)m, (2.10)
where θm is the parameter to the variable of order m for m = 0, 1, 2, ...., M− 1
describ-ing the phase variation [24].
As described in section 1.3, the rays of a transmitted signal will due to multiple effects reach the receiver at different times. Hence, the receiver will sample a sum of multiple signals.Therefore, if h[n] from Eq. (2.6) is described as a polynomial phase signal. A received signal could then be represented as,
r[n] = K−1 X k=0 γke j PM −1 m=0θm,k(n−k)m s[n] + e[n] (2.11)
2.3. Moving linear least squares estimator 17
2.3
Moving linear least squares estimator
To resolve the multipath in the channel, described in section 1.3, the taps of the impulse response have to be separated. From Fig. (1.4) in section 1.5 this is the first step of the estimation according to the proposed method of the thesis. Hence, this section describes one approach to resolve the multipath.
Def 1: Linear least squares estimation
For a linear system Ax=b the least square estimator can be used in order to determine the vector x. The least square solution will then have an error e = b - Ax. Hence, the square error of the residuals are given as
||b − Ax||2. (2.12) This is, in mean square sense the best solution for a linear system. The error is orthogonal to the solution which gives,
AHe = 0 =⇒ AHAx = AHb. (2.13) The solution to the problem then becomes,
x = (AHA)−1AHb [32]. (2.14)
According to the definition the method described is only applicable in a linear case, whereas in a real system that is often not the case [28]. However, considering a non-linear function f (t) which are slowly varying. The function can then be estimated as a set of linear functions as is described in the definition of moving least squares estimation [27].
tp
f(t)
t
18 Theory
Def 2: Moving least squares estimation
Given a non-linear function f (t), an estimation can be written as a sum of linear functions as,
f (t) =
P
X
p=1
fp(tp), (2.15)
where each function is only active an interval tp, where tp is defined for p =
1, ..., P. For visual representation Fig. (2.5) shows an interval for where a function fp(tp) is active. Each of the functions fp(tp) can then be estimated
using the least square approximation [27].
Following the definition of the linear least squares method each sequential function fp(tp) can be solved as
fp(tp) = (AHp Ap)−1AHp bp, (2.16)
where bp is a N × 1 data vector, A a N × K system matrix and fp(t) the
K× 1 estimated function.
If combining Eq. (2.4) and Eq. (2.11) a system with multiple path propagation can be described. The channel impulse responses then distorts the signal as a non-linear complex exponential function. Then the moving linear least squares method could be used to estimate all channel components h(k)[n] in Eq. (2.11). The partial function fp
then represents the sum of h(k)[n] taps at a time instance p. If extracting a small fraction of the signal r[n] from Eq. (2.11), it can be considered constant. All K components h(k)[n] can then be extracted by constructing a convolution matrix P of an M sized window of the pilot sequence p[n]. P then becomes an M × K matrix for each time instance p, where K is the number of components in the channel impulse response and M the window over time where r[n] can be considered constant. For the system to have a solution the condition M > K must be fulfilled. From Eq. (2.4) and Eq. (2.16) estimating the channel response at time instance p then can be described as,
ˆ
h1,p = (PHp Pp)−1PHp rp (2.17)
where (.)H is the Hermitian transpose, ˆh
1,pis the estimation of all taps at time instance
p and rp is the M -sized window of the received signal around time instance p.
During the preamble and post amble of a signal the full signal is known to the receiver compared to the block in between where only the fraction containing the superimposed pilot sequence is known. Therefore, the estimation ˆH1[n] of H[n] is divided into three
2.3. Moving linear least squares estimator 19 ˆ h1,p = (PHp Pp)−1PHp pp+ (P H p Pp)−1PHp ep 1≤ p ≤ Np ˆ h1,p = (PHp Pp)−1PHp αpp + (PHp Pp)−1PpH(1− α)dp Np+ 1 ≤ p ≤ N − Np + (PHp Pp)−1PHp ep ˆ h1,p = (PHp Pp)−1PHp pp+ (P H p Pp)−1PHp ep N − Np+ 1≤ p ≤ N. (2.18)
The vectors ep, pp and dp are M -sized windows of e[n], p[n] and d[n] respectively
around time instance p. For the intervals 1 ≤ p ≤ Np and N − Np + 1 ≤ p ≤ N the
signal is fully known therefore Pp is constructed as the convolution matrix of an M -sized
window in p[n]. However for the interval Np+ 1≤ p ≤ N −Np only a fraction if the signal
is known, dependent on the weight α. The convolution matrix Pp is then constructed
from an M -sized window of αp[n]. according to Eq. (2.4).
The noise (PHp Pp)−1PHp epis considered Gaussian distributed and the unknown symbols
(PHp Pp)−1PHp dp together with the Gaussian noise is seen as noise of an other magnitude
(and maybe other distribution). Eq. 2.18 is therefore rewritten to,
ˆ h1,p = (PHp Pp)−1PHp pp + (PHp Pp)−1PHp e (1) p 1≤ p ≤ Np ˆ h1,p = (PHp Pp)−1PHp αpp+ (P H p Pp)−1PHp e (2) p Np+ 1≤ p ≤ N − Np ˆ h1,p = (PHp Pp)−1PHp pp + (P H p Pp)−1PHp e (1) p N − Np+ 1 ≤ p ≤ N (2.19)
where e(1)p and e(2)p are the two different noise components ep and ep + dp. The full
estimation ˆH1[n] is then described as,
ˆ H1[n] = P X p=0 ˆ h1,p = P X p=0 (PHp Pp)−1PHp rp (2.20)
where the columns in H1[n] are the K time variant taps of the impulse response. H1[n] can therefore be written as,
ˆ H1[n] = [ˆh(1)1 [n]T, ˆh (2) 1 [n]T, . . . , ˆh (K) 1 [n]T]. (2.21)
20 Theory
2.4
Weighted non-linear least squares estimator
The time variation of the channel, described in section 1.3, does separately affect each tap of the impulse response. Therefore an estimation over time for each taps is required in order to resolve the time variation of the full channel. The following section propose an estimation method performed on each tap that was separated in section 2.3 which will describe its variation over time. This procedure is illustrated in Fig. (2.6).
ˆ H1[n] W N LS1 W N LS2 W N LSK ˆ H2[n] ˆ h(1)1 [n] ˆ h(2)1 [n] ˆ h(K)1 [n] ˆ h(1)2 [n] ˆ h(2)2 [n] ˆ h(K)2 [n]
Figure 2.6: The figure shows how the taps separated in the moving least squares estimator are separately estimated using the weighted non-linear least squares (WNLS). These together then describes the full impulse response of the channel.
Given a system model ˆh(k)2 [n] and data ˆh(k)1 [n] from one tap. The model best describing the data is the one minimising the error between them. Therefore, an error function J can be set up as J = N −1 X n=0 (ˆh(k)1 [n]− ˆh(k)2 [n])2, (2.22)
were ˆh(k)2 [n] is a non-linear parametric model. To solve the problem one approach, if possible, is to separate the linear parameters from the non-linear parameters and therefore simplify the complexity of the problem [28]. The model can then be described as
ˆ
2.4. Weighted non-linear least squares estimator 21
model is to be described by Eq. (2.9) of order two the parts in Eq. (2.23) can be written as, θ =[θ1, θ2]T, (2.24) H(θ) =[c01, c1e j(θ1+θ2) , . . . , cN −1e j(θ1(N −1)+θ2(N −1) 2) ]T, (2.25) φ =Ae jθ0 , (2.26)
where in this case the linear parameters are constant. c that is element wise multiplied to the non-linear part of the model is given accordingly,
c =[c0, c1, . . . , cN −1]T, where
ci =1 for i = 0, . . . , Np and i = N− Np + 1, . . . , N
ci =α for i = Np+ 1, . . . , N − Np.
(2.27)
This follows the theory in section 2.1.1 as of how large part of the data consists of superimposed pilots.
In Eq. (2.22) no assumption of a statistic model or weight between samples have been made. If the reliability of one sample is greater than another, it would also be of higher importance. Therefore, specific weights for each sample can be added to Eq. (2.22). Giving the weighted model
J =
N −1
X
n=0
w[n](ˆh(k)1 [n]− ˆh(k)2 [n])2, (2.28) where w[n] is the weight for the model at time instance n. The time variant channel tap ˆh(k)1 [n] represented as a vector can be written as ˆh(k)1 [n] = h. Eq. (2.28) can then be written in matrix form as
J (θ, φ) = (h− H(θ)φ)HW(h− H(θ)φ), (2.29) where the linear and non-linear parts are separated. W is the weighting matrix for the model. Setting the derivative of the error function J (θ, φ) equal to zero, then φ can be estimated accordingly
ˆ
22 Theory
J (θ) = hH(W− WH(HHWH)−1HHW)h (2.31) Knowing the nonlinear parameters in θ the linear or constant parameters in φ are estimated by eliminating the non-linear part from the data in ˆh(k)1 [n]. These parameters are then given according to the following equations [28].
ˆ γ = 1 N N −1 X n=0 ˆ h(k)1 [n]e −j(ˆθ1n+ˆθ2n2) (2.32) ˆ θ0 = ∠ N −1 X n=0 ˆ h(k)1 [n]e −j(ˆθ1n+ˆθ2n2) , (2.33)
Maximum likelihood estimator
If the data in ˆh(k)1 [n] is Gaussian distributed then Eq. (2.31) is the estimation of the mean. Hence, giving the probability density function (pdf) [28],
G(θ, W) = 1 pπ||W||exp − J(θ) . (2.34)
Maximising Eq. (2.34) gives the maximum likelihood estimation of the pdf. The maximum of Eq. (2.34) is found when minimising J (θ). Hence, given that ˆh(k)1 [n] is Gaussian distributed and W being the covariance matrix of ˆh(k)1 [n], minimising Eq. (2.31) also gives the maximum likelihood estimation [33].
2.5
Exhaustive search method
When solving the minimisation problem of Eq. (2.31) the direction of arrival estima-tions are not always optimised simultaneously. This can be solved by using an uniform exhaustive search. If a K-dimension parametric search is considered as,
[ˆθ1, . . . ˆθK] = arg min
θ1,...θK{J(θ)}.
(2.35) The search range is set uniformly around between [−a, a], where a is a positive real value, for each search parameter. The spacing between grid points are set to a step size ∆, spanning the grid according to Fig. (2.5). At each grid point the function J (θ) is evaluated after which the point minimising J (θ) is selected [34]. The method therefore does find the global optima.
2.6. Levenberg-Marquardt method 23
θ2
θ1
∆1
∆2
Figure 2.7: The figure visualises the search grid for a two dimensional exhaustive search.
2.6
Levenberg-Marquardt method
If evaluating all points in the criterion function is not an option an optimisation algo-rithm handling the problem in a more effective manner is needed, where the Levenberg-Marquardt algorithm could be used in such a case.
The Levenberg-Marquardt algorithm is a nonlinear optimisation algorithm used for constrained nonlinear problems [35]. As described in section 2.4 the function to be minimised is the weighted non-linear least square,
min θ J (θ) =||J(θ)|| 2 = N −1 X i=0 (Ji(θ))2, (2.36)
where J (θ) is given according to Eq. (2.31). A problem suited to solve using the Levenberg-Marquardt method.
The Gauss-Newton method follow the same optimisation process as the Newton method with line search. The difference being that the Gauss-Newton method uses the convenient and often efficient approximation [35],
∇2J =∇JH
∇J, (2.37)
of the Hessian. ∇J being the Jacobian of the function J and (.)H the Hermitian
24 Theory
Def 3: The Gauss-Newton method
The function to be minimised in the Gauss-Newton method is,
J (θ) = 1 2 m X i=1 (ri(θ))2, (2.38)
where ri are the residuals. In difference to the Newton method which
uses the true Hessian the Gauss-Newton method uses the approximation ∇JH
∇J, giving the system,
(∇JHk∇Jk)pk=−∇J H
krk, (2.39)
to be solved for the step length pk, then minimising the function J (θ) by iterating,
θk+1 = θk+ pk (2.40)
The approximation can also be used to obtain the Levenberg-Marquardt method, but the method replaces the line search with a trust region strategy. A trust region strategy avoids the weakness of the Gauss-Newton method, namely, the behaviour of the algorithm as the Jacobian ∇J becomes rank-deficient or nearly so [35].
Def 4: The trust region approach
In help of the a quadratic model the trust region method generates a step for its objective function. A region is defined around the current iterate for which the objective function is adequately represented within. The step provides the approximate minimisation of the objective function in the region. If the step is not acceptable the size of the region is reduced and a minimiser of the new area is approximated. The direction of the step is in general changed as the region size is.
From the trust region, for each iteration, the subproblem to be solved is, min dk 1 2||∇Jkdk+ rk|| 2, subject to ||dk|| ≤ ∆k, (2.41)
dk being the step length, rk the residual and ∆k the trust region radius. If the step
length for the Gauss-Newton method lies strictly within the trust region, meaning dk <
∆k. The Levenberg-Marquard can be collapsed to the Gauss-Newton. If not, there is a
λ > 0 such that ||dk|| = ∆k. Then, from Eq. (2.41), the problem to be solved becomes,
(∇J(θk)H∇J(θk) + λkI)dk =−∇J(θk)Hrk. (2.42)
2.6. Levenberg-Marquardt method 25
minimised. Therefore, dk is extracted from the equation and then becomes,
dk =− ∇J(θk)H∇J(θk) + λkI −1 ∇J(θ)HJ (θ) (2.43) If J (θ + dk) < J (θ) then λk+1 = λk and the parameters of θ is updated as,
θk+1= θk+ dk. (2.44)
If J (θ + dk) < J (θ) then λk+1 = λ10k before reevaluating Eq. (2.43). This procedure is
C
HAPTER
3
Method
The work that generated the thesis was performed both on simulated data as well as data collected from field tests performed by FOI. Therefore, methods used to perform cases are described in this chapter. First the parts only performed in the simulated case is described in the section Simulation setup, including the transmission and the simulated channel with AWGN. Where after the section Experimental setup gives a description of the work performed on the field test data, leading up to the data used for the method. As the channel estimation and the parameter estimation was performed both on the simulated data and the real data, both is described in the corresponding sections. The parts described in this chapter were all implemented in separate functions using MATLAB 2019a.
3.1
Simulation setup
3.1.1
Transmission
In line with the theory in section 2.1 a N − 2Np-sample random signal was modulated.
In the same way was a N -sample pilot sequence generated. It made it possible for the signal and pilot sequence to be sent alongside one another according to the theory of superimposed pilots in section 2.1.1. In the transmission process a preamble was placed in the beginning of the signal and a post amble in the end. In between the unknown signal was added along side the superimposed pilot sequence which proportion of the signal energy was varied from one transmission to another. In practice it meant varying weight of α from Eq. (2.4). As the fraction of the signal energy used by the superimposed varied as did the length of the preamble and the post amble. It varied in accordance to Eq. (2.5).
28 Method
3.1.2
Channel
The channel block simulated how the channel affected the transmitted signal, as described in section 2.2.1. In the simulation each impulse response of the channel was constructed as a polynomial phase signal with constant amplitude according to section 2.2.2 and Eq. (2.9). The parameters of the polynomial phased signal for every impulse response was randomly generated inside a constrained range. The channel should vary more than π over a message frame in order for the method to be evaluated, but it was restricted by the assumption that the channel varied slowly. Otherwise the theory in section 2.3 does not apply.
A convolution matrix for the channel response was constructed according to Eq. (2.8). By multiplying the convolution matrix with the signal the channel distorted the signal s[n] according to Eq. (2.7). As an AWGN was assumed complex valued Gaussian noise was added to the distorted signal. The output from the channel was therefore a sum of distorted noisy signals. As the channel response was generated as a polynomial phase signal, the output followed Eq. (2.11).
3.2
Experimental setup
This section describes the transmission of messages during field tests earlier performed by FOI. It also describes the procedure in the receiver done before the channel estimation was performed. This was done using programs previously written by employees at FOI. Before the data from the field test was transmitted all information was saved in order to determine if messages was received.
3.3. Separation of taps in time variant impulse response 29
3.3
Separation of taps in time variant impulse response
According to the assumption described in section 1.5, the channel varied slowly in time and it was assumed that the coherence time in a small region around a single time instance was long enough for the signal to be considered constant. Therefore, a least squares approximation could be performed for a small interval around a time instance n. The performance was divided into three sections, the first during the preamble, the second during the superimposed sequence and the third during the post amble. The method followed the theory in section 2.3 in general and particularly Eq. (2.19). As the length of the channel response in this step was known, so was the dimension K. As one estimation was performed for each time instance the only unknown dimension was therefore the length of the window M , which was set to a value larger then two times K. This created an over determined system and each value ˆh(k)1 [n] was calculated as a mean of at least two received samples. The first section, during the preamble, an M× K convolution matrix was constructed using the MATLAB function convmtx() for each time instance in the preamble.
Solving Eq. (2.19), the response for each channel tap was separated and estimated at the specific time instance n, according to Fig. (2.5). Therefore, the output from each time instance was a K× 1 vector containing the value of the K channel taps at time n. The second section, during the superimposed sequence, only a small weight α of the received signal was known. The convolution matrix P therefore was weighted with α. When weighted, the procedure followed exactly the one of the first section. The second section produced a noisy estimation as the method had multiple noise components in accordance with Eq. (2.19).
As for the third section, during the post amble, exactly as in the first section the full signal was known. Hence, the procedure was the same for the first and the third section.
When separated, a parametric solution could be performed for one tap at a time.
3.4
Parameter estimation
The output from section 3.3 was a matrix of noisy estimations of the time-varying impulse responses. These estimations were to be estimated to a parametric model. The model used was that of polynomial phase signals with constant amplitude according to section 2.2.2 and Eq. (2.9). The estimation method used to match the model to the data was the weighted non-linear least squares estimator described in section 2.4. The optimisation method used to determine the parameters was the Levenberg-Marquardt method from section 2.6.
30 Method
in Eq. (2.31). The minimisation problem was then only dependent of the non-linear parameters in θ. These were estimated using the Levenberg-Marquardt method accord-ing to section 2.6. The Jacobian of the minimisation problem was numerically estimated using symmetric difference quotient for each of the non-linear parameters. A good initial guess was needed for the optimisation algorithm to find the global optima. Therefore, Eq. (2.31) was evaluated in multiple points, making an exhaustive search with a sparse grid. The minimum value from the grid search was set as the initial guess for the op-timisation algorithm. Then by using the MATLAB built in function lsqnonlin() the Levenberg-Marquardt method was used to find the optima of the error function. Know-ing the non-linear parameters, the amplitude and phase were estimated usKnow-ing Eq. (2.32) and Eq. (2.33) respectively.
C
HAPTER
4
Results
This chapter contains the collected results both from simulations and from experimental data. The results from the methods described in chapter 3 are presented. The method was in large extent tested and evaluated using simulated channels before testing on experimental data. Therefore, the chapter starts by presenting the results from the simulation and the methods performance as different parameters of the system varied.
4.1
Simulation results
Simulating the system made it possible to evaluate the performance of the proposed model in multiple aspects. The aspects studied was variance in the length of the message frame, the weights of the superimposed pilot sequence, in the length of the channel response and the rate of phase change over time.
Evaluation of the method
For accuracy validation of the method used in the thesis, it was evaluated by comparison to results known to be exact. The moving least squares estimator were compared to estimate the pilot sequence accurately. Another aspect for the moving least squares estimator was if it estimated the channel well enough for the weighted non-linear least squares estimator to be able to make an accurate estimation.
Fig. (4.1) describes how the moving least squares estimator makes an estimation of one tap in the impulse response. It also visualises how well the weighted non-linear least squares estimator could track the channel impulse response in the noise from the middle of the moving least squares estimation. Both being compared to the true time variety of the tap.
The Levenberg-Marquard method was compared to a grid search of the weighted least square function. This to evaluate the accuracy of the method as an exhaustive search by definition finds the global optima, as described in section 2.5. The result is shown in
32 Results 0 500 1000 1500 2000 2500 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 True channel MLS WNLS
Figure 4.1: The figure shows a comparison of the change in the channel after estimation with the moving least squares and the weighted least squares to the actual channel impulse response.
-60 -40 -20 0 20 40 60 10-2 100 102 104 106 Levenberg-Marquard Exhaustive Search
Figure 4.2: The figure shows the MSE dependent on SNR for the Levenberg-Marquard method and an exhaustive search. The SNR is defined in dB and the MSE is shown in log10-scale.
Fig. (4.2).
When generating the results the following data was used: frame length N = 2500 symbols, channel length K = 11 taps, length of preamble and post amble Np = 250
symbols, the weight of the superimposed sequence was α = 0.4 and the channel varied a maximum of 4π radians over one frame.
A visualisation of the area created by the 2 parametric function J (θ) from Eq. (2.31) is shown in Fig. (4.3). The global optima in this case is placed at coordinates (θ1, θ2) =
4.1. Simulation results 33 -5 500 4 1000 1500 10-7 2 0 2000 10-3 2500 0 -2 5 -4
Figure 4.3: The image shows the surface of the parametric search area for θ.
Superimposed pilot weights
For evaluation of performance using different weights for the superimposed pilots, α was varied. The values assigned were [1, 0.4, 0.25, 0.1], where the value 1 corresponded to the full signal energy, Hence, in that case the whole signal was known. The result from the simulation is presented in Fig. (4.4).
-30 -20 -10 0 10 20 30 10-5 100 105 = 1 = 0.4 = 0.25 = 0.1
Figure 4.4: The figure shows the MSE depending on the SNR when the weight of the superim-posed pilots varied. The SNR is defined in dB and the MSE is shown in log10-scale.
34 Results
symbols and the channel varied a maximum of 4π radians over one frame.
Frame length
The result when the length of the message frame is presented in Fig. (4.5). The figure visualises the performance of the algorithm while the length of the signal frame was N = [900, 2050, 3200, 4350, 5500] symbols. -30 -20 -10 0 10 20 30 10-2 100 102 104 106 N = 900 N = 2050 N = 3200 N = 4350 N = 5500
Figure 4.5: The figure shows the MSE depending on the SNR when varying the frequency rate in the channel. The SNR is defined in dB and the MSE is shown in log10-scale.
When generating the results the following data was used: the weights of the superim-posed sequence α = 0.4, channel length K = 11 taps, length of preamble and post amble Np = 250 symbols and the channel varied a maximum of 4π radians over one frame.
Time variance
The performance when changing how rapid the phase varied in the channel over time is presented in Fig. (4.6). The maximum phase change over one frame was [π, 2π, 3π, 4π, 5π] radians.
When generating the results the following data was used: frame length N = 2500 symbols, channel length K = 11 taps, length of preamble and post amble Np = 250
symbols and the weights of the superimposed sequence α = 0.4.
Channel response length
The channel response length, as dependent on the channel properties, was a variable for evaluation. The number of taps were set to the values [5, 8, 11, 14, 17, 20]. The perfor-mance for different SNR is presented in Fig. (4.7).
4.2. Experimental results 35 -30 -20 -10 0 10 20 30 10-2 100 102 104 106 2 3 4 5
Figure 4.6: The figure shows the MSE depending on the SNR when varying the rate of the time variation in the channel. The SNR is defined in dB and the MSE is shown in log10-scale.
-30 -20 -10 0 10 20 30 10-2 100 102 104 106 K = 5 K = 8 K = 11 K = 14 K = 17 K = 20
Figure 4.7: The figure shows the MSE depending on the SNR when varying the length of the channel. The SNR is defined in dB and the MSE is shown in log10-scale.
and post amble Np = 250 symbols and the weights of the superimposed sequence α = 0.4.
4.2
Experimental results
The results presented in this section is performed on data received from FOI. The data was generated and collected during former field tests. The pre-processing of the data, following the method in section 3.2, was performed using programs written at FOI in MATLAB.
Frame lengths
Fig. (4.8) presents the result from estimating the channel during a field test earlier performed by FOI. The weight used for the superimposed pilot sequence was α = 0.2.
36 Results 0 2000 4000 6000 8000 0.016 0.018 0.02 0.022 0.024 0.026 0.028 0.03
Figure 4.8: The figure presents the MSE of the estimation depending on the length of the message frame using α = 0.2.
The length of the channel impulse response in this case was K = 37.
During the field tests which results are presented in Fig. (4.9) and Fig. (4.10), data was transmitted using the weights α = 0.4 and α = 0.25 for the superimposed pilot sequences. The transmissions were performed during the same field test, therefore also under similar conditions.
When performing the channel estimation presented in Fig. (4.9) was, the block lengths of the message was N = [1010, 2805, 4600, 6395, 8190] symbols, the channel response length K = 21 taps, the preamble length Np = 255 symbols.
4.2. Experimental results 37 0 2000 4000 6000 8000 0.006 0.008 0.01 0.012 0.014 0.016 0.018
Figure 4.9: The figure presents the MSE of the estimation depending on the length of the message frame using α = 0.4.
0 2000 4000 6000 8000 0.01
0.015 0.02 0.025
C
HAPTER
5
Discussion
This chapter will evaluate the work throughout the thesis. This includes the assumptions from section 1.5, how well these held up and the performance of the proposed method. It will discuss if the model could be useful, which parts needed improvements and the appropriate next steps if the work does proceed.
Performance
In previous chapter it was shown the moving least squares estimator was able to separate the taps in the impulse response and weighted non-linear least squares estimator did find a solution comparing to the exact solution, as shown in Fig. (4.1). It validated the combination of the moving least squares method and the weighted non-linear least squares did have the ability to separate the taps and estimate each tap. Together describing the full impulse response.
The weighted non-linear least squares estimator also found the optimal solution for different SNR comparing to the exhaustive search when evaluating Fig. (4.2). The exhaustive search was guarantied, using a fine grid, to find the global optima since it evaluates all points in the region of interest. Hence, the weighted non-linear least squares estimator was validated to perform an accurate estimation of the defined problem.
The frame length of the unknown symbol was varied in both the case of the simulated data and the field test data. From these, Fig. (4.5), (4.8), (4.9), (4.10) determines the method does perform well for different frame lengths. Fig. (4.5) also validates the method for different SNR. In Fig. (4.9) and (4.10) there is not much change in the MSE, to tell from the magnitude of the MSE this could because the channel is stable for all the lengths evaluated. The difference to be of notice is between the different weights of the superimposed pilots. There can be seen that the weight does have an impact towards generating an even better estimation. Fig. (4.8), which had a wider spread in the impulse response was also harder to resolve using the method. This leading to a larger MSE compared to both the simulations and the other field test results. However,
40 Discussion
the pattern in Fig. (4.8) did show expected results as the channel should be harder to estimate with a longer frame as the time between fully known symbols did increase. The accurate estimation even when using long message frames was of great interest since if the gap between training sequences can be increased so can the data rate.
The weight of the superimposed pilot sequence are of greater importance in lower SNR, shown in Fig. (4.4). In high SNR a lower α would be preferred as the SNR for the unknown symbols then becomes higher. Therefore, for channels of low SNR the weight should be increased as it is better to be able to perform an estimation with low SNR for the decoding then to not be able to perform one at all.
The length of the impulse response does affect the result, seen in Fig. (4.7). This since the longer the length of the impulse response the closer the window M used in the moving least squares estimator comes to the coherence length of the channel. If a longer window then the coherence is needed there is, at the time of the writing, no method to resolve such a channel [36].
In the same way as a longer channel response creates a channel more difficult to decode so does a more rapid time variation. Due to the same reason. The result shown in Fig. (4.6) shows that the assumption of being able to estimate the impulse response in a time variant channel holds for a maximum phase change of 5π radians over one frame. This variation is to large for conventional pilots to unwrap and the change can still be considered slow. Hence, the coherence time is long enough for any estimator at all to work.
Together Fig. (4.2)-(4.10) validates the proposed method and its possibilities. This leads to the possibility to increase the frame lengths during transmission even in a highly varying channel, leading to higher bit rate though the channel.
Even though the method performed well for the different SNR tested this factor seemed do represent the larger impact to the performance. A large SNR is of course highly valued and not always possible, this is however using the proposed method the largest factor of performance. From evaluating Fig. (4.3), low SNR could result in a local optima becoming the global one if unfortunate, resulting in an optimised around the wrong parametric values.
Non-linear parametric search
5.1. Future work 41
Superimposed weights
It will be harder to perform the estimation if smaller weights are used for the superim-posed pilots. Noisier and therefore both larger and more local optima will be induced, see Fig. (4.3), which makes it more likely to find a local optima then the global one. On the other hand, if it is possible to perform an accurate estimation with small weights the SNR of the transmission will be much better. This since a larger part of the signal energy is used for the unknown symbols instead of a known sequence.
Computational expense
As shown in Fig. (4.3) multiple optima were present. Therefore, the initial guess was of greatest importance. However, the evaluation of the criterion function J (θ) was compu-tationally expensive and grid search needed to evaluate this function multiple times.
As Eq. (2.31) performs multiple operations in the order of N2 each evaluation is of high cost. As the industry wants to increase the throughput in communication, the length of message frames shall increase resulting in heavier computations. Therefore, the number of evaluations needs to be kept to a minimum for the method to be efficient. In an exhaustive search the evaluation is performed for each grid point. Therefore, searching for an optimal initial guess is computational expensive. As the taps of the response have travelled in the same channel, there might be a chance they are correlated and/or dependent in some extent. If so is the case the search for the initial guess of the strongest tap in the impulse response could be used for the rest, decreasing the number of evaluations of Eq. (2.31).
As mentioned in previous work [16] impulse responses in underwater acoustic com-munication, even though they are long, can often be treated as sparse. Because of the computational expense, if considering the impulse response as sparse it would result in fewer parametric searches and therefore lower the computational expenses.
5.1
Future work
The thesis and the work it included have shown it to be possible to represent the impulse response of a underwater acoustic channel using polynomial phase signals. However, further studies are still needed investigating if the method could be used for real time implementation. The following describes potential improvements to the method.
Analytical determination
42 Discussion
CRLB, if it exists, it would define how well the model could perform. Therefore, this needs to be derived to determine how close to the lower bound the present performance is or how much it can be improved.
Algorithm improvements
Investigation of the dependencies between the taps in the impulse response should be performed. As the taps does propagate in the same channel arriving close in time, some dependence between then might be found. If there is any dependence the initial value grid search for each tap could be weighted between the taps and the number of grid points could then be lowered. This could lower the computational expense of finding a initial value. The impulse response is also likely to be of a sparse structure [16]. One plausible way to make the model more adaptive to the specific channel is to only consider the taps that actually does carry information which would lower the number of parameter estimations needed.
5.2
Ethics
C
HAPTER
6
Conclusion
This thesis has examined the possibility to estimate the impulse response for an underwa-ter acoustic communication channel using polynomial phase signals and estimating the parameters of these signals. The approach was to separate the time variant taps of the impulse response using a moving least squares estimator. Each tap was then estimated to a polynomial phased signal with constant amplitude using the weighted non-linear least squares estimator. The Levenberg-Marquardt method was used for the weighted non-linear least squares estimates.
The thesis shows it to be possible to describe the change in at least some underwater acoustic channels over time using polynomial phase signals with constant amplitude. The method is not markedly sensitive to the frame length of the message, the length of the channel impulse response or a rapid time variance even in low SNR. Overall it could be concluded the method did perform a accurate estimation of an underwater acoustic communication channel. However, due to the complexity of multi-parameter optimisation it is of need of further work before it is ready for real time implementation.
R
EFERENCES
[1] X. Zhang, J. Cui, S. Das, M. Gerla, and M. Chitre, “Underwater wireless com-munications and networks: theory and application: Part 1 [guest editorial],” IEEE Communications Magazine, vol. 53, pp. 40–41, Nov. 2015.
[2] A. Grami, Introduction to Digital Communications. Elsevier Inc., 2016.
[3] C. E. Shannon, “A mathematical theory of communication,” The Bell System Tech-nical Journal, vol. 27, pp. 379–423, Oct. 1948.
[4] U. Madhow, Fundamentals of digital communication. Cambridge University press, 2008.
[5] G. Qiao, Z. Babar, L. Ma, S. Liu, and J. Wu, “Mimo-ofdm underwater acoustic communication systems: A review,” Physical Communication, vol. 23, pp. 56–64, June 2017.
[6] M. Chitre, S. Shahabudeen, L. Freitag, and M. Stojanovic, “Recent advances in underwater acoustic communications & networking,” Oceans 2008, pp. 1–10, Sept. 2008.
[7] M. Stojanovic and J. Preisig, “Underwater acoustic communication channel: Prop-agation models and statistical characterization,” IEEE Communication Magazine, vol. 47, pp. 84 – 89, Feb. 2009.
[8] P. A. van Walree, “Propagation and scattering effects in underwater acoustic com-munication channels,” IEEE Journal of Oceanic Engineering, vol. 38 ), pp. 614 – 631, Oct. 2013.
[9] W. K. A. Quazi, “Underwater acoustic communications,” IEEE Communications Magazine, vol. 20, pp. 24–30, Mar. 1982.
[10] R. E. Williams and H. F. Battestin, “Coherent recombination of acoustic multipath signals propagated in the deep ocean,” The Journal of the Acoustical Society of America, vol. 50, no. 6, pp. 1433 – 1442, 1971.
46
[11] M. Stojanovic, “Efficient processing of acoustic signals for high-rate information transmission over sparse underwater channels,” Physical communications, vol. 1, pp. 146–161, June 2008.
[12] A. K. Jagannatham and B. D. Rao, “Superimposed pilots vs. conventional pilots for channel estimation,” 2006 Fortieth Asilomar Conference on Signals, Systems and Computers, pp. 767–771, Oct. 2006.
[13] H. Kim, J. Chung, M. L. Nordenvaad, and J. Kim, “Superimposed pilots aided estimation of phase varying channels for underwater acoustic communication,” OCEANS 2014 - TAIPEI, pp. 1–4, Apr. 2014.
[14] D. V. Welden, M. Moeneclaey, and H. Steendam, “Parametric versus nonparametric data-aided channel estimation in a multipath fading environment.,” 2006 Symposium on Communications and Vehicular Technology, pp. 25 – 28, Nov. 2006.
[15] M. Stojanovic, J. Catipovic, and J. G. Proakis, “Adaptive multichannel combin-ing and equalization for underwater acoustic communications,” The Journal of the Acoustical Society of America 94, vol. 94, pp. 1621–1631, Aug. 1993.
[16] M. Stojanovic, L. Freitag, and M. Johnson, “Channel-estimation-based adaptive equalization of underwater acoustic signals,” Oceans ’99. MTS/IEEE. Riding the Crest into the 21st Century. Conference and Exhibition. Conference Proceedings (IEEE Cat. No.99CH37008), Sept. 1999.
[17] M. Lopez and A. Singer, “A dfe coefficient placement algorithm for sparse rever-berant channels,” IEEE Transactions on Communications, vol. 49, pp. 1334 – 1338, Aug. 2001.
[18] W. Li and J. C. Preisig, “Estimation of rapidly time-varying sparse channels,” IEEE Journal of Oceanic Engineering, vol. 32, pp. 927 – 939, Oct. 2007.
[19] J. Labat, G. Lapierre, and J. Trubuil, “Iterative equalization for underwater acoustic channels potentiality for the tpident system,” Oceans 2003. Celebrating the Past ... Teaming Toward the Future (IEEE Cat. No.03CH37492), Sept. 2003.
[20] E. Sozer, J. Proakis, and F. Blackmon, “Iterative equalization and decoding tech-niques for shallow water acoustic channels,” MTS/IEEE Oceans 2001. An Ocean Odyssey. Conference Proceedings (IEEE Cat. No.01CH37295), Nov. 2001.
[21] T. Oberg, B. Nilsson, N. Olofsson, M. L. Nordenvaad, and E. Sangfelt, “Underwater communication link with iterative equalization,” OCEANS 2006, Sept. 2006. [22] T. Eggen, A. Baggeroer, and J. Preisig, “Communication over doppler spread
