## 2 System Model
### Linear Interpolation for CPR  
Linear interpolation is a straightforward regression method that is often used as a reference point in research. It serves as the foundation for experiments aimed at mitigating the impact of outliers that degrade phase recovery. The process begins by computing the coarse phase between consecutive pilot positions using the argument of the product between the received signal \(y[n]\), at pilot positions \(k=nL\), where \(k=0,L,\ldots,K\), and the conjugate of the corresponding pilot sequence, \(p[k]\). The pilot sequence is introduced using the \(\pi/4\)-rotated BPSK modulation, as previously described, where \(k\) denotes the pilot positions. Specifically, the estimated phase at each pilot position is given by:  
\[\varphi[k]=\arg(y[Lk]p[k]^{*}), \tag{6}\]  
where \([\cdot]^{*}\) denotes the conjugate operator and \(\varphi[k]\) represents the coarse phase in the pilot position, with pilots spaced in every \(L=32\) samples. This burst length can be modified according to the resource conditions. The estimated phase, \(\varphi[k]\), can be used directly for linear interpolation, eliminating the need for a hard-decision filter. However, this approach may lead to degraded performance under certain conditions. In linear interpolation, it is assumed that the estimated phase \(\hat{q}[n]\) is on the straight line connecting the nearest pilot phases \(p[k]\) and \(p[k+1]\), as illustrated in Figure 3.  
The linear interpolated value \(\hat{q}[n]\) is given by [36]:  
\[\hat{q}[n]=\varphi[k]+\left(\frac{\varphi[k+1]-\varphi[k]}{L}\right)\left([n- Lk]\right),\ \forall Lk<n\leq L(k+1). \tag{7}\]  
In Algorithm 1, \(\varphi[k]\) denotes the value of the coarse phase at the pilot position, and \(n\) represents the position of the received symbols. Modification of linear interpolation using a hard-decision filter involves the use of a threshold \(\varphi_{th}\), which is estimated as \(\varphi_{th}=\pi\sigma_{l}\), where \(\sigma_{l}\) represents the standard deviation of phase noise. The calculation of \(\sigma_{l}\) is based on the burst length plus one \((L+1)\) and the noise variance \(\sigma_{w}^{2}\) of the random phase, and is given by:  
\[\sigma_{l}=\sqrt{(L+1)\sigma_{w}^{2}}. \tag{8}\]  
Figure 3: Example of linear interpolation (blue line) used to obtain the phase tracking of a signal (red line). The red dots represent the pilots. ```
1:\(L\)\(\triangleright\) Set the predetermined burst length
2:\(y[n]\)\(\triangleright\) The received symbols plus phase noise
3:\(p[k]\)\(\triangleright\) The pilot sequence inserted like \(BPSK\) symbols
4:\(K\)\(\triangleright\) Total number of pilots
5:\(N\)\(\triangleright\) Total number of received symbols
6:\(\triangleright\) Coarse Phase Estimation
7:for\(n=0\) to \(N\)do
8:for\(k=0\) to \(K\)do
9:\(\varphi[k]=\arg(y[Lk]p[k]^{*})\)\(\triangleright\) Compute coarse phase for pilots
10:\(\triangleright\) ****************Hard-decision filter for outliers discrimination ****************
11:if\(k\geq 1\)then
12:\(\Delta\varphi[k]=\varphi[k]-\varphi[k-1]\)\(\triangleright\) Phase difference between consecutive pilots
13:if\(\Delta\varphi[k]>\varphi_{th}\)then
14:\(\varphi_{u}[k]=\varphi[k-1]+\varphi_{th}\)\(\triangleright\) Update phase upward for outliers



### Robust Locally Weighted Regression

The robust locally weighted regression method efficiently smooths the received symbols using a burst-based approach, excluding the necessity to process the entire signal or substantial segments of it. This method preferentially assigns greater weights to input symbols characterized by lower phase noise variance, thus reducing the impact of symbols affected by higher phase noise. The robust method effectively identifies and eliminates outliers, mitigating distortions in the smoothed data points. This approach employs an iterative reweighting of the least-squares method [30; 37].