## 2 Materials and Methods
### Bias and Variance Analysis with Conventional Cross-Correlation for TPE
It ensures that computations remain stable and consistent with real-world signal analysis, where measurements always occur over discrete frequency bins [25; 26]:  
\[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]0  
\[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]1  
\[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]2  
\[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]3 \[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]4 \[\hat{\tau_{0}}=\text{argmax}(L(-\tau))\text{ or }\hat{\tau_{0}}=\text{arg}\bigg{(} \frac{d}{d\tau}L(-\tau)=0\bigg{)} \tag{25}\]5\[ \begin{array}{c}\mathbb{E}\left[J_{CCF}(\tau_{0})^{2}\right]=8\pi^{2}\sigma_{N1}^{ 2}\sigma_{N2}^{2}\\ \times\bigg{\{}\frac{f_{0}^{2}}{\Delta^{2}f}\times SNR_{1}.SNR_{2}.\cos\Theta_{ 0}+\frac{f_{0}^{2}}{2B/T}\times(SNR_{1}+SNR_{2})\\ +\frac{B.\Delta f}{12(1/T)^{2}}\bigg{\}}\end{array}  \tag{38}\]  
The previous formulas also present \(\Theta_{0}\) as the argument of the complex number \(D_{1}(f_{0})D_{2}^{*}\) (\(f_{0}\)) and also of \(H_{1}(f_{0})H_{2}^{*}\) (\(f_{0}\)); and the signal-to-noise ratios, which present the ratio between the power spectral of the deterministic part over noise:  
\[SNR_{1}=\frac{\int_{-B}^{B}|D_{1}(f)|^{2}df}{\int_{-B}^{B}|N_{1}(f)|^{2}df}= \frac{|D_{1}(f_{0})|^{2}.\Delta^{2}f}{\sigma_{n_{1}}^{2}};SNR_{2}=\frac{\int_ {-B}^{B}|D_{2}(f)|^{2}df}{\int_{-B}^{B}|N_{2}(f)|^{2}df}=\frac{|D_{2}(f_{0})|^ {2}.\Delta^{2}f}{\sigma_{n_{2}}^{2}} \tag{39}\]  
Therefore, by applying Formulas (4) and (5), the bias and variance of the estimation are derived as:  
\[ \begin{array}{c}\mbox{\it bias }(\hat{\tau}_{0})_{CCF}=-\frac{\mathbb{E}[J_{CCF }(\tau_{0})]}{\mathbb{E}\big{[}J_{CCF}^{\prime}(\tau_{0})\big{]}}=-\frac{1}{2 \pi f_{0}}\times\frac{\mbox{\rm Im }\left\{D_{1}(f_{0})D_{2}^{*}(f)\right\}}{\mbox{\rm Re}\left\{D_{1}(f_{0})D_{2 }^{*}(f_{0})\right\}}=-\frac{\tan\Theta_{0}}{2\pi f_{0}}\\ \\ \mbox{\it var }(\hat{\tau}_{0})_{CCF}=\frac{\mathbb{E}[J_{CCF}(\tau_{0})^{2}]}{ \mathbb{E}[J_{CCF}^{\prime}(\tau_{0})]^{2}}=\\ =\frac{2}{\pi^{2}f_{0}^{2}}\\ \times\Big{\{}\tan^{2}\Theta_{0}+\frac{1}{N_{d}.\cos^{2}\Theta_{0}}\Big{(} \frac{1}{SNR_{1}}+\frac{1}{SNR_{2}}\Big{)}\\ +\frac{B^{2}}{N_{d}.\delta f_{0}^{2}}\times\frac{1}{SNR_{1}.SNR_{2}.\cos\Theta _{0}}\Big{\}}\end{array}  \tag{41}\]  
where the number \(N_{d}=\frac{2B}{\Delta f}\) defines the number of discrete frequency bins \(N_{d}\) used in spectral analysis. This formulation arises from the necessity to sample the continuous frequency spectrum at discrete intervals, ensuring that all relevant signal components within the bandwidth are accounted for in the Fourier-based time period estimation. A finer frequency resolution (\(\Delta f\)) increases \(N_{d}\), leading to improved accuracy in spectral estimation but also higher computational complexity. Conversely, a larger \(\Delta f\) reduces \(N_{d}\), potentially impacting estimation precision. The choice of \(\Delta f\) thus balances the trade-off between resolution and efficiency, ensuring that the cross-correlation-based TPE method effectively captures the required spectral information while maintaining computational feasibility [23; 34]. It is not difficult to see that, in ideal conditions, time constants of \(H_{1}(s)\), \(H_{2}(s)\) coincide in such condition that \(\tau_{1,th}=\tau_{2,th},\tau_{1,\varepsilon}=\tau_{2,\varepsilon}\); \(\arg(H_{1}(f))=\arg(H_{2}(f))\); then, \(\Theta_{0}=\arg(D_{1}(f_{0})D_{2}^{*}(f_{0}))=\arg(H_{1}(f_{0})H_{2}^{*}(f_{0} ))=0\).