Документ взят из кэша поисковой машины. Адрес оригинального документа : http://imaging.cs.msu.su/pub/2000.NFAO.Krylov.en.pdf Дата изменения: Sun Mar 15 23:00:00 2009 Дата индексирования: Sat Apr 9 22:45:08 2016 Кодировка:

NUMERICAL PROJECTION METHOD FOR INVERSE FOURIER TRANSFORM AND ITS APPLICATION

Andrey S. Krylov and Anton V. Liakishev

Faculty of Computational Mathematics and Cybernetics, Moscow State University, Vorobievy Gory, 119899 Moscow, Russia e-mail: kryl@cs.msu.su

Abstract Numerical pro jection method of the Fourier transform inversion from data given on a finite interval is proposed. It is based on an expansion of the solution into a series of eigenfunctions of the Fourier transform. The number of terms of the expansion depends on the length of the data interval. Convergence of the solution of the method is proved. The pro jection method for the case of the sine Fourier transform and the set of the odd Hermite functions being its eigenfunctions are examined and applied to numerical Fourier filtering. Keywords: Fourier transform, Hermite functions, Fourier filtering, pro jection method, convergence theorem, numerical method. AMS Classification (1991) : 42A38, 65R10, 33C45.

1


1

Intro duction

The Hermite functions as the eigenfunctions of the Fourier transform [1, 2, 3] are widely used in the Fourier analysis. Expansions into series of the Hermite functions in Fourier analysis applications can be found in image coding and analysis [4, 5], statistics [6], electrophysiology [7], Monte-Carlo computations [8], diffraction structure investigations [9, 10] and other areas of sciences. The Hermite functions have also attracted considerable interest in applications as real functions minimizing the Gabor uncertainty [11, 12] for the Fourier transform. The Gabor principle by analogy with the Heisenberg uncertainty principle is that the product of localization of a function in real (time) and reciprocal (frequency) spaces has a lower bound and the Gabor functions [13] are complex functions realizing this bound. In many cases we know that the data to be processed by the Fourier transform are not compact supported. Nevertheless, they are given on a finite interval only. The principal problem in this situation is to find a criterion to estimate the number of terms to be used for the expansion of the data into a series of the Hermite functions. The idea to relate the number of terms with the length of the experimental interval and thus to take into consideration the additional physical information on the data was first introduced in [14]. The aim of our paper is to justify this approach and to show some application of the method to the problem of the Fourier filtering. The outline of the paper is the following. We consider in §2 a pro jection method of solution of an inverse problem for linear operator equations including as a special case the problem of the Fourier transform inversion with the data given on a finite interval. A theorem of convergence of the solution for the method is proved. The properties of the Hermite functions are investigated in §3. The feasibility of the proposed method for inversion of the sine Fourier transform using approximation of solution by a series of the Hermite functions is shown. The application of the method of Fourier transform inversion in the Fourier filtering is illustrated in §4. The practical advantages of the proposed

2


pro jection method are discussed in §5.

2

Preliminaries

We consider the problem of finding of a function z from known u Az = u, A : L2 [ 0, ) L2 [ 0, ), (1)

where A is a linear continuous operator and the operator A has a complete orthonormal in L2 [ 0, ) system of eigenfunctions {i }, i=1, 2, . . . such that the absolute values of the eigenvalues {i }, are bounded and are bounded away from zero: > > |i | > > 0. (2)

The inverse problem (1) is well-posed [15] and has an unique solution z for the known u. If instead of the exact right side u(x) of the equation (1) we have a function u (x) and an error u - u then z - z
L2 [0,) 0 L2 [0,)

,

(3)

- 0. Here z is the solution of the equation Az = u .

In many cases the information on u(x) can be obtained on a finite interval [0, x
max

] only and the problem can be posed as Azxmax = u
xmax

,

A : L2 [ 0, ) L2 [ 0, x

max

].

(4)

In general, we have not the uniqueness of the solution zx and thus the problem is ill-posed. Consider now the case when instead of ux [0, x
max
max

max

for this problem

(x) we have u (x) on the interval

]. In this situation it is reasonable to find an approximation of the unique

