CA2359209A1 - Multiexponential signal processing method and apparatus - Google Patents

Abstract

For signal processing of transients such as multiexponential decays, a transform operator (e.g., matrix) is constructed by emphasizing linear resolution, rather than using fitting routines that attempt to find an estimate of an unknown model that fits data. Use of the transform operator t o process multiexponential signals produces outputs that are a better estimate of the unknown model or of some segment of the unknown model.

Description

Title: Multiexponential Signal Processing Method and Apparatus SPECIFICATION
CROSS-REFERENCE TO RELATED APPLICATION
This application claims the benefit of Provisional Application No. 60/116,362 filed January 19, 1999, incorporated herein by reference.
BACKGROUND OF THE INVENTION
This invention is concerned with signal processing of transients, and more particularly, with multiexponential signal processing.
In the recent REVIEW article entitled °Exponential analysis in physical phenomena~~, Review of Scientific Instruments, Vol. 70, No. 2, Feb. 1999, pages 1233-1257 (incorporated herein by reference), authors Andrei'A.
Istratov and Oleg F. Vyvenko consider various aspects of multiexponential analysis. It is apparent from the text of the REVIEW article, and from the over 300 literature citations, that exponential analysis in physical phenomena is a subject of considerable interest in many scientific and technological disciplines, including, inter alia, solid state physics, medicine, biology and biophysics, geophysics, optics, engineering, chemistry and electro-chemistry.
As pointed out in the REVIEW article, many physical phenomena are described by first-order differential equations whose solution is an exponential decay. The amplitude and time constant of the exponential decay carry information about the nature of the phenomenon being studied. As commonly happens, a number of exponential processes take place simultaneously, and equipment employed in analyzing such exponential processes (multiexponential decays) yields a signal which is a sum of a plurality of exponential components.
Obtaining useful information from multiexponential decays is not a simple task. Many methods of multiexponential analysis and many algorithms for this purpose have been proposed, but all of them have limitations.
U.S. Patent No. 5,517,115 issued May 14, 1996 to Manfred G. Prammer (incorporated herein by reference) proposes a method and apparatus for efficient processing of nuclear magnetic resonance (NMR) echo trains in well logging. A priori information about the nature of the expected signals is used in an attempt to obtain an approximation of a model using a set of pre-selected basis functions. A singular value decomposition (SVD) is applied to a matrix incorporating information about the basis functions, and is stored off-line in a memory. During actual measurement, the apparatus estimates a parameter related to the signal-to-noise ratio (SNR) of the received NMR echo trains and uses it to determine a signal approximation model in conjunction with the SVD of the basis function matrix. This approximation is used to determine, in real time, attributes of the earth formation being investigated.
Techniques such as those described in the Prammer patent rely heavily upon a priori information about the nature of the expected signals. In attempting to obtain a reliable estimate of an unknown model, i.e., parameters of a subject being investigated, such techniques rely on one or more numerical algorithms for multiexponential analysis which use so-called "fitting routines" in an attempt to fit data to the model.
Underlying the present invention is the discovery that substantially improved results can be attained without using fitting routines, and, instead, by emphasizing better linear resolution. As pointed in the REVIEW article, fitting routines work well only if the hypothesis of the number of multiexponentiah components is correct and an initial approximization is close to the true solution. The present invention, which does not use fitting routines, avoids the problems that are inherent in the use of such routines.

BRIEF DESCRIPTION OF THE INVENTION
An object of the present invention is to provide improved methods and apparatus for the analysis of transients and for obtaining useful information therefrom.
More particularly, an object of the invention is to provide improved multiexponential signal processing.
One aspect of the invention involves a computer-readable medium containing a set of coefficients that define a transform operator such as a matrix.
Another aspect of the invention involves a method of calculating a transform operator utilizing a plurality of resolution functions.
Another aspect of the invention involves a method of multiexponential signal processing in which multiexponential signals are sampled, and the above-mentioned transform operator is applied to the sampled signals.
Another aspect of the invention involves an apparatus for multiexponential signal processing that comprises a signal processor including the above-mentioned transform operator.
The present invention starts with the construction of appropriate transform operators. Once appropriate transform operators have been constructed, they are incorporated in signal processors of analytical instrumentation for processing data. Generally, such instrumentation includes a computer, as is well known in multiexponential analysis.
Multiexponential data signals are sensed or detected by conventional equipment and are input to the transform operator of the computer for signal processing. The signals may be applied in real time or they may be read out from a suitable storage medium.
Typically, the signals are applied to the transform operator in digital form after conventional sampling and analog-to-digital conversion. For example, digital samples of multiexponential decays may be obtained at equally-spaced instants in time,~beginning at or just after the start of a multiexponential signal.
The invention makes use of linear resolution to obtain a better estimate of an unknown model. The present invention provides the very desirable property of optimal linear resolution of the unknown model when plotted against the log of the time constant of the decay curves.

