ESTIMATION OF A MULTIPATH SIGNAL IN A WIRELESS COMMUNICATION SYSTEM

Abstract

Described herein is a method and system for estimating a time-varying multipath communication channel. The system comprises a receiver for receiving a signal comprising a data packet transmitted over the multipath communication channel. The system also comprises a frequency-domain channel estimator that generates a frequency-domain channel estimate based on the received signal. The system further comprises a time-domain channel modeller that receives the frequency-domain channel estimate and generates an estimate of parameters in a time-domain model of the multipath communication channel. The model characterises the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a gain for the tap. The tap parameters are constant over the data packet.

Description

FIELD OF THE INVENTION

The present invention relates to the field of wireless communications. In particular the present invention relates to joint delay and Doppler frequency estimation from a received signal after the signal has been transmitted through a multipath channel.

BACKGROUND OF THE INVENTION

Wireless communication systems may be represented in terms of a transmitter 100 and receiver 104, separated by a channel 102, as shown in FIG. 1. The transmitter transforms the data into a signal suitable for transmission over the channel. For the purposes of determining the transmitted data, the goal of the receiver 104 is to remove the effects of the channel distortions from the signal and to obtain an estimate of the data.

The channel 102 represents the effects induced by the environment surrounding the wireless communications system. The channel 102 may distort the transmitted signal in some way. Channel distortions may include amplitude distortions, frequency offsets, phase offsets, Doppler effects, distortions resulting from multipath channels, additive noise or interference.

The receiver 104 may include a channel estimator. The channel estimator may observe a received signal that has been distorted by transmission over the channel 102, and generate a channel estimate based upon this observation.

Reference to any prior art in the specification is not, and should not be taken as, an acknowledgement or any form of suggestion that this prior art forms part of the common general knowledge in Australia or any other jurisdiction or that this prior art could reasonably be expected to be ascertained, understood and regarded as relevant by a person skilled in the art.

SUMMARY OF THE INVENTION

According to a first aspect of the present invention there is provided a method of estimating a time-varying multipath communication channel, the method comprising:

• • receiving a signal comprising a data packet transmitted over the multipath communication channel; and, from the received signal
  • generating an estimate of parameters in a model of the multipath communication channel, the model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters for estimation defining a delay, a Doppler and a complex gain for the tap,
  • wherein the estimated tap parameters are constant for the data packet.

The data packet may be an orthogonal frequency division multiplex (OFDM) packet.

According to another aspect of the present invention there is provided a method of estimating a time-varying multipath communication channel, the method comprising:

• • receiving a frequency-domain estimate of the multipath communication channel, the frequency-domain estimate deriving from a data packet transmitted over the multipath communication channel; and, based on the frequency-domain estimate
  • generating an estimate of parameters in a time-domain model of the multipath communication channel, the time-domain model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a complex gain for the tap,
  • wherein the estimated tap parameters are constant for the data packet.

The method of generating estimates of the model parameters may comprise:

• • defining a search window for the delay and Doppler tap parameters;
  • dividing the search window into a plurality of bins;
  • determining a correlator metric for the respective bins; and
  • identifying one or more peaks in the determined correlator metrics; and
  • outputting estimates of the delay and Doppler tap parameters based on the identified peaks.

The method of generating estimates of the model parameters may comprise;

• • receiving a channel estimate;
  • estimating a delay and Doppler tap parameter for a detected multipath component;
  • performing a least-squares estimation to generate an estimate of a complex gain tap parameter for the detected tap; and
  • subtracting the detected multipath component from the channel estimate, using the estimated delay, Doppler and complex gain tap parameters to quantify the detected multipath component.

The method may be repeated to iteratively remove detected multipath components from the channel estimate.

The method of generating estimates of the model parameters may comprise a method of refining estimates, comprising:

• • receiving an initial channel estimate;
  • receiving a set of estimated tap parameters for a plurality of multipath components;
  • subtracting all of the plurality of multipath components from the initial channel estimate except for a currently-considered multipath component;
  • generating a refined estimate of the delay tap parameter and the Doppler tap parameter for the currently-considered multipath component;
  • performing a least-squares estimation to generate a refined estimate of the complex gain tap parameter for the currently-considered multipath component.

The method of refining estimates may be repeated to generate refined estimates for a plurality of multipath components.

Also disclosed is an apparatus to perform the methods described in the preceding paragraphs.

The invention also resides in instructions executable by an apparatus to implement the described methods and in such instructions when stored in a storage medium readable by the apparatus.

According to a further aspect of the present invention there is provided a system for estimating a time-varying multipath communication channel, the system comprising:

• • a receiver for receiving a signal comprising a data packet transmitted over the multipath communication channel;
  • a frequency-domain channel estimator that generates a frequency-domain channel estimate based on the received signal; and
  • a time-domain channel modeller that receives the frequency-domain channel estimate and generates an estimate of parameters in a time-domain model of the multipath communication channel, the model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a complex gain for the tap,
  • wherein the tap parameters are constant over the data packet.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the present invention will now be described with reference to the drawings, in which:

FIG. 1 is a schematic drawing of a communications system;

FIG. 2 illustrates an example of an environment with a two-path channel;

FIG. 3 is an example of a time domain channel having two paths corresponding to the environment of FIG. 2;

FIG. 4 is a schematic representation of an embodiment of a system for processing a received signal and estimating a channel model;

FIG. 5A is a magnitude contour plot of an ambiguity function of an OFDM system with PSK modulation and a signal having 90 symbols;

FIG. 5B is a magnitude contour plot of an ambiguity function of an OFDM system with PSK modulation and a signal having 242 symbols;

FIG. 6A is a flow diagram of a first method for joint delay/Doppler estimation;

FIG. 6B is a flow diagram of an iterative method for joint delay/Doppler estimation via successive multipath cancellation;

FIG. 6C is a flow diagram of an iterative method for refining a set of multipath estimates;

FIG. 7 shows a performance evaluation of the methods of FIG. 6 using Monte Carlo simulations with a single tap; and

FIG. 8 shows a performance evaluation of the methods of FIG. 6 using Monte Carlo simulations with a three tap channel.

DETAILED DESCRIPTION OF THE EMBODIMENTS

Embodiments of an estimator are described for the estimation of multipath components when a signal is transmitted in a mobile environment. The estimator may be used for doubly dispersive channels, where the frequency selectivity of the channel changes with time.

The described techniques have potential application to wireless communications systems, e.g. DVB-T, DVB-H, IEEE 802.11, IEEE 802.16, 3GPP2, Dedicated Short Range Communications (DSRC), Communications Access for Land Mobiles (CALM), and proprietary systems.

An electromagnetic signal transmitted in a mobile environment is scattered by objects in the environment. The receiver 104, which may also possess a certain velocity relative to the transmitter 100, receives the superimposition of delayed and frequency-shifted versions of the transmitted signal caused by reflections from the objects. Each of these versions of the transmitted signal is referred to as a multipath component, a multipath tap, a channel tap, and/or a tap and is represented by P in the equations herein. This type of channel 102 is commonly referred to as a multipath channel. The delay, also referred to as time of arrival (ToA), of each multipath component is caused by the time it takes for the signal to travel the path distance and is represented by τ. The frequency offset, also referred to as Doppler frequency and/or Doppler, is caused by the Doppler effect and is represented by ν.

The estimation of the delay and/or Doppler of each multipath component may be output to a user interface such as a graphical display and/or stored in memory for use in radar applications for target detection and localisation or for channel characterisation in signal communication. The estimation may be provided to a radar application for target detection and localisation, for example in a digital form for use by radar application software. The estimation may be provided to a signal communication application for channel characterisation. The estimation may be used in a signal communication application for the development of efficient and reliable communication schemes.

