Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mrao.cam.ac.uk/projects/OAS/publications/fulltext/Gordon_SPIE_7734-157.pdf Дата изменения: Wed Feb 20 16:47:57 2013 Дата индексирования: Sat Mar 1 03:30:13 2014 Кодировка: Поисковые слова: interferometry

 

















































Bias-free imaging at low light levels
James Gordona , David Buschera and Hrobjartur Thorsteinsson
a a

Astrophysics Group, Cavendish Laboratory, University of Cam bridge, J.J. Thomson Avenue, Cambridge CB3 0HE, UK
ABSTRACT

Measurements from long-baseline interferometry are commo nly analysed in terms of the power spectrum and the bispectrum (or triple-product) of the fringe patterns, as these es timators are invariant in the presence of phase instabiliti es. At low light levels, photon and detector noise give rise to systematic "bias" in the power spectrum and bispectrum. This paper extends previous work on computing the expected biases and variances for these quantities by introducing a general method which can be applied to any fringe-encoding scheme where the measurement equation is linear and to measurements affected by a combination of photon noise and detector noise. We apply our method to a number of interesting practical cases, including systems with unevenly-sampled fringe patterns and in the presence of read noise. Keywords: optical interferometry, fringe decomposition, statistic al bias, bias correction, read noise

1. INTRODUCTION
In both long-baseline stellar interferometry starlight co llected by different parts of an aperture is combined to form an interference pattern. This interference pattern is then analyse d to derive information about one or more Fourier components of the angular brightness distribution of the object being obs erved, and this Fourier data can then be used in the reconstru ction of a high-angular-resolution image of the target. Each Fourier component is a complex quantity with both an amplitude and phase. Under phase-unstable conditions the phase of the measured Fourier component is corrupted in a ran dom and time-dependent fashion. In order to recover the object amplitude and phase information under these conditions, it is common to take multiple short-exposure measureme nts of the fringe pattern and to compute two quantities from each short exposure, the power spectrum and the bispectrum. Thes e quantities are then averaged over multiple exposures -- perh aps thousands or tens of thousands of exposures. While the phases of the short-exposure Fourier components change ran domly from exposure to exposure, the phases of the power spectrum and bispectrum are immune (to first order) to the atm ospheric and instrumental phase perturbations and hence retain information about the object when averaged. The aver aged power spectrum and bispectrum information can therefo re be used to derive reliable images using iterative reconstru ction techniques. The analysis of the short-exposure interference pattern to produce the bispectrum or power spectrum estimate is of necessity a non-linear operation. This non-linearity comb ines with the statistics of the measurement noise, whether due to photon noise or detector noise, to produce biases in the aver aged quantities which become pronounced at low light levels. Goodman1, 2 derived an expression for the bias in the power spectrum and Wirnitzer3, 4 extended Goodman's method to derive the bias of the bispectrum. In both cases the authors used the expressions for the bias to derive bias-free estimat ors for these quantities. The technique used by these authors made two important assumptions. Firstly, Goodman's technique assumes a photoncounting detector and does not address situations where det ector noise is present in addition to photon noise. Secondly it assumes that the amplitudes and phases of the Fourier components incorporated into the power spectrum and bispectrum estimators are to be derived from a continuous Fourier trans form of the interference pattern intensity. This technique can be used with sampled data by substituting a Discrete Fourier Transform (DFT) for the continuous transform, but for the DFT to give an accurate representation of the continuous tra nsform result, the interference pattern needs to be both eve nly sampled and sampled over a number of "pixels" which is an integer number of cycles of the Fourier frequency being measured. Under low-light-level conditions it is often pre ferable to minimise the number of pixels sampling the pattern in order to minimise the effect of detector read noise, as well as to speed up detector rea d-out, which could lead to violations
Further author information: j.gordon@mrao.cam.ac.uk


of the assumed sampling conditions for instrumental reason s. In addition, further processing of the data may be necessary before Fourier transforming to compensate for pixel-to-pi xel gain variations of the detector. Beletic et al5 used a modified version of Goodman's method to derive the bias and variance in the power spectrum, accounting for pixel-to-pixel gain variations in a detector. Le Bouquin et al6 introduced a generalised "pixel-to-visibility" transformation matrix for deriving interference pattern e stimates without the restrictions of a DFT. Using this formalism Tatulli et al7 showed how to compute the bias to the power spectrum under the combined effects of photon noise and detector readout noise, but did not provide equivalent expr essions for the bispectrum. In this paper, we present a method for deriving the biases of t he power spectrum and bispectrum without the restriction to the use of a DFT or only one kind of noise. We show how our results allow improved estimates of the power spectrum and bispectrum under a number of scenarios of practical inte rest.

