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 S/N 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 ${\boldsymbol{s}}$ may be corrupted by additive white noise ${\boldsymbol{v}}$, yielding an observation ${\boldsymbol{x}}={\boldsymbol{s}}+{\boldsymbol{v}}$. To maximize the S/N of this signal by applying a linear filter ${\boldsymbol{h}}$, we can first write the S/N (using the definition of signal power/noise power) as

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\displaystyle \frac{{{\boldsymbol{h}}}^{* }{{\boldsymbol{ss}}}^{* }{\boldsymbol{h}}}{{{\boldsymbol{h}}}^{* }{{\boldsymbol{R}}}_{v}{\boldsymbol{h}}},\end{eqnarray} \tag{ 1 }$$

where * denotes the conjugate transpose and ${{\boldsymbol{R}}}_{v}=E[{{\boldsymbol{vv}}}^{* }]$ is a covariance matrix of the noise ${\boldsymbol{v}}$, where E[ ] is the expectation operator. Under these conditions, the filter ${\boldsymbol{h}}$ which maximizes S/N is

$$\begin{eqnarray}&&{\boldsymbol{h}}=\left[\displaystyle \frac{1}{\sqrt{{{\boldsymbol{s}}}^{* }{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}}}}\right]{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}}=C{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}},\end{eqnarray} \tag{ 2 }$$

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

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