BRIEF DESCRIPTION OF THE DRAWINGS
The invention will be further described in conjunction with the accompanying drawings, which illustrate preferred (best mode) embodiments, and wherein:
Fig. 1 is a diagram showing a matrix with coefficients all . . . a~ for transforming data dl . . . dN to produce parameters ml . . . m~ of an estimate of an unknown model;
Fig. 2 is a flow chart showing a method for calculating the coefficients of a row of a transform matrix in accordance with the invention;
Fig. 3 is a table showing, inter alia, coefficients for three rows of a matrix, for T's of interest, namely 1, 10, 100;
Fig. 4 is a diagram showing data functions, resolution functions, noise response, and point spread functions for a matrix constructed in accordance with the invention where the noise gain (NG) is 1.0000;
Fig. 5 is a similar diagram for another matrix constructed in accordance with the invention for a noise gain of 3.1623;
Fig. 6 is a diagram showing resolution funcaions for four matrices constructed in accordance with the invention with different noise gains;
Fig. 7 is a diagram showing point spread functions associated with four matrices constructed in accordance with the invention for different noise gains;

Fig. 8 is a block diagram showing apparatus for processing data in accordance with the invention;
Fig. 9 is a diagram showing decay curves of an MRI; and Fig. 10 is a diagram showing MRI relaxation distributions for four matrices constructed in accordance with the invention with different noise gains.
DETAILED DESCRIPTION OF THE INVENTION
One of the principal objectives of the present invention is to provide signal processing of transients, such as multiexponential decays, producing outputs that are better estimates of an unknown model to be investigated.
Such estimates permit better interpretation of data, so that a user (researcher, physician, scientist, or engineer, for example) can obtain more accurate information as to the nature of the unknown model. Considered from one point of view, the invention may be looked upon as a better digital lens that provides improved resolution, just as a better optical lens provides improved resolution.
A first step in achieving the objectives of the invention is to construct an improved transform ,operator, conveniently in the form of a matrix. The manner in which a transform matrix may be constructed, pursuant to the invention, will be described in detail later. The actual transform operator will depend upon its intended application, for example, medical imaging or well logging.
For any application, several different transform operators may be constructed, to provide a user with greater flexibility.
The present invention requires an understanding of what is referred to in the art as estimating a solution of a linear inverse problem. The linear inverse problem is one of communicating to an interpreter what is known and what is not known about an unknown model. The almost universal practice in the prior art for estimating the solution of a linear inverse problem is to calculate one or more of an infinite number of estimates of the model which fit the data, i.e., which reproduce the data to within the noise that is present. Whenever possible, a priori information is used to choose which of the estimates to calculate.
However, a priori information may not be sufficiently available or may be suspect.
In applying the present invention to multiexponential decays, an estimate of an unknown relaxation distribution (model) is obtained by linearly resolving each point of the unknown relaxation distribution as precisely as possible within the limits of the noise. In general, such estimates do not reproduce the data. Rather, they are obtained by optimizing the linear resolution in a manner that will be described later.
The term "linear" used herein in connection with "resolution" has the following meaning:

Consider a transform operator A which transforms functions (or vectors) xl and x2 to yl and y2, respectively.
As equations, this is stated:
A x1 = yl and A x2 = y2.
The transform operator A is defined as "linear" if for any two real numbers a and b A ( a xl + b x2 ) _ a yl + b y2 .
A transform operator which is a matrix will always have linear resolution. It is possible to construct non-matrix transform operators which behave similarly to matrices for only a restricted set of functions (or vectors). Moreover, a matrix can be derived from a non-matrix transform operator that expresses this behavior.
The present invention involves the following relationship between a true model mT(y) and data dk (e. g., multiexponential decays):
dk - JmT ~~J) 9x ~~J) dy ( 1 ) I
where I is the appropriate domain of definition and gk(y) represents a data function, more particularly,. one of the N data functions which compose the data kernel of a linear transform. A data function is a mathematical representation of how a model is mapped to a particular data point. Equation (1) is referred to as "the forward problem."

