9512.net
ÌðÃÎÎÄ¿â
µ±Ç°Î»ÖÃ£ºÊ×Ò³ >> >>

# Fast algorithm for directional time-scale analysis using wavelets

Fast algorithm for directional time-scale analysis using wavelets
Rob A. Zuidwijka and Paul M. de Zeeuwa a Center for Mathematics & Computer Science (CWI) P.O. Box 94079, 1090 GB Amsterdam, The Netherlands
Fast algorithms performing time-scale analysis of multivariate functions are discussed. The algorithms employ univariate wavelets and involve a directional parameter, namely the angle of rotation. Both the rotation steps and the wavelet analysis/synthesis steps in the algorithms require a number of computations proportional to the number of data involved. The rotation and wavelet techniques are used for the segregation of wanted and unwanted components in a seismic signal. As an illustration, the rotation and wavelet methods are applied to a synthetic shot record. Keywords: (bi-)orthogonal wavelets, fast wavelet transform, interpolation, rotation, seismic data processing The topic of this paper is the wavelet analysis of functions in several variables using univariate wavelets. In particular, we study fast algorithms to perform such analyses. Observe that the well-known pyramid algorithm for multivariate functions12 , which uses tensor products of univariate wavelets, is such a method. An important aspect of the methods in this paper is the introduction of a directional parameter, namely the angle of rotation. Besides a combination of the usual pyramid scheme and rotation, we also study a method based on the wavelet X-ray transform. The wavelet X-ray transform10,13,21 given by

ABSTRACT

1. INTRODUCTION

t ? b dt a computes a one-dimensional wavelet transform of the restriction of f 2 L2(Rn ) to the line fx + t j t 2 Rg, where 2 Rn is a unit vector and x 2 ? = fy 2 Rn j y ? g. The translation, dilation parameter pair (b; a) is taken from the open upper half plane f(x; y) 2 R2 j y > 0g, just as for the usual 1-D wavelet transform. The wavelet satis es additional conditions so that the function f can be reconstructed from its wavelet X-ray coe cients P f ( ; x; b; a). P f ( ; x; b; a) =
1 f (x + t ) pa R

Z

Properties of the wavelet X-ray transform were discussed in10,13 . The transform was discretized by means of Fourier methods in18 and used to detect linear events in SAR images there. The transform has received further attention in21,22 , where an alternative discretization was proposed. A fast implementation of this discretization simply comes down to the computation of the 1-D fast wavelet transform along parallel gridlines in a rotated grid. Recall that the 1-D fast wavelet transform computes wavelet coe cients of a function starting with function values given at equidistant points. In this paper, we shall not recite the 1-D fast wavelet transform and refer to4,11,12 for details. We mention only the immediate fact that computing fast wavelet transforms along parallel gridlines requires a number of computations proportional to the number of gridpoints. In practice, a function f is usually given on the standard Euclidean grid by means of values at the grid points. These values can be interpreted as point evaluations or as certain averages. The latter interpretation ts in the context of biorthogonal Riesz systems. Some remarks on biorthogonal Riesz systems are given in Section 2. The application of the fast wavelet X-ray transform or a rotated version of the 2-D wavelet transform requires values at the grid points of a rotated (and possibly dilated) version of the Euclidean grid. As a consequence, the following question comes up: Given values of f at the standard Euclidean grid, what are the best values of f at the grid points of the rotated and dilated grid? Clearly, this is an interpolation issue, and it will be dealt with in this paper. Optimal values of f -in least squares sense- at the rotated and dilated grid are described in Sections 3 and
Other author information: R.A.Zuidwijk@cwi.nl, Paul.de.Zeeuw@cwi.nl, dbs.cwi.nl/cwwwi/owa/cwwwi.print
projects?ID=64