$$\begin{eqnarray}&&{T}_{{i}_{0}}=\displaystyle \sum _{i={i}_{0}}^{{i}_{0}+{N}_{s}-1}{x}_{i}{h}_{i-{i}_{0}};\,\,\,\,{i}_{0}\in [0,{n}_{x}-{n}_{s}].\end{eqnarray} \tag{ 3 }$$

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 S/N 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 ${\boldsymbol{s}}$, with a detection threshold set to some multiple of the standard deviation of ${\boldsymbol{T}}$ (e.g., 4σ), or a false positive rate scaled with the variance of ${\boldsymbol{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 data sets). 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 data set can be reshaped to a two-dimensional data set 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⁶ and Measurement Set (MS)⁷ 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_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_sample is identical to output from uvmodel and simobserve.⁸

The principal assumption of a matched filter analysis is that the shape of the signal ${\boldsymbol{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 ${\boldsymbol{f}}$ and the inverted noise covariance matrix ${{\boldsymbol{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:

$$\begin{eqnarray}&&{\boldsymbol{h}}=\left[\displaystyle \frac{1}{\sqrt{{{\boldsymbol{f}}}^{* }{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{f}}}}\right]{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{f}}.\end{eqnarray} \tag{ 4 }$$

Figure 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 Figure 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 Figure 2, generated from the parametric CH₃OH abundance model in Walsh et al. (2016). Details of the model are presented in Appendix D.

Zoom In Zoom Out Reset image size

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 Sections 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 Figure 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.

Figure 3 schematically describes 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 ${\boldsymbol{f}}$ of size (n_uv, n_k). The inverted noise covariance matrix ${{\boldsymbol{R}}}_{v}^{-1}$ is calculated from the data weights ${\boldsymbol{w}}$ and combined with ${\boldsymbol{f}}$ using Equation (4) to produce the full filter kernel ${\boldsymbol{h}}$. This kernel is then cross-correlated with the data ${\boldsymbol{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 ${\boldsymbol{T}}$ is calculated by taking the complex inner product of the windowed data with the kernel:

$$\begin{eqnarray}&&{T}_{{i}_{0}}=\displaystyle \sum _{i={i}_{0}}^{{i}_{0}+{n}_{k}-1}\displaystyle \sum _{j=0}^{{n}_{{uv}}}{V}_{i,j}{h}_{i-{i}_{0},j};\,\,\,{i}_{0}\in [0,{n}_{c}-{n}_{k}].\end{eqnarray} \tag{ 5 }$$

As the data are 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 ${\boldsymbol{T}}$ can be written as:

$$\begin{eqnarray}&&{T}_{{i}_{0}}=\mathrm{Re}\left[\sqrt{2}\displaystyle \sum _{i={i}_{0}}^{{i}_{0}+{n}_{k}-1}\displaystyle \sum _{j=0}^{{n}_{{uv}}}{V}_{i,j}{h}_{i-{i}_{0},j}\right];\ {i}_{0}\in [0,{n}_{c}-{n}_{k}],\end{eqnarray} \tag{ 6 }$$

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

Zoom In Zoom Out Reset image size

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) data set in memory or non-sequential memory access, preventing speed increases through vectorization.⁹ 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 Equation (6), but with the order of the summations switched,

$$\begin{eqnarray}&&{T}_{{i}_{0}}=\mathrm{Re}\left[\sqrt{2}\displaystyle \sum _{j=0}^{{n}_{{uv}}}\displaystyle \sum _{i={i}_{0}}^{{i}_{0}+{n}_{k}-1}{V}_{i,j}{h}_{i-{i}_{0},j}\right];\,\,\,{i}_{0}\in [0,{n}_{c}-{n}_{k}],\end{eqnarray} \tag{ 7 }$$

yielding a final impulse response spectrum identical to that from the sliding window method shown in Figure 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.3 GHz 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 Equation (4) will normalize the filter response such that the spectrum has units of standard deviations (σ) with a rms noise level of unity.

In practice, however, we have found that the calculated data weights are often only accurate to ∼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 Figure 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 Section 3.

A synthetic MS of observations was created from the CH₃OH emission cube described in Section 2.3 by visibility sampling at baselines corresponding to the observations from Walsh et al. (2016) using vis_sample. The complex visibilities were then noise corrupted such that the rms noise was 5 mJy bm⁻¹ across each 0.15 km s⁻¹ channel, equivalent to the Walsh et al. (2016) observations. The noiseless and noisy MSs were imaged in CASA using the CLEAN task, with a CLEAN mask generated from the LIME output emission profile and a circular 1'' FWHM Gaussian taper applied in the Fourier plane to increase the S/N of the images. Only the noiseless MS was CLEANed; the noisy MS 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 Figure 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 ∼3.3 mJy bm⁻¹ km s⁻¹ and the peak integrated flux is ∼13.2 mJy bm⁻¹ km s⁻¹, yielding a S/N of ∼4. A spectrum was extracted from the noisy image cube using an aperture 3'' in diameter, equivalent to the extent of the CH₃OH emission (Figure 4, panel (c)). The spectrum has a peak flux of ∼11.4 mJy and a rms noise of ∼3.2 mJy, yielding an S/N of ∼3.5σ. The rms noise level of the noisy spectrum was estimated from all channels without significant emission (i.e., excluding a velocity range of ±1.5 km s⁻¹ around the systemic velocity of 2.8 km s⁻¹).

Zoom In Zoom Out Reset image size

We cross-correlate the CH₃OH filter kernel with the synthetic observations, as described in Section 2.4, generating the filter response shown in Figure 4, panel (d). The peak value of the filter response, 5.7σ, is the maximum S/N extractable from the data and represents a ∼40% and ∼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 Section 4.1 and Appendices A and B, the level of possible S/N improvements will be higher for data sets that are better-resolved.

Spectral stacking is a common method of S/N 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,

$$\begin{eqnarray}&&{{\boldsymbol{T}}}_{s}=\displaystyle \sum _{i=0}^{{n}_{s}}{{\boldsymbol{T}}}_{i}{w}_{i},\end{eqnarray} \tag{ 8 }$$

where the stacked spectrum ${{\boldsymbol{T}}}_{s}$ is generated by summing n_s individual spectra ${{\boldsymbol{T}}}_{i}$ multiplied by weights w_i, proportional to the S/N of each ${{\boldsymbol{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 S/N 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 Section 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 Figure 5, panels (a)–(c), with peak integrated fluxes of 11.7, 13.2, and 8.5 mJy km s⁻¹ and corresponding S/Ns of 3.5, 4, and 2.6σ, respectively. The individual filter responses are shown in Figure 5, panels (d)–(f), with peak S/Ns of 8.4, 5.7, and 4.9σ, respectively. The filter responses were stacked using a weighted average, yielding the spectrum shown in Figure 5, panel (g), with a peak S/N of 11.6σ. The ratio of the filter responses ∼(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.

Zoom In Zoom Out Reset image size

Compared with traditional line detection methods, the matched filter approach offers an improved S/N. The degree of S/N 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 Section 3, we found that application of a matched filter could increase S/N by ∼60%. The method was found to be similarly effective when applied to real data (∼53% versus ∼60% boost), demonstrating that it is robust to the choice of approximated filter.

Intuitively, the S/N 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 S/N other than the boost from spectral averaging. As emission is spatially resolved, S/N 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 S/N boost scales directly with the spatial resolution of the signal (see Yen et al. 2016 for a detailed image-plane derivation of this S/N boost). Figure 8 illustrates this effect, with an identical simulation to that shown in Figure 4, but with a higher spatial resolution of ∼03. The data were noise corrupted to reach a similar ∼4σ detection in the moment map, although the S/N in the extracted spectrum is now ∼6σ, highlighting how ineffective moment maps are at high spatial resolutions. The filter response is also now much larger (S/N = 23.6σ), yielding an S/N boost over the aperture extracted spectra of ∼400%, compared to ∼60% in Figure 4. Similarly, application of matched filtering to higher-resolution observations of H₂CO (Carney et al. 2017) produced an S/N gain of over 500%, confirming in practice the relationship between S/N boost and spatial resolution.

Zoom In Zoom Out Reset image size

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 S/N (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 data set by Yen et al. (2016). In both cases, an S/N boost of ∼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 (u, v)-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 data set 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 S/N using the signal/rms definition:

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{{\rm{s}}}=\sqrt{\displaystyle \frac{{{\boldsymbol{h}}}^{* }{{\boldsymbol{ss}}}^{* }{\boldsymbol{h}}}{{{\boldsymbol{h}}}^{* }{{\boldsymbol{R}}}_{v}{\boldsymbol{h}}}}.\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{{\rm{y}}}=\sqrt{\displaystyle \frac{{{\boldsymbol{h}}}^{* }\alpha {\boldsymbol{s}}\alpha {{\boldsymbol{s}}}^{* }{\boldsymbol{h}}}{{{\boldsymbol{h}}}^{* }{{\boldsymbol{R}}}_{v}{\boldsymbol{h}}}}\end{eqnarray} \tag{ 10 }$$

which reduces to:

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{{\rm{y}}}=\alpha \sqrt{\displaystyle \frac{{{\boldsymbol{h}}}^{* }{{\boldsymbol{ss}}}^{* }{\boldsymbol{h}}}{{{\boldsymbol{h}}}^{* }{{\boldsymbol{R}}}_{v}{\boldsymbol{h}}}}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\rm{S}}/{{\rm{N}}}_{{\rm{y}}}=\alpha \,{\rm{S}}/{{\rm{N}}}_{{\rm{s}}}.\end{eqnarray} \tag{ 12 }$$

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 data sets at their native spectral resolution is a time-consuming process, often taking many hours or even days. Because much of the information in these data sets 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 data set, 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 data sets 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 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 data set in a small fraction of the time it would take to image that same data set. The resulting full-band impulse response spectrum then provides a convenient first look at the data set, 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.

In general, each visibility V_i in an interferometric data set corresponding to an antenna pair (m, n) will have a characteristic variance ${\sigma }_{{mn}}^{2}$, which is often assumed to be dominated by the system noise (see e.g., Chapter 6 of Thompson et al. 2017). When the system noise dominates, σ_mn can be written in units of Jy as:

$$\begin{eqnarray}&&{\sigma }_{{mn}}(\mathrm{Jy})=\displaystyle \frac{2k}{{n}_{q}{n}_{c}{A}_{\mathrm{eff}}}\sqrt{\displaystyle \frac{{T}_{\mathrm{sys},m}{T}_{\mathrm{sys},n}}{2{\rm{\Delta }}\nu {\rm{\Delta }}t}}\times {10}^{26},\end{eqnarray} \tag{ 28 }$$

where k is Boltzmann's constant, n_q and n_c are the quantization and correlator efficiencies, A_eff is the effective antenna area, T_sys,m and T_sys,n are system temperatures for antennas m and n, respectively, Δν is the effective channel bandwidth, and Δt is the integration time. For identical antennas and system temperatures, this can be simplified as:

$$\begin{eqnarray}&&{\sigma }_{{mn}}(\mathrm{Jy})=\displaystyle \frac{2{{kT}}_{\mathrm{sys}}}{{n}_{q}{n}_{c}{A}_{\mathrm{eff}}}\sqrt{\displaystyle \frac{1}{2{\rm{\Delta }}\nu {\rm{\Delta }}t}}\times {10}^{26}.\end{eqnarray} \tag{ 29 }$$

The data weight for each visibility is then calculated as ${w}_{i}=1/{\sigma }_{i}^{2}$. For instruments which record channelized system temperatures (e.g., ALMA), the weights will also be channelized, and are recorded in CASA as a "weight spectrum" for each visibility.¹⁰

As discussed in Section 2.1, a matched filter can be written as:

$$\begin{eqnarray}&&{\boldsymbol{h}}=\left[\displaystyle \frac{1}{\sqrt{{{\boldsymbol{s}}}^{* }{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}}}}\right]{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}}=C{{\boldsymbol{R}}}_{v}^{-1}{\boldsymbol{s}},\end{eqnarray} \tag{ 30 }$$

where ${{\boldsymbol{R}}}_{v}$ is the noise covariance matrix. When the n_c channels in an interferometric data set are fully independent, ${{\boldsymbol{R}}}_{v}$ can be written for an individual visibility V_i as a ${n}_{c}\times {n}_{c}$ diagonal matrix:

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v}=\left[\begin{array}{cccc}{\sigma }_{1}^{2} & & & \\ & {\sigma }_{2}^{2} & & \\ & & \ddots & \\ & & & {\sigma }_{{n}_{c}}^{2}\end{array}\right]\end{eqnarray} \tag{ 31 }$$

and ${{\boldsymbol{R}}}_{v}^{-1}$ is then:

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v}^{-1}=\left[\begin{array}{cccc}\tfrac{1}{{\sigma }_{1}^{2}} & & & \\ & \tfrac{1}{{\sigma }_{2}^{2}} & & \\ & & \ddots & \\ & & & \tfrac{1}{{\sigma }_{{n}_{c}}^{2}}\end{array}\right]\end{eqnarray} \tag{ 32 }$$

i.e., a diagonal matrix initialized to the channelized data weights. When the weights are not channelized, or can be approximated as equal across channels (w_j ≈ w_i ∀ j), ${R}_{v}^{-1}$ can be written as:

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{v}^{-1}={w}_{i}{{\boldsymbol{I}}}_{{n}_{c}},\end{eqnarray} \tag{ 33 }$$

where ${{\boldsymbol{I}}}_{{n}_{c}}$ is an ${n}_{c}\times {n}_{c}$ identity matrix.

C.2.1. Correlated Channels

In the case of fully independent channels, the filter and normalization prefactor are simple to calculate due to the lack of off-diagonal elements in the noise covariance matrix. In practice, however, channels are often correlated. In particular, astronomical interferometric data sets are generally correlated due to the use of a window function (e.g., Hann, Hamming, etc.) to reduce the ringing effect introduced by the finite maximum lag time of the correlator hardware.¹¹ As the Hann function is a popular choice of window function (and applied directly in the time domain for the ALMA data described in this paper), we derive here the appropriate noise covariances matrices for Hanning-smoothed data. Similar results could be calculated for other choices of window function.

In the frequency/channel domain, the Hann window applied to an observation ${\boldsymbol{x}}$ can be written as:

$$\begin{eqnarray}&&{x}_{i}^{{\prime} }=\displaystyle \frac{1}{4}{x}_{i-1}+\displaystyle \frac{1}{2}{x}_{i}+\displaystyle \frac{1}{4}{x}_{i+1}.\end{eqnarray} \tag{ 34 }$$

For an observation which can be linearly decomposed into signal and additive white Gaussian noise (${\boldsymbol{x}}={\boldsymbol{s}}+{\boldsymbol{v}}$), we can calculate the noise covariance matrix of the smoothed data, ${{\boldsymbol{R}}}_{v}^{{\prime} }$, which now contains off-diagonal elements:

$$\begin{eqnarray}&&{R}_{v}^{{\prime} }[j,k]=E[{v}_{j}^{{\prime} }v{{\prime} }_{k}^{* }],\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{v}_{j}^{{\prime} }=\displaystyle \frac{1}{4}{v}_{j-1}+\displaystyle \frac{1}{2}{v}_{j}+\displaystyle \frac{1}{4}{v}_{j+1}.\end{eqnarray} \tag{ 36 }$$

Noting that in the uncorrelated case

$$\begin{eqnarray}E[{v}_{j}{v}_{k}^{* }]=\left\{\begin{array}{ll}{\sigma }_{j}^{2}, & \mathrm{if}\ j=k\\ 0, & \mathrm{otherwise}\end{array}\right.\end{eqnarray} \tag{ 37 }$$

we find that for diagonal elements of ${R}_{v}^{{\prime} }$, Equation (35) reduces to

$$\begin{eqnarray}&&{R}_{v}^{{\prime} }[j,j]=\displaystyle \frac{{\sigma }_{j-1}^{2}}{16}+\displaystyle \frac{{\sigma }_{j}^{2}}{4}+\displaystyle \frac{{\sigma }_{j+1}^{2}}{16},\end{eqnarray} \tag{ 38 }$$

and for the populated off-diagonal elements it reduces to

$$\begin{eqnarray}&&{R}_{v}^{{\prime} }[j,j\pm 1]=\displaystyle \frac{{\sigma }_{j}^{2}}{8}+\displaystyle \frac{{\sigma }_{j\pm 1}^{2}}{8},\end{eqnarray} \tag{ 39 }$$

and

$$\begin{eqnarray}&&{R}_{v}^{{\prime} }[j,j\pm 2]=\displaystyle \frac{{\sigma }_{j\pm 1}^{2}}{16}.\end{eqnarray} \tag{ 40 }$$

Under the assumption that the channelized weights are all approximately equal for a given visibility (w_j ≈ w_i ∀ j), ${{\boldsymbol{R}}}_{v}^{{\prime} }$ can be written as

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v}^{{\prime} }\approx {\sigma }_{i}^{2}\left[\begin{array}{ccccc}\tfrac{3}{8} & \tfrac{1}{4} & \tfrac{1}{16} & & \\ \tfrac{1}{4} & \tfrac{3}{8} & \ddots & \ddots & \\ \tfrac{1}{16} & \ddots & \ddots & \ddots & \tfrac{1}{16}\\ & \ddots & \ddots & \tfrac{3}{8} & \tfrac{1}{4}\\ & & \tfrac{1}{16} & \tfrac{1}{4} & \tfrac{3}{8}\end{array}\right]\\ \quad \,=\,{\sigma }_{i}^{{\prime} 2}\left[\begin{array}{ccccc}1 & \tfrac{2}{3} & \tfrac{1}{6} & & \\ \tfrac{2}{3} & 1 & \ddots & \ddots & \\ \tfrac{1}{6} & \ddots & \ddots & \ddots & \tfrac{1}{6}\\ & \ddots & \ddots & 1 & \tfrac{2}{3}\\ & & \tfrac{1}{6} & \tfrac{2}{3} & 1\end{array}\right]\,=\,{\sigma }_{i}^{{\prime} 2}\,[{\boldsymbol{M}}]\end{eqnarray} \tag{ 41 }$$

where ${\sigma }_{i}^{{\prime} }=\sqrt{3/8}{\sigma }_{i}$. As tasks such as statwt in CASA do not consider effective channel bandwidth, ${w}_{i}^{{\prime} }=1/{\sigma }_{i}^{{\prime} 2}$ is likely what will be reported as the data weights of the observations.¹² The assumption of equal weights across channels is reasonable when T_sys is stable across channels, which will likely be true unless there were issues with data calibration or there are strong water lines in the data. Under this assumption, the matrix inversion only needs to be computed once:

$$\begin{eqnarray}&&{{\boldsymbol{R}}}_{v}^{{\prime} -1}\approx \displaystyle \frac{1}{{\sigma }_{i}^{{\prime} 2}}\,{[{\boldsymbol{M}}]}^{-1}={w}_{i}^{{\prime} }\,{[{\boldsymbol{M}}]}^{-1}.\end{eqnarray} \tag{ 42 }$$

When the data are initially Hanning smoothed, the edge channels are clipped. Thus for the purposes of preventing boundary effects during inversion, the covariance matrix can be treated as describing a infinitely wide data set. In practice, we extend the matrix to be several times larger than n_c, and then window a n_k sized portion from the center of ${{\boldsymbol{R}}}_{v}^{{\prime} -1}$. In the unbinned case, however, ${{\boldsymbol{R}}}_{v}^{{\prime} }$ is an ill-conditioned matrix (the condition number is >10¹⁰ for a typical ALMA spectral window), making the matrix inversion numerically unstable. The issue can be avoided by channel binning the data by a factor of two, as discussed in the following section.

C.2.2. Averaged Correlated Channels

As Hanning smoothing inflates the effective channel width by a factor of 8/3, it is very common to bin data across channels by a factor of 2; this is the default for ALMA observations. Here we calculate the appropriate covariance matrices for binning factors of 2, 3, and 4. A similar method can be followed for other binning factors, but we note that past a binning factor of 4, the covariance matrix and its inverse rapidly approach the uncorrelated case and the channel correlation can likely be neglected without significant adverse effect.

With a binning factor of 2 applied to data that have been Hanning smoothed, the data (indexed by channel in the averaged data set) can be described as:

$$\begin{eqnarray}&&{x}_{i,\mathrm{bin}\times 2}^{{\prime} }\\ &&=\,\displaystyle \frac{\left(\tfrac{1}{4},{x}_{i-0.5},+,\tfrac{1}{2},{x}_{i},+,\tfrac{1}{4},{x}_{i+0.5},),+,(,\tfrac{1}{4},{x}_{i},+,\tfrac{1}{2},{x}_{i+0.5},+,\tfrac{1}{4},{x}_{i+1}\right)}{2}\\ &&=\,\displaystyle \frac{1}{8}{x}_{i-0.5}+\displaystyle \frac{3}{8}{x}_{i}+\displaystyle \frac{3}{8}{x}_{i+0.5}+\displaystyle \frac{1}{8}{x}_{i+1},\end{eqnarray} \tag{ 43 }$$

and the noise component is therefore:

$$\begin{eqnarray}&&{v}_{i,\mathrm{bin}\times 2}^{{\prime} }=\displaystyle \frac{1}{8}{v}_{i-0.5}+\displaystyle \frac{3}{8}{v}_{i}+\displaystyle \frac{3}{8}{v}_{i+0.5}+\displaystyle \frac{1}{8}{v}_{i+1}.\end{eqnarray} \tag{ 44 }$$

Substituting this into Equation (35) and applying Equation (37), we find that for diagonal elements of ${{\boldsymbol{R}}}_{v}^{{\prime} }$

$$\begin{eqnarray}&&{R}_{v,\mathrm{bin}\times 2}^{{\prime} }[j,j]=\displaystyle \frac{{\sigma }_{j-0.5}^{2}}{64}+\displaystyle \frac{9{\sigma }_{j}^{2}}{64}+\displaystyle \frac{9{\sigma }_{j+0.5}^{2}}{64}+\displaystyle \frac{{\sigma }_{j+1}^{2}}{64},\end{eqnarray} \tag{ 45 }$$

and the only populated off-diagonal elements are

$$\begin{eqnarray}&&{R}_{v,\mathrm{bin}\times 2}^{{\prime} }[j,j\pm 1]=\displaystyle \frac{3{\sigma }_{j\pm 0.5}^{2}}{64}+\displaystyle \frac{3{\sigma }_{j\pm 1}^{2}}{64}.\end{eqnarray} \tag{ 46 }$$

Assuming the weights are roughly equivalent across channels for a given visibility, ${{\boldsymbol{R}}}_{v}^{{\prime} }$ can be approximated as

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 2}^{{\prime} }\approx {\sigma }_{i}^{2}\left[\begin{array}{cccc}\tfrac{5}{16} & \tfrac{3}{32} & & \\ \tfrac{3}{32} & \tfrac{5}{16} & \ddots & \\ & \ddots & \ddots & \tfrac{3}{32}\\ & & \tfrac{3}{32} & \tfrac{5}{16}\end{array}\right]\\ \quad \,=\,{\sigma }_{i}^{{\prime} 2}\left[\begin{array}{cccc}1 & \tfrac{3}{10} & & \\ \tfrac{3}{10} & 1 & \ddots & \\ & \ddots & \ddots & \tfrac{3}{10}\\ & & \tfrac{3}{10} & 1\end{array}\right]\,=\,{\sigma }_{i}^{{\prime} 2}\,[{{\boldsymbol{M}}}_{\mathrm{bin}\times 2}]\end{eqnarray} \tag{ 47 }$$

where ${\sigma }_{i}^{{\prime} }=\sqrt{5/16}{\sigma }_{i}$. Correspondingly,

$$\begin{eqnarray}&&{{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 2}^{{\prime} }}^{-1}\approx \displaystyle \frac{1}{{\sigma }_{i}^{{\prime} 2}}\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 2}]}^{-1}={w}_{i}^{{\prime} }\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 2}]}^{-1}.\end{eqnarray} \tag{ 48 }$$

