Spectral analysis using linear prediction z-transform and autoregression

Spectral analysis using linear prediction z-transform and autoregression

Volume 13 1. number 3 7 November I986 CHEMICAL PHYSICS LETTERS SPECTRAL ANALYSIS USING LINEAR PREDICTION z-TRANSFORM AND AUTOREGRESSION J. TANG a...

375KB Sizes 0 Downloads 37 Views

Volume 13 1. number 3

7 November I986

CHEMICAL PHYSICS LETTERS

SPECTRAL ANALYSIS USING LINEAR PREDICTION z-TRANSFORM

AND AUTOREGRESSION

J. TANG a and J.R. NORRIS ab a Chemistry Dwrsion, Argonne National Laboratory, Argonne, IL 60439, USA b Department of Chemistry, University of Chicago, Chicago, IL 60637, USA Received 11 July 1986; in final form 12 August 1986

LPZ AR, an alternative to fast Fourier transform spectral analysis, is presented with enhanced resolution and sensitivity. By combining the linear prediction z-transform method, and the more efficient Levinson-Durbin algorithm or the Burg algorithm, one can improve computation speed significantly as compared to the Householder or the singular value decomposition. Unlike the maximum entropy method, the LPZ AR method yields phase information and thus can be applied to analyzing signals with arbitrary phase.

