Data Filtering in Instrumental Analyses with Applications to Optical

Sep 7, 2011 - To support such lectures by visual demonstrations and hands-on learning, a windows-based software package, DENOISE.exe, has been develop...
1 downloads 6 Views 3MB Size
ARTICLE pubs.acs.org/jchemeduc

Data Filtering in Instrumental Analyses with Applications to Optical Spectroscopy and Chemical Imaging Frank Vogt* Department of Chemistry, University of Tennessee, Knoxville, Tennessee 37996-1600, United States

bS Supporting Information ABSTRACT:

Most measurement techniques have some limitations imposed by a sensor’s signal-to-noise ratio (SNR). Thus, in analytical chemistry, methods for enhancing the SNR are of crucial importance and can be ensured experimentally or established via pretreatment of digitized data. In many analytical curricula, instrumental techniques are given preference as proper data generation is most important. Nonetheless, the ultimate goal is to utilize computational improvements as well and thus students need to be trained in data pre-treatment tools. This requires teaching a rather extensive math background for which there is a lack of concise teaching materials tuned to analytical chemistry. To overcome this hindrance, three methods for data filtering are presented: Savitzky-Golay smoothing, Fourier filtering, and the more recent and powerful wavelet filtering. By means of 1D signals as acquired with standard instrumentation for optical spectroscopy or chromatography and so forth, the basic principles are discussed and demonstrated. These methods are then expanded towards noise filtering of 2D data sets as obtained in microscopy or the upcoming chemical imaging as well as towards 3D data from spectroscopic imaging. In an extensive Supporting Information section, suggestions on introducing Fourier and wavelet transforms in lectures are presented. To support such lectures by visual demonstrations and hands-on learning, a windows-based software package, DENOISE.exe, has been developed and made available along with the simulated and real-world data used in this article. Because data are imported into DENOISE.exe as ASCII text files, the program can be utilized for any data sets an instructor may want to use. KEYWORDS: Graduate Education/Research, Upper-Division Undergraduate, Analytical Chemistry, Chemoinformatics, Computer-Based Learning, Spectroscopy

N

o experimental measurement is free of uncertainties induced by random processes that are due to very fundamental but unavoidable effects such as quantization of electrical charge (shot noise), Johnson noise induced by thermal movements of charge carriers, or the omnipresent 1/f noise.1 As any chemical analysis’ precision increases with the signal-to-noise ratio (SNR or S/N), it is imperative to improve SNR for any measurement technique. Systematic errors in a measurement impacting a method’s accuracy, however, should always be removed from an experiment prior to acquiring data. Thus, random signal fluctuations (noise) are the main focus of this article rather than systematic measurement errors. There are two general types of data filtering: prior to analog-to-digital conversion or afterward. In the former case, which will not be considered in this article, the filtering is typically done by means Copyright r 2011 American Chemical Society and Division of Chemical Education, Inc.

of electronic circuits. In the latter case, which is the topic of this article, digitized signals are manipulated computationally as a data pretreatment. For teaching purposes, a dedicated Windows-based software package DENOISE.exe has been developed, which is accompanied with simulated and experimental test data. Software, data, and a discussion of the mathematical foundation are concluded in the Supporting Information section and can be found in more detail for instance in references.25 The presentation of noise (or better frequency) filtering techniques begins with 1D, multivariate signals68 such as optical spectra, chromatograms, or any other time series of univariate

Published: September 07, 2011 1672

