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

发布时间：2024-02-09 11:00

1.软件版本

matlab2017b

2.仿真对比分析

1统计同步算法:   

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

基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第1张图片

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

基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第2张图片

 

 基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第3张图片

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

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:

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

基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第4张图片

 基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第5张图片

 

 

    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:

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

基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第6张图片

 基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第7张图片

 

 

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

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

基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第8张图片

 基于matlab频率估计算法对比，包括统计M.Westlund算法，BTDT，CZT，ZOOM-FFT 等的_第9张图片

 

 

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

相关推荐

华为昇腾MindSpore训练过程经验分享——D-Linknet

POM文件详解

Docker实战-部署GPE微服务的监控体系(二)

python中Requests请求的安装与常见用法

心寄开源，合规护航｜2022开放原子全球开源峰会开源合规分论坛即将开幕

java实现求 的近似值_Java 蒙特卡洛算法求圆周率近似值实例详解

微信小程序实现字母索引菜单

C#手动填充DataSet

假期学习小结day1

C语言之题目将一个字符串中的前N个字符复制到一个字符数组中去，不许使用strcpy

全球及中国液化天然气（LNG）行业建设规模现状及十四五产业规划报告2022年版

十大常用web前端UI组件库，赶紧收藏

用python画动态图_数据可视化｜绘制酷炫的GIF动态图（五）

权限提升-MY&MS&ORA等SQL数据库的提权

localhost:8080打不开原因解决方法

Windows启动谷歌浏览器Chrome失败(应用程序无法启动,因为应用程序的并行配置不正确)解决方法

逼真度超越「AI设计师」DALL·E 2！谷歌大脑推出新的文本生成图像模型Imagen

蚓链分红系统有哪些特点？如何进行收益分配？

Serverless 时代下微服务应用全托管解决方案

Javascript中的栈和堆

ItVuer - 免责声明 - 关于我们 - 联系我们

本网站信息来源于互联网，如有侵权请联系：561261067@qq.com

桂ICP备16001015号