4. Moreover, in Section 4, it will be shown that the number of computations required for the rotation procedure is proportional to the number of data. Finally, we consider an application from the eld of seismic data processing. For obvious reasons geophysicist are after geometric information on the strati cation of impedance beneath the earth's surface. Acquisition of data is done by exciting a signal at the surface. Waves are re ected at interfaces (i.e. where the impedances changes rapidly). Seismometers at the surface record the groundmotion as the re ecting waves arrive. As the de ections are measured in time and for an array of seismometers, we can represent the recorded signals as a two-dimensional gridfunction (i.e. a rectangular uniform grid where each gridpoint has been assigned a real number). The information of interest is then constituted by hyperbolic shaped events. Unfortunately, these events are overshadowed by groundroll and surface waves directly stemming from the excitative source. Indeed, these waves often have dominant amplitude and completely blur the picture. It is a major challenge to lter out the surface waves from such given seismic data. We propose a synthetic and strongly simpli ed testproblem that might serve as a benchmark. We sketch a numerical algorithm for the solution of this testproblem in Section 5.

We introduce Riesz systems in Hilbert space and focus on geometric issues related to biorthogonal pairs of Riesz systems. For general reading on Riesz systems, we refer to17 , although many wavelet textbooks such as1,4,5,9 contain material on the subject relevant to wavelets as well. Recall that a Riesz system in a Hilbert space H is a sequence of vectors (xk )1 with two positive constants 0 < A B such that k=1

2. BIORTHOGONAL RIESZ SYSTEMS 2.1. Biorthogonality and Hilbert space geometry

A

N X k=1

jak j

2

N X k=1

ak xk

2

B

N X k=1

jak j

2

H

for all nite sequences a1 ; : : : ; aN . If the constants A; B are chosen optimally for the Riesz system under consideration, then A; B are called the Riesz bounds of the Riesz system. Let (xk )1 be a Riesz system in the Hilbert space H k=1 e k=1 with closed linear span V = span(xk )1 . We rst remark that there exists a unique sequence of vectors (xk )1 in k=1 V such that the biorthogonality condition hxk ; xl iH = k;l e (1) 1 is a Riesz system and we get span(xk )1 = V . Finally, if the Riesz bounds of is satis ed. It turns out that (xk )k=1 e e k=1 e e e k=1 e k=1 (xk )1 are given by A; B , then the Riesz bounds of (xk )1 are given by A = B ?1 ; B = A?1 . A Riesz system (xk )1 k=1 which satis es (1) will be called a dual Riesz system with respect to (xk )1 . Observe that a Riesz system is its own k=1 dual if and only if it is orthonormal. This particular situation corresponds to the case when A = B = 1. If we allow e the sequence (xk )1 to span a subspace V H which is di erent from the subspace V , then the biorthogonality e k=1 condition (1) does not imply that (xk )1 is a Riesz system. However, we can state the following results23 . e k=1 e Theorem 2.1. Let V; V H be subspaces in Hilbert space, and assume that V contains a Riesz basis (xk )1 with k=1 e e e Riesz bounds A; B . Then V contains a Riesz basis (xk )1 , with Riesz bounds A; B , biorthogonal to (xk )1 if and e k=1 k=1 e e only if V V ? = H , i.e., if and only if there exists a bounded projection P onto V along V ? . In that case, the Riesz 1 in V is uniquely determined and the Riesz bounds of the systems are subject to e basis (xk )k=1 e
2

e AB

2

;

2

e AB

2

;

where

