1.

Introduction

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 (${\mu }_s$) and absorption (${\mu }_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.² 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.⁷ It is demonstrated that DC data measured by a CCD can relieve the ill-posedness and improve the resolution of reconstructed images in DOT.⁸ 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.⁹ It is also reported that the reconstructed images from fluorescent optical tomography was improved by large DC datasets.¹⁰ 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.¹¹

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.³ 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.¹⁵ 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.¹⁷ 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 ($\sim$several hundreds).^18,19 E. Gratton et al.²⁰ 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. Carlsson¹⁹ 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.²¹ 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. Gmitro²² 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,19–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.²³ 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.¹⁵ 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.

2.

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 $\phi$, and estimates of them as $\widehat{A}$ and $\widehat{\phi }$, respectively. To describe the estimation process, we first describe the Fourier transform of the noisy data evaluated at the beat frequency $\mathrm{\Omega }$ via

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

For the statistics of $\widehat{A}$ and $\widehat{\phi }$, the statistics of $\widehat{X}$ considering Eqs. (3) and (4) should be investigated first. The mean of $\widehat{X}$ is derived as

It is reasonable to consider $\widehat{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 $\widehat{X}$. Considering Eqs. (6) and (11), the random variable transforming²⁵ from $({\widehat{X}}_R,{\widehat{X}}_I)$ of the multivariate Gaussian to $(\widehat{A},\widehat{\phi })$ generates probability density functions (PDF) for $\widehat{A}$ and $\widehat{\phi }$ as

The mean and variance of $\widehat{A}$ can be analytically calculated from Eq. (12), which are

With the similar mathematical procedure for calculating moments of ${\widehat{X}}_R({\mathrm{\Omega }}_0)$ and ${\widehat{X}}_I({\mathrm{\Omega }}_0)$, the mean and variance for estimated DC $\widehat{D}$ can be derived from that $\mathrm{\Omega }=0$ is substituted to Eqs. (5), (7), and (8), instead of considering $\mathrm{\Omega }={\mathrm{\Omega }}_o$. The condition of $\mathrm{\Omega }=0$ generates $\langle {\widehat{X}}_R(0)\rangle =q{\beta }_{m}T{a}_{\mathrm{DC}}$, $\langle {\widehat{X}}_I(0)\rangle =0$, ${\sigma }_R^2=q^2{\beta }_c^{2}T{\sigma }_{\mathrm{DC}}^2$, and ${\sigma }_I^2=0$, which denotes the PDF of $\widehat{D}$ is univariate Gaussian. Notice that the variance of $\widehat{D}$, $q^2{\beta }_c^{2}T{\sigma }_{\mathrm{DC}}^2$ is the same as ${\sigma }_A^2$ in Eq. (16) only if $\overline{A}=A$. We will observe in simulation that $\overline{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 $\widehat{D}$ is univariate Gaussian, it is not difficult to show that $\overline{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.²⁷ The model of covariance implies that the variance of heterodyne outputs can be expressed as the form of ${\sigma }_p^2(t)$ in Eq. (4) and both ${\sigma }_{\mathrm{DC}}^2$ and ${\sigma }_{\mathrm{AC}}^2$ are affected by a Fano factor $F$²⁴ 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 $a_{\mathrm{AC}}$ and ${\sigma }_{\mathrm{DC}}^2$ in Eq. (14) as²⁷

3.

Simulation Results and Discussion

Figures 1(a) and 1(b) show some examples of ${\mathrm{pr}}_{\widehat{A}}(\widehat{A})$ and ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ for arbitrarily chosen values of $\gamma$, $A$, and $\phi$, respectively, to investigate characteristics of these PDFs. The values of ${\beta }_m$ and ${\beta }_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 ${\beta }_m$, ${\beta }_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 ${\mathrm{pr}}_{\widehat{A}}(\widehat{A})$ and ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ increases as $\gamma$ decreases, which results in biased means. From characteristics of ${\mathrm{pr}}_{\widehat{A}}(\widehat{A})$ shown in Fig. 1(a), it is expected that $\overline{A}$ is higher than $A$ for the smaller value of $\gamma$. Oppositely, the behavior of ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ in Fig. 1(b) shows that $\overline{\phi }$ is smaller than $\phi$ because the value of ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ is increased around $-\pi$ as $\gamma$ decreases. From other simulation results, we confirmed that when $\phi$ is negative, ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ around $\pi$ starts increasing as $\gamma$ decreases. Therefore, it can be stated that $\overline{\phi }$ is always under-estimated regardless of its actual mean. Similarly, we observed that $\overline{A}$ is always over-estimated. It can be also conjectured from Fig. 1 that the amount of bias for both $\widehat{A}$ and $\widehat{\phi }$ is dependent on their actual mean values. Furthermore, it can be expected that smaller $\gamma$ increases the bias of $\overline{A}$ and $\overline{\phi }$ as well as the standard deviations of $\widehat{A}$ and $\widehat{\phi }$.

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 $\gamma$ is decreased.

Since it is definite that $\gamma$ is critical for statistics of $\widehat{A}$ and $\widehat{\phi }$, as observed in Fig. 1, it is worth to briefly discuss the effect of heterodyne measurement parameters on $\gamma$ in Eqs. (14) and (19). It is clear from Fig. 1 that higher $\gamma$ is generally better to correctly estimate $\widehat{A}$ and $\widehat{\phi }$, which is achieved by the larger amount of primary photon, $T{c}_p$ 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, ${\sigma }_{\mathrm{DC}}^2=c_p$ and $a_{\mathrm{AC}}=c_p{m}_p$ in Eq. (14), so $\gamma$ is

We assumed an image intensifier’s gain is large enough to validate Eqs. (17) and (18), so $\mathrm{NF}=4$ is selected for simulations. Figure 2 shows percent differences between actual and estimated means for $\widehat{A}$ and $\widehat{\phi }$, respectively, with parameters of $m_{g=1}=1.3$, and $\sum _{g=1}^{\infty }{m}_g^2=1.5\times m_1^2$ for $\gamma$ in Eq. (19). A value of 1.3 for $m_{g=1}$ is considered from the experimental result²⁷ and $\sum _{g=1}^{\infty }{m}_g^2=1.5\times m_1^2$ indicates the summation of $m_{g\ne 1}^2$ is assumed to be 50% of $m_1^2$. Considering the typical quantum efficiency of a photocathode in image intensifiers, $c_0$ in $\gamma$ is set to 0.3. The larger value among estimated and actual means is used as a denominator to set the percent difference within $\pm 100\%$. Since $T$ is usually much longer than the modulation period of the primary, $T·c_p$ in $x$-axes of Fig. 2 can be considered as the average primary photons during $T$. For frequency-domain diffusive imaging, $m_p$ is largely varied across a phantom exit surface. Thus, $\overline{A}$ and $\overline{\phi }$ in Fig. 2 are simulated for three different $m_p$ values. As shown in Fig. 2, both $\overline{A}$ and $\overline{\phi }$ are severely biased when the number of primary photons from a phantom is small, where $\overline{A}$ and $\overline{\phi }$ show over- and under-estimations, respectively, as expected from the results in Fig. 1. It is verified from other additional simulations that $\overline{\phi }$ always approaches zero as the primary photon number decreases. Furthermore, Fig. 2(b) indicates the amount of bias for $\overline{\phi }$ is strongly dependent on $\phi$, which is consistent with the behavior of ${\mathrm{pr}}_{\widehat{\phi }}(\widehat{\phi })$ in Fig. 1(b). It should be noticed that the typical range of $c_0{\overline{k}}_n^I$ 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 $\widehat{A}$ and $\widehat{\phi }$ may be highly inaccurate.

Fig. 2

Percent differences between actual and estimated means for (a) $\hat{A}$ and (b) $\hat{\phi }$ are plotted as a function of primary photon number of the output beam from a diffusive medium. The terms $m_p$ and $\phi$ 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 (${\mathrm{SNR}}^t$) from the previously known model of Eq. (2) for the same parameters as Fig. 2. For ${\mathrm{SNR}}^t$, $\overline{A}$ and $\overline{\phi }$ are assumed to be correct, $A$ and $\phi$, respectively. Like Fig. 2, the larger one between SNR and ${\mathrm{SNR}}^t$ is used as a denominator for calculating these percent differences. There is a good agreement between SNR and ${\mathrm{SNR}}^t$ when the primary photon number from a phantom is relatively large. However, the agreement starts deviating as primary photon number is decreased, where ${\mathrm{SNR}}_A$ and ${\mathrm{SNR}}_{\phi }$ are larger and smaller than ${\mathrm{SNR}}_A^t$ and ${\mathrm{SNR}}_{\phi }^t$, respectively. For ${\mathrm{SNR}}_{\phi }$, especially, the primary photon number starting the discrepancy and the amount of the percent difference are magnificently increased, compared with the bias of $\overline{\phi }$ 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) $\hat{A}$ and (b) $\hat{\phi }$ are plotted as a function of primary photon number of the output beam from a diffusive medium. $m_p$ and $\phi$ 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 $\widehat{D}$ and $\widehat{A}$ that are random. Averaging these random variables by repeating brute-force simulations can estimate $\overline{D}$ and $\overline{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. $\overline{D}$ and $\overline{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, $\overline{D}$ is unbiased for the entire range of the noise level. However, $\overline{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 tMCimg³⁰ for a $62\times 62\times 30{\mathrm{mm}}^3$ homogeneous phantom of ${\mu }_s=5{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.005{\mathrm{mm}}^{-1}$. 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\times {10}^8$ to $4\times {10}^8$, 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\times {10}^8$ of total photon number is the summation of eight outputs of $0.5\times {10}^8$ source photons. The $40\times 40{\mathrm{mm}}^2$ detector array with $1\times 1{\mathrm{mm}}^2$ 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 $\widehat{D}$, $\widehat{A}$, and $\widehat{\phi }$. Since it is obvious that $\overline{D}$ is unbiased from the developed theory and brute-force simulation, the term ${\overline{m}}_p$ is investigated instead of $\overline{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 (${\overline{m}}_p=\overline{A}/\overline{D}$) and $\overline{\phi }$ from the entire detector pixels at fixed modulation frequencies ($f_m$). 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 $m_p$ and $\phi$ are different for different detector pixels, PDFs of ${\widehat{m}}_p$ and $\widehat{\phi }$ are the same (with different means), as shown in Eqs. (12) and (13), respectively. Therefore, biased estimations would be clearly observed for ${\overline{m}}_p$ and $\overline{\phi }$ averaged from all detector pixels if the theoretical results derived in Sec. 2 are valid. For the comparison, all curves calculated at a fixed $f_m$ are scaled for both ${\overline{m}}_p$ and $\overline{\phi }$ to make the values of the photon number of $0.5\times {10}^8$ are the same. Figures 4(c) and (d) show that ${\overline{m}}_p$ and $\overline{\phi }$ are decreased and increased, respectively, as the number of simulated photons increases, which indicates that ${\overline{m}}_p$ and $\overline{\phi }$ are over- and under-estimated. Furthermore, it is observed that the degree of bias is worse for higher $f_m$, where the number of output photons contributing to the calculation for ${\overline{m}}_p$ and $\overline{\phi }$ is smaller. It is well known that ${\overline{m}}_p$ and $\overline{\phi }$ are decreased and increased, respectively, as $f_m$ 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, $\overline{D}$ and $\overline{A}$ are calculated from sinusoidal signals with white Gaussian noise of increased standard deviations. It is clear that $\overline{D}$ and $\overline{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) $m_p$ and (d) $\hat{\phi }$ 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 $\widehat{A}$ and $\widehat{\phi }$. 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.³¹ Alternatively, increasing the DC primary photon rate $c_p$ 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 $c_p$. The lower gain indicates the noise factor NF of the image intensifier should be increased. Resulting from this, $\gamma$ in Eq. (19) might not be increased enough to diminish the bias even for the increased amount of primary photons.

Percent differences for $\overline{A}$ and $\overline{\phi }$ 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.²⁹ For comparison, results of $\mathrm{NF}=4$ in Fig. 2 are repeatedly presented in Fig. 5. We assume NF linearly increases for decreasing the MCP gain ${\overline{k}}_n^I$ that makes the number of secondary photons constant for the increased primary photons in Fig. 5. This implies ${\sigma }^2\_\mathrm{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 $\overline{A}$ than $\overline{\phi }$. It is observed in Fig. 5(b) that the effect of increased NF is almost ignorable for cases of $\phi$ smaller than $3\pi /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 $4\sim 10$ and $4\sim 16$.

4.

Further Discussion and Conclusion

We observed that modulation amplitude and phase $\widehat{A}$ and $\widehat{\phi }$ are biased when they are estimated by Fourier transforming the noisy heterodyne outputs. Furthermore, the previously known SNR model for $\widehat{A}$ and $\widehat{\phi }$ 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 $\widehat{\phi }$ in diffusive imaging for small volume phantoms^9,16 because a high modulation frequency increases $\phi$ in Eq. (4). However, the theoretical result in this paper shows that high $\phi$ (i.e., the large amount of phase delay of the modulated incident beam to a phantom) deteriorates the reliability of the estimated phase $\widehat{\phi }$. Furthermore, the modulation depth $m_p$ in Eqs. (17) and (19) is dramatically decreased as the modulation frequency increases, which further degrades the estimation for $\widehat{\phi }$. Although recently developed noise models for heterodyne and homodyne measurements indicate that small $m_p$ reduces ${\sigma }_{\mathrm{AC}}^2$ in Eq. (4),²⁷ the decreased ${\sigma }_{\mathrm{AC}}^2$ does not improve the estimation for $\widehat{\phi }$ because $\gamma$ in Eq. (14) depends on only $a_{\mathrm{AC}}$ and ${\sigma }_{\mathrm{DC}}^2$. Therefore, it can be importantly noticed that increasing SNR for $\widehat{\phi }$ 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 $m_p$ as well as the DC photon rate $c_p$ are decreased because of increased scattering and absorbing in the diffusive medium. Therefore, the bias in $\widehat{A}$ and $\widehat{\phi }$ 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 $\widehat{A}$ and $\widehat{\phi }$ 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 $\widehat{\phi }$ 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.²⁰ 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 $\widehat{A}$ and $\widehat{\phi }$ 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 $\widehat{A}$ and $\widehat{\phi }$ is very similar with what is observed in this paper. For example, if the impulse output $a_p^{(i)}(t)$ is measured without amplification, ${\sigma }_{\mathrm{DC}}^2$ in Eq. (4) can be considered as ${\int }_{-T/2}^{T/2}{a}_p^{(i)}(t)/T$ based on that the variance of shot noise induced by a time-variant signal is described by time-averaging the signal.³².

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 approach²² 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.