Optimal Filtering for Deep Sky Lucky Imaging

November 11, 2019

Optimal Filtering for Deep Sky Lucky Imaging

by

Robert E. Majewski Ph.D.

In order to increase the signal to noise ratio of the each short exposure, we use optimal filtering. The optimal or "Wiener" filter was developed by Nobert Wiener in the 1940s.  The shape of this filter adapts to the signal and noise power spectra in order to minimize the least-squares error between the estimated underlying signal and the true signal.  Let |C(f_x,f_y)|² be the signal plus noise power spectral density PSD, |N(f_x,f_y)|² is the noise PSD and f_x and f_y are spatial frequencies. The optimal filter, Φ(f_x,f_y) is given by¹

One of the biggest noise sources for deep sky imaging is shot noise due to the sky radiance.  This is especially a problem at observing sites that are light polluted.  The sky acts like a flood source that is raining photons down on the focal plane in a random fashion.  Since the photons detected from the sky are uncorrelated from one pixel to another, this noise component is spectrally white.  Thus, the average value of this noise power NP is constant, independent of spatial frequency.

NP can be determined from the average level of the sky background Bk.  This is a consequence of Parseval's theorem.  If the FFT region is N x N, the noise spectral power value NP is given by (this depends on the FFT normalization¹ used)

where g_e is the number of photo-electrons per count.

This functional behavior makes sense for this application.  At low spatial frequencies, where the PSD is high compared with this sky shot noise, the filter has a value of one.  At high spatial frequencies, where the PSD is dominated by the sky shot noise, the filter value drops down to zero.  It transitions from one regime to the other in an optimal manner.

Moffat Functions

Moffat functions are widely used to fit stellar profiles, because they can fit the wings of the point spread functions better than can a Gaussian function.  The functional form of the Moffat function is given by

where r is distance from the center and the stellar Full Width at Half Maximum FWHM is given by

When β is equal to one, the profile is called a Lorentzian which is unphysical for this application, because the amount of energy contained is unbounded.  As β increases toward infinity, the Moffat function approaches a Gaussian.  So β controls how sharply the wings tail off.  The Kolmogorov theory of atmospheric turbulence results in a seeing profile that can be best approximated by a Moffat function with value for β = 4.765².  However, in practice β values in the range of 2 to 4 better represent actual data.  For the construction of the optimal filter for deep sky lucky imaging we will use a somewhat more conservative value for β of 1.5.  This value generates a profile with finite total energy and rolls off in the manner expected for real optical systems.  Also, its Hankel transform is a simple exponential function and thus easy to compute. Examination of PSD from actual astro-imagery confirms that this closely matches the data.  This function rolls off in the spatial frequency domain more slowly than the Kolmogorov case and is less likely to generate image artifacts.

Exponential PSD

Since the PSD is mainly driven by atmosphere seeing, we will use star profiles and the FWHM values measured to create the PSD which has the form

The normalization constant A is determined from an average of the PSD computed directly from the data in an annular region near the origin.  The actual PSD does deviate from this form at low spatial frequencies, but this is of no consequence since the value is large and the Wiener filter will take on a value of ~ unity anyway. The reason this model is chosen is that it is important that the optimal filter be smooth in order to avoid creating ringing or other types of image artifacts.  In order to use the data directly would require a great deal of spectral averaging.

This filtering process has little effect of the image sharpness but is effective at reducing high spatial frequency noise in the image.  FWHM values obtained before filtering remain nearly the same afterward.  Figure 1 is an example of what an optimal filter looks like using the exponential PSD.  The smaller the FWHM value is the wider the bandwidth of the optimal filter.

Figure 1  An example of optimal filters for an FWHM value of 2 arcseconds (black), 3 arcseconds (red) and 4 arcseconds (blue) with  0.47 arcsecond pixels.

An exponential roll off in the spatial frequency domain is the 2D Fourier transform (Hankel transform due to the rotational symmetry) pair of a point spread function with the functional form of a Moffat function with beta = 1.5 as

The Hankel Transform of this is proportional to

and the power spectrum will be proportional to the square of this or

Results M57

The Ring Nebula in Lyra was used as a test case.  Figure 2 shows a PSD when the optimal filter was not used.

Figure 2  The PSD without an optimal filter.

Figure 3  The PSD with an optimal filter in red on top of the spectrum from figure 2.

Figure 3 shows the effect of the optimal filter in red on top of the spectrum shown in figure 2.  In figure 2 we see the large flat plateau at high spatial frequencies due to the sky shot noise.  At lower spatial frequencies, we see a downward sloping region that is due to the stars and nebula in the image.  The fact that the falloff is straight in this linear-log plot shows that the PSD is well approximated by an exponential function.

Figure 4 10 minute sub obtained without optimal filtering.

Figure 5. 10 minute sub obtained with optimal filtering.

The 2 images are 10 minute subs which are the result of integrating a large number of 5 second expoures.  The optimal filtering in figure 5 was done on each individual 5 second exposure during the image capture process.  The star FWHM values and background levels were used to create optimal filters which were unique for each 5 second exposure.  The difference is somewhat subtle, but the image in figure 4 appears to be dirty compared figure 5 where the optimal filter is applied.  One may need to enlarge the images to really see this since it is the fine structure of the noise that is being addressed here.  This high spatial frequency noise makes the image look grainy.  Optimal filtering makes the image look smoother. Reducing this high frequeny noise is important if we want to use deconvolution in order to bring out fine details.  Figure 6 below is an example of this of using 6 10 min subs for which optimal filtering has been applied. This methodology only effects high spatial frequency noise.  Lower spatial frequency noise should be reduced by dithering.

Figure 6  Integration with 6 subs with optimal filtering and a bit of deconvolution.

References

1  William H. Press, Saul A. Teukolsky, William T. Vetterling and Brian P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, (Cambridge University Press, New York 1992).

2 I Trujillo,J.A.C. Aquerri, J Cepa and C. M. Gutierrez, The effects of Seeing on Sersic Profiles. II. The Moffat PSF, Mon Not R Astron. Soc. 000 1-7 (2001).