. e e e Observe that in the case when V = V , the theorem reproduces the fact that A = B ?1 and B = A?1 . The bounds in the theorem are not sharp but can be attained in speci c cases, as can be seen even in nite dimensions23 . In the e papers14{16 , the constant as in Theorem 2.1 arises as the reciprocal cosine of the angle between V and V and is used there to obtain an important instance of the following result. This instance is discussed in the next subsection. e Observe that the case when = 1 corresponds to the case when V = V .

Px inf k y = inf kkxkk = 06=y2V kPyk k ; e 06=x2V

Px y = sup kkxkk = sup kPyk k : k 06=y 2V e 06=x2V

e e Theorem 2.2. Let V; V H be subspaces in Hilbert space such that V V ? = H . Let P be the projection onto V e ? and let along V be the orthoprojector onto V . For all g 2 H , we get

kg ? gkH kg ? PgkH
where

kg ? gkH ;

k = sup kPkk kH : k H 06=k2V

Riesz systems consisting of integer translates of a single function ' 2 L2 (R) are of particular interest here. More explicitly, we will assume that ' 2 L2(R) is a function such that the sequence ('( ? k))k2Z is a Riesz system in L2 (R). If V is the closed linear span of this Riesz system, then there exists a unique Riesz system in V biorthogonal to ('( ? k))k2Z which is a basis in V . In particular, for the given function ', one may construct a function ' 2 V , e such that the dual Riesz system of ('( ? k))k2Z is given by ('( ? k))k2Z. We mention that if ' has compact support e and is bounded, then ' has exponential decay at in nity always and in some cases compact support4 . e We shall make further assumptions on the functions ' and '. Indeed, we will assume that these functions e induce multiresolution analyses. Multiresolution analysis forms one of the key issues in wavelet theory; see.4,5,12 A multiresolution analysis (MRA) in L2 (R) is a sequence of subspaces (Vj )j2Z in L2 (R) with the following properties: For each j 2 Z, we have (i) Vj Vj+1 , (ii) f 2 Vj , f (2 ) 2 Vj+1 , and (iii)
T S
j Z j

2.2. Multiresolution analysis

(iv) j2ZVj L2 (R) dense, (v) there exists ' 2 V0 such that ('( ? k))k2Z is a Riesz basis in V0 . A function ' 2 L2 (R) which gives rise to a multiresolution analysis as above will be called a scaling function. Examples of such function are spline functions together with various dual functions1,3,4 . We remark that a biorthogonal pair of scaling functions leads to the de nition of a biorthogonal pair of wavelets and a (fast) wavelet transform. We refer to1,4,11 . Throughout this section, we will assume that ' induces a multiresolution analysis as described in the previous section, i.e., ' is a scaling function. We will also x a dual scaling function ' 2 L2(R). e n , we consider the (rotated and dilated) grid Given an orthonormal basis 1 ; : : : ; n in R
G
;d = f n X r =1

2 V = (0),

3. LEAST SQUARES METHOD AND DUAL METHOD

pr dr r j (p1 ; : : : ; pn ) 2 Zng;

where d1 ; : : : ; dn are positive real numbers. The standard basis in Rn will be denoted by e1 ; : : : ; en and the standard grid by G e accordingly. Observe that G ;d is the image of G e under permutation, rotation and dilation. In fact, G ;d = DRQG e , where Q is an n n permutation matrix, R is an n n rotation matrix and D is an n n diagonal matrix with diagonal d = (d1 ; : : : ; dn )T .

In order to deal with functions on Rn for n 2, we shall construct multivariate functions by means of products of univariate ones. Indeed, given an orthonormal basis 1 ; : : : ; n in Rn , we shall write
;d

(y) =

n Y r =1

d?1=2 '(d?1 hy; r i); almost all y 2 Rn ; r r

where ' is a scaling function. In the same fashion, one de nes e ;d using the dual scaling function '. The next two e lemmas justify that ;d will be called a multivariate scaling function. Lemma 3.1. The system ( ;d( ? p))p2G ;d is a Riesz system in L2(Rn ). The dual Riesz system is given by ( e ;d( ? p))p2G ;d . In particular, h ;d ( ? p); e ;d( ? q)iL2 (Rn) = p;q : The Riesz bounds of the system are given by An ; B n . A multiresolution analysis in L2 (Rn ) associated with a grid G ;d is a sequence of subspaces (V ;d;j )j2Z in L2(Rn ) with the properties: For each j 2 Z, we have (i) V ;d;j V ;d;j+1 , (ii) f 2 V ;d;j , f (2 ) 2 V ;d;j+1 , and (iii)
T S
j Z

(iv) j2ZV ;d;j L2(Rn ) dense, (v) there exists 2 V ;d;0 such that ( ( ? p))p2G ;d is a Riesz basis in V ;d;0.
Lemma 3.2. If ' induces a multiresolution analysis in L2 (R), then
2 n ;d induces a multiresolution analysis in L (R ) associated with the grid G ;d . We shall introduce a multivariate interpolation operator which incorporates the Riesz system induced by the multivariate scaling function. De ne L ;d : 2 (G ;d ) ! L2 (Rn ) by

2 V

;d;j

= (0),

L ;d(ap )p2G ;d (y) =