2. PIXEL TO VISIBILITY TRANSFORM MATRIX
Let us first consider the ideal situation where we have a detec tor system without noise. The continuous fringe pattern created by combining light from multiple telescopes is reco rded by a discrete number of pixels on the detector such that the number of photons, n, falling on each pixel, p, gives the noise-free interferogram p . For the special case where the discrete Fourier transform (D FT) conditions hold, ie evenly sampled pixels, and an integer number of cycles of the Fourier frequency across the detector, we can calculate the complex visibility, C i j , from a given baseline, i j, using ij p H ipj , (1) p ei2 p f = Ci j =
p p

where p is the fringe modulation for temporally sampled fringes and f i j is the fringe frequency from baseline i j. We also introduce the matrix H ipj , which is identical to the DFT under the conditions mentioned above. For instances where the DFT conditions are not met, we modify the matrix H ipj to allow the correct computation of C b . Let us first express Eqn. 1 in matrix form C = H, where C C , R and H C characterise H.
b p bв p

(2)

and b is the number of baselines. In order to solve this equation, we need to

With knowledge of the interferometer, we can construct the s eries of linear equations = W C, (3)

where W C pвb contains the properties of the interferometer. By solving t he inverse problem, we can calculate the elements of H We use the Singular Value Decomposition method to find the optimum (in terms of least squares) pseudoinverse of W .

3. NOISE MODEL
The signal recorded by a real detector system will always have some noise component, especially at low light levels. We rewrite Eqn. 1 in terms of the underlying interferogram, p , embedded in noise -- which we name i p -- to give the noisy complex visibility on a given baseline i p H ipj . (4) ci j =
p

We consider two dominant noise sources; photon noise and rea d noise. Photon counting noise is present in the noisy interferogram , i p . Taking the ideal interferogram, p , such that for a given pixel value, p, the probability of recording N p photons is given by the continuous Poisson distribution


P

P

N p |

p

=
n

p p Np - n exp - p . Np!

N

(5)


Read noise, , is often quoted as a constant across pixels for modern detec tors, however it is sometimes preferable to allow to vary on a pixel by pixel basis as p and we will present results in this form for the sake of genera lity. Assuming the read noise is independent of the expected number of photon counts, the probability of measuring N p photons on a given pixel, i p , is given by the Gaussian distribution Np - ip 2 1 . exp - (6) PG N p - i p | p = 22 2 p 2 p Other Sources of Noise include detector dark current, D p , is also described by the Poisson distribution and could be included in this formalism by replacing p with p + D p . It is noted that even for very low SNR data, dark current prod uces a negligible effect on modern detectors and for the sake of readability will n ot be explicitly included in later derivations. To incorporate the shot noise and read noise into a flexible no ise model, we simulate the Poisson noise model convolved with the fixed Gaussian model. This means the probability of obtaining the noisy data, i p , from the underlying ideal data, p , is given by
i
p

P i p | p , =
0

P

P

N p |

p

PG N p - i p | p d N p .

(7)

3.1 Moments of the Noise Model
We apply the statistics of our noise model by considering the moments of the probability distribution. We require the moment generating functions MP and MG of the noise terms, and combine them using the convolution theorem (applicable since the probability density and moment generating functions are a Laplace transform pair). This give the moment generating function 2 2 (8) exp - + e , M () = MP () MG () = exp 2 for our noise model probability density P i p | p , . We give the first three moments of the noise model required to derive the results given in this paper. Presented are the noisy inte rferograms averaged over all frames, denoted . . . . i i i
2 p 3 p p

d M ( ) d d2 = 2 M ( ) d d3 = 3 M ( ) d =

=0

= = =

p, 2 + p + 2 , p p 3 + 32 + p + 3 p 2 . p p p

(9) (10) (11)

=0

=0

4. BIAS FREE POWER SPECTRUM ESTIMATOR
The noise free power spectrum, S i j , can be expressed in terms of the noise free interferograms, S ij = C
ij 2

=
p1 p2

j p1 p2 H ipj1 H pi2 .

(12)

We require an estimator of S i j that does not suffer from statistical bias introduced by the noisy interferog ram, i p . This estimator was calculated by Tatulli,7 and we present the derivation here as an example. We start by w riting Eqn. 12 in terms of the averaged noisy interferogram ci
j2

=
p1 p2

i p1 i p2 H ipj1 H

ji p2

.

(13)

We split the sum for the cases where p1 = p2 ( p) and for p1 ci
j2

p2 and rewrite Eqn. 13 to give i
p1

=
p

i

2 p

H

ij 2 p

+
p1 p2

i

p2

j H ipj1 H pi2 ,

(14)


noting that i p1 i leads to

p2

=i

p1

i

p2

for p1 ci
j2

p2. Substitution of the moments of the noise model given in Eqns. 10 and 11 2 + p + p
p 2 p

=

H

ij 2 p

+
p1 p2

j p1 p2 H ipj1 H pi2 .

(15)

