10 February 2012 Effect of noise on modulation amplitude and phase in frequency-domain diffusive imaging
Author Affiliations +
J. of Biomedical Optics, 17(1), 016010 (2012). doi:10.1117/1.JBO.17.1.016010
We theoretically investigate the effect of noise on frequency-domain heterodyne and/or homodyne measurements of intensity-modulated beams propagating through diffusive media, such as a photon density wave. We assumed that the attenuated amplitude and delayed phase are estimated by taking the Fourier transform of the noisy, modulated output data. We show that the estimated amplitude and phase are biased when the number of output photons is small. We also show that the use of image intensifiers for photon amplification in heterodyne or homodyne measurements increases the amount of biases. Especially, it turns out that the biased estimation is independent of AC-dependent noise in sinusoidal heterodyne or homodyne outputs. Finally, the developed theory indicates that the previously known variance model of modulation amplitude and phase is not valid in low light situations. Monte-Carlo simulations with varied numbers of input photons verify our theoretical trends of the bias.
Kang and Kupinski: Effect of noise on modulation amplitude and phase in frequency-domain diffusive imaging



Medical imaging in the regime of visual wavelengths, such as diffuse optical tomography (DOT) or fluorescence lifetime imaging (FLI), has emerged as a non-invasive, safe, and cost-effective method to investigate physiologic and/or biological quantities in tissues.1,2 Frequency- or time-domain measurements are widely used data-acquisition methods to improve the system performance in diffusive imaging. These methods measure the frequency-response of diffusive media in addition to the steady-state response (DC). The additional information measured from frequency- or time-domain measurements generally alleviates disadvantages inherent in DC-based diffusive imaging.1,3 For example, the crosstalk between scattering (μs) and absorption (μa) coefficients of reconstructed functional images of tissues in DOT can be dramatically reduced by employing frequency-domain measurement.4,5 Also, the problem of ill-posedness and non-unique solutions typically happens in reconstruction algorithms of tomographic diffusive imaging, which can be effectively managed using frequency- or time-domain measurements.1,6

In addition to acquiring various types of data, such as a frequency-response of tissue, collecting large datasets also enhances system performance in both optical tomography and topography imaging.2 Researchers have considered collecting large datasets using imaging systems with detector arrays, such as charge-coupled device (CCD) or complementary metal-oxide-semiconductor (CMOS) detectors, in diffusive imaging as opposed to conventional fiber-based systems that collect data at relatively a few measurement points.7 It is demonstrated that DC data measured by a CCD can relieve the ill-posedness and improve the resolution of reconstructed images in DOT.8 As another example, a CCD-based DOT system operating at up to 1 GHz frequency of intensity modulation was built for investigating small tissue volume.9 It is also reported that the reconstructed images from fluorescent optical tomography was improved by large DC datasets.10 Recently, the Hotelling observer performance for frequency-domain phased array systems with a detector array is investigated showing that stable and high detectability is achieved at very low modulation frequencies.11

For frequency-domain systems, the intensity of an incident beam is modulated at, typically, hundreds of MHz. The modulation amplitude and phase of an output modulation beam from a tissue are attenuated and delayed due to scattering and absorption characteristics of the tissue. Since the typical modulation frequency in diffusive imaging is usually higher than the achievable sampling rates of current CCD or CMOS detectors, a gain-modulated image intensifier is used to make the modulation amplitude and phase of the output measurable. There are two methods of the gain-modulation measurement: the heterodyne and homodyne methods.3 The heterodyne method makes the output modulate at a beat frequency, that is the difference between source and gain-modulation frequencies. The beat frequency is selected from the purpose of diffusive imaging considering the maximum frame rate of the CCD or CMOS detectors. The homodyne method is conceptually the same as the heterodyne in that a beat frequency output is detected, but the output beam intensity is varied along the beat frequency trajectory by changing the relative phase between gain and source modulations. For this case, the output images are stable in time, so it is easier to average multiple images to reduce measurement noise. The measurement speed of the homodyne method is, however, limited due to instrumentation required to change the modulation phase, which indicates that homodyne measurements might not be appropriate for very fast dynamic diffusive imaging, where heterodyne measurements operating at up to a few kHz beat frequency have the potential to monitor the dynamic variations of biological quantities.12,13

Time-domain systems are conceptually similar to frequency-domain systems. The impulse response function is attenuated and distorted by scattering and absorption in a tissue, which contains the frequency-response information of attenuated amplitude and delayed phase. The time-domain approach has the benefit that one can estimate the frequency-response of tissue for a wide range of frequencies from a single measurement. Furthermore, for very thick tissues, where it is hard in a frequency-domain measurement to extract meaningful data due to very small number of output photons, an impulse response function can be measured by the time correlated single photon counting (TCSPC) method.5,14 However, the instrumentation for time-domain systems is usually more complicated and costly than that of frequency-domain systems. More discussion for the comparison between frequency- and time-domain measurements can be found in Refs. 1 and 3.