X
p G ;d

2

ap

;d

(y ? p)

e for almost all y 2 Rn . The operator L ;d is de ned in the same way using the dual scaling function. Observe that by Lemma 3.1, we get

An k(ap )p2G ;d k22(G ;d ) kL ;d(ap )p2G ;d k2 2 (Rn) B n k(ap )p2G ;d k22 (G ;d ) :  L 
This implies that L ;d : 2 (G ;d ) ! L2 (Rn ) is a bounded injective operator with closed range. It follows that L ;d has a bounded left-inverse. Indeed, the adjoint operator L ;d : L2(Rn ) ! 2 (G ;d ) of L ;d is given by (L ;dg) = hg;
?
;d

( ? p)iL2 (Rn) p2G ;d ; g 2 L2 (Rn );

and now we can identify a left inverse of L ;d. e Theorem 3.3. The operator L ;d : 2 (G ;d ) ! L2(Rn ) has L ;d : L2 (Rn ) ! 2(G ;d ) as a left inverse. Proof. Note that Lemma 3.1 and the continuity of the inner product implies
e e hL ;dL ;d(ap )p2G ;d ; (bq )q2G ;d i2 G ;d = hL ;d(ap )p2G ;d ; L ;d(bq )q2G ;d iL2 Rn =
( ) ( )

h

X

p G ;d

2

ap

;d

( ? p);

X
q G ;d

2

bq e ;d( ? q)iL2 (Rn) =

X
p G ;d

2

ap bp :

e In the particular case when ;d = e ;d, i.e., in the orthonormal case, we obviously get L ;d = L ;d and by Theorem 3.3, we see that in this case L ;dL ;d = I2 (G ;d ) . In the general situation, the self-adjoint operator L ;dL ;d will not be the identity, although the following holds true. Lemma 3.4. The operator L ;dL ;d : 2 (G ;d ) ! 2 ( GG ;d) is a strictly positive, hence boundedly invertible, operator which satis es the estimates (in the sense of self-adjoint operators)

An I L ;dL ;d B n I:
(fp )p2G ;d = (L ;dL ;d)?1 L ;dg makes sense and provides the least-squares solution to the equation Observe that the expression (2) (3) (4)

L ;d(fp )p2G ;d = g
where g 2 L2 (Rn ) is a given function. Another approximate solution to (3) is given by
e (fp )p2G ;d = L ;dg:

L ;d to the right hand side of (2), we get L ;d(L ;dL ;d)?1 L ;dg. The operator ;d = L ;d(L ;dL ;d)?1 L ;d is the orthogonal projection onto ran L ;d. The least squares solution method (2) produces an error kg ? ;dgkL2(Rn) which equals zero whenever g 2 ran L ;d. In general, this error is minimal among all attainable errors, since kg ? ;dgkL2(Rn) equals the distance between g and ran L ;d. e e On the other hand, if we apply L ;d to the right hand side of (4), we get L ;dL ;dg. The operator P ;d = L ;dL ;d is a not necessarily orthogonal projection onto ran L ;d. This dual solution method produces an error kg ? P ;dgkL2(Rn) . Since g ? ;dg ? ;dg ? P ;d g, we get kg ? P ;dgk2 2(Rn) = kP ;dg ? ;dgk2 2(Rn) + kg ? ;d gk2 2(Rn) : L L L
It is immediate that the error caused by the least squares solution method is majorized by the error caused by the dual method. Moreover, the di erence between the two errors can be measured by kP ;dg ? ;dgkL2(Rn) . Theorem 3.5. Let ;d and e ;d be multivariate scaling functions as de ned before which satisfy the biorthogonality condition h ;d( ? p); e ;d( ? q) = p;q ; p; q 2 G ;d ; and let be de ned as in Theorem 2.2. The closed linear span of ( ;d( ? p))p2G ;d is denoted by V ;d and the closed e linear span of ( e ;d( ? p))p2G ;d is denoted by V ;d. The orthoprojector onto V ;d is given by ;d, and the projection e ? reads P ;d . We get for g 2 L2 (Rn ), onto V ;d along V ;d

We will now look at equation (3) and its approximate solutions (2) and (4) more closely. If we apply the operator

