gardner算法算法matlab仿真请教（附程序）

cameion

看了不少gardner位同步的资料，阅读多个gardner算法的matlab仿真程序，始终没有完全掌握。网上找到一个比比较完整的gardner仿真程序，但对不同测试输入数据的仿真结果难以解释，请高手指点。附上gardner的matlab程序。

sunyuanxin

到底有什么样的数据难以解释？你总要一条一条的列出来吧？这么打一个包放在这里，由谁会真正的去下载？

cameion

%%%%%%%%%%原程序中的仿真数据产生方法。这种方法产生的数据，运行后从曲给上看可以收敛。                                                                                                                                                 法                                                         
N=6000;  %符号数
K=4;     %每个符号采4个样点
Ns=K*N;  %总的采样点数

w=[0.5,zeros(1,N-1)];  %环路滤波器输出寄存器，初值设为0.5
n=[0.7 zeros(1,Ns-1)]; %NCO寄存器，初值设为0.9
n_temp=[n(1),zeros(1,Ns-1)];
u=[0.6,zeros(1,2*N-1)];%NCO输出的定时分数间隔寄存器，初值设为0.6
yI=zeros(1,2*N);       %I路内插后的输出数据
yQ=zeros(1,2*N);       %Q路内插后的输出数据
time_error=zeros(1,N); %Gardner提取的时钟误差寄存器

i=1;    %用来表示Ts的时间序号,指示n,n_temp,nco,
k=1;    %用来表示Ti时间序号,指示u,yI,yQ
ms=1;   %用来指示T的时间序号,用来指示a,b以及w
strobe=zeros(1,Ns);
c1=5.41*10^(-3);   c2=3.82*10^(-6);  %环路滤波器系数

%############################
%这段代码是原程序中仿真输入的psk数据
bitstream=randint(1,N,2);              
psk2=pskmod(bitstream,2);
xI=zeros(1,Ns);
xQ=zeros(1,Ns);
xI(1:8:8*N)=real(psk2);%8倍插值？为何要这样设计
xQ(1:8:8*N)=imag(psk2);
%截短后的根升余弦匹配滤波器
h1=rcosfir(0.5,[-8,8],4,1,'sqrt');
hw=kaiser(65,3.97);
hh=h1.*hw.';aI1=conv(xI,h1);
bQ1=conv(xQ,h1);
L=length(aI1);
%仿真输入数据
aI=[aI1(26:2),0,0];%2倍抽取？为何先8倍插值，再2倍抽取？
bQ=[bQ1(26:2),0,0];
%##############################


%下面是我编写的数据产生方法，但产生的数据，程序运行后不收敛？
%程序包下载后，直接运行，就可以查看两种数据模式下的仿真曲线。
N=6000;  %符号数
K=4;     %每个符号采4个样点
Ns=K*N;  %总的采样点数
w=[0.5,zeros(1,N-1)];  %环路滤波器输出寄存器，初值设为0.5
n=[0.7 zeros(1,Ns-1)]; %NCO寄存器，初值设为0.9
n_temp=[n(1),zeros(1,Ns-1)];
u=[0.6,zeros(1,2*N-1)];%NCO输出的定时分数间隔寄存器，初值设为0.6
yI=zeros(1,2*N);       %I路内插后的输出数据
yQ=zeros(1,2*N);       %Q路内插后的输出数据
time_error=zeros(1,N); %Gardner提取的时钟误差寄存器
i=1;    %用来表示Ts的时间序号,指示n,n_temp,nco,
k=1;    %用来表示Ti时间序号,指示u,yI,yQ
ms=1;   %用来指示T的时间序号,用来指示a,b以及w
strobe=zeros(1,Ns);

bitstream=randint(1,N,2);              
psk2=pskmod(bitstream,2);
%直接对psk调制信号进行4倍插值滤波后得到仿真数据
aI=rcosflt(real(psk2),1,4,'sqrt',0.5);
bQ=rcosflt(imag(psk2),1,4,'sqrt',0.5);
%采用这种数据产生方式，从仿真波形看为何不收敛？

cameion

干脆直接将程序代码及仿真波形贴出来吧。请高手诊断。谢谢。
%定时恢复程序：
clear all
N=6000;  %符号数
K=4;     %每个符号采4个样点
Ns=K*N;  %总的采样点数

w=[0.5,zeros(1,N-1)];  %环路滤波器输出寄存器，初值设为0.5
n=[0.7 zeros(1,Ns-1)]; %NCO寄存器，初值设为0.9
n_temp=[n(1),zeros(1,Ns-1)];
u=[0.6,zeros(1,2*N-1)];%NCO输出的定时分数间隔寄存器，初值设为0.6
yI=zeros(1,2*N);       %I路内插后的输出数据
yQ=zeros(1,2*N);       %Q路内插后的输出数据
time_error=zeros(1,N); %Gardner提取的时钟误差寄存器

