Modeling the Echelle Spectra Continuum with Alpha Shapes and Local Regression Fitting

The proposed AFS algorithm is a versatile algorithm that can be used to remove the blaze function whether or not a corresponding laboratory source spectrum is available. In this algorithm, an alpha shape is used to obtain a preliminary estimation of the blaze function's high-level shape. An alpha shape is a generalization of a convex hull, but is not required to be a convex set; it is a region bounded by a set of segments generated from a set of points. To generate an alpha shape, imagine a piece of paper with a plot of spectrum printed on it. Consider a special paper cutter that can only cut out circles with radius α, known as α-balls. An incomplete circle is allowed, but the cutter cannot cut anything from the spectrum. In Figure 1(a), the blue circle is an α-ball that can be cut out from the paper, but we cannot move the α-ball any lower vertically because it would cut the spectrum. Continue to cut as many α-balls as possible, and the remaining paper is called an alpha hull. Figure 1(a) shows the resulting alpha hull with α = 5. By connecting the points where the alpha hull touches the spectrum with straight segments, it becomes an alpha shape, as displayed in Figure 1(b). An alpha shape can capture the general shape of a spectrum, and its upper boundary is used as a starting model for the blaze function estimation.

Zoom In Zoom Out Reset image size

Another technique used in the AFS algorithm is local polynomial regression (Cleveland 1979), which is a nonparametric method for estimating functions given a set of points. At each point, (λ_i, y_i), a p-degree polynomial model is fitted to a subset of neighboring points of λ_i. The p-degree polynomial regression is fitted by weighted least squares, giving more weight to points closer to λ_i and less weight to points further away. Consider a data set $({\lambda }_{1},{y}_{1}),({\lambda }_{2},{y}_{2}),\,...,\,({\lambda }_{n},{y}_{n})$, and a p-degree polynomial function in the neighborhood of λ_i is ${f}_{{\beta }^{({\lambda }_{i})}}({\lambda }_{j})={\beta }_{0}^{({\lambda }_{i})}\,+$ ${\beta }_{1}^{({\lambda }_{i})}({\lambda }_{j}-{\lambda }_{i})\,+$ ${\beta }_{2}^{({\lambda }_{i})}{\left({\lambda }_{j}-{\lambda }_{i}\right)}^{2}+...\,+$ ${\beta }_{p}^{({\lambda }_{i})}{\left({\lambda }_{j}-{\lambda }_{i}\right)}^{p}$, where λ_j is in the neighborhood of λ_i. Then the estimation of ${\beta }^{({\lambda }_{i})}$, ${\hat{\beta }}^{({\lambda }_{i})}=\left({\hat{\beta }}_{0}^{({\lambda }_{i})},{\hat{\beta }}_{1}^{({\lambda }_{i})},\,...,\,{\hat{\beta }}_{p}^{({\lambda }_{i})}\right)$, is obtained by minimizing

