基于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

相关推荐

scss重置浏览器默认样式 normalize.css

分享 9 个不错的 Vue 技巧

TiDB Operator 准入控制器

ES6模块化

LNMP架构搭建

微前端vue项目实战 -- 乾坤qiankun框架 （一）

Centos 7搭建LNMP架构及部署Discuz论坛

程序员必会的核心基础知识：1张导图+10本书

javascript基础#DOM增删改练习 超级详细！基础简单

Docker 容器 ----介绍&安装&基本命令

【GNN】GNN图神经网络工具箱的应用和matlab仿真

千万不要把Request传递到异步线程里面！有坑！

如何在react中使用pixi快速的做一个的炫酷的刹车特效！（全注释，附源码）

嵌入式开发-各种干货

以太坊是什么?｜猿创征文

Spring Cloud Alibaba——Nacos进行多配置文件的加载

使用docker-compose部署Nginx反向代理

MySQL中符号@的作用

DOM增删改

vue-element-ui前后端交互实现分页查询

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

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

桂ICP备16001015号