i=1;    %用来表示Ts的时间序号,指示n,n_temp,nco,
k=1;    %用来表示Ti时间序号,指示u,yI,yQ
ms=1;   %用来指示T的时间序号,用来指示a,b以及w
strobe=zeros(1,Ns);
c1=5.41*10^(-3);   c2=3.82*10^(-6);  %环路滤波器系数

%############################
%这段代码是原程序中仿真输入的psk数据
bitstream=randint(1,N,2);              
psk2=pskmod(bitstream,2);
xI=zeros(1,Ns);
xQ=zeros(1,Ns);
xI(1:8:8*N)=real(psk2);%8倍插值？为何要这样设计
xQ(1:8:8*N)=imag(psk2);
%截短后的根升余弦匹配滤波器
h1=rcosfir(0.5,[-8,8],4,1,'sqrt');
hw=kaiser(65,3.97);
hh=h1.*hw.';aI1=conv(xI,h1);
bQ1=conv(xQ,h1);
L=length(aI1);
%仿真输入数据
aI=[aI1(26:2),0,0];%2倍抽取？为何先8倍插值，再2倍抽取？
bQ=[bQ1(26:2),0,0];
%##############################

ns=length(aI)-2;

while(i<ns)
    n_temp(i+1)=n(i)-w(ms);
    if(n_temp(i+1)>0)
        n(i+1)=n_temp(i+1);
    else
        n(i+1)=mod(n_temp(i+1),1);
        %内插滤波器模块
        FI1=0.5*aI(i+2)-0.5*aI(i+1)-0.5*aI(i)+0.5*aI(i-1);
        FI2=1.5*aI(i+1)-0.5*aI(i+2)-0.5*aI(i)-0.5*aI(i-1);
        FI3=aI(i);
        yI(k)=(FI1*u(k)+FI2)*u(k)+FI3;
        FQ1=0.5*bQ(i+2)-0.5*bQ(i+1)-0.5*bQ(i)+0.5*bQ(i-1);
        FQ2=1.5*bQ(i+1)-0.5*bQ(i+2)-0.5*bQ(i)-0.5*bQ(i-1);
        FQ3=bQ(i);
        yQ(k)=(FQ1*u(k)+FQ2)*u(k)+FQ3;
        strobe(k)=mod(k,2);
        %时钟误差提取模块，采用的是GARDNER算法
        if(strobe(k)==0)
            %每个数据符号计算一次时钟误差
            if(k>2)
               time_error(ms)=yI(k-1)*(yI(k)-yI(k-2))+yQ(k-1)*(yQ(k)-yQ(k-2));
            else
                time_error(ms)=(yI(k-1)*yI(k)+yQ(k-1)*yQ(k));
            end
            %环路滤波器.每个数据符号计算一次环路滤波器输出
            if(ms>1)
                w(ms+1)=w(ms)+c1*(time_error(ms)-time_error(ms-1))+c2*time_error(ms);
            else
                w(ms+1)=w(ms)+c1*time_error(ms)+c2*time_error(ms);
            end
            ms=ms+1;
        end
        k=k+1;
        u(k)=n(i)/w(ms);
    end
    i=i+1;
end

figure(1);
subplot(311);plot(u);
subplot(312);plot(time_error);
subplot(313);plot(w);


%##############################
% bt=0.01;
% c1=8/3*bt;
% c2=c1*4/3*bt;
N=6000;  %符号数
K=4;     %每个符号采4个样点
Ns=K*N;  %总的采样点数
w=[0.5,zeros(1,N-1)];  %环路滤波器输出寄存器，初值设为0.5
n=[0.7 zeros(1,Ns-1)]; %NCO寄存器，初值设为0.9
n_temp=[n(1),zeros(1,Ns-1)];
u=[0.6,zeros(1,2*N-1)];%NCO输出的定时分数间隔寄存器，初值设为0.6
yI=zeros(1,2*N);       %I路内插后的输出数据
yQ=zeros(1,2*N);       %Q路内插后的输出数据
time_error=zeros(1,N); %Gardner提取的时钟误差寄存器
i=1;    %用来表示Ts的时间序号,指示n,n_temp,nco,
k=1;    %用来表示Ti时间序号,指示u,yI,yQ
ms=1;   %用来指示T的时间序号,用来指示a,b以及w
strobe=zeros(1,Ns);

