Detecting Weak Spectral Lines in Interferometric Data through Matched Filtering

The matched filter can be derived in a number of ways. Here, we introduce its derivation by maximizing the SNR of a signal, but it can equivalently be interpreted as a least squares estimator (see e.g. Schwartz & Shaw, 1975; Vio & Andreani, 2016). In general, a signal $\bm{s}$ may be corrupted by additive white noise $\bm{v}$, yielding an observation $\bm{x}=\bm{s}+\bm{v}$. To maximize the SNR of this signal by applying a linear filter $\bm{h}$, we can first write the SNR (using the definition of signal power/noise power) as

where ^∗ denotes the conjugate transpose and $\bm{R_{v}}=E[\bm{v}\bm{v}^{*}]$ is a covariance matrix of the noise $\bm{v}$, where $E[~]$ is the expectation operator. Under these conditions, the filter $\bm{h}$ which maximizes SNR is

i.e., the original signal $\bm{s}$ multiplied by the data weights $\bm{R_{v}}^{-1}$ and a normalization constant $C=1/\sqrt{\bm{s}^{*}\bm{R_{v}}^{-1}\bm{s}}$ (Woodward, 1953; North, 1963; Cumming & Wong, 2005, e.g.).

A simple application of such a filter is locating a signal within a one dimensional dataset such as an emission spectrum. In this case, a short signal $\bm{s}$ of length $n_{s}$ is embedded within a longer noisy observed spectrum $\bm{x}$ of length $n_{x}$, with the location of $\bm{s}$ within $\bm{x}$ unknown. As long as the shape of $\bm{s}$ is known, a filter kernel $\bm{h}$ can be calculated using equation 2 and cross-correlated with $\bm{x}$ to locate $\bm{s}$. This cross-correlation is often thought of as a sliding dot product, and yields a one dimensional impulse response spectrum $\bm{T}$, of length $n_{x}-n_{s}+1$. Each element of $\bm{T}$ at a position $i_{0}$ will then be:

This impulse response spectrum loses physical significance that the original observed spectrum held (it is no longer in units of power or flux). It instead encodes the degree of similarity between the observations and the filter at any given point in the observed spectrum. By projecting the noisy observations in this way, the total noise is decreased and the SNR of the signal is increased. Because the filter is linear, the Gaussian nature of the noise is preserved. The impulse response spectrum can be easily examined for evidence of $\bm{s}$, with a detection threshold set to some multiple of the standard deviation of $\bm{T}$ (e.g. 4$\sigma$), or a false positive rate scaled with the variance of $\bm{T}$.

Each interferometric visibility, $V_{i}$, is measured as the complex product of the output of the first antenna of a baseline pair in the array with the complex conjugate of the output of the second antenna (see e.g. Thompson et al., 2017). The projected baseline distance between the two antennas then defines the location of the visibility on the $(u,v)$ plane. Each visibility $V_{i}$ is associated with a unique weight $w_{i}=1/\sigma_{i}^{2}$, where $\sigma_{i}^{2}$ encodes the variance of $V_{i}$ (see Appendix C for more details on data weights in interferometric datasets). In addition to being measured at discrete spatial frequencies, the visibilities are also measured at a series of spectral frequencies (channels). Cross-correlation in this discretely sampled three dimensional $(u,v,channel)$ space is computationally awkward, but the dataset can be reshaped to a two-dimensional dataset of size $(n_{uv},n_{c})$, where each visibility row corresponds to a unique location on the $(u,v)$ plane in units of distance. Both the complex visibilities and their corresponding real weights are stored this way in the UVFITS¹¹1The UVFITS format definition can be found in AIPS Memo #114 at http://www.aips.nrao.edu/aipsmemo.html and Measurement Set (MS)²²2The MS format definition can be found at https://casa.nrao.edu/Memos/229.html formats of the Common Astronomy Software Applications package (CASA).

