A&A 398, 385-390 (2003)
DOI: 10.1051/0004-6361:20021667

Subtracting the photon noise bias from single-mode optical interferometer visibilities

G. Perrin

LESIA, FRE 2461, Observatoire de Paris, section de Meudon, 5 place Jules Janssen, 92190 Meudon, France

Received 3 July 2002 / Accepted 22 October 2002

Abstract
I present in this paper a method to subtract the bias due to source photon noise from visibilities measured with a single-mode optical interferometer. Properties of the processed noise are demonstrated and examples of subtraction on real data are presented.

Key words: instrumentation: interferometers - techniques: interferometric - methods: data analysis

1 Introduction

The properties of source photon noise are well known. It follows a Poisson distribution whose variance is equal to the average total number of photons. In frequency space, it is a white noise with a flat average power spectral density. In most practical cases where the observables are linearly linked to the number of photons detected, photon noise can be directly averaged out from the data to reduce its variance. For some applications for which the observables are quadratically linked to the number of photons, the data suffer both from photon noise and from a bias linked to the variance of the noise. This is for example the case in speckle imaging techniques where the source spatial intensity distribution is recovered from the power spectral density of a short time exposure (Thiébaut 1994). In astronomical optical interferometry, the observables are the modulus of the visibility and its phase usually expressed as a closure phase quantity. The visibility modulus can be obtained by integrating the modulus of the spectrum of interferograms. However, this estimator is biased by the power spectral density of noises as these add to the power spectral density of the fringe signal. An unbiased estimator of the modulus of the visibility is obtained by forming the squared modulus of the visibility as the power spectral densities of the noises can be independently estimated and subtracted. An example in interferometry is the computation of the fringe squared visibility in the ABCD method where the white light fringe is sampled at four $\lambda/4$ spaced optical path differences. An unbiased single fringe ABCD estimator is obtained by subtracting the source photon noise and the detector noise variances (when the noise has a flat spectrum the power spectral density is constant and equal to the variance) from the fringe power spectral density (Tango & Twiss 1980).

In single-mode interferometers, beams are spatially filtered by single-mode waveguides trading phase fluctuations against intensity fluctuations. A fraction of intensities collected by each aperture can be measured to renormalize interferograms to eliminate the fluctuations due to turbulence. The visibility estimator is no longer directly linked to the power spectral density of the fringe signal. I demonstrate in the following sections that the classical method (explained further in the paper) established by Goodman (1985) can be rigourously extended to such ratios of physical noisy signals under certain assumptions to provide unbiased visibility estimators. Real data reduction cases are presented to illustrate the method.

2 Principles of photon noise bias subtraction

2.1 Description of the signals and assumptions

I refer the reader to Coudé du Foresto et al. (1997) for a full description of the principle to measure fringe amplitudes (also called coherence factors) with a single-mode fiber interferometer and for a detailed description of the fringe signal for coaxial interferometers. Visibility calibration will be addressed in a separate paper (Perrin 2002). I assume here, for sake of simplicity, that no calibration is required and the coherence factor is directly equal to the visibility. Here I will use the more general expression of the interferogram for a two-telescope interferometer:

$$i(x)=P_{\rm A}(x)+P_{\rm B}(x)+2\sqrt{P_{\rm A}(x)P_{\rm B}(x)}~m(x)$$ (1)

where m is an oscillating function containing the fringes. $P_{\rm A}$ and $P_{\rm B}$ are the intensities coupled in the single-mode waveguide at each telescope and are called the photometric signals. i is a function of the optical path difference or of time for coaxial beamcombiners. For multiaxial interferometers, x is a focal plane coordinate. In the following, I will deal with coaxial interferometers only and I will use time t as the variable. The method can be easily adapted to multiaxial interferometers. The photometric signals vary in time with turbulence and the modulation m varies with turbulence and with the optical difference which is varied linearly with time, coding the fringe signal in frequency space. i is measured in photon counts and is defined as the average signal one would obtain if no noise were present.

