Signal-to-noise ratio enhancement in NMR spectroscopy using a vector-space technique

JOURNAL

OF MAGNETIC

RESONANCE

73,423-435

(1987)

Signal-to-Noise Ratio Enhancement in NMR Spectroscopy Using a Vector-Space Technique BEHZAD NOORBEHESHT,* HUA LEE,? AND DIETER R. ENZMANN* *Department of Diagnostic Radiology and Nuclear Medicine, Stanford University SchooI of Medicine, Stanford, California 94305, and tDepartrnent of Electrical and Computer Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801

Received September 22, 1986; revised February 3, I987 A new technique that improves the effective sensitivity of nuclear magnetic resonance by means of post-detection signal processing is presented. This technique will be of value in obtaining localized spectra, allowing a smaller voxel volume for more accurate localization. The technique takes advantage of the a priori knowledge of relaxation times and resonant frequencies of the spectral lines to be measured and uses vector-space methods to determine the corresponding concentrations. With use of this technique, as long as the number of samples is above a certain value, the extent of signal-to-noise ratio improvement is limited only by the spin-spin relaxation time. In particular, under the conditions which usually exist in a typical NMR spectroscopy experiment, this technique results in an eightfold increase in signal-to-noise ratio over that obtainable by conventional algorithms using Fourier transformation. Application of the technique to experimental data confirms the theoretical findings. o 1987 Academic FWS, IIIC. INTRODUCTION

Although NMR spectroscopy has been used for over three decades to characterize molecular structures, the technique has relied upon large sample volumes and long data acquisition times to overcome its low sensitivity and to achieve the signal-tonoise ratio (SNR) necessary for accurate quantification of spectra. It is highly desirable to be able to detect spectra from small volumes in a reasonable length of time and to present this information as a two-dimensional image at high resolution. Due to the low SNR that would result under the above imaging conditions, such an experiment is difficult at the present time. To achieve the two-dimensional image, an improvement in SNR of approximately two to three orders of magnitude over that presently achievable would be required. Factors that influence SNR are numerous and involve almost every aspect of the design of an NMR spectrometer (I). However, they may be divided into two broad categories: those which are primarily determined by hardware design and those which are principally software related. We will address the latter problem. A number of spectrum analysis techniques are available for application to NMR spectra (2). The most widely used technique in present NMR spectrometers is the Fourier transform technique first proposed by Ernst (3). Ernst also proposed apodizing the free induction decay to improve the quality of the detected data. Application of this technique has resulted in improved sensitivity and resolution; however, one usually 423

0022-2364187 $3.00 copyrightQ 1987 by Academic AU r&s

Pms, Inc. of reproduction in any form resend.

424

NOORBEHESHT,

LEE,

AND

ENZMANN

has to trade off resolution for sensitivity and vice versa. In addition, some apodizing functions result in line distortions. Another signal processing technique which has recently received considerable attention is the maximum entropy method, first developed for geophysical applications (4) and later applied to NMR spectroscopy (5). This technique does not require any a priori knowledge about the FID and the results are free of line-broadening artifacts. However, it has been shown that the maximum entropy method does not result in better SNR than that obtainable from the FT techniques (6). The linear prediction (LP) techniques have recently been proposed to overcome some of the drawbacks of the IT algorithms, such as the need for preprocessing the data and truncation artifacts (7, 8). However, LP methods do not result in any significant SNR improvement and indeed their use is limited to applications where SNR is high. In this paper we present a technique for enhancement of SNR by post-detection signal processing. The technique to be presented is a vector-space method (9) and results in improved SNR in an optimum fashion. In terms of SNR improvement, our results are equivalent to those obtained using the maximum-likelihood criterion (10, II). Both of these approaches obviate the need for Fourier transformation. We will first develop the necessary mathematical relationships and then illustrate the use of this technique by applying it to actual phosphorus-3 1 spectral data under varying SNR conditions. THEORY

For mathematical simplicity we will start our analysis using a continuous-time model and later will generalize it to include the discrete case. Let us assume that the signal, originating from a small voxel of volume V, is in the form of a free induction decay and is continuously observed during a time interval TO. Using complex notation, the signal due to the pth spectral line in V may be written as S(t,~p) =~pe-W~p-iWd~, O
VECTOR-SPACE

425

TECHNIQUE

Assuming volume I’ contains N spectral lines, the total signal may be represented bY O=st<

s(t, A) = 5 ApMO,

p=l

where

