Methods for Averaging Spectral Line Data

The ideal spectral averaging method depends on one’s science goals and the available information about one’s data. Including low-quality data in the average can decrease the signal-to-noise ratio (S/N), which may necessitate an optimization method or a consideration of different weighting schemes. Here, we explore a variety of spectral averaging methods. We investigate the use of three weighting schemes during averaging: weighting by the signal divided by the variance (“intensity-noise weighting”), weighting by the inverse of the variance (“noise weighting”), and uniform weighting. Whereas for intensity-noise weighting the S/N is maximized when all spectra are averaged, for noise and uniform weighting we find that averaging the 35%–45% of spectra with the highest S/N results in the highest S/N average spectrum. With this intensity cutoff, the average spectrum with noise or uniform weighting has ∼95% of the intensity of the spectrum created from intensity-noise weighting. We apply our spectral averaging methods to GBT Diffuse Ionized Gas hydrogen radio recombination line data to determine the ionic abundance ratio, y +, and discuss future applications of the methodology.


INTRODUCTION
Averaging spectral line data allows one to increase the signal-to-noise ratio (SNR) of the resultant spectrum.Such averaging is straightforward when the spectra are taken of the same source and have similar noise characteristics.The situation is complicated, however, if the noise or source intensity differ significantly between observations.In such cases, depending on the distribution of peak intensities in the observations and the weighting scheme, averaging spectra can result in a decrease in the of 2-D images that better accounts for noisy data where the SNR of individual observations is difficult to measure.Adaptive smoothing of 2D images, such as using Voronoi Tessellations or Weighted Voronoi Tessellations, can be used to create spatial regions that meet user-specified SNR criteria (Cappellari & Copin 2003;Diehl & Statler 2006).The ideal averaging method may depend on whether the intensities of the spectra to be averaged are uniform or have a large variance.
In this paper, we explore methods for spectral averaging and provide guidance for multiple use-cases.We focus our analytical treatment on radio spectroscopic observations (i.e., we use the variable "T " for intensity), but the method is applicable to any spectral line data set.

THE SIGNAL TO NOISE RATIO
A spectral line has a SNR given by (Lenz & Ayres 1992) SNR = C ∆V ∆λ where C is a constant whose value depends on the line shape, ∆V is the full width at half-maximum (FWHM) line width, ∆λ is the spectral resolution (or width of the smoothing kernel), ⟨T P ⟩ is the average peak line brightness temperature, and ⟨σ⟩ is the average rms noise.For a spectrum that can be modeled as a Gaussian line with white noise, C = 0.7.
Upon averaging n spectra each with (unnormalized) weighting w i , the average intensity at a given spectral channel is (2) At the line center, the peak line intensity ⟨T P ⟩ is therefore The average (uncorrelated) rms spectral noise is The SNR in the average spectrum is then Equation 5 is the fundamental equation that governs the increase in SNR when averaging multiple spectra with weighting w i , assuming uncorrelated noise.We discuss three weighting schemes below.An observer's choice of weighting is dictated by their science goals and the availability of information about their data.If the noise and peak intensity are the same for all spectra such that σ i = σ 0 and T P,i = T P,0 , for the weighting schemes considered here Equation 5 reduces to Equation 6 approximates the SNR when averaging multiple spectra taken of the same source with the same integration times and observing conditions.If the noise is correlated between spectra, then the noise term in Equation 5 will include covariances between the spectra and the exponent of n in Equation 6 will be less than 0.5.One can recover the dependence on n in Equation 6by considering only the number of independent spectra.If the noise is correlated, the average noise decreases slowly when averaging and therefore the exponent in the term n 0.5 decreases and the SNR increases more slowly than in the uncorrelated case.One can recover the expected dependence on n by considering the number of independent samples.

Intensity-Noise Weighting
For "Intensity-Noise Weighting," Using this weighting will bias the average peak line intensity by the highest values of T P,i : If all spectra have the same (uncorrrelated) noise σ 0 , Equation 8 reduces to If all signal strengths T P,0 are the same but the noise is variable, as in the case of averaging data taken of the same source under different observing conditions or integration times, Equation 8becomes

