Revision 1, 2000 May 20
Many of the methods used to process star images assume that the image is well sampled - at least 4 points within the FWHM. The MK IV lens is now giving an FWHM not much greater than one pixel and Tom is liable to improve it further! Therefore, we have to look at the problem of undersampling.
If one samples a band-limited signal sufficiently frequently, one can (in principle) recover all the information. An optical system produces a band-limited signal. A diffraction limited f4 system requires that the image be sampled at interval 4 lambda; for a wavelength, lambda, of 0.6E-6 m, this is 2.4E-6 m. Oops! The CCD pixel size is around 15E-6 m. Therefore, as David McCain among others pointed out, the image would be grossly undersampled if the lens were diffraction limited. Fortunately (or otherwise) our lenses are not, but we still have a problem.
As an example of undersampling, it is instructive to consider a one dimensional gaussian profile. This is not a band-limited signal but the corresponding spectrum, also gaussian and easy to calculate (why do you think I picked on a gaussian profile?), falls off rapidly at high frequency. A gaussian of width parameter s (FWHM = 2sqrt(2ln(2))s = 2.35s), exp(-0.5(x/s)^2), has a spectrum exp(-0.5(2pi s f)^2) where I have left out all those nasty numerical factors with sqrt(pi) in them that I usually get wrong. Putting in the numbers, using my trusty Quattro-Pro spreadsheet, FWHM = 1.5 pixels gives less than 0.5% of the power above the folding frequency of 0.5 cycles/pixel. Not band-limited but quite close. Of course, Tom is hard at work improving things; let's try FWHM = 1.0 pixel. Um! 6% of the signal power above the folding frequency. Rapidly getting not nice.
One major factor has been left out. The CCD pixels are not points.
They occupy a fraction Alpha; it should be quite a good approximation
to take Alpha = 1. The observed signal is the actual image
convolved with the pixel. In frequency space, this multiplies the
observed spectrum by Sin(Alpha pi f)/(Alpha pi f). With Alpha = 1,
FWHM = 1, the power above the folding frequency is reduced by about a
factor of 4. Better than this, with Alpha = 1, Sin(pi f) = 0 for f = 1,
2 ... cycle/pixel. These are the frequencies which alias into f = 0;
all the aliasing errors at f = 0 vanish for any image whatsoever!
[Consider a delta function PSF; a star falls on one pixel or the next.
With Alpha = 1 there are no gaps so one always measures the intensity
correctly.]
This should mean that we can get good magnitudes even when
Tom reduces the lenses to their diffraction limit! Note that this
dispensation applies only for magnitudes, which depend on the
spectrum at f = 0, not to the positions, since these depend on
rate of change of phase of the spectrum around f = 0 and thus include
frequencies for which the aliasing does not vanish identically.
But enough of this spectral stuff! And note that weasel phrase "in principle" back there in the first line. What happens in the real world?
Back to the spreadsheet. I have computed the aperture sum magnitudes for different pixel phasings for useful values of FWHM and Alpha. Alpha = 0.01 is pretty much what happens with point pixels; my spreadsheet blows up if I put in Alpha = 0. CCDs don't work too well with the active area that small either.
Aperture Sum Maximum error in magnitude: 1D case Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.062 0.041 0.026 0.000 1.1 0.029 0.020 0.012 0.000 1.2 0.013 0.009 0.005 0.000 1.4 0.002 0.001 0.001 0.000 1.6 0.000 0.000 0.000 0.000 2 0.000 0.000 0.000 0.000
The error thus goes to zero at Alpha = 1 for any FWHM, as the theory said it should. For Alpha = 0.9 (1.5E-6 m gaps with 15E-6 m pixels; is this a reasonable guess?) the error is going bad very rapidly as the FWHM approaches 1 pixel. The maximum error in two dimensions is just twice that in one so the error for Alpha = 0.9 and FWHM 1 pixel is 0.052 magnitudes.
Of course, this is for the sum over a "sufficiently large" aperture; this does not include those nasty edge effects for small apertures that were discussed elsewhere. In the real world, with noise, one cannot take the aperture large enough to get rid of the edge effects without at the same time bringing in intolerable noise. Ah well - all things in this life involve compromise.
I have computed position errors using a method approximately equivalent to that used by Star. Star uses iterative cubic interpolation to find the median of the marginal sums. I have used the same cubic interpolation but have used a non-iterative method to find the median approximately - I am too lazy to set up over 500 iterative solutions as I can't get Quattro-Pro to do more than one at a time!
Maximum error in position, pixels (1D: approximate) Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.146 0.141 0.138 0.133 1.1 0.121 0.116 0.113 0.109 1.2 0.100 0.097 0.094 0.091 1.4 0.071 0.070 0.069 0.068 1.6 0.055 0.054 0.054 0.053 1.8 0.044 0.044 0.044 0.044 2 0.037 0.037 0.037 0.036
So the theory is confirmed; magnitudes are OK but positions are lousy. In two dimensions the position error is the same as the one dimensional case in both dimensions. The maximum radial error is thus root two times larger than the value in the table.
Note that again, as in the case of the photometry above, the aperture is assumed large and the estimated errors do not include the edge effects of a small aperture.
Now: what happens to my least squares fitting program when it is given undersampled data? I set it up to manufacture synthetic stars of various FWHMs using a 2D Gaussian viewed with various pixel fractions Alpha and let it rip, fitting a circular Gaussian, positions in both axes and amplitude. The first time, I let it fit a different Gaussian for each case. Not a good idea! Fitting of single undersampled sources gets, shall we say, a bit ratty. But the fits did converge and the errors went in all directions: the means were very close to the expected values. So for the results presented here, the width parameter of the fitted Gaussian was held at the expected value and only the two positions and the amplitude were fitted. This is not very different from what I do when fitting a real image - the individual source shape fits are averaged over the whole image to give stable estimates.
Note that the fits, including the adjustable width parameter, did converge. This contradicts the result of Stone (AJ 97, 1227; 1989) who found that his fit became undetermined for an FWHM about 2. I am using the Simplex code, similar to the Amoeba routine in Numerical Methods. This is a very robust method. I also transform my variables so that, for example, FWHM = 0 is transformed to minus infinity. This certainly helps avoid numerical problems with physically silly parameters. Anyway, my code gave reasonable looking results down to FWHM = 1.
Least Squares Fit: Gaussian Maximum error in magnitude: 2D case Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.000 0.042 0.074 0.125 1.1 0.000 0.024 0.041 0.068 1.2 0.000 0.014 0.023 0.038 1.4 0.000 0.004 0.007 0.012 1.6 0.000 0.001 0.002 0.004 1.8 0.000 0.000 0.001 0.001 2 0.000 0.000 0.000 0.000
Not too surprising! The fit is optimised for Alpha = 0 (pure Gaussian) and gets worse as Alpha increases. This is just the opposite of the Aperture Sum above. Bearing in mind that factor of two for going from 1D to 2D, the error for the Aperture Sum at Alpha = 0.01, 0.124 magnitudes at FWHM = 1, is remarkably close to the error in the Least Squares Fit at Alpha = 1, 0.125 magnitudes. And I can do better! Instead of fitting a pure Gaussian I can fit a Gaussian convolved with a pixel having Alpha = 0.9 gives the following results.
Least Squares Fit: Gausian + pixel fraction Alpha0 = 0.9 Maximum error in magnitude: 2D case Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.073 0.031 0.000 0.050 1.1 0.041 0.017 0.000 0.027 1.2 0.023 0.010 0.000 0.015 1.4 0.007 0.003 0.000 0.005 1.6 0.002 0.001 0.000 0.001 1.8 0.001 0.000 0.000 0.000 2 0.000 0.000 0.000 0.000
For comparison here are the Aperture Sum results again; these are the 2D version for direct comparison
Aperture Sum Maximum error in magnitude: 2D case Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.124 0.083 0.051 0.000 1.1 0.059 0.039 0.024 0.000 1.2 0.026 0.017 0.011 0.000 1.4 0.004 0.003 0.002 0.000 1.6 0.000 0.000 0.000 0.000 1.8 0.000 0.000 0.000 0.000 2 0.000 0.000 0.000 0.000
If the PSF is well characterised, therefore, the Least Squares Fit can do rather well. If, for example, Alpha could range over 0.8 to 1.0 and the Least Sqares Fit were done for Alpha = 0.9, the peak error would be 0.050 magnitudes at FWHM = 1, Alpha = 1. The Aperture Sum under the same conditions gives a peak error of 0.083 magnitudes at FWHM = 1, Alpha = 0.8.
Least Squares Fitting is not constrained to use a small fitting box for noisy data and therefore does not suffer from the box edge effects which are a feature of th Aperture Sum method.
Least Squares Fit: Gaussian Maximum radial error (2D), pixels Pixel fraction, Alpha 0.01 0.8 0.9 1 FWHM 1 0.000 0.009 0.013 0.020 1.1 0.000 0.008 0.013 0.021 1.2 0.000 0.007 0.011 0.017 1.4 0.000 0.003 0.006 0.009 1.6 0.000 0.002 0.003 0.004 1.8 0.000 0.001 0.001 0.002 2 0.000 0.000 0.000 0.000
This is already about a factor ten better than the Median method using cubic interpolation. Again, matching the actual value of Alpha gives a further reduction in error.
I conclude that the undersampling resulting from an FWHM as low as 1.0 pixels does not present any problem provided that we use properly chosen algorithms for data reduction. The Aperture Sum estimate of amplitude is essentially unaffected, provided the CCD pixels occupy a substantial fraction of the CCD area, which we believe they do. Least Squares estimation remains accurate provided one has a good estimate of the PSF. The Median method of position estimation is badly affected by undersampling (as is already well known) and should not be used with such a small FWHM. Other methods of position estimation are available which are adequate to our needs; Least Squares estimation (using a good PSF estimate) works well.