In the photometric calibration method, the photometric signals need to be estimated. The estimated signals are filtered, the filtering function being adjusted to reject most of the noise and keep the intensity fluctuations due to turbulence only. Turbulent fluctuations are low frequency fluctuations (limited to frequency ranges of a few tens of Hertz). I note $\overline{P_{\rm A}}$ and  $\overline{P_{\rm B}}$ the estimated photometric signals suitably low-pass filtered. The important property of the filtered photometric signals is that they contain no energy at the fringe frequency and above.

I call g the following gain function:

$$g(t)=\frac{1}{2\sqrt{\overline{P_{\rm A}}(t)\overline{P_{\rm B}}(t)}}$$ (2)

I define the normalized interferogram:

 

                             i_n(t) = g(t).i(t) (3)

$$\frac{P_{\rm A}(t)+P_{\rm B}(t)}{2\sqrt{\overline{P_{\rm A}}(t)\o... ...}(t)P_{\rm B}(t)}}{\sqrt{\overline{P_{\rm A}}(t)\overline{P_{\rm B}}(t)}}~m(t).$$ =

The first term is mainly a low frequency signal whereas the second term is the high frequency signal containing the fringe modulation. I introduce the continuum function:

$$c(t)=\frac{\overline{P_{\rm A}}(t)+\overline{P_{\rm B}}(t)}{2\sqrt{\overline{P_{\rm A}}(t)\overline{P_{\rm B}}(t)}}\cdot$$ (4)

The continuum function is the low frequency part of the normalized interferogram. It is the ratio of the arithmetic and geometric means of the photometric signals. If the two photometric signals are equal then the continuum function is equal to 1. It departs all the more from 1 as the photometric beams get unbalanced. The visibility estimate is computed from the corrected interferogram:

 

$$i_{\rm cor}(t)$$ = i_n(t)-c(t) (5)

$$\frac{\left[P_{\rm A}(t)-\overline{P_{\rm A}}(t)\right]+\left[P_{... ...}(t)P_{\rm B}(t)}}{\sqrt{\overline{P_{\rm A}}(t)\overline{P_{\rm B}}(t)}}~m(t).$$ =

In the corrected interferogram the low frequency components due to turbulence are eliminated. If the photometric signals are perfectly estimated then the average value of the corrected interferogram is equal to zero and the oscillating signal is properly renormalized. In the ideal case where data are noiseless, the corrected interferogram is equal to the oscillating function m which is proportional to the visibility. Although the filtered photometric signals do not contain energy at the fringes frequency, the gain and continuum function may contain some as the residual noise may have been redistributed in the frequency domain by the non-linear combinations of the filtered photometric signals. I make the assumption that this high frequency noise is negligeable by several orders of magnitude compared to the energy of the fringes. This assumption will be validated in Sect. 3 on real data. As a consequence of this assumption, i_n and $i_{\rm cor}$ share the same photon noise. This assumption is crucial for the success of the method. Finding an analytical solution to the problem of bias subtraction without this assumption is certainly a big issue.

2.2 Method

The method to subtract the photon noise bias is a direct generalization of that proposed by Goodman (1985) and I will use the same notations. In the following only photon noise is considered. Methods to subtract additive noises are well established (see for example Coudé du Foresto et al. 1997). Detector and source photon noises being independent the detector noise could be added in the following derivation leading to the classical result on power spectral density bias by detector noise variance. The primary scope of this paper being source photon noise I have decided to keep equations as light as possible and not include detector noise in the equations. Detector noise will nevertheless be considered in the last section on real data.

I call $\tilde{\imath}(t)$ a representation of the interferogram in which photon events are represented by Dirac functions. Assuming that the total number of photons in the interferogram is $\tilde{N}$ for a given realization, then I can write:

$$\tilde{\imath}(t)=\sum_{k=1}^{\tilde{N}}\delta(t-t_{k}).$$ (6)

There are two random variables in this expression: the individual photon detection times t_k and the total number of photons $\tilde{N}$ which is equal to N in average. A model of the normalized interferogram can be built:

$$\tilde{\imath}_{n}(t)=\sum_{k=1}^{\tilde{N}}g(t)\delta(t-t_{k})=\sum_{k=1}^{\tilde{N}}g(t_{k})\delta(t-t_{k}).$$ (7)

Its Fourier transform is therefore:

$$\tilde{I}_{n}(f)=\sum_{k=1}^{\tilde{N}}g(t_{k}){\rm e}^{-2i\pi t_{k}f}.$$ (8)

Thus, the average spectrum of the normalized interferogram is therefore:

$$\langle \tilde{I}_{n}(f)\rangle \left\langle \left\langle \sum_{k=1}^{\tilde{N}}g(t_{k}){\rm e}^{-2i\pi t_{k}f}\right\rangle_{t_{k}}\right\rangle_{\tilde{N}}$$ = (9)

$$\left\langle \sum_{k=1}^{\tilde{N}}\left\langle g(t_{k}){\rm e}^{-2i\pi t_{k}f}\right\rangle_{t_{k}}\right\rangle_{\tilde{N}}\cdot$$ = (10)

In the above expression, the average on t_k does not depend upon k and the average normalized interferogram can be written:

$$\langle \tilde{I}_{n}(f)\rangle=\left\langle \tilde{N} \left\... ... e}^{-2i\pi tf}\right\rangle_{t}\right\rangle_{\tilde{N}}\cdot$$ (11)

The statistics of the photon arrival time t are described by the probability density function $\frac{i(t)}{N}$ with i(t) the average photon flux. Hence the average on the arrival times is equal to:

 

$$\langle g(t){\rm e}^{-2i\pi tf}\rangle_{t} \int_{-\infty}^{+\infty}g(t)\frac{i(t)}{N}{\rm e}^{-2i\pi tf}{\rm d}t$$ = (12)

$$\frac{1}{N}I(f) \star G(f).$$ = (13)

The $\star$ symbol indicates a convolution. Capital letters are used for Fourier transforms. Replacing the average on t by this expression, I obtain for the average normalized interferogram:

$$\langle \tilde{I}_{n}(f)\rangle \langle {\tilde{N}} \rangle_{\tilde{N}} \frac{1}{N}G(f) \star I(f)$$ = (14)

$$G(f) \star I(f).$$ = (15)

As a check, the expression of the average interferogram after applying the inverse Fourier transform to the above expression is:

$$\langle \tilde{i}_{n}(t)\rangle = g(t).i(t).$$ (16)

Let us consider $\tilde{I}_{n}^{(2)}(f)$ the power spectral density of the normalized interferogram with $\tilde{I}_{n}^{(2)}(f)= \vert\tilde{I}_{n}(f)\vert^{2}$. In the ideal case of non noisy data, the integral of the power spectral density is proportional to the squared visibility. The integral of the average power spectral density is:

$$\left\langle \tilde{I}_{n}^{(2)}(f) \right\rangle = {\left\... ... \right\rangle }_{t_{k},t_{l}} \right\rangle}_{\tilde{N}}\cdot$$ (17)

An important assumption on photon events to derive the Poisson statistics is that they are not correlated. Arrival times t_k and t_l are therefore not correlated as long as $k \neq l$. This important property is used to split the above expression in two terms:

$$\left\langle \tilde{I}_{n}^{(2)}(f) \right\rangle = {\left\l... ...{l})f}} \right\rangle}_{t_{k}} \right\rangle}_{\tilde{N}}\cdot$$ (18)

The first average is equal to:

$$\left\langle {\left\langle {\sum_{k=1}^{\tilde{N}} [g(t_{k})]^{2}} \right\rangle}_{t_{k}} \right\rangle_{\tilde{N}} {\left\langle { \sum_{k=1}^{\tilde{N}} {\left\langle [g(t_{k})]^{2} \right\rangle}_{t_{k}} }\right\rangle}_{\tilde{N}}$$ = (19)