Transforming between visibility space and image space requires a gridding and deconvolution routine, such as CLEAN, in one direction and a visibility sampling routine in the other direction, such as uvmodel in MIRIAD, or simobserve in CASA. As using the full simobserve task is relatively slow and uvmodel is not easily interfaced with Python, we have written a Python based visibility sampling routine, vis$\textunderscore$sample, which is able to interface with CASA MS and UVFITS formats. This package builds on an implementation of the sampling algorithm in the DiskJockey package (Czekala et al., 2015; Czekala, 2016) identical to that used in uvmodel and simobserve and uses the spheroidal gridding function approximations described by Schwab (1984). As identical algorithms and gridding functions are used, output from vis$\textunderscore$sample is identical to output from uvmodel and simobserve. ³³3In addition to its utility for filter kernel generation, we note that vis$\textunderscore$sample may be useful for visibility fitting of modern interferometric datasets (e.g. MacGregor et al., 2016; Loomis et al., 2017). vis$\textunderscore$sample is publicly available under the MIT license at https://github.com/AstroChem/vis$_$sample or in the Anaconda Cloud at https://anaconda.org/rloomis/vis$_$sample

The principle assumption of a matched filter analysis is that the shape of the signal $\bm{s}$ is known, or can be reasonably approximated. In traditional applications, such as RADAR, the outbound signal is user-generated and therefore the exact form is known. In astronomical applications, however, the ideal matched filter kernel is unknown and must be approximated. As it is unknown how closely the filter approximates the true signal, any derived detection significance will be a lower limit. The method is relatively robust to choice of filter, however, as long as the filter is a reasonable approximation of the source spatio-kinematic structure.

We suggest two possible approaches: (1) calculating a kernel from a model of the source (model-driven), or (2) calculating a kernel from prior observations of strong emission lines (data-driven). In both cases the kernel is first constructed in the image plane and then Fourier transformed and visibility sampled to match the $(u,v)$ coverage of the observations. The approximated signal $\bm{f}$ and the inverted noise covariance matrix $\bm{R_{v}^{-1}}$ (calculated from the observational data weights, see Appendix C) are then used to compute the full filter kernel, including the normalization prefactor:

Fig. 2 presents three examples of the different filter kernel estimation approaches. First, for objects such as protoplanetary disks or galaxies, the source inclination and position angle are often well-known and a spatio-kinematic model of the gas can be approximated. In the top panels of Fig. 2 we have generated a Keplerian mask for molecular emission from the protoplanetary disk around TW Hya. Alternatively, a more detailed filter kernel can be generated from an astrochemical model of the source, with emission calculated using a radiative transfer code such as RADMC-3D (Dullemond, 2012) or LIME (Brinch & Hogerheijde, 2010). An example is shown in the middle panels of Fig. 2, generated from the parametric CH₃OH abundance model in Walsh et al. (2016). Details of the model are presented in Appendix D.

In the data-driven approach, an assumption is made that an observed molecular transition shares its spatio-kinematic pattern with the desired weak line. The filter will be most effective when the template and target lines have well-matched spatial distributions, e.g. if the two lines are a strong and a weak line, respectively, of the same molecule. In Carney et al. (2017), we used this approach to detect weak H₂CO lines in HD 163296, using a stronger H₂CO line as a data-driven filter (see §4.1 and §4.4). Similarly, lines of a known species can be used as a filter for an undetected but chemically related molecule that is presumed to be co-spatial. The bottom row panels of Fig. 2 present observations of H₂CO around TW Hya (Öberg et al., 2017) which could be used as a filter for CH₃OH emission, due to their linked formation pathways (e.g., Cuppen et al., 2009; Qi et al., 2013; Walsh et al., 2014; Loomis et al., 2015). Details of the observations and kernel generation are presented in Appendix D.

