Wavelet analysis of gamma-ray spectra

Since September 11, 2001, there has been increasing interest in providing first responders with radiation detectors for use in the search for and isotope identification of potentially-smuggled special nuclear material (SNM) or radiological dispersal devices (RDDs). These devices are typically comprised of low-resolution detectors such as NaI, thus limiting their identification abilities. We present a new technique of wavelet analysis of low-resolution spectra for the use in isotope identification. Wavelet analysis has the benefit of excellent feature localization while, unlike with Fourier analysis, maintaining the signal frequency and time characteristics. We will demonstrate this technique with a series of gamma-ray spectra obtained from typical hand-held isotope identifiers, illustrate figures-of-merit to be applied to these results, and discuss future algorithm optimization.


I. Introduction
It has been reported by several sources that the handheld isotope identifiers currently being used by firstresponders have problems with their identification abilities. [1], [2] It can be difficult to implement complicated peak localization and identification algorithms due to limitations of the on-board computer memory, thus limiting their effectiveness in the field of nuclear emergency response. These algorithms are further hindered by problems such as gain drift, which is very common when the devices are moved between temperature extremes such as from an airconditioned vehicle to the warm outdoors. Furthermore, these algorithms become more complicated as the detector resolution increases making the identification of multiplelined sources difficult. A new algorithm is needed that requires minimal computer memory to be used to locate and identify peaks within a medium-or low-resolution gamma-ray spectrometer.
Wavelet analysis was presented in the mid-1980's to solve problems such as image compression, gamma-ray burst analysis, the study of earthquakes, and numerous other problems in signal analysis where the signal is aperiodic, noisy, transient, and so on. More recently, wavelets have been applied to problems such as feature detection and localization in both one-and two-dimensional signals. The concept of wavelet analysis has been applied by other researchers to the analysis of high-resolution gamma-ray spectra where the initial spectrum was denoised using the discrete wavelet transform (DWT) and then the peaks were located by subtracting a fit of the continuum from the signal. [3], [4] We present a new approach where the features of a spectrum are identified directly with a wavelet transform using the modulus maxima technique. [5], [6], [7] This eliminates the requirements of continuum fitting for the detection of peaks and takes full advantage of the capabilities of wavelet analysis without losing any data due to smoothing.

II. Wavelet Transforms and the Modulus Maxima
The technique of wavelet analysis involves observing a signal at a particular time (or, in this case, energy) simultaneously over a broad range of scales, (the wavelet pseudo-equivalent of frequency). Like Fourier analysis, the signal is broken down into frequency components. However, standard Fourier analysis suffers from the fact that it does not maintain the information about the location in a signal of the given frequency. The windowed Fourier transform, or Gabor transform, provides a basic frequency analysis as a function of position. However, the Gabor transform is of fixed time-frequency resolution. In the case of gamma-ray spectroscopy, to analyze a signal over a fixed number of channels across the entire spectrum would be prohibitive since, for most systems, the width of a peak in channels increases with increasing channel number.
Wavelet analysis permits the user to analyze the signal in the time-frequency (or energy-scale) domain simultaneously at different scales, which is referred to as "multiresolution analysis" (MRA). The wavelet transform of a signal x is given by where E is position (in our case, energy), s is the scale, and ψ E,s (t) is the mother wavelet. By adjusting both the scale and the location of the mother wavelet during the convolution, the transform coefficient, T , can be obtained that indicates the degree to which the mother wavelet matches the original signal. The result of the wavelet transform is then a three-dimensional array of position, scale, and wavelet coefficient, most effectively shown as an image called a scalogram where the shading indicates the magnitude of the wavelet coefficients, T (E, s). The exact details of the wavelet transform have been covered extensively by authors elsewhere. [5], [6] There are numerous choices for mother wavelets and the selection of the proper wavelet depends on the application. To be considered a wavelet, a function ψ ∈ L 2 must have finite energy: Also, ifψ(f ) is the Fourier transform of ψ(t), then which implies that the wavelet must have zero mean. ψ is then scaled by s and translated by E as: Thus the wavelet transform can be rewritten as: Several example wavelets are shown in Figure 1 and a sample scalogram of an input signal using the Mexican hat wavelet is illustrated in Figure 2.
As can be observed in Figure 2, the wavelet coefficients in the scalogram create a cone-like appearance whose boundaries converge at small scales on the location of an abrupt change such as a singularity or peak. This position, called the cone of influence, which is defined by: where E 0 is the location of the feature and C is a constant. In order to determine the precise location of E 0 , we employ the wavelet transform modulus maxima technique (WTMM). The WTMM is used to describe any point in the scalogram such that the wavelet coefficients are a local maximum: A modulus maxima line is any line that connects neighboring maxima points in the scalogram. Singularities or peaks   are detected by finding the abscissa where the modulus maxima lines converge at fine scales. Therefore, the wavelet transform can focus on a localized singularity or similar feature with a zooming procedure that gradually reduces the scale.
However, to actually detect a singularity, it is not sufficient to simply follow the modulus maxima lines since small fluctuations such as noise in the signal could create a local maximum, especially at small scales. The singularity itself is characterized by the decay of the wavelet coefficients along maxima lines. This decay quantifies how regular the function is either over all space or at a particular location. One measure of the regularity of a signal at position E 0 is the value of the Lipschitz exponent (also referred to as the Hölder exponent) at that position. If ψ is continuously differentiable with compact support and real values, then for any > 0 , x is uniformly Lipschitz α if and only if there exists a constant, A , such that for where α is the Lipschitz exponent. This is a condition on the asymptotic decay of the wavelet coefficients as the scale approaches zero. If x has a singularity at E 0 , then the Lipschitz exponent characterizes the singular behavior. One must determine α for each line, which is generally calculated by fitting: to the magnitude of the wavelet coefficient along a modulus maxima line. If the slope of this line is calculated at E 0 , this determines the Lipschitz regularity. This is shown in Figure 3 for the cusp in the test signal of Figure 2. In this case, the singularity illustrated was Lipschitz 0.5. Note that large values of α indicate a more regular function whereas smaller values indicate more singular function. An actual singularity is found when 0 < α < 1.

