Gardner定时同步的matlab实现

网络上似乎没有直接开源的代码，并且开源的代码的实现效果也不好（没有画BER-EbN0曲线）
这里把我的放给大家，不知道为什么没达到理想情况，欢迎交流改进。

% BER_timing_syn_qyb
% 需要提升采样频率，4倍-》8倍？
clear
close all;
Norder = 4;
qpskmod = comm.PSKModulator('ModulationOrder',Norder,'PhaseOffset',pi/4,'BitInput',true);
qpskdemod = comm.PSKDemodulator('ModulationOrder',Norder,'PhaseOffset',pi/4,'BitOutput',true);
Ebn0_max = 12;
% data_length
sps = 16;
Data_Len = 2000000;
data = randi([0 1],Data_Len,1);
N_ch = 10000;


% waveform = bpskmodulator(data);
% qpsk调制
data_map = qpskmod(data);
% plot(data_map,'r*')
% check = qpskdemod(data_map);
% plot(real(data_map));
% 信号上采样
data_upsp = upsample(data_map,sps);  

%  t：transmitter 发射端滤波器
span_tx = 16;
RolloffCoef = 0.22;         % Roll Off Coeff.   α设定
% sps = 4;
h_tx = rcosdesign(RolloffCoef,span_tx,sps,'sqrt');
%  r：receiver 接收端滤波器
span_rx = 16;
RolloffCoef = 0.22;         % Roll Off Coeff.   α设定

h_rx = rcosdesign(RolloffCoef,span_rx,sps,'sqrt');

% 对信号进行卷积计算， 此操作名为成型，经过滤波处理掉上采样带来的误差
data_shape1 = conv(data_upsp,h_tx,'same');
% figure;
% plot(real(data_shape1(1:200)));

% ++++++++++++++++++++++++++++++++++++++++++++
% Ebn0 = 10;
Norder = 4;

r_Sample_data = zeros(1,Data_Len);
Rec_out = r_Sample_data;
% filter
% C_2 = inline('1/6*u^3-1/6*u');
% C_1 = inline('-0.5*u^3+0.5*u^2+1*u');
% C0  = inline('0.5*u^3 - 1*u^2 -0.5*u+1');
% C1  = inline('-1/6*u^3+0.5*u^2-1/3*u');
C_2 = @(u) 1/6*u^3-1/6*u;
C_1 = @(u) -0.5*u^3+0.5*u^2+1*u;
C0  = @(u) 0.5*u^3 - 1*u^2 -0.5*u+1;
C1  = @(u) -1/6*u^3+0.5*u^2-1/3*u;

ratio_err = zeros(1,Ebn0_max+1);
for Ebn0 = 0:Ebn0_max
    fprintf('Computing Ebn0 = %d\n',Ebn0)
    SNR = Ebn0 + 10*log10(2) - 10*log10(sps);                   % S/N Ratio for Chanel Simulatioin 
%     E = mean(abs(data_shape1).^2)/log2(Norder);
%     SNR = E/(10^(Ebn0/10)) * sps;
    T_Gain = 0.0005;              % TED Loop Gain

    % 初始化参数值
    timing_error = 0;
    displaysize = Data_Len;
    
    %  ++++++++++channel-------------
%     data_received = awgn(data_shape1,SNR);
    E = mean(abs(data_shape1).^2)/log2(Norder);
    snr = E/(10^(Ebn0/10)) * sps;
    noise = sqrt(snr/2)*(randn(length(data_shape1),1)+1i*randn(length(data_shape1),1));
    data_received = data_shape1 + noise;
    data_r_shape = conv(data_received,h_rx,'same');
%     plot(real(data_r_shape(1:200)))
%     hold on
    for i = 1:2*Data_Len/2
