A TIME-FREQUENCY APPROACH TO BLIND DECONVOLUTION

IN MULTIPATH UNDERWATER CHANNELS

N. Martins, S. Jesus

∗

SiPLAB–FCT, University of Algarve

Campus de Gambelas, 8000 Faro, Portugal

nmartins@ualg.pt, sjesus@ualg.pt

C. Gervaise, A. Quinquis

EIA Dpt., ENSIETA

29806 Brest Cedex 09, France

gervaice@ensieta.fr, quinquis@ensieta.fr

ABSTRACT

Blind deconvolution is presented in the underwater acoustic

channel context, by time-frequency processing. The acous-

tic propagation environment was modelled as a multipath

propagation channel. For noiseless simulated data, source

signature estimation was performed by a model-based me-

thod. The channel estimate was obtained via a time-fre-

quency formulation of the conventional matched-ﬁlter. Si-

mulations used a ray-tracing physical model, initiated with

at-sea recorded environmental data, in order to produce rea-

listic underwater channel conditions. The quality of the

estimates was 0.793 for the source signal, and close to 1

for the resolved amplitudes and time-delays of the impulse

response. Time-frequency processing has proved to over-

come the typical ill-conditioning of single sensor determi-

nistic deconvolution techniques.

1. INTRODUCTION

Two fundamental problems arising in many signal proces-

sing applications are channel and source signature estima-

tion, often referred as deconvolution. Areas such as data

transmission, reverberation cancellation, seismic signal pro-

cessing, image restoration, etc., are some examples where

deconvolution ﬁnds application. One of these areas is ocean

acoustics, where a fundamental problem is the passive cha-

racterization of acoustic transients and the ocean environ-

ment. This paper aims to give an approach to solve the

problem of blind deconvolution in ocean acoustics, i.e., ob-

taining at once estimates of the emitted signal and the me-

dium impulse response (IR) representative of its physical

properties. In view of the typical structure of underwater

channel model-IRs, constituted by a set of close arrivals,

time-frequency (TF) distributions (TFDs) have been chosen

as a mean for separating the various replicas of the typically

transient source signal. A shallow-water realistic scenario

was chosen to support the presentation.

∗

Thanks to MCT, FCT, PRAXIS XXI (BM/19298/99) for funding, and

to Profs. A. Quinquis and C. Gervaise for the favourable reception at EN-

SIETA.

2. MATHEMATICAL BACKGROUND

One of the main mathematical tools essential in the present

work is the Wigner-Ville distribution (WV), deﬁned by[1]

WV

x

(t, f)=

∞

−∞

x

t +

τ

2

x

∗

t −

τ

2

e

−j2πfτ

dτ,

(1)

where

∗

designates complex conjugate. An also important

TFD is the signal-dependent distribution radially Gaussian

kernel distribution (RGK)[2]

RGK

x

(t, f)=IFT2[Φ

RGK,x

(ν, τ )AF

x

(ν, τ )] , (2)

where Φ

RGK,x

(t, f), the kernel of the distribution, is adap-

ted to the signal to be analyzed, AF

x

(ν, τ ) stands for the

ambiguity function of x(t), and IFT2[·] designates the bi-

-dimensional (2D) inverse Fourier transform operator. Sig-

nal reconstruction was performed via an inverse TF trans-

formation. The basis method[3] was employed, where a

linear combination of the basis functions e

k

(t) was used to

obtain the ﬁnal source signal estimate such as

ˆs(t)=

N

b

k=1

α

k

e

k

(t). (3)

The functions e

k

(t) span the signal space to which the esti-

mate is conﬁned. Channel estimation was performed taking

into account the TF formulation of the matched-ﬁlter (MF).

The fact that the WV obeys to Moyal’s formula[4] conduces

to the following equation:

∞

−∞

WV

x

1

(t, f)WV

∗

x

2

(t+τ ),x

2

(t)

(t, f) dt df =

∞

−∞

x

1

(t)x

∗

2

(t + τ) dt

∞

−∞

x

1

(t)x

∗