III. Results on Gamma-Ray Spectra
This paper presents a proof-of-principles experiment using the WTMM and Lipschitz exponent to locate peaks within a spectrum. The data were taken either with a 15 x 15 x 7.5 mm 3 CdZnTe hand-held isotope identifier or a 3 in diameter x 3 in high NaI scintillator. The computations were performed using MATLAB (version 6) with the MATLAB Wavelet Toolbox (version 2) and the Wavelab Toolbox (version 802). [8], [9] The first step was to determine the wavelet transform of the spectrum, using a continuous, real wavelet transform. It may seem counterintuitive to have data that is sampled discretely analyzed with a continuous filter. However, to be clear, what is meant by "continuous" in this case involves the set of scales and positions on which the transform operates. The continuous wavelet transform (CWT) can operate at every scale, from that of the original signal to that of the user's choosing. It also is continuous in terms of shifting the mother wavelet smoothly over the entire signal. The cost in having a very large number of scales or shifted positions is a significant increase in the time required to complete the transform.
In order to calculate the CWT, it was first necessary to select a mother wavelet for the transform. It was observed that the choice of mother wavelet had a significant impact on the WTMM process in terms of the numbers and locations of modulus maxima lines detected. This was not unexpected since different wavelets have different frequency responses and thus different capabilities in the localization of peaks within a signal. Previous work used the Haar wavelet for analysis. [3] This particular wavelet is good at identifying rapid (high-frequency) changes in the signal, such as a those from a high-resolution peak in an HPGe spectrum. However, our goal was to locate peaks within medium-to low-resolution spectra. Therefore, we choose to experiment with two different wavelets that resembled more of the features we were looking for. We decided to use wavelets that resembled the second derivative of a Gaussian function. This decision was made to maintain similarity with other isotope identification algorithms that are based on Mariscotti's technique, which analyzes the second derivative of the signal within a small window or kernel. [10] Unlike Mariscotti's technique, the wavelet approach benefits from the use of MRA, thus permitting the analysis of peaks at different scales instead of being restricted to a constant kernel size that may not be optimal for all peaks. The wavelets chosen for this work both resembled the "Mexican Hat Wavelet" of Figure 1, which is the second derivative of a Gaussian. We used the Wavelab "Sombrero" wavelet: and the "DerGauss" wavelet: which is the second derivative of a Gaussian. Once the CWT of the spectrum was calculated, the modulus maxima lines were found. This algorithm was first performed on a 137 Cs spectrum taken with a CdZnTe detector. The results are shown in Figure 4. As can be seen from the scalograms, the two different wavelets result in very different wavelet coefficients. This does not matter for the subsequent determination of the modulus maxima positions since the algorithm, based on Equation 7, looks only for the places where the first derivative is zero. In reality, this could be either a local maximum or minimum -no further test is applied to determine which it is. What is important is that there are modulus maxima lines that are converging on the location of the 662 keV peak. There are differing numbers of lines converging on 662 keV due to the selection of the wavelet.
There are several other features to notice in these figures. First, in both cases there is a maxima line that crosses the entire scalogram. This is due to edge effects in the scalogram and can be reduced or eliminated by considering a more appropriate range of scales within the scalogram. These particular figures were plotted with a large number of scales to demonstrate the behavior of the CWT at very large scale. Additionally, several maxima lines are observable, particularly at smaller scales, that are not associated with the photopeak. These are due to the noise within the signal and other features such as the Compton edge and backscatter peak within the spectrum.
While the scalogram and plots of the modulus maxima lines yield interesting results on the potential peak locations within the spectrum, there are clearly places where these lines do not represent data of interest, such as the aforementioned noise and other spectral features. As such, it would be desirable to eliminate some of these modulus maxima  lines that do not correspond to a true peak. To create such a filter, the Lipschitz exponent of the various spectral features was calculated for this spectrum. Table I shows the fit slope (i.e. the Lipschitz exponent) of the wavelet coefficient decay along maxima lines for lines associated with different types of spectral features. As is evident from the table, maxima lines associated with a true peaks had significantly larger Lipschitz exponents than those that were associated with the Compton continuum or noise. Based on this observation, it was possible to establish a filter for the modulus maxima lines based on Lipschitz exponent where peaks within the spectrum were identified if they met the requirement of α > 0.65. With this threshold, it is clear that in the cases of both wavelets, the maxima lines corresponding to the photopeak are preserved. In the case of the Sombrero wavelet, some of the edge effect lines will also be preserved as will the backscatter peak. However, this is not true for the DerGauss wavelet, which maintains only the two photopeak maxima lines and the edge effect line.
This filtering technique was then applied to a more complicated spectrum, which was a mixture of 137 Cs and 133 Ba. These results are shown in Figures 5 and 6 with the remaining modulus maxima lines indicating which lines had met the Lipschitz slope filtering requirement. The results of these figures show variations in the sensitivity of the two wavelets. For example, the Sombrero wavelet is capable of detecting the 302, 356, and 383 keV photopeaks of 133 Ba, but it also is more sensitive to the pseudo-peak where the spectrum begins around channel number 25. The DerGauss wavelet, on the other hand, is not as sensitive to this portion, but it fails to detect the 302 keV peak. It also detects the Compton edge from the 662 keV 137 Cs photopeak. This suggests that more work is required on the proper wavelet selection for this application.
In addition to the calculations performed on the CdZnTe spectra, this technique was also applied to NaI spectra to determine its abilities with low-resolution detectors. Due to the lower resolution of these systems, it was expected that α should be larger to indicate that the signal was approaching a more regular function. Table II shows how α varied as a function of spectrum feature. Based on these data, a slope filter was imposed of α > 1.4. Figures 7 and  8 show the results of analyzing a NaI spectrum of 137 Cs using this filter. As is evident from the table and figures, the edge effects were detected due to the fact that they had a significant slope. It is worth considering whether an upper limit on slope should be imposed. The impact of this additional filter will be explored in future work.

