16-QAM的Gardner符号定时

        从网上找的程序修改而来，网上的程序都是4倍采样，而看Gardner的资料介绍都是2倍采样。故改了一个2倍采样的16qam定时程序，第一版。

N = 1e4;
M = 16;
data = randi([0 M-1],N,1);
sig = qammod(data,M,'UnitAveragePower',true);
rrc = rcosdesign(0.25,8,4,'sqrt');
tx = conv(upsample(sig,4),rrc,'same');
tx_off = uk_filter(tx,0.5);               % 构造0.5个符号的延迟
ch = awgn(tx_off,20,'measured');
rx = conv(ch,rrc,'same');
rx_down = resample(rx,1,2);               % 抽取到2倍符号采样,gardner每个符号要求2个采样点

w = 0.5;                                  
y = zeros(length(rx_down)/2,1);
y_cnt = 0;

strobe = 0;      % 最佳观测点指示
w3 = zeros(3,1); % 三个插值输出点
c1 = 5.41e-3;    % 来自https://blog.csdn.net/weixin_44884357/article/details/95619955
c2 = 3.82e-6;    % 来自https://blog.csdn.net/weixin_44884357/article/details/95619955

loopout = 0;     % 环路滤波器输出
loopint = 0;     % 环路滤波器积分
nco = 0;         % NCO
uk = 0;          % 分数间隔

for i=4:length(rx_down)

    y_f = uk_filter(rx_down(i-3:i),uk); % 每次将4个采样点输入插值滤波器

    w3 = [y_f(end);w3(1:2)];            % 每次将插值滤波器的输出存入缓存中,采用点间隔1/2符号,共存3个符号

                                        %                                [y(k) , y(k-1/2) , y(k-1)]
    if strobe == 1                      % 最佳观测点,输出,此时,缓存w3的内容[strobe,  error   , strobe];
        y_cnt = y_cnt + 1;
        y(y_cnt) = y_f(end);            % 输出为最佳观测点

        % uk=yi(k-1/2)*[yi(k)-yi(k-1)] + yq(k-1/2)*[yq(k)-yq(k-1)]    计算TED的公式
        ted = real(w3(2))*(real(w3(1))-real(w3(3))) + imag(w3(2))*(imag(w3(1))-imag(w3(3)));
    else
        ted = 0;                        % 没有中值点,[error,strobe,error],不计算定时误差
    end

    % 环路滤波
    loopout = ted*c1 + loopint;
    loopint = ted*c2 + loopint;

    w = loopout + 0.5;                  % 每个采样点，都将增加0.5,及每2个采样点必将产生一个最佳观测点输出
    if nco < w                          % 0<=uk<1             
        strobe = 1;
        uk = nco/w;                     % 计算最佳观测点的插值分数间隔
    else
        strobe = 0;
    end
    nco = mod(nco - w,1);               % 模1运算,0<=nco<1
end

y0 = y(500:end);

plot(real(rx_down(2:2:end)),imag(rx_down(2:2:end)),'.r');
hold on;
plot(real(y0),imag(y0),'.b');
hold off;

程序运行后的星座图如下，uk为0.7884，猜猜为什么。