clc;
N=1000; % 测量次数X采样时间,就是滤波时间常数
M=100000; % 仿真时间
%F_bias=1e-12; % 晶振最大频偏
F_oxco_error=1e-12; % 晶振短期频率稳定度
F_rub_error=1e-15; % 铷钟频率稳定度
Fr=1e7; % 中心频率
F_ini_error=1e-8; % 初始频差
%vol_adj=1e-10; % 频率调节灵敏度1e-9v/s,要达到1e-15的频率准确度,则要求1e-9/1e-15=1e6,即需要20为的DAC
devia_rub=double(zeros(1,N+1)); % 铷钟的瞬时频率差,初始瞬时相位差为零,即devia_rob(1)=0
devia_oxco=double(zeros(1,N+1)); % 晶振的瞬时频率差,初始瞬时相位差为零,即devia_oxco(1)=0
devia_oxco(1,1)=0; % 初始相位差
adjust_voltage=double(zeros(1,M)); % 压控电压
devia_rub_oxco=double(zeros(1,M)); % 铷钟和晶振相位差
devia_rub_oxco_before=0;
F_error=double(zeros(1,M)); % 两者频差
%sampling_points=1:1e-3:100; % 采样时间间隔为1e-3,采样次数为100,且连续采样
sampling_points=1:M;
F_ini_error2=0;
flag=0;
%ddd=double(zeros(1,M));
ccc=double(zeros(1,M));
%xxx=double(zeros(1,M));
%bb=double(zeros(1,N));
derta_T=1e-3;
for i=1:M
fact11=normrnd(0,F_rub_error,N+1,1);
f_rub_fact=(fact11+1)*Fr;
%f_rub_fact=(normrnd(0,F_rub_error,N+1,1)+1)*Fr; % 铷钟瞬时频率加入高斯频率振荡,标准差为1e-10
fact22=normrnd(0,F_oxco_error,N+1,1);
f_oxco_fact=(fact22+1+F_ini_error)*Fr;
%f_oxco_fact=(1+F_ini_error)*Fr;
%f_oxco_fact=(normrnd(0,F_oxco_error,N+1,1)+1+F_ini_error)*Fr; % 晶振瞬时频率加入高斯频率振荡,标准差为1e-11,且加入了初始频差
devia_rub_sum=0.0; % 铷钟累计相位差
devia_oxco_sum=0.0; % 晶振累计相位差
for j=2:N+1 % 每1e-3秒测量一次,工测量N次,求平均频率差
devia_rub(j)=devia_rub(j-1)+(1e-7-1.0/f_rub_fact(j));%*f_rub_fact(j); %单次测量的频差1s内的相差
%devia_rub_sum=(round(devia_rub(j)*1e11))/1e11+devia_rub_sum; % 1e-10为0.1ns,为系统测量时间时间的精度
devia_oxco(j)=devia_oxco(j-1)+(1e-7-1.0/f_oxco_fact(j));%*f_oxco_fact(j);
%bb(j-1)=1.0/f_oxco_fact(j);
%devia_oxco_sum=devia_oxco_sum+(round(devia_oxco(j)*1e11))/1e11;
end
% N*dater_T时间内产生的相位差
%devia_rub_oxco(1,i)=devia_rub(N+1)-devia_oxco(N+1); % 1ms时间误差的平均差值,相当于数字滤波器,即锁相环中的滤波电路
devia_rub_oxco(1,i)=round(devia_rub(N+1)*1e11)/1e11-round(devia_oxco(N+1)*1e11)/1e11; % 1ms时间误差的平均差值,相当于数字滤波器,即锁相环中的滤波电路
% 频率校正完成,置标志
if((abs(devia_rub_oxco(1,i)-devia_rub_oxco_before)<=1e-11) && (abs(F_ini_error)<=1e-11))
flag=1;
end
if(flag) % 频率校正完成后的频率修改及相位一致性调整
if((abs(devia_rub_oxco(i)-devia_rub_oxco_before)<=1e-11) && (abs(F_ini_error)<=1e-11)) % 判断频率准确度在要求的范围内
%if((abs(F_ini_error)>8e-12)) % 频率准确度临近范围的边界值,不改变VCO频率,跳出至此循环
% continue;
% end
%if((devia_rub_oxco(i)>3e-10) && (F_error(1,i)>8e-12))
% continue;
% end
% if((devia_rub_oxco(i)<-3e-10) && (F_error(1,i)<-8e-12))
% continue;
% end
% 相位一致性要求没有满足,通过微调频率,改变两者之间的相位差
if((devia_rub_oxco(i)>3e-10) && (F_error(1,i)<8e-12))
F_error(1,i)=F_ini_error+1e-12;
elseif((devia_rub_oxco(i)<-3e-10) && (F_error(1,i)>-8e-12))
F_error(1,i)=F_ini_error-1e-12;
elseif((F_error(1,i)~=0) && (abs(devia_rub_oxco(i))<3e-10))
F_error(1,i)=F_ini_error-F_ini_error(1,i)/abs(F_ini_error(1,i))*1e-12;
end
else % 准确度超过要求的方位,改变频率使之满足
if(F_ini_error>=0)
F_error(1,i)=F_ini_error-1e-12;
else
F_error(1,i)=F_ini_error+1e-12;
end
end
F_ini_error=F_error(1,i);
else % 频率校正未完成前的频率调节
% 根据铷钟和晶振的相位差的变化值,来调节压控电压值
adjust_voltage(i)=(round(((devia_rub_oxco(i)-devia_rub_oxco_before)*1e9)*1e4))/1e4;
% 改变VCO的输出频率
if((abs(F_ini_error-adjust_voltage(1,i)*1e-9)<abs(F_ini_error+adjust_voltage(1,i)*1e-9)))
F_error(1,i)=F_ini_error-adjust_voltage(1,i)*1e-9; % VCO改变输出频率
else
F_error(1,i)=F_ini_error+adjust_voltage(1,i)*1e-9; % VCO改变输出频率
end
F_ini_error=F_error(i);
end
% 保存此时的两者相位差值;
devia_rub_oxco_before=devia_rub_oxco(i);
%ddd(i)=devia_oxco_sum/N;
if(i>1)
ccc(i)= devia_rub_oxco(i)- devia_rub_oxco(i-1);
end
% xxx(i)=devia_rub_sum/N;
devia_rub(1,1)=devia_rub(1,N+1);
devia_oxco(1,1)=devia_oxco(1,N+1);
end
subplot(3,1,1);
plot(sampling_points,devia_rub_oxco);
title('deviation rub & oxco');
subplot(3,1,2);
plot(sampling_points,F_error);
title('Frequency error');
subplot(3,1,3);
plot(sampling_points,adjust_voltage);
title('adjust voltage');