For both frequency- and time-domain measurements, the estimated DC, amplitude, and phase are random variables due to measurement noise. The randomness of these quantities detrimentally affects tomographic reconstructions employing these quantities as well as biological parameters estimated from these quantities. Therefore, it is important to investigate the statistical characteristics of the frequency-response quantities with an appropriate noise model. V. Toronov et al.15 theoretically derived variances of modulation amplitude and phase in a frequency-domain DOT by Fourier transforming heterodyne outputs with quantum shot noise. These analytic expressions of amplitude and phase variances have been used to investigate optimized modulation frequencies for small volume tissues and to evaluate a CCD-based phased array system in DOT.11,16 In another study, the influence of system parameters on amplitude and phase signal-to-noise ratios (SNRs) was analyzed.17 The statistical characteristics of the frequency-response quantities also play an important role in FLI, where the fluorescent lifetime of fluorophores, an indicator of biological characteristics of tissues, is estimated. Instead of directly measuring the lifetime, which is typically a few or tens of nanoseconds, the frequency-response of fluorophores is measured using intensify-modulated beams. In FLI, the effect of Poisson noise might be significant because the number of fluorescent output photons is occasionally very small (several hundreds).18,19 E. Gratton et al.20 simulated the disagreement of modulation amplitudes and phases between frequency- and time-domain measurements in accordance with output photon counts in FLI. J. Philip and K. Carlsson19 theoretically studied the effect of photon noise on the fluorescent lifetime SNR estimated from various methods, such as lock-in amplifiers, frequency-domain measurements, and demodulation. A. Esposito et al.21 investigated photon economy and acquisition speed in various FLI detection techniques using Monte-Carlo simulations, where they observed that when the number of detected photons is small, the algorithm for estimating the florescent lifetime converges to incorrect values. Recently, Y. Lin and A. Gmitro22 theoretically found out that fluorescent lifetimes estimated by a lock-in amplifier method are biased in the situation of the low number of detected fluorescent photons because of Poisson photon noise.

The discrete-Fourier transform is widely used to estimate the attenuated amplitude and delayed phase of the measured heterodyne or homodyne or impulse time response outputs, since it is simple and fast.7,9,14,1920.21 In this paper, we derive expressions for means and variances of modulation amplitude and phase that are estimated by Fourier transforming noisy sinusoidal measurements from heterodyne or homodyne methods. Other methods, such as data fitting or maximum-likelihood estimation can be used to extract the frequency-response information, but they generally require more computation time than Fourier analysis.23 When collecting large datasets in diffusive imaging, the computational advantage of Fourier methods over other fitting methods becomes enormous, because the number of detector pixels is very large. Although the kind of methods to estimate modulation amplitudes and phases from measured noisy data might affect the statistics of these quantities, the result of the developed theory in this paper provides insights about the effect of noise on a frequency-domain measurement.

The theoretical derivation of the mean and variance of a modulation amplitude and phase is described in Sec. 2. We first develop the theory by considering that sinusoidal outputs are contaminated by temporally uncorrelated noise. Then the derived noise model of heterodyne outputs is applied to the developed results. We theoretically observed that estimated modulation amplitude and phase are severely biased from the actual values, when the number of output photons is small. Furthermore, the amplitude and phase SNRs show large discrepancies from the previously known SNR model of V. Toronov’s et al.15 Sec. 3 shows computational results based on the developed theory. Monte-Carlo simulations are also presented, which verify the characteristics of biases. Extending the developed theory is briefly discussed in Sec. 4 along with conclusions.


Theory Development

By Fourier transforming noisy heterodyne/homodyne measurements, we can estimate the modulation amplitude and phase. We denote the true amplitude and phase as A and φ, and estimates of them as A^ and φ^, respectively. To describe the estimation process, we first describe the Fourier transform of the noisy data evaluated at the beat frequency Ω via


where i(t) is the noisy photocurrent rate at a detector pixel, which is linearly related to the sinusoidal photon arriving rate ap(t) [used in Eq. (3)]. T is the total measurement time which is usually a multiple of the sinusoidal output period. Increasing T implies that more data are used to estimate A^ and φ^. Actually, the photocurrent i(t) is measured at discrete, short-interval time points, and the discrete-Fourier transform is employed instead of the continuous transform shown in Eq. (1). Although not shown here, we verified that the theoretical results derived using the continuous transform are equivalent to that of the discrete transform. Thus, we discuss the theory with the continuous Fourier transform for simplicity. As mentioned in the previous section, V. Toronov et al. mathematically derived standard deviations of A^=X^R2+X^I2 and φ^=tan1(X^I/X^R) considering a shot noise model for the statistics of i(t) in Eq. (1). They found that the standard deviation of A^ and φ^ to be


where D and A indicate actual DC and AC means of the measured heterodyne output and β is a detector-related parameter which is inversely proportional to T. The detail expression for β is shown in the reference.15 The superscript t is used to differentiate standard deviations derived by Toronov et al. from those derived in this paper. The significance of Eq. (2) is that the variances of modulation amplitude and phase can be given by means. However, only a few series terms are considered for the derivation of Eq. (2) assuming A^=A and φ^=φ, which are not generally true, as shown in this paper.

Assuming temporally uncorrelated noise for i(t) in Eq. (1), the mean and covariance for i(t) can be written as




