matlab2017b

# 2.仿真对比分析

1统计同步算法

统计同步算法的基本思路，主要是通过多次采样测试，然后计算对应的概率分布，来确定其同步时刻。测试信号和频率点为：

其主要思路，其实就是对一段正弦信号进行多次采样，然后计算每次采样后得到的跟踪结论，然后计算对应的概率分布情况，最后根据其概率密度函数得到最后的时钟跟踪点。

``````clc;
clear;
close all;
warning off;
%======================统计算法实现===========================
%测试次数
Sample_Time = 100;
for j = 1:Sample_Time
j
RandStream.setDefaultStream(RandStream('mt19937ar','seed',j));

%产生伪随机数
freqcarrier   = 40e6+round(5000*randn(1,1));
freqSignal    = 10e6;
freqSample    = 640e6;
K             = floor(freqSample/freqSignal);
numSample     = 256;
periodSample  = 1/freqSample;
sampleIndex   = 0:numSample-1;
timeSequence  = sampleIndex/freqSample;

Data1         = 2*(randn(1,numSample)>=0.5)-1;
Data1s        = func_samples(Data1,K);
%模拟伪随机数，即随机数以周期性出现
msg2          = [Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s];

b             = hanning(127);
msg           = filter(b,1,msg2);
msg           = msg/max(msg);
msg(1:1024)   = [];
%调制
ff            = cos(2*pi*freqcarrier.*[0:length(msg)-1]/freqSample);
signalSample  = msg.*ff;
t             = length(signalSample);
[f,sf]        = T2Fv2(t,signalSample);
figure(1);
subplot(311);
plot(msg);
title('测试随机数');
axis([1,length(msg),-1.5,1.5]);
subplot(312);
plot(f,abs(sf));
xlabel('频率 Mhz');
subplot(313);
plot(sf);
xlabel('归一化频率 点数');
%==========================计算Fcourse========================
sf1(1)  = 0;
index   = find(sf== max(sf));
I       = 3;
Fcourse = index;
tic;

Fth = Fcourse;
N   = length(signalSample);

Ntemp         = N*(Fcourse-0.5)/Fth;
numTemp       = round(Ntemp);
[ftemp,sftemp]= T2Fv2(t,signalSample(1:N));
sftemp(1)     = 0;
indexTemp     = find(sftemp== max(sftemp));
if sftemp(indexTemp) > sftemp(indexTemp-1)
Fth = Fth+0.1;
else
Fth = Fth-0.1;
end
Fpresize2(j)=Fth;
t(j) = toc;
end
figure;
hist(Fpresize2,10);
%计算概率分布
[m,n]    = hist(Fpresize2,100);
[V,I]    = max(m);
Fpresize = Fpresize2(I);
fprintf('%12.8f',Fpresize);

figure;
%恢复眼图
delta  = (16320)/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(211)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('标准眼图');
grid on;

%恢复眼图
delta  = Fpresize/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(212)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('同步之后眼图');
grid on;

p=100*abs(Fpresize-16320)/16320;
fprintf('估计精度：');
fprintf('%2.4f',p);
fprintf('%%n');

t2 = mean(t);
fprintf('仿真时间：');
fprintf('%2.4f',t2);
fprintf('sn');
save r4.mat t2 p

``````

BTDT:

BTDT的基本步骤如下所示：

同步方案中，为了重建眼图，需要确定采样步长ΔtΔt是通过Δt = FT/N来确定的，T为信号周期的步长，N是采样点的数目，F是采样数据的包络。通过对这些采样的数据进行傅立叶变换，我们可以利用数值的方法确定包络的数目F’可是离散傅立叶变换的谱分辨率有限，实际的FF’- 0.5 to F’+0.5之间。为重建眼图，我们必须得到精确的F值。

实际的F值可以采用我们最近提出的二进制数据截断法BTDT迅速而有效地得到其工作原理如下：首先采样得到的数据的长度被截断为Nth = N/2F，因而数据的的数目减为N’ =N-Nth。接着对截断的数据进行傅立叶变换来计算其功率谱，并比较P(F’)  P(F’-1)P(f)为频率f的谱功率。如果P(F’-1) > P(F’)，实际的F值应该位于小于F’的范围内。

相反，如果P(F’-1) < P(F’)，如果F值应该高于F’。以相同的方式，我们重复截短数据Nth =N*[Fth-(F’-0.5)]/Fth 倍，Fth为预测的F范围的中心频率。经过m次迭代计算，F的准确度提高了2m倍。采用BTDT方法，采用标准的台式机，眼图重建可以在<0.3s的时间内完成，因此使每次扫描的实时再同步成为可能。

