International Journal of Infrared and Millimeter Waves, Vol. 10, No. 2, 1989

THE MARPLE ALGORITHM FOR THE AUTOREGRESSIVE SPECTRAL ESTIMATES OF THE SMMW FOURIER TRANSFORM SPECTROSCOPY DATA* Zhang Guangzhao and Zhou Guangqun Electronics Department, Zhongshan University Guangzhou People's Republic of China Received November 18, 1988

The Marple algorithm for the autoregressive spectral estimates has been applied to the SMMW Fourier transform spectrum analysis. The experimental results have shown that this method yields AR spectra with three times higher resolution than the FFT method does. The improvements obtained from the Marple algorithm over the maximum entropy algorithm include higher resolution, less bias in the spectral peak frequency estimation and absence of observed spectral line splitting. The effects of the structure of the spectral lines and the noise on the resolution are discussed. Key words: Marple algorithm, autoregressive spectral estimates, Fourier transform spectroscopy. Introduction One of the advantages of the Fourier transform spectroscopy over the conventional spectroscopy is that it has ability to achieve higher resolution. However, the resolution of Fourier transform spectrometers is still restricted by signal to noise ratio, scanning distance of moving mirror and the apodization functions used. In order to get higher *Projects supported by National Natural Science Foundation of China 257 0195-9271/89/0200-0257506.00/09 1989 Plenum Publishing Corporation

258

Zhang and Zhou

resolution, the signal to noise ratio has to be improved by using more sensitive detectors and the scanning distance to be extended. Apparently, it would make the spectrometer more expensive. Theoretically, the interferogram produced by a Michelson interferometer is extended towards doubleside from the zero path difference point infinitely. Buts in fact, the scanning distance of the moving mirror can not be infinite, only a finite-length interferogram can be obtained. The data which are outside the measured interferogram are supposed to be zero when taking Fourier transform. This is equivalent to apply a "window" on the infinite interferogram and it decreases the resolution. If we can adapt a mathematical model to predict the unknown information which are outside the sampling region from the sampling data, and if the model is resonable, it would eliminate the "window" effect and achieve higher resolution. S.M.Key (I) has discussed several common mathematical models, in which the most popular one is autoregressive (AR) model. Among all autoregressive algorithms, the maximum entropy method (Burg algorithm) (2), forward and backward least-square method (LS algorithm) (3,4) and singular value decomposition method (SVD algorithm) (5) are usually used. In 1983, S.Kawata et al (6) applied the Burg algorithm to analyze middleinfrared Fourier transform spectroscopy data first. In 1985, they (7) did the same work by using the SVD algorithm and make the spectrum estimates more reliable. Burg algorithm reduces the numbe# of computational operations to be proportional to the square of AR order p by introducing Levinson recursion constraint condition. It saves computation time, but the further improvement of resolution is restricted and the bias in the spectral peak frequency estimation and the spectral line splitting are caused. The SVD algorithm eliminates the Levinson recursion constraint condition, so the resolution can be further improved, the line bias and splitting are reduced. Because of absence of recursive algorithm, the number of computational operations is proportional to the cube of p. In 1985, Marple (8) introduced a new algorithm for LS method, which not only eliminates the Levinson recursion constraint condition, but also keeps the number of computational operations proportional to square of p. Therefore,

SMMW Fourier Transform Spectroscopy

259

Marple's LS algorithm can give fast computational speed, higher resolution, less line bias and absence of observed spectral line splitting. In this paper, the Marple algorithm is applied to the SMMW Fourier transform spectrum analysis. The experimental results have shown that this method yields AR spectra with three tines higher resolution than the FFT does. The effects of the structure of the spectral lines and the noise on the resolution and prediction accuracy of the spectral lines are discussed. Marple's LS al~orithm Essentially, AR spectrum analysis is a linear prediction method. Let us consider a data sequence (Xo, Xl, ...,xN_1) , in which each sample can be predicted with p past samples by ~n = - ~