kg ?

;d

gkL2(Rn) kg ? P ;dgkL2(Rn)

n

kg ?

;d

gkL2(Rn) :

We shall now apply the considerations in the preceding section to interpolation between dilated and rotated Cartesian grids. The situation is as follows. A collection of data (complex numbers) (aq )q2Ge is given and a function g 2 L2(Rn ) is attached to this collection by means of interpolation. Indeed, we assume that

4. INTERPOLATION BETWEEN ROTATED GRIDS
X
q Ge

g(y) =

2

aq e (y ? q); almost all y 2 Rn :

Note that we have constructed a function g such that

aq = hg; e e ( ? q)i; q 2 G e :
In words, the function g reproduces the original data when inner products are taken with translates along the standard grid G e of the dual multivariate function e e . We now search for a collection of data (fp )p2G ;d associated with the dilated and rotated grid G ;d . This data should reproduce the function g as an interpolant of (fp )p2G ;d . In general, this is not possible. Therefore, we shall write

f (y) =

X
p G ;d

2

fp

;d

(y ? p); almost all y 2 Rn ;

and we shall try to minimize kf ? gkL2(Rn) . In terms of interpolation operators, we have the following situation. Given the data (aq )q2Ge 2 2 (G e ), we construct g = Le (aq )q2G e . The least squares solution to

L ;d(fp )p2G ;d = g = Le (aq )q2Ge provides us with data (fp )p2G ;d on the rotated and dilated grid, such that f = L ;d(fp )p2G ;d minimizes kf ? gkL2(Rn)
among all possible solutions. The least squares solution can be obtained with the least squares method (fp )p2G ;d = (L ;dL ;d)?1 L ;dg = (L ;dL ;d)?1 L ;dLe (aq )q2G e : Here, the least squares method involves the solution of a linear system with a sparse matrix. With complete factorization the number of ( oating point) operations required for the solution is at least proportional to the square of the number of data, see e.g.7 . In order to reduce the number of computations, the dual method
e e (fp )p2G ;d = L ;dg = L ;dLe (aq )q2G e

is proposed. Moreover, the solution of the dual method can be taken arbitrarily close or equal to the least squares e solution. We describe the dual method. The operator L ;dLe : 2 (G e ) ! 2(G ;d ) is applied to the data (aq )q2G e as follows: Observe that for p 2 G ;d ,
e L ;dLe (aq )q2G e Z X
p

=

Z

Rn

Le (aq )q2Ge (y) e ;d(y ? p) dy =
X
q Ge

The in nite matrix M = (Mp;q )p2G ;d ;q2Ge , with matrix elements given by Mp;q = h e ( ? q); e ;d( ? p)iL2 (Rn) provides the solution to the dual method as follows: (fp )p2G ;d = M (aq )q2G e : In order to describe the number of calculations required to perform the dual method, we consider the case when the dataset (aq )q2G e is nite. In particular, we denote the unit cube in Rn by K = ?1; 1]n and assume that aq = 0 for q 62 Q K where Q is a xed positive integer. Further, we will assume that the underlying scaling functions ' and ' have compact support. Indeed, assume that ; e > 0 arenchosen such that supp Pn ? ; ] and supp ' ? e; e]. e ' e P Let (p1 ; : : : ; pn ) 2 Zn and (q1 ; : : : ; qn ) 2 Zn, then q = r=1 qr er 2 G e and p = r=1 pr r dr 2 G ;d . The matrix element Z Y n e r '(hy; er i ? qr )d?1=2 '(d?1 hy; r i ? pr ) dy Mp;q = r is nonzero, only if
Rn r=1

Rn q2G e

aq e (y ? q) e ;d(y ? p) dy =

2

aq h e ( ? q); e ;d( ? p)iL2 (Rn) :

jhy; er i ? qr j

; jd?1 hy; r i ? pr j r

e

; r = 1; : : : ; n:

We get
n X r =1

ky ? qk =
2

n X r =1
2

jhy; er i ? qr j
n X r =1

2

2

n;
e2

ky ? pk =
2

jhy; r i ? dr pr j =
T

d2 jd?1 hy; r i ? pr j2 r r

kdk ;
2