``````clc;
clear;
close all;
warning off;
RandStream.setDefaultStream(RandStream('mt19937ar','seed',1));

%产生伪随机数
freqcarrier   = 40e6;
freqSignal    = 10e6;
freqSample    = 640e6;
K             = freqSample/freqSignal;
numSample     = 256;
periodSample  = 1/freqSample;
sampleIndex   = 0:numSample-1;
timeSequence  = sampleIndex/freqSample;

Data1         = 2*(randn(1,numSample)>=0.5)-1;
Data1s        = func_samples(Data1,K);
%模拟伪随机数，即随机数以周期性出现
msg2          = [Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s];

b             = hanning(127);
msg           = filter(b,1,msg2);
msg           = msg/max(msg);
msg(1:1024)   = [];
%调制
ff            = cos(2*pi*freqcarrier.*[0:length(msg)-1]/freqSample);
signalSample  = msg.*ff;
t             = length(signalSample);
[f,sf]        = T2Fv2(t,signalSample);

figure;
subplot(311);
plot(msg);
title('测试随机数');
axis([1,length(msg),-1.5,1.5]);
subplot(312);
plot(f,abs(sf));
xlabel('频率 Mhz');
subplot(313);
plot(sf);
xlabel('归一化频率 点数');

%==========================计算Fcourse========================
sf1(1)  = 0;
index   = find(sf== max(sf));
I       = 3;
Fcourse = index;
tic;
%==========================利用BTDT计算Fpresize================
Fth = Fcourse;
N   = length(msg);
for k = 1:I
Ntemp = N*(Fcourse-0.5)/Fth;
numTemp = round(Ntemp);
[ftemp,sftemp]=T2Fv2(t,signalSample(1:numTemp));
sftemp(1) = 0;
indexTemp = find(sftemp== max(sftemp));
if sftemp(indexTemp)>sftemp(indexTemp-1)
Fth = Fth+(0.5)^(k+1);
else
Fth = Fth-(0.5)^(k+1);
end
end
format long;
Fpresize = Fth;
fprintf('%6.5fn',Fpresize);
clc;
t=toc;

figure;
%恢复眼图
delta  = (Fcourse)/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(211)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('标准眼图');
grid on;

%恢复眼图
delta  = Fpresize/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(212)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('同步之后眼图');
grid on;

p=100*abs(Fpresize-Fcourse)/Fcourse;
fprintf('估计精度：');
fprintf('%2.4f',p);
fprintf('%%n');

fprintf('仿真时间：');
fprintf('%2.4f',t);
fprintf('sn');
save r2.mat t p

``````

CZT:

CZT的基本步骤如下所示：

http://d.wanfangdata.com.cn/Periodical_zjbgcxyxb201201013.aspx

``````clc;
clear;
close all;
warning off;
RandStream.setDefaultStream(RandStream('mt19937ar','seed',1));

%产生伪随机数
freqcarrier   = 40e6;
freqSignal    = 10e6;
freqSample    = 640e6;
K             = freqSample/freqSignal;
numSample     = 256;
periodSample  = 1/freqSample;
sampleIndex   = 0:numSample-1;
timeSequence  = sampleIndex/freqSample;

Data1         = 2*(randn(1,numSample)>=0.5)-1;
Data1s        = func_samples(Data1,K);
%模拟伪随机数，即随机数以周期性出现
msg2          = [Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s];

b             = hanning(127);
msg           = filter(b,1,msg2);
msg           = msg/max(msg);
msg(1:1024)   = [];
%调制
ff            = cos(2*pi*freqcarrier.*[0:length(msg)-1]/freqSample);
signalSample  = msg.*ff;
t             = length(signalSample);
[f,sf]        = T2Fv2(t,signalSample);

figure;
subplot(311);
plot(msg);
title('测试随机数');
axis([1,length(msg),-1.5,1.5]);
subplot(312);
plot(f,abs(sf));
xlabel('频率 Mhz');
subplot(313);
plot(sf);
xlabel('归一化频率 点数');
%==========================计算Fcourse========================
sf1(1)  = 0;
index   = find(sf== max(sf));
I       = 3;
Fcourse = index;