In Eq. (3), q is the electron charge and βm and βc are coefficients related to the conversion of incident photons to output (i.e., current or voltage) in a detector, where the units of βm and βc are the same.24 The terms ap(t) and σp2(t) in Eq. (4) indicate the mean and variance of the noisy sinusoidal photon rate, respectively, which oscillate at Δω/2π, as shown in Eq. (4). For the heterodyne measurement, Δω/2π is a heterodyne beat frequency and φ indicates the shift of the sinusoidal output from the origin. It is noted that σp2(t) in Eq. (4) becomes constant in time, like the case of thermal noise, when σAC2=0. The terms aDC, aAC, σDC2, and σAC2 in Eq. (4) will be specified later from the known noise model for heterodyne outputs. For mathematical simplicity, exponential functions are considered in the mean and variance of Eq. (4) instead of the co-sinusoidal functions; the final result is unaffected by this consideration.

For the statistics of A^ and φ^, the statistics of X^ considering Eqs. (3) and (4) should be investigated first. The mean of X^ is derived as


where Π(x) and sinc(x) are rect and sinc functions, respectively. For a period function, such as Eq. (4), usually T=n/Ω, where Ω=Ωo=Δω/2π and n is an integer. Resulting from this, Eq. (5) can be simplified, where means of X^R and X^I are


respectively. In order to derive the variances of X^R and X^I, second moments, X^(Ω)X^*(Ω) and X^2(Ω) should be calculated. The mathematical procedure for X^(Ω)X^*(Ω) with Eqs. (1), (3), and (4) is


where the conditions for T and Ω used for deriving means of Eq. (6) are applied to the last step of Eq. (7). Similarly, the result of X^2(Ωo) can be derived as


From Eqs. (7) and (8), the second moments and correlations of X^R and X^I are




respectively. When trigonometric identities are applied to Eqs. (6), (9), and (10), the variances (σR2 and σI2) and cross covariance of X^R and X^I can be derived as


respectively. It is noted that all derived results are for Ωo=Δω/2π, because the heterodyne output i(t) modulates at Δω/2π indicating statistical characteristics of A^ and φ^ at Ωo are our only concerns.

It is reasonable to consider X^ in Eq. (1) is complex multivariate Gaussian from the central limit theorem, because all points of the noisy signal i(t) are weighted and summed to form X^. Considering Eqs. (6) and (11), the random variable transforming25 from (X^R,X^I) of the multivariate Gaussian to (A^,φ^) generates probability density functions (PDF) for A^ and φ^ as




respectively, where


and In(x) and erf(x) indicate a n’th order modified Bessel function of the first kind and an error function, respectively. It is noted that A=XR2+XI2=qβmTaAC is the actual mean of the amplitude of i(t) at Ω0. Separable PDFs of A^ and φ^, as shown in Eqs. (12) and (13) indicate that these random variables are independent each other. Equations (12(13)14) imply that γ plays a critical role for the PDFs of A^ and φ^, which will be investigated in simulation in the next section. Interestingly, γ is not affected by σAC2 and aDC in Eq. (4). Expressions similar to Eqs. (12) and (13) were previously reported by S. Morgan.26 However, he did not specify the detailed form and physical meaning of γ because he derived them from a simple Gaussian noise model.

The mean and variance of A^ can be analytically calculated from Eq. (12), which are




respectively. The amplitude mean and variance of Gaussian speckle are equivalently expressed as Eqs. (15) and (16), although physical meanings of speckle are different from the frequency-domain measurement.25 Equations (11) and (16) show that σA2 becomes the previously known Toronov model of Eq. (2), only when A¯=A, which is not generally acceptable, as indicated from Eq. (15). Equation (15) shows that A¯, the mean of A^ estimated from noisy sinusoidal outputs is biased and the amount of bias depends on σ2 and γ in Eqs. (11) and (14). There are no analytic solutions for the moments of φ^, but they can be numerically calculated from Eq. (13). As shown in the next section, φ¯ calculated from Eq. (13) is also different from φ and the bias depends on σ2 and γ.

With the similar mathematical procedure for calculating moments of X^R(Ω0) and X^I(Ω0), the mean and variance for estimated DC D^ can be derived from that Ω=0 is substituted to Eqs. (5), (7), and (8), instead of considering Ω=Ωo. The condition of Ω=0 generates X^R(0)=qβmTaDC, X^I(0)=0, σR2=q2βc2TσDC2, and σI2=0, which denotes the PDF of D^ is univariate Gaussian. Notice that the variance of D^, q2βc2TσDC2 is the same as σA2 in Eq. (16) only if A¯=A. We will observe in simulation that A¯ is generally greater than A, unless the number of measured photons is very large, so it can be stated that the amplitude variance is generally smaller than the DC variance in frequency-domain measurements. Since D^ is univariate Gaussian, it is not difficult to show that D¯ is unbiased.

Equations (12)(13)(14)(15)–(16) are generally applied for estimating the amplitude and phase of any sinusoidal signal having the noise property of Eq. (3). Recently, the covariance of a heterodyne and/or homodyne output was theoretically derived from the random amplification of a temporal point process.27 The model of covariance implies that the variance of heterodyne outputs can be expressed as the form of σp2(t) in Eq. (4) and both σDC2 and σAC2 are affected by a Fano factor F24 of image intensifiers, which is experimentally verified. A Fano factor is defined by the variance divided by the mean of amplified (secondary) photons per input (primary) photon in image intensifiers. If the gain of an image intensifier in the heterodyne measurement is assumed to be high, and the spatial variation of the amount of primary photons incident on the image intensifier is not abruptly changed, the covariance model specifies aAC and σDC2 in Eq. (14) as27




