Automatic Kalman-Filter-based Wavelet Shrinkage Denoising of 1D Stellar Spectra

TL;DR

This paper introduces an adaptive Kalman-filter-based wavelet denoising method for 1D stellar spectra that improves noise reduction and peak preservation without requiring parameter tuning, outperforming traditional methods.

Contribution

It presents a novel adaptive Kalman thresholding approach combined with wavelet shrinkage, eliminating the need for choosing decomposition levels in denoising stellar spectra.

Findings

Superior noise suppression compared to common methods

Better peak shape preservation in denoised spectra

Effective for synthetic and potentially real survey data

Abstract

We propose a non-parametric method to denoise 1D stellar spectra based on wavelet shrinkage followed by adaptive Kalman thresholding. Wavelet shrinkage denoising involves applying the Discrete Wavelet Transform (DWT) to the input signal, `shrinking' certain frequency components in the transform domain, and then applying inverse DWT to the reduced components. The performance of this procedure is influenced by the choice of base wavelet, the number of decomposition levels, and the thresholding function. Typically, these parameters are chosen by `trial and error', which can be strongly dependent on the properties of the data being denoised. We here introduce an adaptive Kalman-filter-based thresholding method that eliminates the need for choosing the number of decomposition levels. We use the `Haar' wavelet basis, which we found to be the best-suited for 1D stellar spectra. We introduce various levels of Poisson noise into synthetic PHOENIX spectra, and test the performance of several common denoising methods against our own. It proves superior in terms of noise suppression and peak shape preservation. We expect it may also be of use in automatically and accurately filtering low signal-to-noise galaxy and quasar spectra obtained from surveys such as SDSS, Gaia, LSST, PESSTO, VANDELS, LEGA-C, and DESI.

Tables

Table 1. Table 1: Conversion between variables of the Kalman filter used in Appendix A and in § 3 in the main body of this paper.

Appendix	Main Text
${𝐭}_i$	${𝐱}_i$
${𝐧}_i^t$	${𝐰}_i$
${\widehat{𝐄}}_i{𝐭}_{i-1}$	${𝐇}_i$
${𝐦}_i$	${𝐳}_i$
${𝐧}_i^m$	${𝐯}_i$
${\widehat{𝐌}}_i{𝐭}_{i-1}$	${𝐌}_i$
${𝐭}_i^{\mathrm{pred}}$	${\widehat{𝐱}}_i^{-}$
${𝐓}_i^{\mathrm{pred}}$	${𝐏}_i^{-}$
${𝐓}_i^{\mathrm{upd}}$	${𝐏}_i$
${𝐍}_i^t$	${𝐐}_i$
${𝐭}_i^{\mathrm{upd}}$	${\widehat{𝐱}}_i$
${𝐖}_i$	${𝐊}_i$

Table 2. Table 2: The means and standard deviations for the error metrics for different denoising methods applied to the stellar absorption spectrum of Star 1 corrupted with Poisson noise to get a PSNR of 5.0 5.0 5.0 (the results for the other two stars show similar trends). The complete table, along with a detailed caption, is available online. See Figs. 10 – 12 (a) at the x–axis value of 5 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.113\pm 1.5\%$	$0.144\pm 1.3\%$	$0.595\pm 10.2\%$	$0.385\pm 1.8\%$
Moving Mean\parnotew = $(7,7,8,7)$	$0.048\pm 2.8\%$	$0.060\pm 2.1\%$	$0.246\pm 11.9\%$	$0.834\pm 0.8\%$
Exp. Smoothing\parnotespan = $(4,4,5,6)$	$\text{0.067}\pm \text{1.6\%}$	$\text{0.086}\pm \text{1.7\%}$	$\text{0.328}\pm \text{5.6\%}$	$\text{0.690}\pm \text{1.1\%}$
FFT, Gaussian kernel\parnote$\sigma$ = $(2,2,2,3)$	$\text{0.047}\pm \text{2.1\%}$	$\text{0.059}\pm \text{2.5\%}$	$\text{0.219}\pm \text{5.1\%}$	$\text{0.853}\pm \text{0.6\%}$
Wavelet + AKFThresh v1	$0.050\pm 1.7\%$	$0.064\pm 1.7\%$	$0.256\pm 8.6\%$	$0.817\pm 0.6\%$

Table 3. Table 3: Same as Table 2 , but using the tenth–best parameters for each denoising method. See Figs. 10 – 12 (b) at the x–axis value of 5 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.113\pm 1.1\%$	$0.145\pm 1.1\%$	$0.655\pm 14.2\%$	$0.382\pm 1.1\%$
Moving Mean\parnotew = $(13,8,12,4)$	$0.060\pm 1.2\%$	$0.076\pm 1.1\%$	$0.318\pm 3.6\%$	$0.726\pm 1.1\%$
Exp. Smoothing\parnotespan = $(11,11,12,13)$	$0.087\pm 1.4\%$	$0.108\pm 1.2\%$	$0.380\pm 3.4\%$	$0.654\pm 0.5\%$
FFT, Gaussian kernel\parnote$\sigma =(11,11,11,11)$	$0.098\pm 0.6\%$	$0.123\pm 0.4\%$	$\text{0.411}\pm \text{2.3\%}$	$\text{0.614}\pm \text{0.1\%}$
Wavelet + AKFThresh v1	$0.049\pm 2.4\%$	$0.063\pm 2.9\%$	$0.253\pm 6.5\%$	$0.815\pm 1.2\%$

Table 4. Table 4: Same as Table 3 , but using the twentieth–best parameters for the error metrics. See Figs. 10 – 12 (c) at the x–axis value of 5 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.112\pm 0.8\%$	$0.145\pm 0.9\%$	$0.636\pm 8.0\%$	$0.384\pm 0.9\%$
Moving Mean\parnotew = ($21,21,22,22$)	$0.084\pm 0.8\%$	$0.107\pm 0.6\%$	$0.392\pm 4.7\%$	$0.597\pm 0.6\%$
Exp. Smoothing\parnotespan = ($21,21,21,22$)	$0.108\pm 0.7\%$	$0.135\pm 0.5\%$	$0.413\pm 4.1\%$	$0.616\pm 0.7\%$
FFT, Gaussian kernel\parnote$\sigma =(21,21,21,21)$	$\text{0.122}\pm \text{0.4\%}$	$\text{0.150}\pm \text{0.3\%}$	$\text{0.478}\pm \text{1.8\%}$	$\text{0.567}\pm \text{0.1\%}$
Wavelet + AKFThresh v1	$0.050\pm 2.7\%$	$0.064\pm 2.4\%$	$0.267\pm 11.0\%$	$0.814\pm 0.8\%$

Table 5. Table 5: Same as Table 2 , but with PSNR of 10 10 10 . See Figs. 10 – 12 (a) at the x–axis value of 10 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.056\pm 1.4\%$	$0.072\pm 1.3\%$	$0.324\pm 4.3\%$	$0.662\pm 1.0\%$
Moving Mean\parnotew = ($5,5,4,7$)	$0.028\pm 1.6\%$	$0.035\pm 1.7\%$	$0.156\pm 7.9\%$	$0.922\pm 0.3\%$
Exp. Smoothing\parnotespan = ($3,3,3,4$)	$\text{0.042}\pm \text{1.6\%}$	$\text{0.054}\pm \text{1.8\%}$	$\text{0.225}\pm \text{12.4\%}$	$\text{0.827}\pm \text{0.6\%}$
FFT, Gaussian kernel\parnote$\sigma =(2,2,2,2)$	$0.028\pm 2.6\%$	$0.036\pm 2.9\%$	$0.152\pm 22.9\%$	$0.932\pm 0.3\%$
Wavelet + AKFThresh v1	$0.028\pm 2.1\%$	$0.037\pm 2.2\%$	$0.159\pm 11.1\%$	$0.911\pm 0.3\%$