dx.doi.org/10.1021/ed100984c | J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 1. Moving average filters of different window widths W = 5 (left) and W = 41 (right) applied to a medium noise level (#3 in the Supporting Information, N, the number of individual points, = 4096). This and subsequent graphs have been determined by means of DENOISE.exe.

measurement points. These methods are then expanded for treatment of 2D data sets such as images as acquired, for example, in microscopy or chemical imaging. As an outlook, these ideas are further expanded to filtering of 3D data cubes as obtained in spectroscopic imaging.

’ METHODS OF DATA FILTERING Three denoising methods are presented in this article: (A) The Savitzky-Golay algorithm (SavGol),911 which is based on a modified data averaging. Although this method is over 40 years old and has considerable drawbacks, it still is popular and implemented in many instrument software packages. This fact and comparison purposes motivated its inclusion here. (B) Fourier filtering2,4,9 suppresses signal components whose frequencies do not contain relevant information. This method is arguably the most intuitive and interpretable one. (C) Wavelet-based data filtering is the most modern time frequency analysis that determines what frequencies are present where in the signal. Thus, wavelets overcome the limitation of Fourier methods, which completely remove certain frequencies although they might be informative in places. For instance, a Fourier low pass filter interprets all high frequencies as unwanted noise, although they are important whenever the signal changes sharply. A wavelet-based filtering cancels only those high frequencies that are not required to preserve relevant signal features. However, much of the wavelet literature is math driven and the goal here is to present this method from an analytical chemistry perspective (introductory literature, refs 9, 1215; more advanced, refs 4, 1619). Nonetheless, it should be clearly stated to students that a complete separation between relevant signal components and noise is nearly never feasible. This in return means that a data pretreatment intended to remove unwanted compounds most likely also cancels relevant signal components as well; vice versa, if noise filtering is too conservative, artifact removal remains incomplete. For a first demonstration, a synthetic signal with six different noise levels has been generated that consists of Gaussian-shaped peaks of different widths (Figure 1). By means of these peak widths, the filters’ impacts on features of different frequency composition can be studied.

Low Pass Methods Based on Data Averaging

A moving average is the simplest noise reduction filter, which is typically applied when the signal level changes slowly compared to the signal sample rate (= the number of measurement points per time interval). In such situations, it can be assumed in good approximation that signal fluctuations among consecutive data points are due to noise and averaging several measurement points does not considerably change the underlying, true signal but “averages out” fast fluctuations. Such a filter computes the average of a user-selected number of W consecutive measurement points and replaces the signal’s value at this window’s center with this average. Then, this window is moved to the right by one data point and the averaging-and-replacing procedure is repeated. Thus, the window width W is the only required filter parameter and is usually chosen to be an odd number to ensure a unique window center point: W ¼ 2m þ 1

m ¼ 1, 2, :::

ð1Þ

Often, it is helpful to state the window width W in percent of the total signal length N with 5% being a good starting point for optimizing such a filter for a specific application. Figure 1 depicts two examples of how W impacts the smoothing of a synthetic signal. Obviously, a wider window achieves more noise suppression but the needle-shaped (real) signal features on the right end are distorted considerably. Hence, for signal features narrower than W, denoising is traded for lowering and broadening peak heights. A moving average does not introduce such a bias though if the signal is constant or linearly changing in time, a signal is only falsified, if its second derivative is not equal to 0. To reduce such unwanted signal deformations, Savitzky-Golay (SavGol) filters911 have been introduced. Instead of averaging the data inside the moving window W, eq 1, they are fitted to a polynomial of order k g 1 pðxÞ ¼ a0 þ a1 x þ a2 x2 þ ::: þ ak xk

ð2Þ

that follows the data features more flexibly. If no x axis is available, the x positions in eq 2 are simply the integer numbers ranging from m e x e m. This also simplifies determining the polynomial’s value at the window center p(x = 0) = a0; a0 then replaces the original measurement point at this position. For SavGol, the user has to select W and k, whereas W > k + 1 is 1673

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 2. Impact of the polynomial order, eq 2, on denoising and preservation of narrow signal features; k = 3 (left) and k = 5 (right); W = 41 was chosen for both.

Figure 3. Examples for filter functions F(ω), eqs 47; as the filter function is symmetric with respect to the y axis, only positive frequencies are shown here for clarity purposes.

required in eq 2 to introduce least-squares fitting. Furthermore, at both ends of the signal, m measurement points are lost because they never can be center of W. It is also noteworthy, that for k = 0, SavGol collapses to the standard moving average.

Figure 2 compares two polynomial orders (k = 3 and k = 5) for the same W and obviously the impact of k is less than that of W (cf. Figure 1). Lower-order polynomials produce a slightly smoother signal than higher ones because the latter has a tendency to follow 1674

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 4. Fourier low pass filtering: h(t), eq 8, shown in the top panels in black is Fourier transformed (f H(ω) black bottom row); H(ω) is then low pass filtered by means of the red F(ω), eq 4, with s1 = 0.25 and ω0,1 = 12% 3 ωmax (ωmax = 512 Hz), which preserved the 5 Hz and the 20 Hz components ~ but not the 100 Hz. The iFT of H(ω) resulted then in the red signal (top row panels).

Figure 5. Fourier high pass filtering: function, eq 8 has been Fourier transformed (black, left panel), high pass filtered, eq 5 (gray, left panel) with s2 = 0.25 and ω0,2 = 12% 3 ωmax and inverse Fourier transformed (black, right panel)—now the low frequency components at 5 Hz and 20 Hz have been removed and only the high-frequency components at 10 0 Hz remains.

coincidental trends in the noise. On the other hand, for very narrow signal features (see right end of signal), a higher-order polynomial

introduces less distortions because it can better follow sharp signal changes than a lower order. In general, the user should not choose 1675

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 6. Fourier band-pass filtering, eq 6, with s1 = s2 = 2, ω0,1 = 5% 3 ωmax, ω0,2 = 2% 3 ωmax; only the medium frequencies (here: 20 Hz) have been retained and all others (here: 5 Hz and 100 Hz) were suppressed.

Figure 7. Fourier band block filtering, eq 7 with s1 = s2 = 2, ω0,1 = 5% 3 ωmax, ω0,2 = 2% 3 ωmax removes medium frequencies (here: 20 Hz) and retains all others (here: 5 and 100 Hz).

k too large (i.e., e5) or the polynomial describes too much of the signal artifacts and thus defeats the purpose of noise reduction. On a side note, the Savitzky-Golay algorithm is also commonly applied to compute derivative signals that help, for instance, to reduce baseline drifts in spectroscopy or to emphasize small shoulders in signals911,20 (see the Supporting Information for examples). Filtering Methods Based on Fourier Transforms

Fourier filtering comprises three steps: (i) A time-dependent measurement signal h(t) is Fourier transformed (see the

Supporting Information for an overview) into the frequency domain and represented by H(ω). (ii) This frequencydependent function H(ω) is then weighed (eq 3) with the filter function F(ω), which suppresses certain, user-selectable frequencies and preserves others. Frequencies ω for which F(ω) = 0 are removed from the signal and frequencies for which F(ω) = 1 remain unchanged; 0 < F(ω) < 1 reduces these frequencies contributions to H(ω). (iii) The filtered ~ function in Fourier domain, H(ω), is then inverse Fourier transformed back into the original time domain resulting in 1676

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 8. Results on noise level #3 (cf. the Supporting Information) from different Fourier filters; (top left) denoising via low pass filtering with s1 = 1, ω0,1 = 25% 3 ωmax; (top right) ditto but with s1 = 1, ω0,1 = 5% 3 ωmax; (middle left) example for high pass filtering with s2 = 1, ω0,2 = 1% 3 ωmax as an example for suppression of broad signal features; the bending toward negative values and the oscillations on the right end are due to removal of frequencies that would have compensated this resulting in a flat baseline; (middle right) impact of low pass filtering on noise-free data (noise level #0, s1 = 1, ω0,1 = 5% 3 ωmax) demonstrating the removal of real signal components; (bottom left) band-pass with s1 = s2 = 1, ω0,1 = 15% 3 ωmax, and ω0,2 = 1% 3 ωmax; (bottom right) band block with same parameters as band-pass. Note the impact of the filters on the Gaussians of different widths.

the filtered function ~h(t). filtering

FT iFT ~ hðtÞ f HðωÞ f HðωÞ ¼ FðωÞ 3 HðωÞ f ~hðtÞ

ð3Þ

A low pass filter is typically applied to remove high-frequency noise. A high pass filter suppresses slow signal drifts assuming that

the dca and low-frequency signal components do not contain relevant information. A band-pass filter is applied when only signal components at a specific frequency (range) contain relevant information; all others can or should be suppressed. In photoacoustic spectroscopy, for instance, a chopper generates a concentration dependent ac signal at a well-defined modulation 1677

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education frequency; other frequencies cannot be related to chemical information. A band block (also: band stop or notch filter) quenches a disturbance of known frequency such as the 60 Hz (or 50 Hz) line frequency. Filter design is a wide field and in this introduction only the basic ideas can be covered without being rigorous from a mathematical or engineering perspective.4,5,9 Filters utilized in this tutorial (see eqs 47 and Figure 3) incorporate two parameters that “come as a package”: The first is a steepness s, which describes in how many frequency units the filter transmission F(ω) changes from F(ω) = 0 and F(ω) = 1 are. The other, ω0, defines the center frequency ω of the filter at which F(ω = ω0) = 0.5. The maximum frequency ωmax contained in a signal is defined by the data acquisition processb and in order to state ω0 independent from a specific signal, it will be given in percent of ωmax: ω0 = X% 3 ωmax with ωmax; the accompanying software DENOISE.exe makes this conversion from hertz to percent automatically. In the remainder, a subscript 1 on filter parameters refers to a low pass filter; the subscript 2 refers to a high pass filter.

ARTICLE

Low pass filter in frequency domain (Figure 3, top left): FðωÞ ¼ 1 

1 1 þ expfs1 ðω  ω0, 1 Þg

ð4Þ

High pass filter (= one minus low pass, Figure 3, top right): FðωÞ ¼

1 1 þ expfs2 ðω  ω0, 2 Þg

ð5Þ

Band pass filter, that is, a product of low pass, eq 4, times high pass, eq 5, (Figure 3, bottom left): ! FðωÞ ¼

1 1 1 þ expfs1 ðω  ω0, 1 Þg

! 1 1 þ expfs2 ðω  ω0, 2 Þg

ð6Þ Band block filter, that is, an inverted band-pass (= one minus band-pass, eq 6 (Figure 3, bottom right)): FðωÞ ¼

! ! 1 1 1 1 1 þ expfs1 ðω  ω0, 1 Þg 1 þ expfs2 ðω  ω0, 2 Þg

ð7Þ

Figure 9. Band pass filtering for removal of very low-frequency baseline drifts and high-frequency noise from a Raman spectrum22 with s1 = s2 = 1, ω0,1 = 15% 3 ωmax, ω0,2 = 0.075% 3 ωmax .

Note for a band-pass as defined in eq 6, the transmission range of the low pass, eq 4, and the high pass, eq 5, must overlap or otherwise all frequencies would be suppressed and thus ω0,1 > ω0,2 must be chosen (especially for steep flanks s1,2 g 1). The same requirement applies to band blocks or all frequencies would be transmitted. Most signals in the Fourier domain are complex-valued functions while the filters F(ω), eqs 47, used in this article are real-valued. Thus, the same F(ω) is applied to the real and the ~ imaginary part of the Fourier transformed signal: H(ω) = F(ω) 3 ReH(ω) + i 3 F(ω) 3 ImH(ω). Because the real and the imaginary part of H(ω) are multiplied with the same number, only the amplitudes r(ω) are changed but not the phases ϕ(ω) (see the Supporting Information for more details). For demonstration purposes (Figure 47), the following function has

Figure 10. Wavelet noise removal performed on the gray signal in the left panel; after applying the threshold, eq 10, for relevant coefficients cτ,s (right panel: gray dots f black squares), an iWT derives the smoothed black signal (left). 1678

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 11. Impact of denoising threshold th, eq 10, and wavelet type on the filtering results (gray original, black smoothed signal): (top left) Daub12 and th = 5%; this example demonstrates that choosing the threshold to high removes too many wavelet coefficients and causes signal artifacts with shapes similar to the used wavelet type; (top right) Daub12 and th = 1% is an example for a well-balanced filtering procedure as all relevant signal features are preserved and only a little noise remains; (bottom left) Daub12 and th = 0.5% is a typical borderline case when the threshold was too low as apparently some wavelet coefficients are just below and others just above the threshold, eq 10, resulting in this partially smooth and partially noisy signal; (bottom right) the particular unsmooth nature of the Daub2 (= Haar) wavelet (th = 1%) introduces steps in the filtered signal due to omission of small wavelet coefficients that would partially have filled these steps; such situations are best avoided by selecting a smoother wavelet type.

been used: hðtÞ ¼ 7:0 3 3 3 sinð2π 3 5 Hz 3 tÞ þ 10:0 3 3 3 sinð2π 3 20 Hz 3 tÞ þ 0:5 3 sinð2π 3 3 3 100 Hz 3 tÞ

ð8Þ

The impact of ω0,1 in a low pass is demonstrated in Figure 8 top row. The further the cutoff frequency is shifted toward lower frequencies, the lower is the resulting noise level; however, real signal components that are high-frequency such as the needle-shaped spikes at the right end of the signal are canceled as well. In such situations, selecting s1 , 1 might help to reduce detrimental impacts as the filter function can then reach into the high frequency regime. However, low pass filtering will usually require a compromise between noise suppression and blurring of relevant but high-frequency components. This is further demonstrated by low pass filtering a noise free signal (Figure 8, middle right); high frequencies are suppressed resulting in reduced peak heights for sharp features. Pointing this out to students is particularly important as it clearly shows that high frequencies are not pure noise, a common misconception!

Figure 8 (middle left) demonstrates the outcome of high pass filtering and shows that the needle-shaped features are in the same frequency region as what apparently is mostly noise in other parts of the signal. The bottom row of this graph demonstrates the impact of a band-pass that removes the broadest and narrowest peaks but (partially) transmits features of medium width. The band block suppresses features of medium frequency composition. Figure 9 presents a real-world application of a band-pass filter to Raman spectrum that is impacted by an unwanted fluorescent background. As such fluorescence is often rather broad compared to the Raman bands of interest, a band-pass has been used to suppress the unwanted dc component and the lowest frequencies as well as high-frequency noise.22 SNR Enhancement Methods Based on Wavelet Transforms

Applications of Fourier transforms (FTs) and Wavelet transforms (WTs)9,1624 in signal analyses have a common goal, namely, preserving the wanted signal frequencies and suppressing unwanted components. However, whereas a FT decomposes a 1679

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 12. Comparing an original color image (top left) and its 2D Fourier low pass filtered counterpart (top right, s1,x = s1,y = 0.1, ω0,1,x = ω0,1,y = 10% 3 ωmax) depicts the blurring effect; for such color images, filters are applied to all three color channels separately. (Bottom left) Microscope picture of microalgae cells dispensed onto a hemocytometer whose grid lines define areas/volume for cell counting; (bottom right) after 2D Fourier high pass filtering (s2,x = s2,y = 1, ω0,2,x = ω0,2,y = 10% 3 ωmax) only retains the high-frequency components at the cells’ outline when the signal levels considerably changes within a very short distance.

signal h(t) into sines and cosines, a WT utilizes “wavelets” Ψτ,s(t) instead. Resulting from this, the major difference between FT and WT is their time-frequency localization: A FT analyzes what frequencies are contained in a signal without specifying at what point of time they occur. In contrast to that, a WT determines in which signal windows (=x axis region) what frequency bands are present. The latter can be achieved because the wavelets, unlike “true waves”, are unequal to zero only in a restricted time window (the wavelet's “support”). Thus, only the portion of the signal h(t) inside this window contributes to the corresponding wavelet coefficientc:   Z 1 ∞ tτ hðtÞ 3 ψ dt cτ, s ¼ pffiffiffi a ∞ s Z ∞ ¼ hðtÞ 3 ψτ, s ðtÞ dt ð9Þ ∞

Aside from a normalization constant a, a wavelet Ψτ,s(t) contains two parameters τ and s, which are varied during a WT: τ shifts the wavelet along the t-axis over the signal h(t) and the “scale” parameterd s defines the window width in which it is unequal to zero (= its “support”). Translating a wavelet ψ[(tτ)/s] via τ and dilating it via s is best understood by a simple example: Assume a function 8 < 1 for : 0 e x e 1 jðxÞ ¼ : 0 otherwise and determine for which x is j(x2) = 1? For 2 e x e 3, because for these x-values the function’s argument x  2 falls into the range 0 e x e 1 (translation). For which values of x is j (x/2)=1 ? For: 0 e x e 2

(scaling); and finally, for which values of x is j [(x/4) +3] = 1? For: 12 e x e 8 (translation and scaling). A WT starts with the narrowest wavelet and scans it over the entire signal and at each position τ a wavelet coefficient cτ,s, eq 9 is computed. Then, the wavelet is made wider by adapting s and the scanning procedure is performed again; this procedure is repeated multiple times. As demonstrated in the Supporting Information, the impact of a wavelet on a signal is best understood when considering the signal in the frequency domain. Wavelets of increasing widths are bandpass filters with transmission functions shifted more and more toward dc. The wavelets’ shift in time translates into a phase shift in the Fourier domain. These band passes with different phase shifts establish the aforementioned time-frequency analysis of the signal. Wavelet Noise Removal12,15,18 consists of three steps that follow an idea very similar to Fourier filtering, eq 3: (i) wavelet transform a signal; (ii) replace small wavelet coefficients with zeros; (iii) inverse wavelet transform the modified signal (Figure 10). For canceling small wavelet coefficients, the user has to define a threshold th discriminating relevant from irrelevant wavelet coefficients (Figure 10): jcτ, s j < th w cτ, s ¼ 0

ð10Þ

In DENOISE.exe, th = X% max|cτ,s| has been defined with a user-selectable parameter X denoting how many percent of the largest absolute wavelet coefficient defines the cutoff.X% = 1% has being found empirically to be a good starting point for application-dependent adaptations. Furthermore, the user also has to select a wavelet type for denoising, which also has a considerable impact on the filter outcome (cf. Figure 11). Although DENOISE.exe has been restricted to wavelets from the 1680

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

Figure 13. Impact of threshold th, eq 10, and 2D-wavelet type on the 2D WT image filtering; the dark, grid-shaped light pattern was generated for an unrelated purpose23 but helps to assess filtering results; (top left) original, grainy image; (top right) Daub12 with th = 1% introduced considerable blurring possibly due to a too high threshold; (bottom left) Daub12 with th = 0.1% has removed much of the image’s graininess without introducing blurring; (bottom right) Daub2 (= Haar) with th = 1%; especially for such a large th, this wavelet type introduces a boxlike structure in the wavelet-filtered image similar to the results shown in Figure 11 for 1D-signals (bottom, right) and thus should be avoided.

Daubechies family,9,17 namely, Daubechies 2 (Daub2), ..., Daubechies 20 (Daub20), many more can be applied accordingly. Daub8 and Daub12 have been found in several tests to be good general-purpose wavelets; Daub2, on the other hand, is prone to introduce unwanted, step-shaped artifacts. The advantage of wavelet denoising over Fourier filtering is that the user does not need to set an artificial cutoff frequency above which no relevant signal components are expected and thus all features are suppressed by the filter. Instead, the size of a wavelet coefficient decides whether this particular timefrequency component is important. This is a crucial advantage leading to less signal blurring compared to FTs as all relevant signal components are preserved and all irrelevant ones are removed independent from their frequency. Figure 11 presents some selected examples for wavelet denoising. Expanding Data Filtering to Images and Data Cubes

Data filtering method for 1D-signal can be expanded in a straightforward way toward 2D-data sets (matrices or images) and 3D-data cubes. Averaging methods incorporate a certain area of pixels; multidimensional FTs and WTs consecutively apply 1D-FTs and 1D-WTs, respectively, in different directions. Averaging Based Methods. For images and data cubes, a common method for noise reduction is pixel binning, which averages neighboring pixels (2D) or voxels (3D) and replaces an original pixel (voxel) with its neighbors’ average. The drawback of this approach is a reduced spatial resolution as fewer pixels cover the same field of view; averaging voxels in spectroscopic imaging additionally reduces wavelength resolution. Furthermore, in color images, color distortions can be introduced if neighboring pixels of different colors are averaged

Expanding Fourier filtering from 1D-signals to multidimensional data sets is straightforward as the FTs in different directions are independent from each other and can be performed consecutively. A matrix for instance is 2D transformed by consecutively performing 1D FFTs on each row of the matrix followed by running 1D FFTs on each column of the row-transformed matrix. A 3D FFT expands this concept by performing 2D FFTs for each XY slice of a data cube followed by transforming each vector in Z direction of the XY 2D transformed data cube. The multidimensional inverse FFTs perform the same type of operations based on 1D iFFTs. Because a large number of 1D FFTs has to be performed, N-dimensional FFT algorithms becomes time-consuming and data handling was speed optimized.9 Fourier filters for multiple dimensions need to have one frequency axis per dimension. Such N-dimensional filters are then multiplied element-wise with the N-dimensional Fourier transformed data sets. Images (matrices) for instance require filter matrices with the x axis being the frequency axis containing information on how fast the signal level changes in horizontal direction (filter parameters s1,x and s2,x as well as ω0,1,x and ω0,2,x) and the y axis comprising information of frequency components in vertical direction (filter parameters s1,y and s2,y as well as ω0,1,y and ω0,2,y). The two filter functions in x  and y  direction are then multiplied. To visualize the impact of 2D Fourier filtering on image data, Figure 12 (top row) compares a digital photo to its low pass filtered version. As this image contains homogeneous areas such as the sky and heterogeneous ones (bushes), the impact of different Fourier filters can be studied. Low-pass filtering of rows and columns suppresses fast signal changes, that is, changes of color or light intensity within a few pixels. This results in image details such as branches of the bushes being “washed out” and the image becomes blurred. Figure 12 (bottom row) shows an 1681

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education example of high-pass filtering. This microscope image was acquired in transmission mode microalgae cells which are contained as black dots as they block the light; the straight lines are borders of hemocytometer wells into which cells dispensed in an aqueous medium have been placed. Areas without cells feature a homogeneous, high light level. After high-pass filtering, large negative signal values (cf. Figure 8, middle left) are encoded into dark pixels and high signal levels in bright pixels; medium gray pixels represent signal intensities close to zero. Thus, only at image locations where the original signal changed drastically, a nonzero filtered signal was obtained. The aforementioned hemocytometer’s grid lines equally blocked the light; however, they did not change in either x or y direction (≈ dc) and thus the high pass canceled them. This is an example of how filtering can suppress unwanted information that may not necessarily be random noise. Applications of high-pass filtering might be suppression of constant, noninformative signal contributions, which only limit the resolution of the information carrying signal components. Multidimensional WTs24,25 sequentially combine 1D WTs in all dimensions of the data set. Images, for instance, are 2D wavelet transformed by applying a 1D WT to each row of the matrix followed by performing 1D WTs of the columns of the row-transformed matrix. Then, noise suppression is performed in an equivalent way as described for 1D signals in eq 10 and related discussion. The interpretation of 2D WT can be found in ref 9 but this is not relevant to wavelet denoising as the wavelet coefficients are simply checked whether they are large enough to be required (eq 10). The advantages of WTs over FTs also apply to multidimensional transforms. In addition, WTs require less floating point operations and thus can be performed faster, an advantage that becomes more important the larger the dimensionality of the data sets is. Figure 13 gives three examples for 2D WT noise filtering applied to image processing in microscopy.

’ CONCLUSIONS Enhancing the signal-to-noise ratio is of central importance in analytical chemistry and improvements in experimental setups are not always sufficient or feasible. For such cases, three computational methods for noise reduction were presented that operate on digitized data sets. Signal averaging based on the Savitzky-Golay method has been compared to Fourier and wavelet filtering; their advantages and disadvantages were discussed and compared. Algorithms are introduced and assessed first by means of 1D signals such as typically encountered in spectroscopy, chromatography, and so forth. In a subsequent step, these computation tools were then expanded and applied to 2D data sets like images and finally to 3D data cubes. The learning goal is to utilize the mathematical knowledge upper-level undergraduates and first-year graduate students have to advance real-world applications in analytical chemistry. Another aim is to excite students through a synergetic combination of two fields, applied math and instrumental analyses, they might have perceived as unrelated. To facilitate hands-on teaching including immediate visualization of these data filtering techniques, a Windows-based program and test data were provided. ’ ASSOCIATED CONTENT

bS

Supporting Information Suggestions on introducing Fourier and wavelet transforms in lectures; a windows-based software package, DENOISE.exe, to support lectures by visual demonstrations and hands-on learning;

ARTICLE

the simulated and real-world data used in this article. This material is available via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ADDITIONAL NOTE a In the remainder, dc (direct current, “current” although this applies to all kind of signals) refers to signal (components) that are constant over time; ac (alternating current) refers to time dependent signal (components) in the widest meaning. b

The maximum frequency is determined by the band limitation of the amplifier electronics and the analog-to-digital conversion frequency. c

Wavelet coefficients play the same role as the amplitudes of a certain sine or cosine in Fourier transforms. This becomes clear when comparing eq 9) to the definition of the FT (cf. the Supporting Information).

d

This scale parameter s is not related to the steepness s used in Fourier filters, eqs 47.

’ REFERENCES (1) Skoog, D., Holler, F., Crouch, S. Principles of Instrumental Analysis, 6th ed.; Thomson, Brooks, Cole: Belmont, 2007; Chapter 5. (2) Hayes, M. Digital Signal Processing; Schaum’s Outline Series, McGraw-Hill: New York, 1999. (3) Fountain, A. W. J. Chem. Educ. 2001, 78, 271. (4) Strang, G, Nguyen, T. Wavelets and Filter Banks; WellesleyCambridge Press: Wellesley, MA, 1997. (5) Bremaud, P. Mathematical Principles of Signal Processing—Fourier and Wavelet Analysis; Springer: New York, 2002. (6) Sen, A.; Srivastava, M. Regression Analysis—Theory, Methods, Applications; Springer: New York, 1990. (7) Draper, N.; Smith, H. Applied Regression Analysis, 3rd ed.; Wiley: New York, 1998. (8) Weisberg, S. Applied Linear Regression, 3rd ed.; John Wiley & Sons: Hoboken, 2005. (9) Press, W, Teukolsky, S, Vetterling, W, Flannery, B. Numerical Recipes in C++, 2nd ed.; Cambridge University Press: New York, 2002. (10) Savitzky, A.; Golay, M. Anal. Chem. 1964, 36, 1627. (11) Steiner, J.; Termonia, Y.; Deltour, J. Anal. Chem. 1972, 44, 1906. (12) Walczak, B; Massart, D. Trends Anal. Chem. 1997, 16, 451–463. (13) Barclay, V.; Bonner, R.; Hamilton, I. Anal. Chem. 1997, 69, 78–90. (14) Walczak, B. Wavelets in Chemistry: Elsevier Science: New York, 2000. (15) Ehrentreich, F. Anal. Bioanal. Chem. 2002, 372, 115–121. (16) Mallat, S. IEEE Trans. Pattern Anal. Machine Intell. 1989, 11, 674–693. (17) Daubechies, I. Ten Lectures on Wavelets; CBMS-NSF regional conference series in applied mathematics (61): Philadelphia, PA, 1992. (18) Donoho, D. IEEE Trans. Inf. Theory 1995, 41 (3), 613–627. (19) Cai, C.; Harrington, P. J. Chem. Inf. Comput. Sci. 1998, 38, 1161–1170. (20) Vogt, F. Spectrophotometry: Derivative Techniques. In Encyclopedia of Analytical Science, 2nd ed.; Worsfold, P.; Townshend, A.; Poole, C., Eds.; Elsevier: Oxford, 2005; Vol. 8, pp 335343. 1682

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683

Journal of Chemical Education

ARTICLE

(21) Strang, G. Linear Algebra and its Applications, 3rd. ed.; Harcourt Brace Jonanovich College Publishers: Fort Worth, 1988. (22) Bu, D.; Huffman, S.; Seelenbinder, J.; Brown, C. Appl. Spectrosc. 2005, 59, 575. (23) Gilbert, M.; Vogt, F. Anal. Chem. 2007, 79, 5424–5428. (24) Vogt, F.; Booksh, K. J. Chemom. 2005, 19, 575–581. (25) Vogt, F.; Cramer, J.; Booksh, K. J. Chemom. 2005, 19, 510–520.

1683

dx.doi.org/10.1021/ed100984c |J. Chem. Educ. 2011, 88, 1672–1683