solution of the well-posed problem (1) instead of finding of an approximation of one of the solutions of the problem (4). The lack of the information for x > x
necessary to find an approximate solution zxmax that zx max

can be considered as an "error" in the right side of the equation (1). So, it is
max

-z

where z is the solution of (1). 3

xmax 0

-

0,


Definition. Set of functions {i } is called norm-ordered with index k if
i L2 [ x,)

<

i+j L2 [ x,)

,

i = 1, 2, . . . ,

j = k , k + 1, . . . ,

(5)

for all x, where x > x1 0 and k 1 is a constant. If for some k > 0, > x1 0 we can reindex set of the eigenfunctions of our operator A so that the set i becomes norm-ordered with index k , then the algorithm of solution of our problem can be based on the following theorem. Theorem 1. Let z where d,i xmax c,i = xmax , i
xmax xmax nm

(y ) =
i=1

d,i i (y ), xmax

c

,i xmax

=
0

u

xmax

(x)i (x) dx,

i = 1, 2, . . . , nm ,

and nm is the maximum n, satisfying the condition
i L2 [ x
max

,)



n
1/2+

,

i = 1, 2, . . . , n,

(6)

> 0, > 0 ­ constants. Then
zx

max

-z

L2 [ 0,)

xmax 0

-

0.

Pro of. The second term in the right side of the inequality z
xmax

-z

L2 [ 0,)

zx

max

-z



L2 [ 0,)

+ z - z

L2 [ 0,)

tends to zero for 0. For the first term the following estimate has place
zx
max

zxmax 2

(y ) - z (y )) dy

1 /2

- z

L2 [ 0,)