The equivalent split form for the noise-free power spectrum given in Eqn. 12 is S ij =
p

2 H p

ij 2 p

+
p1 p2

j p1 p2 H ipj1 H pi2 .

(16)

Substituting Eqn. 16 into Eqn. 15 gives S i j = ci
j2

-
p

ip +

2 p

H

ij 2 p

.

(17)

The estimator, S i j , is then given by S i j = ci as presented by Tatulli.
7 j2

-
p

ip +

2 p

H

ij 2 p

(18)

5. BIAS FREE BISPECTRUM ESTIMATOR
We can calculate the bias free bispectrum estimator followi ng the same methodology as for the power spectrum. Defining the bispectrum j (19) Bi jk = C i jC jk C ki = p1 p2 p3 H ipj1 H pk H ki3 , 2p
p1 p2 p3

where the index i j relates to the baseline created between apertures i and j. The bias free estimator, Bi jk is given by Bi
jk

= ci j c jk cki - ci j ip +
p

2 p 2 p

j H pk H

ki p ki p

-

c

jk p

ip + ip +
p

H ipj H H ipj H

-c +

ki

2 p

jk p

2i p + 3
p

2 p

j H ipj H pk H ki . p

(20)

5.1 Comparison with pervious bispectrum estimators
j If we assume DFT conditions, H pk H ki = H ji , and that read noise is negligible ( p = 0), we reduce Bi jk to the estimator p presented by Wirnitzer3

Bi jk Wirnitzer

= ci j c jk ck - ci
j2

i jk 2

-c

- ck

i2

+ 2 N,

(21)

where N =

(as this goes as N ), and we reduce our estimator further to the estimator used b y Tatulli,8 Bi jk Tatulli = ci j c jk cki .

p ip 3

is the mean number of photons in the interferogram. For high l ight levels, the ci j c jk cki term dominates


0.10 a Closure Phase Bias, bias [ radians]

0.04 Closure Phase Bias, bias [ radians]
Our Estimator Tatulli Estimator

b

Our Estimator Tatulli Estimator

0.05

0.02

0.00

0.00

-0.05

-0.02

-0.10 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

-0.04 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

Figure 1. The bias in closure phase estimated for 40 and 100 photons per aperture per exposure, in frames a and b respectively, for an average of 105 exposures under DFT conditions and zero read noise. The trace for th e closure phase from the Wirnitzer estimator is not shown, as it it identical to that of our estimator in the zero read noise limit. Note that the bias is much larger than the variance for this number of exposures.

6. CLOSURE PHASE BIAS FROM BISPECTRUM ESTIMATORS
As an example of the effects of statistical bias, we modeled interferograms with re alistic noise and the deviations from DFT conditions. We present the closure phase data extracted usi ng our new estimator, the Wirnitzer estimator, and the Tatul li estimator. We see that the closure phase can be significantly effected by statistical bias at low light levels.

6.1 Simulation model
We simulated an ideal fringe pattern from a three aperture in terferometer creating a noise-free fringe pattern, I p , such that
N
ap

N

ap

Ip =
p

N + 2N V

pq i< j

cos 2 p f i j +

ij

,

(22)

is the number of apertures, N the flux per aperture, V i j = 1/Na p the visibility of the baseline between j, p the sampling modulation, f i j fringe frequency on baseline i j, and i j the object phase on baseline i j j = C i j . Poisson shot noise is added to each pixel of the interferogram for all of the simulations. We then of adding read noise, 2 , and adjusting p to deviate from DFT conditions. We simulate objects with closure phase, true , between - and such that true = i j + jk + ki where two of these phases are set at random and the third phase component select ed to get the desired true . where Na p = 3 telescopes i and i such that V i j ei have the option The expected bispectrum is calculated using our estimator, closure phase is extracted from the bispectrum estimate in e ach Im = atan Re the Wirnitzer estimator and the Tatulli estimator, and the case using Bi jk . (23) i jk B

We present the closure phase bias, bias = - true , after applying an algorithm to remove phase wrapping.

6.2 DFT conditions and zero read noise
In the under shows Wirnit true situation where photon noise dominates the noise model (ie read noise is zero) and the interferometer is working DFT conditions, our bispectrum estimator is identica l to the Wirnitzer estimator as shown in Eqn. 5.1. Fig. 1 the returned closure phase estimate for 40 and 100 phot ons per aperture per exposure. Our estimator and the zer estimator show no bias with respect to object phas e, however the Tatulli estimator shows a bias term, peaking a t -/2, /2. This is a well known result, and stimulated the creation of the Wirnitzer estimator.


0.4 Closure Phase Bias, bias [ radians]

a Closure Phase Bias, bias [ radians]

0.04

b

Our Estimator Wirnitzer Estimator Tatulli Estimator

0.2

0.02

0.0

0.00