2

(t) dt

∗

. (4)

This shows that temporal correlation can be performed in

the TF domain.

II - 12250-7803-7402-9/02/$17.00 ©2002 IEEE

3. SIMULATION RESULTS

The simulated data set corresponds to a canonical two-

-layered shallow-water waveguide whose environment is the

real data acquisition scenario of the INTIMATE ’96 sea

trial[5]. The acoustic source is positioned at 90-m depth,

and the receiver is located at 5.6-km range and 115-m depth,

as illustrated in Fig. 1. The 135-m water column is superim-

Depth

Range

90 m

Source

115 m

Sediment

5.6 km

Receiver

1508

1521

(m/s)

Fig. 1. Scenario of the INTIMATE ’96 sea trial.

posed to a perfectly rigid-modelled bottom. The considered

sound speed proﬁle is shown within Fig. 1. Application of

the ray/beam propagation model BELLHOP[6], taking as

input the environment in Fig. 1, allowed to obtain the refe-

rence channel IR h(t)–Fig. 2. A total number of M =45

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0

0.2

0.4

0.6

0.8

1

Time (s)

Normalized amplitude

Resolved arrivals

Fig. 2. Reference IR h(t) of the canonical scenario. The dashed

line separates unresolved from resolved arrivals.

arrivals have been predicted by the model, spanning the time

interval [0.31, 0.86] s. The amplitudes a

m

and time-delays

τ

m

of the arrivals are the channel parameters to estimate,

grouped into the vectors a and τ , respectively. The source

signal s(t) consists of a T =0.0625-s duration linear fre-

quency modulation (LFM) pulse, spanning the [300, 800]-

Hz instantaneous frequency (IF) range. The IR is split into

2 packets of arrivals: one packet contains the ﬁrst 8 arrivals,

and the other contains the remaining 37 arrivals, which can

be seen as unresolved and resolved arrivals, respectively.

Unresolved arrivals are adjoint arrivals whose arrival times

differ by less than the inverse of the bandwidth of s(t). The

received signal is modelled by

r(t)=x(t)+ξ(t), (5)

where x(t) and ξ(t) are the signal and noise terms, respec-

tively. Next two sub-sections describe results for noiseless

data, whereas Sec. 4 addresses deconvolution robustness to

noise.

3.1. Source Signature Estimation

Considering the multi-component structure of the noiseless

received signal, source signature estimation was performed

by application of the basis method, whose output is given by

(3), and whose model function was deﬁned by the product

of WV

r

(t, f) by a 2D function, centered on the IF estimate

of s(t). This procedure is an approach to solve a mul-

ti-component signal separation problem, conditioned on the

overlap between the very early arrivals.

The information of the purely phase modulated emitted

signal s(t)=exp [jϕ

i

(t)] is contained in the IF f

i

(t), apart

from a phase term, since the signal’s amplitude modulation

component is constant. This implies that most part of the

signal energy concentrates around the line deﬁned by the

points [t, f

i

(t)], in the TF space[7]. This characteristic mo-

tivated source signal estimation, by a ﬁrst estimation of the

IF of s(t). Due to the typical presence of one or a few strong

arrivals in h(t), the IF estimation consisted of the global

maximization with respect to t, of the signal-dependent dis-

tribution RGK of the received signal, RGK

r

(t, f):

[t,

ˆ

f

i

(t)] = { (t, f ):t=arg max

t

RGK

r

(t, f),f∈B}. (6)

In RGK

r

(t, f), for an unitary kernel volume, the ﬁrst 8 ar-

rivals are grouped in a support of large energy, as shown in

Fig. 3. The remaining 37 arrivals, “clustered” in quadruplets,

Fig. 3. Contour plot of the RGK signal-dependent distribution of

the received signal, RGK

r

(t, f).

are clearly distinguished. Maximization of RGK

r

(t, f)

II - 1226

with respect to t, gave rise to an accurate IF estimate[8].

The IF is estimated via maximization of a signal-dependent

TFD, since here, the interference terms (ITs) are quite re-