= (
0

1 dy
/2





n

m

(c,i xmax
0 i=1

1 - c ) i (y ) i i

2

1 dy

/2

+



i=nm +1

0

1 ci i (y ) i

2

.

4


According to (5),(6), n

m

when xmax . Therefore the second term

tends to zero when xmax . At the same time,
0 i=1

n

m

(

c,i xmax

-

c i

1 ) i (y ) i

2

1 dy

/2

= 2 1
/2

=





n


m

x

max




i=1 0

u ( )i ( ) d -
0

0

1 u ( )i ( ) d i (y ) dy i
2

=

=



nm



1

/2

0

i=1x max

1 u ( )i ( ) d i (y ) dy i
1 /2

2 1
/2

1




n
m







1 2 ( ) d i

/2

2 u ( ) d

i (y ) dy



0

i=1

xmax nm

xmax 2

1 u 1/2+ nm Q.E.D.

dy

1 /2

i (y )
0 i=1



u - 0, n nm m

Theorem 2. For the coefficients d estimate has place

,i xmax

,

i = 1, 2, . . . , nm ,

the fol lowing

|d,i - di | + u xmax where 1 di = i Pro of.
xmax

L2 [x

max

,)

n

1/2+ m

,

u(x)i (x) dx.
0



|d

,i xmax

- di | =
0

(u(x) - u (x))i (x) dx +
x
max

u(x)i (x) dx i (x)

u(x)-u (x)



L2 [0,xmax ]

i (x)

L2 [ 0,x

max

]

+ u(x) n

L2 [x

max

, )

L2 [xmax ,)

+ u(x) Q.E.D.


L2 [xmax ,) 1/2+ m

,

5


3

Hermite functions

The Hermite functions are the eigenfunction of the Fourier transform [2]. The properties of this set of functions enables us to construct an algorithm of inversion of the Fourier transform described above in the general case. The Hermite polynomials are defined as [16, 17] Hn (x) = (-1)n e and satisfy the following relations Hn (x) = 2nH
n-1 x
2

dn e-x , dxn

2

(7)

(x),

(8)

Hn

+1

(x) = 2xHn (x) - 2nH

n-1

(x).

(9)

The Hermite functions are defined as follows n (x) = e They satisfy the recursion relation
n+1 -
x2 2

Hn (x).

(10)

(x) = 2xn (x) - 2nn

-1

(x).

(11)

The same definitions and relations for the orthonormal in L2 [0, ) Hermite polynomials and functions are ^ Hn (x) = 2n-1 n! (-1)
n

e

x2

dn e-x , dxn n ^ H n+1

2

(12)

^ H

n+1

(x) = x

2 ^ Hn (x) - n+1

n-1

(x),

(13)

^ n (x) = e

-

x2 2

^ Hn (x) =

1 2
n-1

- e n!

x2 2

Hn (x),

^

n L2 [ 0,)

= 1,

(14)

^

n+1

(x) = x

2^ n (x) - n+1 6

n^ n-1 (x). n+1

(15)


Lemma 1. The norms of the Hermite functions n satisfy the recursion relation n
2 L2 [0,R]

= 2n n

2 -1 L2 [0,R]

- n (R)

n-1

(R).

(16)

Pro of. Applying (8) and integrating by parts we obtain
R 2 n L2 [0,R]

=
0

e

-x2

1 Hn (x)Hn (x)dx = 2(n + 1)
R 0
2

R

e
0

-x2

Hn (x)dHn+1 (x) =

2 1 e-x Hn (x)H = 2(n + 1)

R

n+1

(x)
-R

-
0

H

n+1

(x) d e

-x

2

Hn (x) =

1 e = 2(n + 1) 1 - 2(n + 1)
R

Hn (r)H
-x2

n+1

(R)-
2

H
0

n+1

(x) -2xe

Hn (x)dx + e-x dHn (x) .

According to (9) we have
R R R

2xHn (x)H
0

n+1

(x)e-x dx =
0 R

2

H

2 n+1

(x)e
R

-x

2

dx+
0

2nH

n-1

(x)H

n+1

(x)e-x dx =

2

=
0

H

2 n+1

(x)e

-x

2

dx +
0

Hn

+1

(x)e

-x

2

dHn (x).

Therefore
2 n L2 [0,R]

1 = e 2(n + 1)

-R

2

Hn (R)H

n+1

1 (R ) + 2(n + 1)

R

H
0

2 n+1

(x)e-x dx.

2

Thus, for the Hermite functions we have n
2 L2 [0,R]

=

1 n (r) 2(n + 1)

n+1

(R ) +

2 n+1 L2 [0,R]

,

and the required result follows.

7


^ For the orthonormal Hermite functions n the recursion relation may be written as follows ^ n
2 L2 [0,R]

^ =

2 n-1 L2 [0,R]

1^ ^ - n (R)n-1 (R). 2n

(17)

We will show the idea of the method of the Fourier transform inversion with the case of the sine Fourier transform on a half-line. The odd Hermite functions (see Fig. 1) ^ n (x) = 2
n+1

(x)

(18)

are the eigenfunctions of the sine Fourier transform


n (x) sin xy dx = (-1)n
0

n (y ). 2

(19)

They satisfy the following recursion relation (x) = 2x2 - 4n - 3 (2n + 2)(2n + 3) n (x) - 2n(2n + 1) (2n + 2)(2n + 3) (x). (20)

n+1

n-1

Lemma 1a. The norms of the orthonormal Hermite functions n satisfy the recursion relation - 1 2R 2 ( R ) + n
2 n L2 [0,R]

= n

2 -1 L2 [0,R]

-
2 n-1

4n + 1 2n(2n + 1)

n (R)n-1 (R) +
2 L2 [0,x]

(R) .

(21)

The general behavior of the norms n
2 n+1 L2 [0,x]

as functions of x is shown
2 n L2 [0,x]

in Fig. 2. According to (21) there are intersections between . Nevertheless, the following theorem has place [18].

and

Theorem 3. The set of functions {n } is norm-ordered with index 2 for x1 = 1, that is
2 n+i L2 [x,)

>

2 n-1 L2 [x,)

,

i, n 1, x 1.

(22)

8


4

Pro jection scheme of Fourier filtering

The general scheme of the Fourier filtering for function smoothing looks as follows
F F ^ ^ ~ f (x) - f (k ) - f (k ) = f (k ) · G(k ) - f (x),
-1

(23)

where F is the Fourier transform and G(k ) is a "window" function. We will consider the case of the sine Fourier filtering F (f (x)) = Fs (f (x)) = 2


f (x) sin xk dx.
0

An example of a window function is the "natural" window: G0 (k ) = 1 for 0 max

, and G0 (k ) = 0 for k > k

max

. Its use blur the initial function

f (x), x 0 and add false oscillations to the result. These disadvatages are typical for all window functions used in practice [19]. The pro jection method enables us to construct an effective algorithm of the Fourier filtering to avoid these problems and to store the filtered function in a compact form. ^ Instead of the scheme (23) the filtering is performed by approximating of f (k ) with a finite series of the odd Hermite functions n concentrated on [0, k ~ f (x) =
nm kmax max

]

ci i (x),
i=0

ci = (-1)

i 0

^ f0 (k )i (k ) dk ,

(24)

nm = max where 1 > > 0 is a constant.

n|

i

2 L2 [0,k

max

]

1- ,

(25)

^ By this method the multiplication of f (k ) by a window function G(k ) in the ^ scheme (23) is replaced by expanding of f (k ) into the finite series of n Hermite functions concentrated on the selected interval. The n numerically found using Lemma 1a. Figs. 3-5 show an example of use of the pro jection method in the Fourier filtering. We added uniformly distributed error to the initial function f (x). The resulting function f (x) was treated by the Fourier filtering with the natural 9
m m

odd

value can be


window and with the pro jection method for kmax = 8. Thus, the sine Fourier ^ transform f (k ) of f (x) was cut at k = 8 (see Fig. 4) and was continued by 0 for k > k
max

. Then the inverse sine Fourier transform was calculated by the
m

^ pro jection method (the obtained approximation fp (k ) with = 0.1, n

= 17

is plotted with dotted line in Fig. 4) and by the explicit formula of the sine ~ Fourier transform inversion. The results by the pro jection method fp (x) and ~ by the natural window application fs (x) are shown in Fig. 3 along with initial functions f (x) and f (x). Figs. 5a, 5b show more detailed parts of Fig. 3.

5

Conclusion

The methods usually used for numerical sine Fourier transform inversion from data given on a finite interval are based either on the explicit formula for the inverse Fourier transform or on the solution of an integral equation of the first kind with a fixed upper limit of integration. The standard method based on the use of the explicit formula of the sine Fourier transform is equivalent to the expansion of the data into an infinite series of the odd Hermite functions with the assumption that the data is zero outside the experimental interval. In this case the use of the Hermite functions which are not concentrated on the experimental interval in practice adds noise to the result, because the coefficients of the expansion for these functions depend critically on the unavailable data outside our interval. The inverse problem for the integral equation of the first kind being an illposed problem is solved, for example, in [20], using the regularizing method of Tikhonov [15]. The idea to replace the initial problem which has a continuos dependence of the result on the initial data, to an unstable problem and to suppress then this unstability by a regularizing method does not look natural. It is preferable to use methods, which take into account the properties of the sine Fourier transform while no important information about the initial problem is rejected. The pro jection method was illustrated for the case of the sine Fourier trans10


form inversion. This method can be also used for the complex Fourier transform [21]. The pro jection method in the case of the Hankel transform [22] is based on the system of Laguerre functions. Use of Gauss-Hermite quadrature to find the coefficients of the expansion [3] allows us to construct fast algorithms for the proposed method. Finally, some of the features of the pro jection method can be mentioned: (i) Additional physical information is used to find the solution (length of the data interval), (ii) Smoothing of the result is performed without loss of the physical information, (iii) It is possible to test if the data quality is sufficient for the given length of the data interval, basing on the quality of the fitting of the data obtained, (iv) The approximation of the data and its Fourier transform are stored in a compact form. The authors are grateful to Professor A.M. Denisov for many valuable discussions and helpful suggestions.

REFERENCES
[1] N. Wiener, The Fourier Integral and Certain of its Applications, Cambridge, 1933. [2] E.C. Titchmarsh, Introduction to the Theory of Fourier Integrals, Oxford University Press, N.Y., 1937. [3] W.F. Eberlein, A New Method for Numerical Evaluation of the Fourier Transform, J. of Math. Analysis and Appl., 65, 1978, 80-84. [4] J.-B. Martens, The Hermite Transform - Theory, IEEE Trans. Acoust., Speech, Signal Processing, 38, 1990, 1595-1606. [5] J.-B. Martens, The Hermite Transform - Applications, IEEE Trans. Acoust., Speech, Signal Processing, 38, 1990, 1607-1618. 11


[6] R. Deutsch, Nonlinear Transformations of Random Processes, Prentice-Hall Inc., Englewood Cliffs, N.J., 1962. [7] L.R. Lo Conte, R. Merletti and G.V. Sandri, Hermite Expansions of Compact Support Waveforms: Applications to Myoelectric Signals, IEEE Trans. Biomed. Eng., 41, 1994, 1147-1159. [8] A.J. Chorin, Hermite Expansions in Monte-Carlo Computation, J. Comp. Physics, 8, 1971, 472-482. [9] A.S. Krylov and B.M. Shchedrin, Numerical Method for Radial Distribution Function Calculation, Crystallographya, 34, 1989, 1088-1092 (in Russian). [10] D.I. Svergun, A Direct Indirect Method of Small-Angle Scattering Data Treatment, J. Appl. Cryst., 26, 1993, 258-267. [11] D. Gabor, Theory of Communication, J. Inst. Electr. Eng., 93, 1946, 429457. [12] J.A. Bloom and T.R. Reed, An Uncertainty Analysis of Some Real Functions for Image Processing Applications, Proceedings of the IEEE Intern. Conf. on Image Processing, 3, 1997, 670-673. [13] H.G. Feichtinger and T. Strohmer, eds., Gabor Analysis and Algorithms. Theory and Applications, Birkhauser, 1997. [14] A.S. Krylov and A.V. Vvedenskii, Software Package for Radial Distribution Function Calculation, J. of Non-Cryst. Solids, 192&193, 1995, 683-687. [15] A.N. Tikhonov and V.Y. Arsenin, Solution of Ill-posed Problems, V.H. Winston, Washington, DC, 1977. [16] M. Abramowitz and I.A. Stegan, eds., Handbook of Mathematical Functions, McGraw-Hill, N.Y., 1969.

12


[17] P.K. Suetin, Classical Orthogonal Polinomials, Nauka, Moscow, 1976 (in Russian). [18] A.S. Krylov and A.V. Liakishev, An Inequality for the Norms of Hermite Functions on a Finite Interval, Vestnik Moskovskogo Universiteta, Seriya 15 - Vychislitel'naya Matematika i Kibernetika, No.1, 1999, 17-19 (in Russian). [19] J. Max, Mґ des et Techniques de Traitement du Signal et Applications etho aux Mesures Physiques, Tome I, MASSON, 1981. [20] N.V. Ershov, A.L. Ageev, A.V. Serikov, Regular Method in X-Ray Scattering Analysis of Amorphous and Liquid Metals, Phys. Stat. Sol., 121, 1984, 451-460. [21] A.S. Krylov, J.F. Poliakoff and M.Stockenhuber, J. of Physics: Cond. Matter, to appear. [22] A.S. Krylov and M.A. Spiridonov, Compact Presentation of Structure Dependencies of the Surface Layer of Liquid Copper, Rasplavi, N.4, 1993, 81-84 (in Russian).

13


Figure 1: Odd Hermite functions

Figure 2: Norms of the first 20 odd Hermite functions in consequent order (from the left to the right) 14


Figure 3: Sine Fourier filtering results

Figure 4: Sine Fourier filtering results

15


Figure 5: Details of Fig. 3

16