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 and , where white Gaussian noise is added. Fourier transforming the noise-contaminated signals determines and that are random. Averaging these random variables by repeating brute-force simulations can estimate and . 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. and 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, is unbiased for the entire range of the noise level. However, 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 homogeneous phantom of and . 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 to , 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 of total photon number is the summation of eight outputs of source photons. The detector array with 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 , , and . Since it is obvious that is unbiased from the developed theory and brute-force simulation, the term is investigated instead of 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 () and from the entire detector pixels at fixed modulation frequencies (). 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 and are different for different detector pixels, PDFs of and are the same (with different means), as shown in Eqs. (12) and (13), respectively. Therefore, biased estimations would be clearly observed for 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 are scaled for both and to make the values of the photon number of are the same. Figures 4(c) and (d) show that and are decreased and increased, respectively, as the number of simulated photons increases, which indicates that and are over- and under-estimated. Furthermore, it is observed that the degree of bias is worse for higher , where the number of output photons contributing to the calculation for and is smaller. It is well known that and are decreased and increased, respectively, as 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.