%         plot(2*i-1,real(data_r_shape(2*i-1)),'r*');
%         1000个symbol，因此2000个sample √
       r_Sample_data(i) =  data_r_shape(8*i-6);
    end
    Rec_out(1:3) = r_Sample_data(1:3);
    mk = 0;
    uk = 0;
    Interp_Pos = 1;
    count = 0;
    timing_error = 0;
    fprintf('Detection & Interpolation\n')
    TE_out = zeros(1,Data_Len/2);
    for i = 3:2:Data_Len-3
        count = count + 1;
        %++++++++++++++Gardner TED+++++++++++++++++++定时误差检测器
        % temp= ( ted_data1 - ted_data3 ) * ted_data2;
        %++++++++++++++++++++++++++++++++++++++++++++
        ted_data1 = Rec_out(i-2);
        ted_data2 = Rec_out(i-1);
        ted_data3 = Rec_out(i);
        if(real(ted_data1) * real(ted_data3) < 0)||(imag(ted_data1)*imag(ted_data3) < 0)
            temp = -( ted_data1 - ted_data3 ) * conj(ted_data2);
            timing_error = real(temp) * T_Gain;
        else
            timing_error = 0;
        end
        TE_out(count) = timing_error;
        %      更新插值位置
        Interp_Pos = Interp_Pos + 2 + timing_error;
        
        % +++++++++++++Interpolation_1+++++++++++++++++++插值
        % out = C_2 * in_i(4) + C_1 * in_i(3) + C0 * in_i(2) + C1 * in_i(1);
        %+++++++++++++++++++++++++++++++++++++++++++++
        mk = floor(Interp_Pos);
        uk = Interp_Pos-mk;
        
        C_2u = C_2(uk);
        C_1u = C_1(uk);
        C0u  = C0(uk);
        C1u  = C1(uk);
        
        %     3
        data_2 = r_Sample_data(i-2);
        data_1 = r_Sample_data(i-1);
        data0 = r_Sample_data(i);
        data1 = r_Sample_data(i+1);
        
        Interp_data = C_2u * data_2 + C_1u * data_1 + C0u * data0 + C1u * data1;
        Rec_out(i) = Interp_data;
        %     4
        data_2 = r_Sample_data(i-1);
        data_1 = r_Sample_data(i);
        data0 = r_Sample_data(i+1);
        data1 = r_Sample_data(i+2);
    
        Interp_data = C_2u * data_2 + C_1u * data_1 + C0u * data0 + C1u * data1;
        Rec_out(i+1) = Interp_data;
        %     5
        data_2 = r_Sample_data(i);
        data_1 = r_Sample_data(i+1);
        data0 = r_Sample_data(i+2);
        data1 = r_Sample_data(i+3);
    
        Interp_data = C_2u * data_2 + C_1u * data_1 + C0u * data0 + C1u * data1;
        Rec_out(i+2) = Interp_data;
    end
    
    r_Sample_show = zeros(1,Data_Len/2);
    Rec_out_show = zeros(1,Data_Len/2);
    for i = 1:Data_Len/2
        r_Sample_show(i) = r_Sample_data(2*i-1);
        Rec_out_show(i) = Rec_out(2*i-1); 
    end
    Data_BER = zeros(Data_Len/2,1);
    for i = 1:Data_Len/2
        if(real(r_Sample_show(i))>0)
        if(imag(r_Sample_show(i))>0)
            Data_BER(i) = exp(0.25*pi*1j);
        else
            Data_BER(i) = exp(-0.25*pi*1j);
        end
        else
        if(imag(r_Sample_show(i))>0)
            Data_BER(i) = exp(0.75*pi*1j);
        else
            Data_BER(i) = exp(-0.75*pi*1j);
        end
        end
    end
%     data_downsp = downsample(data_r_shape,sps);  
% %     plot(data_downsp,'r+');
% %     grid on;
    data_get = qpskdemod(Data_BER);
    
    ratio_err(Ebn0 +1) = sum(data((N_ch+1):end) ~= data_get((N_ch+1):end))/(Data_Len-N_ch);
end

EbN0_dB = 0:Ebn0_max;
figure;

ber_Theory = berawgn(EbN0_dB, 'psk',4, 'nondiff');
semilogy(EbN0_dB,ratio_err,'-ro')
hold on;
semilogy(EbN0_dB,ber_Theory,'b')
title('BER')
ylabel('BER')
xlabel('Ebn0')

% figure;
% plot(TE_out(1:50000))