mlse_cpm_zl.asv

cpm连续相位调制系统的matlab仿真程序

% MLSE CPM demodulation

clear;

% *******************
% 参数控制
% *******************
% modulator index array
h=[4,5]/16;
M=4;

% fc=carrier frequency
% T=symbol period
% Ts=sample period

fc=100*(10^3);
fsym=160*10^3;
num_fs=4;
fs=num_fs*fsym;
T=1/fsym;
Ts=1/fs;

state=(0:31)*pi/16;

% length of wave displayed
display_sym=20;

% test sequence
sym_seq_4ary=(floor(rand(1,10000)*4)-2)*2+1;
sym_seq_2ary=1-2*round(rand(1,10000));

if M==4 
    sym_seq=sym_seq_4ary(1:1000);
    D=[3 1 -1 -3];
else
    sym_seq=sym_seq_2ary(1:1000);
    D=[1 -1];
end     
    
% step-test VB
% sym_seq=repmat([ 3 1 -1 -3 1 -3 3 1 -1 -3],1,25);

% CPM phase trellis table

% h1_state_pre[32*4]:Table with modulaton index h1=4/16
% h2_state_pre[32*4]:Table with modulaton index h2=5/16
% for every state:symbol value with branch is 
%                 [ +3 +1 -1 -3]

h1_state_pre = ...
[   21    29     5    13; ...
    22    30     6    14; ...
    23    31     7    15; ...
    24    32     8    16; ...
    25     1     9    17; ...
    26     2    10    18; ...
    27     3    11    19; ...
    28     4    12    20; ...
    29     5    13    21; ...
    30     6    14    22; ...
    31     7    15    23; ...
    32     8    16    24; ...
     1     9    17    25; ...
     2    10    18    26; ...
     3    11    19    27; ...
     4    12    20    28; ...
     5    13    21    29; ...
     6    14    22    30; ...
     7    15    23    31; ...
     8    16    24    32; ...
     9    17    25     1; ...
    10    18    26     2; ...
    11    19    27     3; ...
    12    20    28     4; ...
    13    21    29     5; ...
    14    22    30     6; ...
    15    23    31     7; ...
    16    24    32     8; ...
    17    25     1     9; ...
    18    26     2    10; ...
    19    27     3    11; ...
    20    28     4    12 ];
    
h2_state_pre = ...
[   18    28     6    16; ...
    19    29     7    17; ...
    20    30     8    18; ...
    21    31     9    19; ...
    22    32    10    20; ...
    23     1    11    21; ...
    24     2    12    22; ...
    25     3    13    23; ...
    26     4    14    24; ...
    27     5    15    25; ...
    28     6    16    26; ...
    29     7    17    27; ...
    30     8    18    28; ...
    31     9    19    29; ...
    32    10    20    30; ...
     1    11    21    31; ...
     2    12    22    32; ...
     3    13    23     1; ...
     4    14    24     2; ...
     5    15    25     3; ...
     6    16    26     4; ...
     7    17    27     5; ...
     8    18    28     6; ...
     9    19    29     7; ...
    10    20    30     8; ...
    11    21    31     9; ...
    12    22    32    10; ...
    13    23     1    11; ...
    14    24     2    12; ...
    15    25     3    13; ...
    16    26     4    14; ...
    17    27     5    15];
    
% ****************************
% Send part:
%          CPM modulation main program
% ****************************

%***********************
% CPM modulator process
%***********************

% initialized
% k      = sample index
% n      = symbol index
% theta =start phase of each symbol
% fai   =instant phase of sample
k=1;
theta=0;
fai=0;

for n=1:length(sym_seq)
   In=sym_seq(n);
   hn=h(mod(n-1,2)+1);
   for i=1:num_fs
   fai=mod(theta+2*pi*In*hn*((k-1)*Ts-(n-1)*T)/(2*T),2*pi);
   sd(k)=cos(2*pi*fc*(k-1)*Ts+fai);
   k=k+1;    
   end    
   theta=mod(theta+pi*hn*In,2*pi);
end    

% display wave in time domain
plot(sd(1:num_fs*display_sym));
title('CPM time wave');
xlabel('n');
ylabel('s(n)');
grid on;
zoom on;

debug=1;


% *****************************
% channel
%        AWGN
% *****************************