Fig. 3 schematically diagrams how these filter kernels would be applied to the data to produce an impulse response spectrum. First, the image plane kernel is Fourier transformed and visibility sampled to produce a complex $(u,v)$ plane kernel $\bm{f}$ of size $(n_{uv},n_{k})$. The inverted noise covariance matrix $\bm{R_{v}^{-1}}$ is calculated from the data weights $\bm{w}$ and combined with $\bm{f}$ using Eq. 4 to produce the full filter kernel $\bm{h}$. This kernel is then cross-correlated with the data $\bm{V}$, of size $(n_{uv},n_{c})$. The kernel and the data both have the same number of visibilities, $n_{uv}$, but different numbers of channels, and the kernel slides through the data along the spectral axis. At each channel, the filter impulse response spectrum $\bm{T}$ is calculated by taking the complex inner product of the windowed data with the kernel:

As the data is complex, the filter response spectrum is also complex with a normalized total noise power. Signal power will leak from the real to the imaginary portion of the response, however, if there is a phase misalignment between the sky locations of the filter and the source. Thus if the filter has been properly phase shifted to be aligned with the source, the resultant impulse response spectrum $\bm{T}$ can be written as:

with the factor of $\sqrt{2}$ introduced to normalize the noise power in the real portion of the spectrum.

This method of calculating the cross-correlation is conceptually simple, but computationally inefficient. Computing inner products of the windowed data requires either manipulation of the (very large) dataset in memory or non-sequential memory access, preventing speed increases through vectorization.⁴⁴4Using FFT cross-correlation is even slower for typical interferometric dataset sizes. There is no restriction, however, on the order of operations in which the inner products are internally calculated. We use this to our advantage and treat the partial two-dimensional cross-correlation as a series of $n_{uv}$ one dimensional cross-correlations along the spectral axis, yielding $n_{uv}$ individual impulse response curves. The UVFITS and MS data formats store visibilities in a row-major order such that these one-dimensional cross-correlations quickly access data sequentially in memory. The resulting impulse response curves are then summed along the spatial frequency dimension, identical to Eq. 6, but with the order of the summations switched,

yielding a final impulse response spectrum identical to that from the sliding window method shown in Fig. 3. The $n_{uv}$ one-dimensional cross-correlations are independent, making parallelization trivial. Using this approach, the full bandwidth of a typical ALMA spectral window (e.g. 3840 channels over a 1 hr integration with 43 antennas) can be filtered very quickly on a desktop (e.g. a few seconds on a quad-core 3.3GHz processor).

Assessing the probability of a line detection from the filter impulse response requires understanding the noise properties of the response spectrum, which no longer holds the same physical significance as an emission spectrum. The response spectrum at a given frequency now represents how closely the data correspond to the filter, rather than the flux at that frequency. As the filter is linear, uncorrelated thermal noise in the visibilities manifests as Gaussian noise in the filter response. If the data weights are properly calibrated (see Appendix C), the prefactor in Eq. 4 will normalize the filter response such that the spectrum has units of standard deviations ($\sigma$) with a RMS noise level of unity.

In practice, however, we have found that the calculated data weights are often only accurate to $\sim$20% compared to the actual variance of the visibilities. Thus we strongly recommend comparing the data weights and visibility scatter and recalculating the weights using a task such as statwt in CASA if there is a discrepancy. Alternatively, the filter response itself can be manually normalized by dividing by the standard deviation of the spectrum (excluding any channels with obvious signal).

Additionally, the linear nature of the filter means that any unsubtracted continuum emission will result in a constant offset in the response spectrum. A model of the continuum can either be subtracted in the $(u,v)$ plane as a pre-processing step prior to filtering using a task such as uvcontsub in CASA, or be subtracted after filtering by subtracting the mean of signal-free channels in the impulse response spectrum. Subtraction after filtering removes the possibility of small inaccuracies in the $(u,v)$ subtracted model, to which the filter will be sensitive. Conversely, subtraction prior to filtering may be more convenient when a large bandwidth is covered for a source with a non-zero spectral index.

Once the response spectrum is normalized and any offsets are removed, peaks can then be evaluated against a detection threshold, set at some number of standard deviations corresponding to an acceptable false alarm rate. We stress, however, that the detection significance is a lower limit, as it is unknown how closely the filter approximates the ideal matched filter.

