## 4 Numerical Implementation and Analysis of the Phase-Locked Loop Statistics  
In this section, numerical implementation and analysis of the first-order phase-locked loop statistics are presented by solving the normalized FP equation. Computer codes were developed using the Matlab(r) software. According to Equation (9), the initial condition is supposed to be a Dirac delta function. Since the delta function is a continuous-time signal and therefore, for numerically resolved purpose, the value \(1/\Delta\phi\) is used for approximating the delta function. In the following illustrative examples, the initial value of phase error \(\phi_{0}=0\) is utilized for discussion. First of all, the time-evolution of PDFs is investigated to reveal detailed evolution of the time-dependent PDFs of the PLL. Typical results for the time-evolution of nonmodulo-\(2\pi\) transient statistics for SNR \(\alpha=1\) and detuning factor \(\gamma=0\) are shown in Figures 2-4. The grid spacing for the phase-error domain is chosen as \(\Delta\phi=\pi/50\) rad; the grid spacing in the time marching direction is set as 0.01. Figure 2 shows the transition phase-error density function \(g(\phi,\tau)\) and the PDF \(p(\phi,\tau)\). The numbers beside individual curves indicate the normalized times, switching \(\tau\) from 0.625 to 20. The corresponding cumulative probability distribution of the phase-error is depicted in Figure 3. Typical time-evolution of \(p(\phi,\tau)\), as shown in Figure 2b, revealed that as P diffused with time, the probability that the phase would jump out of the primary interval \(-\pi<\phi<\pi\) was the same as the probability that it would jump into that interval, implying that the possibility of larger phase error increased. In the case of \(\gamma=0\), both of the error density functions \(g(\phi,\tau)\) and \(p(\phi,\tau)\) diffuse with time while remaining symmetric with respect to the phase-error axis, and thus the initial global peak values of the density functions \(g\) and \(p\) at zero phase error (\(\phi=0\)) decreased as time progresses. As functions of \(\phi\), the local peak positions of \(p\) appeared at a period of \(2\pi\) rad. The amplitudes of local peaks were monotonically nonincreasing when moving outwards from their mean values with respect to the phase-error axis. The figures and later ones as well only show part of the results. In Figures 2 and 3, the results at a phase error between \(-20\) and \(20\) rad were the area of interest. Figure 4 provides the three-dimensional plots for nonmodulo-\(2\pi\) transition \(g\) and \(p\), respectively, for more profound insights into the statistics information, especially the relation between \(g/p\) and phase error, as time progresses. Both the front view and back view are shown such that better understanding on the behavior of PLL statistics is accessible. The result enabled us to mathematically quantify the difference among different PDFs, providing a meaningful conceptual link between the stochastic process and information geometry, and could be utilized for optimizing various desired outcomes. Figures 5-7 provide information basically similar to that of Figures 2-4, and thus the parameter values were identical except that the detuning factor was now changed to \(\gamma=\sin(\pi/4)\). The identical numerical values including the normalized time (\(\tau\)), the grid size for the phase-error domain (\(\Delta\phi\)) and the grid size in the time marching direction (\(\Delta\tau\)) were utilized for comparison. It can be seen that in this case the density functions diffused to one direction as time progresses, given a certain SNR (e.g,\(\alpha=1\) as demonstrated here). The error probabilities \(g\) and \(p\) diffused virtually in opposite directions. In this example, \(g\) virtually diffused leftwards while \(p\) diffused rightwards. As can be seen, \(p\) had a right-skewed (or positively skewed) distribution, and most data fell to the right (positive direction) of the graph's peak.