respectively. Terms cp and mp are the DC photon rate and the modulation depth of modulating primary photons incident onto an image intensifier. The modulated gain in the image intensifier is generally non-sinusoidal, which can be expanded to Fourier cosine terms of amplitude c0mg.28 The subscript g in Eqs. (17) and (18) is an integer and c0mg=1 indicates the amplitude for the main frequency of the modulated gain. The term k¯nI is the gain achieved on the micro channel plate (MCP) of image intensifiers, which separately functions from the DC cathode gain c0. Therefore, c0k¯nI indicates the overall image intensifier’s gain. Substituting aAC and σDC2 of Eqs. (17) and (18) to Eq. (14) produces γ for the case of heterodyne measurements, which is


The term Tcp in Eq. (19) indicates the number of primary photons incident into the image intensifier during T, which determines the unit of γ; it is unitless.


Simulation Results and Discussion

Figures 1(a) and 1(b) show some examples of prA^(A^) and prφ^(φ^) for arbitrarily chosen values of γ, A, and φ, respectively, to investigate characteristics of these PDFs. The values of βm and βc in Eq. (3) are also arbitrarily chosen as 0.5 and 0.6, respectively, which will be maintained for all simulations in this paper. Having βm, βc<1 indicates the conversion loss from photons to photoelectrons without an electronic gain in the CCD or CMOS detectors. As expected from the definition of PDF, all PDFs plotted in Fig. 1 are non-negative functions. Although all PDFs appear Gaussian shaped, they are actually not, as mathematically verified from Eqs. (12) and (13). Figure 1 indicates that the asymmetry in both prA^(A^) and prφ^(φ^) increases as γ decreases, which results in biased means. From characteristics of prA^(A^) shown in Fig. 1(a), it is expected that A¯ is higher than A for the smaller value of γ. Oppositely, the behavior of prφ^(φ^) in Fig. 1(b) shows that φ¯ is smaller than φ because the value of prφ^(φ^) is increased around π as γ decreases. From other simulation results, we confirmed that when φ is negative, prφ^(φ^) around π starts increasing as γ decreases. Therefore, it can be stated that φ¯ is always under-estimated regardless of its actual mean. Similarly, we observed that A¯ is always over-estimated. It can be also conjectured from Fig. 1 that the amount of bias for both A^ and φ^ is dependent on their actual mean values. Furthermore, it can be expected that smaller γ increases the bias of A¯ and φ¯ as well as the standard deviations of A^ and φ^.

Fig. 1

Some examples of PDFs of the modulation amplitude and phase are plotted with arbitrary parameters. Notice that these PDFs deviate from the symmetry as γ is decreased.


Since it is definite that γ is critical for statistics of A^ and φ^, as observed in Fig. 1, it is worth to briefly discuss the effect of heterodyne measurement parameters on γ in Eqs. (14) and (19). It is clear from Fig. 1 that higher γ is generally better to correctly estimate A^ and φ^, which is achieved by the larger amount of primary photon, Tcp in Eq. (19). If modulating primary photons are assumed to be directly detected without heterodyne and/or homodyne processes, it is reasonable to consider that the main noise for the primary is Poisson. For Poisson photon noise, σDC2=cp and aAC=cpmp in Eq. (14), so γ is


It is straightforward that γ in Eq. (20) is always larger than that in Eq. (19). This indicates that even though the image intensifier in heterodyne processes amplifies primaries to increase the number of detected secondary photons, the estimation for A^ and φ^ is worse because the Poisson noise is also amplified. For Eq. (19), γ can be increased as mp and/or mg=1 increase. As the gain-modulation deviates from a perfect sinusoidal function, g=1mg2 in the denominator is more rapidly increased than mg=12 in the numerator, hence reducing γ. The term F/c0k¯nI in Eq. (19) is equivalent to the Noise Factor NF, when primary photons are assumed as Poisson.27 Therefore, NF, commonly used in noise analyses for photon amplifiers, affects γ profoundly. It is known that NF can be minimized with a very high image intensifier gain, which is 3.5 to 4.2 and 1.6 to 2.2 for second- and third-generation image intensifiers, respectively.29 Since NF is increased as the gain decreases, γ becomes small at low-gain situations.

We assumed an image intensifier’s gain is large enough to validate Eqs. (17) and (18), so NF=4 is selected for simulations. Figure 2 shows percent differences between actual and estimated means for A^ and φ^, respectively, with parameters of mg=1=1.3, and g=1mg2=1.5×m12 for γ in Eq. (19). A value of 1.3 for mg=1 is considered from the experimental result27 and g=1mg2=1.5×m12 indicates the summation of mg12 is assumed to be 50% of m12. Considering the typical quantum efficiency of a photocathode in image intensifiers, c0 in γ is set to 0.3. The larger value among estimated and actual means is used as a denominator to set the percent difference within ±100%. Since T is usually much longer than the modulation period of the primary, T·cp in x-axes of Fig. 2 can be considered as the average primary photons during T. For frequency-domain diffusive imaging, mp is largely varied across a phantom exit surface. Thus, A¯ and φ¯ in Fig. 2 are simulated for three different mp values. As shown in Fig. 2, both A¯ and φ¯ are severely biased when the number of primary photons from a phantom is small, where A¯ and φ¯ show over- and under-estimations, respectively, as expected from the results in Fig. 1. It is verified from other additional simulations that φ¯ always approaches zero as the primary photon number decreases. Furthermore, Fig. 2(b) indicates the amount of bias for φ¯ is strongly dependent on φ, which is consistent with the behavior of prφ^(φ^) in Fig. 1(b). It should be noticed that the typical range of c0k¯nI is from hundreds up to ten thousands, so the number of secondary photons for the primaries in Fig. 2 is actually very large. Figure 2 implies that for the diffusive imaging, where the number of output photons from a phantom is small, tomographic functional images or biological quantities estimated from A^ and φ^ may be highly inaccurate.