To,

4 P(t) = e-(Wp-i2~f,U

PI [31

and A=[AL,A2,.

. . ,AN]=.

141

The received data, r(t), consists of the above signal corrupted by random noise, n(t), r(t) = s(t, A) + n(t),

O
To.

[51

We will assume n(t) to be white and Gaussian with zero mean, variance u2, and power spectral density 7,~/2.Our goal is to determine A+, from a knowledge of r(t) from 0 to To, and the values of T,s and f,s. Addition of noise converts the unknown deterministic constants, A,, to random processes whose expectations and variances are Ap and 2, respectively. Our task then becomes the estimation of the expectations of these random processes. We will obtain these estimates by first defining a complex vector-space S whose elements are &p(t). We will then determine a vector A in S which is closest to the vector of actual complex amplitudes, A. By the projection theorem (9), the elements of A will be the optimal estimates of the elements of [A], whose magnitudes are the desired spectral line concentrations. We define the inner product of two complex functions, a(t) and b(t) as (4) (a(t),

b(t))

=

="a(t)b*(t)dt. s0

Applying the above operator to Eq. [5] gives (r(t), 4XQ) = (dt, Ah 4,(t)) + (40, &At)).

171

In Eq. [6] *denotes the complex conjugate operator. Substituting for s(t, A), we get (r(t), t&At)) = p+&#t&),

&At)) + (MO, 4&>).