duced, comparatively to non-signal-dependent TFDs, as the

WV. In that sense, non-signal-dependent TFDs are inappli-

cable in (this) case of signals composed by a sum of repli-

cas, besides the large variance of the estimate, in presence

of noise.

It is likely that the proposed approach will yield a good

IF estimate mainly for the class of purely phase modulated

source signals, for which WV

s

(t, f) is usually concentrated

on a continuous region centered on the line deﬁned by the

pair [t, f

i

(t)]. However, for a non-linear frequency modu-

lation source signal, there are always ITs in its auto-WV,

what may require a smaller value for the kernel volume of

RGK

r

(t, f), if some ITs of that auto-WV have greater am-

plitude than the signal terms. A quadratic IF has been accu-

rately estimated in [9]. The extension to the case of incons-

tant amplitude modulation may hypothetically be treated in

the same manner, provided that the instantaneous amplitude

is a sufﬁciently smooth function, what weakly increases the

spread of the auto-WV along the IF.

For application of the basis method, the considered sig-

nal space was the subspace of band-limited signals, in the

band [0,f

s

/2], which is a natural choice when dealing with

discrete-time signals. The obtained source estimate is de-

picted in Fig. 4, and corresponds to a cross-correlation coef-

ﬁcient of 0.793 with s(t). The band limitation signal cons-

0.25 0.3 0.35 0.4 0.45

−1

−0.5

0

0.5

1

Time (s)

Normalized amplitude

Fig. 4. Source estimate –real part of ˆs(t).

traint is reﬂected not only in the increased signal duration,

but also in an increase of the IF range, since the limitation

to [0,f

s

/2] gave the freedom to the estimate’s IF to extend

in the whole [0, 850] Hz interval.

3.2. Channel Estimation

This sub-section describes the estimation of the parameters

{a

m

,τ

m

; m =1, ..., M} that characterize the channel’s

IR, by use of the TF formulation of the MF mentioned in

Sec. 2. For the case of an LFM signal at emission, and

for the particular cases x

1

(t)=r(t) and x

2

(t)=s(t), (4)

is well approximated by the integration of WV

r

(t, f) along

f

i

(t)[8]. Thus, the IF estimate obtained in Sec. 3.1 was used

to do this “coherent integration” of WV

r

(t, f). The obtai-

ned channel estimate is shown in Fig. 5. For the problem at

0 0.1 0.2 0.3 0.4 0.5 0.6

0

0.2

0.4

0.6

0.8

1

Time (s)

Normalized amplitude

Fig. 5. Channel estimate, obtained by coherent integration of

WV

r

(t, f).

hand, it can be veriﬁed that the obtained resolution is com-

parable to that of the MF result–Fig. 6. It can be shown

0 0.1 0.2 0.3 0.4 0.5 0.6

0

0.2

0.4

0.6

0.8

1

Time (s)

Normalized amplitude

Fig. 6. Channel estimate obtained by classic matched-ﬁltering.

that the relation between the two resolutions depends on the

relative values of the IF range and the inverse of T [8]. In

practice, the ﬁrst 8 arrivals of h(t) are not resolved, and the

quality of the channel estimate will refer only to the remai-

ning 37 arrivals. One normalized measure of the quality can

be given as

ρ

a

=1−

||a −

ˆ

a||

2

; ρ

τ

=1−

||τ −

ˆ

τ ||

2

, (7)

where the vectors must ﬁrst be unitarily normalized. The

quality of the channel estimate is given by ρ

a

=0.9664 and

ρ

τ

=0.9996. The present approach can be applied also to

RGK

r

(t, f).

Alternatively to the procedure of TF coherent integra-

tion, the channel estimate can be obtained by cross-

-correlation between ˆs(t) and r(t), once ˆs(t) be available.

This approach allowed to obtain an estimate very similar to

the previous estimates[8].

4. PERFORMANCE IN NOISE

In this section, the received signal in (5) contains non-zero

white noise ξ(t). Looking at the expected value of

II - 1227

WV

r