Fig. 2

Percent differences between actual and estimated means for (a) A^ and (b) φ^ are plotted as a function of primary photon number of the output beam from a diffusive medium. The terms mp and φ indicate an actual modulation depth and phase of the output beam, respectively.


Figure 3 shows percent differences between SNR calculated from the derived theory in this paper and SNR (SNRt) from the previously known model of Eq. (2) for the same parameters as Fig. 2. For SNRt, A¯ and φ¯ are assumed to be correct, A and φ, respectively. Like Fig. 2, the larger one between SNR and SNRt is used as a denominator for calculating these percent differences. There is a good agreement between SNR and SNRt when the primary photon number from a phantom is relatively large. However, the agreement starts deviating as primary photon number is decreased, where SNRA and SNRφ are larger and smaller than SNRAt and SNRφt, respectively. For SNRφ, especially, the primary photon number starting the discrepancy and the amount of the percent difference are magnificently increased, compared with the bias of φ¯ in Fig. 2. Figure 3 shows that it is not appropriate to estimate system performance based on the previously known noise model for frequency-domain diffusive imaging systems with the small output photon numbers, such as the system observing dynamic variation of biological quantities or FLI.

Fig. 3

Percent differences between SNR calculated from the theory and SNR from the previously known model for (a) A^ and (b) φ^ are plotted as a function of primary photon number of the output beam from a diffusive medium. mp and φ indicate an actual modulation depth and phase of the output beam, respectively.


We examine the estimated means and their characteristics using brute-force and Monte-Carlo simulations. For brute-force simulations, we generate sinusoidal signals of one period with known D and A, where white Gaussian noise is added. Fourier transforming the noise-contaminated signals determines D^ and A^ that are random. Averaging these random variables by repeating brute-force simulations can estimate D¯ and A¯. Figure 4(a) shows one example of the noisy sinusoidal signal that mimics a heterodyne or homodyne output. The standard deviation of the noisy data is 10 and the solid line indicates the mean sinusoidal signal. D¯ and A¯ estimated from the ensemble of the noisy signals are shown in Fig. 4(b) as increasing the standard deviation from 0 to 30. As theoretically investigated, D¯ is unbiased for the entire range of the noise level. However, A¯ shows the bias estimation, the amount of which is increased as the noise level increases—over-estimated. In order to investigate biased estimations in the situation similar to diffusive imaging, Monte-Carlo simulations are conducted using tMCimg30 for a 62×62×30mm3 homogeneous phantom of μs=5mm1 and μa=0.005mm1. The tMCimg simulator is widely used in diffusive imaging Monte-Carlo simulations, where the Henyey-Greenstein phase function is assumed for scattering. The anisotropy parameter for the Henyey-Greenstein phase function is set to 0.8 in this simulation. The source of a 2 mm diameter injects photons to the phantom, of which number is varied from 0.5×108 to 4×108, as shown in Fig. 4. The degree of noise on output photons are decreased as the total photon number increases, because the increased photon number is achieved by summing independently simulated outputs. For example, the output from 4×108 of total photon number is the summation of eight outputs of 0.5×108 source photons. The 40×40mm2 detector array with 1×1mm2 pixels on the exit side of the phantom collects simulated output photons. Different from brute-force simulations in Figs. 4(a) and (b), temporal functions of output photons on each detector pixel are Fourier transformed to calculate D^, A^, and φ^. Since it is obvious that D¯ is unbiased from the developed theory and brute-force simulation, the term m¯p is investigated instead of A¯ to minimize the numerical error of the Monte-Carlo simulations, such as photon number mismatching. Figures 4(c) and 4(d) show the averaged modulation depth (m¯p=A¯/D¯) and φ¯ from the entire detector pixels at fixed modulation frequencies (fm). Averaging from the entire pixels avoids a heavy computational load because enormous number of simulations is required to calculate means for one detector pixel. Although actual means mp and φ are different for different detector pixels, PDFs of m^p and φ^ are the same (with different means), as shown in Eqs. (12) and (13), respectively. Therefore, biased estimations would be clearly observed for m¯p and φ¯ averaged from all detector pixels if the theoretical results derived in Sec. 2 are valid. For the comparison, all curves calculated at a fixed fm are scaled for both m¯p and φ¯ to make the values of the photon number of 0.5×108 are the same. Figures 4(c) and (d) show that m¯p and φ¯ are decreased and increased, respectively, as the number of simulated photons increases, which indicates that m¯p and φ¯ are over- and under-estimated. Furthermore, it is observed that the degree of bias is worse for higher fm, where the number of output photons contributing to the calculation for m¯p and φ¯ is smaller. It is well known that m¯p and φ¯ are decreased and increased, respectively, as fm increases, which shows the validity of the Monte-Carlo simulations. Although other systematic factors in heterodyne or homodyne measurements, such as a photon amplifier, signal and noise amplification, and photon losses in these brute-force and Monte-Carlo simulations, Fig. 4 clearly shows that characteristics of estimated quantities from simulations are coincident with what expected from the developed theory.