Noise Weighting
For "Noise weighting," This weighting will bias the average peak line intensity toward spectra with lower noise: If all spectra have the same (uncorrelated) noise σ 0 , Equation 12 reduces to If all signal strengths T P,0 are the same, but the noise is variable, we again find Equation 10.

Uniform Weighting
For "Uniform weighting," If all spectra have the same (uncorrelated) noise σ 0 , we again find Equation 13.If all signal strengths T P,0 are the same but the noise is variable, we again find Equation 10.

Maximizing the Signal to Noise Ratio
For a given weighting scheme and averaging method, the optimal value for n is often found when the SNR reaches a maximum value SNR max , or when For intensity-noise weighting, or in the case that all spectra have the same values of T P,i and σ 0 , averaging all available spectra will result in the highest SNR (cf.Equation 6).For noise and uniform weighting, if the values of T P,i or σ i are different, the ideal number of spectra to average may be less than the total number of spectra available.
To determine the ideal number of spectra to average for noise and uniform weighting, our method requires that the spectra be ordered by decreasing SNR.For individual spectra, To determine n and SNR max , one therefore must: 1. compute or estimate the peak line intensity, T P,i , and the rms spectral noise, σ i , for all spectra; 2. order the spectra in terms of SNR i (using Equation 17); 3. determine when the average SNR is maximized (SNR max ), either theoretically using Equation 5or by fitting the average spectra with a model.
Below, we use this method to estimate SNR max and n for simulated distributions of T P .

SIMULATED SIGNAL AND NOISE DISTRIBUTIONS
We perform Monte Carlo simulations to assess the effects of different intensity and noise distributions, as well as weighting schemes.

Distributions for T P
We investigate two characteristic distributions for T P : half-normal and power law.We plot the distributions in Figure 1, for a range of half-normal standard deviations (see Section 3.1.1)and power law indices (see Section 3.1.2).The half-normal distribution is what is measured from a compact source and a Gaussian telescope response, whereas the power law distributions are meant to model diffuse (low power law indices) and compact (high power law indices) sources.For both distributions, we assume that the distribution of noise values is Gaussian, characterized by a mean value of σ 0 and a standard deviation of s σ (measured in units of the index).

Half-normal distribution for TP
We explore how a half-normal distribution of T P affects the derived values of n and SNR.A half-normal distribution can be a good approximation for data sets where the brightest spectra have a much higher signal strength than the mean.If T P follows a half-normal distribution, where T P,max is the maximum line height in the dataset and the standard deviation in the distribution of T P is s T (measured in units of the index).From Equation 5, the SNR is then To illustrate the basic functional dependencies, we can assume that the noise is uncorrelated, is the same in all spectra to be averaged, and is equal to σ 0 .In this case, for intensity-noise weighting, we have .
(20) For noise and uniform weighting, we have .
(21) As can be seen in Equations 20 and 21, in the case of constant noise the SNR depends on the ratio of the maximum line intensity divided by the noise, rather than the individual value of either quantity.We use this ratio to parameterize the simulations.
We create 100 simulated peak signal and noise distributions for values of T P,max /σ 0 of 0.01, 0.05, 0.1, 0.5, 1, 2, and 5 in two noise distributions: "constant" noise (all spectra have the same noise value) and Gaussian noise with s σ = 0.1 (each spectrum has a noise value drawn randomly from a normal distribution).All signal distributions have a standard deviation s T = 1.5.For the Gaussian noise trials, we split the analysis into two categories: 1) the estimation of the signal strength is unaffected by the noise; and 2) the estimation of the signal strength is modified by the normal distribution of standard deviation s σ .The former case represents the theoretical situation when noise does not affect the estimation of the signal; the latter case is more realistic.We analyze both cases to determine how noise affects the analysis.
For each set of distributions, we estimate SNR max using Equation 5.For noise and uniform weighting trials, we additionally compute the signal at the maximum SNR compared to the maximum signal, T P,i=n /T P,max .We give our results from all three weighting schemes in Table 1 and  Intensity-noise weighting leads to an increase in the SNR without bound and therefore sets SNR max for any averaging method.For constant-noise half-normal signal distributions and noise or uniform weighting we find: • SNR max is obtained when averaging all spectra satisfying T P ≳ 0.4T P,max ; • SNR max is ∼ 5% less than that from averaging all spectra using intensity-noise weighting, assuming T P can be reliably estimated.
• The SNR can decrease by up to 30% relative to SNR max when averaging spectra down to T P /T P,max ≃ 0.01.
• SNR max for noise and uniform weighting is ∼ 95% that found for intensity-noise weighting.
These relationships also hold for the variable noise distributions when the noise and signal strength are uncorrelated.The above are theoretical best-case scenarios.
If noise affects the estimation of the signal strength, as it does for actual data, the inability to reliably order the highest SNR spectra affects the SNR; these effects are larger if the noise is comparable to the signal.

