View - APC

15 Pages·2001·0.96 MB·English

PHYSICALREVIEWD,VOLUME63,042004
Adaptive filtering techniques for gravitational wave interferometric data:
Removing long-term sinusoidal disturbances and oscillatory transients
E. Chassande-Mottin
Max-Planck-Institutfu¨rGravitationsphysik,Albert-Einstein-Institut,AmMu¨hlenberg1,D-14476Golm,Germany
S. V. Dhurandhar
Max-Planck-Institutfu¨rGravitationsphysik,Albert-Einstein-Institut,AmMu¨hlenberg1,D-14476Golm,Germany
andIUCAA,Postbag4,Ganeshkhind,Pune411007,India
~Received27March2000;published29January2001!
Itisknownbytheexperiencegainedfromthegravitationalwavedetectorprototypesthattheinterferometric
outputsignalwillbecorruptedbyasignificantamountofnon-Gaussiannoise,alargepartofitbeingessen-
tially composed of long-term sinusoids with a slowly varying envelope ~such as violin resonances in the
suspensions,ormainpowerharmonics!andshort-termringdownnoise~whichmayemanatefromservocontrol
systems, electronics in a nonlinear state, etc.!. Since non-Gaussian noise components make the detection and
estimationofthegravitationalwavesignaturemoredifficult,adenoisingalgorithmbasedonadaptivefiltering
techniques ~LMS methods! is proposed to separate and extract them from the stationary and Gaussian back-
groundnoise.Thestrengthofthemethodisthatitdoesnotrequireanyprecisemodelontheobserveddata:the
signals are distinguished on the basis of their autocorrelation time. We believe that the robustness and sim-
plicityofthismethodmakeitusefulfordatapreparationandfortheunderstandingofthefirstinterferometric
data.Wepresentthedetailedstructureofthealgorithmanditsapplicationtobothsimulateddataandrealdata
fromtheLIGO40mprototype.
DOI:10.1103/PhysRevD.63.042004 PACSnumber~s!: 04.80.Nn,07.05.Kf,07.50.Hp,07.60.Ly
I.INTRODUCTION instruments! consisting mainly of noise from a host of
sources both environmental and intrinsic to the apparatus.
Over the next decade, several large-scale interferometric Buried in this noise will be the gravitational wave signature.
gravitational wave detectors will come on-line. These in- Sophisticateddataanalysistechniqueswillhavetobedevel-
clude LIGO, composed of two Laser Interferometer oped to optimally extract astrophysical data. Many of the
Gravitational-wave Observatories situated in the U.S. @1#, techniques developed so far @7–9# are based on matched fil-
VIRGO, a French-Italian project located near Pisa @2#, tering and assume stationary Gaussian noise.
GEO600, a German-British interferometer under construc- However, the real data stream from the detectors is not
tion near Hannover @3#, TAMA in Japan, a medium-scale expected to satisfy the stationary and Gaussian assumptions.
laser interferometer @4#, and with funding approval In fact, the data from the Caltech 40 meter proto-type inter-
AIGO500, the proposed 500 meter project sponsored by ferometer has the expected broadband noise spectrum, but
ACIGA. There are also separate proposals for space-based superposed on this are several other noise features @7#; such
detectors which could be operational twenty-five years from aslong-termsinusoidaldisturbancesemanatingfromsuspen-
[email protected].,theLaserInterferometerSpaceAntenna~LISA!,a sions and electric main harmonics and also transients occur-
cornerstone project of the European Space Agency @5##. In ringoccasionally,typicallyduetoservo-controlsinstabilities
the meantime, a number of existing resonant bar detectors or mechanical relaxation in suspension system etc. While no
will have had their sensitivities further enhanced. precise a priori model can be given for this noise until the
Thekeytogravitationalwavedetectionistheveryprecise detectoriscompletedandfullytested,matchedfilteringtech-
measurementofsmallchangesindistance.Forlaserinterfer- niques cannot be used to locate or remove these noisy sig-
ometers,thisisthedistancebetweenpairsofmirrorshanging nals.
at either end of two long, mutually perpendicular vacuum This disparity between standard Gaussian assumptions
chambers. Gravitational waves passing through the instru- and real data characteristics poses a major problem to the
ment will shorten one arm while lengthening the other. By direct application of matched filtering techniques. This is
using an interferometer design, the relative change in length true when searching for burst sources such as blackhole bi-
of the two arms can be measured, thus signaling the passage nary quasinormal ringings @10#. This is also the case for the
ofagravitationalwaveatthedetectorsite.Longarmlengths, inspiral searches in Caltech 40 m data, where one has to
high laser power, and extremely well-controlled laser stabil- introduce a veto @7# on the decision taken with the matched
ity are essential to reach the requisite sensitivity, since the filtertoensurethatthedetectedsignalisactuallytheonewe
gravitational waves will be faint and will modify only are looking for.
weaklythestructureofspace-timeinthedetector’sarms~see It is possible that in the future, improved experimental
e.g., @6#!. techniques and greater experience, will reduce or even com-
Gravitational wave detectors produce an enormous vol- pletely eliminate some of these nonstationary and non-
umeofoutput~e.g.,oftheorderof16MB/secfortheLIGO Gaussian features. Nevertheless, it will take probably some
0556-2821/2001/63~4!/042004~15!/$15.00 63042004-1 ©2001TheAmericanPhysicalSociety
E.CHASSANDE-MOTTINANDS.V.DHURANDHAR PHYSICALREVIEWD 63042004
time to reach such acceptable and high quality of data. define our variables and therefore stands for ‘‘equal by defi-
Therefore, it is necessary and desirable to somehow combat nition.’’ Finally, the vectors xt, ¯x and x†[¯xt denote respec-
this noise. Since such noise features defy modeling, a novel tivelythetransposeofthevector,thecomplexconjugateand
approach to the problem is called for. the Hermitian transpose of the vector x.
We propose a denoising method based on least mean
square (LMS) adaptive linear prediction techniques which
C.Meansquarelinearprediction
does not require any precise a priori information about the
noise characteristics. Although our method does not pretend First, we recall some standard principles to design an op-
to optimality, we believe that its simplicity makes it useful timallinearpredictor.Thequestiontoaddressistooptimally
for data preparation and for the understanding of the first predict the data sample xk with a collection of past samples
data. taken between the indexes k2d2(N21) and k2d. These
In the following, we present the principles of LMS adap- past samples may be placed in a column vector xk whose N
tive denoising ~Sec. II!, a characterization of its behavior on components are given by
asinglemodelofthenoisefromtheinterferometer~Sec.III!,
theprecisestructureofthedenoisingalgorithm~Sec.IV!and xk(m)[xk2d2m, m50,1,...,N21. ~1!
results ~Sec. V! obtained with simulated data and also with
For now, the delay d>1 also referred to as prediction
real data taken from the Caltech 40 m prototype interferom-
depth is fixed arbitrarily ~an explanation on how it may be
eter @11#.
chosen is given later in Sec. III!. The prediction is obtained
This work here is preliminary; its goal is to explore how
by linearly combining these data samples weighted by the N
effectivelyadaptivefilteringtechniquesperformontheprob-
correspondingcoefficientsw(m), formingthetap-weight~col-
lem we address. It is a first step towards a more complete
umn! vector w[(w(m),m50,1,...,N21)t. Therefore, the
statistical evaluation of the algorithm.
prediction y of x reads
k k
y [wtx . ~2!
II.METHODS k k
A.Fromhypothesistomethod The predictor is optimal in the mean square sense when
the variance of the prediction error e 5x 2y is minimum.
WeassumethatthenoiseconsistsofbroadbandGaussian k k k
Therefore, the problem is to find the set of weight coeffi-
noise plus large amplitude oscillating interference signals.
cients which minimizes
The model does not include any a priori knowledge of the
signal such as its exact frequency or shape of the envelope. Jk~w![E@ek2#5E@~xk2wtxk!2#, ~3!
Theonlyassumptionwemakeisthatitsautocorrelationover
a small time-lag d — the time-lag d chosen greater than the where E@(cid:149)# denotes the expectation value operator.
decorrelation time scale of the broadband noise — is appre- This leads to the minimization of the following quadratic
ciable, while for the broadband noise it is essentially zero. form
This difference can be used to advantage to discriminate be-
tween the narrow band interferences and the broadband Jk~w!5sk222wtpk1wtRkw, ~4!
noise.
The idea is to predict the current signal sample given a wheresk2[E@xk2#, pk[E@xkxk# andRk[E@xktxk#. Thereex-
collectionofpastsamplesofthedata,thesetwoobjectbeing ists only one solution wk*, obtained when the gradient of Jk
separated by the time-delay d. The prediction is effective, vanishes. This situation is realized when
only if the target sample shares enough information with
~i.e., is sufficiently correlated to! the previous samples. In Rkwk*5pk. ~5!
otherwords,theonlypredictablepartofthesignalistheone
When the signal is stationary, R 5R and p 5p are con-
whose correlation length is sufficiently large ~i.e., long-term k k
stant ~independent of k). In this case, R defines the autocor-
sinusoids or ringdowns!. Conversely, broadband noise can-
relation matrix of the signal x and the solution of Eq. ~5! is
not be predicted, as it is not possible to guess the next value k
referred to as the Wiener filter.
in this way. It is this crucial underlying idea we use to dis-
criminate between the two noise signals.
D.LinearpredictionandLMSmethod
B.Notation Equation ~5! requires the computationally expensive in-
version of the matrix R . An alternative and more efficient
As a general rule, we will denote scalar quantities by k
solutionforfindingtheminimumofJ (w) inEq.~4!consists
plainitalics,e.g.,x ; vectorsbyboldfaceletters,e.g.x ; and k
1 1 in starting from an arbitrary initial value w , and iterating
matricesbyboldfacecapitals,e.g.,X . Wewillrepresentthe k,0
1 the tap-weight vector along the steepest descent direction,
components of vectors and matrices with superscripts within
brorawckaentds,net.hg.coXlu1(mmnn) doefstihgenamteastritxheXe.leSmuebnstcrliopctast~endoibnramckth- wk,n115wk,n2m(cid:147)wJk~wk,n! ~6!
1
ets!denotethetimeindex;e.g.w denotesthevectorwatthe ~the new index n50,1,... counts the number of iteration
k
timeindexk.Thesymbol[ willbeusedinthefollowingto being implemented! given by the gradient
042004-2
ADAPTIVEFILTERINGTECHNIQUESFOR... PHYSICALREVIEWD63042004
J ~w!52~R w2p !. ~7!
(cid:147)w k k k
For a sufficiently small gain m, the weight vectors will
eventually converge to the optimal predictor filter w*. This
k
procedurerequiresthesecondorderstatistics~namelyR and
k
p ) ofthesignal.Inourcase,thisinformationisnotavailable
k
and one has therefore to estimate these quantities. Instead of
estimating R and p directly and combining them with Eq.
k k
~7!, a more efficient solution is to estimate the gradient. FIG. 1. The figure illustrates the principle of the underlying
From the derivation of Eq. ~3!, one can rewrite the gradient methodonwhichthealgorithmweproposeisbased.Thealgorithm
is designed to discriminate the nonstationary and non-Gaussian
as
noise features from the broadband background noise in interfero-
(cid:147)wJk~w!522E@ekxk#. ~8! amdeatpritcivgerlaivnietaetniohnaanlcwemavenetdaantda.itTshoibsjemcetitvheodisitsorceofemrrpeadretotheassiLgMnaSl
anditslinearprediction,thepredictorcoefficientsbeingadjustedby
A simple and natural way to obtain an estimator of this
afeedbackloopcontrolledbythepredictionerror.
quantity is to omit the expectation operator:
or the ringdown signals! and we send the rest of the signal
ˆJ 522e x . ~9!
(cid:147)w k k k for further analysis for detection.
Because the noise perturbs this estimate, the algorithm
mayiterateinadirectionwhichdoesnotliealongtheactual E.PropertiesoftheLMSmethod
~noisefree!directionofsteepestdescent,thuspreventingthe
Themethodwedescribedaboveisreferredtoasadaptive
filter from converging to the Wiener filter. For this purpose,
line enhancer (ALE). It is a special case of the LMS algo-
we stabilize the estimation above by setting the algorithm
rithm. Both, ALE and LMS algorithms have been first intro-
iteration index n equal to signal time index k in the Eq. ~6!. duced by Widrow and Hoff @12# in the 1960s.
Thefinalevolutionequationforthetap-weightvectorfinally The acronym LMS ~least mean square! designates a gen-
reads
eral scheme to design signal processing methods where a
minimization ~in a statistical sense! of a definite positive
wk115wk12mekxk. ~10! quadratic cost function ~usually related to some mean qua-
dratic error! is needed. Its central idea is the use of the esti-
This amounts to ~i! doing just one iteration at given time mate of the gradient of this function given in Eq. ~9!. The
k @i.e., applying Eq. ~6! once with n50] and ~ii! setting the
LMS technique has been extensively used for the last 30
result as the starting point for the next time index ~i.e., set-
years in communications problems such as echo cancella-
ting wk11,05wk,1). tion,channelequalization,antennaprocessing,etc.Themain
At a fixed time k, the weight vector evolves along the
advantages to be gained by applying the LMS technique are
crude estimate of the steepest descent direction. But on a ~i! adaptivity, ~ii! robustness and ~iii! simplicity.
longer duration, the direction followed by the tap-weight
In this context, the term ‘‘adaptivity’’ has two different
vector is governed by the sum of the successive gradient
meanings. First, it means that the LMS technique will auto-
estimates obtained with different noise samples. In other
matically adjust its parameters so as to reach the best setup
words,wehavereplacedanensembleaverageinEq.~8!bya
for a problem which has not been initially precisely defined.
time average. It also implies that we have implicitly called
Second,itisalsoabletofollowchangesinthecharacteristics
for further assumptions on the signal x : first its local sta-
k ofthedatabeingprocessedintheeventthattheyoccur.The
tionarity ~more precisely, the second order statistics are sup-
latter property also shows that the method is robust. In fact,
posed to be constant during the convergence time of the
this method has been proved to be robust according to spe-
algorithm! and second, its ergodicity. cific statistical criterion such as the minimax criterion @13#.
Obviously, the early-time estimates ~i.e., before the con-
The ALE is an adaptive prediction algorithm using the
vergenceiscompleted!maybefarfromoptimal,butaccord-
LMS technique. We have seen that the signal is predicted
ing to the previous argument, as the data accumulate, the
from a reference signal which is the signal itself. In some
estimation is expected to converge to the Wiener filter. Pre-
other applications, although the same principles are applied,
cise proofs of the statistical convergence of the algorithm in
the reference signal can be another signal, e.g., echo cancel-
Eq. ~10! can be found in @15#.
lation or denoising. In such cases, the quantity of interest
Summarizing, the method we propose consists in linearly
mightnotbethepredictionoutputbutthelinearfilterusedto
filtering the data to extract the part of the signal with a long
compute it, e.g., deconvolution.
correlationtime.AsillustratedwiththeblockdiagraminFig.
1,thefiniteimpulseresponsefilter~givenbyw ) ismodified
k III.ADAPTINGALEFILTERTOCANCELINGNOISE
at each iteration according to the relation ~10! with the final
INGWDATA
goal to minimize the mean square error. Once the filter has
converged ~i.e., w is stable in time!, we reject the predicted In this section we essentially describe a model for under-
k
part of the signal ~corresponding to the long-term sinusoidal standing the behavior of the ALE algorithm. The model we
042004-3
E.CHASSANDE-MOTTINANDS.V.DHURANDHAR PHYSICALREVIEWD 63042004
assume consists of a high amplitude narrowband signal su-
2
perposed on broadband noise. For simplicity, we assume the w*(m)5 cos~m1d!df, m50,1,...,N21,
broadband noise to be white and Gaussian and the narrow- N14s2
band signals are sinusoids of constant envelope. The results ~15!
we obtain hold for more realistic signals when the evolution
where we have chosen the length of the filter to be half-
of their amplitude and/or instantaneous frequency occurs
integral number of cycles for reasons of simplicity, i.e.,
adiabatically, i.e., the change is small over the period of the
Ndf5lp, where l is an integer.
sinusoid.
In other words, the optimum linear predictor is nothing
The assumption of white noise is not too restrictive be-
butacopyoftheexpectedsignalitself.ThefilterinEq.~15!
causethisisequivalenttochoosingthenoisecorrelationtime
is also referred to as the matched filter. In our situation, in
to be zero and therefore we are free to choose the prediction
depth~i.e.,thetimedelaybetweenthecurrentpredicteddata practice, N@s2 and the term 4s2 can be omitted from the
sample and the reference signal to the LMS filter! to be ar- amplitude of w*.
For the reasons detailed before, we propose to use the
bitrarilysmall.Inarealsituation,wemustfixthedelaytobe
ALEalgorithminordertofindagoodapproximationofw*.
greater than the correlation time of the broadband noise. We
Startingfromanarbitraryinitialtap-weightvector,weiterate
first analyze the case of the sinusoid because it is easier to
the weights w according to Eq. ~10! to converge to w*.
investigate and provides invaluable insights into the work- k
Oncethefilteris‘‘close’’enoughtotheoptimalsolution~the
ings of the LMS algorithm.
word ‘‘close’’ will be defined later in the text!, we then say
It may be remarked that the denoising of sinusoids in
that the filter has locked on to the signal.
whitenoisehasbeentreatedintheliteraturewithgreatdetail
~see @13,14# for a review!. We give here only pertinent re-
sults ~with a short proof! for introducing the structure of the B.Approachtolocking
algorithm, which we present later in the text.
a. Continuous time approximation of the locking trajec-
tory. We may analyze the approach to locking by deriving a
A.Optimalfilter differenceequationfortheaveragedevolutionoftheweights
and then investigating this equation. It is impossible to ob-
We consider the data to be of the following form
tain the average evolution of the weights by using the stan-
x [cos~2pf t 1F!1n , ~11! dard definition of the expectation operator E because of the
k 0 k k
nonlinearity and the recursive scheme involved in evolving
where t [kd, d[1/f being the sampling interval and F is the weights. We therefore adopt the time-average over suc-
k s
arandomphase~attheorigin!withuniformprobabilityden- cessive data points as the operational definition of E.
sity function between 2p and p. The sinusoid has fre- Shifting the origin to w* by defining v [w 2w*, we
k k
quency f and the units are so chosen that it is of unit am- may write the LMS evolution equation ~10! in the following
0
plitude. The additive white noise n with variance s2 form @15#:
k
satisfies the relation,
vk112vk522m~xkxkt!vk12mek*xk, ~16!
E@nknm#5s2dkm, ~12!
where e*[x 2w*tx is the prediction error produced when
k k k
where d is the Kronecker delta. using the optimal filter.
km
As mentioned in Sec. IIC, the problem is to optimally During the locking phase, the filter is far apart from the
predict x from a reference signal x which is just the de- optimallocation~i.e.,v hasalargemodulus!.Thehomoge-
k k k
layeddatabytheamountdd~orequivalently,anumberofd neous term dominates the forcing term in the difference
samples!. More precisely, the reference is a column vector equation ~16! which then can be approximated by
whose components are equal to xk(m)[xk2d2n, n
50,1,...,N21. vk112vk522m~xkxkt!vk. ~17!
The best linear prediction filter can be obtained by evalu-
atingtheautocorrelationmatrixRandthevectorpinEq.~4!. Inthesituationwherethestepgainparametermischosen
In case of signal ~11!, their components are given by to be very small so that the weight coefficients are almost
constantoveragiventimeinterval,therecursivityeventually
R(mn)51/2cos~m2n!df1s2d , ~13! acts as an averaging operation on both sides of the equation
mn
above.Thisleadstothedifferenceequationwhichweuseto
p(m)51/2cos~m1d!df, ~14! describe the tap-weight trajectory in the space of weight co-
efficients we denote W:
where (m,n)50,1,...,N21 and df52pf d. Note that
0
we have dropped the index k because the autocorrelation R vk112vk522mRvk. ~18!
doesnotdependuponk,sincewearedealingwithastation-
ary signal. Let Q be the transformation which diagonalizes R. The
From the above expressions of R and p and solving Eq. above difference equation is best analyzed by changing the
~5!, we obtain the optimum prediction ~Wiener! filter frame in W to the principal axis
042004-4
ADAPTIVEFILTERINGTECHNIQUESFOR... PHYSICALREVIEWD63042004
˜vk112˜vk522mR˜˜vk, ~19! oftThheiseiegxepnrveaslsuieoneqhuaastitoonbSevc5omlpva,reledawdiinthgttohetwseocosonlduttieornms
where R˜[QRQ215diag(l(0), ...,l(N21)) and ˜v [Qv . for a, namely, a5g/2 and a5(g2p)/2. These yield the
Equation ~19! gives decoupled difference equationsk for thke eigenvectors v6 and the corresponding eigenvalues l6:
components ˜vk(m),m50,...,N21 of ˜vk which can be v15nexp~2ig/2!1¯nexp~ig/2!, ~29!
solved given the initial weight vector˜v(m)
0
v25inexp~2ig/2!2i¯nexp~ig/2!,
˜v(m)5˜v(m)~122ml(m)!k. ~20!
k 0
l65~N6b!/2. ~30!
b. Eigenvalues and eigenvectors of R. We need to com-
pute the eigenvalues l(m) of R. This can be conveniently If we choose N large enough and Ndf5mp, where m is
implementedbysplittingRintoanoisepart,whichisjusts2 an integer then the analysis becomes simpler. This amounts
times the identity, plus the signal part which we denote by tochoosingthelengthofthefiltertohavehalf-integralnum-
S/2, and thus ber of cycles: we have b50 and l65N/2. ~Geometrically,
this means that the eigenvalue problem is degenerate with
R5s2I1S/2, ~21!
respect to the two signal eigenvectors: there is a two dimen-
sional eigenspace belonging to the eigenvalue N/2. The
where S(mn)[cos(m2n)df. It is easily verified that the
weights thus evolve non-preferentially with respect to the
eigenvectors of R and S are identical and the eigenvalues of
signal eigendirections.! Since typical cases imply generally
RareobtainedfromthoseofSbyfirsthalvingthemandthen
N@b, we will assume this simplification in the rest of the
addings2 totheresult.Itremains,therefore,tocomputethe
paper.
eigenvalues and eigenvectors of S. We do this by observing
In this situation, the spectrum of R
that we can write S as follows:
sp~R!5$l(0)5l(1)5N/41s2 and
S5~nn†1¯n¯n†!/2, ~22!
l(m)5s2,m52,...,N21%, ~31!
where n[@1,exp(idf),exp(2idf),...,exp (N21)idf#t.
(cid:132) (cid:133)
Since the matrix S is real and is essentially made out of consist of two sets of eigenvalues : the first two correspond
two external products of nand ¯n, its rank equals 2 (S has to directions in the signal space associated to ‘‘signal
N22 degenerateeigendirectionsinWwitheigenvaluezero! 1noise’’ ~or ‘‘signal,’’ for short! whereas the remaining N
andhastwononzerorealeigenvalues.Letv beaneigenvec- 22 characterize ‘‘noise’’ directions.
torassociatedtooneofthenontrivialeigenvalues.According According to Eq. ~20!, the weight vector will converge
tothestructureofS, thevectorv canbewrittenwithoutloss more rapidly in directions associated with the largest eigen-
of generality as the following linear combination: values, which are the signal eigenvalues. The other noise
eigenvaluesareunimportantinthisconsideration.Theeigen-
v5nexp~2ia!1¯nexp~ia!, ~23! vectors pertaining to the signal provide preferred directions
in W: it is along these directions that the slope of the per-
where the coefficients have been chosen arbitrarily to have formance surface is steep and hence promotes faster conver-
unit modulus. gence.
Using the two scalar products n†n5N and
C.Steadystateevaluation
ntn511exp~2idf!1 1exp 2~N21!idf ~24!
(cid:149)(cid:149)(cid:149) (cid:132) (cid:133)
If the step gain factor is sufficiently small, the tap-weight
[bexpig, ~25! coefficients eventually converge and stabilize in a neighbor-
hood of the optimal value. At this stage, the assumptions
where the geometric series can be summed up and modulus made in obtaining the approximated evolution equation ~18!
and phase ascertained do not hold anymore. In contrast with the case of ‘‘the ap-
proach to locking,’’ the right hand side of the difference
sin~Ndf! equation ~16! is now dominated by the forcing term:
b5 , ~26!
sindf
vk112vk52mek*xk. ~32!
g5~N21!df, ~27!
Roughly speaking, the trajectory of the vector w during
k
the steady state can be viewed as a random walk centered
weobtaintheeffectofthematrixSonthevectorv, givenby
aroundw* lyingwithinaregionofWspacewhoseextentis
Sv5nexp~2ia!@N1bexp~2ig12ia!#/21c.c., determined by two factors, namely, mand the intrinsic ge-
~28! ometry of W in the vicinity of w*.
The misalignment between the actual ALE filter w and
k
where c.c. denotes complex conjugate. the optimal one w* creates an additional error in the output.
042004-5
E.CHASSANDE-MOTTINANDS.V.DHURANDHAR PHYSICALREVIEWD 63042004
In fact, a direct calculation from Eq. ~4! shows that the total and vectors!, the second, j(s) with the signal. Because the
k
mean square error may decomposed as @12# signal eigenvalues are much larger than those of the noise,
thesuminEq.~34!isessentiallydominatedinthebeginning
Jk~wk!5jmin1jk, ~33! ~for small k) by j(s). These two errors decrease during the
k
where ~i! j [J (w*) is the minimum mean square error locking phase until reaching a steady state value. The lock-
min k ing time ~i.e., time at which the steady state is reached! is
arising from the fraction of the input noise which still re-
definedtobethat,whenj(s) isoftheorderofthetotalEMSE
mainsintheoutput,assumingthattheALEfilterhasreached k
exact optimality and ~ii! j[vtRv is the excess mean expected in the steady state.
square error ~EMSE! due tokthekmiskalignment between the From Eqs. ~20!, ~30! and ~34! we obtain
ALE filter and the Wiener filter.
j(s)[l(0)~˜v(0)!21l(1)~˜v(1)!2 ~38!
One can verify that the EMSE vanishes when reaching k k k
optimality i.e., when v 50. In other words, this term quan- S D
tifies the nonoptimalitykof the current filter in use. We can N mN 2k
’ ~˜v(0)!21~˜v(1)!2 12 ,
imaginej asthesquareofanaturaldistanceinWandRas 4 (cid:132) 0 0 (cid:133) 2
k
a intrinsic metric over W. ~39!
A good approximation of j can be found for large
min
number of weight coefficients for the specific case of sinu- where we have assumed N@b.
soidal signals with high SNRs. Using Eqs. ~4! and ~5!, we We set the starting point w0 in W to be 0. It corresponds
maywritejmin5E@xk2#2ptw*. WhenN!‘, adirectcalcu- to the initial value ˜v0 in the eigenspace which is given by
lationshowsthatthesecondtermptw* tendstotheenergyof ˜v 52Qw*. The first two coordinates of˜v can be directly
0 0
the sinusoid, which means that the remaining energy is that obtained since the first two row vectors of Q are just the
of the noise: jmin’s2. normalized signal eigenvectors of the R matrix, leading to,
We complete the characterization of the mean square er- for half integral wavelength filters and N@s2,
ror ~33! with the evaluation of the average value of the
EMSE, which we denote by j(st). First, noticing that the ˜v(0)5A2/Ncos~ddf1g/2!, ~40!
EMSE is invariant under the principal axis transformation 0
A
N21 ˜v(1)5 2/Nsin~ddf1g/2!. ~41!
( 0
j5˜vtR˜˜v 5 l(m)~˜v(m)!2, ~34!
k k k k
m50 These considerations yield
S D
and secondly, using the approximation E@˜vk˜vkt#’mjminI 1 mN 2k
proposedin@12#toobtainthetypicalvaluefor(˜v(m))2 yields j(s)’ 12 . ~42!
k k 2 2
N21
(
j(st)’m l(m)j . ~35! This error must now be compared with the averaged
min
m50 EMSEinEq.~36!inordertofindthetimet atwhichj(s)
lock
and j(st) are equal
Since the signal is of much larger amplitude than the
broadband noise, the trace of R˜ is essentially due to the sig- ln~mNs2!
nal eigenvalues @see Eq. ~31!#. Combining with the expres- t ’d . ~43!
lock 2ln~12mN/2!
sion of j above, this leads to
min
j(st)’mNs2/2. ~36! It is important to mention that, when the product mN/2
tendsto1, theconvergencetimedivergestoinfinitymeaning
A better estimate of j(st) can be obtained starting with that the weights do not converge toward w* anymore. In
more realistic hypotheses and using more sophisticated ap- order to ensure the stability of the algorithm, the parameters
proximations @14# : need to satisfy the stability condition 0,mN/2,1. How-
ever, we have observed in our simulations that when 1/2
N(21 ml(m) <mN/2,1, theconvergenceissloweddown.Thisisdueto
j(st)’m50 12ml(m)jmin. ~37! thepresenceofoscillatorytermsinthegradientwhichdonot
average to zero anymore. In practice, it is advisable to
choose the parameters so that mN/2,1/2.
Inthelimitofsmallstepsize,thisapproximationtendsto
the simpler one in Eq. ~36!. ForasinusoidofamplitudeAinsteadofunityaswehave
considered before, the condition for stability can be simply
obtained by replacing the parameter mN/2 by r[mNA2/2
D.Convergencetime
leading to 0,r,1.
In the expression of the EMSE in Eq. ~34!, we separate We illustrate in Figs. 2 and 3 with an example the results
the sum into two parts: the first, j(n) associated with the of this section pertaining to the approach to locking and
k
noise ~i.e., consisting of terms involving noise eigenvalues steady state analysis.
042004-6
ADAPTIVEFILTERINGTECHNIQUESFOR... PHYSICALREVIEWD63042004
FIG.2. ApplyingtheALEtoasinusoidalsignal:approachtolockingandsteadystate.ThefigureillustrateshowtheALEperformson
a test signal composed of a sinusoid (f 550 Hz sampled at f 51 kHz! corrupted by white noise @SNR5A2/(2s2)550(17 dB)#. We
0 s
initialize the N540 tap-weight coefficients to 0, set the prediction depth d55 ~ms! and the step gain parameter m50.003. ~a! The filter
output signal y ~solid line! converges rapidly towards the actual noise free sinusoidal signal ~dashed line!. ~b! This is confirmed by
k
observingthatthe‘‘signal’’EMSEjs definedinEq.~38!whichdecreaseswithtimeuntilitreachesitssteadystatevalue.Forcomparison,
k
thehorizontaldashedlineindicatesthetheoreticalmeanvalueofthetotalEMSEjst inthesteadystate@seeEq.~36!#.Wecanverifythat
thetheoreticalvalueobtainedinEq.~43!fortheconvergencetimet correspondseffectivelytothetimeinstantatwhichjs andjst are
lock k
ofthesameorder.Finally,thetwocontourplotsofthebottomlinedisplaythetrajectoryfollowedbytheadaptivefiltercoefficientsinthe
eigenweightspace:in~c!,theaxisarethefirsttwoeigenvectorsofRnamely˜v(0) and˜v(1) ~i.e.,the‘‘signal’’eigenvectors!whereasin~d!
k k
the diagram plane is given by ˜v(1) and ˜v(2) ~i.e., a ‘‘signal’’ direction vs a ‘‘noise’’ direction!. As proved in Sec. IIIB, the weight
k k
coefficientsconvergemorerapidlyalongthetwodirections˜v(0) and˜v(1) givenbyeigenvectorsassociatedwiththelargesteigenvalues.
k k
IV.THEALEINPRACTICE frequency remains large. Therefore, the model introduced
previously,namelythatoflargeamplitudesinusoidalsignals
Intheprevioussectionswehavecharacterizedthebehav-
embeddedinbroadbandnoise,isareasonableapproximation
ioroftheALEincasesofinterest.Wewillnowelaborateon
within the relevant small bandwidth of frequencies.
how this algorithm can be adapted to the interferometric
Forthisreason,wedividethefrequencyaxisinpdisjoint
data.
frequency subbands of the same size. The p signals lying in
In the scheme we present here, we first decompose the
each of the subbands are heterodyned and decimated to the
signal in p frequency subbands to which we apply the ALE
sampling frequency fband[f /p.
twice with different sets of parameters. In the first stage, the s s
Thetilinghastheadvantagethat,ifpissufficientlylarge,
parameters are tuned to best remove long-term sinusoidal
we can consider the interferometer background noise almost
components of the noise; whereas in the second stage, the
white within a subband, which implies that the noise has
target consists of shorter oscillatory transients.
vanishingcorrelationtime.Thepredictiondepthdwhichhas
to be larger than the correlation time, can be then simply
A.Subbanddecomposition fixedtoanyvaluegreaterthan1 sampleperiodineachofthe
subbands.
Interferences such as mains power and violin mode har-
monics are distributed over a large dynamic scale ~the first
B.Long-termsinusoidremoval
harmonics are of much larger amplitude than those of high
order!. But, since the interferometer noise curve also de- Certain parts of the spectrum may not contain any long-
creases at low frequencies, their relative amplitude as com- term periodic interferences. We apply a preliminary test to
pared to the background noise power spectrum at the same exclude subbands which may not require the first denoising
042004-7
E.CHASSANDE-MOTTINANDS.V.DHURANDHAR PHYSICALREVIEWD 63042004
FIG.3. ApplyingtheALEtoasinusoidalsignal:convergencetime.Thefigureshowsthecomparisonbetweentheconvergencetimet
lock
obtainedbysimulations~solidlinewith‘‘1’’associatedwithonenoiserealization!anditstheoreticalvalue~solidline!giveninEq.~43!.
ThetestsignalisasinusoidinwhiteGaussiannoise@seeEq.~11!withA51, f 550Hz,sampledat f 5200Hz#.Theconvergencetime
0 s
isshownasafunctionofs2 in~a!~wehavefixedtheremainingALEparameterstoN5200andmNA2/250.01) andin~b!asafunction
ofmNA2/2(N5200ands250.1). SimulationsgloballyconfirmtheresultsobtainedinEq.~43!exceptwhen1/2<mNA2/2,1. Thereason
forthisdiscrepancyisthatthedifferencebetweentheactualgradient~8!anditsestimates~9!canbeshowntobeanoscillatingtermwhich
doesnotaveragetozeroanylongerwhenthestepgainparameterapproachesthecriticalvalueforstability.Therefore,inpractice,wechoose
parameterssothatmNA2/2,1/2.
step. The test is crudely done by estimating the amplitude A The L2 norm iw*i25(2/N2)@N1bcos(g12ddf)# is ob-
2
of the sinusoid from the largest peak of the power spectrum tainedbysquaringandsummingtheEq.~15!fortheoptimal
~Welch estimate! and comparing it to the variance s2 of the filter. Since in typical cases b!N, this leads to simpler ex-
broadband noise ~also estimated from the power spectrum!. pression iw*i2’2/N.
If it is found that A.s, we decide that there exists a long- Consequent2ly, the number of tap-weight coefficients N
term sinusoidal signal of sufficient amplitude in the band has to be chosen so that
whichneedstoberemoved,otherwiseweproceeddirectlyto
the second step. N>2/h . ~46!
noise
We apply the ALE in each of the selected subbands
choosing parameters as follows. Step gain parameter m
Number of tap-weight coefficients N
Wefixthestepgainparameterbyimposingthecondition
The number of tap-weight coefficients is fixed by pre-
that on the average, the distance of the ALE filter from op-
scribing an upperbound 0,h ,1 to the ratio between
noise timality in the steady state is smaller than a given threshold.
the noise power corrupting the filtered output y of the opti-
k As we have seen in Sec. IIIC, this can be done naturally by
mal filter w* and the input noise power. Let the column imposing an upperbound 0,h ,1 on the ratio of the ex-
vector nk[(nk2d2m,m50,1,...,N21)t be a collection of cess square mean error to the ssiiggnal power E 5A2/2:
noise samples, then the above condition reads s
E@~w*tnk!2#<hnoiseE@~nk!2#, ~44! j(st)<hsigEs. ~47!
which,withthestationarityandwhitenessofthebackground Using the expression obtained in the steady state analysis
noise n , results in a bound on the optimal filter gain: in Eq. ~36! for the EMSE, this condition reduces to
k
w*tw*<h . ~45! m<h /~Ns2!. ~48!
noise sig
042004-8
ADAPTIVEFILTERINGTECHNIQUESFOR... PHYSICALREVIEWD63042004
Generally, this equation leads to small values of mwhich on each data segment. Accordingly, we have added a super-
prevent the convergence of the ALE filter from its initial vision test which decides whether or not the denoising algo-
state ~i.e., all tap-weight coefficients are fixed to 0) in a rithm should be applied to a given data segment. The test
reasonable time ~convergence faster than a tenth of second, consists of observing the Gaussianity of the filtered output
which is the duration of the chunk of data!. We solve this y 5wtx . Iftheinputsignalx isazero-meanwhiteGauss-
k k k k
problembyfirstapplyingtheALEonasequenceoftraining ian process of variance s2, then the output of the filter y
k
data, the step gain parameter being set at the beginning to a shares the same characteristics, except that the variance gets
largevalue~forfastconvergence!anddecreasedgraduallyto multiplied by the filter gain: vary 5iw i2s2. Furthermore,
the value given in Eq. ~48!. The filter obtained after the under this hypothesis, the envelopke Y k52iH(y) i2 @H(y)
completion of this training is close to the objective ~i.e., the denotes the ~complex valued! analyticksignal1 asksociated tok
Wiener filter!. We then start the longterm sinusoid removal the real signal y ] follows by definition a chi-square distri-
k
using this prepared filter. bution with 1 degree of freedom.
We remark here that although mis small, it is nonzero This implies that, up to an arbitrary probability P , the
0
thus giving the ALE filter some flexibility of adapting to envelope Y does not exceed the threshold given below:
changes~nonstationarities!inthesignalsuchasslowdriftsin k
frequency and amplitude modulation. This property, how-
Y ,k~P !2fbandiw i2s2, ~51!
ever, needs to be investigated more in detail. k 0 s k 2
where k( ) is the inverse function of the ~unit variance! x2
C.Ringdownremoval (cid:149)
cumulative distribution function ~cdf!.
The aim of the second step of the algorithm is to remove IfEq.~51!issatisfied,weconcludethatthefilteredoutput
oscillatory transients ~ringdowns! of large amplitude. These is essentially due to a Gaussian background noise and we
transients are either frequency bands excited from time to leave the input signal as it is. Otherwise, we conclude that
time ~caused by dysfunctions in the interferometer! or relics the filtered output carries a ringdown signal and decide to
from the previous step ~when the envelope of a long-term remove it from the input data.
sinusoid possesses fast variations to which the algorithm The functioning of the second step of the denoising algo-
cannotadaptorconvergetoduringthefirststepofremoval!. rithm could be interpreted as follows : it removes from the
ThecleaningprocedureconsistsinapplyingALEthesec- input data, regions in the time-frequency plane presumably
ondtimetoeachofthesubbandsbutnow,theparametersare associatedwithtransients,whosesupportisdefinedalongthe
so adjusted that, ~i! they select features with a larger band- frequency axis by the ALE filter, and along the time axis by
width than in the previous step, and ~ii! converge rapidly the supervision criterion ~51!.
onto an oscillatory noisy signal that may appear. After completing these two steps, we recombine the sig-
Number of tap-weight coefficients N nal in all the subbands together to retrieve a single strain
The impulse response duration and frequency selectivity signal.
~i.e., the filter bandwidth Df) of the transfer function are
dual in character. This follows from the uncertainty relation.
The rough approximate relation between these quantities is V.NUMERICALRESULTS
given by
A.Simulateddata:testoftheringdownremoval
N5f /~pDf!, ~49! In this section, the goal is to test how effectively the sec-
s
ond stage of the denoising algorithm ~i.e., the ringdown re-
where fs is the sampling frequency. We choose the number moval!describedinSec.IVoperatesonasimplesignal.The
of tap-weight coefficients N by imposing a minimum band- test signal is composed of three ringdown signals ~of fixed
width Df to the filter and using the above equation. amplitude and frequency! occurring successively in the data
min
Step gain parameter m stream and embedded in a additive Gaussian white noise.
Assuming that the ringdown can be locally approximated This model may be used to represent ringdown disturbances
by a sinusoid, we choose the step gain parameter by impos- originating from the same underlying physical mechanism.
ing a convergence time of the order of a typical transient Each of these ringdown signals is a sinusoidal waveform,
duration ~i.e., t ’Nd). More concretely, setting r similar to Eq. ~11! ~with A51, f 550 Hz and sampling
lock 0
[mNA2/2 in the unnormalized form of Eq. ~43! ~i.e., for frequency f 5200Hz!,whosesupportislimitedintimebya
s
arbitrary ringdown amplitude A), we solve for Gaussian envelope:
ln~2s2r! r [Aexp 2p~t 2t !2/T2 cos~2pf t 1F!, ~52!
5N. ~50! k (cid:132) k c (cid:133) 0 k
2ln~12r!
UsingthecrudeestimateA2/2’ixki22/N fortheringdown 1The analytic signal y 5H(x) associated to the signal x is es-
amplitude, the step gain parameter is finally obtained as m n n n
sentiallyobtainedbycancellingitsnegativefrequencies;morepre-
5r/ixki22. cisely, Y(f)[2U(f)X(f), with U(f)51 when fP@0,1/2# and 0
Since the ringdown signals are of short duration and can when fP]21/2,0@ andwhereX(f) @andY(f)] denotestheFourier
occurwithlargetimegaps,theALEdoesnotneedtooperate transformofthecorrespondingsignalX(f)[(N x e22pinf.
n50 n
042004-9
E.CHASSANDE-MOTTINANDS.V.DHURANDHAR PHYSICALREVIEWD 63042004
FIG.4. ApplyingtheALEtooscillatingtransients:testingtheringdownremovalalgorithm.Thefiguredepictstheresultsofthetransient
removalalgorithmpresentedinSec.IVCtoasignalin~a!composedofthreesuccessiveoscillatingbursts~seetextfordetails!embedded
in Gaussian white noise @SNR5A2/(2s2)58(9 dB)#. In addition to the ALE, we measure the deviations of the filtered output from
Gaussianity: when its envelope ~c! exceeds some threshold @dashed line, see Eq. ~51!#, we decide that the filtered output is not normally
distributedand,thereforecontainsatransientwhichhastoberemoved~thisisindicatedbydotsatthetopofthegraph!.Thefinalneteffect
oftheoperationisthatoftimevariantfilteringoftheinputdata.Thecorresponding~timevarying!transferfunctionisrepresentedin~d!:the
regions indicated with dark colors are parts of the time-frequency plane where the data are selected and removed from the input ~i.e., it
correspondstothetime-frequency‘‘band’’passofthefilter!.Comparingthespectrograms@21#@seeEq.~53!foradefinition#in~b!and~f!
respectivelyoftheinputandoutput~e!signals,weobservethatthethreetransientsareprogressivelyremovedfromtheinput~darkregions
representlargevaluesofthetime-frequencyenergydensity!.
where three different reference times t are given and the Notice that real time and frequency coordinates can be
c
equivalent time duration is T5200 ms ~giving a frequency retrieved through the relations t5n/f and f5mf .
s s
bandwidth of Df’1/T55 Hz and Q[f T’10 cycles!.
0
Figure 4 describes the application of the denoising algo-
B.ResultsonCaltech40mprototypedata
rithm configured with d55 sampling periods ~equal to 25
ms! Df 53 Hz, and P 50.01. It can be seen that the Here we have applied the algorithm to the Caltech 40 m
algorithmminoperates better o0n the transient encountered later prototype data taken in October 1994 @11#. This data was
in the data train than its predecessor. The explanation is that recorded with a sampling frequency of fs59.86 kHz. We
a transient duration is too short for the filter to reach the have used the calibrated strain signal @17# ~relative arm
steady state but, when it encounters the next transient, the length measurement! for applying our algorithm.
filter benefits from the distance to w* previously covered, Wetilethecompletespectruminto p532frequencysub-
thus improving the convergence towards optimality. bandsofapproximately154Hzeach.Eachsubbandencoun-
This can be verified with a time-frequency representation ters typically one or two long-term sinusoidal interferences.
@16# of the output signal such as Fig. 4, where we have Wehavechosenthepredictiondepthtobed55 sampling
chosenthespectrogramSh@n,m#[uFh@n,m#u2 definedasthe periods, which corresponds to a delay of pd/fs’16 ms in
x x
squared modulus of the short-time Fourier transform : real time. The correlation time of the broadband noise is
effectively smaller in each subband except at the extremities
of the spectrum where the steep slope of the spectrum does
(
Fxh@n,m#[ xnhk2ne22pinm, ~53! not allow us to assume the background noise to be locally
k white.Itonlyaffectsthefirstandlastsubbandswhicharenot
too important for detection purposes.
where nP@1,2,...,N#, mP]21/2...1/2] and h is an ar- In the first stage, we have chosen h 50.01 @giving
k noise
bitrary window ~a Gaussian window here!. N5200accordingtoEq.~46!#andh 50.01. Inthesecond
sig
042004-10