Fig. 4

Estimated mean DC and AC, D¯ and A¯ are calculated from sinusoidal signals with white Gaussian noise of increased standard deviations. It is clear that D¯ and A¯ are unbiased and over-estimated, respectively, in (b). In (c) and (d), Monte-Carlo simulated means are shown as a function of simulated photons, where the relative amount of noise is decreased as the number of incident photons is increased. Estimated means of (c) mp and (d) φ^ from simulations are over- and under-estimated, respectively.


As shown from the theory and simulation, increasing the output photon number from a phantom can definitely remedy the bias estimation and make the simplified noise model of Eq. (2) acceptable. This is definitely achieved by increasing the total measurement time T which indicates that quite a number of heterodyne and/or homodyne output periods are acquired to estimate A^ and φ^. However, increasing measurement time is not allowed or difficult in some situations, where, for example, the imaging purpose is to observe the dynamic variation of biological quantities, such as hemodynamic changes.31 Alternatively, increasing the DC primary photon rate cp decreases the amount of bias estimations with a fixed T. Because the photoelectron capacity of CCD or CMOS detector pixels is usually limited, however, the gain of image intensifier should be decreased to avoid the saturation of detector outputs for this increased cp. The lower gain indicates the noise factor NF of the image intensifier should be increased. Resulting from this, γ in Eq. (19) might not be increased enough to diminish the bias even for the increased amount of primary photons.

Percent differences for A¯ and φ¯ with increased NFs are shown in Figs. 5(a) and 5(b), respectively, where NF ranges 4 to 10 and 4 to 16. It was experimentally shown that these NF ranges are easily achieved in common image intensifiers of decreased gain.29 For comparison, results of NF=4 in Fig. 2 are repeatedly presented in Fig. 5. We assume NF linearly increases for decreasing the MCP gain k¯nI that makes the number of secondary photons constant for the increased primary photons in Fig. 5. This implies σ2_DC of Eq. (18) is increased by the increased NF. Other simulation parameters in Fig. 5 are the same as Fig. 2. As shown in Fig. 5, the degree of bias is worse than the results in Fig. 2, which is caused by increased NF. Interestingly, the detrimental effect of increased NF is much more significant for A¯ than φ¯. It is observed in Fig. 5(b) that the effect of increased NF is almost ignorable for cases of φ smaller than 3π/4. Higher NF indicates smaller SNR of secondary photons for given SNR of primary photons, which physically means the secondary photons are more random. Figure 5 indicates that estimating modulation amplitude is much more sensitive to the randomness of output photons than phase.

Fig. 5

Percent differences between actual and estimated means for (a) amplitudes and (b) phases are plotted. For the increasing primary photon number, NF is linearly increased with ranges of 410 and 416.



Further Discussion and Conclusion

We observed that modulation amplitude and phase A^ and φ^ are biased when they are estimated by Fourier transforming the noisy heterodyne outputs. Furthermore, the previously known SNR model for A^ and φ^ is not valid when the number of output photons from a phantom is small. It is known that a very high modulation frequency is required to achieve high SNRs for φ^ in diffusive imaging for small volume phantoms9,16 because a high modulation frequency increases φ in Eq. (4). However, the theoretical result in this paper shows that high φ (i.e., the large amount of phase delay of the modulated incident beam to a phantom) deteriorates the reliability of the estimated phase φ^. Furthermore, the modulation depth mp in Eqs. (17) and (19) is dramatically decreased as the modulation frequency increases, which further degrades the estimation for φ^. Although recently developed noise models for heterodyne and homodyne measurements indicate that small mp reduces σAC2 in Eq. (4),27 the decreased σAC2 does not improve the estimation for φ^ because γ in Eq. (14) depends on only aAC and σDC2. Therefore, it can be importantly noticed that increasing SNR for φ^ does not always reduce the degree of bias estimation.

In frequency-domain diffusive imaging with a detector array, increasing detecting points usually accompanies decreasing the detector pixel size, which reduces the measured primary photons per pixel. Furthermore, for detector pixels far away from the incident source, the modulation depth mp as well as the DC photon rate cp are decreased because of increased scattering and absorbing in the diffusive medium. Therefore, the bias in A^ and φ^ is substantially increased for high-dimensional frequency-domain diffusive imaging. Additionally, the typical loss of photons in heterodyne and/or homodyne frequency-domain measurements makes the photon number necessary for correctly estimating A^ and φ^ increased. For example, in some types of systems, where the phosphor of an image intensifier is imaged onto a CCD by an imaging system, the amount of photon loss due to the imaging system is usually significant. This loss of photons is equivalent to decreasing primary photons.

Although the theory in this paper is derived from Fourier analysis of a noisy sinusoidal signal, it can be expanded for time-domain measurements. In fact, it is reported elsewhere that φ^ estimated by Fourier transforming simulated impulse outputs from TCSPC in FLI are under-estimated as compared to the exact values, where the amount of under-estimation becomes larger for a higher modulation frequency.20 Similarly, the Monte-Carlo simulation results in Fig. 4 are calculated by Fourier transforming impulse response functions of output photons at each detector pixel. As long as A^ and φ^ are estimated by Fourier transforming a temporal impulse function that is contaminated by temporally uncorrelated noise, the theoretical analysis and characteristics for the statistics of A^ and φ^ is very similar with what is observed in this paper. For example, if the impulse output ap(i)(t) is measured without amplification, σDC2 in Eq. (4) can be considered as T/2T/2ap(i)(t)/T based on that the variance of shot noise induced by a time-variant signal is described by time-averaging the signal.32.

