Gardner环 matlab仿真 BPSK系统 小数倍采样

内插用farrow

clc;clear all;close all;


%% Tx
% Fs = 1.7e6; 
Rs = 1e6;
N_symbol = 1e4;               % 仿真时间内的符号点数
% N_noise = randi([0 10] , 1 , 1);        % 随机噪声头/噪声尾长度               

osr = 8;
rolloff = 0.35;
% fir_hCos = rcosfir(rolloff , [-3 3] , osr , 1 , 'sqrt');
fir_hCos = rcosdesign(rolloff,6,osr,'sqrt');

data = randi([0 1] , 1 , N_symbol);
data_ini = pskmod(  data,2); % bpsk
data_upsample = upsample(data_ini , osr);
data_baseshape = conv(data_upsample , fir_hCos , 'same'); 

% 重采样
osr_new = 1.7;
data_base_resample = Farrow ( osr*Rs , osr_new*Rs , data_baseshape );

% 载波多普勒
data_send = data_base_resample;


%% Rx
EbN0 = 20;
SNR = EbN0 + 10 * log10(2)  - 10 * log10(osr_new);
err_number=zeros(1,length(SNR));
bit_err_ratio=zeros(1,length(SNR));

data_rece = awgn(data_send , SNR , 'measured');  %信道加噪

fir_hCos_new = Farrow ( osr*Rs , osr_new*Rs , fir_hCos );
data_mf = conv(data_rece , fir_hCos_new , 'same');
    
    %对比图
    % data_mf_down = downsample(data_mf,osr_new,0);
    % scatterplot(data_mf_down(1000:end))
    %% Gardner
    NCO = ones(1,length(data_mf));
    u = zeros();
    NCO_cnt = 0;
    d_out = 0;
    flag_interpolation = 0; % 每个符号进行一次误差纠正
    FTW_gardner = 1 / osr_new;% NCO每周期溢出两次，根据DDS公式计算码
    % FTW：2 * Rs = fs/2^N * FTW_gardner => FTW_gardner = 2 * 2^N * Rs / Fs = 2^N / osr  * 2
    index_out = 1;
    y_0 = 0;
    y_0p5 = 0;
    y_1 = 0;
    data_best_sample = [];
    err_detect = zeros();
    
    err_loop = zeros();
    Branch_A = 0;
    Branch_B = 0;
    
    for i = 4 : length(data_mf)-2
    
        NCO(i+1) = NCO(i) - FTW_gardner;
        
        if(NCO(i+1) <= 0)
            NCO_cnt = NCO_cnt+1;
    
            %数控振荡器寄存器值
            NCO(i+1) = mod(NCO(i+1), 1);
    
            %farrow 内插 整点
            %确定索引
            mk(NCO_cnt) = i;%基本指针   %从4开始
            uk(NCO_cnt) = NCO(i)/FTW_gardner ;
    
            D(4)  = data_mf( mk(NCO_cnt) +2 );
            D(3)  = data_mf( mk(NCO_cnt) +1 );
            D(2) = data_mf( mk(NCO_cnt)    );
            D(1) = data_mf( mk(NCO_cnt) -1 );
    
            f1 = 0.5*D(4) - 0.5*D(3) - 0.5*D(2) + 0.5*D(1);
            f2 = -0.5*D(4) + 1.5*D(3) - 0.5*D(2) - 0.5*D(1);
            f3 = D(2);
            d_out = uk(NCO_cnt)^2*f1 + uk(NCO_cnt)*f2 + f3;
    
            y_0 =  y_1;
            y_1 = d_out;
            
    
            %farrow 内插 0.5点
            %确定0.5点索引
            mk_0p5 = fix( mk(NCO_cnt) + uk(NCO_cnt) - 0.5/FTW_gardner );% Gardner I (4)   %向0取整
            uk_0p5 = mk(NCO_cnt) + uk(NCO_cnt) - 0.5/FTW_gardner - mk_0p5 ;% Garder I (5)
            
            D(4)  = data_mf( mk_0p5 +2 );
            D(3)  = data_mf( mk_0p5 +1 );
            D(2) = data_mf( mk_0p5    );
            D(1) = data_mf( mk_0p5 -1 );
    
            f1 = 0.5*D(4) - 0.5*D(3) - 0.5*D(2) + 0.5*D(1);
            f2 = -0.5*D(4) + 1.5*D(3) - 0.5*D(2) - 0.5*D(1);
            f3 = D(2);
            d_out = uk_0p5^2*f1 + uk_0p5*f2 + f3;
    
            y_0p5 = d_out;
            %%
            data_best_sample(NCO_cnt) = y_1;
            % 误差检测
            err_detect(NCO_cnt) = real(y_0p5) * (real(y_1) - real(y_0)) + imag(y_0p5) * (imag(y_1) - imag(y_0));
    
            % 环路滤波
            c1 = 1 / 2^5;
            c2 = 1 / 2^12;
            Branch_A = c1 * err_detect(NCO_cnt);
            Branch_B = Branch_B + c2 * err_detect(NCO_cnt);
            err_loop(NCO_cnt) = Branch_A + Branch_B;
    
            NCO(i+1) = NCO(i+1) - err_loop(NCO_cnt);

        end 
    end
    

figure;
plot(err_detect / 2^32);
title('Gardner环路误差')

figure;
plot(NCO ,'r-*');
title('NCO寄存器变化过程')

figure;
plot(uk/2^32);
title('miu')

scatterplot(data_best_sample(1000:2:end))

figure;
plot(real(data_best_sample),'b*');hold on;plot(imag(data_best_sample),'r*');

figure;
theoryBer = berawgn(EbN0,'psk',2,'nondiff');
figure;
subplot(2,2,1);
plot(EbN0,10*log10(bit_err_ratio),'-*');
hold on;
plot(EbN0,10*log10(theoryBer),'-+');
xlabel('Eb/N0');ylabel('误码率');
title('不同信噪比下误码率仿真曲线');
legend('实验曲线','理论曲线');
%%
function data_farrow = Farrow ( Fs_in , Fs_out , data_in )
    i=1;
    while( i * Fs_in / Fs_out <= length ( data_in ))
        m_k = floor( i * Fs_in / Fs_out );
        mu_k = i * Fs_in / Fs_out - m_k;
        if ( m_k > 3)
            x_1 = 0.5 * data_in (m_k) - 0.5 * data_in (m_k-1) - 0.5 * data_in (m_k-2) + 0.5 * data_in (m_k-3) ;
            x_2 = 1.5 * data_in (m_k-1) - 0.5 * data_in (m_k) - 0.5 * data_in (m_k-2)- 0.5 * data_in (m_k-3) ;
            x_3 = data_in (m_k-2);
            data_farrow(i) = x_1 * ( mu_k^2 ) + x_2 * mu_k + x_3;
        else
            data_farrow(i) = 0;
        end
        i=i+1;
    end
end