IV. Summary and Future Work
We have demonstrated a promising new technique involving the localization of peaks within a gamma-ray spectrum of low-to medium-resolution using the concepts of the wavelet transform with the modulus maxima technique. The modulus maxima lines can be used as a first-pass indication on the potential location of a peak. It is clear from the results that the Lipschitz exponent provides useful information on the nature of the modulus maxima lines to determine whether or not they actually are a peak. It was observed that there was a significant difference between the values of the Lipschitz exponent between different spectral features. It was shown that a filter could be established by requiring a minimum value of the slope. Using this filter, we were able to discriminate between photopeaks, Compton edges, and noise within the spectrum. In some cases, we were also able to distinguish the difference between a photopeak and a backscatter peak, although more work is required in this area.
There are still important areas to consider with this technique. First, while the wavelets chosen for this analysis were reasonably mathematically similar (both approximated the second derivative of a Gaussian), their results were rather different. This suggests that a detailed comparison of the results using different wavelets and what parameters (detector resolution, photopeak shape, etc.) are important for optimal wavelet selection are necessary. Second, in keeping with the eventual goal, a fitting algorithm needs to be employed to determine precisely where the modulus maxima lines corresponding to a peak intersect the abscissa. This will provide the means to identify the energy of the photopeak for comparison with a library to be used in isotope identification. Finally, it is desirable to use wavelets or some other technique to fit the continuum under each peak to provide estimates of the net area for each peak, which could be used to obtain crude isotopic ratio information.