where kdk is the Euclidean norm of d = (d1 ; : : : ; dn ) . This provides As a result, we see that

0 = k(y ? q) + (q ? p) + (p ? y)k kp ? qk ? ky ? pk ? ky ? qk kp ? qk ? ( n + ekdk):

p

kp ? qk ( pn + ekdk): p e We may conclude that if kp ? qk > ( n +pkdk), then Mp;q = 0. This means that the matrix M has a band-like structure with band width majorized by ( n + ekdk). The number of data N , i.e., the number of gridpoints of
G
;d

in Q K , is proportional to (2Q)n ?1 ,p where = d1 : : : dn . The number of computations required to calculate (fp )p2G ;d is majorized by the number N ( n + ekdk) which is of the same order as the number of data involved.

In this section we consider an application of the theory described in the previous section. We propose a testproblem and sketch the outline of a numerical algorithm to solve it.

5. AN APPLICATION IN SEISMIC DATA PROCESSING

5.1. An idealized testproblem

We propose a simple synthetic shotrecord that is depicted by the picture on the left of Figure 1. It shows a simple

50

100

100

200

150

300

200

400

250

500

300

600
50 100 150 200 250 300

100

200

300

400

500

600

Figure 1. Idealized testproblem (frequency= 4)
synthetic dataset on a grid of dimensions 300 300. The dataset consists of low-amplitude hyperbolic events on which a high-amplitude function g has been superimposed. The picture is a caricature of a seismic shotrecord. The vertical axis corresponds to time, the horizontal axis to a line at the earth's surface. The hyperbolic events would correspond to measurements of seismometers recording waves excited at the earth's surface and re ected at interfaces (i.e. where the impedances changes rapidly). The challenge consists of segregating the high-amplitude distortion from the low-amplitude events of interest. The problem might serve as a benchmark for testing numerical algorithms and is therefore made available19 at WWW. The dimensions of the grid, the angle and frequency of the distortion (and its width), and other relevant parameters can all be controlled.

5.2. A segregation template

We investigate the applicability of wavelet transforms as an alternative to the Fourier transform. Results of earlier research in this respect are reported in2,6,8 . What's new in our method is determined by: rotation/interpolation, the transform (of wavelet-type) and the ltering. We devise numerical algorithms for the solution of problems as described in Section 5.1. The methods we discuss can be characterized by the following template with seven subsequent steps:

step 1 step 2 step 3 step 4 step 5 step 6 step 7

Flatten data. Rotate data. Transform data. Adapt (limit, mute, interpolate) appropriate coe cients, i.e. ltering. Apply backtransform. De-rotate data. De- atten data.

Apart from step 2 (and its counterpart step 6) this is quite a conventional framework. Most or all present methods t in this framework where often fast Fourier' is used as transform (step 3). Step 1 (and its counterpart step 7) may be a simple method to compensate for the exponential decay in time of the signal, but more sophisticated preprocessing at this point is quite common. Step 2 (and its counterpart step 6) is new. For reasons to be explained later, we rotate the data such that they more or less align with either horizontal or vertical gridlines. The actual ltering takes place at step 4, it is the heart of the algorithm. In a way, the combination of steps 1{3 can be seen as one big transform with the aim of making the undesired components explicit. The purpose of steps 1 { 3 is to transform a signal step by step in such a manner that the coe cients with respect to the new basis of representation can be divided into a set corresponding to the desired components and another set corresponding to the undesired components. Ideally, we then could simply mute the coe cients going with the undesired components. By numerical experience our ways of adapting (step 4) evolved from A to C :

A c := max(min(c; upperbound); lowerbound); B if c satis es mute-criterion then c := 0 (mute c); C if c satis es mute-criterion then replace c by interpolation between nearby coe cients that do not satisfy the
mute-criterion. Here A is called limiting, B is called muting, C is called interpolation. The values lowerbound' and upperbound' are de ned in a problem-dependent way. Often A turns out to be too crude: the groundroll remains too large. This is improved by B which is more radical by muting coe cients completely. An important di erence of C with both A and B is that the coe cients which need to be adapted are replaced by interpolation between remaining nearby coe cients that (are supposed to) correspond to the wanted components. The mute-criterion in B, C still needs to be lled in. It can be simple like: mute-criterion = (c > upperbound or c < lowerbound) be it that the values lowerbound' and `upperbound' still need to be determined. More details and results can be found in20 .