$${\left\langle \sum_{k=1}^{\tilde{N}} { \int_{-\infty}^{+\infty}[g(t)]^{2}\frac{i(t)}{N}{\rm d}t }\right\rangle}_{\tilde{N}}$$ = (20)

$${\left\langle \tilde{N} { \int_{-\infty}^{+\infty}[g(t)]^{2}\frac{i(t)}{N}{\rm d}t }\right\rangle}_{\tilde{N}}$$ = (21)

$$\int_{-\infty}^{+\infty}[g(t)]^{2}i(t){\rm d}t$$ = (22)

where I have used once again the probability density of the photon statistics $\frac{i(t)}{N}$. The second average of Eq. (18) can be written as the sum of factors yielding:

$${\left\langle {\left\langle { 2\sum_{k=1}^{\tilde{N}}\sum_{l<... ...ngle }_{t_{k},t_{l}} \right\rangle}_{\tilde{N}} = \nonumber$$

$${\left\langle { 2\sum_{k=1}^{\tilde{N}}\sum_{l<k} \left\lang... ... t_{l}f}\right\rangle_{t_{l}}} \right\rangle}_{\tilde{N}}\cdot$$ (23)

The averages on t_l and t_k can be substituted by the expression of Eq. (13) yielding:

$${\left\langle {\left\langle { 2\sum_{k=1}^{\tilde{N}}\sum_{l<... ...ngle }_{t_{k},t_{l}} \right\rangle}_{\tilde{N}} = \nonumber$$

$${\left\langle 2\sum_{k=1}^{\tilde{N}}\sum_{l<k} \frac{1}{N^{2}} \vert I(f) \star G(f)\vert^{2} \right\rangle}_{\tilde{N}}$$ (24)

and:

$${\left\langle 2\sum_{k=1}^{\tilde{N}}\sum_{l<k} \frac{1}{N^{2}} \vert I(f) \star G(f)\vert^{2} \right\rangle}_{\tilde{N}} =$$

$$\left\langle\tilde{N}(\tilde{N}-1)\right\rangle_{\tilde{N}}\f... ...rt I(f) \star G(f)\vert^{2} = \vert I(f) \star G(f)\vert^{2}.$$ (25)

$\vert I(f) \star G(f)\vert^{2}$ is the power spectral density of the normalized interferogram. From the equation above I therefore derive an unbiased estimate of this quantity:

$$\vert I(f) \star G(f)\vert^{2}=\left\langle \tilde{I}_{n}^{(2... ...ight\rangle - \int_{-\infty}^{+\infty}[g(t)]^{2}i(t){\rm d}t .$$ (26)

Since the normalized and the corrected interferograms share the same noise bias, the integral term of the above equation is also the bias of the power spectral density of the corrected interferogram.

Figure 1: Examples of raw signals acquired for Mira in October 2000 with FLUOR. a) and b) The two photometric signals $P_{\rm A}(t)$ and $P_{\rm B}(t)$. c) Interferometric signal i(t). d) Power spectral density of the interferometric signal. Energy units are photon counts. The power spectral density is expressed in squared photon counts. The fringe peak is located between 350 and 400 Hz. Other peaks are due to the noise and disappear when averaging the power spectral densities. The source is resolved hence the fringe amplitude is small and hardly shows up in the interferogram signal. The fringe peak in the power spectral density is proportional to the squared modulus of the visibility which is obtained by integrating the peak.

Open with DEXTER

	Figure 2: a) Gain function g(t). b) Continuum function c(t). c) Power spectral density of the gain function. d) Power spectral density of the continuum function. The gain function is homogeneous to the reciprocal of a number of photon counts. The continuum function has no unit.
Open with DEXTER

2.3 Comments

The expression of the photon noise bias is very intuitive. i(t) is the average number of photons detected at time t. Being a Poisson statistics, it is also the variance of the photon noise. When the signal is multiplied by the gain g(t), the variance at time t becomes [g(t)]²i(t). Photon events at different times being uncorrelated, the total variance of the photon noise is therefore equal to the integral of the local variance.

