BPSK调制gardner算法的MATLAB实现

什么插值滤波器、定时误差检测器，环路滤波器，定时误差控制器这里就不说了。
代码如下：

%BPSK gardner  2022年8月30日 北京理工大学

clc;clear;close all;

%参数设定
symbol_rate = 1000;     %每秒symbol_rate个符号
sps = 4;                %每个符号sps个采样点
fs = sps*symbol_rate;   %每秒就是sps*symbol_rate个采样点
fc = 4*fs;              %载波频率
P = 40000;              %P是resample函数的分子
Q = 40002;              %Q是resample函数的分母
N = 1E5;                %设置了N个符号

%信号源
msg_source = randi([0 1],1,N);          %生成1行N列范围是[0,1]的随机序列
bipolar_msg_source = 2*msg_source -1;   %将[0,1]映射为[-1,1]，便于成型卷积

%根升余弦滤波器
rolloff_factor = 0.5;                   %滚降系数
span = 6;                               %横跨6个符号的长度
rcos_fir = rcosdesign(rolloff_factor,span,sps,'sqrt');

%sps倍上采样
up4_bipolar_msg_source = upsample(bipolar_msg_source,sps);
%成型滤波
rcos_msg_source = conv(up4_bipolar_msg_source,rcos_fir,'same');
%载波调制
nn = 0:length(rcos_msg_source)-1;
rcos_msg_source_module = rcos_msg_source.*cos(2.*pi.*nn.*fc./fs);

%环路滤波器系数
k1 = 5.41*10^(-4);
k2 = 3.82*10^(-6);

EbN0 = 0:8;
SNR = EbN0 - 10*log10(sps) + 10*log10(2);
err_number=zeros(1,length(SNR));
bit_err_ratio=zeros(1,length(SNR));    
for i = 1:length(SNR)
    %信道加噪
    rcos_msg_source_addnoise = awgn(rcos_msg_source_module,SNR(i),'measured');
    %解调
    rcos_msg_source_addnoise_demodulation = rcos_msg_source_addnoise.*cos(2*pi*nn*fc/fs);
    %采样时钟出现了偏差
    data_re = resample(rcos_msg_source_addnoise_demodulation,P,Q);
    %匹配滤波
    rcos_msg_source_MF = conv(data_re,rcos_fir,'same');

    len_rcos_msg_source_MF = length(rcos_msg_source_MF);
    
    %gardner环参数设置
    w = [0.5,zeros(1,N-1)]; %环路滤波器输出寄存器
    n = [0.9,zeros(1,N-1)]; %NCO寄存器
    
    %nco输出的定时分数间隔寄存器
    u = [0.6,zeros(1,len_rcos_msg_source_MF - 1)];
    yI = zeros(1,2*N);%I路内插之后的输出数据
    time_error = zeros(1,N);%时钟误差寄存器
    
    ii = 1;%用来表示Ts的时间序号，指示n,n_temp,nco
    k = 1; %用来表示Ti的时间序号，指示u,yI,yQ
    ms = 1;%用来指示T的时间序号，用来指示a,b以及w
    
    for ii = 2:(len_rcos_msg_source_MF - 3)
        n(ii+1) = n(ii) - w(ms);
        if(n(ii+1)<0)
            u(k) = n(ii)/w(ms);
            n(ii+1) = mod(n(ii+1),1);
            
            %内插滤波器
            F1=0.5*rcos_msg_source_MF(ii+2)-0.5*rcos_msg_source_MF(ii+1)-0.5*rcos_msg_source_MF(ii)+0.5*rcos_msg_source_MF(ii-1);
            F2=1.5*rcos_msg_source_MF(ii+1)-0.5*rcos_msg_source_MF(ii+2)-0.5*rcos_msg_source_MF(ii)-0.5*rcos_msg_source_MF(ii-1);
            F3=rcos_msg_source_MF(ii);
            yI(k)=(F1*u(k)+F2)*u(k)+F3;
            
            % 时钟误差提取模块
            if(mod(k,2)==1)
                %每一个数据符号计算一次误差
                if(k>2)
                    time_error(ms)=yI(k-1)*(yI(k)-yI(k-2));
                else
                    time_error(ms)=0;
                end
                if(ms>1)
                    w(ms+1)=w(ms)+k1*(time_error(ms)-time_error(ms-1))+k2*time_error(ms);
                else
                    w(ms+1)=w(ms)+k1*time_error(ms)+k2*time_error(ms);
                end
                ms = ms +1;
            end
            k=k+1;
        end
    end
    
    %抽样判决
    yI_decision = (sign(yI(1:2:end))+1)/2;
    [err_number(i),bit_err_ratio(i)] = biterr(msg_source(1:N),yI_decision(1:N));
end

%仿真结果
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('实验曲线','理论曲线');
grid on; 
subplot(2,2,2);plot(u);xlabel('运算点数');ylabel('分数间隔');
subplot(2,2,3);plot(time_error);xlabel('运算点数');ylabel('定时误差');
subplot(2,2,4);plot(w);xlabel('运算点数');ylabel('环路滤波器输出');