[81

We assume that for long enough To, n(t) and & are poorly correlated. Therefore, the second term on the right-hand side of Eq. [S] may be neglected and that equation may be rewritten as (m44w)

q= 1,2,. . . ,N

=pi+@Jm#dt~)

where, because of the above approximation,

[91

we have replaced the actual amplitudes,

A,, by their estimates, a,. We now define the column vector R and matrix CDas R=B1,R2,.

where

. . ,&vIT,

4 = (r(t), &dO)

1101 1111

426

@'I1 a= :[.'**@Nl .:1

NOORBEHESHT,

and

LEE,

AND

1121

@NN

+lN

where

ENZMANN

-i2rf,)To %Q

= (dJp(O,

&7(O)

=

l ---;2*f

w

iI31

with am= 1/T,+

1/T,

f,=f*-f,.

t141

WI

Matrix @ is sometimes referred to as the Gram matrix. Equation [9] may now be written in matrix form as R=CD& [I61 where w is the vector of estimates, a,. Equation [ 163 is often called the normal Eq. [9]. The estimates may be determined by the solution of Eq. [ 161, i.e., A=@-‘R.

1171

The magnitude of &, gives the concentration of the pth line, whereas its phase is the phase of the @h FID. The above estimates are equivalent to those derived using the maximum-likelihood criterion (10, II) and it can be shown (II) that they are unbiased and efficient estimates of [A], that is, their variances are minimum. Therefore, they reach the Cramer-Rao lower bound. The variances of&s are given by the diagonal elements of the Fisher information matrix, W, (II) defined as ly =p-1. In deriving Eq. [ 181, we have taken advantage of the fact that the autocorrelation white noise is a delta function, i.e., g(7) = ; 6(T).

iI81 of

1191

In the above derivation no restrictions were placed on the relative values of Tp and fP, consequently this technique is applicable even if there is considerable overlap between the spectral peaks. However, in situations where there is only minimal overlap the above results may be simplified further because matrices W and CDmay be approximated by their diagonal elements. The diagonal elements of ly are then given by qpp=L PO1 2@PP substituting for aPPfrom Eq. [ 131, we get

Pll

VECTOR-SPACE

427

TECHNIQUE

In the above derivation we have assumed r(t) to be continuously detected. Because, in many applications, data are sampled discretely, it would be desirable to have relationships that are applicable to discrete data as well. Although, when the number of samples is large, Eqs. [ 131 and [2 l] would be approximately true for discrete data, exact relationships may be derived by transforming the integrals in those equations into summations and evaluating the resulting sums. For example, the discrete version of Eq. [21] may be derived as follows:

P24

%P = s0 Tol&,,(t)12d~g”i’ ]&,12AT n=O = AT

M-l 2

1 _ e-2MATlTp e-WTpWT)

= AT 1 _ e-2ATl

n=O

Tp

Wbl with To = (M-

l)AT,

~231 where M and AT represent the number of samples and the sampling interval, respectively. Assuming r(t) to be band-limited, we can rewrite Eq. [22b] as

where 2F is the noise bandwidth, and 2 is the noise power. The above results are valid only if the noise samples are statistically independent. This will be strictly true only if samples are taken at the zeros of the autocorrelation function of the noise. In our case, because we have assumed n(t) to be white and band-limited, the autocorrelation function is a sine function given by sin(2?rFT) g(7)= a2 2?rFT * Therefore, for samples to be independent, dition

the sampling interval must satisfy the con-

AT=’

2F’

WI

This condition not only assures independent sampling but also satisfies the Nyquist criterion for alias-free sampling. Substituting Eq. [26] into Eq. [24] we get

Equation (27) relates variance of each estimate to variance of the noise. Because the power of a random process is given by its variance, C? and qPP represent the noise powers before and after processing, respectively. Defining SNR as the ratio of initial amplitude of the FID to standard deviation of the noise, i.e.,

428

NOORBEHESHT,

LEE, AND ENZMANN

we may rewrite Eq. [27] in terms of the corresponding SNRs as (SNR), = (SNR), [ ll$z]li2,

1291

where (SNR), and (SNR), represent SNRs of the processed and unprocessed data, respectively. RESULTS

The expected improvement in SNR, as given by Eq. [29], may be determined by substituting for the various parameters in that equation values that are normally found in a typical NMR spectroscopy experiment. Figures 1 and 2 show the resulting plots of Eq. [29] as a function of TO and Tp. As expected, SNR improves as observation time is increased. However, because data decays exponentially with time, no significant further improvement in SNR is achieved as To is extended beyond two to three time constants. Signal-to-noise ratio also improves with increasing Tp. Again this is to be expected because, for a given To, the relative amplitude of each signal sample increases as Tp increases. The increase in SNR with T, asymptotically approaches its maximum as T, tends to infinity. The maximum gain in SNR for large Tp is proportional to the square root of the number of samples, M, because for large Tp the exponential decay is negligible and all samples contribute equally to the signal level. The range of values of T, in Fig. 2 has been chosen to cover most of the values that are likely to be encountered in P-3 1 spectroscopy.

6.0 2 0

1

2

3

4

To/T, FIG. 1. Expected enhancement in SNR as a function of observation time, T,, normalized by the spinspin relaxation time, TP. TP = 15 ms and AT = 0.125 ms.

VECTOR-SPACE

429

TECHNIQUE

16 -

11 -

12 -

0

20

40 T,

FIG. 2. Expected AT = 0.125 ms.

enhancement

in SNR

as a function

60

b-1

of spin-spin

relaxation

time,

Tn. To = 125 ms and

Equation [29] indicates that the extent of SNR improvement is a function of AT, for large M, such as is found in most spectroscopic experiments, the dependence on A4 (or TO) vanishes and the numerator in Eq. [29] approaches unity. It can be shown that if A4 satisfies the lower bound Tp, and A4 (or To). However,

ML

60000

[301

AT

3 F

0

2

I

~I~1III~IIl,,~,I,,,,r,~~~, 20 10 CHEMICAL

YV’ 0

-10 SHIFT

-20

(6)

-30

-40

ppm

FIG. 3. “P spectrum obtained by Dm at high SNR. 1 = j3AF, 2 = CUATP, 3 = TATP, (PCr); 5 = phosphodiesters (PD); 6 = inorganic phosphate (Pi); 7 = sugar phosphates.

4 = phosphocreatine

430

NOORBEHESHT,

LEE, AND ENZMANN TABLE 1

Experimental Parameters Magnetic field strength

Surface coil diameter

(T)

AT (ms)

km)

2.0

3.0

0.125

To (ms)

M

125

1000

the improvement in SNR is 93% of its maximum value (see Fig. 1). Therefore, if M is larger than the above minimum value, the improvement in SNR for each spectral line is governed only by the ratio of the sampling interval to the corresponding spinspin relaxation time. However, the value of AT is dictated by Eq. [26] and is not, in general, under experimental control. This leaves the spin-spin relaxation time as the only determinant of SNR improvement. To evaluate the above technique we applied it to experimental data. Figure 3 shows the frequency spectrum of the data we used. It represents the phosphorus-3 1 spectrum from a rabbit brain acquired at 2 T. Table 1 summarizes pertinent experimental parameters. There are eight peaks discernable in the spectrum, seven narrow peaks and one broad peak. The broad component is due to the phosphorus in bone and was neglected in our calculations. The other seven peaks are due to adenosinetriphosphate (ATP) (three peaks), phosphocreatine, sugar phosphates, inorganic phosphate, and phosphodiesters. Relaxation times ( Tp) and resonant frequencies (f,) for each spectral line were calculated and are given in Table 2. The data in Fig. 3 were obtained using a conventional fast Fourier transform algorithm and represent the average of 256 acquisitions resulting in very good SNR. In our experiments we neglected the small amount of noise present in the data and assumed it to consist of signal only. We then used it as a reference to determine the efficiency of our noise-reduction algorithm. Figures 4-9 show the results of progressively decreasing SNR and the corresponding processed data. Signal-to-noise ratio was reduced by addition of computer-generated Gaussian white noise to the original data.

TABLE 2 Known Signal Parameters

Spectral line

Relaxation time (ms)

Chemical shift (ppm) (relative to PCr)

/3ATP olATP yATP Phosphocreatine (PCr) Phosphodiesters (PD) Inorganic phosphate (Pi) Sugar phosphates (SP)

8.0 1.4 10.8 15.4 5.5 14.9 13.0

18.1 8.5 2.1 0.0 -3.1 -5.8 -7.3

VECTOR-SPACE

431

TECHNIQUE

c 60000 z 2 J 2

40000

2 E x %I

20000

Y F: 2 z

0

20

0 -10 10 CHEMICAL SHIFT (6)

-20 ppm

-30

-40

FIG. 4. Spectrum of noisy data, SNR = 10 (20 dB).

Figure 4 shows the resulting spectrum when SNR is 10 (for the smallest peak). It can be seen that compared to the original noiseless data the spectrum is somewhat degraded, while the corresponding processed spectrum (Fig. 5) is essentially equivalent to the original data. As SNR is further decreased to 5 (Fig. 6), the spectrum is degraded further. The heights and widths of the true peaks now deviate substantially from their actual values. The corresponding processed spectrum, shown in Fig. 7, still quite closely resembles the superimposed original data. Finally, as SNR is reduced to 2 (Fig. 8) the original seven peaks are almost completely buried in noise but there is still reasonably good agreement between the processed and original data (Fig. 9).

60000

L

20

10 0 -10 CHEMICAL SHIFT (6)

-20 ppm

-30

-40

FIG. 5. Spectrum of noisy data after processing. Solid curve = processed data, solid with circles = original data. SNR = 10 (20 dB).

432

NOORBEHESHT,

LEE, AND ENZMANN

60000

n

- 20

10

0

CHEMICAL

-10 SHIFT

-20 (6)

-30

-40

pprn

FIG. 6. Spectrum of noisy data, SNR = 5 (14 dB).

It should be noted that due to the proximity of the peaks in the above spectra, matrix W was not diagonal; therefore, calculations of the estimates involved a matrix inversion operation as indicated by Eq. [ 171. However, in situations where the peaks do not overlap significantly nondiagonal elements of VI are negligible and its diagonal elements are given by Eq. [27]; therefore, no matrix inversion operation would be involved in calculation of the estimates. Also, as pointed out earlier, A,s are complex quantities whose magnitudes, the areas under the corresponding peaks, are usually the unknowns of interest. There is, therefore, no need to obtain the spectrum explicitly. Consequently no Fourier transform algorithms are necessary.

60000

0

20

10 CHEMICAL

0

-10 SHIFT

-20 (6)

-30

-40

ppm

FIG. 7. Spectrum of noisy data after processing. Solid curve = processed data, solid with circles = original data. SNR = 5 (14 dB).

VECTOR-SPACE

433

TECHNIQUE

100000 E k

60000

3 4

60000

p: G E

40000

Y F:

20000

2 2

n 20

10

0 CHEMICAL

FIG. 8. Spectrum

-10 SHIFT

of noisy

-20 (a)

data, SNR

-30

ppm

= 2 (6 dB).

DISCUSSION

The technique presented here is based on vector-space optimization and results in an SNR improvement equivalent to that obtainable from maximum-likelihood estimation technique (10). Although the final results of these two techniques are similar, the approaches are different. For mathematical simplicity the necessary relationships were derived using continuous-time representation and were extended to the discrete case by simple mathematical manipulations of the final result. Due to the simpler form of the resulting mathematical relationships, this approach allows perhaps better

60000

CHEMICAL

FIG. 9. Spectrum of noisy data. SNR = 2 (6 dB).

data after

processing.

SHIFT

Solid curve

(a)

ppm

= processed

data, solid with circles

= original

434

NOORBEHESHT,

LEE,

AND

ENZMANN

understanding by the nonspecialist and as a result makes its implementation easier. Taking advantage of the functional form of the elements of the Gram matrix, CD, closed form expressions for its elements were derived which made the technique computationally more efficient. By using actual values for resonant frequencies and relaxation times, we were able to determine the extent of SNR improvement achievable by our technique. Under conditions which usually exist in most NMR spectroscopy experiments, the extent of SNR improvement for each spectral line was found to be determined solely by the spin-spin relaxation time, either machine-dependent or object-dependent, whichever is shorter. Thus, the longer the relaxation time the greater the improvement in SNR. For Tp’s in the range of lo-20 ms, the resulting improvement in SNR was found to be on the order of a factor of 8. This improvement in SNR can be translated into either shorter experiment time or smaller voxel size. The former is important in NMR spectroscopy where shorter experiment time allows study of materials which are otherwise too unstable for longer data acquisition times. The latter would have application for in viva spectroscopy with smaller volumes. In conventional in vivo spectroscopy it allows detection of low concentration lines. This technique has application in two-dimensional spectroscopic imaging by affording better spatial resolution. In its present form this technique is directly applicable only to situations where the relaxation times and resonant frequencies are known a priori. Such conditions may exist, for example, in imaging applications where one is interested in measuring relative concentrations of only a few lines with well-known chemical shifts in a controlled and stable magnetic environment where the relaxation times are fairly constant. Another favorable situation may arise in studies where one has to run a large number of experiments on the same sample to study the relative changes in the concentration of various lines over a period of time (for example, before and after therapeutic intervention in an in viva study). It is emphasized that while absolute resonant frequencies may not be constant and may shift with perturbations in various chemical, physiological, and electromagnetic parameters, their relative positions along the frequency axis remain more or less constant. Therefore, using the signal from a strong line (e.g., water), their absolute values may be obtained (12). It has been suggested (13) that the minimum SNR for acceptable image quality in clinical applications of MRI is on the order of 20. Therefore, application of this technique when used in conjunction with an appropriate spatial localization algorithm would result in clinically useful images even when the SNR is on the order of 3. However, much lower SNR than what is given above may be tolerated when the spectra are presented two-dimensionally in an image format (5). The reason for this is that with such a presentation a priori knowledge of anatomy and physiology can be used by the observer to offset the effect of apparently poor SNR to extract useful clinical information. It should be added that the inherent spatial integrating ability of the human eye further improves a two-dimensional image’s effective SNR. We have not addressed the problem of spatial localization in a two-dimensional spectral imaging system. However, with appropriate modification of the results presented in this paper, similar SNR improvement could be expected when this technique is applied to NMR signals containing spatial information.

VECTOR-SPACE

TECHNIQUE

435

ACKNOWLEDGMENTS The authors thank General Electric Company’s Fremont, California, division and Dr. Mike Moseley of University of California at San Francisco for kindly providing the raw experimental spectral data. REFERENCES 1.

2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.

D. SHAW, “Fourier Transform NMR Spectroscopy,” Elsevier Scientific, New York, 1976. S. M. KAY AND S. L. MARPLE, Proc. IEEE 69, 1380 (1981). R. R. ERNST, in “Advances in Magnetic Resonance” Vol. 2, Academic Press, New York, 1966. J. P. BURG, “Maximum Entropy Spectral Analysis,” Thesis, Stanford University, 1975. S. SIBISI,Nature (London) 301, 134 (1983). J. F. MARTIN, J. Magn. Reson. 65,291 (1985). H. BARKUSIJSEN,R. DE BEER, W. M. M. J. BOVEE, AND D. VAN ORMONDT, J. Magn. Reson. 61,465 (1985). J. TANG, C. P. LIN, M. K. BOWMAN, AND J. R. NORRIS, J. Magn. Reson. 62, 167 (1985). D. G. LUENBURGER, “Optimization By Vector Space Methods,” Wiley, New York, 1969. A. MACOVSKI AND D. SPIELMAN, Magn. Reson. Med. 3,97 (1986). A. D. WHALEN, “Detection of Signals in Noise,” Academic Press, New York, 197 1. A. MACOVSKI, Department of Radiology, Stanford University, private communication. W. A. EDELSTEIN, “Society of Magnetic Resonance in Medicine, Third Annual Meeting, New York, August 13-17, 1984.”