map_cpm.m

连续相位调制CPM的最大后验概率MAP解调

function [decoder_output]=MAP_cpm(channel_output,M,h0,h1,L,Ns,x_data,Eb_n0)
tic;
%   仍然是4倍采样，L=3,h=1/2;是K-lag解调，可以假设K=3或4。
%   分几个步骤：设初值、第一部分进行迭代，第二部分进行迭代、相加进行判决
%   首先所求得概率值，注意不能超过单位1，这是最基本的。
%   首先分析状态网格，进行基本的运算。
% L=3;h=1/2;Ns=4;M=2;number_of_states=16;
h=h0/h1;
if rem(h0,2)==1
    theta_of_states=2*h1;
else
    theta_of_states=h1;
end
number_of_states=M^(L-1)*theta_of_states;
for j=0:number_of_states-1
  for l=0:M-1     %   考虑M进制
    [next_state,branch_output]=nxt_stat_cpm(j,l,M,h0,h1,L,Ns);     %   四个参数:状态序号和输入值，后面两个是状态参数  
    %   已验证 input 矩阵是对的。
    input(j+1,next_state+1)=l;
    %   已验证 nextstate 矩阵是对的。
    nextstate(j+1,l+1)=next_state;
    output(j+1,l+1,1:Ns)=branch_output;
  end
end
%   网格深度（数据长度）
depth_of_trellis=length(channel_output)/Ns;  
%   信道输出矩阵，即译码输入值
channel_output_matrix=reshape(channel_output,Ns,depth_of_trellis);   

%   初始化两个矩阵SkAdd1Yk,DkSkYk，这两个矩阵必须同时更新和计算。
%   最好编写两个函数，这样在循环中不停的调用,实现第一部迭代
SkAdd1Yk=zeros(number_of_states,depth_of_trellis+1);      %16*1001
SkAdd1Yk(:,1)=ones(number_of_states,1)/number_of_states;
DkSkYk=zeros(number_of_states,M,depth_of_trellis);   %   二进制,16*2*1000
for i=1:depth_of_trellis
    [temp_DkSkYk,DkSkYk(:,:,i)]=DkSk_Yk(SkAdd1Yk(:,i),channel_output_matrix(:,i),output,Eb_n0);
    %   中间还有一些处理。
    SkAdd1Yk(:,i+1)=Skadd1_Yk(DkSkYk(:,:,i),nextstate);
end
%   下面实现第二部迭代，首先是DkYkSkAdd1矩阵，为16*2*1000，为下一步迭代计算作准备。
DkYkSkAdd1=Dk_YkSkAdd1(DkSkYk,SkAdd1Yk,nextstate);
%   只需加3次，因为第一次是条件独立的。
%   我们选取K-lag的K等于4，有待验证。
%   计算p(Dk-4|Yk,Dk,Sk+1)，因为条件独立p(Dk-4|Yk,Dk,Sk+1)=p(Dk-4|Yk-1,Sk)，
%   我们认为p(Dk-4|Yk-1,Sk)是DkSub3YkSkAdd1，便于表达
%   生成DkSub3YkSkAdd1为16*2*996矩阵（p(d1|y4,s5)——p(d996|y999,s1000)）
%   中间会有DkSub2YkSkAdd1，DkSub1YkSkAdd1矩阵，为16*2
%   反正加了3次，上面3个矩阵表达上是反过来的。
%   好像也不等于4，尽管连续4个1会使得状态达到一致，但是好像是6，得经过6个才能达到状态一致，产生最小欧氏距离。
%   根据DkYkSkAdd1这个矩阵一步步计算得到所需要的后验判决值。
%   没有噪声，性能太好，将NaN处理掉，要处理的函数有：DkYkSkAdd1和DkSub3YkSkAdd。加噪后就不需要特殊的除法程序了。
for i=1:depth_of_trellis-4
    temp_DkSub1YkSkAdd1=zeros(number_of_states,M);
    temp_DkSub2YkSkAdd1=zeros(number_of_states,M);
    temp_DkSub3YkSkAdd1=zeros(number_of_states,M);
    for j=0:number_of_states-1  %   加第一次
        [line_temp,column_temp]=location_state(nextstate,j);
        for k=1:size(line_temp,2);
            temp_DkSub1YkSkAdd1(j+1,:)=temp_DkSub1YkSkAdd1(j+1,:)+DkSkYk(line_temp(k),column_temp(k),i+1).*DkYkSkAdd1(line_temp(k),:,i);
        end
        %   下面2行可放在for循环外。下面同
    end
    temp_SkAdd1Yk=SkAdd1Yk(:,i+2);
%     加噪后
%     temp_DkSub1YkSkAdd1=temp_DkSub1YkSkAdd1./temp_SkAdd1Yk(:,ones(1,2));
    temp_DkSub1YkSkAdd1=divider16_2(temp_DkSub1YkSkAdd1,temp_SkAdd1Yk(:,ones(1,M)));
    for j=0:number_of_states-1  %   加第二次
        [line_temp,column_temp]=location_state(nextstate,j);
        for k=1:size(line_temp,2);
            temp_DkSub2YkSkAdd1(j+1,:)=temp_DkSub2YkSkAdd1(j+1,:)+DkSkYk(line_temp(k),column_temp(k),i+2).*temp_DkSub1YkSkAdd1(line_temp(k),:);
        end
    end
    temp_SkAdd1Yk=SkAdd1Yk(:,i+3);
%     加噪后
%     temp_DkSub2YkSkAdd1=temp_DkSub2YkSkAdd1./temp_SkAdd1Yk(:,ones(1,2));
    temp_DkSub2YkSkAdd1=divider16_2(temp_DkSub2YkSkAdd1,temp_SkAdd1Yk(:,ones(1,M)));
    for j=0:number_of_states-1  %   加第三次
        [line_temp,column_temp]=location_state(nextstate,j);
        for k=1:size(line_temp,2);
            temp_DkSub3YkSkAdd1(j+1,:)=temp_DkSub3YkSkAdd1(j+1,:)+DkSkYk(line_temp(k),column_temp(k),i+3).*temp_DkSub2YkSkAdd1(line_temp(k),:);
        end
    end
    temp_SkAdd1Yk=SkAdd1Yk(:,i+4);
%     加噪后
%     temp_DkSub3YkSkAdd1=temp_DkSub3YkSkAdd1./temp_SkAdd1Yk(:,ones(1,2));
    temp_DkSub3YkSkAdd1=divider16_2(temp_DkSub3YkSkAdd1,temp_SkAdd1Yk(:,ones(1,M)));
    %   最后生成16*2*996的DkSub3YkSkAdd1矩阵
    DkSub3YkSkAdd1(:,:,i)=temp_DkSub3YkSkAdd1;
end
%   最后生成p(Dk-4|Yk)（从p(d1|y5)—P(d996|y1000)）概率矩阵DkSub4Yk（其中K-lag的K等于4），为2*996矩阵。
for i=1:depth_of_trellis-4
    temp_all_DkSub4Yk=zeros(1,M);
%     temp0_DkSub4Yk=0;temp1_DkSub4Yk=0;
    for j=1:M
        temp_DkSkYk=0;
        for k=1:M
            temp_DkSkYk=temp_DkSkYk+DkSkYk(:,k,i+4);
        end
        temp_all_DkSub4Yk(k)=sum(DkSub3YkSkAdd1(:,j,i).*temp_DkSkYk);
        DkSub4Yk(j,i)=temp_all_DkSub4Yk(k);
    end
end
%    下面开始判决，取两者概率的相对大值。解调输出
[temp_Max,I]=max(DkSub4Yk);
decoder_output=I-1;






toc;