Decaying systems generally have a point, usually defined as t=0 (where t is time), at which a system begins decaying. The system decays to a constant value, often zero, as t approaches infinity. In nuclear magnetic resonance (NMR), for example, the t=0 point is when an excitation pulse energizes a sample.
The multiexponential forward problem commonly used in data analysis and inverse theory has the form dx - ~T~e t~lTr ( 2 ) In this equation dk represents the data; mi represents the unknown model, and ~e'~k/n represents the data function, where T designates the time constant of the exponential decay.
Equation (2) can be expressed in the continuous form dk = ~~~a~T _'.~~ a t,,l'r l3) Jo i where 8~~ is the Dirac delta function.
Incidentally, as used herein, the range of indices such as i, j and k are determined by the equation in which they are used. For example in the equation R~1/) _ ~ bi9i(EI) ~ 4 ) the index i would be assumed to range over all the data functions 9;(y~ . which are normally numbered with positive integers starting at 1. But in the equation NG; _ ~~i (5) i the index i ranges over the rows of a transform matrix and the index j ranges over the columns of the transform matrix, which is the same range as the data functions in the previous equation.
It is desirable to plot the relaxation distribution (unknown model) versus the logarithm of T, and to determine the resolution of a relaxation distribution on a log scale.
Applying a change of variables y=ln(T) to equation (3), without simplification, yields d~ _ ~+~~~3(ev _ ev:)e ckar(evd~) -°°
Applying the Dirac delta function identity a(f (x)) _ ~ f,~~~~ a(~ - xi) .f ~' each, j(x;) = U ( ~ ) to equation (6) along with standard simplification yields d~ = J~~ maa(U - ~Js)e tk' "dy ( g ).
-°°
Substituting mT (U) _ ~ mia(J - ~Ji) ~ 9 ) s into equation (8) gives i~oo dx = J mT ~1I) a t''' ~ d~- ( 10 ) where mT(y) represents the unknown model, and e'~
represents the data function, which may be expressed as:
9~~t1) = e-t'~.' ". ( 11 ) Hereinafter, the form of the multiexponential transform in equation (10) will be termed the forward problem.
It should be noted that the data function approaches 1 as T approaches infinity for all values of tk. Pursuant to the invention, the data function of equation (11) has been found to be very useful in multiexponential analysis, but other data functions may be appropriate for other applications of the invention.
As stated earlier, a first step in achieving the objectives of the invention is to construct a 'transform operator, such as a transform matrix, which maps data (e. g.
decay signals) to unknown model parameters. Fig. 1 is a diagram showing a matrix with coefficients all . . . aiR for transforming data dl. . . dH to produce parameters ml . . . mk of an unknown model. The transform matrix is typically a matrix in which each row corresponds to a T of interest.
Selection of r's of interest and data points for initial coefficients will be guided by available information in the particular field in which the matrix is to be used.

In accordance with the invention, it is desired to obtain an estimate of an unknown model with linear resolution, and since the data functions used are linear functions of the model, the estimate of the unknown model can be calculated by multiplying the data by a matrix. The matrix is chosen so that each point of the estimate of the unknown model linearly resolves the corresponding point of the model as well as possible with an acceptable noise gain i.e., with optimal linear resolution. Since matrix multiplication is a linear operation, it yields an estimate which does not necessarily reproduce the data, but which does have linear resolution. Linear resolution is a desirable property of an estimate, because each point of the estimate resolves the corresponding point of the model in the same way, independent of any particular model.
A goal in constructing a transform matrix in accordance with the invention is to calculate linear combinations of the data functions which yield a resolution function that resolves as small a region of the unknown model as possible.
Accordingly, it is preferred that each resolution function corresponding to a row.of the matrix be characterized by optimal linear resolution. Preferred criteria for selecting coefficients of a transform matrix which yield optimal linear resolution are described later. Together, the resolution function and its noise gain give a concise formulation of the ambiguity of a point in an estimate of an unknown model from information given by data and corresponding data functions.
Calculation of a transform matrix proceeds row by row.
It is convenient to use the same noise gain for all rows of a transform matrix, but this is not required. As noted above, each row normally corresponds to a T of interest.
The rows can be ordered by increasing T from top to bottom of the matrix. A spacing of 16 z values per decade has been found to work well. As an example, T values may be calculated between O.ltm;a and lOt~X, where tm;n is the smallest time at which a decay signal is to be sampled and t~,~ is the largest time. Each row of the matrix corresponds to a resolution function that is centered on a z of interest.
It is presently believed that the best way to calculate the coefficients of a particular transform matrix (which, incidentally, may have only one row) is to solve a constrained minimization problem.
The information needed before the constrained minimization can be performed are (1) the data function, g;(y), for each data point, d;, (2) the standard deviation of the expected noise, Q;, of each data point, (3) the desired noise gain, NG, and (4) the z of interest, Tk, on which the resolution function is to be centered.
The constrained minimization problem to be solved is to minimize I in the equation:
1 = ~ ~~I!)d~i ( 12 ) WO 00!43906 PCT/IB00/00212 by varying bi, where R(3~1~ _ ~; bi9e(1I~ ( 4 i and where b; are trial coefficients of a transform matrix row and R(y) is a trial resolution function. More information on resolution functions can be found in Parker, Robert L. 1994; Geophysical Inverse Theory; Princeton, New Jersey; Princeton University Press, incorporated herein by reference.
The trial resolution function preferably complies substantially with the following constraints:
(i) R(y) > 0 [Nonnegative constraint]
(ii) R(yk) > 1 [Peak constraint]
(iii) R(y)monotonically decreases from yk [Monotonicity constraint]
N
( iv) ~ bi2ai < NG I +'~ ~ bt9t(y)dy [Noise gain constraint i-__L1. ~1 °° i=1 Constraint (iii) may be satisfied automatically for the data function used for the multiexponential problem. Thus, while it is stated as a separate constraint, it is to be understood that it may, in fact, be satisfied inherently.
To arrive at the final coefficients, air, for row i of the transform matrix, once the trial coefficients bi which minimize I have been found, it is highly desirable that the trial coefficients be scaled by a constant so that the area of their associated resolution function is unity on a log scale. This is accomplished by N
0.j = b~~ ~ ~ bk9k (y)d~l ( 13 ) -°° m Making the area of the resolution function unity insures that each point in the estimate of the unknown model is a local average.
Incidentally, for a data point dj, sampled on a decay curve at time t~, the data function previously described can be modified by a balancing function B(y), so that the data function becomes 9k(~) - B(U) ( 14 ) This modification of the data function has the effect of making the resulting transform matrix generate an estimate of the unknown model multiplied by B(y) but with the required noise gain. While many balancing functions are possible, a useful balancing function is B(y~ = e'~ . The choice w = 1 may be useful since it increases the amplitudes of the relaxation components at late time constants by a factor of r.
Constructing a transform matrix in which each resolution function corresponding to a row of the matrix is characterized by optimal linear resolution preferably involves a process termed constrained optimization. In the preferred embodiment of the invention, this process includes an outer loop, a middle loop and an inner loop. The outer loop converts the constrained optimization problem to a series of unconstrained multidimensional minimization problems using the simple penalty method described by Fletcher (Fletcher, R. 1987 Practical Methods of Optimization; Toronto; John Wiley & Sons) incorporated herein by reference. The middle loop converts the unconstrained multidimensional minimization problems into a series of one dimensional (1D) minimization problems via the conjugate gradient method described by Press et al. (Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P.
1992; Numerical Recipes in C: The art of scientific computing; 2nd Ed.; Cambridge UK; Cambridge University Press) incorporated herein by reference. The 1 dimensional minimization problems are solved using the golden section search or parabolic interpolation described by Press et al.
A desirable condition for each resolution function is that its area should be substantially minimized.
Calculating the area of the trial resolution function I to be minimized, given the trial coefficients, b;, can be accomplished using the equation:
m 5~
where Aiis the area of the data function, gi(y), given by the equation 2 5 ~ _ /'Yi ( 16 ) ~~ 9i(tJ)d~!

The integral's upper limit of infinity has been replaced by the value Y1, which is a large finite number. Experience has shown that Y1 = ln(1000*t~x), is a good choice.
The simple penalty method converts the constrained optimization problem to an unconstrained one via residues, rm. A residue is a measure of how much a particular constraint is violated. If a particular residue is greater than zero, then the corresponding constraint is considered satisfied. The more negative the residue, the more the constraint is violated. The constrained optimization problem then becomes the unconstrained one of minimizing IP
in the equation:
I p = ~ b~Ac + P ~ min(rn, p)2 ( T 7 ) s m where min(x,y) returns the minimum and x and y. The value of P starts small, and is then increased by a small factor, for example, 0.1, after each unconstrained minimization. The solution to one constrained minimization is used as a starting point for the next constrained minimization. The value of P is increased until IP has stabilized. The test for stability is whether the value of Ip changes by an accumulative factor of 10-6 over four consecutive increases in P.
Each of the non-negative, peak, monotonicity and noise gain constraints will usually have associated residues.
Before the residues can be calculated, the data functions and the resolution functions are discretized.

The data functions and the resolution functions are discretized by sampling them at e.g., 16 points per decade of T starting at T= 0.01 time units and continuing to z=
1000*t~ time units. The discretized data functions, g;(y), are represented as gil. The total number of points at which the functions are sampled is denoted L. The discretized trial resolution function is represented by R1.
Thus .Rt _ ~ bi9it ( 18 ) The residues of the nonnegative constraint, rl'~'s, are defined to be r =
(19) giving L residues.
The residues of the peak constraint, rP, are defined to be r~-j~p_1 (20) where lP is the index of the discretized trial resolution function at the z of interest. This yields one residue.
The residues of the monotonicity constraint are defined to be ~.ML=~_~_1 2>1>_lr (21) and rM~ _ .Ri - Ri+i lp >_ 1 > L -1 ( 2 2 ) These two equations together yield another L - 1 residues.
The residue of the noise gain constraint is defined to be (23) rNC = lltG ~ Asb~ - ~ b;

This equation yields one additional residue.
In total there are 2*L+1 residues to express the four constraints.
The conjugate gradient technique is explained by Press et al. The termination conditions used may be an accumulative factor change in Ip of no more than 10-' over 25 consecutive iterations. The total number of iterations can be limited to 3000*N, where N is the number of data points.
The one dimensional minimization problems are solved using the golden section search or parabolic interpolation described by Press et al. The 1D optimization is terminated if the functional value does not decrease by a certain fraction each step. The fraction can be set to 10-8, for example. A step limit of 1000 on each 1D optimization can be set to prevent infinite or unproductive near infinite loops.
The linear transform for the trial row coefficients to the trial resolution function would usually be computed a large number of times during each one dimension optimization and would be computationally demanding. However, since the relationship between the trial row coefficients and the trial resolution function is always linear, this linearity can be used to greatly reduce the number of times the trial resolution functions have to be calculated directly from the trial coefficients.
In one dimensional optimization, a scalar parameter c is varied along a line with direction dl; the direction being provided by the conjugate gradient method. Therefore, for the trial coefficients b;, optimization occurs along a line defined by bn = b° + cdi . ( 2 4 ) where b° is also supplied by the conjugate gradient method.
Because of the linearity relationship between b; and R1, a similar equation can be written for the corresponding trial resolution function R1, Ri = R~° ~- cdR ( 2 5 ) where d~9~t ( 2 6 ) and dR - ~ ~9~t. ( 2 7 ) Therefore, along each line of optimization, the linear transform from the trial coefficients to the trial resolution function only needs to be calculated twice, to calculate R1 and dl, instead of once for every value of c considered. This results in a many fold reduction in the one dimensional optimization time.
The effectiveness of the minimization can be checked by how close the value of the peak of the resolution function is to unity, since the minimum area should yield a peak value of exactly unity.
An upper limit on the noise gain, NG, of a particular point in an estimate can be applied by the noise gain constraint stated earlier, where the area of the resolution function is included because the area of the resolution function must be normalized (set equal to 1) before the noise gain is calculated.
Figure 2 is a flow chart for calculating the coefficients of a row of a transform matrix pursuant to the foregoing description, using a digital computer conventionally. Coefficients for successive rows of the matrix can be calculated in the same manner for each T of interest and for each noise gain deemed appropriate. For calculating the coefficients of any row, the coefficients of the preceding row can be used as starting points:., Constrained optimization methods may perform better if all the data have noise with a standard deviation of 1.
Fortunately, it is quite simple to modify the data functions so this is the case. Given the real data functions, gi(y), the adjusted data function, giA(y) is ~(11~ = 9t(ll)~Qt ( 2 8 ) The adjusted data functions are then supplied to the constrained optimization method along with the assumption that the standard deviations of noise at all the data points are 1. The final coefficients, a;j, produced by the constrained optimization based on the adjusted data functions, need to be corrected for the adjustment. The correction to the final coefficients is air = a~Q.i ( 2 9 ) If a data set has 100's, 1000's, or even more evenly spaced points on a decay curve, it is much more efficient computationally for calculating the coefficients of the transform matrix (and also for applying it to data) to average adjacent data points together to create a new group data point. The corresponding data functions must also be averaged together to get the grouped data function corresponding to the new grouped data points. The size of the groups is important. While, in general, it is best to have larger groups at later sample times, groups which are too large will reduce the resolution of the resolution functions. Groups which are too small are inefficient. A
logarithmic group size appears to be a good choice.
The equation G" = max(1, int(C 10"~D). ( 3 0 ) appears to produce good group sizes for C = 0.1, D = 10, where 1 is a positive integer starting at 1 and increasing to whatever value is appropriate. The function max(x,y) returns the maximum of x and y, and the function int(z) rounds z down to the nearest integer. For 512 data points the above equation gives group sizes of Occasionally, instead of measuring data at a point in time on a decay curve, data will be measured by averaging over a window between two points in time. In this case, the data function can be approximated by averaging together a larger number of sample points over the window.
Modification of the data functions will also be required if a trigger that initiates a decay curve is not a single point in time. For example, in time resolved fluorescence spectroscopy, the flash of light triggering the decay will last a finite length of time. If the intensity versus time for the flash is L(t) then each data point will be convolved with this function. The corresponding data functions will also have to be convolved with the same function.
Fig. 3 illustrates a transform matrix constructed in accordance with the invention. Three sets of coefficients for T=1, 10 and 100 are shown as columns so that they fit on the page. Each set has 48 coefficients calculated from 48 data functions. The first 32 data functions correspond to 32 sampling points that are evenly spaced by one time unit at times 1 to 32. The next 16 data functions are evenly spaced by 30 time units between time 62 and time 512. The standard deviation of the noise for the first 32 points is assumed to be 1, and the standard deviation of the noise of the next 16 is assumed to be 0.18257. The standard deviation of the noise of the last 16 is a factor of 1/sqrt(30) less than the first 32. This drop in the standard deviation of the noise would result if the cut-off frequency of a low pass filter before the analog-to-digital converter were dropped by a factor of 3.0 before the point was acquired. For the sake of simplicity, the matrix was constructed without the use of balancing functions or data function groups.
Fig. 4 illustrates data functions, resolution functions, noise response, and point spread functions for a matrix constructed in accordance with the invention where the noise gain is 1.000. The data function curves are for time samples at t=1, t=2 . . . t=N from left to right. Each resolution function corresponds to a set of data functions.
The abscissa in each diagram is in units of T on a log scale [ln(T)], and each resolution function is localized about a particular T value. As a result, each point spread. function tends to be localized about a particular T value.
Fig. 5 is a diagram similar to Fig. 4 but for a matrix constructed with a noise gain of 3.1623. It should be noted that in each of Figs. 3, 4, and 5 there are 48 data points (time samples).

A comparison of resolution functions produced by matrices with different noise gains is shown in Fig. 6.
Higher resolution is achieved with higher noise gains, but there is a trade-off between resolution and noise gain.
Greater noise tends to make the results achieved less reliable.
In addition to resolution functions, performance of a transform matrix can be judged using the PSF's and noise response. To judge the performance of a transform matrix, information is required as to corresponding time points on the decay curve at which data should be acquired as well as assumed noise at each data point.
Fig. 7 shows point spread functions associated with four matrices constructed in accordance with the invention for different noise gains. PSF's can be calculated using simulated data at required time points for monoexponential decays at T = 1,2,5,10,20,50,100, ... time units, assuming the initial amplitude of the decay curve is 1. No noise is added to the simulated data. It should be noted that in each of Figs. 6 and 7, and also on Figs. 9 and 10 to be described, there are only 32 equally spaced data points (with the same noise).
To calculate the noise response, several realizations of a decay curve which are purely noise are generated. The noise may be assumed to be Gaussian. The standard deviation of the generated noise is scaled so that the standard deviation of the noise of the first data point is 1. Then the sampled noise decay curves are multiplied by the transform matrix to obtain the relaxation distribution.
Several realizations of the relaxation distribution of the noise can be plotted to obtain a "feel" for the distribution.
An important characteristic at each point in the estimate of the unknown model is the standard deviation due to the noise in the data. If the noise in the data is uncorrelated, has a mean of zero, and all points have the same finite standard deviation, it can then be characterized by a single standard deviation. The noise gain for each point can be defined to be the standard deviation of the point divided by the standard deviation of the noise in the data and can be calculated directly from the coefficients Q
NG= ~ai i (31) W
Once a transform matrix has been constructed, experiments can be performed in which the spacing of the data acquisition points, the number of data acquisition points, and the grouping of data acquisition points are changed. The results of these experiments can b~ compared by observing the resolution functions and/or the PSFs that are produced. As a rule of thumb, the linear resolution is proportional to the height of the resolution function provided that the resolution function has unity area on a log scale.

Each trial resolution function is localized about a point of interest in the unknown model by requiring the peak of the resolution function to have a value greater than unity at a z of interest and then by minimizing the area of the resolution function. The effectiveness of the minimization can be checked by how close the value of the peak of the resolution function is to unity, since the minimum area should yield a peak value of unity.
From the foregoing description, it is evident that each transform matrix is constructed so that each point of the estimate of the unknown model linearly resolves the corresponding point of the unknown model as well as possible within an acceptable noise gain. Since matrix multiplication is a linear operation, it yields an estimate which does not necessarily reproduce the data but does have linear resolution. With linear resolution, each point of the estimate resolves the corresponding point of the unknown model in the same way independent of any particular unknown model. A highly desirable property of linear resolution in providing an estimate of an unknown model is that it enables a human interpreter to obtain an intuitive "feel!' of what the data reveals and does not reveal about the unknown model.
Each resolution function gives a concise mathematical description of the linear resolution at each point of the unknown model and is independent of the unknown model and the data. Viewing the transform matrix as a digital lens, linear combinations of the data functions are calculated which yield a resolution function that resolves as small a region of the unknown model as possible.
After a transform matrix is constructed and selected, it is incorporated in a signal processor of a computer, as software or hardware, for example, as indicated in Fig. 8.
In this figure the INPUT represents a source of sampled digital signals such as multiexponential decays, e.g., NMR
decays obtained from well-logging. The DATA PROCESSOR and STORAGE are components of a conventional digital computer.
The OUTPUT MODULE may have conventional DISPLAY, NETWORK
INTERFACE, and PRINTER COMPONENTS, for example.
An estimate of an unknown model is calculated by multiplying data, e.g., multiexponential decays, by the transform matrix. Several transform matrices for different noise gains, for example, can be incorporated in the signal processor and accessed selectively to provide a user with greater flexibility.
As stated earlier, an object of the present invention is to provided a better estimate of an unknown model, from which useful information can be obtained. In the aforementioned Prammer patent, Figure 8 of the patent is a mapping of estimated NMR signal decay times into pore sizes of an investigated earth formation. The curve of Fig. 8 is an estimated relaxation distribution. The present invention provides a better estimate of the relaxation distribution, and also a better signal-to-noise ratio (SNR). Resolution functions produced by matrices constructed in accordance with the invention are superior, in terms of peak value and localization to resolution functions produced by matrices of the Prammer patent for the same data and noise.
In general, the present invention can be used to provide better estimates of an unknown model than the prior art, estimates from which more reliable and useful information can be obtained. These improved results are achieved by emphasizing linear resolution irrespective of whether data fits a model, and, in fact, without any attempt to fit data to a model. The invention is particularly useful in multiexponential signal processing, such as in the analysis of multiexponential decays.
As mentioned earlier, with particular reference to the REVIEW article, multiexponential analysis is useful in a host of scientific and technological applications. One of those applications, NMR, such as NMR in well-logging, has already been discussed. Another application is MRI
(magnetic resonance imaging). Fig. 9 shows decay curves in an MRI application. The decay curves represent pixels taken from a series of 32 magnetic resonance images of the brain of a multiple sclerosis patient. The TZ relaxation data were acquired with a lOms sample time out to 320ms. The strength and decay of the signals contain valuable information about the tissues. Some of the major tissue types of interest in multiple sclerosis are normal appearing white matter (NAWM), cerebral spinal fluid (CSF) and lesions, which can be classified as chronic or acute. Fig.
shows relaxation distributions yielded by applying to the decay curves of Fig. 9 transform matrices of the invention with various noise gains. The larger the noise gain of the 5 estimate, the wider the range of relaxation estimates available. This is a characteristic which shows up in the resolution functions and is due to the fact that the larger the noise gain the more likely there is a feasible resolution function available.

10 Estimating the noise in the estimates shown in Fig. 10 can be accomplished in several ways. The first is to estimate the noise in the data and multiply by the noise gain for each transform matrix. The ideal way to measure the noise in the data is to repeat the measurements a large number of times and calculate mean, standard deviation and covariance. These statistics can then be propagated through transform matrices using standard statistical procedures.
Unfortunately, the measurement of the decay curves takes about 20 minutes to complete on a patient, so large numbers of repetitions are impractical.
Another way to measure the noise in the data i$ to estimate the standard deviation from previous measurements using the same instrument. This is not always reliable because the noise can vary from sample to sample. For the data in Fig. 9, previous measurements of noise gives the standard deviation of about 10 scanner units. After multiplying the standard deviation by the noise gain, the predicted standard deviation for the noise in Fig. 10 is 100.
A third way to estimate the noise is to consider that while the noise gain increases by a factor of 10 in Fig. 10, the resolution gain increases by only 1.4 for relaxation rates around five sample times. In~Fig. 10, a portion of the signal between 30 and 200ms increases proportionately to the noise gain. This strongly suggests that the portion is due to random uncorrelated additive noises in the data. It is possible then to measure the standard deviation of the noise between 30 and 200ms and work back to the noise in the data.
Further applications of the invention will be apparent to persons of ordinary skill in the art. For example, decaying sinusoids are a common problem in inverse theory.
A decaying sinusoid can be handled if it is band pass filtered at the particular bandwidth of interest and then the magnitude of the decay curve taken. The relaxation distribution of the magnitude decay curve can then be calculated.
Any row of a transform matrix, since it corresponds to a resolution function, can be applied to a time series in the same way as digital filters. The "relaxation" digital filter will be useful for applications such as ultrasound and radar. Using quadrature detection, it would be possible to measure the magnitude of a reflection off of an interface. If the signal from an interface of interest oscillated for a while after the initial sound wave had passed, the decay time of the oscillation would reveal information about the interface. Applying various relaxation digital filters to the ultrasound time series would allow characterization of the interfaces. If it were desirable to detect a particular range of relaxation time, selection of the coefficients of a relaxation digital filter would give a resolution function with a more boxcar-like shape.
Other applications of the invention include photoluminescence and time-resolved fluorescence spectroscopy of biological and other types of samples, electrical signals radiating from ore bodies in geophysical exploration and acoustic, electrical and electromagnetic decay processes. It is possible to design a transform matrix which handles data integrated over a window or unevenly spaced in time by modifying the data functions.
Principles of the invention can be applied to problems which data functions other than decay curves. If the data functions are cosine functions, resolution functions can be generated for low pass, band pass and high pass .filters as well as windows for discrete Fourier transforms. A low pass filter can be designed by requiring maximum area below the cutoff frequency, minimum area above the cutoff and an optional requirement of monotonicity to eliminate wiggles..
The bounds on all parts of the filter would be 0 and 1. A

limit on broadband noise gain could also be imposed to improve the robustness of the filter.
While preferred forms of the invention have been shown and described,' these forms are intended to be exemplary, not restrictive, and it will be apparent to those skilled in the art that modifications can be made without departing from the principles and spirit of the invention, the scope of which is set forth in the appended claims.

Claims (21)

WHAT IS CLAIMED IS:

1. A computer-readable medium containing a transform operator constructed to provide an estimate of an unknown model with substantially optimal linear resolution.

2. A computer-readable medium according to Claim 1, wherein the transform operator comprises a matrix having at least one row of coefficients corresponding to a resolution function.

3. A computer-readable medium according to Claim 2, wherein the resolution function has substantially no negative values.

4. A computer-readable medium according to Claim 2, wherein the resolution function has a peak value of at least substantially 1 substantially at the point of interest.

5. A computer-readable medium according to Claim 2, wherein the resolution function monotonically decreases from its peak value.

6. A computer-readable medium according to Claim 2, wherein the resolution function has substantially no negative values, has a peak value of at least substantially 1 substantially at the point of interest, and decreases substantially monotonically from its peak value.

7. A computer-readable medium according to Claim 6, wherein the resolution function complies substantially with the following noise gain constraint:

8. A computer-readable medium according to Claim 2, wherein the matrix has a plurality of rows of coefficients, each row corresponding to a different resolution function.

9. A computer-readable medium according to Claim 8, wherein the matrix is constructed for use in multiexponential decay signal processing and wherein each resolution function is substantially centered on a different time constant.

10. A computer-readable medium according to Claim 9, wherein each coefficient has a corresponding data function in the form e-tke-v or derivations of this form, where y is the natural logarithm of a time constant and tk is a sampling point on the decay signal.

11. A computer-readable medium according to Claim 2, wherein the coefficients are selected so that the corresponding resolution function has a peak value of at least substantially 1 substantially at the point of interest and wherein the area of the resolution function is substantially minimized.

12. A computer-readable medium according to Claim 11, wherein the area of the resolution function is set to be substantially 1 on a logarithmic scale.

13. A method of constructing a transform operator in which a plurality of coefficients are selected to produce a corresponding resolution function that provides substantially optimal linear resolution.

14. A method according to Claim 13, wherein the resolution function has substantially no negative values, has a peak value of at least substantially 1 substantially at the point of interest, decreases substantially monotonically from the peak value, and has a gain within a predetermined range.

15. A method according to Claim 13, wherein the area of the resolution function is set to be substantially 1 on a logarithmic scale.

16. A method according to Claim 13, wherein the transform operator is constructed to provide substantially optimal linear resolution between data and an unknown model without any consideration of whether the data fits the unknown model.

17. A method of multiexponential signal processing, which comprises:
sampling a multiexponential signal and applying the transform operator of Claim 1 to the sampled signal.

18. Apparatus for multiexponential signal processing, which comprises a signal processor that has the transform operator of Claim 1.

19. A method of constructing a transform operator, in which a plurality of coefficients are calculated using data functions that emphasize linear resolution irrespective of data.

20. A method according to Claim 19, wherein each data function is in the form e -tke-v or derivations of this form, where y is the natural logarithm of a time constant and tk is a sampling point on a decay signal.

21. A method of exponential signal processing, which comprises:
providing a sampled multiexponential signal; and applying the sampled multiexponential signal to a transform operator constructed to provide substantially optimal linear resolution between outputs of the transform operator and an unknown model.