The theoretical procedure in this paper is based on a heterodyne measurement, but the similar results can be derived for a homodyne measurement. The lock-in amplifier method is widely used for estimating fluorescent lifetimes in FLI, where two (or more) co-sinusoidal signals are multiplied to measured modulated fluorescent outputs. The relative phases between modulation outputs and multiplied signals are set to be different like in a homodyne method. Therefore, the concept of lock-in amplifier is very similar with that of heterodyne measurement except the photon amplification process in a homodyne method, which is optional in lock-in amplifiers. We expect that this similarity is the reason that fluorescent lifetimes estimated from a lock-in amplifier show bias, when the detected fluorescent photons are small.19,22 Although the theoretical approach22 to show the biased fluorescent lifetime is different from what is presented in this paper, we consider the origin of the phenomenon is conceptually similar with what we have observed in this paper.

The effect of noise on modulation amplitude and phase in frequency-domain diffusive imaging was investigated. It was theoretically derived that the Fourier analysis for a noisy sinusoidal signal induces the estimated amplitude and phase to be biased from actual values when the detected photon number is small. By combining the developed theory with the known noise model of heterodyne and/or homodyne outputs, the degree of bias was measured as a function of the number of diffused photons exiting the phantom. Furthermore, physical insight including the effect of heterodyne and/or homodyne measurement parameters on the bias was elucidated. Especially, we showed that the amplitude and phase biases are deteriorated by the increased noise factor of image intensifiers, which is common for heterodyne and homodyne measurements with low-gain image intensifiers. It is interesting that the noise depending on output signal variation (i.e., AC variance on noisy sinusoidal signals) does not affect the estimation of modulation amplitude and phase. Finally, we pointed out that the previously known noise model for the modulation amplitude and phase is not valid for the situation of small output photons. Therefore, considering the known noise model should be careful because, for many frequency-domain diffusive imaging situations, output photon counts are low. In the near future, experimentally measuring the bias for the estimated amplitude and phase will be investigated. Also, it is worth assessing the performance of frequency-domain diffusive imaging systems in photon-limited situations based on statistical properties of the estimated amplitude and phase developed in this paper.


1. A. P. GibsonJ. C. HebdenS. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005).PHMBA70031-9155 http://dx.doi.org/10.1088/0031-9155/50/4/R01 Google Scholar

2. B. W. Pogueet al., “Image analysis methods for diffuse optical tomography,” J. Biomed. Opt. 11, 033001 (2006).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2209908 Google Scholar

3. B. Chanceet al., “Phase measurement of light absorption and scatter in human tissue,” Rev. Sci. Instrum. 69, 3457–3481 (1998).RSINAK0034-6748 http://dx.doi.org/10.1063/1.1149123 Google Scholar

4. T. O. McBrideet al., “Separation of absorption and scattering heterogeneities in NIR tomographic imaging of tissue,” in Biomedical Topical Meetings, OSA Technical Digest, Optical Society of America, Washington, D. C., pp. 339–341 (2000). Google Scholar

5. G. XuD. PiaoC. F. BuntingH. Dehghani, “Direct-current-based image reconstruction versus direct-current included or excluded frequency-domain reconstruction in diffuse optical tomography,” Appl. Opt. 49, 3059–3070 (2010).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.49.003059 Google Scholar

6. S. R. ArridgeW. R. B. Lionheart, “Nonuniqueness in diffusion-based optical tomography,” Opt. Lett. 23, 882–884 (1998).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.23.000882 Google Scholar

7. A. B. ThompsonE. M. Sevick-Murace, “Near-infrared fluorescence contrast-enhanced imaging with intensified charge-coupled device homodyne detection: measurement precision and accuracy,” J. Biomed. Opt. 8, 111–120 (2003).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1528205 Google Scholar

8. Z. M. Wanget al., “Experimental demonstration of an analytic method for image reconstruction in optical diffusion tomography with large data sets,” Opt. Lett. 30, 3338–3340 (2005).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.30.003338 Google Scholar

9. U. J. NetzJ. BeuthanA. H. Hielscher, “Multipixel system for gigahertz frequency-domain optical imaging of finger joints,” Rev. Sci. Instrum. 79, 034301 (2008).RSINAK0034-6748 http://dx.doi.org/10.1063/1.2840344 Google Scholar

10. G. Y. Panasyuket al., “Fluorescent optical tomography with large data sets,” Opt. Lett. 33, 1744–1746 (2008).OPLEDP0146-9592 http://dx.doi.org/10.1364/OL.33.001744 Google Scholar

11. D. KangM. A. Kupinski, “Signal detectability in diffusive media using phased arrays in conjunction with detector arrays,” Opt. Express 19, 12261–12274 (2011).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.19.012261 Google Scholar

12. H. M. Watzmanet al., “Arterial and venous contributions to near-infrared cerebral oximetry,” Anesthesiology 93, 947–953 (2000).ANESAV0003-3022 http://dx.doi.org/10.1097/00000542-200010000-00012 Google Scholar

