mlse_1h_cpm.m

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

% MLSE CPM demodulation

clear;

% *******************
% 参数控制
% *******************
% modulator index array
h=5/16;
M=2;

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

fc=12*(10^3);
fsym=16*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:100);
    D=[3 1 -1 -3];
else
    sym_seq=sym_seq_2ary(1:100);
    D=[1 -1];
end     
    
% step-test VB
%sym_seq=repmat([1 -1  1 -1],1,25);

% CPM phase trellis table

% state_pre[32*4]:Table with modulaton index h2=5/16
% for every state:symbol value with branch is 
%                 [ +3 +1 -1 -3]
   
state_pre_tab = ...
[   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];

state_pre_tab2 = ...
[   28     6; ...
    29     7; ...
    30     8; ...
    31     9; ...
    32    10; ...
     1    11; ...
     2    12; ...
     3    13; ...
     4    14; ...
     5    15; ...
     6    16; ...
     7    17; ...
     8    18; ...
     9    19; ...
    10    20; ...
    11    21; ...
    12    22; ...
    13    23; ...
    14    24; ...
    15    25; ...
    16    26; ...
    17    27; ...
    18    28; ...
    19    29; ...
    20    30; ...
    21    31; ...
    22    32; ...
    23     1; ...
    24     2; ...
    25     3; ...
    26     4; ...
    27     5];
    
% ****************************
% 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);
   for i=1:num_fs
   fai=rem(theta+2*pi*In*h*((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=rem(theta+pi*h*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
% *****************************

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=11;
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*h*bn*(t-n*T)/T),sin(pi*h*bn*(t-n*T)/T)} 
%              bn=+3 +1 -1 -3
% ***********************

% symbol index
  n=floor((k-1)/num_fs)+1;
   
% h determine
  angle_sym=rem(2*pi*h*((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
% *************************

% *************************    
% search minimum value inter subset for every state
% i:index of state
    for i=1:STATE_NUM
   
       % **************
       % Branch metric
       prestate=state_pre_tab2(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,128)==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;