Equation (26) demonstrates that the noise of the corrected interferogram remains a white noise whose mean power spectral density is constant. This property is also the result of the independence of photon events.

With the above result, the computation of the unbiased estimate of the visibility is easy. The bias is simply obtained by co-adding the individual counts of the normalized interferogram.

	Figure 3: The corrected interferogram and its power spectral density. The corrected interferogram being the ratio of photon counts has no unit.
Open with DEXTER

Figure 4: a) Average corrected interferograms power spectral density. b) Average processed dark current scans power spectral density. c) Difference of the two to suppress the bias due to detector noise. d) Difference of the two with photon noise bias subtracted. The energy peak at low frequency in the processed dark signal power spectral density is due to the multiplication of the dark signals by the gain function. This peak causes a drop when the dark signal power spectral density is subtracted from that of the fringe signal.

Open with DEXTER

3 Example of bias subtraction on real signals

The purpose of this section is to illustrate the theoretical results of the previous section and to show that the assumptions on the noise of the continuum and gain functions are correct.

I have selected a series of scans of Mira observed with FLUOR in October 2000. The observations were carried out with a baseline long enough that the fringe contrast is of a few percent only. This is an interesting case because the source is very bright and the fringe contrast is small, hence the photon noise bias is relatively important and accounts for a large fraction of the visibility if not corrected.

The raw photometric and interferometric signals are displayed in Fig. 1. The intensities are in photon counts. The fringes are hardly visible and their amplitude is smaller than that of the photometric fluctuations. The fringes peak is visible in the power spectral density between 350 and 400 Hz. In Fig. 2 are presented the gain and continuum functions as well as their power spectral densities on a log scale. The noise level of the continuum at the fringes frequency range can be directly compared to that of the corrected interferogram in Fig. 3. A residual photometric fluctuation is visible on the corrected interferogram hence the low frequency component. In the data reduction procedure, the detector dark current signals undergo the same correction process as the source signals. In case the dark currents would combine both detector dark signal and sky background if measured with a chopping technique, the dark signals would contain both background photon noise and detector noise. But both noises can be treated as additive noise, in particular the background photon noise does not need to be processed like the source photon noise. In the case of the FLUOR signals, background photon noise is totally negligible and the dark signals only contain the detector contribution. In Fig. 4 I have represented the average of power spectral densities of both the corrected detector dark current and corrected interferograms to reduce the noise on the power spectral densities. On the same figure, the difference of the two is plotted. At this stage, the bias due to source photon noise is not removed to make its magnitude obvious to the reader. The fringe squared visibility being obtained by integrating the power spectral density at the fringe frequency, the non-zero mean level under the fringe peak causes the bias. This is the bias due to source photon noise. This graph also shows an increasing residual noise at high frequency due to detector instabilities. The value of the bias has been computed with the method presented in this paper and subtracted as illustrated by the last graph of this figure. The power spectral density background level is now equal to zero in average showing that the bias has been correctly subtracted.

4 Conclusion

A method to subtract the photon noise bias from visibility data acquired with a single-mode interferometer has been presented. An analytical expression for the bias can be established under verified assumptions. Other non-analytical methods based on the fit of the average level of power spectral densities may suffer from confusion with fringe signal and from detector instabilities and the analytical method should be preferred.

Acknowledgements

I wish to thank P. Bordé for his careful reading of the paper and for his precious comments that helped me write the final version.

References

 
• Coudé du Foresto, V., Ridgway, S. T., & Mariotti, J.-M. 1997, A&AS, 121, 379 In the text NASA ADS  
• Goodman, J. 1985, Statistical Optics (J. Wiley & Sons) , 43 In the text  
• Perrin, G. 2002, submitted In the text  
• Tango, W. J., & Twiss, R. Q. 1980, Progress in Optics 17, ed. E. Wolf, 239 In the text  
• Thiébaut, É. 1994, Ph.D. Thesis, Université Paris VII dissertation In the text