where the weights ${w}_{i}^{{\prime} }$ were calculated from the binned data using a task such as statwt. In contrast to the unbinned case, ${{\boldsymbol{R}}}_{v}^{{\prime} }$ is now tridiagonal and well-conditioned for inversion.

Repeating these calculations for binning factors of 3 and 4, we find that:

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 3}^{{\prime} }\approx {\sigma }_{i}^{2}\left[\begin{array}{cccc}\tfrac{1}{4} & \tfrac{1}{24} & & \\ \tfrac{1}{24} & \tfrac{1}{4} & \ddots & \\ & \ddots & \ddots & \tfrac{1}{24}\\ & & \tfrac{1}{24} & \tfrac{1}{4}\end{array}\right]\\ \quad =\,{\sigma }_{i}^{{\prime} 2}\left[\begin{array}{cccc}1 & \tfrac{1}{6} & & \\ \tfrac{1}{6} & 1 & \ddots & \\ & \ddots & \ddots & \tfrac{1}{6}\\ & & \tfrac{1}{6} & 1\end{array}\right]\,=\,{\sigma }_{i}^{{\prime} 2}\,[{{\boldsymbol{M}}}_{\mathrm{bin}\times 3}]\end{eqnarray} \tag{ 49 }$$

where ${\sigma }_{i}^{{\prime} }=\sqrt{1/4}{\sigma }_{i}$,