Power Law distribution for TP
We perform a similar analysis assuming a power law distribution for T P : where the maximum value is T P,max and α is the power law index.The relevant SNR equation is then In the case of constant noise, for intensity-noise weighting, we have For noise and uniform weighting, we have Once again we see that the SNR is linearly proportional to the ratio of the maximum peak to the standard deviation.
We investigate the effect of different power law distributions for α = 0.5 to 4.5 in increments of 0.5 with T P , max/σ 0 = 1.0, and for T P , max/σ 0 = 0.1, 0.5, 1.0, 5.0, 10 and 50 with α = 2.0.We show these results in Figures 3 and in Table 3.We do not consider variable noise (and so set s σ = 0), which we assume has a minor effect, as it does for the half-normal distributions.For a constant-noise power law signal distribution with noise or uniform weighting, we find: • SNR max decreases with increasing power law index and decreasing values of T P,max /σ 0 ; • SNR max is obtained when averaging all spectra satisfying T P ≳ 0.35T P , max for the values of α and T P , max/σ 0 investigated; • As for the half-normal signal distribution, SNR max is ∼ 5% less than that from averaging all spectra using intensity-noise weighting.
• The SNR can decrease by up to 30% relative to SNR max when averaging spectra down to T P /T P,max ≃ 0.01.
• SNR max for noise and uniform weighting is ∼ 95% that found for intensity-noise weighting.size of 30 ′′ , and a spectral resolution of 0.5 km s −1 .The rms spectral noise per spaxel is ∼ 10 mK.We test the above spectral averaging methods using GDIGS data to constrain the ionic 4 He + / H + abundance ratio by number, y + .Measurements of elemental abundances provide key constraints for our understanding of Galactic chemical evolution.We define y + as where T P is the peak line intensity and ∆V is the FWHM line width.The uncertainty on y + is therefore where σ denotes parameter uncertainties.If the source is optically thin, y + measures the 4 He + / H + abundance ratio directly.
Because the mass of helium is greater than that of hydrogen, its RRL velocity is shifted by ∼ −122 km s −1 from that of hydrogen.Both lines therefore fall within the same GDIGS bandpass and are subject to the same systematic effects.
To spectrally average GDIGS RRL data using intensity-noise weighting, we: • align the spectra in velocity using the velocity centroids from the Automatic Gaussian Decomposition (AGD) results described in  • remove a fifth-order polynomial baseline and determine the SNR in the average spectrum using Gaussian fits to the hydrogen RRLs.
To spectrally average GDIGS RRL data using noise or uniform weighting, we: • determine the SNR and peak intensity for each spaxel using the results from the AGD analysis; • align the spectra in velocity using the velocity centroids from the AGD analysis; • average spectra, starting with the highest SNR spectrum; • remove a fifth-order polynomial baseline from linefree portions of the spectrum; • determine the SNR in the average spectrum using Gaussian fits to the hydrogen RRLs; • and cease averaging when the average spectrum SNR stops increasing, with a buffer of 100 spectra (once a peak in SNR is reached, continue averaging the next 100 to determine if the SNR peak is local).
For all weighting schemes, we only use spaxels fit by a single Gaussian component in the AGD.We determine y + for all average spectra by fitting the helium line using velocities from −135 km s −1 to −110 km s −1 and the hydrogen line using velocities from −20 km s −1 to +20 km s −1 .
We perform this analysis on the GDIGS H n α data in a 100 ′ × 60 ′ zone centered on the massive star forming region W43 that was first analyzed in Luisi et al. (2020).The GDIGS data of this zone has 24,000 spectra, of which 19,849 are fit in the AGD with a single hydrogen line.This zone has numerous H II regions and also diffuse ionized gas (see Luisi et al. 2020).We show the distribution of AGD-derived peak line intensities in Figure 4 for the 1000 highest-intensity values in the field.We also separate this distribution into those derived from spaxels falling within H II regions defined by the WISE Catalog of Galactic H II Regions (Anderson et al. 2014, hereafter the "WISE Catalog"), and those that  and those that are not ("DIG"; diffuse ionized gas).The "H II" spectra are more numerous in this field at all intensities studied, and all intensities TP > 0.4 K are cospatial with H II regions.A half-normal distribution fits the lower intensity values well, but drastically over-predicts the high values, indicating that the intensity distribution is more complicated than the simple models considered here.
do not fall within H II regions.The peak line intensities approximately follow a half-normal distribution of 4000 values with s T = 400, T P,max = 10 K, and that is scaled so the minimum value is 0.14 K (a power law with α = 4 also fits fairly well).The exception to this good fit is at intensities ≳ 0.6 K where the model over-predicts the data.Thus, the signal distribution is more complicated than the simulated distributions considered here.Most of the values, and a greater fraction of the high-intensity values, are associated with H II regions.
We create five different average spectra and compute y + for each: intensity-noise weighting all spectra, noise weighting with SNR maximization, uniform weighting with SNR maximization, noise weighting all spectra and uniform weighting all spectra.The noise and uniform SNR maximizations use 569 and 1102 of the ∼ 20,000 spectra, respectively, corresponding to approximate intensity values of T P > 0.15 K and T P > 0.14 K.
We show the five average spectra in Figure 5.Each spectrum is independently normalized.All five spectra have the same basic shape, although the SNR maximization spectra have the smallest deviations from a single Gaussian line.In Table 3 we summarize the H and He line height (T P ) and FWHM line width (∆V ) for the H and He RRLs, as well as their 1σ fit uncertainties, y + (Equation 26) and its 1σ uncertainty (Equation 27), and the spectral rms σ.The derived values of y + differ depending on the averaging method and the weighting scheme.Differences in y + are not accounted for by the uncertainties in σ y + .As expected, uniformly weighting all spectra results in the largest rms spectral noise; the other spectral noise values are similar.In this paper we explored methods for averaging spectra.Intensity-noise weighting leads to the highest possible SNR.For noise and uniform weighting, averaging the 35 − 45% highest intensity individual spectra (assuming similar noise characteristics for each) results in the maximum SNR average spectrum, in agreement with the results of Zhang & McElvain (1999).This average spectrum created from the 35 − 45% highest intensity individual spectra has ∼ 95% the SNR of the intensitynoise weighted average spectrum.Our results are largely independent of the intensity distribution; other peaked signal distributions should have similar results.
We apply our averaging methods to Green Bank Telescope (GBT) Diffuse Ionized Gas (GDIGS) H n α data (Anderson et al. 2021) to determine the ionic abundance ratio, y + .The different averaging methods give values of y + that differ by ∼ 25%.
Differences in the derived values of y + can be explained by which locations are weighted more heavily during averaging.Intensity-noise weighting obviously preferences the spectra with the highest peak intensity.For GDIGS, the highest intensities are found toward discrete H II regions; the highest intensity diffuse regions are found just outside of the discrete H II regions (see Luisi et al. 2020).Noise weighting preferences the spectra with the lowest noise, whereas uniform weighting weights all spectra evenly.Since H II regions have bright radio continuum emission, noise weighting can preference the diffuse regions.The SNR maximization method only averages the highest SNR spectra, which means that only the brightest regions may appear in the average, regardless of their noise levels.
That the value derived for y + depends on the weighting scheme employed indicates that there are differences in y + in the GDIGS field studied; if y + were invariant, all averaging techniques would produce the same result.This piece of evidence is not as apparent without averaging, as the He RRL signal that goes into the y + computation is weak and can only be seen in a fraction of the GDIGS spectra.We caution that studies of y + that include a range of intensity values (i.e., from both H II regions and from diffuse ionized gas, as in our example) will be biased depending on the weighting scheme.In future research with the GDIGS data, we will investigate and model y + over the survey area with these considerations in mind.
The SNR maximizing procedure allows for the creation of more sensitive spectra, and therefore a more accurate determination of y + , but the derived y + values in all average GDIGS spectra are low relative to those found previously for Galactic H II regions.For comparison, an analysis of the 80 high-quality RRL spectra toward H II regions in Quireza et al. (2006) by Wenger et al. (2013) found y + = 0.075 ± 0.024.Wenger et al. (2013) found y + = 0.068 ± 0.023 in a sample of 54 high-quality RRL spectra towards Galactic H II regions.For the H II region W43, which is in the studied field, y + = 0.068 ± 0.0052 Bania et al. (1997Bania et al. ( , 2007)).It may be that the inclusion of the diffuse ionized gas outside of H II regions has caused the discrepancy with values derived for H II regions; we will investigate the cause of the low y + values in a subsequent paper.

