mlse_cpm-zl.asv

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

% MLSE CPM demodulation

% 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]

clear;

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
% ****************************

% modulator index array
h=[4,5]/16;
M=4;

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

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

state=(0:31)*pi/16;%%%%%%%%%%%%%%%原来是(0:31)*pi/32,改！！！！！

% length of wave displayed
display_sym=20;

% test sequence
sym_seq_4ary=(floor(rand(1,10000)*4)-2)*2+1;
sym_seq_2ary=round(rand(1,10000)-1);%%%%%%%%%%%%%%%%%%%%%并不是二进制的数据，是－1、0组成的数据,改?

if M==4 
    sym_seq=sym_seq_4ary(1:20);
    D=[3 1 -1 -3];%%%%%%%%%%%%%%%可能发送的数据
else
    sym_seq=sym_seq_2ary(1:20);
    D=[1 -1];
end     
    
% step-test VB
sym_seq=1.*ones(1,100)'; %%%%%%%%%%%%%%%%%%假设发送一组3组成的数据

%***********************
% 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=rem(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=rem(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
% *****************************

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=10;
STATE_NUM=32;

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

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=2*pi*fc_r*(k-1)*Ts;
   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=2*pi*hn*((k-1)*Ts-(n-1)*T)/(2*T)*D;%%%%%%%%%%匹配信号的相位偏移,前面的信号已经解调到了基带,因此匹配信号就不需要用载波了吧。改！！！
  
%   angle_fc=angle_fc.*ones(1,M);
  angle_fc=0.*ones(1,M);
  angle_sym=angle_fc+angle_sym;%所有可能的匹配相位,加上了载波相位,但是没有加初始相位值  
%   rd_fcdem_i=rd(k);
%   rd_fcdem_q=rd(k);%%%%%%%%%%改，将前面所得到的结果覆盖了
      
  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,num_fs);
      rd_mf_acc_q=zeros(1,num_fs); 
      rd_mf_acc_qi=0;
      rd_mf_acc_iq=0;

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

% **********
% 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    %%%%%%%%%%根据不同的n的值产生初始相位的表
% ***********       

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