As a proof of concept, we apply the method to synthetic observations of CH₃OH emission in TW Hya and compare to an aperture-based spectral extraction in the image plane. The modeled emission from the middle panels of Fig. 2 is used to generate both the observations (with noise added) and the filter kernel. As this is a true matched filter, it provides a useful benchmark for comparison with the approximate filter results as presented in §3.

A synthetic measurement set of observations was created from the CH₃OH emission cube described in §2.3 by visibility sampling at baselines corresponding to the observations from Walsh et al. (2016) using vis$\textunderscore$sample. The complex visibilities were then noise corrupted such that the rms noise was 5 mJy bm^-1 across each 0.15 km/s channel, equivalent to the Walsh et al. (2016) observations. The noiseless and noisy measurement sets were imaged in CASA using the CLEAN task, with a CLEAN mask generated from the LIME output emission profile and a circular 1$\arcsec$ FWHM Gaussian taper applied in the Fourier plane to increase the SNR of the images. Only the noiseless measurement set was CLEANed; the noisy measurement set was dirty imaged to prevent bias from over-CLEANing, as the emission is practically at the noise limit in any given channel. Integrated intensity (moment-0) maps of the noiseless and noisy 3₁₂-3₀₃ transitions are shown in Fig. 4, panels a and b, respectively, and were generated by integrating all channels with emission. No clipping threshold was used. In the noisy case, the moment-0 map rms is $\sim$3.3 mJy bm^-1 km s^-1 and the peak integrated flux is $\sim$13.2 mJy bm^-1 km s^-1, yielding a SNR of $\sim$4. A spectrum was extracted from the noisy image cube using an aperture 3$\arcsec$ in diameter, equivalent to the extent of the CH₃OH emission (Fig. 4, panel c). The spectrum has a peak flux of $\sim$11.4 mJy and a rms noise of $\sim$3.2 mJy, yielding a SNR of $\sim$3.5$\sigma$. The rms noise level of the noisy spectrum was estimated from all channels without significant emission (i.e., excluding a velocity range of $\pm$ 1.5 km s^-1 around the systemic velocity of 2.8 km s^-1).

We cross-correlate the CH₃OH filter kernel with the synthetic observations, as described in §2.4, generating the filter response shown in Fig. 4 panel d. The peak value of the filter response, 5.7$\sigma$, is the maximum SNR extractable from the data and represents a $\sim$40$\%$ and $\sim$60$\%$ improvement over the moment-0 and spectral detections, respectively. This already corresponds to a factor of 2–3 increase in effective observing time, but as discussed in both §4.1 and Appendices A & B, the level of possible SNR improvements will be higher for data sets that are better-resolved.

Spectral stacking is a common method of SNR improvement for observations of multiple transitions of the same molecule (e.g. Langston & Turner, 2007; Kalenskii & Johansson, 2010; Loomis et al., 2016; Walsh et al., 2016). If the excitation conditions of two or more transitions are similar and their rest frequencies are well known, then the signals can be combined through weighted averaging,

where the stacked spectrum $\bm{T_{s}}$ is generated by summing $n_{s}$ individual spectra $\bm{T_{i}}$ multiplied by weights $w_{i}$, proportional to the SNR of each $\bm{T_{i}}$. Knowledge of the relative strengths of each transition is therefore important to gain the most signal improvement. Application of a matched filter results in an estimated SNR for each transition, which can be used as a proxy for their relative strengths. The resultant impulse response spectra are then easily stacked to generate an appropriately weighted stacked spectrum.

