Design of optimized finite impulse response digital filters for use with

1760

Anal. Chem. 1990, 62, 1768-1777

wavelengths shorter than that available from the Hg arc lamp. A demonstration of the use of this detector for the analysis of a mixture of iodinated and chlorinated hydrocarbons in a real air sample is shown in Figure 9. The PDM-ECD simultaneously provides the normal ECD response to this sample shown in Figure 9A and the PD-modulated response shown in Figure 9B. While responses to all of the halogenated hydrocarbons are provided by the normal ECD response in Figure 9A, responses to only the two iodinated compounds are provided by the PD-modulated response in Figure 9B. It is significant to note that ethyl iodide emerges from the GC along with a 1O-fold excess of chloroform. Nevertheless, the response to ethyl iodide in Figure 9B is unaffected by the simultaneous presence of chloroform in the detector and, therefore, provides an unambiguous means of determining the concentration of ethvl iodide in this samde. With the aDparatus described here, the detection limitof the PDM-EC'D for CH,I in air is about 0.01 ppb. Anticipated improvements in the electronic components of this instrument should lead to additional reductions in detection limits.

(4) Mock, R. S.; Zook, D. R.; Grimsrud, E. P. Int. J . Mass Spectrom. Ion Processes 1989. 9 1 . 327-338. (5) Mock, R. S.; Grimsrk, E. P. Int. J . Mass Spectrom. Ion Processes 1989, 9 4 , 293-303. (6) Wentworth, W. E.; Chen, E. C. M. In Nectron Capture, Theory and Practice in Chromatography; Ziatkis, A., Pooie, C. F., Eds.; Elsevier: New York, 1981; pp 27-68. (7) Zook. D. R.; Knighton, W. B.; Grimsrud, E. P. Int. J . Mass Spectrom. Ion Processes, in press. (8)Warman, J. M.; Sauer, M. C. J . Chem. Phys. 1975, 62, 1971-1981. (9) Alajajian, S. H.; Bernius, M. T.; Chutjian, A. J . Phys. B : At., Mol. Opt. Phys. 1988, 27, 4021-4033. IO) Aige, E.; Adams. N. G.; Smith, D. J . Phys. B: At. Mol. Phys. 1984, 17, 3827-3833. 11) Smith, D.; Adams, N. G.; Alge, E. J . Phys. B: At. Mol. Phys. 1984, 17 . . , 461-472 . - . . . -.

12) Orient, 0.J.; Chutjiin, A.; Crompton, R. W.; Cheung, B. Phys. Rev. A : Gen . Phys . 1989, 3 9 , 4494-450 1. 13) Patterson, P. L. J . Chromatogr. 1977, 734, 25. 14) Grimsrud, E. P. In Detectors for Capillary Gas Chromatography; Hill, H. H., McMinn, D., Eds.; Wiley & Sons, Inc.: New York, in press. 151 Adams. N. G.: Smith. D.: Herd. C. R. Int. J . Mass Soectrom. Ion Processes 1988, 8 4 , 243-253. (16) Grimsrud, E, p,; Knighton, W. B, Ana,, C b m , 1982, 5 4 , 565-570. (17) Grimsrud, E. P. Anal. Chem. 1984, 56, 1797-1803. (18) Smith, D.; Adams, N. G. J . Phys. 8: At. Mol. Phys. 1987, 2 0 , .^^^

4YUY-4Y 18.

(19) Mandl, A.; Hyman, H. A. Phys. Rev. Lett. 1973, 3 7 0 , 417-419. 13n\

ACKNOWLEDGMENT The authors express their appreciation to Douglas Zook for providing the APIMS spectra reported in Figure 6.

\-"I

Lhnrll m.,Gane"t,

A

Dhwc

n. "a,'.

Dnw

A.

llrl. n .

P2-

w a r .

Dhwr

,a','.

iIa 7e S,",

1 1 QAC QAQ 1 7 , V-rY-Y-r".

