基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOMFFT等的

matlab2017b

1统计同步算法:   

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

最后得到的信号的眼图为：

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

clc;
clear;
close all;
warning off;
addpath 'func\'
%======================统计算法实现===========================
%测试次数
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('s\n');
save r4.mat t2 p

BTDT:

这个部分的仿真结果如下所示：

    BTDT的基本步骤如下所示：

    软同步方案中，为了重建眼图，需要确定采样步长Δt，Δt是通过Δt = FT/N，来确定的，T为信号周期的步长，N是采样点的数目，F是采样数据的包络。通过对这些采样的数据进行傅立叶变换，我们可以利用数值的方法确定包络的数目F’，可是离散傅立叶变换的谱分辨率有限，实际的F在F’- 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;
addpath 'func\'
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.5f\n',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('s\n');
save r2.mat t p

CZT:

这个部分的仿真结果如下所示：

CZT的基本步骤如下所示：

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

clc;
clear;
close all;
warning off;
addpath 'func\'
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('s\n');
save r3.mat t p

ZOOM_FFT

这个部分的仿真结果如下所示：

clc;
clear;
close all;
warning off;
addpath 'func\'
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('s\n');
save r4.mat t p

A15-03