-0.2
Our Estimator Wirnitzer Estimator Tatulli Estimator

-0.02

-0.4 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

-0.04 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

Figure 2. The same simulation model as in Fig. 1 but with the introduction of re alistic read noise = 3e- . Frame a and b are again for 40 and 100 photons per aperture per exposure respectively.
0.2 a Closure Phase Bias, bias [ radians]
Our Estimator Wirnitzer Estimator Tatulli Estimator

0.3 b Closure Phase Bias, bias [ radians]

Our Estimator Wirnitzer Estimator Tatulli Estimator

0.1

0.2

0.0

0.1

-0.1

0.0

-0.2 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

-0.1 -1.0

-0.5 0.0 0.5 Object Closure Phase, true [ radians]

1.0

Figure 3. The same simulation model as in Fig. 1 but with the introduction of non-DFT conditions. We simulate non-linear fringe modulation as described in Eqn. 24. Frames a and b are again for 40 and 100 photons per aperture per exposure respec tively. We offset the closure phase bias of the Tatulli estimator in frame b by 0.15 for clarity (the zero being shown again by the dashed line).

6.3 DFT conditions but non-zero read noise
At low light levels the read noise on the detector becomes important. Fig. 2 shows simulations under DFT conditions, as in Fig. 1, but with = 3e- read noise. This level of read noise is quite low compared to m any of the optical detectors in use today. We see that failing to account for read noise introduc es a significant bias in closure phase estimates. For very low light levels, for example 40 photons per exposure per apertu re as seen in Fig. 2 frame a, the Wirnitzer estimator introduces more bias than the Tatulli estimator.

6.4 non-DFT conditions, zero read noise
We simulate a small deviation from DFT conditions by adjusting p , the fringe modulation, which results in unevenly sampled fringes for temporal fringe encoding schemes. The i deal case is for p to be a triangle wave described by n = 1, 3, . . . , if DFT conditions, (-1)(n-1)/2 8 n p p = 2 (24) s in n = 1, 3, 5 L n n2 if non-DFT conditions where L is the period. It is clear however that no interferome ter is able to create this ideal modulation, and that this summation will always be cropped at some cut-off frequency. We convolve Eqn. 24 with low pass filter to remove a ll terms above n = 5 to simulate a more probable fringe modulation. Fig. 3 shows the closure phase returned by the three estimato rs. In this case we simulate our estimator using the correct version of the matrix H and the Wirnitzer and Tatulli estimators assuming DFT conditions. Our estimator and the


Wirnitzer estimator both show no bias, however the variance of the Wirnitzer estimator is significantly larger than that of our estimator. The Tatulli estimator shows the same bias as f or DFT conditions, but also the large variance seen in the Wirnitzer estimator. We note here that using the Wirnitzer a nd Tatulli estimators with the correct version of the matrix H gives results identical to using the estimators under DFT conditions.

7. CONCLUSIONS
We have presented bias free estimators for both the power spe ctrum and the bispectrum that work in the presence of both read noise and photon noise under general, non-DFT conditions. The estimators reduce to the previous estimators under special cases and are given in a form that is simple to implement into existing data-reduction pipelines. The bias correc ted bispectrum estimator presented in this paper will allow rob ust and reliable image reconstruction -- essential when work ing with sparse uv coverage as is always the case in optical interferometry.

REFERENCES
[1 ] Goodman, J. W., Belsher, J. F., and LAB, S. U. C. I. S., "Photon limited images and their restoration," (1976). [2 ] Goodman, J. W., Belsher, J. F., and CA, S. U., "Precompensation and postcompensation of photon limited degraded images," (1976). [3 ] Wirnitzer, B., "Bispectral analysis at low light levels and astronomical speckle masking," Journal of the Optical Society of America A 2, 14­21 (Jan. 1985). [4 ] Karbelkar, S. N. and Nityananda, R., "Atmospheric noise on the bispectrum in optical speckle interferometry," Journal of Astrophysics and Astronomy 8, 271 (Sept. 1987). [5 ] Beletic, J. W., "Deterministic photon bias in speckle imaging," Optics Communications 71, 337­340 (June 1989). [6 ] Le Bouquin, J. B. and Tatulli, E., "Pupil plane optimization for single-mode multi-axial optical interferometry wi th a large number of telescopes," Monthly notices of the royal astronomical society 372(2), 639­645 (2006). [7 ] Tatulli, E. and Duvert, G., "AMBER data reduction," New Astronomy Reviews 51(8-9), 682­696 (2007). [8 ] Tatulli, E., Millour, F., Chelli, A., Duvert, G., Acke, B., Hofmann, K. H., Kraus, S., Malbet, F., Mege, P., and Petrov, R. G., "Interferometric data reduction with AMBER/VLTI. principle, estimators and illustration.," Arxiv preprint astroph/0603046 (2006).