Table 6. Table 6: Same as Table 3 , but with PSNR of 10 10 10 . See Figs. 10 – 12 (b) at the x–axis value of 10 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.056\pm 1.3\%$	$0.073\pm 1.1\%$	$0.330\pm 7.6\%$	$0.661\pm 0.7\%$
Moving Mean\parnotew = ($11,11,11,12$)	$0.044\pm 1.0\%$	$0.056\pm 0.8\%$	$0.235\pm 7.8\%$	$0.799\pm 0.2\%$
Exp. Smoothing\parnotespan = ($11,11,11,12$)	$0.081\pm 0.5\%$	$0.102\pm 0.5\%$	$0.328\pm 3.4\%$	$0.713\pm 0.3\%$
FFT, Gaussian kernel\parnote$\sigma =(11,11,11,11)$	$\text{0.096}\pm \text{0.4\%}$	$\text{0.121}\pm \text{0.2\%}$	$\text{0.395}\pm \text{1.5\%}$	$\text{0.616}\pm \text{0.1\%}$
Wavelet + AKFThresh v1	$0.028\pm 2.0\%$	$0.036\pm 1.9\%$	$0.157\pm 9.4\%$	$0.911\pm 0.3\%$

Table 7. Table 7: Same as Table 4 , but with PSNR of 10 10 10 . See Figs. 10 – 12 (c) at the x–axis value of 10 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.056\pm 0.8\%$	$0.072\pm 0.9\%$	$0.313\pm 8.7\%$	$0.664\pm 0.7\%$
Moving Mean\parnotew = ($21,21,21,21$)	$0.081\pm 0.4\%$	$0.104\pm 0.2\%$	$0.370\pm 5.7\%$	$0.614\pm 0.2\%$
Exp. Smoothing\parnotespan = ($21,21,21,21$)	$0.107\pm 0.5\%$	$0.131\pm 0.4\%$	$0.398\pm 2.4\%$	$0.632\pm 0.4\%$
FFT, Gaussian kernel\parnote$\sigma =(21,21,21,21)$	$\text{0.121}\pm \text{0.1\%}$	$\text{0.149}\pm \text{0.1\%}$	$\text{0.475}\pm \text{1.3\%}$	$\text{0.568}\pm \text{0.0\%}$
Wavelet + AKFThresh v1	$0.028\pm 1.5\%$	$0.036\pm 1.9\%$	$0.154\pm 7.3\%$	$0.913\pm 0.3\%$

Table 8. Table 8: Same as Tables 2 and 5 , but with PSNR of 20 20 20 . See Figs. 10 – 12 (a) at the x–axis value of 20 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.028\pm 1.2\%$	$0.036\pm 1.2\%$	$0.151\pm 10.8\%$	$0.869\pm 0.3\%$
Moving Mean\parnotew = ($3,5,3,5$)	$0.017\pm 1.2\%$	$0.021\pm 1.2\%$	$0.089\pm 4.2\%$	$0.969\pm 0.1\%$
Exp. Smoothing\parnotespan = ($2,2,2,2$)	$\text{0.024}\pm \text{1.0\%}$	$\text{0.031}\pm \text{0.9\%}$	$0.124\pm 8.7\%$	$\text{0.922}\pm \text{0.2\%}$
FFT, Gaussian kernel\parnote$\sigma =(2,2,2,2)$	$0.021\pm 1.2\%$	$0.027\pm 1.2\%$	$\text{0.132}\pm \text{14.5\%}$	$0.958\pm 0.1\%$
Wavelet + AKFThresh v1	$0.016\pm 2.4\%$	$0.020\pm 2.1\%$	$0.083\pm 6.9\%$	$0.966\pm 0.1\%$

Table 9. Table 9: Same as Tables 3 and 6 , but with PSNR of 20 20 20 . See Figs. 10 – 12 (b) at the x–axis value of 20 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.028\pm 1.1\%$	$0.036\pm 1.0\%$	$0.152\pm 7.9\%$	$0.868\pm 0.4\%$
Moving Mean\parnotew = ($11,11,11,11$)	$0.041\pm 0.6\%$	$0.053\pm 0.4\%$	$0.237\pm 5.1\%$	$0.842\pm 0.1\%$
Exp. Smoothing\parnotespan = ($11,11,11,11$)	$0.080\pm 0.5\%$	$0.100\pm 0.4\%$	$0.310\pm 1.9\%$	$0.727\pm 0.2\%$
FFT, Gaussian kernel\parnote$\sigma =(11,11,11,11)$	$\text{0.096}\pm \text{0.2\%}$	$\text{0.121}\pm \text{0.1\%}$	$\text{0.397}\pm \text{1.0\%}$	$\text{0.616}\pm \text{0.0\%}$
Wavelet + AKFThresh v1	$0.016\pm 1.4\%$	$0.020\pm 1.6\%$	$0.091\pm 9.0\%$	$0.965\pm 0.1\%$

Table 10. Table 10: Same as Tables 4 and 7 , but with PSNR of noisy signal 20 20 20 . See Figs. 10 – 12 (c) at the x–axis value of 20 for a graphical illustration of this table.

\parnotereset	${𝑳}^{\mathrm{𝟏}}$	${𝑳}^{\mathrm{𝟐}}$	${𝑳}^{\infty }$	$𝑺𝑺𝑰𝑴$
Noisy data	$0.028\pm 1.0\%$	$0.001\pm 0.9\%$	$0.165\pm 6.7\%$	$0.868\pm 0.2\%$
Moving Mean\parnotew = ($21,21,21,21$)	$0.080\pm 0.2\%$	$0.103\pm 0.1\%$	$0.354\pm 2.8\%$	$0.615\pm 0.1\%$
Exp. Smoothing\parnotespan = ($21,21,21,21$)	$0.106\pm 0.3\%$	$0.131\pm 0.2\%$	$0.380\pm 0.9\%$	$0.637\pm 0.2\%$
FFT, Gaussian kernel\parnote$\sigma =(21,21,21,21)$	$\text{0.121}\pm \text{0.1\%}$	$\text{0.149}\pm \text{0.0\%}$	$\text{0.476}\pm \text{0.8\%}$	$\text{0.568}\pm \text{0.0\%}$
Wavelet + AKFThresh v1	$0.016\pm 1.6\%$	$\text{0.020}\pm \text{1.5\%}$	$\text{0.082}\pm \text{7.2\%}$	$0.967\pm 0.2\%$

Equations