To illustrate this process, we have repeated the simulations and filtering described in §2.6 for three CH₃OH transitions: 2₁₁-2₀₂, 3₁₂-3₀₃, and 4₁₃-4₀₄, with relative strengths of 1.8:1.3:1.0. Moment-0 maps of the emission from each of these transitions are shown in Fig. 5, panels a, b, and c, with peak integrated fluxes of 11.7, 13.2, and 8.5 mJy km s^-1 and corresponding SNRs of 3.5, 4, and 2.6$\sigma$, respectively. The individual filter responses are shown in Fig. 5, panels d, e, and f, with peak SNRs of 8.4, 5.7, and 4.9$\sigma$, respectively. The filter responses were stacked using a weighted average, yielding the spectrum shown in Fig. 5, panel g, with a peak SNR of 11.6$\sigma$. The ratio of the filter responses $\sim$(1.7:1.2:1) recovers the flux ratio of the input models (1.8:1.3:1.0) fairly well, even though the 2₁₁-2₀₂ transition appears weaker than would be expected in the imaged data (likely due to random noise fluctuations in the inherently more noisy moment-0 maps). This highlights one of the advantages of applying the matched filter in the Fourier plane.

Compared with traditional line detection methods, the matched filter approach offers an improved SNR. The degree of SNR boost depends on both the accuracy of the approximated kernel as well as the specific properties of the data (particularly the spatial resolution). In the synthetic and real data test cases presented in §3, we found that application of a matched filter could increase SNR by $\sim$60%. The method was found to be similarly effective when applied to real data ($\sim$53% vs $\sim$60% boost), demonstrating that it is robust to the choice of approximated filter.

Intuitively, the SNR boost and spatial resolution of the emission should be coupled. By definition, a spatially unresolved signal encodes no spatio-kinematic information, and in this limit the matched filter technique will provide no increase in SNR other than the boost from spectral averaging. As emission is spatially resolved, SNR will decrease roughly with the square of the degree of spatial resolution (source width / beam size), with additional losses due to spatial filtering (see e.g. Crane & Napier, 1986). With appropriate knowledge of the velocity structure, a matched filter essentially negates this effect, and thus the SNR boost scales directly with the spatial resolution of the signal (see Yen et al., 2016, for a detailed image-plane derivation of this SNR boost). Fig. 8 illustrates this effect, with an identical simulation to that shown in Fig. 4, but with a higher spatial resolution of $\sim$0$\farcs$3. The data were noise corrupted to reach a similar $\sim$4$\sigma$ detection in the moment map, although the SNR in the extracted spectrum is now $\sim$6$\sigma$, highlighting how ineffective moment maps are at high spatial resolutions. The filter response is also now much larger (SNR=23.6$\sigma$), yielding a SNR boost over the aperture extracted spectra of $\sim$400%, compared to $\sim$60% in Fig. 4. Similarly, application of matched filtering to higher resolution observations of H₂CO (Carney et al., 2017) produced a SNR gain of over 500%, confirming in practice the relationship between SNR boost and spatial resolution.

Recently Matrà et al. (2015) and Marino et al. (2016) introduced an image-plane line detection technique (also independently introduced and formalized by Yen et al., 2016) that provides some similar benefits to the matched filter technique. In their approaches, pixels from a dirty image are adjusted for an assumed velocity offset (from a source model), and the velocity corrected spectra are then stacked. In many ways, this can be seen as an image-plane analog to a Fourier plane matched filter, and it should yield comparable increases in SNR (see Appendix B for more details). This is confirmed by comparing the results of the matched filtering technique to detect H₂CO in HD 163296 (Carney et al., 2017) with those obtained on the same dataset by Yen et al. (2016). In both cases, a SNR boost of $\sim$500% is achieved.

Several subtle differences between matched filtering and pixel stacking, however, may motivate their use in a synergistic fashion. First, application of a matched filter in the uv-plane requires no imaging of the data, and is therefore much faster and more robust than image-plane spectral stacking. Second, the matched filter technique allows for a more accurate emission model than simple Keplerian rotation to be applied to the data (i.e., the spatial distribution of molecular signal can be properly used for weighting). On the other hand, stacking pixels in the image plane makes extracting a flux measurement for the line much simpler and allows for a radial profile to be estimated (Yen et al., 2016). The two techniques could therefore be used sequentially to exploit the benefits of each, with a matched filter first used to quickly identify and confirm line detections in a dataset and pixel stacking then used to better characterize the lines.