(211 Zook. D. R.: Grimsrud. E. P. J . Phvs. Chem. 1988. 92. 6374-6379. (22) Arshadi, M.; Yamdagni, R.; Kebark, P. J . Phys. Chem. 1970, 7 4 , 1475. (23) Hiraoka, K.; Mizuse, S.; Yamabe, S. J . Pbys. Chem. 1988, 9 2 , 3943-3952.

LITERATURE CITED (1) Ziatkis, A., Poole, C. F., Eds. ~ k c b o nCapture, TheoryandPractfce in Chromatography; Elsevier: New York, 1981. (2) Mock, R. S.;Grimsrud. E. P. Anal. Chem. 1988, 6 0 , 1684-1694. (3) Mock, R. S . ; Grimsrud, E. P. J . A m . Chem. SOC. 1989, 7 7 1 , 2861-2870.

RECEIVEDfor review March 9, 1990. Accepted May 18, 1990. This work was supported by grants from the National Science Foundation (Grant CHE-8711618) and from the Idaho National Engineering Laboratory, Idaho Falls, ID.

Design of Optimized Finite Impulse Response Digital Filters for Use with Passive Fourier Transform Infrared Interferograms Gary W. Small* and Amy C. Harms' Department of Chemistry, The University of Iowa, Iowa City, Iowa 52242 Robert T. Kroutil, John T. Ditillo, and William R. Loerop U S . Army Chemical Research, Development, and Engineering Center, Aberdeen Proving Ground, Maryland 21010

A deslgn technlque for flnlte Impulse response (FIR) dlgltal filters Is descrlbed that allows a narrow fllter bandpass to be achleved wlth a small number of fllter coeffklents. Based on multiple regression analysts, this approach produces filters that are tallored to the specific type of slgnal to which they will be applled. I n the present appllcatlon, fllters are developed for use In extractlng lnformatlon pertalnlng to single Infrared spectral bands dhectiy from Fourier transform Infrared (FTIR) Interferograms. The regresston-based technlque Is descrlbed and compared to two standard approaches to FIR fllter design. Interferograms collected by a passive FTIR remote sensor are used to evaluate the performance of the developed filters In terms of slgnal-to-noise ratios and computatlonal efflclency. I n addlon, the bandpass characterlstics of the computed filters are examined. On the basts of these evaluation criteria, the regresslon approach is judged to provlde superior performance In terms of Its ability to produce fitters that are computatlonally effklent and that achieve high slgnal-to-noise ratlos. Present address: Department of Chemistry, T h e Pennsylvania State University, University Park, PA 16802.

INTRODUCTION Passive Fourier transform infrared (FTIR)remote sensors represent an emerging technology for use in a variety of environmental monitoring applications. These devices consist of a conventional FTIR spectrometer in which the infrared source has been removed and a specialized set of input optics added. The infrared background emission in the field-of-view of the spectrometer is collected and analyzed for the presence of spectral bands corresponding to specific target analytes. Passive infrared spectrometers have come into recent use for remote sensing of atmospheric pollutants (1-4). Specific application environments for these sensors include monitoring a t hazardous waste sites, leak detection a t chemical plants, and regulatory monitoring of smokestack emissions. In these applications, the interferometer can be positioned in a stationary configuration or mounted in a ground or airborne vehicle. In either case, a premium is placed on the degree to which the sensor is compact, rugged, and reliable. The most fragile component of a typical FTIR remote sensor is the interferometer drive mechanism, which must allow the collection of a stable interferogram of 1024 or 2048 points. The required interferogram length is dictated by the spectral resolution required to detect the target analyte(s) of

0003-2700/90/0362-1768$02.50/00 1990 American Chemical Society

ANALYTICAL CHEMISTRY, VOL. 62, NO. 17, SEPTEMBER 1, 1990

interest. This relationship between interferogram length and spectral resolution derives from an inherent characteristic of the fast Fourier transform (FFT), the data processing tool used to extract infrared spectra from the collected interferograms. The FFT assumes that the interferogram is an infinitely long waveform that contains zeroes for all points not explicitly collected. This has the effect of adding (sin x ) / x components to the computed spectrum, resulting in spectral broadening. As the number of collected interferogram points is increased, the (sin x ) / x contribution is decreased. Similarly, spectra computed from very short interferograms are severely distorted due to these effects. Recently, work in our laboratories has focused on design concepts for an FTIR remote sensor based on a simplified “short-scan” interferometer design. The drive mechanism for such a system would allow the collection of only a 100-200point interferogram segment. Conceptually, this system would be much more rugged than a conventional design, as the moving mirror of the interferometer would need to maintain optical alignment for only a very short distance. The drawback to such a system is that a conventional spectral-based analysis could not be performed, due to the characteristics of the FFT noted above. The key to overcoming this problem is the ability to analyze a short interferogram segment directly for information pertaining to a target analyte. In an initial study, we have demonstrated the feasibility of this approach through the use of bandpass digital filters to extract from short interferogram segments information pertaining to characteristic frequencies of an analyte (5). This approach overcomes the limitations of the FFT but still allows the data analysis to be based on selected infrared frequencies. The success of this approach is determined by the ability to generate digital filters with an optimum bandpass. For the detection of a specific spectral band, the filter bandpass must be matched to that band. A wider than optimum bandpass will result in decreased sensitivity for the determination, as interferogram information other than that pertaining to the desired frequencies will be passed by the filter. Unfortunately, the number of computations required to apply a digital filter is related to the width of the bandpass. Narrow-bandpass filters require more computations than do wide-bandpass filters. For a typical infrared spectral band, the application of the optimum digital filter to a short interferogram segment may require more computations than the application of the FFT to the entire interferogram. In the work presented here, this drawback is overcome through the introduction of a new strategy for the construction of bandpass digital filters. Based on multiple regression analysis, this technique produces filters that require fewer computations than are typically required with filters generated by conventional design methods. T o evaluate this design strategy, filters are developed and compared to conventional filters in terms of bandpass characteristics, signal-to-noise ratios, and computational speed. EXPERIMENTAL SECTION Data used to test the filter design methodology were collected with a passive infrared interferometer constructed by Honeywell Corp. to specifications provided by the U S . Army Chemical Research, Development,and Engineering Center, Edgewood, MD. The system consisted of a flex-pivot “porch swing” Michelson interferometer based on the design of Walker and Rex (6). The spectrometer used a liquid nitrogen cooled Hg:Cd:Te detector to collect spectral background radiation in the region 800-1200 cm-’. For the work reported here, 1024-point interferograms were collected, with a corresponding nominal spectral resolution of 4 cm-‘.

1768

The data were collected in a mobile environment. The spectrometer was mounted on a shock-absorbingplatform and placed in both a helicopter and a ground vehicle. Interferograms were collected at a variety of speeds and distances from a ground source of SF6(MG Industries, Jessup, MD, 97% purity). sF6was selected as a test analyte because of its use as a standard test compound in pollution monitoring. It has a single strong absorption at 940 cm-’. All computer software used in this work was written in Fortran 77 and implemented on a Prime 9955 interactive computer system operating at the Gerard P. Weeg Computing Center at The University of Iowa. Some computations made use of subroutines from the IMSL library (InternationalMathematical and Statistical Library, Houston, TX). Plots were generated by use of the TELLAGRAF interactive graphics system (Integrated Software Systems Corp., San Diego, CA). A Hewlett-Packard7475A digital plotter was used as the output device. RESULTS AND DISCUSSION Application of Digital Filters to Direct Interferogram Analysis. The top left and top center plots in Figure 1 are single-beam spectra obtained by applying the FFT to two single-scan interferograms collected by the FTIR remote sensor described above. The top left spectrum corresponds to an interferogram acquired during the helicopter data collection, while the top center spectrum results from the ground-based data collection. Both spectra contain the spectral band of SF6 a t 940 cm-’. Due to the relative differences between the temperature of the analyte and that of the infrared background, the spectral band arises as an absorption in the top left spectrum and as an emission in the top center spectrum. The helicopter-based data collection produces an absorption band because the analyte is observed against the strong infrared emission from the surface of the earth. These spectra were ratioed to single-beam background spectra collected when the SF6source was out of the fieldof-view of the spectrometer, producing the corresponding transmittance spectra depicted as the bottom left and bottom center plots in the figure. In the transmittance spectra, the respective SF6absorption and emission bands are seen clearly. The uneven spectral baseline results from the change in the infrared background emission with the change in the fieldof-view of the spectrometer. The goal of any detection algorithm for the analysis of passive FTIR data is to determine if the spectral bands corresponding to the target analyte are present. The application of a bandpass digital filter helps to accomplish this goal by extracting the information in the region of a target spectral band from the remainder of the information in the spectrum. The bandpass of a filter is determined by its frequency response. The frequency response of a digital filter centered on the SF6 band is depicted in the top right plot of Figure 1 and expanded in the bottom right plot. The full width at half-maximum (fwhm) of this frequency response is 45.4 cm-’. A variety of functional forms can be used to define the frequency response. The frequency response depicted in Figure 1 is a simple Gaussian curve. While the Gaussian frequency response does not necessarily represent an optimum shape, we have found it to outperform other similar functions in our digital filtering work. For this reason, the remainder of the discussion here will focus on Gaussian-shaped frequency response functions. The action of a bandpass digital filter in both the spectral and interferogram domains is depicted in Figure 2. In the spectral domain, the filter is applied by multiplying the frequency response by the spectrum to be filtered. In effect, this operation superimposes the spectral features within the filter bandpass on the frequency response curve. In Figure 2, the top left, center left, and bottom left spectra result from the

1770


J

,."e"

i \ 100

IO00

1100

loo0

FREQUENCY (cm-')

FREQUENCY (cm-')

1.10

1.05

1.05

Y

z

FREQUENCY (cm-')

I

st

,

l

4

e 0.95

0.90

8

900

,000

FREQUENCY (cm-')

1

BOO

100

loo0

I

0

FREQUENCY (cm-')

Flgure 1. Left column: single-beam spectrum derived from the helicopter-based data collection (top) and Corresponding transmittance spectrum showing SFe absorption band at 940 cm-' (bottom). Center column: single-beam spectrum derived from the ground-based data collection (top) and corresponding transmittance spectrum showing SFe emission band at 940 cm-' (bottom). Right column: Gaussian frequency response function used to extract the SF, band from the single-beam spectrum (top) and expanded plot of this function (bottom).

multiplication of the frequency response depicted in Figure 1 by three synthetic transmittance spectra. From top to bottom, the synthetic transmittance spectra correspond to Lorentzian absorption bands with peak transmittances of 100% (i.e., no band), 8570,and 70%. In the resulting filtered spectra, the spectral band is observed superimposed on the frequency response of the filter. Synthetic spectra were used in this illustration, as the filtering concepts can be demonstrated more clearly without the presence of noise in the data. The filtered spectra were inverse Fourier transformed to obtain the corresponding filtered interferograms. Two segments of these interferograms are plotted in the corresponding center and right columns in Figure 2. The center column contains plots of interferogram points 1-80, while the right column contains plots of points 175-250. In this example and all others to follow, point number 1 in the interferogram is the maximum-intensity point. In an FTIR interferogram, this point is often termed the centerburst. An inspection of the interferogram data reveals two fundamental characteristics of interferograms that allow interferogram-based digital filters to be used in a detection application. First, a visual separation of information exists in the interferogram according to the width of the corresponding spectral feature. The three spectra produce interferograms that are virtually identical over points 1-80. The information that is visually dominant here corresponds to the shape of the frequency response of the filter. Since the same frequency response was used in each case, the interferogram segments are virtually identical. Points 175-250 are very different across the three interferograms, however. It is in this region that the information pertaining to the Lorentzian absorption band is apparent, as the signal due to the frequency response has

largely damped to zero. The second key point to note is that the magnitude of the interferogram segment is directly related to the magnitude of the Lorentzian absorption band. In fact, for the three interferograms, the square root of the sum of squares of points 175-250 has a correlation coefficient of 1.000 with the integrated area of the corresponding Lorentzian bands. This indicates that the magnitude of the interferogram segment can be used to obtain the area of the spectral band. The ability to use the interferograms in Figure 2 for a direct analysis was keyed by two operations. First, digital filtering was used to eliminate information outside of the frequency range of interest. This was essential, as each spectral frequency contributes to each interferogram point. The elimination of those frequencies outside of the region of 940 cm-' prevents other spectral features from contributing to the filtered interferogram, thereby allowing just the information pertaining to the targeted spectral band to be observed. Second, the interferogram was "windowed" to remove a segment in which the information pertaining to the shape of the filter bandpass had damped to zero. Performed together, the filtering and windowing operations allow qualitative and quantitative information to be extracted directly from a short interferogram segment. The drawback to the above procedure is that it began with a spectrum computed from a full-size interferogram. We wish to collect only a short interferogram segment, however. To overcome this problem, the action of multiplying the spectrum by the filter frequency response must be duplicated in the interferogram domain. Development of Interferogram-Based Digital Filtering Methodology. Mathematical Background. The action of multiplying a spectrum by the frequency response function


1771

0.03

I

1.0-

0.01-

0.1-

4.0-

* !4l t

Ez

>

0.6-

t

-

z

0.4

g

0.01-

1.0-

A

0.00-

Ez

-0.01-

-1.0-

-

-4.0

0.2

~

-0.02-

-6.0-

-0.03-

0 7

-0.O-l

a,ol

0.01-

0.1-

4.0

e v)

-

z

IW

z

>

0.6-

0.4

0.1

5

-

0.01-

:I

OA0-0.01-

-0.01-

-0.03-

0, 100

900

loo0

I100

,

,

,

210

230

-1.0 20

40

60

80

00

IBO

1 150

INTERFEROCRAM POINT 1.0 0.01

0.8-

I

-

a.0

-

4.0-

$

Ez

*

0.6-

0.4

t v) z

0.01-

*

t!

-

z

1.0-

t v) z

0.00-

0.0-

t! E

-0.01-

-2.0-4.0-

0.1-

-0.01-

,.

-6.0-

.

-

I

I - 1 3 . 1 800

900

FREQUENCY

loo0

I100

(em-')

IO 40 60 INTERFEROGRAM POINT

I O

no

190

110

130

150

INTERFEROGRAM POINT

Flgure 2. Left column: Gaussian frequency response function with a superimposed Lorentzian absorption band corresponding to 100% (top), 85 % (middle), and 70 % (bottom)transmittance. Center cdmnn: points 1-80 in the interferogram corresponding to the left column. Right column: points 175-250 in the interferogram corresponding to the left column.

of a digital fiiter has a corresponding form in the interferogram domain that is expressed by the convolution integral

y(t)= l m h ( k )x ( t - k) dk = P [ H ( f ,X O ] -m

(1)

where y ( t ) is the filtered interferogram, the product H ( f ) X o , is the filtered spectrum, H ( f )is the frequency response of the filter, X ( f )is the transformed spectrum, x ( t ) is the raw interferogram, and h(t)is the inverse Fourier transform of the frequency response, termed the impulse response of the filter. The terms H ( f ) and X ( f ) are functions of frequency, f, while the terms y ( t ) ,x ( t ) , and h(k)are functions of the time variables t or k. The inverse Fourier transform, P',relates the functions in the two domains. The above equation pertains to continuous functions. For digitized data, the discrete convolution integral has the form m

where yt is the intensity of point t in the filtered interferogram, and hk and xt-k define discrete points on the impulse response function and raw interferogram, respectively. The values yt and x t define the corresponding points on the filtered and unfiltered interferograms, respectively. Given that the interferogram is sampled over only a finite range, the discrete form of the convolution integral must be approximated by use of a finite series of terms. This approximation most often takes the form y t = h,,xt

+ ... + hnxt-,,

(3)

where n + 1 terms are summed to estimate each filtered interferogram point yt. The terms hk can be considered weighting coefficients that determine the frequency dependence of the filter. Each fiitered interferogram point is formed from a linear combination of the corresponding raw data point and a set of preceding raw data points. This type of filter, defined by the h k , is termed a finite impulse response (FIR) filter. From an inspection of eq 3, it can be seen that filtering a short interferogram segment requires only the n interferogram points preceding each point being filtered. Conventional FIR Filter Design. The hk terms must be found that best approximate the action of the filter frequency response function. A variety of techniques can be used to accomplish this task. For example, the frequency response can be inverse Fourier transformed to generate its corresponding form in the interferogram domain. Once obtained, this impulse response function can be truncated to the desired number of terms for use in eq 3. Alternatively, a polynomial series approximation to the frequency response can be generated. By approximating the frequency response with a f i i t e series of terms, the hk can be computed by direct calculation. The Remez exchange method (7, 8) is perhaps the most commonly used algorithm for generating a finite series approximation of the frequency response function. Figure 3 displays six frequency response functions obtained by use of these conventional techniques. The top three plots display the frequency responses corresponding to truncated versions of the Gaussian function depicted in Figure 1. The depicted Gaussian was inverse Fourier transformed and truncated to 60, 40,and 20 points. These points define the filter coefficients (hk)that would be used in eq 3. To obtain


1772

1-

' 7 0.8

-

0 W

25

0.6-

4

I w

VI

2

0.4-

Ya

IO

I

1

0.1

0.8

0 w

z5

2z

0.6 l

0

0.6-

I

w

W

VI

g

0.80 W

z In

0.4

0.b-

0 VI (L

Y a

a

Y

0.2

0.2

-

0

0

$00

IO00

1500

2000

0

FREQUENCY (cm-')

Flgure 3. Top row: frequency response functions obtained by use of truncated Gaussian functions containing 60 (left), 40 (center) and 20 (right) points. Bottom row: frequency response functions derived through application of the Remez exchange algorithm. These frequency responses are produced by 60 (left), 40 (center), and 20 (right) filter terms.

the corresponding frequency responses, the truncated functions were then forward transformed to produce the top left, top center, and top right plots in Figure 3. From left to right, the fwhm values for these frequency response functions are 84.0, 122.4, and 234.0 cm-'. The lower three frequency responses were produced by the Remez exchange algorithm. In each case, the targeted frequency response was the Gaussian frequency response displayed in Figure 1. Input to the Remez exchange method is a series of stop-band and pass-band specifications that encode the desired frequency response. T o model the Gaussian function, the regions of zero response were specified as the stop bands (0-844.8, 1037.5-1974.75 cm-l), and the region at the top of the Gaussian curve was specified as the pass band (939.2-943.1 cm-'). While the Remez exchange algorithm was not designed specifically to produce a Gaussian-shaped response, the focus here is simply to obtain a response of the desired width. If that width is achieved, the obtained response will be shaped similar to a Gaussian function. The three frequency responses plotted differ in the number of terms used in the finite series approximation. From left to right, the frequency responses corespond to 60,40, and 20 terms. The fwhm values for the curves are 102.4, 124.4, and 194.4 cm-'. These two sets of frequency responses confirm the assertion made above that increasing numbers of filter terms are required to achieve a narrower filter bandpass. The signal-tonoise ratio of the frequency response also improves with an increase in the number of filter terms. The generated filters were tested by applying them to points 175-250 in a series of interferograms collected during two data runs. Before application of the bandpass filter, each interferogram was filtered with a four-term high-pass filter to remove low-frequency noise and subsequently normalized on the basis of the sum of squares of points 175-250. This

two-step normalization process is the best method found for correcting for differences in overall infrared energy between interferograms. The high-pass fiiter is required, as the power supply in the instrument used here introduces low-frequency noise components into the interferogram. The presence of low-frequency noise interferes with the normalization calculation. After bandpass filtering, the sum of squares was computed across the range of filtered points. A plot of the computed sums of squares vs interferogram number provides an effective means for evaluating the performance of a filter. Figure 4 depicts a series of six such plots for a set of interferograms collected with the spectrometer mounted on a helicopter flying at 80 knots and at an altitude of lo00 ft. Five passes were made past a ground source of SFG.The six plots derive from the application of the six filters whose corresponding frequency responses are depicted in Figure 3. As motivated by the discussion of Figure 2, the peaks in each sum-of-squares plot (i.e., large values of the sum of squares) correspond to interferograms containing SF, information. The same scale is used for each plot. The expected trend is observed of an increase in signal strength with an increase in the number of filter terms. This trend is observed in both sets of filters. Figure 5 is an analogous series of plots for a set of interferograms collected when the spectrometer was mounted on a ground vehicle moving at 10 mph past an SF6source. One pass was made past the SF6source, as indicated by the single peak in each plot. The same scale is used for each plot. Again, the strength increases with an increasing number of filter terms, although the baseline noise in the plot also increases as the number of filter terms increases. In comparing Figures 4 and 5, much larger SF6 signals are always observed in the helicopter-based data, as the spectrometer is looking a t the target gas against the earth, a strong infrared radiator.


2 g

1773

0.3-

0 0.1-

1

m

0.1-

,

0.04

n

n

n

I

h

INTERFEROGRAM NUMBER

0.7

0.1

0

z 0.5

c K

t

B

t

0.1

4

D

4

D

4

0.1

0.1-

0 3

0

&

0.4-

4

&

0.1

I

0.1-

r

I VI

VI

\

0.1-

0.1

,

o

I

IOO

0.03

,

zoo 100 400 a00 INTERFEROGRAM NUMBER

roo

i -

700

J-:


I

0.0

1


Flgure 4. Plots of interferogram sum of squares after filtering vs interferogram number for the helicopter-based data set. The plots were produced by use of the filters whose frequency responses are depicted in the corresponding plots of Figure 3. 0.1,

0 24

z 5

z 0.11

5

2

K

K

b-

Y

4

4

Y

I 0

0

0

z

5 '1

0.16

K

t m

Y

o.12

0

0.12.

K 3

2

4

0.11-

I&

m

!A k.

;

; k.

0.0s

0.06

3

YI 3

0.00

0 00 100

100 400 aoo INTERFEROGRAM NUMBER

100

IOO

0

600

zoo 100 400 aoo INTERFEROGRAM NUMBER

600

io0

o

100

zoo

100 400 aoo INTERFEROGRAM NUMBER

600

0.24

0

0

z

3

8

2

g

0.16

2

0.11

B

K

c

L

4 0.12

K Y

3

3

3 4

y)

L?

YI

;L ;

; I&

0.06

0.06

3

0.00 IO0

200

100

400

500


600

:

0.00 0

100

ZOO INTERFEROGRAM 100 400NUMBER 500

I00

- 0 0 . 0


Figure 5. Plots of interferogram sum of squares after filtering vs interferogram number for the ground-based data set. The plots were produced by use of the filters whose frequency responses are depicted in the corresponding plots of Figure 3.

1774

ANALYTICAL CHEMISTRY, VOL. 62,

NO. 17, SEPTEMBER

1, 1990

The results described above serve to motivate our work in the development of new filter design techniques. To extract small analyte signals, a filter bandpass of optimum width and signal-to-noise must be obtained. For the example used here, the targeted fwhm value (45.4 cm-') for the filter frequency response could not be achieved, even with 60 filter terms. Given application environments in which multiple analytes must be detected, it may be necessary to apply several filters. Clearly, for the interferogram-based method to be viable, filters must be developed that achieve acceptable bandpass and signal-to-noise characteristics with as few terms as possible. Regression Approach to FIR Filter Generation. As noted, eq 3 is a linear model that relates a filtered interferogram point, yt, to a series of n + 1raw interferogram points, xt+ Regression analysis is the branch of statistics that focuses on techniques for the construction and evaluation of such models. In regression terms, the yt define a dependent variable, while the Xt-k defiie a set of independent variables. The hk are regression coefficients. In our initial work (5),it was judged potentially fruitful to investigate the application of regression techniques to the construction of FIR filters. Two motivating factors led to this conclusion. First, in any regression model, some independent variables are always more statistically significant than others. These significant variables contribute most to the explanation of variance in the dependent variable. It was hypothesized that if regression procedures could be used to estimate the hk,it may be possible to delete some of the Xt-k terms as being insignificant. The resulting filters would thus have fewer terms. Second, an inspection of eq 1 reveals that, in a conventional FIR filter, the hk depend only on the frequency response function. In one sense, this is an advantage in that a set of hk can be used to filter any type of signal and any segment of that signal. For the present application, however, the only signal encountered is the interferogram produced by the FTIR interferometer and detector. Moreover, as demonstrated previously, it should be possible to base a detection algorithm on a short segment of that interferogram. I t was hypothesized that fewer filter terms would be required if a regression approach could be used to derive a set of hk that are dedicated to filtering specific segments of FTIR interferograms. The generation of an FIR filter via regression methodology begins as before with the definition of the frequency response of the filter. For example, to generate an SF6filter, a frequency response such as that displayed in Figure 1would be used. As noted previously, if this Gaussian function is multiplied by a sample spectrum, a filtered spectrum results. The inverse Fourier transform of this filtered spectrum is a corresponding filtered interferogram. If the desired interferogram-domain filter were available, the action of the filter on the sample interferogram would produce a filtered interferogram identical with that obtained through the Fourier transform/Gaussian multiplication/inverse transform step outlined above. Thus, points in this generated filtered interferogram define the dependent variable for the regression analysis. The specific points used would be those defining the interferogram segment for which the dedicated filter is desired. A standard multiple linear regression analysis can thus be performed to derive the desired hk. In this calculation, a pool of potential independent variables can be used at k = 0, ..., nmax.These variables would be derived from points in the same raw interferogram whose transformed spectrum was used in the generation of the dependent variable. Thus, in regression terms, the hk are derived by solving

b = (XTX)-'XTY (4) where X is a matrix whose columns are the independent

variables, XT is the transpose of X, Y is a vector that defines the dependent variable, and b is a vector that contains the regression coefficients (i.e., the hk). As an example, if the filter were to be generated for points 175-177 and independent variables from k = 0 to k = 2 were used, X would be defined as

2 1:: x175

x=[

X174

:.:1 X173

(5)

where X I 7 5 is point 175 in the raw interferogram, X I 7 6 is point 176, etc. Similarly, the dependent variable, Y,would be defined as where ~ 1 7 is 5 point 175 in the filtered interferogram obtained through the Fourier transform/ Gaussian multiplication/ inverse Fourier transform step. The other elements in y are defined similarly. Thus, the independent variable corresponding to k = Z(X173, xlT5)represents the collection of raw interferogram points displaced two points behind each of the points being filtered. Stated differently, the first row of X represents the collection of points to be used in attempting to model the fitered interferogram point whose value is the first element of Y. In our initial work (5), we added consecutive terms to the model, beginning with k = 0 and continuing until no further variance in the dependent variable was explained by an added term. While this procedure produced useful filters, it was clear that many of the added terms were not statistically significant. In a second study (9), a stepwise regression algorithm (IO)was used to select only those 7xt-k that contributed strongly to the model. This algorithm begins by selecting the single independent variable that has the highest correlation with the dependent variable. A one-term model is then formed based on this selected variable. Consecutive terms are added to the model in a stepwise manner. At each step, the variables remaining in the pool are evaluated for their correlation with the variance in the dependent variable that has not been explained by the terms previously selected. The variable chosen through this procedure is added to the model, and the process is repeated for a fixed number of iterations or until no remaining variables in the pool meet a minimum standard of correlation. FIR Matrix Filters. An inspection of eq 1 reveals that each point on the filtered interferogram is derived through the separate evaluation of the convolution integral for that point. As noted previously, a filter that operates on the interferogram must approximate this convolution integral at each point. In a conventional FIR filter, however, the same filter coefficients and the same surrounding points are used to define each filtered point. To achieve the best approximation of the convolution integral, however, it can be argued that a different fiiter model should be used at each point. For each point, such a filter would have a different set of filter coefficients and be based on a different set of surrounding points. By deriving separate f i t e n to be used on each interferogram point, it was hypothesized that the overall accuracy of the filtering process could be increased. It is also true that such a procedure may be faster computationally. Since each filter is now perorming a simpler task, fewer terms in each filter may be needed. The derivation of a set of filters is a direct extension of the stepwise regression methodology discussed above. Separate stepwise regressions are performed, one for each interferogram point treated. The principal difference in this case, however, is that a set of interferograms must be employed to define the

ANALYTICAL CHEMISTRY, VOL. 62, NO. 17, SEPTEMBER 1, 1990 RAW INTERFEROGRAM

1775

FILTERED INTERFEROGRAM

0.6-

1

0

z

W E

3

0.5-

e W E

0.4-

6

/

fi E

6 0.3-

3

a wl

DEPENDENT VARIABLE FOR POINT 180 FILTER

REGRESSION ANALYSIS

0.2-

I 3 VI

0.1 -

Figure 6. Schematic of the assembly of the dependent and independent variables used in the computation of a filter optimized for a single interferogram point.

L 0

200

100

400

300

500

I


dependent and independent variables for the regression analysis. The dependent variable would be formed as before through the Fourier transformation of each interferogram in the set, multiplication by the frequency response of the filter, and inverse Fourier transformation back to a set of filtered interferograms. For example, the dependent variable needed to generate a filter for point 180 would be formed by assembling each point 180 across the set of filtered interferograms. Analogously, the k = 1 independent variable would be formed from each point 179 across the set of raw interferograms. The number of observations in each variable would be equal to the number of interferograms in the set. Figure 6 provides a graphical illustration of the formation of the dependent and independent variables. The stepwise regressions produce a set of filters stored as a matrix of filter coefficients. The columns of this matrix contain the individual filter coefficients for a given interferogram point. On the basis of this storage format, we have termed the resultant filter set a “matrix filter”. This filter design strategy was tested by developing a matrix filter for SF6 based on the same frequency response and interferogram segment used above. A set of 3200 interferograms was used in this calculation. This interferogram set contained interferograms derived from both ground-based and helicopter-based data runs. The interferogram normalization procedure described previously was used. The pool of independent variables used at each point was defined as k = 0-100, and the stepwise regression procedure was used to compute 30 coefficients for each element of the matrix filter. The mean of the R2values for the 76 stepwise regressions was 97.8%. The 0-100 range and 30-coefficient limit were used after investigating several ranges and filter sizes. Filters with greater numbers of coefficients performed somewhat better in testing, but the 30-term filters were judged superior in terms of performance vs speed of application. The matrix filter was applied to the helicopter and ground-based data sets discussed previously. The interferograms used in the generation of the filter were not part of either of these data sets. The four-term high-pass filter was applied to each interferogram, along with the normalization procedure based on points 175-250. Following application of the matrix filter, the sum of squares was computed across points 175-250 as before. The top of Figure 7 is a plot of the computed sum of squares vs interferogram number for the helicopter data run, while the bottom of the figure depicts the corresponding plot for the ground-based data. Quantitative Comparison of SF6Filter Performance. The filtering results presented above are best compared by evaluating three criteria: (1)filter bandpass characteristics, (2) signal-to-noise ratio of the sum of squares plots, and (3) computational efficiency.

i Y W

L c

6 0 W

0.12-

w

2a VI

0.06 wl

0

100

200

300

400

500

600

;


Figure 7. Plots of interferogram sum of squares after filtering vs interferogram number for the helicopter-based (top) and ground-based (bottom) data sets. The matrix FIR filter was used to produce both plots. The plot scaling is identical with that used for the corresponding plots in Figures 4 and 5.

Comparison of Bandpass Characteristics. The frequency response curves depicted in Figure 3 can be evaluated in terms of their fwhm value, as well as their ability to attenuate frequencies outside of the filter bandpass. The matrix filter possesses a different set of filter coefficients for each point, thus complicating the calculation of an overall frequency response function. This problem is overcome by applying the fiter to a series of synthetic single-frequencycosine waveforms. At a given frequency f , a cosine waveform of unit amplitude is defined as A , = COS (27rft) (7) where A , is the amplitude of the waveform at time t. Values off range from 0 to 0.5, where 0.5 corresponds to one-half of the data sampling rate (i.e., the maximum frequency observable as defined by the Nyquist criterion). The frequency range was divided into 512 increments, corresponding to the same frequency spacing used in Figures 1-3. A t frequency f, the response magnitude of the filter R, is computed as

R, =

[e i=l

( Y * ~- j j * ) 2 / ?

i=l

(yi - j92]1/2

(8)

where the y*i are the intensities of the n points in the filtered segment of the synthetic cosine waveform, jj* is the mean intensity in the filtered segment, and yi and j j are the corre-

1776


0

0

20

40

60

80

100

120

140

NUMBER OF FILTER COEFFICIENTS

Figure 8. Frequency response produced by the 30-term matrix filter. The singWxam spectrum produced by the passive FTIR spectrometer is depicted with the dashed line. 250

-

'E 2.

200

I

3

E

I50

4

20

2

100

60

80

100

I20

0

Signal-to-noise ratios computed from the sum-of-squares plots presented for the helicopter-based (top) and ground-based (bottom)data sets. The Gaussian filters are indicated by solid circles, while the filters generated with the Remez exchange algorithm are indicated by open circles. The signal-to-noisevalues produced by the matrix filter are indicated by horlzontal dashed lines, while the value of 30 coefficients for the matrix filter is indicated by the vertical dashed lines. Figure 10.

c

I .b

2

40

NUMBER OF FILTER COEFFICIENTS

2 L

50

0

0

,

,

I

20

40

60

,

80

,

100

#

120

Figure 9. Full width at half-maximum for the computed filters. The Gaussian filters are indicated by s o l i circles, while the filters produced by the Remez exchange algorithm are indicated by open circles. The analogous value for the matrix filter is indicated by the horizontal dashed line. The vertical dashed line indicates that the matrix filter is based on 30 Coefficients.

sponding values in the raw (i.e., unfiltered) waveform. Equation 8 computes the ratio of the root-mean-square (rms) powers of the filtered and unfiltered signals. Those frequencies attenuated by the filter will produce R, values near 0, while those frequencies passed by the filter will produce R, values near 1. In this calculation, the matrix filter is applied to points 175-250 in each cosine waveform. Thus, n in eq 8 equals 76. Figure 8 depicts the computed frequency response for the 30-term matrix filter. Superimposed on the plot is a single beam spectrum produced by Fourier transforming one of the passive FTIR interferograms. The frequency response of the matrix filter is clearly tied to the infrared detector response. The filter sacrifices performance in frequency regions in which the specific infrared detector used dictates that no signal will be present to filter. In this scheme, the infrared detector is also used as a filter. The regression-based computation automatically builds into the filter information regarding which frequencies will be present. In this way, the number of filter coefficients can be minimized. Figure 9 compares the fwhm value of the matrix filter to those obtained by use of the truncated Gaussian (solid circles) and Remez exchange (open circles) algorithms. In this comparison, conventional filters consisting of 20,40,60, 80, 100, and 120 coefficients were used. The fwhm value for the matrix filter (82.1 ern-') is indicated as a horizontal dashed line in

the plot. The vertical dashed line markes the value of 30 coefficients for the matrix filter. The intersection of the dashed lines with the plotted curves allows comparisons to be made between the different filter design strategies. An examination of Figure 9 reveals that the fwhm value for the matrix filter corresponds to conventional filters requiring twice the number of coefficients. Conversely, both conventional filters with 30 coefficients would not perform nearly as well as the matrix filter. The ability to tailor the filter coefficients to each point and to the signal being filtered clearly allows the number of filter coefficients to be reduced. Evaluation of Sum-of-Squares Plots. On the basis of inspections of transformed spectra, SF6-activeregions were established in both the helicopter- and ground-based data runs. The SF6-activeregions in the helicopter-based data set correspond to interferograms 231-246, 320-333, 419-427, 532-535, and 599-605, while the single SF,-active region in the ground-based data set corresponds to interferograms 358-387. A baseline set of 200 interferograms was selected in each data set from among those interferograms judged to contain no SF,. For each plot, the sum of squares values for these baseline interferograms was used to compute a standard rms noise level. In addition, simple linear regression was used to define a baseline model for each plot. For each interferogram in the SF,-active regions, the baseline model was used to compute a predicted baseline value, and the difference was computed between the sum of squares value for the interferogram and the baseline value. The ratio of this difference to the rms noise was taken as the signal-to-noise ratio. This computation corrects for sloping baselines. In a manner analogous to that of Figure 9, Figure 10 com-


pares the computed signal-to-noise ratios for the sum-ofsquares plots produced by the computed filters. The upper plot in the figure describes the application of the filters to the helicopter-based data, while the lower plot presents the corresponding information for application of the filters to the ground-based data. The Gaussian filters are again indicated by solid circles, while the filters computed by use of the Remez exchange algorithm are indicated by open circles. The results depicted for the helicopter-based data represent mean values computed over the five SF6-active regions in the data set. The signal-to-noise ratio produced by the matrix filter is indicated by a horizontal dashed line, while the value of 30 filter coefficients is again indicated by a vertical dashed line. The signal-to-noise ratios produced by the conventional filters reach a maximum and then decrease in three of the four cases. This phenomenon indicates that the obtained signal-to-noise is clearly related to the filter width. The optimum width appears to be different for the two types of data. The curve for the filters produced by the Remez exchange algorithm does not invert for the helicopter-based data, as none of these filters reaches the width of the filter producing the maximum signal-to-noise value, the 100-term Gaussian filter. For the helicopter-based data, the 30-term matrix filter performs in a manner analogous to conventional filters with 80 terms or greater. For the weaker ground-based data, the matrix filter performs in a manner that is superior to each of the conventional FIR filters. The observed performance of the matrix filter is better than what would be expected on the basis of fwhm values in Figure 9. The significantly enhanced performance observed with the lower signal-to-noise data indicates to us that the removal of statisticaly insignificant terms from the convolution approximation greatly enhances the performance of the filter in low signal-to-noise cases. The conventional filters suffer in performance due to the presence of these insignificant (i.e., noise-dominated) terms. While the 100-term Gaussian filter and the 120-term Remez exchange filter do achieve better performance for the high signal-to-noise helicopter-based data, we consider improvements in signal-to-noise above a value of 200 to be effectively not worth the computational cost of the additional filter coefficients. The key to evaluating the filters lies in the consideration of the filter performance for the low signal-tonoise data. The fact that the matrix filter achieves superior performance with a 62.5% reduction in the number of filter coefficients represents an extremely significant advancement. Evaluation of Computational Efficiency. As described by eq 3, the number of mathematical operations required to implement an FIR filter is linearly related to the number of filter coefficients. Thus, a matrix filter requiring 30 filter coefficients can be applied in half the time that is required for a conventional FIR filter containing 60 coefficients or in one-third the time required for a 90-term filter. For comparison to conventional spectral-based processing, the application of the FFT to a 1024-point interferogram requires 20 480 multiplications. A 30-term matrix FIR filter applied to a 76-point interferogram segment requires 2280 multiplications. Thus, if an equivalent detection can be made by analyzing a short-filtered interferogram segment, an 88.9% savings in computational time is realized vs the spectral-based analysis. Furthermore, this comparison underestimates the computational cost of the spectral-based analysis, as no consideration is given to the extraction of the analyte information from the transformed spectrum. As noted in the discussion of Figure 2, this extraction of the analyte information is automatically a part of the interferogram-based filtering procedure.

1777

CONCLUSIONS On the basis of the work described here, the matrix FIR filter design method offers substantive advantages over conventional filter design techniques. The matrix filter approach offers a means to obtain narrow-bandpass filters that require a minimum number of computations to apply. The principal drawback to the regression-based method is that it produces filters that are highly specific to a given type of signal (e.g., a specific segment of an FTIR interferogram). In the majority of applications, however, this drawback is not significant. Additionally, the increased computational time required to generate separate sets of filter coefficients is not significant, as these computations are only performed once. For passive FTIR remote-sensing applications, we consider the most significant achievement of this work to be the ability to perform a detection by use of only a short interferogram segment. In addition, the interferogram-based analysis allows analyte information to be extracted without the necessity of having a separate infrared background measurement. The lack of a stable infrared background in the remote-sensing application severely complicates the use of conventional spectral-based analysis strategies with passive FTIR data. The significant savings in computational time provided by the matrix filter is an added benefit. While the matrix filter design strategy was evaluated here in terms of applications to passive remote sensing, the approach is similarly applicable to other types of signals. As an example, interferograms collected with conventional laboratory FTIR spectrometers differ in no substantive way from the interferograms employed here. This suggests that an analysis based on filtered interferogram segments may be useful in other applications in which FTIR monitoring is used (e.g., process monitoring/control, chromatographic experiments). Many of these applications have similar requirements to remote sensing in terms of the ruggedness needed in the spectrometer hardware. While many process monitoring problems require the analysis of complex condensed-phase mixture spectra, the ability to extract from the interferogram information pertaining to specific spectral bands may still be very useful. In principle, a multivariate analysis could be built around the application of several different bandpass filters to an interferogram segment. A multivariate analysis of this type might allow the effects of overlapping spectral signals to be overcome, while still allowing the potential simplification of the spectrometer hardware made possible by an analysis requiring only a short interferogram segment. LITERATURE CITED Hanst, P. L. I n Advances in Environmental Science and Technology; Pitts, J. N., Metcalf, R. L., Eds.; Wiley: New York, 1971. Chan, S.H.; Lin, C. C.; Low, M. J. D. fnviron. Sci. Technol. 1973, 7 , 424-428. Walter, H.; Flanigan, D. Appl. Opt. 1975, 14, 1423-1428. Low, M. J. D.; Clancy, F. K. Environ. Sci. Techno/. 1967, 1 , 73-74. Small, G. W.; Kroutil, R. T.; Ditillo, J. T.; Loerop, W. R. Anal. Chem. 1988, 6 0 , 264-269. Walker, R. A.; Rex, J. D. Proc. SOC. Photo-Opt. Instrum. f n g . 1979, 191, 88.

McClellan, J. H.; Parks, T. W. I€€€ Trans. Circuit Theory 1973, CT20, 697-701. McClellan, J. H.; Parks, T. W.; Rabiner, L. R. I f f € Trans. Audio Electroacoust. 1973, AU-21, 506-525. Bjerga, J. M.; Small, 0 . W. Anal. Chem. 1989, 6 1 , 1073-1079. Draper, N. R.; Smith, H. Applied Regression Analysis, 2nd ed.;WileyInterscience: New York. 1981.

RECEIVED for review March 28,1990. Accepted May 18,1990. This work was supported by the U S . Army Chemical Research, Development, and Engineering Center, Edgewood, MD, under Contracts DAAA15-86-C-0034 and DAAA15-89c-0010.

Design of optimized finite impulse response digital filters for use with

Recommend Documents