If the coherence bandwidth of the multipath is small compared with the bandwidth of the transmitted signal, the resulting channel is frequency selective. Orthogonal frequency division multiplexing (OFDM) is a transmission strategy specifically designed to combat frequency-selective channels with relatively low receiver complexity. In OFDM, the signal bandwidth is divided into several non-overlapping (hence orthogonal) narrowband subcarriers, such that each subchannel is frequency non-selective.

FIG. 2 shows an example of an environment with an inflector 200 inducing a two path channel between the transmitter 100 and receiver 104, where:

• • T is a point representing the position of the transmitter 100;
  • R is a point representing the position of the receiver 104; and
  • Q is a point representing the position of the signal inflector 200.

FIG. 3 shows an example of a channel 102 in the time domain (with normalised power delay profile) corresponding to the environment of FIG. 2. The direct path TR corresponds to channel tap h₁ 700 at time t₁. The inflected path via inflector 200 corresponds to channel tap h₂ 702 at delay t₂. In this example h₂ 702 has lower power relative to tap h₁ 700 due to increased propagation loss (as the inflected path is longer than the direct path) and attenuation at the point of inflection 200. The time difference between the two channel taps is Δt₁₂=t₂−t₁. The instantaneous phase, and rate of change of phase, of taps h₁ 700 and h₂ 702 may also differ.

As seen in FIG. 4, a transmitted signal 106 is subjected to the channel 102 induced by the environment. The received signal is collected at the receive antenna 108, and input to the receiver 104. The receiver 104 processes the received signal to determine the transmitted data. FIG. 4 shows receiver processing in the context of the channel estimation and may omit details pertaining to the common operation of a receiver 104 apparent to those skilled in the art. For example, the receiver 104 may also make the transmitted data available to other components of a communication system.

The receiver 104 outputs receiver information, for example received signal samples, to a frequency-domain channel estimator 110 that provides an initial channel. estimate 112 in the frequency domain as the channel varies over the duration of a packet. Our previous commonly-assigned International (PCT) Applications, PCT/2006/AU001201, PCT/2007/AU000231 and PCT/2007/AU001506 published under WIPO publication numbers WO2007022564 WO2007095697, WO2008040088 (respectively), the contents of which are incorporated herein by reference, disclose systems and methods for the channel estimator 110.

The frequency domain channel estimator 110 is shown as a separate functional module to the receiver 104, but in practice the channel estimator 110 may be a component of the receiver 104.

The frequency-domain channel estimate 112 is provided to a time-domain channel modeller 114 that estimates time domain parameters including the delay, complex gain and Doppler for each multipath tap. The methods used in the time-domain channel modeller are described in more detail below.

The time-domain channel modeller 114 may be a functional module incorporated in the receiver 104. Alternatively, the modeller 114 may be implemented as a separate unit in data communication with the receiver 104. The modeller 114 may also be run off-line, that is, using data that was received earlier by receiver 104.

The functional modules described herein, including the frequency-domain channel estimator 110 and the time-domain channel modeller 114, may be implemented in hardware, for example application-specific integrated circuits (ASICs). Other hardware implementations include, but are not limited to, field-programmable gate arrays (FPGAs), structured ASICs, digital signal processors and discrete logic. Alternatively, the functional modules may be implemented as software, such as one or more application programs executable within a computer system. The software may be stored in a computer-readable medium and be loaded into a computer system from the computer readable medium for execution by the computer system. A computer readable medium having a computer program recorded on it is a computer program product. Examples of such media include, but are not limited to CD-ROMs, hard disk drives, a ROM or integrated circuit. Program code may also be transmitted via computer-readable transmission media, for example a radio transmission channel or a networked connection to another computer or networked device.

The methods used in the modeller 114 use the system model described in the following section.

1. System Model

Consider a K subcarrier OFDM system, where packets of length L OFDM symbols are transmitted. Let X_l,k ∈χ denote the symbol of subcarrier k of OFDM symbol l, where χ⊂ is an arbitrary complex constellation. Each symbol is drawn at random, independently and uniformly from χ and has unit average energy, i.e. [|X_l,k|²]=1 and [X_l,kX*_n,m]=0 for any n≠l or m≠k. The receiver is assumed to have complete knowledge of the transmitted symbols X_l,k, e.g. a pilot signal, or from the feedback of decoder decisions.

We let χ(t)=Σ_t=0^L−1 χ₁(t) denote the continuous-time transmitted OFDM signal, where

$\begin{array}{cc}{x}_l\left(t\right)=\frac{1}{\sqrt{\mathrm{KL}}}\sum _{k=-\lfloor K/2\rfloor }^{\lceil K/2\rceil -1}X_{l,k}{}^{\mathrm{j2\pi }k\left(t-T_{\mathrm{cp}}\right)/T}w\left(t-{\mathrm{lT}}_d\right),& \left(1\right)\end{array}$

is the lth OFDM symbol, T_d is the OFDM symbol duration (seconds). 1/T is the subcarrier spacing (Hz). T_cp=T_d−T_cp is the cyclic prefix duration (seconds). and w(t) is a windowing function such that