5.2.1. Rotation

Why rotate the data? As with every transform in the context of seismic data processing the aim is to make a better separation between the unwanted components and the events of interest. For the two-dimensional waveletdecomposition, it has been shown2 that if the angle between groundroll and the gridlines is within well-de ned bounds a bias exists for distribution of the energy towards either horizontal detail or vertical detail wavelet coe cients. These ndings supported our research to investigate the e ects of rotation. Figure 1 shows how a simple synthetic dataset on a 300 300-grid is rotated (interpolated) onto a 600 600-grid (both grids are uniform). The \right" angle of rotation is chosen visually, in practice we need to rely on a numerical approach to detect directional bias in a dataset. In order to add a touch of realism, the function g is not completely aligned with the grid. Note that the rotated grid wants a smaller meshwidth than the original if we do not allow for loss of information. Generally, we need to be concerned about costs and accuracy of the rotation, see Section 4 and20 .

5.2.2. A numerical result

We apply the above segregation template to the idealized testproblem of Section 5.1. Steps 1 and 7 are omitted. We choose a compactly supported two-dimensional separable biorthogonal spline wavelet and use a corresponding decomposition of 6 levels. Muted coe cients at step 4 are replaced by interpolation (C ). Figure 2 shows the results before and after de-rotation (step 6). If we compare Figure 2 with Figure 1 we observe a reasonable result: the
Corrected, level = 6

100

50

200

100

300

150

400

200

500

250

600 100 200 300 400 500 600

300

50

100

150

200

250

300

Figure 2. Numerical segregation result for the idealized testproblem (frequency= 4)
high-amplitude distortion has been removed e ectively and (nearly) all of the hyperbolic events retain coherence. The notion of a biorthogonal pair of Riesz systems which induce multiresolution analyses can be de ned in the multivariate case. This leads to a fast rotation algorithm which uses biorthogonal (spline) wavelets. By employing rotation and wavelet decomposition we were able to devise and demonstrate a numerical algorithm that can segregate a high-amplitude distortion with directional bias from low-amplitude hyperbolic events. This work was supported nancially by the Technology Foundation (STW), project no. CWI44.3403.

6. CONCLUSIONS

ACKNOWLEDGMENTS