apk Xn_ k ,

n=p, p+1, ..., N-1

(1)

where ~

is prediction value of x , a is AR paran n pk meter and p is AR order. The Eq.(1) is also called forward prediction model, in which each current output sample x n is only relative to past output "

samples. From Eq.(1), are given by

the forward prediction errors

fP n=xn-~n= ~=o ~ apk Xn-k '

(2)

n=p, p+l, ..., N-l, where apo=l. Similarly, the backward prediction errors are given by p

bpn=~o

a*pk Xn-p+k '

(3)

n=p, p+1, ..., N-I. According to the Eqs.(2) and (3), the sum of the forward and backward prediction error energy is given by e

=

fpn 2 +

~

bp n

(4)

In the forward and backward LS algorithm, the total prediction error energy is required to be minimized to all the AR parameters, that is

260

Zhang and Zhou

ep api

2~ a rp(i,J) j=O pj

O

(5)

i= 1,2, ... , p where rp(i,j) -

N-p-1 ~ (x

-

~+P-Jx*k+p-i

k=O

+

x~+j ~+i)

(6)

By Eqs.(2)-(5), the minimum prediction error energy is found to be P e ~ a . rp(O,j) (7) pmin = j=O P0 Eqs.(5)-(7) can be combined into a single (p+l)(p+l) matrix equation: Rp Ap = Zp where

(8)

rp(O,O) .... rp(O,p)~ @

Rp =

" rp i p,O) .... rp(p,p) " epminl

(9)

0

Ap =

,

Ep

@

= @

app Solving the matrix equation (8),

AR parameters apk

and minimum prediction error energy epmin can be obtained. by P(f)

The power spectral density P(f) is given

~2 At =

............

P

.....

(10)

2

where ~ p = epmin/N is called the minimum prediction error power, f is frequency and 4t is the autooorelation lag interval in seconds. Solving Eq.(8) by Gaussianmethod requires on the order of the cube of p operations. Marple found

SMMW Fourier Transform Spectroscopy

261

that the matrix Rp, though, does not have the Toeplitz symmetry, it has structure composed of the sum of two Toeplitz data matrix products. This structure allows a recursive algorithm to solve matrix equation (8) and the number of computational operations is proportional to the square of p. For the detailed algorithm, computational flowchart and the FORTRAN program, the readers are recommended to reference (8). Experimental results We have applied the Marple LS spectral estimation to the interferograms obtained by SMMW Fourier transform spectrometer. In our previous works (9-ql), the procedure used to realize superresolution of an absorption spectrum is two-step spectrum analysis. That is, the measured background and sample interferograms are Fourier transformed individually, the absorptance spectrum is formed by dividing sample spectrum by background spectrum, the spectral region of interest is truncated from the spectrum, and the background-free sample interferogram is reformed by the inverse Fourier transform, then, the AR spectrum analysis is applied to this interferogram. This procedure, essentially, is divided into two steps: first, produce the background-free sample interferogram of the interesting spectral region from the measured interferograms by Fourier transform and inverse transform, and then, apply the AR spectrum analysis to the background-free sample interferogram~ The advantage of this procedure is that only a software package is needed without any change in experimental set up. However, the Fourier transform and inverse transform may cause information loss and distortion, which can effect the resolution and prediction accuracy of spectral line positioning. In order to overcome the drawbacks of the twostep spectrum analysis procedure, we introduce a new method to handle experimental interferograms. We measure the background and sample interferograms with optical filter respectively, them subtract sample interferogram from background interferogram to produce the background-free sample interferogram, finally, apply Marple LS algorithm to the backgroundfree sample interferogram. This procedure is shown in Fig.1.

262

Zhang and Zhou

I

M~aaure sample interferogramIs l and backgroundinterferogramIb

1

j.......Background-freesample interferogram=I~-I~

1

Calculate the power spectr~ d e n s i t y P(~) by AR spectral estimates