$${\tilde{w}}_{j,i}^{SOFT}=\{\begin{array}{ll}& \text{sgn}(w_{j,i})(|w_{j,i}|-{\lambda }_j),\mathrm{for}|w_{j,i}|\ge {\lambda }_j\\ & 0,\mathrm{for}|w_{j,i}|<{\lambda }_j.\end{array}$$

$${\tilde{w}}_{j,i}^{HARD}=\{\begin{array}{ll}& w_{j,i},\mathrm{for}|w_{j,i}|\ge {\lambda }_j\\ & 0,\mathrm{for}|w_{j,i}|<{\lambda }_j.\end{array}$$

$${\text{x}}_i$$

$${\text{z}}_i$$

$${\text{xˆ}}_i^{-}$$

$${\text{P}}_i^{-}$$

$${\text{K}}_i$$

$${\text{xˆ}}_i$$

$${\text{P}}_i$$

$${\mathbf{t}}_i={\hat{\mathbf{E}}}_i{\mathbf{t}}_{i-1}+{\mathbf{n}}_i^t.$$

$${\mathbf{m}}_i={\hat{\mathbf{M}}}_i{\mathbf{t}}_i+{\mathbf{n}}_i^m.$$

$${\mathbf{t}}_i^{\mathrm{pred}}={\hat{\mathbf{E}}}_i{\mathbf{t}}_{i-1}^{\mathrm{pred}}.$$

$${\mathbf{T}}_i^{\mathrm{pred}}={\hat{\mathbf{E}}}_i{\mathbf{T}}_{i-1}^{\mathrm{upd}}{\hat{\mathbf{E}}}_i+{\mathbf{N}}_i^t.$$

$${\mathbf{t}}_i^{\mathrm{upd}}={\mathbf{t}}_i^{\mathrm{pred}}+{\mathbf{W}}_i{\mathbf{\Delta }}_i^{m-p}.$$

$${\mathbf{\Delta }}_i^{\mathrm{m}-\mathrm{p}}\equiv {\mathbf{m}}_i-{\hat{\mathbf{M}}}_i{\mathbf{t}}_i^{\mathrm{pred}}.$$

$${\mathbf{W}}_i^{\mathrm{noiseless}m,i.e.e.sign.}\to {\hat{\mathbf{M}}}_i^{-1}.$$

$${\mathbf{M}}_i^{\mathrm{pred}}={\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T;$$

$${\mathbf{M}}_i={\mathbf{M}}_i^{\mathrm{pred}}+{\mathbf{N}}_i^m.$$

$${\mathbf{W}}_i\sim {\mathbf{M}}_i^{-1}={\left({\mathbf{M}}_i^{\mathrm{pred}}+{\mathbf{N}}_i^m\right)}^{-1}={\left({\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T+{\mathbf{N}}_i^m\right)}^{-1},$$

$${\mathbf{W}}_i^{\mathrm{noiseless}m}={\hat{\mathbf{M}}}_i^{-1}=\mathbf{A}{\left({\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T\right)}^{-1}.$$

$${\hat{\mathbf{M}}}_i^{-1}=\mathbf{A}\mathbf{B},$$

$${\mathbf{B}}^{-1}=\left({\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T\right),$$

$$\mathbf{A}={\hat{\mathbf{M}}}_i^{-1}\left({\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T\right).$$

$$\mathbf{A}={\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T.$$

$${\mathbf{W}}_i={\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T{\left({\hat{\mathbf{M}}}_i{\mathbf{T}}_i^{\mathrm{pred}}{\hat{\mathbf{M}}}_i^T\right)}^{-1}.$$

Full text

††pubyear: 2019††pagerange: Automatic Kalman-Filter-based Wavelet Shrinkage Denoising of 1D Stellar Spectra–17

1 Introduction

Most commonly-used methods for spectral denoising are parametric and thus require laborious user intervention for selection of appropriate parameters. The most popular among them is the Fourier transform (FT). Unlike other time-based techniques for noise removal such as the moving average (MA) and the exponential moving average (EMA), the FT provides a vision of a signal in the complementary domain (i.e., frequency domain, if the original signal is a time series). Because of the periodicity of the sinusoids used as the basis functions in the FT, it is well-suited to processing stationary signals. However, it has several drawbacks. It is inefficient when dealing with a non-stationary signal with varying shape and width, has zero resolution in the original (time) domain, and requires a user-specified cut-off frequency to estimate the noise in the signal. While the time-windowed FT offers an improvement by introducing a moving window that represents a compromise between time and frequency resolution, it is still limited by use of a fixed window size (Kaiser, 2010).

Wavelet transform (WT), on the other hand, is able to achieve an excellent trade-off between resolutions in time and frequency domains, by providing high (low) resolution in the frequency (time) domain for small frequency values, and high (low) resolution in the time (frequency) domain for large frequency values. This provides a clear advantage over FT and short time-windowed FT for exploring features of interest in the input signal in both time and frequency domains. However, traditional denoising approaches built using wavelet transform still require trial-and-error by the user to determine and suppress the noise in the signal, as we discuss further in §2 and §4.

In this work, we present a new denoising approach, based on wavelet decomposition (Graps, 1995) and adaptive Kalman-filter-based thresholding of the resultant components, which significantly improves denoising performance. It also provides a simpler, quicker implementation relative to previous methods, including FT, by being effectively non-parametric and hence independent of user input. For characterization of an object (e.g. a star, galaxy, or quasar) based on its spectrum, the precision and accuracy of the derived parameters depend strongly on two factors—the spectral resolution, and the signal-to-noise ratio (SNR). While the former is completely dependent on hardware (spectrograph), the latter can be controlled to some extent with software (algorithms). This dependence of object characterization on spectral SNR becomes even more prominent when dealing with low SNR spectra, as is common in large galaxy surveys. Our proposed method recursively decomposes the input stellar absorption spectrum—corrupted with Poisson (shot) noise—into high-frequency and low-frequency signals via a method called the wavelet transform. These high-frequency signals are then denoised (i.e. filtered) using Kalman filters. Kalman filters use a heuristic to derive noise statistics, and by employing a predictor-corrector scheme iteratively filter the signals. The denoised high-frequency components are then re-combined with the low-frequency signal at the last level, and an inverse wavelet transform used to recover the denoised stellar spectrum. We demonstrate that this method outperforms the conventional signal filtering methods when optimal method parameters are unavailable for the latter. This is indeed likely to be the case in practice due to an absence of a realistic ‘reference’ signal and noise to employ for obtaining optimal parameters for the conventional, parametric methods.

This paper is organized as follows.In §2 we discuss the wavelet transform and the wavelet shrinkage method. In §3 we briefly describe the Kalman filter. §4 describes our choice of wavelets and selection of decomposition levels to denoise, and demonstrates the effectiveness of our proposed method in addressing these issues. In §5, we introduce the synthetic PHOENIX spectra (Husser et al., 2013) used for comparing various filtering methods, and provide an in–depth explanation of our proposed algorithm. In §6 we use the PHOENIX data to compare the performance of our proposed method with commonly-used filtering methods from the literature. Finally, in §7 we comment on our findings and briefly discuss future extensions to the work.

2 Wavelet-based Signal Denoising

2.1 Discrete Wavelet Transform

The wavelet transform is similar to the Fourier transform (and more so to the windowed Fourier transform) in that both transform the input data into a different domain. However, the main difference between them is this—the Fourier transform decomposes the signal into sines and cosines, i.e. functions localized in Fourier space; in contrast the wavelet transform uses functions (wavelets) that are somewhat localized in both real and Fourier space. The sine and cosine functions used in the Fourier transform are ‘infinite’—in that they never go to zero in the time domain—and hence the signal is deconstructed into waves that are infinitely long. This introduces the ‘resolution problem’—if we have a high resolution in the frequency domain (i.e. focusing on one frequency in the signal) it is hard to isolate it in time, as each frequency exists across all time. Being uncertain of the time when focusing on the frequency is the flip-side of being uncertain of the frequency when focusing on time.

To overcome this resolution problem a wavelet transform is used to deconstruct the signal. The time-limited quality of wavelets is useful as it provides more resolution in the time domain. Instead of modeling with a wave that has infinite support, one uses a compactly-supported wave convolved with the input signal. To handle different frequencies, the wavelet transform employs varying scales for the wavelet. Wavelets in the WT are akin to sinusoidal waves in the FT: just as a sinusoid gets compressed (stretched) for high (low) frequencies, so do wavelets. For the WT, the signal is deconstructed using the same wavelet with different scales, rather than the same sinusoid at different frequencies. As there are hundreds of different wavelets, hundreds of different transforms are possible, each with an associated domain. However each domain has ‘scale’ on the horizontal axis rather than ‘frequency’, and provides better resolution in the time domain by being finite.

In this work, we use a specific type of wavelet transform—the Discrete Wavelet Transform (DWT) (Bultheel et al. 1995 gives an in-depth review of wavelet transforms and their properties). DWT is an implementation of the wavelet transform using a discrete set of wavelet scales and translations that obey several rules. DWT is particularly suitable for signal compression and denoising. There are several implementations of the DWT algorithm. The oldest and best-known one—and the one used in this work—is the Mallat (1989) (pyramidal) algorithm. With this algorithm, the data length is restricted to a power of $2$, requiring padding of the signal vector before applying the DWT. The noisy signal vector is simultaneously passed through a low-pass filter to obtain large scales (analogous to low frequencies in an FT), and a high-pass filter to obtain small scales (analogous to high frequencies in an FT). From the low-pass filter, we get a vector of approximation coefficients A0 that represents an estimation of the original signal with half resolution. From the high-pass filter, we obtain a vector of detail coefficients D0 that contains the details of the signal. The vector A0 can be further decomposed to form a new vector of approximation coefficients A1 and a new vector of detail coefficients D1.

With increasing decomposition level, less information will be included in the approximation coefficients. The lost information between approximation coefficients of two successive decompositions is encoded in the detail coefficients. This process can be iterated to level $\log_{2}N$, where the length of the original input is $N$. We must have $N$ a power of $2$ because the DWT is a recursive splitting. As a result, a vector of approximation coefficients and a series of vectors of detail coefficients are obtained that forms the DWT coefficients. The signal can be reconstructed by inverse DWT. Fig. 1 presents a block diagram demonstrating this concept. Fig. 2 demonstrates the DWT process when using a ‘Haar’ wavelet on a sawtooth signal. As is demonstrated there, at each successive level (scale), the mother and father wavelets are dilated by a factor of 2, and convolved with the Approximation coefficient from the previous level. Since the father wavelet is non–negative, this is equivalent to adding successive terms of the input signal (i.e., the Approximation coefficient from the previous level) in the time domain to get the Approximation coefficient for the next level. Similarly, the action of the mother wavelet is to find differences between pairs of points for all points in the time domain; this is how we obtain the Detail coefficient for the next level. These dual actions of scaling up (equivalent to increasing the window-size of the kernel in Fourier transform) and convolving with the mother and father wavelets help DWT achieve finite resolutions in both time and frequency domains—a significant advantage over the Fourier transform.

2.2 Shrinkage

DWT can be used for easy and fast denoising of a noisy signal. If we take only a limited number of highest approximation coefficients of the discrete wavelet transform spectrum and perform an inverse transform, we can obtain a partially denoised signal; this process is known as wavelet shrinkage. A more nuanced approached involves thresholding the values of the approximate coefficients that lie above a certain level (the threshold), instead of rejecting coefficients entirely (Johnstone & Silverman, 2005).

There are several different thresholding schemes in the literature, e.g. Johnstone & Silverman (2005), Li et al. (2007). The most common two are soft and hard thresholding. In soft thresholding, the $i^{th}$ component of the $j^{th}$ coefficient, $w_{j,i}$, will be set to zero if it is smaller than the threshold for that level $\lambda_{j}$. We note that $j$ runs from [math] up to the maximal decomposition level $k$, and we use just one threshold value at each level $j$. Otherwise, the $w_{j,i}$ are reduced by $\lambda_{j}$. With $\widetilde{w}_{j,i}$ denoting the thresholded coefficient, we have

[TABLE]

In hard thresholding, the coefficients that are smaller than the threshold value are set to zero; those equal to it or larger are not altered. We have

[TABLE]

The noise threshold $\lambda_{j}$ for a given Detail component can be obtained either using the universal level–independent threshold (Donoho & Johnstone, 1994) $\lambda=\sigma^{\rm Noise}\sqrt{2\log{N}}$, with $\sigma^{\rm Noise}$ an estimate of the noise and $N$ the length of the input signal, or a level-dependent threshold such as in Johnstone & Silverman (1997), which sets $\lambda_{j}=\sigma_{j}^{\rm Noise}\sqrt{2\log_{2}N_{j}}$ , where $N_{j}$ is the length of the $j^{th}$ Detail component, and $\sigma_{j}^{\rm Noise}$ is an estimate of the noise level.

Hard thresholding is better suited for the case where there is a Detail wavelet coefficient is either a signal or a noise coefficient. On the other hand, soft thresholding performs better when a Detail wavelet coefficient contains both signal and noise. While there are several thresholding functions in the literature that can be used (Chang et al., 2000; Zhao & Cui, 2015; Zhang & Bao, 2003), in this work, we have used the universal, hard, and soft thresholding schemes to compare to our proposed algorithm. Fig. 3 summarizes the steps involved in the denoising process, while Fig. 4 displays a block diagram of the same.

A detailed explanation and description of wavelet-based methods can be found in Graps (1995); Taswell (2000); Li et al. (2005); Tary et al. (2018).

3 Kalman Filter

The Kalman filter is a method for updating predictions in the presence of noisy measurements. Named for its inventor Rudolph E. Kalman (Kalman, 1960), it uses a predictor-corrector system to recursively improve the accuracy of predictions by accounting for the noise properties of both the underlying process and the measurements. It has been extensively applied in many areas including control systems, tracking, navigation, and hydrology (Andrews & Grewal, 2008; Yeh & Huang, 2005). Typically, the Kalman filter (and its advanced variants, such as the Unscented Kalman filter) are used to estimate the state of a time-varying process—such as tracking the location of an airplane using radar data and other relevant observables.

However, in this work we have used Kalman filtering to de-noise stellar absorption spectra by re-interpreting the time domain as the wavelength domain. Just as with time, points in wavelength follow a fixed ordering. However, unlike with time, wavelength does not have a fixed overall direction for correlations. Whereas by causality, a given timestep of a process should correlate only with previous timesteps, the value of a spectrum’s flux on a bin in wavelength is likely correlated with both lower-wavelength and higher-wavelength bins, as spectra are generally roughly predictable. Our problem is thus a more symmetric one than the usual cases where Kalman filtering is applied. To preserve this symmetry, we apply Kalman filtering in both backward and forward wavelength directions while denoising the spectra, as detailed in §5. For an intuitive derivation of the linear Kalman filter, we refer the readers to Appendix A.

3.1 Modeling the Process

Let us denote the measured, noisy signal at wavelength bin i by $\textbf{z}_{i}$, and the true, unknown signal by $\textbf{x}_{i}$. Let the process noise, i.e. the intrinsic noise describing the stochasticity of the system itself, independent of measurement, be $\textbf{w}_{i}$. Let the measurement noise—arising due to the in-practice imperfect process of measurement—be denoted by $\textbf{v}_{i}$.

The Kalman filter implementation used in this paper assumes $\textbf{w}_{i}$ and $\textbf{v}_{i}$ are normally-distributed random variables with independent respective covariance matrices $\textbf{Q}_{i}$ and $\textbf{R}_{i}$. Since the Kalman filter is a recursive filter, the true (clean) signal at wavelength i depends on the signal at wavelength $i-1$; let the evolution operator from $i-1$ to $i$ be denoted by $\textbf{M}_{i}$. Finally, let the operator relating the measured signal $\textbf{z}_{i}$ to the true signal $\textbf{x}_{i}$ be denoted by $\textbf{H}_{i}$. Then for any wavelength point i the true signal (‘state’, hereafter) and the measurement signal satisfy the process equation (3) and the measurement equation (4), respectively:

[TABLE]

Now that we have introduced the main elements of the Kalman filter, we proceed to briefly describe the set of recursive operations that estimate the clean signal at each wavelength i. This is achieved in two steps (Welch & Bishop, 1995): ‘prediction’ and ‘update’, as detailed in §3.2 and §3.3. We outline the Kalman filtering process assuming we are moving in the direction of increasing wavelength; however as discussed earlier (§3), we also apply our process backwards. The appropriate equations can simply be obtained by switching $i-1$ and $i$.

3.2 Predicting the Process

In the ‘predict’ phase, the Kalman filter uses information from the previous state to estimate the a priori current state. This estimate is then updated using information from a measurement to produce an a posteriori state estimate: this is termed the ‘update’ phase.

In more detail, first we initialize a mean ($\textbf{\^{x}}_{-1}$) and a covariance matrix ($\textbf{P}_{-1}$) for the state of the system, which are advanced to obtain an a priori estimate of the state (mean $\textbf{\^{x}}_{0}$ and covariance $\textbf{P}_{0}$). The covariance part is necessary to handle the uncertainty in the estimate arising due to the failure of the linear model used for the system. This is the ‘prediction’ part of the Kalman filtering process—where we use the estimate of the system state from the previous run of the process to predict the a priori for the next run; i.e., use the flux predicted for wavelength i-1 to get an a priori estimate for the flux at wavelength i. This is encoded in equations (5) and (6):

[TABLE]

3.3 Updating the Process

Next, we use the measured signal at wavelength i ($\textbf{z}_{0}$, since $i=0$ in the first iteration) to update the a priori state prediction from the previous step (i.e., $\textbf{\^{x}}_{0}^{-}$ and $\textbf{P}_{0}^{-}$). This measurement comes from an instrument (a spectrograph in our case); its uncertainty is the measurement noise ($\textbf{v}_{i}$ at wavelength i, with assumed mean [math] and covariance matrix $\textbf{R}_{i}$). This is induced by inaccuracies in the measurement process (e.g. an imperfect CCD chip or an unstable spectrograph). In addition, we class within this measurement error any errors arising due to imperfect modeling of the system from our assumption that its evolution is linear; we highlight in §5 our rationale for this choice.This phase is known as the ‘correction’ or the ‘update’ phase; it is encoded in equations (7), (8) and (9):

[TABLE]

In equation (7), $\textbf{K}_{i}$ refers to the Kalman gain, which encodes the filter’s confidence in the prediction $\textbf{\^{x}}_{i}^{-}$ versus in the measurement $\textbf{z}_{i}$. As one might expect, the Kalman gain depends on the error covariance matrices associated with both these quantities. In the limiting case where the measurement error at a wavelength point i is much smaller than the prediction error (i.e., $\textbf{R}_{i}\ll\textbf{P}_{i}^{-}$), $\textbf{K}_{i}\to\textbf{H}_{i}^{-1}$ and the a posteriori state matches the measurement as acted upon by the inverse of the measurement matrix (i.e, $\textbf{\^{x}}_{i}\to\textbf{H}_{i}^{-1}\textbf{z}_{i}$ and $\textbf{P}_{i}\to 0$). Since the present work deals with denoising spectra as measured, which are obtained by spectrographs, both the true state x̂ and the measurement z refer to spectra. Hence all $\textbf{H}s$ have been assigned the identity matrix, and the a posteriori state matches the measurement exactly ($\textbf{\^{x}}_{i}\to\textbf{z}_{i}$). This makes intuitive sense—when we are extremely confident about a measurement, we can assume that to be the real state of the system. On the other hand, when the measurement error dominates (i.e., $\textbf{R}_{i}\gg\textbf{P}_{i}^{-}$), $\textbf{K}_{i}\to 0$ and the a posteriori state matches the a priori state (i.e, $\textbf{\^{x}}_{i}\to\textbf{\^{x}}_{i}^{-}$ and $\textbf{P}_{i}\to\textbf{P}_{i}^{-}$). When we lack confidence in a measurement, it is prudent not to update the original a priori state estimate. In equation (9), I refers to the identity matrix.

This entire process is repeated recursively until we have traversed the length of the input signal. As noted earlier, to encode the symmetry under traversing the wavelength range forwards or backwards, we run our spectral denoising algorithm for each traversal direction.

4 Limitations of Current Methods

The standard soft– and hard–thresholding functions, along with other commonly used threshold-selection schemes, suffer from a few major drawbacks that preclude their use for automated signal denoising:

The choice of $k$, the number of wavelet transform decomposition levels to threshold, is data-dependent and needs to be carefully chosen by the user. 2. 2.

The choice of $\sigma_{j}^{\rm Noise}$, the noise estimate at the $j^{\rm th}$ decomposition level, greatly influences the noise threshold $\lambda_{j}$, but there is no definitive way to estimate it. Different noise estimates yield different noise thresholds (Donoho & Johnstone, 1994, 1995; Johnstone & Silverman, 1997; Srivastava et al., 2016). 3. 3.

Assuming white Gaussian noise (WGN), a single noise threshold is selected and applied to the magnitudes of both negative and positive Detail coefficients. However for WGN, the assumption of no distribution bias (i.e. a distribution symmetric about zero) may often be flawed (Srivastava et al., 2016). 4. 4.

The methods used result in fixed thresholds that are not adjusted. Flexibility in adjusting thresholds is extremely important to optimally threshold for signals, especially when the wavelet coefficients of weak signals are close to the maximum magnitude of noise (Srivastava et al., 2016).

The method proposed in the current work is considerably more ‘non-parametric’ than both the traditional wavelet denoising schemes discussed above, and traditional smoothing or filtering schemes like the Fourier transform, moving mean, etc. The only ‘parameter’ of consequence that our proposed algorithm requires from the user is a suitable wavelet function. The criteria for choosing this basis function are: maximization of its correlation with the signal, and minimization of its correlation with the noise. Clearly, we do not want to choose an arbitrary wavelet function. The selection and development of such data-specific basis wavelets is well-studied (e.g. Rafiee et al. 2009). In the current work, we use the simplest and computationally cheapest wavelet (‘Haar’).

5 Automatic Stellar Denoising

In §5.1, we describe the nature and origin of the stellar spectra used to compare the performance of the proposed method with the major denoising methods from the literature. In §5.2, we explain in detail the implementation of two versions of our proposed algorithm. We then introduce error metrics used to compare the various filtering algorithms, and end by demonstrating the importance of the non-parametric nature of our proposed algorithm.

5.1 Test data

We use synthetic 1D PHOENIX spectra (Husser et al., 2013) for three stars of different spectral types. Star 1 has T${}_{\rm eff}=4,000$K, [Fe/H] = $1.0$ dex, $\log g=2.0$ dex, and spectral type = KIII. Star 2 has T${}_{\rm eff}=6,000$K, [Fe/H] = $1.0$ dex, $\log g=3.5$ dex, and spectral type = GV. Star 3 has T${}_{\rm eff}=7,000$K, [Fe/H] = $1.0$ dex, $\log g=2.5$ dex, and spectral type = KVI. The spectra are convolved from the native PHOENIX resolution of $R=\lambda/\Delta\lambda=500,000$ down to the three lower resolutions of $2,000$, $5,000$ and $10,000$, which are then re-sampled onto the optical wavelength interval of 5,000-5,495.9$$\mathrm{\SIUnitSymbolAngstrom}, with a sampling of 0.1$$\mathrm{\SIUnitSymbolAngstrom}. The wavelength range and the sampling chosen are for illustration purposes and have no bearing on the generalizability of our proposed algorithm. Our sampling yields vectors of length $4,096=2^{12}$. Each of the nine spectra were disturbed by twenty-five different levels of Poisson noise, to obtain final PSNRs (peak signal–to–noise ratios) going in steps of unity from $1.0-25.0$. For each spectrum and PSNR, we produced $100$ realizations; this point is further discussed in §6. Both the reference spectra and the corresponding noisy spectra were normalized by the peak flux in the reference spectra to make them all lie on a comparable flux scale. In total we deal with $9\times 25\times 100=22,500$ unique vectors to denoise, each consisting of $4,096$ data points.

5.2 Adaptive Kalman-Filter-based Thresholding (AKFThresh)

As described in Section 4, denoising performance with conventional wavelet shrinkage methods is a strong function of the number of decomposition levels and the thresholding scheme. Our proposed algorithm deals with this dependence by treating the individual decomposed Detail components of the noisy input signal as separate signals to denoise, and utilizes distinct Kalman filters (see equations 5-9) for each.

In brief, our proposed thresholding scheme works as follows. With every Detail level $j$, from $0\leq j\leq k$, we associate a Kalman filter. At this level, the noisy signal $\textbf{Z}_{Dj}$ is the Detail component, the denoised signal is $\textbf{X}_{Dj}$, and $\textbf{H}_{Dj}$ is set to I (the identity matrix). We are thus tasked with obtaining $\textbf{Q}_{Dj}$, $\textbf{R}_{Dj}$, and $\textbf{M}_{Dj}$ (see equations 5-9).

Since the $0^{\rm th}$ level Detail component consists of the highest-frequency part of the spectral signal, we assume it to be completely noise; specifically, we set the process noise covariance vector for the $0^{\rm th}$ Detail component $\textbf{Q}_{D0}$ to 0, and use the algorithm suggested by Meng et al. (2016) to obtain the measurement noise covariance matrix $\textbf{R}_{D0}$. Fig. 6 illustrates the importance of choosing $\textbf{R}_{D0}$ in this manner. Starting from the $1^{\rm st}$ level Detail component, all the measurement noise covariance matrices $\textbf{R}_{Dj}$ are set equal to $\textbf{R}_{D0}$, while the process noise covariance matrices $\textbf{Q}_{Dj}$ are obtained by subtracting the variance of the $0^{\rm th}$ Detail component from the covariance of the $0^{\rm th}$ and the $j^{th}$ Detail components. Even though we are dealing with synthetically-generated spectra in the absence of a measurement device, the measurement noise matrices $\textbf{R}_{Dj}$ are not 0 since they capture the errors emanating from the imperfect assumption of linearity in our Kalman filters.

At the same time, we estimate the state transition matrix ${\textbf{M}_{Dj}}$ by following the Support Vector Regression (SVR) procedure outlined in Salti & Di Stefano (2013) rather than setting it to unity. It is precisely this heuristic of obtaining $\textbf{R}_{Dj}$, $\textbf{Q}_{Dj}$, and $\textbf{M}_{Dj}$ for the Kalman filters at the various Detail levels that allows us to make the process effectively non-parametric, instead of finding the optimal number of decompositions as is done traditionally. The pseudo-code for the proposed method is given in Fig. 7. All the variances and covariances are obtained in a robust fashion.111https://scikit-learn.org/stable/modules/generated/sklearn.covariance.MinCovDet.html

We use two versions of our denoising algorithm—AKFThresh v1 refers to the algorithm outlined above (and in Fig. 7), whereas AKFThresh v2 uses the denoised output vector from AKFThresh v1 as input, takes the natural logarithm of the values (after adding an offset to make them strictly positive, if required), applies the AKFThresh v1 algorithm again, and exponentiates the resulting vector (followed by removing any offset that was added). The Kalman filter assumes the process and measurement noises to be additive (see equations 5 and 6). Since Poisson noise—unlike Gaussian noise—is not strictly additive, this additional process is done to convert any residual multiplicative noise to additive form, and then remove it (Zhang et al., 2014).

We compare the denoising performance of the above two versions of the proposed algorithm with those of popular denoising algorithms in the field. These are: the moving mean,222https://pandas.pydata.org/pandas-docs/stable/generated/pandas.Series.rolling.html the moving median,††footnotemark: exponential smoothing,333https://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.ewm.html Fast Fourier Transform with a Gaussian kernel,444https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.gaussian_filter.html Wiener filter,555https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.wiener.html the Savitzky–Golay (SG) filter,666https://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.signal.savgol_filter.html wavelet transform followed by soft thresholding, and wavelet transform followed by hard thresholding. Since these commonly used methods are all parametric, any such comparison must be preceded by the selection of an appropriate loss function to minimize, followed by a search in their respective parameter spaces to optimize this loss function. Following Köhler & Lorenz (2005) we calculate three loss functions—mean absolute error, or the $L^{1}$ norm, least squares error, or the $L^{2}$ norm, and maximum absolute error, or the $L^{\infty}$ norm. We also consider a fourth function, the structural similarity index measure or SSIM, which is a measure of the similarity between two images (Wang et al., 2004; Wang & Bovik, 2009).

Each error metric probes the differences between given signals in a different way. We illustrate this with an example in Fig. 5, by highlighting the differences between the mean squared error (MSE)—a ‘global’ error metric—and the structural similarity index measure (SSIM)—a ‘local’ error metric. We assume a sinusoidal form for the ‘clean’ signal, and create two corrupted versions of it. The MSEs for the two versions of the ‘noisy’ signal (dotted blue and dashed red curves) with respect to the real signal (solid black curve) are the same, while the SSIM is significantly higher for the dotted blue curve. This is also apparent by visual inspection—the dotted blue curve, while having a larger amplitude than the original signal, is clearly more ‘similar’ to it, in terms of capturing the peak and the valley. In contrast, the dashed red curve, being identically zero, is devoid of any information about the original signal.

We should also optimize the parameters within each method with respect to the measure of quality we want to use. For instance, in a use case, if one wished to achieve minimum least-squares error on a given data set, clearly one would optimize one’s parameters (if possible) under that error metric. We analyzed each norm to fully assess the strengths and weaknesses of the different denoising methods under each.

We define appropriately expansive parameter spaces for the different algorithms and do a thorough search to find the optimal parameters to minimize the $L^{2}$ norm between the denoised signal and the original reference signal. We repeat this analysis two more times, with the difference that now we look for parameters that give the tenth-lowest and the twentieth-lowest $L^{2}$ norms. This is done to simulate the situation where, in the absence of a reference signal, manual parameter search by trial-and-error yields sub–optimal parameters for a given denoising algorithm. We then compare the denoised signals from both versions of our proposed algorithm (AKFThresh v1 and AKFThresh v2) to the denoised signals from the optimal and sub-optimal versions of the other algorithms; partial tables displaying these results are available in Appendix B (complete tables are available online only777https://github.com/astrogilda/stellar_denoising). For robustness, this entire procedure was repeated for $100$ different realizations of the noise to create synthetic noisy signals for each of the three stars, and the results averaged. A high-level pictographic overview of our proposed method is presented in Fig. 8, where we show its denoising performance in action for the absorption spectrum for Star 1 corrupted by Poisson noise to bring the PSNR down to $10$.

Our code is written entirely in python; to implement Kalman filters and wavelet decompositions, we have used the libraries filterpy (Labbe, 2014) and pywavelets (Lee et al., 2006), respectively. The complete codebase and datasets used may be found online.888https://github.com/astrogilda/stellar_denoising

6 Results

We discuss here the relative denoising performances of the methods introduced in §5.2. In the interest of saving space, we only present results for synthetic spectra of Star 2 (T${}_{\rm eff}=4,000$K, [Fe/H] = $1.0$ dex, $\log g=2.0$ dex, spectral type = KIII, further details in §5.1); the results for the other two stars are similar.

Tables 2 through 10999https://github.com/astrogilda/stellar_denoising present the values of the different comparative metrics, at three different noise levels, for the optimal and sub-optimal parameter cases, at a resolution of $5,000$ (again, the results are qualitatively similar at the other two resolutions of $2,000$ and $10,000$). Each table shows the mean value and the standard deviation of the results for the $100$ different runs of the respective method. The best and the worst values for each column are highlighted—the best in boldface and the worst in italics. In addition, the median parameter values used for the parametric denoising methods have been specified.

In Figs. 10, 10, 12 and 12, we present the values of the four performance metrics for the best four denoising algorithms—FFT with Gaussian kernel, Savitzky-Golay, Wavelet + Soft Thresholding, and Wavelet + AKFThresh v2—as a function of input signal PSNR and quality of parameters for the different denoising methods. The highest quality parameter is that which gives the lowest $L^{2}$ for a given method, the second-highest quality parameter that which gives the second lowest $L^{2}$, and so on. As the standard deviations are relatively small, they have been omitted for clarity.

In Fig. 13, Fig. 16, and Fig. 17 we compare the denoising performance of the best version of the proposed algorithm (v2) with that of Fourier transform with optimal and sub-optimal parameters. We have only shown a small section of the input spectrum to better emphasize the differences between the different spectra. At PSNRs $\geq 5.0$, our algorithm is able to extract the peaks and valleys of the original input spectrum from the noisy spectrum much better than the Fourier transform. Its superiority in filtering out the Poisson noise is especially highlighted when the width of the Gaussian kernel for the FFT is not optimal for the specific noisy spectrum. This is indeed likely to be the case when the user is responsible for picking it without access to the original spectrum.

Finally, in Figs. 14 and 15, we show the performance of both versions of our denoising methods as a function of the decomposition level. Performance as measured by all four error metrics more or less saturates once decomposition exceeds the $6^{\rm th}$ layer. This saturation justifies our decision, and recommendation, to decompose the input noisy signal to half the maximum possible number of decomposition levels. In this work, the maximum number of levels possible is $12$ because we had spectra of $4096=2^{12}$ wavelength bins).

Based on these results, we can draw several conclusions:

At almost all noise levels and all choices of parameters, exponential smoothing has the worst performance. 2. 2.

DWT with both soft and hard thresholding gives excellent results across all four metrics for the optimal choice of parameters (i.e. maximal decomposition level, thresholding), but the performance quickly deteriorates for an incorrect choice of parameters. We emphasize that one does not know the correct parameters a priori, so this issue is a significant one. 3. 3.

At extremely low PSNRs ($<5$), both versions of our algorithm provide poor denoising results, even when compared to the performance of other methods with sub-optimal parameters. 4. 4.

At PSNRs $>5$, both versions of our algorithm start outperforming other methods, especially when sub-optimal parameters have been chosen for the other methods. This scenario is the more likely one in practice since one does not always have a reference spectrum with which to work. Indeed, calibrating a parametric denoising method on synthetic spectra before applying to real data comes with its own set of issues, primarily interpolation errors when creating the most closely-matched synthetic spectrum, as well as errors originating from instrument peculiarities’ not being captured in the synthetic spectra. 5. 5.

Savitzky–Golay is the only algorithm that compares favorably to the algorithm proposed in this paper with respect to the four error metrics we explore (Figs. 10–12). However, we note that during parameter selection for the commonly used parametric methods including Savitzky-Golay, we explored a rather wide range of parameters (e.g. for the choice of the window-length ‘w’, see footnotes for Tables 2 through 10). This was computationally practical for a small number of spectral types (we considered just three different types of star in this work). However for a large survey dataset with millions to tens of millions of objects (e.g. DESI will have $\sim 10^{7}$ stellar spectra), with many different possible spectral types, such an exploration would likely be computationally demanding. In particular, one might, for each expected spectral type, optimize the parameters using synthetic spectra of that type. But then one would still need to build in a step where, for a given data spectrum, the correct type was automatically identified to allow the optimal parameters for that particular type for a given algorithm to be used. Having a method that simply does not require optimal parameter selection at all seems more efficient.

7 Conclusions and Future Directions

In this work, we have proposed a new method for denoising 1D stellar spectra. Our proposed method involves two stages—a discrete wavelet transform, followed by adaptive Kalman-filter-based thresholding. By using the ‘Haar’ wavelet basis, the only free parameter for the user is the number of levels of decomposition for the wavelet transform. By setting this to half the maximum possible value (where the improvement of adding more levels essentially saturates), the number of free parameters can in fact be reduced to zero. Our proposed method has been validated on several types of synthetic PHOENIX spectra with Poisson noise injected. Our tests show that the method can improve denoising better than conventionally employed tools (up to $3\times$ if the conventional methods have even slightly sub-optimal parameters; for instance as shown in Figs. 10–12). The method performs well all the way down to a PSNR of roughly 5.

Our method may also be useful denoising galaxy and quasar spectra, which often have low SNRs ($<20$). There are many large current and upcoming such datasets that offer fruitful use-cases: the Sloan Digital Sky Survey (SDSS), Gaia, the Large Synoptic Survey Telescope (LSST), the Public ESO Spectroscopic Survey for Transient Objects (PESSTO), the ‘deep VIMOS survey of the CANDELS UDS and CDFS fields’ (VANDELS), the Large Early Galaxy Astrophysics Census survey (LEGA-C), and the Dark Energy Spectroscopic Instrument survey (DESI). DESI in particular will have low SNR, low-resolution spectra of roughly 30 million galaxies and quasars, and denoising these might aid both in redshift fitting for large-scale structure clustering analyses as well as in ancillary science. Denoising spectra of the Lyman-$\alpha$ forest as measured using quasars as “skewers” might also enhance the science possible e.g. regarding warm dark matter constraints (for instance, Viel et al. 2013, Iršič et al. 2017) or intergalactic medium physics (e.g. Gurvich et al. 2017).

We note that although our proposed algorithm has been tested and validated on stellar absorption spectra, we believe it would perform just as well on emission spectra, which can be thought of as ‘flipped’ absorption spectra. In future work, we plan to evaluate the performance of AKFThresh v2 on spectra from the public data releases of the surveys listed above.

The current work is just a first step in developing better denoising algorithms based on sophisticated statistical methods, drawing on advances in the fields of digital signal processing (Oktar et al., 2016), autonomous vehicle control (Kasper & Schmidt, 2008), and financial timeseries analysis (Renaud et al., 2005), among others. There is consequently significant scope for improvement and development. For instance, we are currently updating the linear Kalman filter employed here to the Unscented Kalman filter capable of modelling non-linear state-space systems (Wan & Van Der Merwe, 2000; Deng et al., 2013). We are also exploring the impact of developing and using a custom mother wavelet instead of the ‘Haar’ wavelet. Another possible line of research would be to apply Kalman filtering to reconstructed decomposition levels, estimating the noise and process covariances via an ARMA process (Huang, 2015). One might also investigate deviations from the assumption of Gaussianity for the process or measurement noise, or of the assumption of zero covariance between them (Liu et al., 2016). These ideas might enable denoising down to lower PSNRs, $<5$.

In addition to wavelet shrinkage as demonstrated here, one might also explore the application of Kalman filtering in conjunction with other DWT-based methods like wavelet packet thresholding (Oktar et al., 2016; Zhao & Zhu, 2017), wavelet coefficient modeling (Chen & Bui, 2003; Jiang et al., 2010; Kumar & Saini, 2012), and the wavelet transform modulus maxima (WTMM) method (Liu & Cheng, 2007). Finally, we see no reason that this method cannot be extended to denoising astronomical time series, though this would require analysis of the data’s properties to develop an appropriate basis wavelet.

Acknowledgements

We thank Stephen Eikenberry and Paul Torrey for productive discussions, and Peter Barnes and Kai Polsterer for a careful review of the original draft. We particularly acknowledge Jian Ge for his kind guidance and mentoring over the past several years, and many helpful insights on making this paper more applicable to current science needs.

Appendix A Derivation of the Linear Kalman Filter

Here we present an intuitive derivation of the linear Kalman filter used in this work. We have presented the derivation in terms of filtering a temporal signal, as this is the traditional application of a Kalman filter. However we emphasize that in this work the Kalman filter has been used to denoise spectral data in wavelength domain, leading to applying it forwards and backwards in wavelength (as discussed in §5.2). We have also used different, easier–to–understand symbols than in the main text above; the latter have been used to match the practice in signal processing literature. The transformations between the symbols used in the main text above and those used in the derivation here are specified in Table 1.

Kalman filtering answers the question of how one should optimally update predictions in the presence of noisy measurements. The prediction is smooth but likely, eventually (in a time series) to become biased; the measurements are, one hopes, unbiased (noise averages to zero), but noisy. How does one employ both?

Let us denote the true signal (also called the state vector) at timestep $i$ by ${\bf t}_{i}$ and the measured signal by ${\bf m}_{i}$. These are vectors, denoted by boldface. We denote the intrinsic noise, i.e. describing stochasticity of the system itself, independent of measurement, as ${\bf n}_{i}^{t},$ superscript $t$ for “true.” We denote the measurement noise (associated with the process of measurement) ${\bf n}_{i}^{m},$ superscript $m$ for “measurement.” We also define an evolution matrix ${\bf\hat{E}}_{i}$ that takes the true signal from time $i$ forward to time $i+1$. The hat represents that this matrix is an operator, and operates on a state.

We also define a measurement matrix $\hat{{\bf M}}_{i}$, where the hat denotes that it is an operator (and is useful in distinguishing it later from another matrix ${\bf M}$ for the covariance matrix of the measurement). We also note that the hat implies the operator may not commute. However, that is not problematic—in general matrix multiplication does not commute, so we are not implying an additional property beyond what should already hold.

Finally, we make the following assumptions (the derivation and operation of the Kalman filter when these assumptions do not hold is beyond the scope of this paper; interested readers are referred to the references in §5.2):

• •

${\bf n}_{i}^{m}$ is 0-mean white noise: $E\left({\bf n}_{i}^{m}\right)=0$, ${\rm cov}\left({\bf n}_{i}^{m},{\bf n}_{j}^{m}\right)={\bf N}_{i}^{m}\delta_{ij}^{K}$

• •

${\bf n}_{i}^{t}$ is 0-mean white noise: $E\left({\bf n}_{i}^{t}\right)=0$, ${\rm cov}\left({\bf n}_{i}^{t},{\bf n}_{j}^{t}\right)={\bf N}_{i}^{t}\delta_{ij}^{K}$

• •

intrinsic noise and measurement noise are uncorrelated: ${\rm cov}\left({\bf n}_{i}^{t},{\bf n}_{i}^{m}\right)=0$

• •

true signal at timestep $i=0$ (denoted by ${\bf t}_{0}$) is uncorrelated with the other two noise vectors. $E\left({\bf t}_{0}\right)=\overline{t_{0}}$, ${\rm cov}\left({\bf t}_{0},{\bf t}_{0}\right)={\bf T}_{0}$.

Above, $E$ means the expectation value (also called the first moment), $\rm{cov}$ is covariance matrix, and $\delta_{ij}^{K}$ is the Kronecker delta ($\delta_{ij}^{K}=1$ if $i=j$, [math] otherwise).

With this notation in hand we may now describe the system.

Description of system

[TABLE]

The true signal at time (really, wavelength bin) $i$ is given by the true signal at bin $i-1$ evolved by the evolution matrix at bin $i$ and with the intrinsic (“process”) noise ${\bf n}_{i}^{t}$ at that time added in. Equation 10 is known as the state transition equation.

Similarly, the measurement ${\bf m}_{i}$ at bin $i$ is given by the measurement equation

[TABLE]

We take the true signal at bin $i,$ apply the measurement matrix $\hat{{\bf M}}_{i}$, and then add in the measurement noise at that time ${\bf n}_{i}^{m}$.

Note that both of the above equations are linear functions of the state (i.e., true signal), and have additive noise. This is because in the current work we have used the linear Kalman filter; for descriptions of non–linear filters such as the Extended Kalman Filter, we offer references in §5.2. We now may also describe our predictions.

Prediction step

Our predicted signal at bin $i$ is

[TABLE]

To predict, we simply evolve our prediction from the previous bin using the evolution matrix. We also would like the covariance matrix of the predicted true signal, which former we denote ${\bf T}_{i}^{{\rm pred}}.$ We have

[TABLE]

The right–hand side says that we take our updated covariance matrix from the $i-1^{th}$ bin, evolve it (we need two evolution matrices since one can think of the true signal matrix as the outer product of two true signal vectors, each of which must be evolved by $\hat{{\bf E}}$), and then add the covariance matrix of the true signal noise, i.e. ${\bf N}^{t}$ is the covariance matrix of ${\bf n}^{t}.$ We have now completed our prediction step. But where does the updated signal covariance matrix ${\bf T}_{i-1}^{{\rm upd}}$ come from? How did we update this prediction from a previous one, and how will we update this prediction to move to the next one?

Update step

The update step is the non-trivial part of Kalman filtering: it balances the prediction with whatever measurement we have made. The updated signal is

[TABLE]

Our updated true signal (${\bf t}_{i}^{{\rm upd}}$) is the predicted one (${\bf t}_{i}^{{\rm upd}}$) plus a weight (${\bf W}_{i}$) applied to the difference between our prediction for the measured signal and the actual signal we measured. We denote this difference as

[TABLE]

The superscript $m-p$ means “measured minus predicted.” The weight ${\bf W}_{i}$ is commonly referred to as the “Kalman Gain”.

In other words, our actual signal measured at bin $i$ was ${\bf m}_{i},$ but our prediction for the signal we would measure was our measurement matrix operating on our prediction ${\bf t}_{i}^{{\rm pred}}$ for the true signal. Equation 14 makes intuitive sense: to update, we begin with our predicted signal, and then correct it by adding in the weighted difference between our actually measured and predicted-to-be-measured signals (the weight can also take care of converting a measured signal to a true signal).

The correct form of ${\bf W}_{i}$ is the only non–obvious piece of Kalman filtering. Clearly, if we had no measurement noise and if the covariance matrix for our predicted signal were the identity matrix (i.e. all bits of the signal are independent from each other and have equal error bars), then all we would need to do would be to convert ${\bf\Delta}$ from the “measured” space to the “true signal” space. In this limit, all we would need to do would be to invert our measurement matrix. We do not need anything additional to inverse–covariance–weight since we have said that our intrinsic true signal has covariance matrix as an identity matrix; inverse–covariance weighting would thus also be an identity matrix. So we have

[TABLE]

$``{\rm i.e.e.\;sign}.$” stands for “independent, equal error bars signal.”

However, in practice we have measurement noise, and we also want to inverse–covariance weight the update, which is less trivial if the true signal has intrinsic covariance matrix different from the identity matrix. In this case, the predicted $measurement$ covariance matrix is

[TABLE]

we just convert the predicted covariance matrix of the true signal into one for the measurement by applying two copies of the measurement matrix. We use two copies for the same reason as discussed earlier, regarding evolving the predicted true signal covariance matrix.

We note that the matrix ${\bf M}_{i}^{{\rm pred}}$ describes the intrinsic covariance of a noiseless measurement—it occurs because our underlying true signal has some intrinsic covariance, and we are measuring that signal. However, the full covariance of the measurement will also include a contribution from the measurement noise covariance matrix, denoted ${\bf N}_{i}^{m}$. So our total measurement covariance matrix is

[TABLE]

Thus our inverse covariance weight so far looks like

[TABLE]

where to obtain the second equality we simply substituted equation 16 for ${\bf M}_{i}^{{\rm pred}}$. However, in the limit of no measurement noise, this weight should simply reduce to the inverse of the measurement matrix. So we have

[TABLE]

The first equality comes from demanding that the weight, in the limit of no measurement noise, should just be the inverse of the measurement matrix. We then take ${\bf A}$ to be an unknown matrix that multiplies our guess for the weight in the noiseless limit (with ${\bf N}_{i}^{m}\to 0$).

In the second equation above, ${\bf A}$ is an unknown matrix for which we wish to solve. We have a matrix equation of the form

[TABLE]

so we apply ${\bf B}^{-1}$ from the right on each side. We have

[TABLE]

meaning

[TABLE]

This simplifies to

[TABLE]

So our full weight is now

[TABLE]

We now insert the weight from equation 24 into equation 14 to update our prediction.

Appendix B Tables

Incomplete tables are displayed here. Complete tables are available online.101010https://github.com/astrogilda/stellar_denoising.

Appendix C Denoising Comparison

We present the comparative denoising results of our proposed algorithms with state–of–the–art algorithms, for peak signal–to–noise ratios of 10.0 and 20.0. See also Fig. 13 for comparison.