$$\begin{eqnarray}&&{{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 3}^{{\prime} }}^{-1}\approx \displaystyle \frac{1}{{\sigma }_{i}^{{\prime} 2}}\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 3}]}^{-1}=w{{\prime} }_{i}\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 3}]}^{-1},\end{eqnarray} \tag{ 50 }$$

and

$$\begin{eqnarray}{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 4}^{{\prime} }\approx {\sigma }_{i}^{2}\left[\begin{array}{cccc}\tfrac{13}{64} & \tfrac{3}{128} & & \\ \tfrac{3}{128} & \tfrac{13}{64} & \ddots & \\ & \ddots & \ddots & \tfrac{3}{128}\\ & & \tfrac{3}{128} & \tfrac{13}{64}\end{array}\right]\\ \quad =\,{\sigma }_{i}^{{\prime} 2}\left[\begin{array}{cccc}1 & \tfrac{3}{26} & & \\ \tfrac{3}{26} & 1 & \ddots & \\ & \ddots & \ddots & \tfrac{3}{26}\\ & & \tfrac{3}{26} & 1\end{array}\right]\,=\,{\sigma }_{i}^{{\prime} 2}\,[{{\boldsymbol{M}}}_{\mathrm{bin}\times 4}]\end{eqnarray} \tag{ 51 }$$

where ${\sigma }_{i}^{{\prime} }=\sqrt{13/64}{\sigma }_{i}$,

$$\begin{eqnarray}&&{{{\boldsymbol{R}}}_{v,\mathrm{bin}\times 4}^{{\prime} }}^{-1}\approx \displaystyle \frac{1}{{\sigma }_{i}^{{\prime} 2}}\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 4}]}^{-1}=w{{\prime} }_{i}\,{[{{\boldsymbol{M}}}_{\mathrm{bin}\times 4}]}^{-1}.\end{eqnarray} \tag{ 52 }$$