1. C.K. Chui, An introduction to wavelets, Academic Press, Boston, 1992. 2. J.C. Cohen and T. Chen, \Fundamentals of the Discrete Wavelet Transform for Seismic Data Processing," Preprint, 1993, www.mathsoft.com/wavelets.html 3. A. Cohen, I. Daubechies, J.-C. Feauveau, \Biorthogonal bases of compactly supported wavelets," Comm. Pure Appl. Math. 45: 485{560, 1992. 4. A. Cohen, R.D. Ryan, Wavelets and multiscale signal processing, Chapmann & Hall, London, 1995. 5. I. Daubechies, Ten lectures on wavelets, SIAM-CBMS Series 61, Philadelphia, PA, 1992. 6. Bingwen Du, Larry Lines, \Wavelet ltering of tube waves in the Glenn Pool crosswell seismic survey," M.U.S.I.C. Annual report on wavelet application, Memorial University of Newfoundland, 1998. 7. I.S. Du , A.M. Erisman, J.K. Reid, Direct Methods for Sparse Matrices, Monographs on Numerical Analysis, Clarendon Press, Oxford, 1986. 8. Liu Faqi, M.M. Nurul Kabir and D.J. Verschuur, \Seismic processing using the wavelet and the Radon transform," J. of Seismic Exploration 4, 375{390, 1995. 9. M. Holschneider, Wavelets: An analysis tool, Clarendon Press, Oxford, 1995. 10. G. Kaiser, R.F. Streater, \Windowed Radon transforms, analytic signals, and the wave equation," Wavelets: A tutorial in theory and applications, C.K. Chui (ed.), pp. 399{441, Academic Press, 1992. 11. A.K. Louis, P. Maass, A. Rieder, Wavelets: Theory and Applications, John Wiley & Sons, Chichester, 1997. 12. S.G. Mallat, \A theory for multiresolution signal decomposition,"IEEE Tr. Pattern Anal. Machine Intelligence 11: 674{693, 1989. 13. T. Takiguchi, \On invertibility of the windowed Radon transform," J. Math. Sci. Univ. Tokyo 2: 621{636, 1995. 14. M. Unser, \Quasi-orthogonality and quasi-projections," Appl. Comp. Harm. Anal. 3: 201{214, 1996. 15. M. Unser, \Approximation power of biorthogonal wavelet expansions," IEEE Tr. Signal Proc. 44(3): 519{527, 1996. 16. M. Unser, A. Aldroubi, \ A general sampling theory for nonideal acquisition devices," IEEE Tr. Signal Proc. 42 (11): 2915{2925, 1994. 17. R.M. Young, An introduction to nonharmonic Fourier series, Academic Press, New York, 1980. 18. A.L. Warrick, P.A. Delaney, \A wavelet localized Radon transform," SPIE Proceedings 2569, Wavelet Applications in Signal and Image Processing III, pp. 632{643, 1995. 19. P.M. de Zeeuw, \A benchmark problem representing a strongly simpli ed and synthetic shotrecord", 20. P.M. de Zeeuw, R.A. Zuidwijk, \Numerical methods for decomposition of 2D signals by rotation and wavelet techniques", CWI-report PNA-R98xx (in preparation), 1998. 21. R.A. Zuidwijk, \The wavelet X-ray transform", CWI-report PNA-R9703, 1997. 22. R.A. Zuidwijk, \The discrete and continuous wavelet X-ray transform," SPIE proceedings 3169, Wavelet Applications in Signal and Image Processing V, pp. 357{366, 1997. 23. R.A. Zuidwijk, P.M. de Zeeuw, \The fast wavelet X-ray transform," CWI-report PNA-R98xx (in preparation), 1998.
www.cwi.nl/ftp/pauldz/Demos/SimpleSet

REFERENCES

ÔÞÖúÉÌÁ´½Ó

¸ü¶àÏà¹ØÎÄÕÂ£º
Ó¢ÎÄÎÄÏ×¾äÐÍÕûÀí
The parameters u and s are scale factor and ...Heisenberg box, it is usually used to analyze ...far there has been no satisfactory fast algorithm...
Overview of Fourier Transformation
wavelet multi-scale analysis of signals thought ...algorithm and using its fast characteristics, for ...It is an ideal tool for signal time frequency ...
Ð¡²¨ÄÜÁ¿ìØµÄÑÐ¾¿¼°ÆäÔÚµçÁ¦ÏµÍ³Ð³²¨¼ì²âÖÐµÄÓ¦ÓÃÍâÎÄ...
fast algorithm of wavelet analysis, that is ...scale energy leakage and produce frequency aliasing...used to locate a system under certain conditions,...
Êý×ÖÍ¼Ïñ´¦ÀíÓ¢ÎÄÎÄÏ×·­Òë²Î¿¼
we used two different types of gray-scale ...Wavelet Transform based algorithm and get the ...time is about 2 seconds, operation is very fast...
Finite element analysis system development present ...
scales smaller but to use nimble, the price is...uses the non-linear finite element algorithm to ...time finite element analysis software research key ...
Medium and small-scale analysis of financial data(ÖÐÐ¡¹æÄ£µÄ...
Medium and small-scale analysis of financial data(ÖÐÐ¡¹æÄ£µÄ½ðÈÚÊý¾Ý·ÖÎö)-ÍâÎÄ·­Òë_¼ÆËã»úÈí¼þ¼°Ó¦ÓÃ_IT/¼ÆËã»ú_×¨Òµ×ÊÁÏ¡£Medium and small-scale analysis of...
Analysis of the use of computer image technology i
Analysis of the use of computer image technology ... With the increasing scale of chemical industry, ...operation is to achieve a fast and good ...
lecture 3 notes - time and frequency analysis
analyze systems in both the time and frequency ...Scaling and shifting is described by the systems ...Use partial fraction expansion to find the inverse...