$\begin{array}{cc}w\left(t\right)=\{\begin{array}{cc}\tilde{w}\left(t\right)& 0\le t\le T_d\\ 0& \mathrm{otherwise},\end{array}& \left(2\right)\end{array}$

and ∫₀^Td|w(t)|²dt=1. Hence the transmitted signal has unit energy, i.e. [∫_−∞^∞|χ(t)|² dt]=1. The most simplest windowing function is when {tilde over (w)}(t)=1/√{square root over (T_d)}. Note that the requirement that w(t)=0 for t∉(0, T_d) is not necessarily required, but we assume this for simplicity. In practice, (1) is implemented via the inverse discrete Fourier transform (IDFT).

We assume transmit and receive filter responses g_Tx(t) and g_Rx(t) respectively, and a doubly dispersive channel model with P multipath components,

$\begin{array}{cc}c\left(t,\tau \right)=\sum _{p=1}^Pa_p{}^{-j2\pi v_pt}\delta \left(\tau -{\tau }_p\right).& \left(3\right)\end{array}$

where α_p, ν_p and τ_p denote the pth complex gain, Doppler and delay respectively. Note that these parameters are assumed to he constant over the L symbols of the OFDM packet. Hence, we write the overall channel response as

$\begin{array}{cc}h\left(t,\tau \right)=\left(g*c\right)\left(t,\tau \right)=\sum _{p=1}^{P-1}a_p{}^{-j2\pi v_pt}g\left(\tau -{\tau }_p\right),& \left(4\right)\end{array}$

where g(t)=(g_Tx*g_Rx) is the combined transmit/receive filter response and * denotes the convolution.

An ambiguity function for the transmitted signal x(t) is defined as:

A _χ(τ, ν)=∫_−∞^∝χ(t)χ*(t−τ)c^−j2πνt dt.  (5)

For the OFDM signal defined in (1), the ambiguity function can be written as follows.

$\begin{array}{cc}{A}_x\left(\tau ,v\right)=\frac{1}{\mathrm{KL}}\sum _{l,k,l^{\prime },k^{\prime }}X_{l,k}X_{l^{\prime },k^{\prime }}^{*}{}^{-j2\pi \left(k-k^{\prime }\right)T_{\mathrm{cp}}/T}{}^{j2\pi k^{\prime }\tau /T}{}^{-j2\pi vT_Al}A_w\left(\tau +\left(l^{\prime }-l\right)T_d,v+\left(k^{\prime }-k\right)/T\right).\mathrm{where}& \left(6\right)\\ A_w\left(\tau ,v\right)={\int }_{-\infty }^{\infty }w\left(t\right)w^{*}\left(t-\tau \right){}^{-j2\pi \mathrm{vt}}t& \left(7\right)\end{array}$

is the ambiguity function of w(t). In practice, w(t) is usually designed such that A_w(τ,ν)≈0 for |τ|>T_cp and |ν|>1/T, which implies

$\begin{array}{cc}{A}_x\left(\tau ,v\right)=\frac{A_w\left(\tau ,v\right)}{\mathrm{KL}}\sum _{k,l}{X_{t,k}}^2{}^{j2\pi \mathrm{kr}/T}{}^{-j2\pi {\mathrm{vT}}_dl}.& \left(8\right)\end{array}$

For example, for the special case of ŵ(t)=1/√{square root over (T_d)}, then we have

$\begin{array}{cc}{A}_w\left(\tau ,v\right)=\{\begin{array}{cc}\left(1-\frac{\tau }{T_d}\right){}^{-j\pi v\left(T_d+\tau \right)}\sin c\left(\pi v\left(T_d-\tau \right)\right)& 0\le \tau <T_d\\ \left(1+\frac{\tau }{T_d}\right){}^{-j\pi v\left(T_d+\tau \right)}\sin c\left(\pi v\left(T_d+\tau \right)\right)& -T_d\le \tau \le 0\\ 0& \mathrm{otherwise}.\end{array}& \left(9\right)\end{array}$

In practice OFDM systems are designed such that the cyclic prefix is greater than the delay spread of the channel. i.e. we are primarily concerned with 0<τ<<T_d. In addition, Doppler frequencies induced by the mobile channel are typically |ν|<<1/T. Under these conditions. A_w(τ,ν)≈1 and the ambiguity function is approximately

$\begin{array}{cc}{A}_x\left(\tau ,v\right)\approx \frac{1}{\mathrm{KL}}\sum _{k,l}{X_{l,k}}^2{}^{j2\pi k\tau /T}{}^{-j2\pi {\mathrm{vT}}_dl}.& \left(10\right)\end{array}$

For the special case when |X_l,k|=1, i.e. PSK modulation, then the ambiguity function further simplifies to

$\begin{array}{cc}{A}_x\left(\tau ,v\right)\approx \frac{1}{\mathrm{KL}}{}^{-j2\pi \left[K/2\right]\tau /T}{}^{-j\pi \left(L-1\right)vT_d)}\frac{\sin \left(\pi \tau K/T\right)}{\sin \left(\pi \tau /T\right)}\frac{\sin \left(\pi {\mathrm{vT}}_dL\right)}{\sin \left(\pi {\mathrm{vT}}_d\right)}.& \left(11\right)\end{array}$

From this expression it is apparent that the sin/sin terms will introduce sidelobes, which are de-coupled in time and frequency. As an example, the ambiguity function of the IEEE 802.11 standard is plotted in FIG. 5A, showing a 64 subcarrier system with 52 active subcarriers and zero DC subcarrier, T_d=8 μsec and T=6.4 μsec. In this example, K=52 and the summation terms in (10) range from k=−26, . . . −1 to k=1, . . . , 26. For OFDM systems where there is no DC subcarrier, (11) becomes,

$\begin{array}{cc}{A}_x\left(\tau ,v\right)\approx \frac{2{}^{-j\pi \left(L-1\right)l/T_d}}{\mathrm{KL}}\cos \left[\pi \left(K/2+1\right)\left(\tau /T\right)\right]\frac{\sin \left(\mathrm{\pi \tau }K/\left(2T\right)\right)}{\sin \left(\mathrm{\pi \tau }/T\right)}\frac{\sin \left(\pi {\mathrm{vT}}_dL\right)}{\sin \left(\pi {\mathrm{vT}}_d\right)}.& \left(12\right)\end{array}$

The vertical axis represents delay and the horizontal axis represents Doppler. FIG. 5B shows the ambiguity function when the number of symbols, L, is increased to 242. It may be seen that increasing L improves the Doppler resolution but not the delay resolution, which is dependent only on K and the subcarrier spacing (1/T).

In practice the OFDM signal of (1) is filtered and upconverted in transmitter 100, transmitted through multipath channel 102, down-converted and receive filtered in receiver 104. Assuming the up/down conversion process is done perfectly, the received continuous-time baseband signal may be written as:

$\begin{array}{cc}\begin{array}{c}y\left(t\right)=\sum _{l^{\prime }}x_{l^{\prime }}\left(t\right)*h\left(t,\tau \right)+\delta \left(t\right)*g\left(t\right)\\ =\sum _{l^{\prime }}{\int }_{-\infty }^{\infty }x_{l^{\prime }}\left(t-\tau \right)h\left(t,\tau \right)\tau +z\left(t\right)\\ =\frac{1}{\sqrt{{\mathrm{KLT}}_d}}\sum _{l^{\prime },k^{\prime }}X_{l^{\prime },k^{\prime }}{\int }_{-\infty }^{\infty }{}^{j\pi k^{\prime }\left(t-\tau -T_{\mathrm{cp}}\right)/T}\\ w\left(t-\tau -l^{\prime }T_d\right)h\left(t,\tau \right)\tau +z\left(t\right)\\ =\frac{1}{\sqrt{{\mathrm{KLT}}_d}}\sum _{l^{\prime },k^{\prime },p}a_pX_{l^{\prime },k^{\prime }}{\int }_{-\infty }^{\infty }{}^{j\pi k^{\prime }\left(t-\tau -T_{\mathrm{cp}}\right)/T}\\ w\left(t-\tau -l^{\prime }T_d\right){}^{-j2\pi v_pt}g\left(\tau -{\tau }_p\right)\tau +z\left(t\right)\\ =\frac{1}{\sqrt{{\mathrm{KLT}}_d}}\sum _{l^{\prime },k^{\prime },p}a_pX_{l^{\prime },k^{\prime }}{\Phi }_{l^{\prime },k^{\prime }}\left(t,{\tau }_p\right)\\ {}^{j2\pi k^{\prime }\left(t-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pt}+z\left(t\right)\end{array}& \left(13\right)\end{array}$

where {tilde over (z)}(t) is an AWGN process and

φ_l,k(t,τ_p)=∫_−∞^∞g(τ−τ_p)w(t−τ−lT_d)e^−j2πkτ/T dτ  (14)

To simplify (14) let

G(f)={g(t)}∫_−∞^∞g(t)e^−j2πf1dt  (15)

W(f)={w(t)}∫_−∞^∞w(t)e^−j2πf2dt  (16)

The Fourier transform of convolution in (14) is the same as the multiplication of the constituent Fourier transform. i.e.

$\begin{array}{cc}\begin{array}{c}ℱ\left\{{\Phi }_{l,k}\left(t,{\tau }_p\right)\right\}=ℱ\left\{g\left(t-{\tau }_p\right){}^{-j2\pi \mathrm{kt}/T}\right\}ℱ\left\{w\left(t-{\mathrm{lT}}_d\right)\right\}\\ ={}^{-j2\pi {\tau }_p\left(f+k/T\right)}{}^{-j2\pi {\mathrm{lT}}_df}G\left(f+k/T\right)W\left(f\right).\end{array}& \left(17\right)\end{array}$

In practical systems, G(f) has a bandwidth that exceeds the bandwidth of the signal, i.e. K/T. In addition W(f) has a bandwidth of approximately 1/T_d. Assuming ideal filtering, then G(f+k/T)=1, −F_s/2−k/T<f<F_s/2−k/T and zero otherwise. e.g. g(t)=sinc(πtF_s), where F_s>K/T is the bandwidth of the filter. Hence we may write

{φ_l,k(t,τ_p)}≈e^{−j2πTp(f+k/T)e−j2π/Tdf}W(f).  (18)

and therefore

φ_l,k(t,τ_p)≈e^{−j2πkτp/T}w(t−τ_p−lT_d).  (19)

The received signal can now be written as

$\begin{array}{cc}y\left(t\right)=\frac{1}{\sqrt{\mathrm{KL}}}\sum _{l^{\prime },k^{\prime },p}a_pX_{p,k^{\prime }}{}^{j2\pi k^{\prime }\left(t-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pl}{}^{-j2\pi k^{\prime }{\tau }_p/T}w\left(t-{\tau }_p-l^{\prime }T_d\right)+z\left(t\right)& \left(20\right)\end{array}$

The receiver performs the matched filter to the transmitted sinusoids (less the cyclic prefix),

$\begin{array}{cc}{Y}_{l,k}=\frac{1}{\sqrt{\mathrm{KL}}}{\int }_{T_{\mathrm{cp}}+{\mathrm{lT}}_d}^{\left(l+1\right)T_d}y\left(t\right)w^{*}\left(t-{\mathrm{lT}}_d\right){}^{-j2\pi k\left(t-T_{\mathrm{cp}}\right)/T}t=\frac{1}{\mathrm{KL}}\sum _{l^{\prime },k^{\prime },p}a_pX_{l^{\prime },k^{\prime }}{\int }_{T_{\mathrm{cp}}+\mathrm{lTd}}^{\left(l+1\right)T_d}{}^{j2\pi k^{\prime }(t-T_{\mathrm{cp}}/T}{}^{-j2\pi v_pt}{}^{-j2\pi k^{\prime }{\tau }_p/T}\times w\left(t-{\tau }_p-l^{\prime }T_d\right)w^{*}\left(t-{\mathrm{lT}}_d\right){}^{-j2\pi k\left(t-T_{\mathrm{cp}}\right)/T}t+Z_{l,k}=\frac{1}{\mathrm{KL}}\sum _{k^{\prime },p}^{}a_pX_{l,k^{\prime }}{\int }_{T_{\mathrm{cp}}+{\mathrm{lT}}_d}^{\left(l+1\right)T_d}{}^{j2\pi k^{\prime }\left(t-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pt}{}^{-j2\pi k^{\prime }{\tau }_p/T}\times w\left(t-{\tau }_p-{\mathrm{lT}}_d\right)w^{*}\left(t-{\mathrm{lT}}_d\right){}^{-j2\pi k\left(t-T_{\mathrm{cp}}\right)/T}t+Z_{l,k}=\frac{1}{\mathrm{KL}}\sum _{k^{\prime },p}a_pX_{l,k^{\prime }}{}^{-j2\pi k^{\prime }{\tau }_p/T}\times {\int }_{T_{\mathrm{cp}}+{\mathrm{lT}}_d}^{\left(l+1\right)T_d}w\left(t-{\tau }_p-{\mathrm{lT}}_d\right)w^{*}\left(t-{\mathrm{lT}}_d\right){}^{j2\pi \left(k^{\prime }-k\right)\left(t-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pt}t+Z_{l,k}=\frac{1}{\mathrm{KL}}\sum _{k^{\prime },p}a_pX_{l,k}{}^{-j2\pi k^{\prime }{\tau }_p/T}\times {\int }_{T_{\mathrm{cp}}}^{T_d}w\left(t-{\tau }_p\right)w^{*}\left(t\right){}^{j2\pi \left(k^{\prime }-k\right)\left(t+{\mathrm{lT}}_d-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_p\left(t+{\mathrm{lT}}_d\right)}t+Z_{l,k}=\frac{1}{\mathrm{KL}}\sum _{k^{\prime },p}a_pX_{l,k^{\prime }}{}^{-j2\pi k^{\prime }{\tau }_p/T}{}^{j2\pi \left(k^{\prime }-k\right)\left({\mathrm{lT}}_d-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pT_dl}\times {\int }_{T_{\mathrm{cp}}}^{T_d}w\left(t-{\tau }_p\right)w^{*}\left(t\right){}^{j2\pi \left(k^{\prime }-k\right)t/T}{}^{-j2\pi v_pt}t+Z_{l,k}=\frac{1}{\mathrm{KL}}\sum _{k^{\prime },p}X_{l,k^{\prime }}a_p{}^{-j2\pi k^{\prime }{\tau }_p/T}{}^{j2\pi \left(k^{\prime }-k\right)\left({\mathrm{lT}}_d-T_{\mathrm{cp}}\right)/T}{}^{-j2\pi v_pT_dl}{\hat{A}}_w\left({\tau }_p,v_p+\left(k-k^{\prime }\right)/T\right)+Z_{l,k},\mathrm{where}& \left(21\right)\\ {\hat{A}}_w\left(\tau ,v\right)={\int }_{T_{\mathrm{cp}}}^{T_d}w\left(t-\tau \right)w^{*}\left(t\right){}^{-j2\pi \mathrm{vt}}t,& \left(22\right)\end{array}$

which resembles the ambiguity function of w(t). The second line of (21) follows since the integral is approximately zero for any l′ not equal to l. Typically, w(t) has a bandwidth of 1/T_d, and hence Â_w(τ_p, ν_p+(k−k′)/T)≈0 for k′≠k. Therefore,

$\begin{array}{cc}{Y}_{l,k}=\frac{1}{\mathrm{KL}}\sum _pa_pX_{l,k}{}^{-j2\pi k{\tau }_p/T}{}^{-j2\pi v_pT_dl}{\hat{A}}_w\left({\tau }_p,v_p\right)+Z_{l,k}.& \left(23\right)\end{array}$

For the special case of ŵ(t)=1/√{square root over (T_d)}, i.e. a rectangular window function, then

$\begin{array}{cc}{\hat{A}}_w\left(\tau ,v\right)=\frac{T}{T_d}{}^{-j\pi v\left(T_d+T_{\mathrm{cp}}\right)}\sin c\left(\pi \mathrm{vT}\right)0<\tau <T_{\mathrm{cp}}& \left(24\right)\end{array}$

In practice, 0<max_p τ_p<T_cp and max_p ν_p<<1/T. Therefore

${\hat{A}}_w\left({\tau }_p,v_p\right)\approx \frac{T}{T_d}{}^{-j\pi v\left(T_d+T_{\mathrm{cp}}\right)}$

and hence (23) simplifies to

$\begin{array}{cc}{Y}_{l,k}=X_{l,k}\sum _{p=1}^P{\tilde{a}}_p{}^{-j2\pi v_pT_dl}{}^{-j2\pi k{\tau }_p/T}+Z_{l,k},\mathrm{where}{\tilde{a}}_p=\frac{1}{\mathrm{KL}}\frac{T}{T_d}{}^{-j\pi v_p\left(T_d+T_{\mathrm{cp}}\right)}a_p.& \left(25\right)\end{array}$

2. Joint Delay/Doppler Estimation

It is assumed that the receiver 104 has perfect knowledge of the transmitted symbols X_l,k. This may, for example, result from the use of pilot symbols or the feedback of correct decoder decisions in the receiver. As a result, we may assume the modeller 114 has perfect knowledge of the channel gains H_l,k:

Y _l,k =H _l,k X _l,k +Z _l,k.  (26)

where from (25),

$\begin{array}{cc}{H}_{l,k}=\sum _{p=1}^P{\tilde{a}}_p{}^{-j2\pi v_pT_dl}{}^{-j2\pi k{\tau }_p/T}& \left(27\right)\end{array}$

Hence the problem is to estimate P, α_p, τ_p and ν_p for p=1, . . . P from H_l,k, or its noisy estimate.

Equation (27) is simply the superimposition of 2-dimensional exponential signals.

Two estimation algorithms are described below: the first is based on a conventional correlator approach; the second is based on a matching pursuit or basis pursuit algorithm which takes advantage of the sparse nature of the channel taps.

2.1. Method 600—Correlator Approach

A first method 600 for the time domain channel modeller is shown in FIG. 6A. Method 600 is based on a correlator approach and may for example be implemented as software instructions running on the modeller 114.

The input is the frequency domain channel estimate 112. In step 602 an initial delay search window and initial Doppler search window are defined. Then, in step 604 a grid-based approach is used to divide the (τ,ν) plane into bins. An iterative sequence follows. In step 606 G and S matrices (associated with the delay and Doppler effects respectively) are constructed as described in more detail below with reference to equations (30) to (32). In step 608 a correlator output metric is determined using equation (32). In step 610 the channel modeller identifies the P largest peaks in the correlator metric.

Let Ĥ be a noisy estimate obtained from Y given knowledge of X, i.e.

Ĥ=H+{circumflex over (Z)}.  (28)

where {circumflex over (Z)} is additive noise, not necessarily Gaussian, e.g. if Ĥ=Y⊙X is the zero-forcing estimate, then {circumflex over (Z)}=Z⊙X* is not necessarily Gaussian. The correlator approach computes the following metric,

ψ(τ, ν)=s^t(ν)Ĥg^t(τ),  (29)

where s(ν) is a L×1 column vector, whose lth element is given by s_l(ν)=e^−j2πνTdl, l=0, . . . , L−1, and g(τ) is a l×K column vector whose kth element is given by g_k(τ)=e^−j2πkτ/T.

To implement the correlator estimator we evaluate (29) over a region of the (τ, ν)-plane and take the P largest peaks. In particular, we use a grid based approach, i.e. we divide the region into discrete bins. Let M and N denote the number of delay and Doppler bins respectively. Let (τ_min, τ_max) and (ν_min, ν_max) denote the delay and Doppler region. Let {circumflex over (τ)}_m and {circumflex over (ν)}_n denote bin m, n of the delay-Doppler region, i.e.

$\begin{array}{cc}{\hat{\tau }}_m={\tau }_{\min }+m\frac{{\tau }_{\max }-{\tau }_{\min }}{M}& \left(30\right)\\ {\hat{v}}_n=v_{\min }+n\frac{v_{\max }-v_{\min }}{N}.& \left(31\right)\end{array}$

for m=0, . . . , M−1 and n=0, . . . , N−1. Let Ŝ and Ĝ denote the matrices constructed using the above bins, i.e. Ŝ_l,n=^−j2πνnTd and Ĝ_m,k=e^−j2πktm/T. The correlator output can be computed via the matrix multiplication,

ψ=Ŝ^tĤĜ^t.  (32)

The estimator then picks the {circumflex over (P)} largest peaks. Note that there is some similarity between (32) and the 2D DFT. Complexity reduction may be possible in computing (32) by taking advantage of the fast Fourier transform (FFT).

To investigate the correlator estimator from (27) we have,

$\begin{array}{cc}\begin{array}{c}\varphi \left(\tau ,v\right)=\sum _{p,k,l}a_p{}^{-j2\pi (v_p-v)}{}^{-j2\pi \left({\tau }_p-\tau \right)}+s^{†}\left(v\right)\hat{Z}g^{†}\left(\tau \right)\\ =\sum _pa_pA_x\left(\tau -{\tau }_p,v_p-v\right)+s^{†}\left(v\right)\hat{Z}g^{†}\left(\tau \right),\end{array}& \left(33\right)\end{array}$

i.e. the output correlator metric is the sum of complex scaled versions of the transmit signal ambiguity function. Hence we see that correlator estimator performance of method 600 is affected by sidelobes of the ambiguity function, which may cause false and miss detections. Although method 600 has relatively low complexity, it may suffer from missed and false channel detection due to the sidelobes in the ambiguity function. In particular, for strong line of sight (LOS) channels, the peaks corresponding to weaker channel taps may be hidden by the sidelobes of the LOS component.

2.2. Method 630—Estimation via successive multipath cancellation

From (27) we see that H_l,k is the sum of a sparse number of points. This suggests that the use of sparse estimation techniques may be beneficial in terms of reducing complexity with relatively low loss in performance. One approach is method 630 shown in FIG. 6B. This uses an orthogonal matching pursuit (OMP) algorithm to provide an iterative algorithm that successively removes previously-estimated multipath components from the original channel estimates as they are detected. A basis pursuit algorithm can also be used in the iterative process of removing previously-estimated multipath components. Method 630 takes advantage of the simplified signal derivations given above in equations (1) to (25) and allows the detection of components that would otherwise be hidden by erroneous sidelobes.

Let {circumflex over (τ)}_p, {circumflex over (ν)}_p be 1×p vectors of delay and Doppler estimates respectively.

Let S({circumflex over (ν)}_p)=(s({circumflex over (ν)}₁) . . . s({circumflex over (ν)}_p)) and G({circumflex over (Σ)}_p)=(g^t({circumflex over (τ)}₁) . . . g^t({circumflex over (τ)}_p))^t denote the corresponding matrices

The estimation of the complex gain vector α_p involves finding the solution to

$\begin{array}{cc}{\hat{a}}_p=\arg \underset{a_p}{\min }{\hat{H}-S\left({\hat{v}}_p\right)\mathrm{diag}\left(a_p\right)G\left({\hat{\tau }}_p\right)}^2,& \left(34\right)\end{array}$

which can be shown to be

{circumflex over (α)}_p=R_p⁻¹w_p.  (35)

where w_p is a p×1 column vector that is computed recursively using,

$\begin{array}{cc}{w}_p=\left(\begin{array}{c}{w}_{p-1}\\ \sum _{l=0}^{L-1}\sum _{k=0}^{K-1}{\hat{H}}_{l,k}g_k^{*}\left({\hat{\tau }}_p\right)s_l^{*}\left({\hat{v}}_p\right)\end{array}\right)=\left(\begin{array}{c}{w}_{p-1}\\ s^{†}\left({\hat{v}}_p\right)\hat{H}g^{†}\left({\hat{\tau }}_p\right)\end{array}\right).& \left(36\right)\end{array}$

The p×p matrix R_p is also computed recursively, i.e.

$\begin{array}{cc}{R}_p=\left(\begin{array}{cc}{R}_{p-1}& r^{†}\\ r& \mathrm{KL}\end{array}\right).& \left(37\right)\end{array}$

where r is a 1×p−1 row vector whose nth element, n=1, . . . , p−1, is given by

$\begin{array}{cc}\begin{array}{c}{\tau }_n=\sum _{l=0}^{L-1}\sum _{k=0}^{K-1}g_k^{*}\left({\hat{\tau }}_p\right)g_k\left({\hat{\tau }}_n\right)s_l^{*}\left({\hat{v}}_p\right)s_l\left({\hat{v}}_n\right)\\ =g\left({\hat{\tau }}_n\right)g^{†}\left({\hat{\tau }}_p\right){\left(s^{†}\left({\hat{v}}_n\right)s\left({\hat{v}}_p\right)\right)}^{*}\\ ={}^{-j\pi \left(L-1\right)\left({\hat{v}}_n-{\hat{v}}_p\right)T_d}{}^{-\mathrm{j\pi }\left(K-1\right)\left({\hat{\tau }}_n-{\hat{\tau }}_p\right)/T}\left(39\right)\\ \frac{\sin \left(\pi L\left({\hat{v}}_n-{\hat{v}}_p\right)T_d\right)}{\sin \left(\pi \left({\hat{v}}_n-{\hat{v}}_p\right)T_d\right)}\frac{\sin \left(\pi K\left({\hat{\tau }}_n-{\hat{\tau }}_p\right)/T\right)}{\sin \left(\pi ({\hat{\tau }}_n-{\hat{\tau }}_p)/T\right)}.\end{array}& \left(38\right)\end{array}$

A detailed description of method 630 is given in algorithm 1.

Algorithm 1 Estimation via successive multipath cancellation.
1:  Δ0 = Ĥ
2:  p = 1
3:  while ∥Δp−1∥2 > a and p ≦ Pmax do
4:  ({circumflex over (τ)}p, {circumflex over (v)}p) = arg maxτ,y |s† (v)Δp−1g†(τ)|
5:  if p > 1 then
6:  Compute Rp and wp recursively using (37) and (36)
7:  else
8:  R1 = KL
9:  w1 = s†({circumflex over (v)}1)Ĥg†({circumflex over (τ)}1)
10: end if
11: âp = Rp−1wp.
12: Δp = Δp−1 − aps({circumflex over (v)}p)g({circumflex over (τ)}p)
13: p = p + 1
14: end while

As summarised in FIG. 6B, the input to method 630 is a channel estimate H. Method 630, commencing in step 632, performs an iteration over detected multipath components p, up to a maximum value of P_max.

For the current multipath component, a correlation metric is calculated to provide an estimate of the delay and Doppler in step 634 (see point 4. of algorithm 1). This may be determined using the 2-dimensional bisection approach described below in algorithm 3.

Then, in step 636 the time-domain channel modeller 114 constructs an R matrix and w vector for a least squares optimisation on the current multipath component p (see points 5. to 9. of algorithm 1). In step 638 the modeller 114 determines the estimate for the complex gain of the current multipath component p using the constructed R and w. In step 640 the detected multipath component is removed from the residual signal.

3. Iterative Refinement

The inventors have observed that the estimates produced by method 630 may be affected by components that have not been removed. This may lead to some degradation in the delay/Doppler resolution. Method 650 of FIG. 6C may be used to refine a set of estimated multipath components by removing all previously estimated multipath components except for a component of interest.

The input to method 650 is a set of p multipath estimates, for example a set of estimates 652 generated by the modeller 114 using method 630. In step 654 a first iteration is commenced that continues while the criterion set out in point 5 of algorithm 2 below is satisfied. For each of these outer iterations, a further set of iterations commences in step 656, stepping through each of the multipath component estimates in the set p. In step 658 all multipath components are removed except for the component of current interest. Then in step 660 the correlation approach is used .to refine the delay and Doppler estimates for the component of current interest. Then a least-squares optimisation refines the estimate of the complex gain for the component of current interest in step 662. Process flow returns to step 656 to step through the remaining estimates.

Once the loop of steps 656-662 is complete, the iteration criterion is updated in step 664 and process flow returns to step 654 to continue the iterative process if appropriate.

More detail on method 650 is shown in algorithm 2 as follows:

Note that in line 7 of algorithm 2, q denotes all elements excluding element q.

Algorithm 2 Estimate refinement algorithm.
1:  Given p >1 estimates, {circumflex over (τ)}, {circumflex over (v)} and â.
2:  Set Δ1 = Ĥ − S({circumflex over (v)})diag(â)G({circumflex over (τ)}).
3:  Set ∥Δ0∥2 = ∞
4:  i = 1
5:  while |∥Δi∥2 − ∥Δ1−1∥2| > ε and i ≦ Niter do
6:  for q = 1 to p do
7:  Δq = Ĥ − S({circumflex over (v)}q)diag(âq)G({circumflex over (τ)}q).
8:  ({circumflex over (τ)}q, {circumflex over (v)}q) = arg maxτ,v |s†(v)Δqg†(τ)|.
9:  Recursively create w and R.
10: â = R−1w
11: end for
12: i = i + 1
13: Δi = Ĥ − S({circumflex over (v)})diag(â)G({circumflex over (τ)})
14: end while

To perform the maximisation step in line 4 of algorithm 1 and line 8 of algorithm 2, the time domain channel modeller 114 may use the 2-D bisection approach as described in algorithm 3 below.

In its most naive implementation, algorithm 3 has a complexity of ML(N+K) complex multiply-accumulates per iteration. The complexity of line 14 may be reduced by exploiting the radix-2 algorithm.

Algorithm 3 2D Bisection Algorithm.
1:  Given the initial delay search window: τmin, τmax.
2:  Given the initial Doppler search window: νmin, νmax.
3:  Given number of delay and Doppler bins: M > 1 and N > 1 respectively.
4:  Given the search window overlap factor: β > 0.
5:  Given the error difference threshold: ε > 0.
6:  Given the maximum number of iterations: NIter.
7:  Set i = 0, μ−1 = ∞ and μo = 0.
8:  while |μi−1 − μi| > ε and i < NIter do
9:  Determine                                                                                                                                      delay                                                                                                                                      resolution                                        ,                                          Δτ                      =                                                                                                                                  τ                              max                                                        -                                                          τ                              min                                                                                M                                                .
10: Determine                                                                                                                                      Doppler                                                                                                                                      resolution                                        ,                                          Δν                      =                                                                                                                                  ν                              max                                                        -                                                          ν                              min                                                                                N                                                .
11: Construct delay bins, {circumflex over (τ)}m = τmin + mΔτ, for m = 0, . . . , M − 1.
12: Construct Doppler bins, {circumflex over (ν)}n = νmin + nΔν, for n = 0, . . . , N − 1.
13: Construct matrices G({circumflex over (τ)}) and S({circumflex over (ν)}).
14: Compute Ψ = GY†S.
15: Find largest magnitude squared value of Ψ and its indices, i.e.
    μi+1 = maxm,n|Ψm,n|2 and (m*, n*) = arg maxm,n|Ψm,n|2.
16: Update delay window, i.e. τmin = {circumflex over (τ)}m* − βΔτ and τmax = {circumflex over (τ)}m* + βΔτ.
17: Update Doppler window, i.e. νmin = {circumflex over (ν)}n* − βΔν and νmax = {circumflex over (ν)}n* + βΔν.
18: i = i + 1
19: end while

4. Performance Evaluation

To examine the performance of our estimation algorithm, we conducted a comprehensive set of Monte Carlo simulations. The assumption is made that P is known a-priori and that the multipath components are sufficiently separated to minimise false/miss detects due to the delay/Doppler resolution.

In all simulations, we assume an OFDM signal similar to the IEEE802.11 standard, with K=64 and 52 active subcarriers, a subcarrier spacing of 156:250 kHz and an OFDM symbol period of 8 μsec.

In these simulation results we assume the number of taps P is known a-priori. Only P iterations of algorithm 1 are performed and the detected taps are subsequently refined using algorithm 2. For each trial the delay/Doppler are drawn randomly, but not independently and identically distributed (i.i.d.). For each trial we chose a vector of P delays by continually drawing a vector P delays from an i.i.d. uniform distribution on the interval (τ_min·τ_max) until the minimum distance between the delays is greater than Δτ. The delays are then re-arranged in ascending order. Similarly for the Doppler offsets, however no sorting is performed. This allows us to generate the delay and Doppler offset of each tap randomly whilst maintaining a minimum delay/Doppler resolution. Note that Δτ≦(τ_max−τ_min)/P and similarly Δν≦(ν_max−ν_min)/P.

To ensure that the estimator performance is not dominated by the delay/Doppler bin resolution we first consider the P=1 case and compare the root mean squared (RMS) error with the Cramer-Rao lower bound (CRLB) of a 2-D sinusoid in AWGN, i.e. let r(m, n)=A exp [−j(w₁l+w₂k)]+z(m,n), for l=0, . . . L−1 and k=0, . . . K−1, where z(m,n)˜(0,σ²). Then

$\begin{array}{cc}\mathrm{var}\left[{\hat{w}}_1\right]\ge \frac{6}{\mathrm{KL}\left(L^2-1\right)}\frac{{\sigma }^2}{A^2}& \left(40\right)\\ \mathrm{var}\left[{\hat{w}}_2\right]\ge \frac{6}{\mathrm{KL}\left(K^2-1\right)}\frac{{\sigma }^2}{A^2}.& \left(41\right)\end{array}$

For P=1, algorithm 1 is the ML estimator of a 2-D sinusoid in AWGN. To observe this we set (τ_min, τ_max)=(0,200) nsec and (ν_min, ν_max)=(−500, 500) Hz. Results of the simulations are shown in FIG. 7, in which a single-tap estimation root-mean-square (RMS) error for the present estimation techniques (shown in solid circles) is compared to the CR lower bound shown as a dashed line. K=52. The upper curves are for L=128, the middle curves are for L=256 and the bottom curves are for L=512. For the single tap channel the RMS error is very low, e.g. at an SNR of 20 dB with L=256 the RMS delay/Doppler errors are approximately 0.4 nsec and 3.5 Hz respectively.

Now we consider a P=3 tap channel with power delay profile |α₁|²=−0.54, |α₂|²=−10.54 and |α₃|²=15.54 dB. The phase of each tap was chosen uniformly random on (0.2π).

FIG. 8 shows the results of a simulation for the 3-tap estimation. Again, the dashed lines show the CRLB. The solid lines with squares indicate the performance of algorithm 1 only.

Solid lines with circles show the performance of algorithm 1 followed by an application of the refinement algorithm 2 after 20 iterations. It can be seen that the refinement provides a noticeable improvement in accuracy.

For each trial, RMS error statistic's were only collected when all taps are detected, i.e. each tap is closest (in Euclidean distance) to a single estimate. It is possible for a particular tap to be closest to two or more estimated taps. These events are counted as miss-detections, but are not included in the RMS error statistics.

It will be understood that the invention disclosed and defined in this specification extends to all alternative combinations of two or more of the individual features mentioned or evident from the text or the drawings. All these different combinations constitute various alternative aspects of the invention.

The parameter estimation described herein has many applications. In one application the estimates of path delays and Dopplers provides a means of characterising a reflective environment of the communication channel. This approach may, for example, be used in Dedicated Short Range Communications (DSRC) systems. Methods for environment estimation in a wireless communication system are described, for example, in a commonly assigned PCT application PCT/AU2010/000768 claiming priority from Australian provisional application 2009902848 filed 19 Jun. 2009.

It will also be understood that the term “comprises” and its grammatical variants as used in this specification is equivalent to the term “includes” and should not be taken as excluding the presence of other elements or features.

Appendix: Derivation of Gain Vector Solution

The following derivation is for the equation used in line 10 of algorithm 2 and line 11 of algorithm 1 (Eq. 35).

We require the solution to

$\begin{array}{cc}{\hat{a}}_p=\arg \underset{a_p}{\min }{\hat{H}-S\left({\hat{v}}_p\right)\mathrm{diag}\left(a_p\right)G\left({\hat{\tau }}_p\right)}^2.& \left(42\right)\end{array}$

First, we rewrite the problem as a least-squares optimisation problem,

$\begin{array}{cc}{\hat{a}}_p=\arg \underset{a_p}{\min }{\mathrm{vec}\left(\hat{H}\right)-{\Omega }_pa_p}^2.& \left(43\right)\end{array}$

where vec(M) denotes the vectorize operation, i.e. stacking the columns of the matrix M, and Ω_p is the KL×p matrix,

$\begin{array}{cc}{\Omega }_p=\left(\begin{array}{ccc}{g}_0\left({\hat{\tau }}_1\right)s_0\left({\hat{v}}_1\right)& \dots & g_0\left({\hat{\tau }}_p\right)s_0\left({\hat{v}}_p\right)\\ g_1\left({\hat{\tau }}_1\right)s_0\left({\hat{v}}_1\right)& \dots & g_1\left({\hat{\tau }}_p\right)s_0\left({\hat{v}}_p\right)\\ ⋮& \dots & ⋮\\ g_{K-1}\left({\hat{\tau }}_1\right)s_0\left({\hat{v}}_1\right)& \dots & g_{K-1}\left({\hat{\tau }}_p\right)s_0\left({\hat{v}}_p\right)\\ ⋮& \dots & ⋮\\ g_0\left({\hat{\tau }}_1\right)s_1\left({\hat{v}}_1\right)& \dots & g_0\left({\hat{\tau }}_p\right)s_1\left({\hat{v}}_p\right)\\ g_1\left({\hat{\tau }}_1\right)s_1\left({\hat{v}}_1\right)& \dots & g_1\left({\hat{\tau }}_p\right)s_1\left({\hat{v}}_p\right)\\ ⋮& \dots & ⋮\\ g_{K-1}\left({\hat{\tau }}_1\right)s_1\left({\hat{v}}_1\right)& \dots & g_{K-1}\left({\hat{\tau }}_p\right)s_1\left({\hat{v}}_p\right)\\ ⋮& \dots & ⋮\\ g_{K-1}\left({\hat{\tau }}_1\right)s_{L-1}\left({\hat{v}}_1\right)& \dots & g_{K-1}\left({\hat{\tau }}_p\right)s_{L-1}\left({\hat{v}}_p\right)\end{array}\right).& \left(44\right)\end{array}$

The analytic solution to the least-squares problem (43) is well known [19],

{circumflex over (α)}_p=(Ω_p⁵Ω_p)⁻¹Ω^tvec({circumflex over (H)}).

We now let R_p=Ω_p⁵Ω_p, which is a p×p matrix whose m, nth element is given by

$\begin{array}{cc}\begin{array}{c}{\left(R_p\right)}_{m,n}=\sum _{l=0}^{L-1}\sum _{k=0}^{K-1}g_k^{*}\left({\hat{\tau }}_m\right)s_l^{*}\left({\hat{v}}_m\right)g_k\left({\hat{\tau }}_n\right)s_l\left({\hat{v}}_n\right)\\ ={}^{-j\pi \left(L-1\right)\left({\hat{v}}_n-{\hat{v}}_m\right)T_d}{}^{-\mathrm{j\pi }\left(K-1\right)\left({\hat{\tau }}_n-{\hat{\tau }}_m\right)/T}\left(47\right)\\ \frac{\sin \left(\pi L\left({\hat{v}}_n-{\hat{v}}_m\right)T_d\right)}{\sin \left(\pi \left({\hat{v}}_n-{\hat{v}}_m\right)T_d\right)}\frac{\sin \left(\pi K\left({\hat{\tau }}_n-{\hat{\tau }}_m\right)/T\right)}{\sin \left(\pi ({\hat{\tau }}_n-{\hat{\tau }}_m)/T\right)}\\ \approx \mathrm{KL}{}^{-j\pi \left(L-1\right)\left({\hat{v}}_n-{\hat{v}}_m\right)T_d}{}^{-\mathrm{j\pi }\left(K-1\right)\left({\hat{\tau }}_n-{\hat{\tau }}_m\right)/T}\left(48\right)\\ \sin c\left(\pi L\left({\hat{v}}_n-{\hat{v}}_m\right)T_d\right)\sin c\left(\pi K\left({\hat{\tau }}_n-{\hat{\tau }}_m\right)/T\right).\end{array}& \left(46\right)\end{array}$

and ω_p=Ω^tvec(Ĥ) is a column vector whose mth component is given by

$\begin{array}{cc}{\left(w_p\right)}_m=\sum _{l=0}^{L-1}\sum _{k=0}^{K-1}{\hat{H}}_{l,k}g_k^{*}\left({\hat{\tau }}_m\right)s_l^{*}\left({\hat{v}}_m\right).& \left(49\right)\end{array}$

Claims

1. A method of estimating a time-varying multipath communication channel, the method comprising: receiving a signal comprising a data packet transmitted over the multipath communication channel; and, from the received signalgenerating an estimate of tap parameters in a model of the multipath communication channel, the model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters for estimation defining a delay, a Doppler and a complex gain for the tap,wherein the estimated tap parameters are constant for the data packet.

2. The method of claim 1 wherein the data packet is an orthogonal frequency division multiplex (OFDM) packet.

3. A method of estimating a time-varying multipath communication channel, the method comprising: receiving a frequency-domain channel estimate of the multipath communication channel, the frequency-domain estimate deriving from the data packet transmitted over the multipath communication channel; and, based on the frequency-domain estimategenerating an estimate of tap parameters in a time-domain model of the multipath communication channel, the time-domain model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a complex gain for the tap,wherein the estimated tap parameters are constant for the data packet.

4. The method of claim 3 wherein said generating the estimate of tap parameters comprises: defining a search window for the delay and Doppler tap parameters;dividing the search window into a plurality of bins;determining a correlator metric for the respective bins;identifying one or more peaks in the determined correlator metrics; andoutputting estimates of the delay and Doppler tap parameters based on the identified peaks.

5. The method of claim 3 wherein said generating the estimate of tap parameters comprises: estimating the delay and Doppler tap parameters for a detected multipath component;performing a least-squares estimation to generate an estimate of the gain tap parameter for the detected multipath component; andsubtracting the detected multipath component from the channel estimate, using the estimated delay, Doppler and gain tap parameters to quantify the detected multipath component.

6. The method of claim 5 wherein said estimating, least-squares estimation and subtraction are repeated to iteratively remove detected multipath components from the channel estimate.

7. The method of claim 3 further comprising a method of refining an estimate of the time-domain model of the multipath communication channel, said refining comprising: subtracting all of the plurality of multipath components from the estimate of the multipath communication channel except for a currently-considered multipath component;generating a refined estimate of the delay tap parameter and the Doppler tap parameter for the currently-considered multipath component; andperforming a least-squares estimation to generate a refined estimate of the gain tap parameter for the currently-considered multipath component.

8. The method of claim 7 wherein said refining is repeated to generate refined estimates for the plurality of multipath components.

9. The method according to claim 1 wherein generating the estimate of tap parameters uses an ambiguity function.

10. The method of claim 9 wherein the ambiguity function is described by the following equation:

11. The method of claim 1 further comprising: providing the channel estimate to a radar application for target detection and/or localisation.

12. The method of claim 1 further comprising: providing the channel estimate to a signal communication processor for channel characterisation and/or processing of the received signal.

13. The method of claim 7 further comprising: providing the refined estimates to a radar application for target detection and/or localisation.

14. The method of claim 7 further comprising: providing the refined estimates to a signal communication processor for channel characterisation and/or processing of the received signal.

15. The method of claim 1 further comprising: characterising a reflective environment of the communication channels based on the time-domain model of the multipath communication channel

16. The method of claim 3 further comprising determining the frequency domain estimate, said determining comprising: estimating data from the data packet; andusing the data as pilot symbols to obtain the frequency domain estimate from a signal received over the multipath communication channel.

17. A computer program product comprising machine-readable program code recorded on a non-transitory machine-readable recording medium, for controlling the operation of a data processing apparatus on which the program code executes to perform the method of claim 1.

18. A non-transitory computer readable medium including instructions for: receiving a frequency-domain channel estimate of the multipath communication channel, the frequency-domain estimate deriving from the data packet transmitted over the multipath communication channel; and, based on the frequency-domain estimategenerating an estimate of tap parameters in a time-domain model of the multipath communication channel, the time-domain model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a gain for the tap,wherein the estimated tap parameters are constant for the data packet.

19. The non-transitory computer readable medium of claim 18 further including instructions for: defining a search window for the delay and Doppler tap parameters;dividing the search window into a plurality of bins;determining a correlator metric for the respective bins;identifying one or more peaks in the determined correlator metrics; andoutputting estimates of the delay and Doppler tap parameters based on the identified peaks.

20. The non-transitory computer readable medium of claim 18 further including instructions for: estimating a delay and Doppler tap parameter for a detected multipath component;performing a least-squares estimation to generate an estimate of the gain tap parameter for the detected tap; andsubtracting the detected multipath component from the channel estimate, using the estimated delay, Doppler and gain tap parameters to quantify the detected multipath component.

21. The non-transitory computer readable medium of claim 20 further including instructions for: repeating the estimating, least-squares estimation and subtraction to iteratively remove respective multipath components from the channel estimate.

22. The non-transitory computer readable medium of claim 18 further including instructions for: subtracting all of the plurality of multipath components from the estimate of the multipath communication channel except for a currently-considered multipath component;generating a refined estimate of the delay tap parameter and the Doppler tap parameter for the currently-considered multipath component; andperforming a least-squares estimation to generate a refined estimate of the gain tap parameter for the currently-considered multipath component.

23. The non-transitory computer readable medium of claim 22 further including instructions for: repeating the steps of subtracting, generating a refined estimate and performing a least-squares estimate as defined in claim 16 to generate refined estimates for the plurality of multipath components.

24. A system for estimating a time-varying multipath communication channel, the system comprising: a receiver for receiving a signal comprising a data packet transmitted over the multipath communication channel;a frequency-domain channel estimator that generates a frequency-domain channel estimate based on the received signal; anda time-domain channel modeller that receives the frequency-domain channel estimate and generates an estimate of parameters in a time-domain model of the multipath communication channel, the model characterising the multipath communication channel as a delay line with a plurality of taps, each tap having tap parameters defining a delay, a Doppler and a gain for the tap,wherein the tap parameters are constant over the data packet.

25. An apparatus for estimating a time-varying multipath communication channel based on a signal comprising a plurality of multipath components, the signal associated with a data packet, and based on a frequency-domain channel estimate deriving from the data packet, the apparatus comprising: a time-domain channel modeller that receives the channel estimate and generates an estimate of tap parameters in a time-domain model of the multipath communication channel, the model characterising the multipath components according to the tap parameters, the tap parameters defining a delay, a Doppler and a gain for the multipath components,wherein the tap parameters are constant over the data packet.