Figure 1 .
Figure 1.Simulated intensity distributions (see Sections 3.1.1and 3.1.2).The left panel shows the distributions themselves, while the right panel shows a histogram of the values.Half-normal distributions have a larger number of extremely low intensity values.Low power law indices have more high intensity values whereas higher power law indices have more low intensity values.
show the noise-weighting analysis in Fig-ure 2 (uniform weighting produces nearly identical results).

4.Figure 2 .
Figure 2. SNR analysis for half-normal-distributed values of TP , with noise weighting.The SNR in intensity-noise weighting (not shown) increases without bound whereas uniform weighting (also not shown) produces nearly identical results to those of noise weighting.Panels in the left column show TP,max/σ0 = 0.1 and those of the right column show TP,max/σ0 = 1.0.The top row of panels has sσ = 0 (constant noise) and the bottom row has sσ = 0.1σ0.In all panels, solid lines show individual values and dashed lines show integrated values.The blue curves show the SNR (and use the left y-axis), the red curves show TP , and the green curves show the noise (both use the right y-axis).The shaded regions in the lower panels show the standard deviations from the Monte Carlo simulations.The vertical gray lines indicate the peaks of the SNR distributions; the gray shaded areas show the range within one standard deviation of the SNR distributions.

Figure 3 .
Figure 3. SNR analysis for power law-distributed values of TP , with noise weighting.Uniform weighting (not shown) produces nearly identical results..All panels have TP , max/σ0 = 1.0, and clockwise from the top-left panel α ranges from 1 to 4 in integer steps.In all panels, solid lines show individual values and dotted lines show integrated values.The blue curves show the SNR, the red curves show TP , and the green curves show the noise.The vertical gray lines indicates the peaks of the SNR distribution.

Figure 4 .
Figure 4. GDIGS data towards W43.Shown are the 1000 highest intensity fitted line height values (TP ), separated into those that are spatially coincident with H II regions ("H II")and those that are not ("DIG"; diffuse ionized gas).The "H II" spectra are more numerous in this field at all intensities studied, and all intensities TP > 0.4 K are cospatial with H II regions.A half-normal distribution fits the lower intensity values well, but drastically over-predicts the high values, indicating that the intensity distribution is more complicated than the simple models considered here.

Figure 5 .
Figure5.Spectra created using the averaging and weighting methods discussed in the text.The data are the same in both panels, but we adjust the y-axis in the lower panel to better show the helium line centered at −122 km s −1 .Line-free regions used to fit the baseline are shaded in diagonal gray lines.Dotted black lines show the boundaries used in the fit.Light green vertical lines show the expected velocities of the hydrogen, helium, and carbon RRLs.

Table 1 .
SNR analysis for half-normal intensity distributions

Table 2 .
SNR analysis for power law intensity distributions

Table 3 .
Analysis of GDIGS data a All intensity values are normalized such that the hydrogen line intensity has a value of 1.00.