bitstream=randint(1,N,2);              
psk2=pskmod(bitstream,2);
%直接对psk调制信号进行4倍插值滤波后得到仿真数据
aI=rcosflt(real(psk2),1,4,'sqrt',0.5);
bQ=rcosflt(imag(psk2),1,4,'sqrt',0.5);
%采用这种数据产生方式，从仿真波形看为何不收敛？

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ns=length(aI)-16;
while(i<ns)
    n_temp(i+1)=n(i)-w(ms);
    if(n_temp(i+1)>0)
        n(i+1)=n_temp(i+1);
    else
        n(i+1)=mod(n_temp(i+1),1);
        %内插滤波器模块
        FI1=0.5*aI(i+2)-0.5*aI(i+1)-0.5*aI(i)+0.5*aI(i-1);
        FI2=1.5*aI(i+1)-0.5*aI(i+2)-0.5*aI(i)-0.5*aI(i-1);
        FI3=aI(i);
        yI(k)=(FI1*u(k)+FI2)*u(k)+FI3;
        FQ1=0.5*bQ(i+2)-0.5*bQ(i+1)-0.5*bQ(i)+0.5*bQ(i-1);
        FQ2=1.5*bQ(i+1)-0.5*bQ(i+2)-0.5*bQ(i)-0.5*bQ(i-1);
        FQ3=bQ(i);
        yQ(k)=(FQ1*u(k)+FQ2)*u(k)+FQ3;
        strobe(k)=mod(k,2);
        %时钟误差提取模块，采用的是GARDNER算法
        if(strobe(k)==0)
            %每个数据符号计算一次时钟误差
            if(k>2)
               time_error(ms)=yI(k-1)*(yI(k)-yI(k-2))+yQ(k-1)*(yQ(k)-yQ(k-2));
            else
                time_error(ms)=(yI(k-1)*yI(k)+yQ(k-1)*yQ(k));
            end
            %环路滤波器.每个数据符号计算一次环路滤波器输出
            if(ms>1)
                w(ms+1)=w(ms)+c1*(time_error(ms)-time_error(ms-1))+c2*time_error(ms);
            else
                w(ms+1)=w(ms)+c1*time_error(ms)+c2*time_error(ms);
            end
            ms=ms+1;
        end
        k=k+1;
        u(k)=n(i)/w(ms);
    end
    i=i+1;
end

figure(2);
subplot(311);plot(u);
subplot(312);plot(time_error);
subplot(313);plot(w);

cameion

swas

既然你诚心贴出来，我就帮你一次吧。

你只要在两条语句中间插一条语句就对了。
psk2=pskmod(bitstream,2);  % 你的原程序

psk2 = upsample(psk2,4);   % 就插这么一句，下面还是你的程序，
%直接对psk调制信号进行4倍插值滤波后得到仿真数据
aI=rcosflt(real(psk2),1,4,'sqrt',0.5);

cameion

我添加了psk2 = upsample(psk2,4);后，运行结果如下：

从运行结果看，u值仍然在0与1之间振荡。这是为何？
我查了下matlab的帮助，upsample函数的目的在对输入数据内插0值。upsample(psk,4)相当于在每个psk数据之间内插3个0值。如果紧接其后，采用低通滤波器滤波，相当于对对psk2信号的采样频率提高4倍。而rcosflt(real(psk2),1,4,'sqrt',0.5)本身就具备内插4倍且完成低通成形滤波的功能。因此，在程序中添加upsample(psk,4)后，相当于在原数据基础上进行了4倍内插，再加上rcosflt(real(psk2),1,4,'sqrt',0.5)，就成了16倍内插了。

还请指正！谢谢。

cameion

添加  psk2 = upsample(psk2,4);   后，程序运行结果。

swas

实际上你也可以2倍上采样后，再下采样，你取奇数点和偶数点，结果是不同的，包括原来的程序，如果取奇数点，则是你的效果，如果取偶数点，则是原图效果，你比较一下。另外原来的程序也是有问题的，如h1=rcosfir(0.5,[-8,8],4,1,'sqrt');实际上应为：h1=rcosfir(0.5,[-8,8],8,1,'sqrt');

cameion

“实际上你也可以2倍上采样后，再下采样，你取奇数点和偶数点，结果是不同的”这句话提醒了我！1.程序中的gardner位定时算法是正确的；
2.采样时，不同的起始点采样，u值的最后收敛值是不同的，但都在0~1范围内。
3.当每个符号采样4个点时，程序中我编写的代码aI=rcosflt(real(psk2),1,4,'sqrt',0.5)，表示对一个符号从起始点到最末点的时间段时每隔T/4时间采样一次，4个采样点时刻依次为0,T/4,2T/4,3T/4;这样，gardner算法收敛后，mk为2T/4时，uk=0,计算出的插值为2T/4处的值（眼图张开最大时的值），mk为T/4,uk=1，计算出的插值仍为2T/4处的值（眼图张开最大时的值）。因此，出现算法收敛后uk在0与1来回振荡的现象；
4. 当进行上采样，再进行下采样，取不同的下采样序号，根据gardner算法原理，只会出一固定的值，因此不会出现振荡现象。
5.“”另外原来的程序也是有问题的，如h1=rcosfir(0.5,[-8,8],4,1,'sqrt');实际上应为：h1=rcosfir(0.5,[-8,8],8,1,'sqrt');“ 的说法是正确的；
6. “psk2 = upsample(psk2,4);   % 就插这么一句，下面还是你的程序”只加这么一句是不可以的，后续可以再进行下采样，则可以得到正确的结果。

谢谢！非常感谢！