Fast Fourier transformation (FFT) is the most commonly used form of spectral analysis of discretely sampled signals. However, limitations in FFT occur when data are truncated, for example when data observation time is limited by spectrometer inadequacies. Also for noisy or poorly resolved signals, FFT apodization routines can improve results but are subject to a trade-off between resolution and sensitivity [l-2]. Consequently, the maximum entropy method (MEM) [3-lo], by maximizing the Jaynes configuration entropy, has been employed in an attempt to overcome some of the deficiencies of FFT limitations. However, the current MEM cannot be applied to signals with arbitrary phase. As an alternative we recently presented the linear prediction z-transform (LPZ) method [11,12] of spectral analysis with improved resolution and sensitivity while retaining phase information. A discrete signal of ftite duration can be extended by linear prediction to an infinite data sequence and consequently truncation artifacts are avoided. The LP coefficients can be determined by linear least-squares fit using the Householder triangularization decomposition (QRD) [ 1 l- 151 or the singular value decomposition (SVD) [16,17]. However, the requirements or large computer memory and extensive computational time impose a limitation on the data size of about 1K. To 252

overcome these restrictions we present here a different approach, LPZ AR, by combining the LPZ method and autoregression method (AR) [18]. With LPZ AR the computational time can be shortened by about two orders of magnitude for 1K data as well as process a much larger data set (32K) such as a high resolution NMR FID. After an early beginning by Prony in 1795 [ 191, linear prediction (LP) theory has not been applied extensively until just recently. Using the linear prediction filtering method [ 181, for a signal containing M damped sinusoids each data point yn can be approximated by a linear combination of M previous data points (known as forward prediction) as M

Y,,=-~?~b,y,_,,

n=M,M+

l,..., N-

1,

(1) where N is the number of data points, the b, are the forward LP coefficients and M is the filter length. Instead of using QRD or SVD to solve the LP coeffiCientsasinLPZQRD [11,12],onemayusethe more efficient Levinson-Durbin algorithm or the Burg algorithm commonly used in the autoregression method [ 181 for a stationary process. Using the forward LP relation in eq. (1) and the ad0 009-2614/86/$ 03.50 0 Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

Volume I3

CHEMICAL PHYSICS LETTERS

1,number3

ditional assumption of a stationaryprocess, a similar LP relation known as the Yule-Walker equation [ 181 exists for the auto-correlation function as given by

M R(k)=--c b,R(k-n), for k>O, ?l=l

M R(O)=-nzl lqwo+

u2,

(2)

7 November 1986

Since only a magnitude spectrum is calculated using eq. (S), phase information for any individual peak is unknown and thus the MEM method is not eompletely satisfactory for many spectroscopies. To retain phase information, instead of the MEM, we proposed the LPZ method by extending unobserved data by LP [I 11. For data sequence yn of finite length, the finite discrete Fourier transform is given by

where a2 is the noise power. R(m), the auto-correlation function with lag I)Z,is defined by N-m-l

(6)

In the stationary process, the auto-correlation function depends only on the relative time lag, but not on the absolute initial time. The matrix equation (2) with Toeplitz form [18] can be solved efficiently with the recursive Levinson-Durbin algorithm [ 181. The algorithm is initialized by

where z = exp(--io7) and lf~ is the data sampling rate. Using the LP relation, a f&te data record can be extended to an infinite sequence to avoid truncation errors. The infinite Fourier series can be reduced to a finite summation with a closed form using LPZ [I I 1. The corresponding LPZ spectral function S(z) given below is equivalent to the result derived previously [11,12], except that forward prediction relation is used here instead of backward prediction,

b,, = -Nl)/R(O),

S(z) = H(z)IG(z),

4 = (1 - I b,, 12W~,

G(z)=

R(m) =

j$

ncO Yn+mYn*. =

(3)

m$o b,z-, I

and with the recursion for k = 2,3, .... M given by

M-l H(z) =

22 c,z-rn,

m=O

bki = b,_,,, + b&‘&-i,

i
u$=(1 - lbkkf2)@;-1.

(4)

blMl,blM2, .... biMM(or defined later as b,, b2, .... b,) are the final LP coefficients for falter length M. The best choice of the filter length M is not generally known a priori. Too low a guess for M results in a broader and highly smooth spectrum. Too high a filter length yields less improvement in se~itivity. Several criteria can be used for a less subjective selection of the AR model order [ 181. Once the LP coefficients are determined, the magnitude spectrum can be calculated using the Burg rn~~ entropy method [I813 -1

iS(o)I a 1 + a61 b, expfinwr) I

. I

(5)

where cm =Z,m=ob,ym_,,andbo~l. A fitted FID can be used in generating cm of the above equation for better sensitivity improvement. One may use the following iterative LP relation to estimate the fitted FID,

M Y*=-gl&;Yn+m,

n=-l,-2,...,-M+1

(8)

n=O, I,..., M- 1.

(9)

and &4 y,, =-f$lbmyn_m,

In eq. (8), one first caiculates the signal at negative time by backward prediction. For a stationary process, the backward LP coefficients are simply complex conjugate of the forward LP coefficients [ 181. The fitted 253

Volume 131, number 3

FID for the first M points can then be calculated using the forward prediction relation in eq. (9). There is one drawback of the LPZ AR as compared to the LPZ QRD, i.e. the assumption of a stationary process, which may not be a good approximation in some cases. For a highly non-stationary process, using eqs. (8) and (9) to generate the fitted FID may cause spectral distortion. Weaker peaks tend to be lower in intensity than the actual values. In this case, one may simply use the raw FID to generate the c, in eq. (7). The resultant spectrum may be noisier than using the fitted FID, but less distorted in intensity and shape. For better results, one may first use the routines in eqs. (8) and (9) to estimate the FID and its residues of fitting. Despite the large content of noise, the residues generally contain perceptible signals ignored by the first estimation. A second AR procedure can then be applied to recover the signals. In most cases, reasonably good results can be achieved in less than five cycles of the above procedures. Quantitative results of several testing cases using LPZ QRD, LPZ AR and FFT will be presented elsewhere [20]. Although the spectral function in eq. (7) is defined by a continuous variable z, in practice, a finite and discrete spectrum is calculated. The summation in H(z) and G(z) can be carried out using the FFT algorithm. The z transform [2 l] is a combination of discrete Fourier transform and discrete Laplace transform theory, where z = exp(-iwr - /CT).The ordinary Fourier spectrum is a special case of the z transform with k = 0. A two-dimensional LPZ spectrum can be constructed as a function of w and k, i.e. of the frequency and the decay rate. Spectra with reduced linewidth and thus enhanced resolution can be obtained by slicing through the 2-D spectrum [12]. This procedure has some similarities to the 2-D reconstruction using MEM reported by Sibisi [5], and is superior to apodization using a positive exponential function prior to FFT spectral analysis, since such apodization procedure in FFT can produce noise amplification and truncation artifacts while these problems are minimized by LPZ. There are some similarities between the LPZ AR and the Burg MEM method [ 18,221 in which the Levinson-Durbin algorithm was also used and modified. The conventional MEM method using eq. (5) only yields a power spectrum, and consequently phase information for each peak is lost. However, 254

7 November 1986

CHEMICAL PHYSICS LETTERS

I,

A

I

-500

8

-300

I

-100

1

I

100

300

0

500

HZ

Fig. 1. Frequencydomain spectra of simulated data with 13 damped sinusoidal signals with various intensity, frequency, phase, and Iinewidth. (A) 1K FFT spectrum without noise. (B) 1K FFT spectrum with noise. (C) 1K LPZ AR spectrum using the same noisy data as in (B) (in this example a filter 1engthM of 200 was used). The spectrum obtained by LPZ AR has better sensitivity than FFT.

using the LPZ spectral function defined in eq. (7), one may construct a spectrum with phase as in the conventional FFT spectrum. A few spectra, shown in figs. 1-3, illustrate the advantages of the LPZ AR over FFT or MEM. The LPZ AR method is particularly useful for truncated signals and in many cases yields an accurate spectrum with minimal artifacts, whereas the FFT spectrum is subject to degraded resolution and distortion in phase and amplitude. Applications of the LPZ AR method to 2-D NMR spectroscopy are particularly important because the 2-D signals are commonly subject to data truncation and in many cases have low sensitivity. In conclusion, we have presented the LPZ AR method, which exhibits better resolution and sensi-

CHEMICAL PHYSICS LETTERS

Volume 13 1, number 3

7 November 1986

tivity than FFT and has computational advantages over LPZ QRD, LPSVD, and LPQRD. Unlike the MEM method, the LPZ AR can be applied to signals with arbitrary phase. Thus, LPZ AR is a practical and powerfool tool for spectral analysis in magnetic resonance, FTIR, FT mass spectroscopies, as well as realtime speech analysis and coding. The authors thank Drs. R.E. Botto and P.H. Neil1 for providing us NMR data for the analysis. This work was supported by the Division of Chemical Sciences, Office of Basic Energy Sciences of the US Department of Energy under contract W-31.109-Eng38.

References I

-500

I

-300

1

-100

I

I

1

100

300

500

Hz Fig. 2. Frequencydomain magnitude spectra using the same data as in fii. 1B. (A) LPZ AR magnitude spectrum. (B) MEM magnitude spectrum with the same filter length. LPZ AR method yields better sensitivity than MEM. In addition, LPZ AR can yield phase information, which is not available using MEM.

A A

B

I

I

0

220

I

440

HZ Fig. 3. (A) A part of proton NMR spectra of cellobiose obtained by conventional FFT. The whole spectrum is calculated using 2K complex FFT of 2K raw data without zero-filling. (B) The same part of the spectrum obtained by 16K complex FFT with zero-padding of the data from 2K to 16K. (C) The same part of the spectrum using 16K complex FFT with apodiaation by a cosine waveform prior to zero-filling. (D) The LPZ AR spectrum of the same 2K raw data with filter 1engthMof 900. The LPZ spectrum avoids the truncation artifacts which appear in (B) and yields a spectrum with better resolution and less noise than the spectra in (A) and CC),

111 R.R. Ernst, Advan. Msgn. Reson. 2 (1966) 1. 121 J.C. Lindon and A.G. Ferrige, Progr. Nucl. Magn. Reson. Spectry. 14 (1980) 27. 131 S.F. Gull and G.J. Daniell, Nature 272 (1978) 686. [41 D.G. Collins, Nature 298 (1982) 49. PI S. Sibisi, Nature 301 (1983) 134. 161 S. Sibisi, J. Skillhig, R.G. Brereton, E.D. Laue and J. Staunton, Nature 311 (1984) 446. [71 E.D. Laue, J. Skilling, J. Staunton, S. Sibisi and R.G. Brereton, J. Magn. Reson. 62 (1985) 437. 181 P.J. Hore, J. Magn. Resort. 62 (1985) 561. PI E.D. Laue, J. Skill& and J. Staunton, J. Magn. Reson. 63 (1985) 418. WI J.F. Martin, J. Magn. Reson. 65 (1985) 291. 1111 J. Tang and J.R. Norris, J. Chem. Phys. 84 (1986) 5210. WI J. Tang and J.R. Norris, J. Magn. Reson. 69 (1986) 180. v31 J. Tang, C.P. Lin, M.K. Bowman and J.R. Norris, J. Magn. Resort. 62 (1985) 167. [ 141 C.L. Lawson and R.J. Hanson, Solving least squares problems (Prentice-Hall, Englewood Cliffs, 1975). [15] J.J. Dongarra, C.B. Moler, J.R. Bunch and G.W. Steward, LINPACK user’s guide (SIAM, Philadelphia, 1975) ch. 9. [16] R. Kumaresan and D.W. Tufts, IEEE Trans. ASSP-30 (1982) 833. [17] H. Barkhuijsen, J. de Beer, W.M.M.J. Bovee and D. van Ormondt, J. Magn. Reson. 61 (1985) 465. [18] S.M. Kay and S.L. Marple, Proc. IEEE 69 (1981) 1380. [19] G.R.B. Prony, J. L’Ecole Polytechnique 1 (1795) 22. [20] J. Tang and J.R. Norris, to be published. (2 1 ] A.V. Oppenheim and R.W. Schafer, Digital signal processing (Prentice-Hall, Englewood Cliffs, 1975). [22] J.P. Burg, Thesis, Stanford University (1975).

255