The main utility of matched filtering when applied to interferometric spectral line data is in the rapid detection of weak lines, rather than their detailed characterization. Once a line is identified, it might be further characterized through careful imaging, spectral stacking, or model-fitting to the visibilities. For the weakest lines, however, detailed characterization will likely require additional observations. Matched filtering provides useful predictive utility when planning these observations, robustly confirming weak lines which might be desirable targets.

In particular, after a weak line is identified, the matched filter method can be used to estimate a line flux if the emission is too weak to be seen directly in the image cube. When a data-driven approach is taken, the responses of the target and template lines to the filter can be compared. If the two lines have a similar emission morphology, the ratio of their responses will be similar to the ratio of their fluxes, with the flux of the strong line being easy to measure.

This can be proven by considering a modified version of equation 1, writing down the SNR using the signal/rms definition:

We can treat the two lines as signals $\bm{y}$ and $\bm{s}$, where $\bm{y}$ differs only from $\bm{s}$ by an arbitrary constant $\alpha$, i.e. $\bm{y}=\alpha\bm{s}$. The SNR of $\bm{y}$ after filter application is:

which reduces to:

Therefore we can use the filter impulse response ratio to roughly estimate the flux, with the accuracy dependent on the closeness of the filter kernel to the true emission distribution. It is important to note that this estimate will always be a lower limit. An upper limit can additionally be derived from the imaged weak line, bounding the flux measurement. This approach was used by Carney et al. (2017) to determine the flux ratio of multiple detected H₂CO lines, enabling them to constrain the H₂CO excitation temperature.

In addition to aiding the detection of specific known weak lines, interferometric matched filtering provides substantial benefits for the processing of spectral line surveys where source locations and approximate spatio-kinematic structures are known. Imaging the full bandwidth of these large datasets at their native spectral resolution is a time consuming process, often taking many hours or even days. Because much of the information in these datasets is contained in spatio-kinematic patterns of the spectral lines, decreasing spectral resolution through channel averaging is typically not a viable option and can result in signal loss. A choice must therefore be made between imaging only small targeted windows of the broadband dataset, or spending time and computing resources on imaging the full bandwidth. For sparsely populated line surveys (e.g., of protoplanetary disks), imaging the entire data set is inherently inefficient, since most of the channels do not contain signal. Conversely, selective imaging reduces the likelihood of serendipitously detecting weak species, and conflicts with the motivations of an unbiased survey.

Numerous tools have been developed to aid in identifying spectral line emission in broadband datasets from current and future instruments such as ALMA, ASKAP, VLA, SKA, and the ngVLA (e.g., Koribalski, 2012; Whiting, 2012; Whiting & Humphreys, 2012; Friedel et al., 2015; Serra et al., 2015), but these methods often rely on a fully imaged datacube as input. In instances where the locations of the sources being targeted are known a-priori, our presently described method of matched filtering can help streamline this process by quickly and robustly identifying lines in the native visibilities. Then only these lines need be imaged and analyzed. In sources with a single dominant velocity pattern, a strong line could be imaged, converted to a filter kernel, and cross-correlated through the entire dataset in a small fraction of the time it would take to image that same dataset. The resulting full-band impulse response spectrum then provides a convenient first look at the dataset, guiding the observer as to which sections of the data are worth windowing out for further imaging and analysis. In particular, matched filtering will highlight weak lines that the observer would have missed even in a careful CLEAN of the data.

We note that our $(u,v)$-plane method could likely be extended to full 3D searches in blind surveys, where source locations are not known a-priori, but such an implementation is beyond the scope of this paper. It is not immediately clear whether the large speed benefits of $(u,v)$-plane analysis over full survey imaging would be maintained in such a 3D search space, and we encourage further research in this area.