Signal Processing 11 (1986) 133-143 North-Holland
MULTI-COMPONENT F. G L A N G E A U D
133
FILTERING
OF
MAGNETIC
SIGNALS
a n d C. L A T O M B E
Centre d'Etude des Ph~nom~nes Aldatoires et Grophysiques, UA 346, INPG / ENSIEG, B.P. 46, F-38402 Saint-Martin-d'H~res Codex, France
Received 24 October 1984 Revised 23 November 1985 and 6 March 1986
Abstract. This paper presents two treatments of real signals recorded by an immersed passive array during an underwater magnetic detection experiment. In this situation, the signal-to-noise ratio is very low, and the noise characteristics are close to the characteristics of the signal that is to be detected. Methods that estimate the signal by minimizing the mean-square distance between the signal and its estimate (using a Wiener filter) are optimal, but some a priori knowledge must first be available before these methods can be used. In this specific application, this a priori knowledge can only be acquired by analyzing the data, and a straightforward application of the optimal method may be misleading. Two filtering methods are therefore implemented to estimate a vector valued signal in a noisy observation. These filtering methods are simple and robust, and they give good results on real data. One method is an adaptation of Wiener filtering. The other consists of global filtering using the eigen-structure of the spectral density matrix of the received signals. Results for each method are presented and compared. Zusammenfassung. Fiir ein Mehrkomponentensignal, das durch additives Rauschen gestrrt ist, werden zwei Methoden der Filterung beschrieben. Hierbei wird zum einen eine Kettenschaltung mehrerer Wiener-Filter angewendet; zum anderen wird das Signal global mit Hilfe der Eigenwerte und -vektoren der Spektralmatrix gefiltert. Die theoretischen Grundlagen beider Methoden werden dargelegt; auBerdem werden praktische Ergebnisse fiir Unterwasseraufzeichnungen magnetiseher Signale diskutiert. Rrsumr. Cot article prrsente deux traitements appliqurs h des signaux rrels qui ont 6t6 enregistrrs par une antenne passive immergre dans une exprrience de drtection de signaux magn&iques en mer. Dans cette situation prrcise, le rapport signal bruit est trrs faible et la forme du signal bruit pout ~tre trrs proche de cello du signal h estimer et ~ drtecter. Los mrthodes estimant le signal en minimisant la distance quadratique entre le signal et son estimre (filtrage de Wiener simple) sont optimales mais nrcessitent une connaissance a priori sur le signal. Dans l'application particuli~re traitre, cette connasissance a priori ne pout qu'~tre estimre h partir des donnres et une application brutale de la mrthode optimale pout alors contuire une solution erronnre. Nous implantons ici deux mrthodes de filtrage pour estimer un signal vectoriel dans une observation bruitre. Cos mrthodes sont simples et robustes et donnent des rrsultats intrressants sur cos donnres rrelles. Ce sont, d'une part une adaptation du filtrage de Wiener, et d'autre part un filtrage global utilisant los 616ments propres de la matrice spectrale des signaux re~us. Los rrsultats de chaque mrthode sont prrsentrs et comparrs. Keywords. N-dimensional signal, multicomponent filtering, Wiener filtering, spectral matrix.
can vary using diving wings. Three parameters
1. Introduction
describe the spatial position of the towed vehicle: v e h i c l e c a n b e u s e d to
t h e a n g l e s o f list a n d t r i m , a n d the i m m e r s i o n
r e c o r d t h e s u r r o u n d i n g m a g n e t i c field. T h e v e h i c l e h a s a d e l t a - s h a p e d tail w i n g t h a t c a n b e c o n s i d e r e d
( d i s t a n c e to t h e s u r f a c e ) . A c q u i s i t i o n a n d p r e p r o c e s s i n g a r e d e s c r i b e d in m o r e d e t a i l in [5].
as p l a n a r . A t t h e e x t r e m i t i e s o f t h e tail w i n g , t w o
From the observed data, a vector valued signal
A towed
underwater
magnetometers
measure
the
horizontal
com-
p o n e n t s H1 a n d H2 o f t h e m a g n e t i c field. D e p t h
x is c o n s t r u c t e d ; its c o m p o n e n t s are, in o r d e r : [ H i - / - / 2 , H1, list, t r i m , i m m e r s i o n ] .
0165-1684/86/$3.50 © 1986, Elsevier Science Publishers B.V. (North-Holland)
F. Glangeaud, C. LaWmbe / Filtering of magnetic signals
134
mqgnei~me~rs
densities of the desired signal and of the noise are known. An adaptation of this method is presented in the case where the spectral densities are not known but are estimated from a second observation.
2.1, Wiener filtering using two scalar observations f
~
divingwings
~
X ~ : Trim y ~ z : List
Fig. 1. Immersed vehicle; z is the vertical going through the
gravity center O, y is perpendicular to the wing plane, and x is perpendicular to the vehicle axis, in the vertical plane. Notice that only the first two components are field dependent. The experiment is designed to obtain a map of the magnetic field, thereby obtaining the best possible estimates of H1 and H2. When a magnetic body is present, its magnetic signature is a characteristic of this body. However, movements of the immersed vehicle induce parasitic noises on the magnetic records. These noises closely resemble the desired magnetic 'signature'. In addition to this ambiguity between the useful signal and parasites, the signal-to-noise ratio is low, and uncorrelated additive noises may occur on each component. Assuming that each component is a scalar signal that can be considered as the output of a linear and time-invariant filter excited by the signal represented by an other component, this study attempts to identify these filters and to use them to select the useful part of the observations (here the real magnetic signal). The proposed techniques are: - Wiener filtering of pairs of scalar signals either in parallel or in cascade form, - global vector filtering using the projection on the eigenvector associated with the largest eigenvalue of the spectral matrix of the recorded signals.
2. Wiener-filtering o f an N-component vector signal
Two observations x~(t) and x2(t) are considered: Xl(t) = sl(t) + p l ( t ) ,
s2(t)+p2(t), where signals sl(t) and SE(t ) are uncorrelated with the parasitic signals pi(t). XE(t)
A relationship may exist either between p~ and P2, or between sl and s2. This relationship will be assumed to admit a linear time-invariant filter model ~. In practice, either the frequency response F ( f ) of the filter ~, or the signals si(t), or the parasitic p~(t) are of interest. For reasons of simplicity, two particular cases are commonly considered: xl(t) = sl(t)+pl(t), x2( t) = s2( t),
(1)
where s2(t) is the output of filter ~ with input sl(t), x2(t) is called a signal reference. x](t) = s](t) +pl(t), x2(t) =
p2(t),
(2)
where p](t) is the output of filter ~ with input p2(t), x2(t) is called a noise reference. The magnetic signals considered hereafter are of type (2): the useful signal only emerges from magnetic sensors whereas the position sensors (list, trim, and immersion) are sensitive to parasites only. To estimate the linear time-invariant filter from the noisy observations x](t) and x2(t), a linear time invariant flter (~ is constructed so that, when its input is x2(t), the output is as 'near' as possible to the parasitic signal pl(t). The filter (g has an impulse response g(t) and a frequency response
G(f) (gG~)~ G(f)). Classical Wiener filtering [7, p. 336] extracts the desired signal from a noisy observation if spectral Signal Processing
Consider the diagram in Fig. 2, where * indicates convolution.
135
F. Olangeaud, C. Latombe/ Filteringof magnetic signals one would have
X I I f ) = sl{f)+pl(f) U ( f } = X 1(t) - X 2 { t ) ~ g ( t )
G ( f ) = F(f), p,(t) =x2(t) * g(t) =p2(t) * g(t),
x2{t) =p2(t) u(t) = xl(t)
Fig. 2.
- x2(t) * g(t)
= sl ( t ) .
An estimate of the frequency response G(f) is The equivalent expression, in the frequency domain is
^ ....
d(s)= oUJL T A
U(f) = X~(f) - X2(f) G(f), where Xi(f) denotes the Fourier transform of the scalar signal xi(t). U(f) is the difference between the observation X](f) and an element of the one-dimensional vector space spanned by X2(f). The space of stochastic stationary processes, with the inner product defined as (X, Y)= E{XY*} (E denotes expectation and y* the complex conjugate of y) is a Hilbert space. According to the projection theorem [10] in Hilbert spaces, a minimum norm element U(f) exists, that is orthogonal to the vector space spanned by X2(f), i.e., E{ U(f)X*2 (f)} = O,
E{ X , ( f ) X * ( f ) } = G(f)E{ X2(f)X*2 (f)}. Hence (cf. [10]), the frequency response of the noncausal Wiener filter is given as
F(f)
-
E{X,(f)X*2(f)}
T~2(f)
~,{x~(f)x~(/)} v~:(S)
['yl,(J)l ''~
= c(f)L~
1
,
where ~/.(f), T22(f), and 712(f) are the spectral (respectively cross spectral) densities of observations 1 and 2, and C ( f ) is the complex coherence function. In this study, the parasitic pi are uncorrelated with the signals s~, whence 3',2(f) = E { X , ( f ) X * ( f ) }
= E{P,(f)P*2(f)}. In the case of perfect modelling and estimation,
F _ A, [ ~ ' ~ ' l l / 2 / T11kJ 1/
A
/
,
A
where C(f), y]~(f), and Y22(f) are estimated from the data. However, some difficulties arise in the A implementation of the noncausal filter G(f), mainly due to the difficulty of obtaining good estimates of the needed spectral functions C(f), y , ( f ) , and Y22(f). These difficulties will be outlined below, and an explanation will be given showing how they can be overcome partially. 2.2 Implementation di~iculties The computed Wiener filter (~(f) explicitly depends on the spectral densities and on the coherence function. These quantities are statistical ones and are defined in terms of ensemble averages. These ensemble averages (under the hypotheses of stationarity and of ergodicity) can be estimated from a single sample with a smoothed periodogram (i.e., convolving the sample periodogram with a spectral window W(f)). The characteristic parameter of the method is the product of the equivalent bandwidth of the window by the time of integration; this product is denoted by BwT and can be expressed in terms of the Fourier coefficients of the window (see [6, 7]). Nevertheless, this smoothing can be misleading when narrow band spectral densities are present: in this case, spectral lines are spread over a frequency band as wide as the equivalent bandwidth of the window. Hence, one proceeds by steps: first, narrow spectral windows are chosen in order to detect possible spectral lines. Wider windows are then used. When the coherency coefficiency is low, its estimation is highly biased (see [6]). In [3], it is shown Vol. I1, No. 2, September 1986
F. Glangeaud, C. Latombe / Filtering of magnetic signals
136
that the expectation of the error on G ( f ) satisfies
p ( f ) defined by p ( f ) = Co(f) * W ( f ) , where W ( f ) is the spectral window used in the estimation, and
E { I G ( f ) - G(f)l} IG(T)I 2 1~
1 ( 1
Co(f)={~
- n - ~ \ l f ( f ) l =- /" Hence, if the coherency coefficient is not close to 1, a good estimate of G ( f ) cannot be obtained, and in this case, the estimated filter G ( f ) may degrade the useful part of the signal. To prevent this damaging effect, Comon and Lacoume [3] set a limit L to the coherency coefficient: when Ic(f)l is lower than L, it has to be considered as being zero. This limit is computed by evaluating the ratio of the signal to the extra noise due to estimation errors of the filter and it depends explicitly on the product BwT, L=
if l C ( f ) l > L, otherwise.
In Fig. 3(a), the coherency coefficient is presented, as well as the lowest limit L (0.47) and the correcting function p ( f ) . Fig. 3(b) compares the moduli of Wiener filters estimated with or without the correction: at frequency 0.037 Hz, the modulus of the frequency response G ( f ) (large because T22(f) is low) is reduced when the correction p ( f ) is used.
2.3. Cascade Wiener filtering For our specific underwater detection and estimation case, the observation can be modelled as follows:
1
I+BwT"
xl(t) = sl(t)+pl(t),
Nevertheless, if all the values of IC(f)l that are lower than the limit are set to zero, sharp cut-off filters are introduced. When the level is exceeded only occasionally, the filter will be selective at such cross-overs, too. To reduce this effect, in this study, the coherency is convolved by a correcting function
x2(t) = s2(t) + p2(t),
Xa(t) = p3(t), x~(t) = p,(t),
x5(t) = ps(t),
i
0.47
0
0.3Hz
(a)
1.04
/[,Wiener fiLher
0.69 ~
/.I
B,,T : 7
[on'ecfecl W filter
/'~k.~ _--=A_/~ 0.017
0.3ttz
A FREQUENEY (b) Fig. 3. Wiener filter spoils the signal when coherency is lower than the threshold 0.47, thus a correction p(F) is applied, p(F) has no sharp cut-off. Signal Processing
F. Glangeaud, C. Latombe / Filteringof magnetic signals where the parasitic noises pl(t) and p2(t) are mixtures of filtered versions of p3(t), p4(t), and ps(t). Moreover, some additive uncorrelated noises hi(t) may be recorded on each sensor. To take advantage of all the noise references, the following steps were chosen, in this specific problem: - All the coherence functions between xi(t) and xt(t) ( i = 1, 2, j = 3 , 4, 5) are computed. - The most coherent pair [(x~, Xto), io = 1, 2,jo = 3, 4, 5] is chosen and the estimation of the parasitic /~o(t) is determined using a noise-reference Wiener filter. - Coherence functions between x ~ ( t ) - ~ ( t ) and xj(t) (j =3, 4, 5 and j ~jo) are computed and a Wiener filtering is performed for the most coherent pair. If the estimate /~(t) is close to the parastitic p~(t), so is the estimate ~ ( t ) to the useful signal sg(t). - A last Wiener filtering of x~(t) = s~(t)+p~(t) (i = 1, 2 and i # io) with ~io(t) as a signal reference provides the desired estimate for s~(t).
3. Filtering by projection onto the dominant eigenvector o f the spectral matrix
Filtering by using the eigen-elements of the spectral matrix has become a well-known theoretical topic [1, 4, 8, 14]. This method of filtering has been applied on real data in underwater acoustics [1] and in geophysics [4, 9]. The general backgrounds will be briefly reviewed below. Consider a passive array with N sensors that receives an incoming wave-field. This field is generated at a few specific points in space by P stochastic excitations, called sources, that are generally assumed to be uncorrelated. It is assumed that the medium of propagation is linear and time invariant. The sources are characterized by: aJ(t) : amplitude of the j t h source signal that is assumed to be ergodic and stationary; its Fourier transform is denoted
AJ(f),
137
s~(t) : stochastic signal radiated from the j t h source which could be recorded at the kth sensor if the jth source were alone, hJk(f) : the transfer function between source j and sensor k. The vector whose components are the Fourier transforms of s~(t) is denoted SJ(f). The vector whose components are h~(f) is denoted H i ( f ) . The vector sJ(f)= Ht(f)/llHt(f)ll has a norm equal to 1 (1[. II denotes the Euclidean norm of a vector). It describes the wave-front in the vicinity of the sensors and can be related to the source itself only if the properties of the medium between source and receiver are perfectly known. Nevertheless, it is called the 'directional source vector'. The preceding assumptions imply S t ( f ) = ,4~(f)st(f),
.~J(f)=AJ(f)/tlHJ(f)ll
where is the apparent amplitude (in frequency domain) of the jth source: for reason of simplicity we still use the notation AJ(f) in the following. The usual hypotheses are made: - T h e sources are assumed to be statistically uncorrelated (E{At,(f)A*'2(f)}=O; Jl #j2, E denoting expectation) and their directional vectors are assumed to be linearly indep,endant (the subspace spanned by { s l ( f ) , . . . , s P ( f ) } is of dimension P). - The medium of propagation between the sensors is assumed to be linear, time-invariant, and homogeneous. - T h e recorded data are corrupted by additive noise sources that are assumed to generate white zero-mean noises. The white zero-mean noises are uncorrelated with each other and uncorrelated with the signal. They are isotropic and of variance No. These noises are denoted by the vector notation b(f). The vector recorded signal may be written as P
X ( f ) = E S t ( f ) + b(f) j=l P
= ~, At(f)sJ(f) + bt(f). i=1 Vol. 11, No. 2, September 1986
F. Glangeaud, C. Latombe / Filtering of magnetic signals
138
Under these hypotheses, the N x N spectral matrix of the data is F ( f ) = E{X(f)X+(f)} P
= ~ ~(f)sJ(f)sJ(f)++Nol, j=l
where dj(f)= E{IA(T)I 2} is the apparent spectral energy of the j t h source and ÷ denotes the conjugate transposed. If )tl(f) t> a2(f) ~>" • • i> AN(f) are the eigenvalues of F ( f ) and { V ~ ( f ) , . . . , VN(f)} are their associated vectors, a known property [1, 4, 8, 14] is that Ap+~= Ap+2. . . . . AN = No. Hence, the number of sources is obtained from the multiplicity of the minimal eigenvalue of the true spectral matrix. However, in practice, only the estimate /~(f) can be obtained from the data and the number of sources must be estimated from the eigenvalues Ai(f) of /~(f). It has been shown [8] that this estimation highly depends on the "equivalent bandwidth x time of integration product" (denoted BwT) of the method used to estimate the spectral matrix of the records. Nevertheless, some estimation criteria are now well known and can be classified into two groups: -criteria comparing the eigenvalues to a fixed level [1, 8, 11], -statistical criteria obtained from information theoretic principles such as Akaike's or Rissanen's criteria [11, 13]. When a single source is detected, we know [1, 4, 9] that its normalized directional vector s~(f) is equal to the eigenvector Vt(f) associated with the dominant eigenvalue A~(f):
X (f) = S l ( f ) -+'b(f) : A ~(f) Vl(f) + b(f). If the noise level is low enough, the amplitude of this signal source is
A~(f) = (X(f), Vm(f)), where (., .) denotes the scalar product of vectors and where the residual noise in the direction of Vm(f), denoted (b(f), Vl(f)), is neglected. Thus, all the complex components of the source vector S~(f) = Am(f) V~(f) can be identified. Signal Processing
In this single source case, the first eigenvector must be considered as describing the wave-front surrounding the sensors at each frequency. The modulus of the ith component gives the relative influence of the source on sensor i (a modulus equal to 1 means that the source reaches this sensor only) whereas the phase of component i gives the phase difference between sensor i and sensor 1 (taken here as the phase reference).
Remark Prior to the processing, the recorded components have to be normalized so that the noises have nearly the same energies. One possible procedure is to normalize each record by the inverse of its standard deviation. The treatment can be performed, and finally the signal obtained after projection can be remultiplied. Moreover, if the data have different physical natures, projections keep the same nature.
4. Results for the magnetic underwater recordings The present data (Fig. 4) were recorded at sea, 60 meters below sea level. Of the five records, the last three ones (list, trim, and immersion) are noise references. The magnetic signal is present on the first two records only but it is very weak. The record /-/1 -/-/2 will be denoted by grad H.
4.1. Cascade Wiener filtering When the cascade Wiener filtering is applied, the strategy presented in Section 2.3 leads to the following choices, after coherencies are computed: (a) the pair (grad H, trim) is selected, with trim as noise reference, to estimate the parisite related to trim and contained in grad H. The signal obtained by subtracting this parasite from grad H is denoted grad~trim. (b) The last signal is filtered using immersion as a noise reference, to estimate the parasite related to immersion and contained in grad/trim. The final signal obtained by subtracting this estimate from grad/trim is denoted grad/trim/immersion. The
F Glangeaud, C. Latombe / Filtering of magnetic signals
3~-| (]RAO H
It can be seen in Fig. 4 that the high intensity and low frequency perturbation that occurs on noise references near the middle of the time duration influences the signal plus noise records. One way to evaluate results of the cascade Wiener filtering, is to consider the standard deviation (denoted or) of each signal. Assuming that each record is a sample of a stationary and ergodic stochastic process, the best estimate of the low level magnetic signal would be the one having the weakest standard deviation. We have found (see Fig. 5):
A
Signal 4-
Noise
Noise
7.?°
TR,M
Ii
/
.
IMMERSION
0
TIME
for grad H:
t7 = 0.66,
for grad/trim:
tr -- 0.24,
for grad/trim/immersion: tr -- 0.20.
References 6.3m -
139
~52s
Fig. 4. D a t a w a v e forms.
most logical procedure would be to extract the part correlated to trim from immersion, and to use it as a noise reference to filter grad/trim. However, in order to avoid an additional filtering, this procedure is not followed because the signal-to-noise ratio is already high enough. (c) The magnetic signal H, is filtered using grad/trim/immersion as a signal reference to obtain the part of the magnetic signal related to the gradient component and from which perturbations related to movement are cleared. Here again the most logical procedure would be to estimate the part of H, rid of the parasites related to trim and immersion and to filter it with grad/trim/immersion as a signal reference. To avoid two more filterings and because the present results are acceptable, this procedure is not followed.
Hence, the first filtering (step (a)) has significantly improved the estimate of the magnetic signal contained in grad H. Step (b), although useful, is less significant. Fig. 6 outlines results of step (a): the shape of the modulus of the estimated Wiener filter (denoted WF) is compared to the coherency coefficient and the spectral density of grad/trim is seen to be 12 dB lower in the low-frequency band than the spectral density of grad H. Finally, Fig. 7 compares the signals obtained after the three steps of the cascade Wiener filtering to the signals obtained after a classical Wiener filtering between H, and grad H. The standard deviations of the resulting signals are, respectively, tr -- 0.19 and (7 = 0.78, hence proving the improvement that was performed. In fact, the data presented are for a case where no magnetic signature occurs. The magnetic signal resembles a noise.
4.2. Filtering using the first eigenvector The a priori hypothesis that a relationship between all the components exists is made after inspecting the presented records. If the hypothesis is justified, the method using the eigenelements of the spectral matrix may be fruitful. Vol. 11, NO. 2, September 1986
F. Glangeaud,C. Latombe/ Filteringof magneticsignals
140
Chained Wiener filleting
To estimate the number of excitations, three criteria were applied on these data. Their results are presented in Fig. 8. The first criterion [11], denoted CSS, is a simple criterion that compares the eigenvalues to the fixed level
0-=0.66
ORAD H
1
1
N
fj Ni=l
S=--max-NOISE REF:TRIM
0-=I.53 1
I I q ORAD/TRIM
o- =0.2/, ---1
NOISE REF:IMM
~-
o-- =172
TIME
B52E
Fig. 5. Steps(a) and (b) of the cascade Wienerfilteringproviding the final signal referencegrad/trim/immersion.
First, the records are normalized following the remark of Section 3. Each element of the spectral matrix of the data is obtained with the smoothed periodogram method [6, 7], using a smoothing spectral window that has the product Bw T equal to 10. Eigenvalues of this spectral matrix are computed with classical subroutines that implement tridiagonalisation and Rutishauser methods. Results are presented in the low frequency band, from 0 to 0.15 Hz. Signal Processing
N
Y~ Ai(fj)
(Fig. 8(b)). The second criterion uses explicitly [12] the chi-square distribution of the likelihood ratio of observations (Fig. 8(c)). The last one [11, 13] is the Rissanen (or MDL) criterion (Fig. 8(d)). Using the second detection criterion, we say that a single excitation exists in the frequency bands [0.01 ; 0.03 Hz] and [0.05 ; 0.07 Hz] (of. Fig. 8). This single excitation or 'source' is hence described by the eigenvector associated with the dominant eigenvalue (Fig. 9). It can be seen on the moduli of the various components that: - the source energy is nearly equally distributed among components in the frequency band [0.01 ; 0.03 Hz] meaning that, in this frequency band, the parasite also affects the magnetic records, justifying our a priori hypothesis, - t h e modulus of the third component of the 'source' vector is close to one in the frequency band [0.05;0.07Hz]. This means that the 'source' is in fact a noise source that affects the trim record and has no influence on the other components. The 'source' vector S l ( f ) = (X, VOsl(f) is identified only for low frequencies (0.01 <~f<~ 0.03 Hz) (of. Fig. 10, traces (b)), and is considered as the estimated parasite. It coincides rather well with the data. On each component, differences between data and estimated parasite are computed. They are to be compared with signals in Fig. 7. It can be noted that, if the results for the grad H component are similar, the result obtained in estimating the Ht component seems better when using the Wiener cascade filtering. This is due to the fact that the one-source decision in this frequency band is not correct and that in reality a low level second source exists whose magnetic components are significant.
F. Glangeaud, C. Latombe / Filtering of magnetic signals
J ORAD
~
--~ ~ -
cr =0.66
I ~ i
-/-~-M".~ \,",,',-,,-.,.
!
141
IClf)l Ck/~
~1 /" /i~''~
i"~JI\v,.~i/~-I ~ !.,r~ it~~, +~
I
0 PARTOFGRAD RELATEDTOTRIH
IE-"~._L~,.,
cr =0 59
/\
SPErrRALD~NSrrlES'I
®
I . ~, .I~TRI G HyR %A D
o- =0.2~+
0
FREQUENCY
0,15Hz
-'~.~-~,~
0
TIHE
~52.,
Fig. 6. Part of the gradient signal related with trim is filtered. The result is called gradient/trim. Wiener filter WF is presented in modulus en phase.
Single Wiener filtering I H1
-
{ step ¢ }
o-=1.0/+/
.,
, l~_,~, SISNALREF:GRAB
ReferenceJ,/
£oscode Wiener filtering
- / o-=0.66 J
,~
~
HI
-
-
o- =I .04
~
StONALREF:SRAD/TRIH/IHM
o- =0.20
,i
Result
Result
o- =019
[_
0-
'
TI.E
852~
o
T~ME
8~2~
Fig. 7. Estimations of the magnetic signal on H~. Vol. 11, No. 2, September 1986
142
F. Glangeaud, C. Latombe / Filtering of magnetic signals
Eigen vector Orndient (o)
_ Odb
~ ....... L.._L_~'~z-~=.. Jr-~_ _ __1 ~..~_.~_~
CSS
-8db (b)
'
L
CHI SQUARETEST (c)
I o
,
0
Trim
mE
MDL,2
'°'IV2,0 0
0.03Hz
0.15Hz FREQUENCY
Fig. 8. (a) Estimated eigenvalues. (b) Estimated eigenvalues (logarithmic scale) and the fixed level S. (c) Estimated number of sources using the 'chi-square test'. (d) Estimated number of sources using the Rissanen (MDL) criterion.
5. Conclusion Two procedures' for vector-valued signals have been presented. The procedures were conducted on magnetic underwater signals. Cascade Wiener filtering gives good results in the specific application presented here. A preprocessing must be carded out during which the coherencies for each pair of signals are evaluated. Global filtering, which identifies the useful signal after projecting the recorded signal onto the first eigenvector of the spectral matrix, also performs well when a single source is detected. However, its Signal Processing
0.1 O.3 0.5 0.7
0.15
FREQUENCY
Fig. 9. Moduli of the fivecomponents of the eigenvectorassociated with the dominant eigenvalue.
application is limited to the single source case, when the assumption of planar waves is not fulfilled. The proposed Wiener cascade filtering that uses the same subroutine several times is simple to implement. In this specific application a reasonable estimate is obtained in three iterations. Other treatments o f similar records have shown that, due to the introduction of a minimum level for the coherency, the procedure always converges towards a solution for which estimation errors can be evaluated, whereas this is not the case in the classical multicomponent Wiener filtering [3].
F. Glangeaud, C. Latombe / Filtering of magnetic signals
143
References
(a}
I
LOW_FREQ.ESTIMATEDPARASITE o-=05B
(d-(b)
a-=O.1B
Jlc) ~
a---1.0~ I
r!
A
,
a) LOW. FREQ.ESTIMATED PARASITE o-:0.66
[a)-(b)
o-=0.62
Fig. 10. Estimates of the low frequency magnetic signal using global filtering.
Acknowledgment The authors are indebted to the DCAN/ BREST/GESMA, which provided the pretreated records. They also would like to thank Mr. Pretet who proposed this study and provided useful advices and P. Lorenzino who helped in the computing work.
[1] G. Bienvenu and L. Kopp, "Optimality of high resolution array processing using the eigensystem approach", IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-31, No. 5, October 1983, pp. 1235-1247. [2] D.R. Brillinger, Time Series Data Analysis and Theory, Holden-Day, San Francisco, CA, 1981, Chap. 8. [3] P. Comon and J.L. Lacoume, "Noise reduction for an estimated filter using noise references", IEEE Trans. Information •,Theory, Vol. IT-29, March 1986, to appear. [4] F. Glangeaud and C. Latombe, "Identification of electromagnetic sources", Ann. Geophysicae, Vol. 1, No. 3, May/June 1983, pp. 245-252. [5] P. Granjean and G. Pretet, "Drtection magnrtique en mer: Etude de l'amrlioration des mrthodes de drtection et de localisation", Proc. Collque GRETSI, Nice, June 1975, pp. 531-538. [6] G.M. Jenkins and D.G. Watts, Spectral Analysis and Its Applications, Hoiden-Day, San Francisco, CA, 1968. [7] J.L. Lacoume and C. Latombe, "Analyse interspectrale", Proc. Colloque GRETSI, Nice, June 1981, pp. 311-318. [8] C. Latombe, "Nonconventional array treatment using the eigensystem of the spectral matrix", in: Schussler, ed., Proc. EUSIPCO '83, Erlangen, FRG, September 1983, pp. 499-502. [9] C. Latombe, F. Glangeaud and C. Lathuillere, "Traitement matriciel de signaux & N composantes et 6tude par couple", Proc. Colloque GRETSI, Vol. 1, Nice, May 1983, pp. 253-258. [10] A. Papoulis, Signal Analysis, McGraw-Hill, New York, 1977. [11] I. Tas and C. Latombe, "Performance de divers algorithmes pour la drtection, dans les probl~mes d'imagerie, du nombre d'excitations non corrrlres", Proc. Colloque GRETSI, Vol. 3, Nice, May 1985, pp. 363-367. [12] I. Tas and C. Latombe, "Drtection multiple par les valeurs propres de la matrice spectrale", Traitement du Signal, September 1986, to appear. [13] M. Wax and T. Kailath, "Determining the number of signals by information theoretic criteria", ASSP Spectrum Estimation WORKSHOP II, Tampa, FL, November 1983, pp. 192-196. [14] M. Wax, T. Shan and T. Kailath, "Spatio-temporal spectral analysis by eigenstructure methods", IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-32, No. 4, August 1984, pp. 817-827.
Vol. 11, No. 2, September 1986