1.算法仿真效果
matlab2022a仿真结果如下：
RBS算法的累积误差在50us大于80%；ERBS算法的累积误差在25us大于80%；本文所提出的MRF-MAP-Gardner ERBS算法的累积误差在10us大于80%。因此，在降低同步误差方面，MRF-MAP-Gardner ERBS算法性能优于RBS算法和ERBS算法。
RBS算法的同步误差在40us到100us所在比重最大，ERBS算法的的同步误差则主要集中在20us以内，而MRF-MAP-Gardner ERBS算法的同步误差集中在18us之内.
2.算法涉及理论知识概要
ERBS是在RBS算法的基础上改进得到的，但其性能上依旧存在如下几个方面的缺陷：在非确定性相位偏差估计方面，EBRS算法假定了相位偏差测量值服从均值为0，方差为的高斯分布，此外，其非确定性相位偏差估计是基于有限个相位偏差测量值来获得的，因此这种方法当WSN网络环境发生较大程度的改变的时候，则会产生较大的估计误差。在时钟漂移估计方面，EBRS算法将其节点时钟漂移模型简化为线性模型进行估计，而实际中，时钟漂移是一个动态变化的过程，对时钟漂移的估计需要进行实时的跟踪和估计才能获得一个较为精确的估计结果。针对这两个问题，提出了一种新型的高精度能量有效无线传感器网络时间同步算法。
算法通过一种基于马尔科夫随机场的最大后验估计方法对无线传感器网络中不相邻的两个接收节点在多个参考广播消息条件下的相位偏差进行估计，然后通过Gardner环对时钟漂移进行实时的估计和跟踪，实现时间同步过程。
在非确定性相位偏差估计方面，本文采用了一种基于马尔科夫随机场的最大后验估计算法。定义WSN中任意两个不相邻的的节点i和节点j交换m个参考消息包的相位偏差为：
在时钟漂移估计方面，由于WSN中各个节点的时钟漂移是由于晶振的频率漂移导致的，假设一个晶振的频率为100Mhz，其频率漂移为20ppm，那么晶振的实际输出频率为100Mhz2Khz。因此，接收节点将存在接收信号与本节点频率不同步的问题。同时，即便是同一型号的晶振，其频率漂移也是存在偏差的，因此需要通过一种实时跟踪和估计的方法来获得各种晶振的时钟漂移值。
3.MATLAB核心程序
```function yout=func_gardner(data0);
nsamp = 2;
data = zeros(1,nsamp*length(data0));
for i = 1:length(data)-2
data(i) = data0(floor(i/nsamp)+1);
end
interplen = length(data);
w = 0.5;%lf输出
q0 = 0.9;
q = zeros(1,interplen);
m = 5;
j = 3;
C1 = 0.001;
C2 = C1 0.0005;
datarcosI = data(9:end) ;
datarcosQ = data(9:end) ;
s0 = 2;%插值乘法器值
Detector_out = zeros(1,interplen);%误差检测输出
interp_outI = zeros(1,interplen 2);%插值filter输出;
interp_outQ = zeros(1,interplen * 2);
interp_outI(1) = datarcosI(2);
interp_outI(2) = datarcosI(4);
interp_outQ(1) = datarcosQ(2);
interp_outQ(2) = datarcosQ(4);
for i = 2 : interplen - 12
for k = 1 : 1
y_temp = q0 - w;%=============================
q(m) = q0;
if y_temp > 0
q0 = y_temp;
else
q0 = mod(y_temp,1);
mk = m;
uk = s0 q(m);
uu(j) = uk;
data1 = datarcosI(mk);
data2 = datarcosI(mk + 1);
interp_outI(j) = uk data2 + (1 - uk) data1;
if mod(j,2) ~= 0
qoutI((j+1)/2) = interp_outI(j);
end
data1 = datarcosQ(mk);
data2 = datarcosQ(mk + 1);
interp_outQ(j) = uk data2 + (1 - uk) data1;
if mod(j,2) ~= 0
qoutQ((j+1)/2) = interp_outQ(j);
end
j = j + 1;
end
m = m + 1;
end
Detector_out(i) = (interp_outI((i - 2) 2 + 2)) ((interp_outI((i - 1) 2 + 1)) - (interp_outI((i - 2)* 2 + 1)))...
+ (interp_outQ((i - 2)* 2 + 2)) * ((interp_outQ((i - 1)* 2 + 1)) - (interp_outQ((i - 2)* 2 + 1)));
w = w + (C1 (Detector_out(i) - Detector_out(i - 1)) + C2 Detector_out(i));
Wm(i) = w;
end
yout = (interp_outI + interp_outQ)/2/nsamp;
```