(t, f), one ﬁnds that

E[WV

r

(t, f)] = WV

x

(t, f)+σ

2

ξ

, (8)

where σ

2

ξ

is the noise variance. It can readily be seen that

the blind deconvolution method can be performed the same

way as in the noiseless data case of Sec. 3, provided that

WV

r

(t, f) is replaced by its expected value. Regarding

RGK

r

(t, f), one can verify that

E[AF

r

(ν, τ )] = AF

x

(ν, τ )+σ

2

ξ

δ(ν)δ(τ ), (9)

hence not representing a signiﬁcant change on the calculus

of the optimal kernel, relatively to the noiseless case. The

signal-to-noise ratio (SNR) is deﬁned as:

SNR

dB

=10log

1

N

N−1

n=0

s

2

(n)

σ

2

ξ

, (10)

where N is the number of time samples. The empirical mea-

sures below were estimated from 100 trials for each SNR

value. Source estimation performance in noise is illustrated

in Fig. 7, by a plot of the normalized correlation coefﬁcient

as a function of the SNR. Channel estimation performance

−20 −15 −10 −5 0 5 10

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

SNR (dB)

Normalized correlation coefficient

Fig. 7. Source estimation performance in noise.

in noise was measured by the normalized performance mea-

sures (7), for values of SNR ranging from -5 to 10 dB –

Fig. 8. For values of SNR below -5 dB, it was not techni-

cally possible to clearly distinguish the 37 resolved peaks,

with a ﬁxed threshold.

−5 0 5 10

0.992

0.993

0.994

0.995

0.996

0.997

0.998

0.999

1

SNR (dB)

Performance measure

:

amplitudes

time−delays

:

Fig. 8. Channel estimation performance in noise.

5. CONCLUSIONS AND PERSPECTIVES

This work has shown the feasibility of single sensor blind

deconvolution by TF processing, for a realistic multiple time

delay-attenuation channel, and a deterministic non-stationary

source signal. The simulated data set was based on actual

environmental data recorded at sea during a shallow wa-

ter experiment, leading to a particularly difﬁcult test case,

where there are a high number of closely spaced multipaths.

The deconvolved source signal showed a correlation with

the emitted signal close to 0.8, while the channel parameter

estimation was shown to be effective and robust for SNR as

low as -5 dB. Application of the proposed methods on at-sea

recorded data is being performed at this time.

6. REFERENCES

[1] J. Ville, “Th

´

eorie et applications de la notion de signal

analytique,” C

ˆ

ables et Transmissions, vol. 2A, pp. 61–

74, 1948.

[2] R.G. Baraniuk and D.L. Jones, “A radially Gaussian,

signal dependent time-frequency representation,” in

IEEE Int. Conf. Acoust., Speech and Sig. Proc., 1991,

vol. May, pp. 3181–4.

[3] F. Hlawatsch and W. Krattenthaler, “Bilinear signal

synthesis,” IEEE Trans. Sig. Proc., vol. Feb., pp. 352–

63, 1992.

[4] P. Flandrin, “A time-frequency formulation of optimum

detection,” IEEE Trans. Acoust. Speech, Sig. Proc., vol.

36, No. 9, 1988.

[5] X. D

´

emoulin, Y. St

´

ephan, S. Jesus, E. Coelho, and M.B.

Porter, “Intimate96: a shallow water tomography expe-

riment devoted to the study of internal tides,” in Proc.

SWAC’97, 1997.

[6] M. Porter, The Kraken Normal Mode Program, Saclant

Undersea Research Centre, 1991.

[7] L. Cohen, Time-Frequency Analysis, Prentice Hall

PTR, 1995.

[8] N.E. Martins, A time-frequency approach to blind de-

convolution in multipath underwater channels, M.Sc.

thesis, University of Algarve, 2001.

[9] C. Gervaise, A. Quinquis, and N. Martins, “Time-

-frequency approach to the study of underwater acous-

tic channel estimation and source reconstruction,” in

Physics in signal and image processing, Marseille,

France, January 2001.

II - 1228