1

""'Theabsorptio~spectrum

J

s(~)=Cp(~))

I

J

i]

J

J

Fig.l. Experiment.-and computational flowchart of LS spectral estimates Only one step spectrum analysis is used. It can reduce the information loss and distortion. Therefore, it can improve the resolution and prediction accuracy of spectral llne positioning. Particularly, this method is sutable for the measurements of gas samples. We have measured the 14NH 3 gas ground-state

in-

verse spectrum with a far-infrared Fourier transform spectrometer. The sample cell is 5cm long with 9mm total thick crystalline quartz windows as an optical filter. The beamsplitter is qO~m thick mylar. The working spectral region of this spectrometer is 50130cm -I. The gas pressure is about 30 torr. For convenience of comparing the LS spectrum estimation with the Fourier transform spectrum, we measure 2048 points double-side background and sample interferograms (sampling interval is qo.q248~m) and only took 600 points around the zero path difference point for the LS spectrum analysis. The results are shown in Fig.2. Fig.2(a) and (b) are 600 points background and sample interferograms, respectively. Fig.2(c) is the background-free sample interferogram generated from background interferogram minus sample interferogram. Fig.2(d) is the LS prediction spectrum with

SMMW Fourier Transform Spectroscopy

263

-~I I~---'I i'-~--'-! ['~-~ ('~

e9.~5

i,i,,,J 1/

!1

(e)

i84

I~9,719

Q

,;~, ?,3

(o)

optical path difference

(r

48.Z3 Wo'r

~,81 (ci~ -1 )

1Z~.39

Fig.2. 14NH_ FFT spectrum and LS spectral estimates N=6OO and p=248~ Fig.2(e) and (f) are 512 points and 2048 points FFT spectra, respectively. Comparing Fig.2(d) with Figo2(f), we found that the 600 points LS prediction spectrum is comparable with the 2048 points FFT spectrum. That is, the resolution of the LS spectral estimates is about three times higher than that of the FFT method. Table I shows the LS prediction values with N=6OO and p=248 and the FFT values with N=2048 of the 14NH 3 spectral line wavenumbers.

For convenience of

comparison, Tab.1 also shows the calculational values which are average results from Ref~176 From Tab.1 1 we found that the deviation of the LS prediction values are very small.

Zhang and Zhou

264

Tab.1. The comparison of the calculational values of the 14NH 3 spectral line wavenumbers with the LS spectral estimates and the FFT results No. of

C~,

LS Speotral e s t i - ...... FFT method

spectra] value= mates(N=660~p=248) lines from predic%ed deviat ibn

1 2 3 4 5 6 7 8

Ref.|J values 58.83 58.55 60.35 60.30 78.65 78.73 80.12 8o.13 98.40 98.54 99.86 99.84

118.11 t18.12 119.47119.51

0.28 0.05 0.08 O.Ol

o.14 0.02

o.01 0.04

I~=2048 cal. values 58.70 60.33 78.59 79.96 98.31 99.88 118.1o 119.1o

deviation's"

o,13 0.02 0,06 o,16 o,o9 o.o2 O,Ol o.37

Discussions We shall discuss the effects of the structure of the spectral lines and the noise on resolution and the prediction accuracy of the spectral line positioning with simulated interferograms given by Yn=

k - (n-N/2)2 ~_~ e J=l +en ,

x2 (11)

where N is the number of sampling points, A x is the sampling interval, aj and O"~ are the line-width and the wavenumber of the jth GaussianJ spectral line, respectively, and e is the Gaussian white noise, n

By changing noise level, number and interval of the spectral lines, we can study their effects on the resolution and the prediction accuracy. (I) The effects of the noise level and the structure of the spectral lines on the resolution From Eq.(11), supposing N=128, ~x=O.OO4cm and a~=O, the limitations of resolution by LS algorithm with different noise levels and different structures of spectral lines can be calculated and the results are shown in Fig.3. We can see that the LS algorithm has the highest resolution for the spectrum which consists of two lines only, and the resolution is decreased with the increasing of number of spectral lines. Furthermore,

265

SMMW Fourier Transform Spectroscopy

desh lines: the interval between groups is 10 resolution units

~1.5 m

o

1.0

/

/ /

/

5 6 /

/

2

/

.~" ~roups

~ o.~ B

2 r~l~tive

3' noise

levsl

4 (%)

Fig.3. The effects of noise levels and spectral line constructions on the resolution the resolution for the spectrum which consists of several close equi-interval lines (solid lines in Fig.3) is decreased faster than that for the spectrum which consists of several double lines (desh lines in Fig.3). Therefore, in order to get higher resolution, we should try to reduce the number of spectral lines by using optical or electronical filters. From Fig.3, we can also see that noise has great influence upon the resolution. (2) The effects of the noise level and the interval of the spectral lines on the prediction accuracy of spectral line positioning From Eq.(S1), supposing N=128, ~x=O.OO4cm, aj=O , k=3 and relative noise level equal to 1%, calculate the prediction deviations of the spectral line positioning with different spectral line intervals, the results are shown in Fig.4. Supposing the spectral line interval equal to 1 resolution unit, calculate the prediction deviations of the spectral line positioning with different

266

Zhang and Zhou

O.OG

.~

0.05

0.04

0.03

0.02

0.01

0 0.6

0.7

0.8

0.9

1~

interval of z~cc,~s~ive ~!ectral lines

Fig.4. The effects of spectral line intervals o n prediction deviations of spectral line positioning relative noise levels from Eq.(11), the results are shown in Fig.5.

0.06

.~ 0.05

0.04

:

0.03

0.02

0.01

rolative noise lev~l (%)

Fig.5. The effects of noise levels on prediction deviations of spectral line positioning

SMMW Fourier Transform Spectroscopy

267

From Figs.4 and 5, we can see that the larger the spectral line interval is, the smaller the prediction deviation is; and the smaller the noise is, also the smaller the prediction deviation is. Conclusions Marple LS algorithm for the AR spectral estimates of the SMMW Fourier transform spectrum can improve the resolution effectively. The experimental results with 14NH 3 gas spectrum have shown that the resolution of the LS spectral estimates is about three times higher than that of the FFT method. Furthermore, the simpler the structure of the spectral lines and the lower the noise levels are, the higher the resolution and the smaller the prediction deviations are. References I. S.M.Kay and S.L.Marple, Proc. IEEE, 6__9, (1981) 1380 2. J.P.Burg, "Maximum entropy spectral analysis", in Proc. 37th Meeting Society of Exploration Geophysicists, Oct. 31, 1967 3. T.J.Ulrych and R.W.Clayton, Phys. Earth Planetary Interior, 1.~2, (1976) 188 4. A.H.Nuttall,"Spectral analysis of an univariate process with bad data points, via maximum entropy, and linear predictive techniques", Naval Underwater System Center, Tech. Rep., 5303, New London, CT, May 19, 1976 5. D.W.Tufts and R.Kumaresan, IEEE Trans.ASSP, ASSP-~, (1982) 671 6. S.Kawata, K.MINAMI and S.Minami, Appl. Opt., 22, (1983) 3593 7. K.Minami, S.Kawata and S.Minami, Appl. Opt., 2_~4, (1985) 162 8. L.Marple, IEEE Trans.ASSP, ASSP-28, (1980) 441 9. G.Zhang and J.Fu, Chin. J. of Microw. & Radiofreq. S p e c t . , 3 (4), (1986) 403 10. G.Zhang, Z.Xie and J.Fu, Chin. J. of Microw. & Radlo-freq. Spect., ~ (3), (1987) 275 11. Z.Xie and G.Zhang, J. of Zhongshan Univ.,

No.2, (1988) 55 12. S.Ueban et al, J. of Molec. Spect., 88, (1981) 274