tic;
%==========================CZT实现===========================
F1        = 2e7;%细化频率段起点
F2        = 6e7;%细化频率段终点
M         = 2^16;%细化频段的频点数，（这里其实就是细化精度）
w         = exp(-j*2*pi*(F2-F1)/(freqSample*M));%细化频段的跨度（步长）
a         = exp(j*2*pi*F1/freqSample);%细化频段的起始点，这里需要运算一下才能代入czt函数
xk        = czt(signalSample,M,w,a);
h         = 0:1:M-1;%细化频点序列
deltaFreq = freqSample/length(signalSample);
index1    = F1/deltaFreq;
index2    = F2/deltaFreq;
f0        =(F2-F1)/M*h+F1;%细化的频率值
f00       =(index2-index1)/M*h+index1;%细化的频率值
Fpresize  = f00(find(xk==max(xk))-1);%%%%不对
fprintf('%6.5f',Fpresize);
t=toc;

figure;
%恢复眼图
delta  = (Fcourse)/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(211)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('标准眼图');
grid on;

%恢复眼图
delta  = Fpresize/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(212)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('同步之后眼图');
grid on;

p=100*abs(Fpresize-Fcourse)/Fcourse;
fprintf('估计精度：');
fprintf('%2.4f',p);
fprintf('%%n');

fprintf('仿真时间：');
fprintf('%2.4f',t);
fprintf('sn');
save r3.mat t p

``````

ZOOM_FFT

``````clc;
clear;
close all;
warning off;
RandStream.setDefaultStream(RandStream('mt19937ar','seed',1));

%产生伪随机数
freqcarrier   = 40e6;
freqSignal    = 10e6;
freqSample    = 640e6;
K             = freqSample/freqSignal;
numSample     = 256;
periodSample  = 1/freqSample;
sampleIndex   = 0:numSample-1;
timeSequence  = sampleIndex/freqSample;

Data1         = 2*(randn(1,numSample)>=0.5)-1;
Data1s        = func_samples(Data1,K);
%模拟伪随机数，即随机数以周期性出现
msg2          = [Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s,Data1s];

b             = hanning(127);
msg           = filter(b,1,msg2);
msg           = msg/max(msg);
msg(1:1024)   = [];
%调制
ff            = cos(2*pi*freqcarrier.*[0:length(msg)-1]/freqSample);
signalSample  = msg.*ff;
t             = length(signalSample);
[f,sf]        = T2Fv2(t,signalSample);

figure;
subplot(311);
plot(msg);
title('测试随机数');
axis([1,length(msg),-1.5,1.5]);
subplot(312);
plot(f,abs(sf));
xlabel('频率 Mhz');
subplot(313);
plot(sf);
xlabel('归一化频率 点数');
%==========================计算Fcourse========================
sf1(1)  = 0;
index   = find(sf== max(sf));
I       = 3;
Fcourse = index;
tic;
%==========================ZOOMFFT===========================
F1        = 2e7;%细化频率段起点
F2        = 6e7;%细化频率段终点
M         = 2^16;%细化频段的频点数，（这里其实就是细化精度）

fi        = 1e8;   %最小细化截止频率
np        = 32;    %放大倍数
nfft      = M;
clc;
xk        = zfft_m(signalSample,fi,freqSample,nfft,np);
t=toc;
h         = 0:1:M-1;%细化频点序列
deltaFreq = freqSample/length(signalSample);
index1    = F1/deltaFreq;
index2    = F2/deltaFreq;
f0        = (F2-F1)/M*h+F1;%细化的频率值
f00       = (index2-index1)/M*h+index1;%细化的频率值
Fpresize  = f00(find(xk==max(xk))+1);
fprintf('%6.5f',Fpresize);

figure;
%恢复眼图
delta  = (Fcourse)/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(211)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('标准眼图');
grid on;

%恢复眼图
delta  = Fpresize/length(msg);
Xpoint = mod((1:length(msg))*delta,64);
Ypoint = msg;
subplot(212)
tip = ceil((1:length(msg))*delta);
for k = 1:max(tip)
indexSS = find(tip==k);
plot(Xpoint(indexSS),Ypoint(indexSS));
hold on;
end
axis([2.3,18.4,-0.96,0.96]);
title('同步之后眼图');
grid on;

p=100*abs(Fpresize-Fcourse)/Fcourse;
fprintf('估计精度：');
fprintf('%2.4f',p);
fprintf('%%n');

fprintf('仿真时间：');
fprintf('%2.4f',t);
fprintf('sn');
save r4.mat t p``````

A15-03

THE END