13. M. Atlanet al., “Cortical blood flow assessment with frequency-domain laser Doppler microscopy,” J. Biomed. Opt. 12, 024019 (2007).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.2715184 Google Scholar

14. H. E. GreccoP. Roda-NavarroP. J. Verveer, “Global analysis of time correlated single photon counting FRET-FLIM data,” Opt. Express 17, 6493–6508 (2009).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.17.006493 Google Scholar

15. V. Toronovet al., “Optimization of the signal-to-noise ratio of frequency-domain instrumentation for near-infrared spectro-imaging of the human brain,” Opt. Express 11, 2717–2729 (2003).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.11.002717 Google Scholar

16. H. K. Kimet al., “Optimal source-modulation frequencies for transport-theory-based optical tomography of small-tissue volumes,” Opt. Express 16, 18082–18101 (2008).OPEXFF1094-4087 http://dx.doi.org/10.1364/OE.16.018082 Google Scholar

17. T. Tuet al., “Analysis on performance and optimization of frequency-domain near-infrared instruments,” J. Biomed. Opt. 7, 643–649 (2002).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1501562 Google Scholar

18. M. KollnerJ. Wolfrum, “How many photons are necessary for fluorescence-lifetime measurements?,” Chem. Phys. Lett. 200, 199–204 (1992).CHPLBC0009-2614 http://dx.doi.org/10.1016/0009-2614(92)87068-Z Google Scholar

19. J. PhilipK. Carlsson, “Theoretical investigation of the signal-to-noise ratio in fluorescence lifetime imaging,” J. Opt. Soc. Am. A 20, 368–379 (2003).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.20.000368 Google Scholar

20. E. Grattonet al., “Fluorescence lifetime imaging for the two-photon microscope: time-domain and frequency-domain methods,” J. Biomed. Opt. 8, 381–390 (2003).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.1586704 Google Scholar

21. A. EspositoH. C. GerritsenF. S. Wouters, “Optimizing frequency-domain fluorescence lifetime sensing for high-throughput applications: photon economy and acquisition speed,” J. Opt. Soc. Am. A 24, 3261–3273 (2007).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.24.003261 Google Scholar

22. Y. LinA. F. Gmitro, “Statistical analysis and optimization of frequency- domain fluorescence lifetime imaging microscopy using homodyne lock-in detection,” J. Opt. Soc. Am. A 27, 1145–1155 (2010).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.27.001145 Google Scholar

23. A. ElderS. SchlachterC. F. Kaminski, “Theoretical investigation of the photon efficiency in frequency-domain fluorescence lifetime imaging microscopy,” J. Opt. Soc. Am. A 25, 452–462 (2008).JOAOD60740-3232 http://dx.doi.org/10.1364/JOSAA.25.000452 Google Scholar

24. H. H. BarrettK. J. Myers, Foundations of Image Science, John Wiley & Sons, Inc., Hoboken, New Jersey (2004). Google Scholar

25. J. W. Goodman, Chap. 2 in Speckle Phenomena in Optics, Theory and Applications, Roberts & Company, Englewood, Colorado (2006). Google Scholar

26. S. P. Morgan, “Detection performance of a diffusive wave phased array,” Appl. Opt. 43, 2071–2078 (2004).APOPAI0003-6935 http://dx.doi.org/10.1364/AO.43.002071 Google Scholar

27. D. KangM. A. Kupinski, “Noise characteristics of heterodyne/homodyne frequency-domain measurements,” J. Biomed. Opt. 17(1), 015002 (2012).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.JBO.17.1.015002 Google Scholar

28. B. Q. SpringR. M. Clegg, “Image analysis for denoising full-field frequency-domain fluorescence lifetime images,” J. Microscopy 235, 221–237 (2009).JMICAR0022-2720 http://dx.doi.org/10.1111/jmi.2009.235.issue-2 Google Scholar

29. S. E. Moranet al., “Intensified CCD (ICCD) dynamic range and noise performance,” Proc. SPIE 3173, 430–457 (1997).PSISDG0277-786X http://dx.doi.org/10.1117/12.294535 Google Scholar

30. D. A. Boaset al., “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10, 159–170 (2002).OPEXFF1094-4087 Google Scholar

31. T. Durduranet al., “Optical measurement of cerebral hemodynamics and oxygen metabolism in neonates with congenital heart defects,” J. Biomed. Opt. 15, 037004 (2010).JBOPFO1083-3668 http://dx.doi.org/10.1117/1.3425884 Google Scholar

32. P. J. Winzer, “Shot-noise formula for time-varying photon rates: a general derivation,” J. Opt. Soc. Am. B 14, 2424–2428 (1997).JOBPDE0740-3224 http://dx.doi.org/10.1364/JOSAB.14.002424 Google Scholar

© 2012 Society of Photo-Optical Instrumentation Engineers (SPIE)
Dongyel Kang, Matthew A. Kupinski, "Effect of noise on modulation amplitude and phase in frequency-domain diffusive imaging," Journal of Biomedical Optics 17(1), 016010 (10 February 2012). https://doi.org/10.1117/1.JBO.17.1.016010


Phase shift keying


Signal to noise ratio

Monte Carlo methods

Homodyne detection

Image intensifiers

Back to Top