snr=10^(snr_in_db/10);
sgma=sqrt(



rd=sd;%%%%%%%%%%%%%不加噪声,假设接受信号即为发送信号,并且载波同步
fc_r=fc;

% ***************************************************************
% CPM receive: 
% Principle:  IEEE trans. on comm., April,1987
%            "Reciver structure for Multi-h signaling formats"
% ***************************************************************
% argument control
PUN_LEN=8;
STATE_NUM=32;

rd_mf_acc_i=zeros(1,M);
rd_mf_acc_q=zeros(1,M);
rd_mf_acc_iq=zeros(1,M);
rd_mf_acc_qi=zeros(1,M);

metric_acc=zeros(1,STATE_NUM);
metric_acc1=zeros(1,STATE_NUM);
inf_bank0=zeros(STATE_NUM,PUN_LEN);
inf_bank1=zeros(STATE_NUM,PUN_LEN);

det_cntr=0;
det_flag=0;
bank_sel=0;

for k=1:length(rd)

% **********************
% Carrier demodulation
%      rd*{cos(2*pi*fc_r*t),sin(2*pi*fc_r*t)}             
% **********************
   
   angle_fc=rem(2*pi*fc_r*(k-1)*Ts,2*pi);
   rd_fcdem_i=rd(k)*cos(angle_fc);
   rd_fcdem_q=rd(k)*sin(angle_fc);
 
% ***********************
% Matched filter: 
%     *{cos(pi*hn*bn*(t-n*T)/T),sin(pi*hn*bn*(t-n*T)/T)} 
%              bn=+3 +1 -1 -3
% ***********************

% symbol index
  n=floor((k-1)/num_fs)+1;
   
% h determine
  hn=h(mod((n-1),2)+1);

  angle_sym=mod(2*pi*hn*((k-1)*Ts-(n-1)*T)/(2*T)*D,2*pi);
      
  rd_mf_acc_i=rd_mf_acc_i+rd_fcdem_i*cos(angle_sym);
  rd_mf_acc_q=rd_mf_acc_q+rd_fcdem_q*sin(angle_sym);
  rd_mf_acc_qi=rd_mf_acc_qi+rd_fcdem_q*cos(angle_sym);
  rd_mf_acc_iq=rd_mf_acc_iq+rd_fcdem_i*sin(angle_sym);

  if mod(k,num_fs)==0
      rd_mf_i=rd_mf_acc_i/num_fs;
      rd_mf_q=rd_mf_acc_q/num_fs;
      rd_mf_qi=rd_mf_acc_qi/num_fs;
      rd_mf_iq=rd_mf_acc_iq/num_fs;
      
      rd_mf_acc_i=zeros(1,M);
      rd_mf_acc_q=zeros(1,M); 
      rd_mf_acc_qi=zeros(1,M);
      rd_mf_acc_iq=zeros(1,M);

% *************************
% MLSE :Viterbi algorithm
% Branch metric:
%   cos(S(i))*rd_mf_i-cos(S(i))*rd_mf_q-sin(S(i))*rd_mf_qi-sin(S(i))*rd_mf_iq
% *************************

% **********
% determine state_pre_tab by hn

       if mod((n-1),2)==0
          state_pre_tab=h1_state_pre;
       else
          state_pre_tab=h2_state_pre;           
       end
% ***********       

% *************************    
% search minimum value inter subset for every state
% i:index of state
    for i=1:STATE_NUM
   
       % **************
       % Branch metric
       prestate=state_pre_tab(i,:);       
       angle_si=state(prestate); 
       BM=cos(angle_si).*rd_mf_i-cos(angle_si).*rd_mf_q-sin(angle_si).*rd_mf_qi-sin(angle_si).*rd_mf_iq;
       % **************       
       
       %****************
       % survival branch
       BM_up=metric_acc(prestate)+BM;
       [max_val,index]=max(BM_up);
       metric_acc1(i)=max_val;
       surv_sym=D(index);
       surv_state=prestate(index);
       %****************       
       
       % ***************
       % survval sym sequence memorize
        
        if bank_sel==0
            inf_bank1(i,2:PUN_LEN)=inf_bank0(surv_state,1:PUN_LEN-1);
            inf_bank1(i,1)=surv_sym;                                
        else
            inf_bank0(i,2:PUN_LEN)=inf_bank1(surv_state,1:PUN_LEN-1);
            inf_bank0(i,1)=surv_sym;
        end    
        % *****************
    
    end    % end "for i"

    metric_acc=metric_acc1;%状态更新    

% ********************
% bank alter
    if bank_sel==0
        bank_sel=1;
    else
        bank_sel=0;
    end
% ********************    

% *************************    
%       detect out 
%**************************
 if det_flag==1

    [max_val,max_pos]=max(metric_acc);
       
     if bank_sel==0
         det_out(n-PUN_LEN+1)=inf_bank0(max_pos,PUN_LEN);
     else
         det_out(n-PUN_LEN+1)=inf_bank1(max_pos,PUN_LEN);
     end

debug=1;
      
 elseif n>=PUN_LEN-1
        det_flag=1;%%%%%%%%%当接收到至少十个数据之后便开始判决
 end
% 

 if mod(n,64)==0
    metric_acc=metric_acc/2;
 end    
 
end        %end  "if mod(k,num_fs)==0"
 
end  % end for k

% **************************
%        error rate

err_cntr=0;
k=1;
test_len=length(sym_seq)-PUN_LEN;
for i=1:test_len
    if det_out(i)~=sym_seq(i)
        err_cntr=err_cntr+1;
        err_rec(k)=i;
        k=k+1;
    end
end

BER=err_cntr/test_len;
debug=1;