目录

1.算法描述

2.仿真效果预览

3.MATLAB核心程序

4.完整MATLAB

1.算法描述

        随着数字通信技术的飞速发展，数字通信的调制技术也得到了快速发展，其中连续相位调制(Continuous Phase Modulation，CPM)作为一种全新的通信调制方法得到了越来越多的关注。CPM和传统调制方法不同之处在于，其具有带外辐射小，抗干扰能力强以及频率利用效率高等优势。这是由于传统的FSK和PSK等调制方式，其存在调制相位不连续的情况，而这种不连续的现象会导致相位突变问题，即在信号的传输过程中，信号会从一个相位特性突然切换到另外一个相位特性使得信号在频域存在较大的旁瓣，那么此时需要一个较大的频带来进行传输，导致频谱利用率低下，而CPM调制方式则可以有效避免这种现象发生。CPM是一种通过恒包络调制，从而降低了CPM信号的发送功率。       

       CPM(连续相位调制)是近年来发展很快的一种高效的调制技术。与其它调制技术相比,它具有较高的频谱利用率和功率效率,这也令它在通信资源日益紧张的今天得到了越来越多的关注。CPM信号包含大量的信号形式,其共同的特点是信号包络恒定、相位连续。CPM调制体系是一种信道编码和调制相结合的方式，通过生成相位状态格子序列控制下一时刻的状态转移，因此信息符号的调制直接具有编码效应，从而无须更多的冗余符号.CPM调制体系可以通过最大似然序列估计（MLSE）对接收信号进行软判决，实现低信噪比条件下的高检测概率.这些优良特性使得CPM技术在幅度衰落信道中实现高数据传输速率的通信系统中具有广阔的应用前景.

        CPM即使连续相位调制(Continue Phase Modulation），是一种相位调制技术，它具有相位连续的特点，频谱特性优良，相比PSK调制方式，具有更高的频带利用率。并且，CPM调制体系是一种信道编码和调制相结合的方式，通过生成相位状态格子序列控制下一时刻的状态转移，因此信息符号的调制直接具有编码效应，从而无须更多的冗余符号.CPM调制体系可以通过最大似然序列估计（MLSE）对接收信号进行软判决，实现低信噪比条件下的高检测概率.这些优良特性使得CPM技术在幅度衰落信道中实现高数据传输速率的通信系统中具有广阔的应用前景.

      CPM基带信号可以表示为
      

其瞬时相位信息可以表示为：

CPM信号最大的优势在于其频谱利用率较高，当传输的符号为等概率情况下，CPM信号的功率谱表达式为 ：

CPM解调系统如下： 

        在CPM解调过程中，接收信号通过载波同步模块，得到解调之后的基带信号，然后通过低通滤波器滤除存在的高频信息，然后将得到的基带信号进行定时同步操作得到定时同步的基带数据，从而完成CPM解调。  目前, CPM信号的载波同步方法主要分为数据辅助同步和非数据辅助同步。数据辅助同步方法的同步性能较好, 但需占用额外的传输带宽和功率; 非数据辅助同步方法不需额外的带宽和功率, 但同步精度较差, 且仅适用于信噪比较高的情况, 很难实现有效同步。

    传统的CPM信号使用Viterbi算法实现相干解调, 不便于软信息的提取, 为此, 采用MAP算法对CPM信号进行软解调, 以获取载波相位同步所需的软解调信息。

由Rimoldi分解可知, CPM信号隐含一个具有记忆性的网格编码器, 类比卷积码的网格图, 可定义 

时刻CPM信号的状态为: 时刻的状态即为: 相邻时刻两个相位状态之间的关系为:

CPM信号的网格图如图1所示, 其中 

和 分别表示起始状态和结束状态, 和 分别表示输入信息符号和输出码字符号, 表示状态集合, 表示输入符号集合, 表示任意两个网格状态之间的所有可能状态转移。

2.仿真效果预览

matlab2022a仿真结果如下：

3.MATLAB核心程序

。。。。。。。。。。。。。。。。。。。。。。。。。
 
OVS=8;
N_fft=OVS*2*4;   
M_fft=OVS*2;   
L_fft=N_fft/M_fft;
 
N_data=100;
preamble=[1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0];
tao0=0.2*T/OVS;
 
t=[0:T/OVS:L*T];
gt=( 1-cos( 2*pi*(t-tao0)/(L*T) ) )/(2*L*T);
gt=gt(1:(length(gt)-1));
 
 
databit=rand(1,N_data)>0.5;
modbit=[preamble,databit];
length_modbit=length(modbit);
length_modsignal=(length_modbit+2)*OVS;
 
mod_gt=zeros(1,length_modsignal);
for ii=1:length_modbit
    if modbit(ii)==1    %1is modulated to -1
        mod_gt((ii-1)*OVS+1:(ii-1)*OVS+L*OVS)=mod_gt((ii-1)*OVS+1:(ii-1)*OVS+L*OVS)-gt;
    else
        mod_gt((ii-1)*OVS+1:(ii-1)*OVS+L*OVS)=mod_gt((ii-1)*OVS+1:(ii-1)*OVS+L*OVS)+gt;
    end
end 
 
phai=zeros(1,length_modsignal);
for ii=1:length_modsignal
 
 phai(ii)=mod(2*pi*h*sum(mod_gt(1:ii)),2*pi);
 
end 
modsignal=cos(phai);
 
I=cos(phai);
Q=sin(phai);
 
%假设已经鉴频完毕
freq_d=[zeros(1,16*OVS),mod_gt+2000*2*pi];
freq_demod=awgn(freq_d,10);   
 
 
tmp_fft_L=zeros(1,length(freq_demod)-N_fft+1); %tone component
tmp_deltaf=zeros(1,length(freq_demod)-N_fft+1);%carrier offset
tmp_symbole=zeros(1,length(freq_demod)-N_fft+1);%symbol timing error
for ii=N_fft:(length(freq_demod))
    tmp_deltaf(ii-N_fft+1)=sum(freq_demod(ii-N_fft+1:ii))/(2*pi*N_fft);
    nn=0;
    for jj=ii-N_fft+1:ii
        tmp_fft_L(ii-N_fft+1)=tmp_fft_L(ii-N_fft+1)+freq_demod(jj)*exp(-i*(2*pi/M_fft)*nn);%第L个
        nn=nn+1;
    end
    tmp_symbole(ii-N_fft+1)=-((T/ts)/pi)*angle(tmp_fft_L(ii-N_fft+1));
end 
 
tmp_symbole=abs(tmp_symbole);
 
 
tao=tmp_symbole(150:278);
  for i=1:length(tao)/(2*OVS)
             for k=1:2*OVS
                 mtao(i,k)=tao((i-1)*2*OVS+k);
             end
  end
            ktao=min(mtao,[],2);
            ktao=ktao';
mse=0;
for ii=1:length(ktao)
        mse=mse+ (  ktao(ii)-tao0*OVS/T )^2;   
end
01_173m

EbN0=[5:2.5:30]; 
     for yy=1:length(tao0)
%         for yy=1:1
        t=[0:T/CY:L*T];
        gt=1/(2*L*T).*ones(1,length(t));
        qt=(  t-tao0(yy) )/(2*L*T); 
%     gt=(1-cos(2*pi*(t)/(L*T)))/(2*L*T);%gt脉冲函数LRC
%     qt=-1/4*(-2*pi*(t-tao0(yy))+L*T*sin(2*pi*(t-tao0(yy))/L/T))/L/T/pi;

        pha=zeros(1,length(data));
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% modulation
         for i=L:length(data)
        
                k=i-L+1;
                if(k>1)
                pha(k)=mod(pha(k-1)+pi*h*data(k-1),2*pi);
                end
                pha_t((k-1)*CY+1:k*CY)=mod(2*pi*h* (data(i)*qt(1:CY))+pha(k) ,2*pi);
                %s((k-1)*CY+1:k*CY)=cos(pha_t((k-1)*CY+1:k*CY));
        end
% plot(pha)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        mo_exp=exp(pha_t*j);

            for tt=1:length(EbN0)
%          for tt=1:1
                %snr(tt) = EbN0(tt)- 10 * log10(CY) + 10 * log10(log2(m));
                mo_exp_res=awgn(mo_exp,EbN0(tt));

%                     mo_exp_res=mo_exp;
                cnt=1;
                ll=length(mo_exp_res)/CY;
                x_multi=zeros(1,CY,ll);
                x_multi_temp=zeros(1,CY);
                 all_ee_sum=zeros(1,CY);
                  ee_sum=zeros(1,CY);
              for i=1:m  
                  cnt=i+1;
                   x_multi=zeros(1,CY,ll);
                    for k=(i+1):ll
                        xk1=mo_exp_res( (k-1-i)*CY+1:(k-i)*CY);
                        xk2=mo_exp_res( (k-1)*CY+1:(k)*CY);
                        x_multi(1,:,cnt)=( conj(xk1).*xk2 ).^2/ll;
                        cnt=cnt+1;
                    end
                   ee_sum=zeros(1,CY);
                    for a=1:CY
                        for it=1:cnt-1
                            ee_sum(a)=ee_sum(a)+x_multi(1,a,it);
                            ee_sum(a)=ee_sum(a);
                        end
                    end
                    abs(ee_sum);
                all_ee_sum=all_ee_sum+abs(ee_sum);

              end
                
                rm_sum=0;
                for a=1:CY
                       rm_sum=rm_sum+abs(all_ee_sum(a))*exp(-j*2*pi*(a-1)/CY)*1/4;
                       rm_sum=rm_sum;
                end
                tao(yy,tt)=-T*angle(rm_sum)/(2*pi);
            end


    end

4.完整MATLAB

V