$$\begin{eqnarray}&&\displaystyle \sum _{{\lambda }_{j}\in {N}_{{m}_{0}}({\lambda }_{i})}{\omega }_{i}({\lambda }_{j}){\left({y}_{j}-{f}_{{\beta }^{({\lambda }_{i})}}({\lambda }_{j})\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where ${N}_{{m}_{0}}({\lambda }_{i})$ is the neighboring set of λ_i containing the nearest $\lfloor {m}_{0}n\rfloor$ pixels to λ_i, m₀ is a smoothing parameter, and ${\omega }_{i}({\lambda }_{j})=K\left(\tfrac{| {\lambda }_{j}-{\lambda }_{i}| }{\mathop{\max }\limits_{{\lambda }_{{j}^{* }}\in {N}_{{m}_{0}}({\lambda }_{i})}| {\lambda }_{{j}^{* }}-{\lambda }_{i}| }\right)$, with $K(x)={\left(1-{x}^{3}\right)}^{3}$ as a common weighting function. The local polynomial estimate at λ_i is ${\hat{\beta }}_{0}^{({\lambda }_{i})}$. An advantage of local polynomial regression over ordinary polynomial regression is its ability to adapt to local characteristics of a data set rather than fitting all the data points using one single model. This flexibility makes it useful for modeling complex data sets where regular polynomial regression fails.

Let ${\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}$ be an observed spectrum, where λ_i is the wavelength of pixel i, and y_i is the intensity of pixel i. Our method is summarized in Algorithm 1 and illustrated in Figure 2 with the details of the AFS algorithm presented next.

Zoom In Zoom Out Reset image size

In step 1, the intensity vector, $y=({y}_{1},{y}_{2},\,...,{y}_{n})$, is rescaled by multiplying by a value $u=\tfrac{\max (\lambda )-\min (\lambda )}{10\times \max (y)}$. This u corresponds to the α value for the alpha shape. Because the construction of an alpha shape depends on coverage areas of small circles, the relative range of λ and y affects results: if the range of y is too large or too small compared to the range of λ, larger alpha values should be used for the alpha shape. For a common blaze function, $u=\tfrac{\max (\lambda )-\min (\lambda )}{10\times \max (y)}$ works well by scaling the range of y and λ to be 1: 10, in coordination with a recommended value for α later. In step 2, the alpha shape, AS_α, is constructed with radius α, which is an infinite point set containing all the points within the boundary of the alpha shape (see Figure 2(a)). In AS_α, for each λ_i, there are infinite ${y}_{i}^{* }$ such that $({\lambda }_{i},{y}_{i}^{* })\in {\mathrm{AS}}_{\alpha }$. The upper boundary of AS_α is defined as ${\widetilde{\mathrm{AS}}}_{\alpha }$, which is a finite point set only including the largest ${y}_{i}^{* }$ such that $({\lambda }_{i},{y}_{i}^{* })\in {\mathrm{AS}}_{\alpha }$ for each λ_i, as displayed in Figure 2(a). In step 3, we fit a local polynomial regression using all the points in ${\widetilde{\mathrm{AS}}}_{\alpha }$, denoted as ${\hat{B}}_{1}$, such that ${\hat{B}}_{1}$ is a smoothed version of ${\widetilde{\mathrm{AS}}}_{\alpha }$ and is the initial model of the blaze function (see Figure 2(b)). ${\hat{B}}_{1}$ is not an accurate estimate of the blaze function, but rather an approximation of its shape. Next, y is divided by ${\hat{B}}_{1}$ to get the first estimate of the flattened spectrum, denoted as ${\hat{y}}^{(1)}$, displayed in Figure 2(c).

In step 4, we denote the intersection of ${\widetilde{\mathrm{AS}}}_{\alpha }$ and the spectrum ${\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}$ as W_α, shown in Figures 2(b) and (c). Notice that W_α contains points mostly near the continuum. The blue line in Figure 2(b), which connects the points in ${\widetilde{\mathrm{AS}}}_{\alpha }$, can be thought of as a collection of segments, and points in W_α are vertices of those segments. Each point in W_α is a vertex of a window where splits are defined—the spectrum is cut into small windows by the points in W_α in order to get the local quantiles. The purpose of the next steps is to select a subset of points that do not fall into an absorption feature. This is accomplished by selecting ${\hat{y}}_{i}^{(1)}$ in the upper quantiles of these windows. In particular, in the jth window, we select the points whose ${\hat{y}}^{(1)}$ values are larger than the q quantile of the ${\hat{y}}^{(1)}$ in the window, and the selected points make up the set ${S}_{j,\alpha ,q}$. For $j=1,2,\,...,\,| {W}_{\alpha }| -1$, where $| {W}_{\alpha }|$ is the number of points in W_α, the combined set of ${S}_{j,\alpha ,q}$ for different j's is defined as S_α,q, displayed in Figure 2(c). S_α,q contains points that are locally in the upper 1 − q quantile, which will be used to estimate the blaze function because these points generally do not fall into an absorption line. In step 5, we run a local polynomial regression on S_α,q and fit it to the whole spectrum. This regression is our final estimate of the blaze function, denoted as ${\hat{B}}_{2}$. The red line in Figure 2(d) shows the final estimate of the blaze function. Then, in step 6, y is divided by ${\hat{B}}_{2}$ to get the blaze-removed spectrum, denoted as ${\hat{y}}^{(2)}$ and shown in Figure 2(e).

Algorithm 1. AFS Algorithm

Step 0: Let ${\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}$ be an observed spectrum.

Step 1: Let $u=\tfrac{\max (\lambda )-\min (\lambda )}{10\times \max (y)}$. Multiply y by u.

Step 2: Let ${\mathrm{AS}}_{\alpha }={alpha}\,{shape}(\left\{({\lambda }_{i},{y}_{i}),i\,=\,1,...,n\right\})$ with radius value α. Then ${\widetilde{\mathrm{AS}}}_{\alpha }=\left\{({\lambda }_{i},\tilde{y}({\lambda }_{i})):{\lambda }_{i}\in \left\{{\lambda }_{i},i\,=\,1,...,n\right\},\,\tilde{y}({\lambda }_{i})=\mathop{\max }\limits_{\forall ({\lambda }_{i},{y}_{i}^{* })\in {\mathrm{AS}}_{\alpha }}{y}_{i}^{* }\right\}$.

Step 3: Run a local polynomial regression on ${\widetilde{\mathrm{AS}}}_{\alpha }$ with smoothing parameter m0, denoting the fit model as ${\hat{B}}_{1}$. Calculate ${\hat{y}}^{(1)}=\tfrac{y}{{\hat{B}}_{1}}$.

Step 4: Let ${W}_{\alpha }={\widetilde{\mathrm{AS}}}_{\alpha }\cap {\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}=\left\{({\lambda }_{i},{y}_{i}),i={w}_{1},{w}_{2},...,{w}_{| {W}_{\alpha }| }\right\}$. Let ${S}_{j,\alpha ,q}=\left\{{w}_{j}\leqslant i\leqslant {w}_{j+1}:\tfrac{{\sum }_{k={w}_{j}}^{{w}_{j+1}}{\mathbb{1}}({\hat{y}}_{i}^{(1)}\geqslant {\hat{y}}_{k}^{(1)})}{{w}_{j+1}-{w}_{j}+1}\geqslant q\right\}$. The ${S}_{\alpha ,q}=\bigcup _{j=1,...,| {W}_{\alpha }| -1}{S}_{j,\alpha ,q}$.

Step 5: Run a local polynomial regression on set ${\left\{({\lambda }_{i},{y}_{i})\right\}}_{i\in {S}_{\alpha ,q}}$ with m0 and fit to the whole spectrum, denoted as ${\hat{B}}_{2}$.

Step 6: Calculate ${\hat{y}}^{(2)}=\tfrac{y}{{\hat{B}}_{2}}$. Output ${\left\{({\lambda }_{i},{\hat{y}}_{i}^{(2)})\right\}}_{i=1}^{n}$.

Download table as:  ASCIITypeset image

The ALSFS algorithm is a method for removing the blaze function when a laboratory source spectrum, such as an LED or quartz lamp spectrum, is available as a reference. Generally, when a reliable reference spectrum is available, it can be used as the preliminary estimate of the blaze function shape. The use of a reference spectrum can be particularly advantageous in situations when the science spectrum contains wide absorption lines, which often make estimation of blaze function shape especially challenging.

The AFS algorithm in Section 2.1 can be adapted to take advantage of this additional information in the reference spectrum. To a first approximation, a laboratory source continuum spectrum should trace the instrumental blaze function for each order. However, the calibration source will typically have some effective blackbody temperature—that is, its intrinsic intensity will peak at a specific wavelength. Over the limited wavelength range covered by a single spectral order, this effective blackbody function is approximately linear. Therefore, an observation of a laboratory source spectrum can be modeled as a linear transformation of the instrumental blaze function, and the difference between the blaze function and the corresponding laboratory source spectrum is only a location-scale transformation. Thus, the blaze estimation problem can be translated into an optimization problem to find the best intercept, scale, and slope of the laboratory source.

2.2.1. ALSFS Algorithm

The ALSFS algorithm is initialized in the same manner as the AFS algorithm in steps 1 and 2, but it differs slightly in steps 3 and 4 and greatly in step 5. Let the reference laboratory source spectrum be ${\left\{({\lambda }_{i},{l}_{i})\right\}}_{i=1}^{n}$. The ALSFS algorithm is summarized in Algorithm 2.

Steps 1 and 2 are the same as the AFS algorithm. Because prior knowledge of the shape of blaze function is available, fewer points are needed for the second local polynomial fitting. Step 3 is similar to step 3 of the AFS algorithm, but we make set S_α,q more accurate by also calculating the 2q − 1 quantile of ${\hat{y}}^{(1)}$, denoted as ${Q}_{2q-1}$. ${Q}_{2q-1}={Q}_{1-2(1-q)}$, which means the upper 2(1 − q) quantile of ${\hat{y}}^{(1)}$, displayed in Figure 3(a). Compared to the upper 1 − q quantile of ${\hat{y}}^{(1)}$ within each window, a smaller quantile is used for the global quantile; otherwise, too many points will be excluded if the same quantile is used. Step 4 also departs from the AFS algorithm: for the jth window, we select the points λ_i where both ${\hat{y}}_{i}^{(1)}\geqslant$ the q quantile of the window and ${\hat{y}}_{i}^{(1)}\geqslant {Q}_{2q-1}$ into set S_α,q, shown in Figure 3(a). In step 5, we use the laboratory source spectrum's intensity curve as a reference model, displayed in Figure 3(b), to find the best linear coefficients for blaze function estimation using the points selected in the previous step. We apply a linear transformation on the laboratory source spectrum: ${\hat{l}}_{i}(a,b,c)=a+{{bl}}_{i}+c{\lambda }_{i}$, where a, b, and c are intercept, scale, and slope parameters, respectively. Because our ultimate goal is to have a flat spectrum without the blaze function, we use an objective function ${\sum }_{i\in {S}_{\alpha ,q}}{\left(\tfrac{{l}_{i}}{{\hat{l}}_{i}(a,b,c)}-1\right)}^{2}$, which measures the total distances from the removed spectrum to the constant 1. We minimize this objective function on the set S_α,q to get estimates for a, b, and c. Then the modified laboratory source spectrum $\hat{a}+\hat{b}{l}_{i}+\hat{c}{\lambda }_{i}$ is our final estimate for blaze function, displayed in Figure 3(b). Step 6 is the same as the AFS algorithm, y is divided by ${\hat{B}}_{2}$ to get the blaze-removed spectrum.

Zoom In Zoom Out Reset image size

Algorithm 2. ALSFS Algorithm

Step 0: Let ${\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}$ be an observed spectrum, and ${\left\{({\lambda }_{i},{l}_{i})\right\}}_{i=1}^{n}$ be the corresponding laboratory source.

Step 1: Let $u=\tfrac{\max (\lambda )-{\min }(\lambda )}{10\times \max (y)}$. Multiply y by u.

Step 2: Let ${\mathrm{AS}}_{\alpha }={alpha}\,{shape}(\left\{({\lambda }_{i},{y}_{i}),i\,=\,1,...,n\right\})$ with radius value α. and ${\widetilde{\mathrm{AS}}}_{\alpha }=\left\{({\lambda }_{i},\tilde{y}({\lambda }_{i})):{\lambda }_{i}\in \left\{{\lambda }_{i},i\,=\,1,...,n\right\},\,\tilde{y}({\lambda }_{i})=\mathop{\max }\limits_{\forall ({\lambda }_{i},{y}_{i}^{* })\in {\mathrm{AS}}_{\alpha }}{y}_{i}^{* }\right\}$.

Step 3: Run a local polynomial regression on ${\widetilde{\mathrm{AS}}}_{\alpha }$ with m0, denoted as ${\hat{B}}_{1}$. Calculate ${\hat{y}}^{(1)}=\tfrac{y}{{\hat{B}}_{1}}$. Denote ${Q}_{2q-1}={quantile}({\hat{y}}^{(1)},2q-1)$.

Step 4: Let ${W}_{\alpha }={\widetilde{\mathrm{AS}}}_{\alpha }\cap {\left\{({\lambda }_{i},{y}_{i})\right\}}_{i=1}^{n}=\left\{({\lambda }_{i},{y}_{i}),i={w}_{1},{w}_{2},...,{w}_{| {W}_{\alpha }| }\right\}$. Let ${S}_{j,\alpha ,q}=\left\{{w}_{j}\leqslant i\leqslant {w}_{j+1}:\tfrac{{\sum }_{k={w}_{j}}^{{w}_{j+1}}{\mathbb{1}}({\hat{y}}_{i}^{(1)}\geqslant {\hat{y}}_{k}^{(1)})}{{w}_{j+1}-{w}_{j}+1}\geqslant q\ \,{and}\,{\hat{y}}_{i}^{(1)}\geqslant {Q}_{2q-1}\right\}$. ${S}_{\alpha ,q}=\bigcup _{j\,=\,1,...,| {W}_{\alpha }| -1}{S}_{j,\alpha ,q}$.

Step 5: Consider a linear transformation: ${\hat{l}}_{i}(a,b,c)=a+{{bl}}_{i}+c{\lambda }_{i}$, $i=1,...,n$. $(\hat{a},\hat{b},\hat{c})={\mathrm{argmin}}_{a,b,c}{\sum }_{i\in {S}_{\alpha ,q}}{\left(\tfrac{{l}_{i}}{{\hat{l}}_{i}(a,b,c)}-1\right)}^{2}$. ${\hat{B}}_{2}=\hat{a}+\hat{b}{l}_{i}+\hat{c}{\lambda }_{i}$, $i=1,...,n$.

Step 6: Calculate ${\hat{y}}^{(2)}=\tfrac{y}{{\hat{B}}_{2}}$. Output ${\left\{({\lambda }_{i},{\hat{y}}_{i}^{(2)})\right\}}_{i=1}^{n}$.

Download table as:  ASCIITypeset image

2.2.2. Laboratory Source Smoothing by AFS Algorithm

The process of flat-fielding spectra is necessary for removing pixel-to-pixel QE variations in charge-coupled devices that are used as detectors in astronomical spectrographs. In the case of fiber-fed echelle spectrographs, flat-fielding can be carried out by extracting a featureless calibration spectrum and dividing the extracted science spectrum order-by-order. However, the flat-field source is generally not perfectly uniform in intensity over a large wavelength range and, therefore, will typically have an effective blackbody temperature that does not match the stellar effective temperature. As a result, this division leaves behind residual trends. By using the AFS algorithm, a model can be fitted to each order of the flat-field calibration spectrum. Division of the flat-field echelle orders by this fitted model will then produce a normalized spectrum that can be used to divide out the QE variation. Figure 4 shows an example of how this process was used to create a normalized flat in red orders that exhibit fringing from interference of red wavelengths in thinned silicon detectors. This fringing can be removed by dividing stellar spectra with this normalized flat.

Zoom In Zoom Out Reset image size

In some orders, a cosmic ray or a pixel with very low QE will create an upward or downward spike that is confined to one or a few pixels. In this case, the AFS algorithm is slightly modified to iteratively reject these pixels using outlier rejection. Let the original laboratory source spectrum be ${\left\{({\lambda }_{i},{L}_{i})\right\}}_{i=1}^{n}$, and ${\rm{\Delta }}L=\left\{| {L}_{i}-{L}_{i-1}| ,i=2,\,...,\,n\right\}$. Let ${Q}_{{q}_{s}}$ be the q_s quantile of ΔL, where q_s could be a number between 0.95 and 0.99. In the beginning, the 0.99 quantile of ΔL is denoted as ${Q}_{0.99}^{(0)}$. In the jth iteration, remove pixels where $| {L}_{i}-{L}_{i-1}| \gt {Q}_{0.99}^{(j-1)}$ and calculate the 0.99 quantile of the new ΔL, denoted as ${Q}_{0.99}^{(j)}$. Continue the iteration until ${Q}_{0.99}^{(j)}\lt {Q}_{{q}_{s}}$. The remaining pixels are used for smoothing by the AFS algorithm, displayed in Figure 5.

Zoom In Zoom Out Reset image size

In the proposed algorithms, there are several parameters that need to be selected by users. In the AFS algorithm, there are three parameters: α for the alpha shape, quantile q for the point selection, and m₀ in two local polynomial regressions. The α determines how many windows a spectrum is cut into, because the number of windows is determined by smoothness of the alpha shape, which is controlled by α. Using these windows, q determines how many points are selected in each window. After selecting the points, m₀ determines the smoothness of local polynomial fitting on these selected points. Practically, the three parameters are robust within appropriate ranges. For the example spectrum in Figure 6(a), we recommend a set of standard parameters $(\alpha =\tfrac{1}{6}\times$ wavelength range, q = 0.95, m₀ = 0.25) as a default. In Figures 6–8, we show the blaze estimates using extreme parameter choices (top panels) and its resulting normalized spectra compared with the standard parameter choice (bottom panels).

Zoom In Zoom Out Reset image size

First, the α can be selected according to the shape of a blaze function. For example, in Figure 2(a), we find that the blaze function increases first and then decreases. Also, it is convex first, then becomes concave, and turns to be convex in the end. We want to select an α that can capture the convex parts (not too large), but will not go too deep into absorption lines (not too small). The choice of α is mainly determined by the shape, e.g., curvature and concavity, of a blaze function. However, when there is a wide absorption line, a larger α may be needed than the shape would suggest. For echelle spectra orders, because each convex portion generally takes up about one-sixth of an order, we recommend selecting an α that is one-sixth of the wavelength range of the order, but have found empirically that α is rather robust from $\tfrac{1}{12}$ of the wavelength range to one-third of the wavelength range. Figure 6(a) shows an example of a very large α equal to the order's entire wavelength range. This α value is too large for the α-balls to capture the shape of the spectrum near the boundaries; in the blaze-removed spectrum, the left part drops downward like a wide absorption feature because points in that portion of the spectrum were not selected in S_α,q. Figure 6(b) shows an example of a small α of $\tfrac{1}{50}$ of the wavelength range. With such a small α, the α-balls fit into absorption lines so that points inside an absorption are selected into set S_α,q. In the removed spectrum, there are regions that are above the reference line y = 1, compared with the cyan spectrum where $\alpha =\tfrac{1}{6}\times$ wavelength range.

The parameter q depends on the signal-to-noise ratio (S/N) of a spectrum and the amount of absorption. The goal when selecting q is to find points on the spectrum that do not drop into absorption lines, but instead are on the continuum. After flattening the spectrum using a preliminary estimate of its shape, we select points in the local upper 1−q quantile for the set S_α,q to be used in the final estimation. If we happened to know the true blaze, then after dividing the spectrum by the blaze we would expect to see points randomly scattered around 1. In the proposed algorithms, these points are approximated by set S_α,q. If the S/N is high or there is a large amount of absorption, a larger q is needed to select points in S_α,q so that these points do not fall in absorption lines. Conversely, if S/N is low or there is minimal absorption, a smaller q is needed to get enough points into set S_α,q.⁴ For an order with a similar amount of absorption as the one displayed in Figure 6(a), a q from 0.95 to 0.99 works for S/N 300, a q from 0.85 to 0.95 works for S/N 150, and a q from 0.5 to 0.85 works for S/N lower than 150. Figure 7(a) shows the effects of a small q such as 0.7: too many points are selected into S_α,q so that the blaze function estimate is dragged downward by the points in absorption lines, and in the removed spectrum, the left and right boundary regions are above the reference line y = 1. In contrast, Figure 7(b) shows a large q of 0.999. In this example S_α,q contains too few points, which makes the estimate sit almost completely above the spectrum. In the removed spectrum, almost all the pixels are under the reference line.

Zoom In Zoom Out Reset image size

The smoothing parameter m₀ depends on the distribution of points in S_α,q along the spectrum, which is determined by the amount of absorption. If there are many absorption lines or any absorption lines that are wide, the set S_α,q has large gaps between pixels, and, thus, a large m₀ is needed to get a good estimate. If there are few absorption lines or absorption lines that are narrow, a small m₀ is needed so that the estimation better adapts to local regions. For an order with a similar amount of absorption as the one displayed in Figure 6(a), an m₀ value from 0.15 to 0.3 has worked well empirically. In Figure 8(a), m₀ is set to be 0.5 (too large) and the estimate is off on the left part of the spectrum, where the blaze-removed spectrum rises to above 1.5. Figure 8(b) shows the results of a too small m₀ value of 0.1: the blaze estimate has some small bumps, but the blaze-removed spectrum looks reasonable to the eye. However, because the true blaze function does not have small bumps, this blaze estimate is not as good a fit as it might appear at first glance.

Zoom In Zoom Out Reset image size

In general, the blaze function estimate is more sensitive to small changes in q than to α and m₀. For an echelle spectrum order with a similar amount of absorption as the one displayed in Figure 6(a), we can start with an α equal to one-sixth of the wavelength range of the order, an m₀ equal to 0.25, and tune the parameter q within the range according to its S/N and amount of absorption. A more detailed set of recommendations for parameters in different situations is provided in the Appendix.

The ALSFS algorithm has the same parameters: α, q, and m₀. Because the final estimate is the reference laboratory source spectrum with linear modification, m₀ has less influence on the results than for the AFS algorithm. In the laboratory source spectrum smoothing process, α, q, and m₀ operate the same as in the AFS algorithm. The quantile parameter q_s in ${Q}_{{q}_{s}}$ depends on the particular appearances of spikes. If spikes are long, use a smaller q_s such as 0.95; if spikes are very small, use a larger q_s such as 0.98 or 0.99.

We have found that the proposed algorithms work well in most cases. However, there are several special cases that can result in poorer estimates of the blaze function, for which we have developed corrections to mitigate these issues. Because the AFS algorithm relies only on the science spectrum itself, it is more susceptible to complications than the ALSFS algorithm.

2.4.1. Boundary Correction

An order is normalized by dividing by its estimated blaze function. Because the blaze function approaches zero near the edges of the order, small errors in the blaze shape will be magnified in the divided spectrum. The ALSFS algorithm is less susceptible to this problem because of the strong constraints provided by the laboratory source, but other blaze estimation methods, including the AFS algorithm, can be strongly affected.

The AFS algorithm also has difficulty in this region in cases where the edge of an order splits an absorption line. Fortunately, neighboring orders often have some region of overlap that can be used to correct the boundaries. A weighted average of the blaze-removed spectrum of the two orders can be used as an estimate of the blaze function in the overlapping region. For example, Figure 9 shows two neighboring orders that share an overlapping region. Figure 9(a) shows the right boundary of the left order, which looks good as a blaze-removed spectrum using AFS algorithm. Figure 9(b) shows the left boundary of the right order, which spuriously rises above 1. We correct the overlapping region using the following:

$$\begin{eqnarray*}&&{y}_{\mathrm{corrected},{\rm{l}}}={w}_{l}{y}_{1,l}+(1-{w}_{l}){y}_{2,l},\,\,l=1,\,...,\,m,\end{eqnarray*}$$

where y₁ and y₂ are the intensities for overlapping regions from the left and right orders, respectively, y_corrected is the array of intensities for the corrected overlapping region, m is the number of pixels in the overlapping region, and ${w}_{l}=1-\tfrac{l}{m}$, for $l=1,\,...,\,m$. The result of the correction is shown in Figure 9(c). The boundary-corrected spectrum is much better than the original estimate, and we can achieve further improvement by changing the definition of w_l. For example, if it is known that one of the two orders has a better estimate on its boundary, we can assign more weight toward the better order. This might be the case for a pair of orders with a broad spectral feature that is cut off on only one of the orders.

Zoom In Zoom Out Reset image size

2.4.2. Wide Absorption Lines

Sometimes an order has wide absorption lines that influence the performance of the blaze function estimation. Wide absorption lines can significantly influence the AFS algorithm performance, but only slightly influence the ALSFS algorithm. Figure 10(b) shows the order containing the two deep Na D lines, which are each so broad that the spectrum does not fully return to the nominal continuum level between them. When attempting to fit the blaze function, the alpha shape will dip into wide features like these and pull the final blaze function estimate downward. In Figure 10(a), the spectrum (same as in Figure 10(b)) before the blaze function removal is displayed. Despite failing to return to the proper continuum level, this attempt at normalization looks reasonable to the eye. Without prior information, the proposed data-driven AFS algorithm cannot consistently determine the continuum level over regions of greatly extended absorption like this one. If we know there is a wide absorption feature before applying the algorithm, the regions can be masked in step 4 of the AFS algorithm, excluding those points from the estimate.

Zoom In Zoom Out Reset image size

2.4.3. Continuous Opacity

The iterative method is commonly used to remove blaze function from a spectrum. A polynomial model is fit to a spectrum order, and the fit is considered as the starting estimation for the blaze function. Next, the method iterates. In each iteration, there is a threshold curve obtained by:

$$\begin{eqnarray}&&\mathrm{threshold}({\lambda }_{i},t)={\mathrm{fit}}_{\mathrm{polynomial}}^{(t)}({\lambda }_{i})-\displaystyle \frac{0.5}{t+1},\end{eqnarray} \tag{ 2 }$$

where ${\mathrm{fit}}_{\mathrm{polynomial}}^{(t)}({\lambda }_{i})$ is the seventh-order polynomial regression estimate at λ_i and t is the iteration time. As t increases, the threshold curve increases. The method builds a subset M_t containing wavelength λ_i's with intensity larger than threshold(λ_i, t). Then in the next iteration, it only uses spectrum points whose wavelength values are in M_t to fit a polynomial model and the fit is a new estimation for continuum. To stop the iteration, a stopping time, denoted as T, can be set at a particular iteration. Another option is to set a standard deviation value (sd), such that the algorithm stops when the standard deviation of the blaze-removed spectrum is smaller than sd. The sd cannot be too small (much smaller than the standard deviation of the true spectrum without the blaze function), otherwise the iteration will not stop. In this work, we stop the iteration if either it arrives at the Tth iteration or the standard deviation of the blaze-removed spectrum is smaller than sd. Empirically, an sd value around 0.05 works well, but the performance of the iterative method is influenced by the choice of T. A higher S/N requires a larger T and, in general, we have found a T value equal to S/N seems to well.

The first order (blue) has wavelengths ranging from 4478 to 4528 Å, called "order B," and the second order (red) has wavelengths ranging from 6154 to 6221 Å, called "order R." For this simulation, we require a realistic blaze function to add to the NSO spectrum. To do this, we estimate a blaze function of a B-star, HR 5501, observed with EXPRES (Jurgenson et al. 2016). We first use the ALSFS algorithm, with the corresponding LED spectrum (after using the AFS algorithm on it) as the laboratory source, on the B-star spectrum to get the estimate of the blaze function. Then we apply this estimate as the blaze function to the two orders. The raw spectrum of the B-star spectrum is then used as the laboratory source reference for the ALSFS algorithm. The AFS algorithm, ALSFS algorithm, and iterative method are applied to estimate blaze functions of the two orders and residuals are calculated. The residuals of order B and R with S/N = 300 are displayed in Figures 11(a) and (b), respectively. Overall, ALSFS has the smallest residuals that are consistent across the whole order. AFS has larger residuals than ASLFS, which increase in the boundary regions, but the iterative method has larger residuals than AFS in both boundary regions and middle regions.

Zoom In Zoom Out Reset image size

Figure 11 displays a single realization of the noisy spectrum. We repeat the procedure by adding different realizations of noise 1000 times. For each of the 1000 realizations, the AFS algorithm, ALSFS algorithm, and the iterative method are applied to the two orders. This is carried out for an S/N 300, 150, and 50. The results are displayed in Figure 12 and Table 1. The rms error is calculated as $\sqrt{\tfrac{1}{n}{\sum }_{i=1}^{n}{r}_{i}^{2}}$, where r_i is the residual at pixel i. The ALSFS algorithm has the smallest median rms error for three scenarios: S/N 300-order B, S/N 300-order R, and S/N 150-order B. For the other three scenarios, the AFS algorithm has the smallest median rms error. Except for S/N 50-order R, the iterative method has the largest median rms error.

Zoom In Zoom Out Reset image size

An example of continuous opacity is illustrated using the same simulation setup: we apply a blaze function obtained from a B-star spectrum to the NSO spectrum, and then we attempt to recover the underlying spectrum. Here we examine five consecutive artificial orders. A region of continuous opacity spans about 80 Å, which is captured in parts of the third, fourth, and fifth orders, displayed in Figure 13(a). After applying the ALSFS algorithm, the resulting blaze-removed spectra are displayed in Figure 13(b). Although the blaze-removed spectra are flat, the ALSFS algorithm does not capture the continuous opacity correctly on its own. We know that the first two orders are not affected by the continuous opacity, and so they are used as the reference to correct the other orders: each order is linearly adjusted in intercept and slope to be align with its left neighboring order using the points within the overlapping region whose normalized intensities are in the top $\tilde{q}$ quantile. Ideally, the combined segment should recover the continuous opacity in the normalized intensities after the blaze function removal. However, error in the slope estimation of the first order is amplified: the combined spectrum in Figure 13(c) is not perfectly horizontal and it goes upward from left to right. We fit another linear regression using the points in the combined spectrum whose normalized intensities are in the top $\tilde{q}$ quantile to remove the extra slope. The resulting spectrum recovers the continuous opacity well, as displayed in Figure 13(d).

Zoom In Zoom Out Reset image size

The proposed methods can also be used for spectra with cosmic rays. In particular, the modified AFS algorithm for laboratory source smoothing, described in Section 2.2.2, can be used directly to deal with the presence of cosmic rays. To demonstrate this, an order with simulated cosmic rays is displayed in Figure 14(a), which contains two upward spikes. The blaze function estimate from the AFS algorithm is shown in Figure 14(a), and the blaze-removed spectrum is displayed in Figure 14(a) in comparison to the true spectrum. The ALSFS algorithm can be modified similarly to deal with cosmic rays.

Zoom In Zoom Out Reset image size