siso_cpm.m

连续相位调制CPM的软输入软输出SISO解调

function [ decoder_output , output_MAP ]=SISO_cpm(channel_output,M,h0,h1,L,Ns,x_data,Eb_n0, aa_state)
tic;
%   仍然是4倍采样，L=3,h=1/2;是K-lag解调，可以假设K=3或4。
%   分几个步骤：设初值、第一部分进行迭代，第二部分进行迭代、相加进行判决
%   首先所求得概率值，注意不能超过单位1，这是最基本的。
%   首先分析状态网格，进行基本的运算。
% L=3;Ns=4;M=4;h0=1;h1=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);   
%   这个程序可分为三部分，前向度量、后向度量和最后一步整合。
%   afa和beta为 numannel_output_matrix(:,i))'*real(temp_output));
            temp_metric_imag=Eb_n0*2/Ns*(imag(channel_output_matrix(:,i))'*imag(temp_output));
            temp_metric=exp(temp_metric_real+temp_metric_imag);
            beta(j+1,i)=beta(j+1,i)+beta(nextstate(j+1,k)+1,i+1)*temp_metric*apP(k,i);
        end
    end
    beta(:,i)=beta(:,i)./sum(beta(:,i));
end
%   last
%   输出M*depth_of_trellis矩阵output_MAP
output_MAP=zeros(M,depth_of_trellis);
for i=2:depth_of_trellis+1
    for m=1:M
        for j=0:number_of_states-1
            for index_output=1:size(output,3)
                temp_output(index_output,1)=output(j+1,m,index_output);
            end
            temp_metric_real=Eb_n0*2/Ns*(real(channel_output_matrix(:,i-1))'*real(temp_output));
            temp_metric_imag=Eb_n0*2/Ns*(imag(channel_output_matrix(:,i-1))'*imag(temp_output));
            temp_metric=exp(temp_metric_real+temp_metric_imag)*apP(m,i-1);
            temp_metric=afa(j+1,i-1)*temp_metric*beta(nextstate(j+1,m)+1,i);
            output_MAP(m,i-1)=output_MAP(m,i-1)+temp_metric;
        end
    end
%     归一化
    output_MAP(:,i-1)=output_MAP(:,i-1)./sum(output_MAP(:,i-1));
    
    [temp,decoder_output(i-1)]=max(output_MAP(:,i-1));
end
decoder_output=decoder_output-1;



toc;