viterbi_detector_v2.m

小区初搜为GSM系统中的一个关键过程

function [ rx_burst_hard , rx_burst_soft ] = viterbi_detector(SYMBOLS,NEXT,PREVIOUS,START,STOPS,Y,Rhh,Y1)
%
% VITERBI_DETECTOR:
%           This matlab code does the actual detection of the
%           received sequence. As indicated by the name the algorithm
%           is the viterbi algorithm, which is a MLSE. At this time 
%           the approch is to use Ungerboecks modified algorithm, and
%           to return hard output only.
%
% SYNTAX:   [ rx_burst ] 
%           = 
%           viterbi_detector(SYMBOLS,NEXT,PREVIOUS,START,STOPS,Y,Rhh)
%
% INPUT:    SYMBOLS:  The table of symbols corresponding the the state-
%                     numbers. Format as made by make_symbols.m
%           NEXT:     A transition table containing the next legal
%                     states, as it is generated by the code make_next.m
%           PREVIOUS: The transition table describing the legal previous
%		      states as generated by make_previous.m
%           START:    The start state of the algorithm.
%           STOPS:    The legal stop states.
%           Y:        Complex baseband representation of the matched 
%                     filtered and down converted received signal, as it
%                     is returned by mf.m
%           Rhh:      The autocorrelation as estimated by mf.m
%
% OUTPUT:   rx_burst: The most likely sequence of symbols.
%
% SUB_FUNC: make_increment
%
% WARNINGS: None.
%
% TEST(S):  Tested with no noise, perfect syncronization, and channel 
%           estimation/filtering. (Refer to viterbi_ill.m)
%
% AUTOR:    Jan H. Mikkelsen / Arne Norre Ekstr鴐
% EMAIL:    [e-mail] / [e-mail]
%
% $Id: viterbi_detector.m,v 1.7 1997/11/18 12:41:26 aneks Exp $

% % KNOWLEDGE OF Lh AND M IS NEEDED FOR THE ALGORITHM TO OPERATE

[ M , Lh ] = size(SYMBOLS);

% THE NUMBER OF STEPS IN THE VITERBI
%
STEPS=length(Y);

% INITIALIZE TABLES (THIS YIELDS A SLIGHT SPEEDUP).
%
METRIC = zeros(M,STEPS);
SURVIVOR = zeros(M,STEPS);
LLR = zeros(1,STEPS);  %  20070926

% DETERMINE PRECALCULATABLE PART OF METRIC
%
INCREMENT=make_increment(SYMBOLS,NEXT,Rhh);

% THE FIRST THING TO DO IS TO ROLL INTO THE ALGORITHM BY SPREADING OUT 
% FROM 	THE START STATE TO ALL THE LEGAL STATES. 
%
PS=START;

% NOTE THAT THE START STATE IS REFERRED TO AS STATE TO TIME 0
% AND THAT IT HAS NO METRIC.
%
S=NEXT(START,1);
METRIC(S,1)=real(conj(SYMBOLS(S,1))*Y(1))-INCREMENT(PS,S);
SURVIVOR(S,1)=START;

S=NEXT(START,2);
METRIC(S,1)=real(conj(SYMBOLS(S,1))*Y(1))-INCREMENT(PS,S);
SURVIVOR(S,1)=START;

PREVIOUS_STATES=NEXT(START,:);

% MARK THE NEXT STATES AS REAL. N.B: COMPLEX INDICATES THE POLARITY
% OF THE NEXT STATE, E.G. STATE 2 IS REAL.
%
COMPLEX=0;

for N = 2:Lh,
  if COMPLEX,
    COMPLEX=0;
  else
    COMPLEX=1;
  end
  STATE_CNTR=0;
  for PS = PREVIOUS_STATES,
    STATE_CNTR=STATE_CNTR+1;
    S=NEXT(PS,1);
    METRIC(S,N)=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N))-INCREMENT(PS,S);
    SURVIVOR(S,N)=PS;
    USED(STATE_CNTR)=S;
    STATE_CNTR=STATE_CNTR+1;
    S=NEXT(PS,2);
    METRIC(S,N)=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N))-INCREMENT(PS,S);
    SURVIVOR(S,N)=PS;    
    USED(STATE_CNTR)=S;
  end
  PREVIOUS_STATES=USED;
end
% AT ANY RATE WE WILL HAVE PROCESSED Lh STATES AT THIS TIME
%
PROCESSED=Lh;

% WE WANT AN EQUAL NUMBER OF STATES TO BE REMAINING. THE NEXT LINES ENSURE
% THIS.
%

if ~COMPLEX,
  COMPLEX=1;
  PROCESSED=PROCESSED+1;
  N=PROCESSED;
  for S = 2:2:M,
    PS=PREVIOUS(S,1);
    M1=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S));
    PS=PREVIOUS(S,2);
    M2=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S)); 
    if M1 > M2,
      METRIC(S,N)=M1;
      SURVIVOR(S,N)=PREVIOUS(S,1);
    else
      METRIC(S,N)=M2;
      SURVIVOR(S,N)=PREVIOUS(S,2);
    end
  end
end

% NOW THAT WE HAVE MADE THE RUN-IN THE REST OF THE METRICS ARE
% CALCULATED IN THE STRAIGHT FORWARD MANNER. OBSERVE THAT ONLY
% THE RELEVANT STATES ARE CALCULATED, THAT IS REAL FOLLOWS COMPLEX
% AND VICE VERSA.
%
N=PROCESSED+1;
while N <= STEPS,
  for S = 1:2:M-1,
    PS=PREVIOUS(S,1);
    M1=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S));
    PS=PREVIOUS(S,2);
    M2=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S)); 
    if M1 > M2,
      METRIC(S,N)=M1;
      SURVIVOR(S,N)=PREVIOUS(S,1);     
    else
      METRIC(S,N)=M2;
      SURVIVOR(S,N)=PREVIOUS(S,2);
    end
  end
  N=N+1;
  for S = 2:2:M,
    PS=PREVIOUS(S,1);
    M1=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S));
    PS=PREVIOUS(S,2);
    M2=METRIC(PS,N-1)+real(conj(SYMBOLS(S,1))*Y(N)-INCREMENT(PS,S)); 
    if M1 > M2,
      METRIC(S,N)=M1;
      SURVIVOR(S,N)=PREVIOUS(S,1);     
    else
      METRIC(S,N)=M2;
      SURVIVOR(S,N)=PREVIOUS(S,2);
    end
  end
  N=N+1;
end  


% HAVING CALCULATED THE METRICS, THE MOST PROBABLE STATESEQUENCE IS
% INITIALIZED BY CHOOSING THE HIGHEST METRIC AMONG THE LEGAL STOP 
% STATES.
%
BEST_LEGAL=0;
for FINAL = STOPS,
  if METRIC(FINAL,STEPS) > BEST_LEGAL,
    S=FINAL;
    BEST_LEGAL=METRIC(FINAL,STEPS);
  end
end

% UNCOMMENT FOR TEST OF METRIC
%
% METRIC
% BEST_LEGAL
% S
% pause

% HAVING FOUND THE FINAL STATE, THE MSK SYMBOL SEQUENCE IS ESTABLISHED
%
IEST(STEPS)=SYMBOLS(S,1);
N=STEPS-1;
while N > 0,
  S=SURVIVOR(S,N+1);
  IEST(N)=SYMBOLS(S,1);
  N=N-1;
end

% THE ESTIMATE IS NOW FOUND FROM THE FORMULA:
% IEST(n)=j*rx_burst(n)*rx_burst(n-1)*IEST(n-1)
% THE FORMULA IS REWRITTEN AS:
% rx_burst(n)=IEST(n)/(j*rx_burst(n-1)*IEST(n-1))
% FOR INITIALIZATION THE FOLLOWING IS USED:
% IEST(0)=1 og rx_burst(0)=1
%
rx_burst(1)=IEST(1)/(j*1*1);
for n = 2:STEPS,
  rx_burst(n)=IEST(n)/(j*rx_burst(n-1)*IEST(n-1));
end

% rx_burst IS POLAR (-1 AND 1), THIS TRANSFORMS IT TO
% BINARY FORM (0 AND 1).
%
rx_burst_hard = (rx_burst+1)./2;

%generator the log likelihood ratio(LLR) of a symbol Sk. %  20070926
Y1 = [Y1 zeros(1,Lh)];   %in order to compensate the indices. 
IEST = [[-i 1 i -1] IEST [-i 1 i -1]]; %IEST = [mskmap(zeros(1,Lh)) IEST zeros(1,Lh)];

IEST = IEST((4-Lh)+1:length(IEST)-(4-Lh));

for stepk = 1:STEPS
    YHk = 0;              %define a temporary variable!
    for stepi = 1:Lh      %0:L-1 <=>1:L
        HSk = 0;          %define a temporary variable!
        for stepj = 1:Lh   
            if stepi ~= stepj
                HSk = HSk + Rhh(stepj)*IEST(Lh +stepk+stepi-stepj);%so the indices must be compensated by Lh.                
            end
        end
        YHk = YHk - abs( Y1(stepk+stepi-1)-HSk-Rhh(stepi)*IEST(stepk+Lh) )^2+...
            abs( Y1(stepk+stepi-1)-HSk-Rhh(stepi)*(-IEST(stepk+Lh) ))^2;  %IEST(stepk) = Al
    end
    LLR(stepk) = YHk;%LLR(stepk) = Lh*log(2*pi)/4*YHk; using sigma. 20070927
end
% LLR = abs(LLR);            %  keeping LLR to be positive. 20070927
rx_burst_soft = (-rx_burst).*LLR;%  20070926  -1->1,1->0



% %method 2,DERTA!
% % KNOWLEDGE OF Lh AND M IS NEEDED FOR THE ALGORITHM TO OPERATE
% 
% [ M , Lh ] = size(SYMBOLS);
% 
% % THE NUMBER OF STEPS IN THE VITERBI
% 
% STEPS = length(Y);
% 
% % INITIALIZE TABLES (THIS YIELDS A SLIGHT SPEEDUP).
% 
% METRIC = zeros(M,STEPS);
% SURVIVOR = zeros(M,STEPS);
% DERTA = zeros(M,STEPS)+10^2; %10^2=NaN   20070921
% SOFT_VALUE = zeros(1,STEPS); %  20070921
% 
% % DETERMINE PRECALCULATABLE PART OF METRIC
% 
% INCREMENT = make_increment(SYMBOLS,NEXT,Rhh);
% 
% % THE FIRST THING TO DO IS TO ROLL INTO THE ALGORITHM BY SPREADING OUT 
% % FROM 	THE START STATE TO ALL THE LEGAL STATES. 
% 
% PS = START;
% 
% % NOTE THAT THE START STATE IS REFERRED TO AS STATE TO TIME 0
% % AND THAT IT HAS NO METRIC.
% 
% S=NEXT(START,1);
% METRIC(S,1)=real(SYMBOLS(S,1)*Y(1))-INCREMENT(PS,S);
% SURVIVOR(S,1)=START;
% 
% S=NEXT(START,2);
% METRIC(S,1)=real(SYMBOLS(S,1)*Y(1))-INCREMENT(PS,S);
% SURVIVOR(S,1)=START;
% 
% PREVIOUS_STATES=NEXT(START,:);
% 
% % MARK THE NEXT STATES AS REAL. N.B: COMPLEX INDICATES THE POLARITY
% % OF THE NEXT STATE, E.G. STATE 2 IS REAL.
% 
% for N = 2:Lh  
%   STATE_CNTR=0;
%   
%   for PS = PREVIOUS_STATES,
%     STATE_CNTR=STATE_CNTR+1;
%     S=NEXT(PS,1);
%     METRIC(S,N)=METRIC(PS,N-1)+real(SYMBOLS(S,1)*Y(N))-INCREMENT(PS,S);
%     SURVIVOR(S,N)=PS;
%     USED(STATE_CNTR)=S;
%     
%     STATE_CNTR=STATE_CNTR+1;
%     S=NEXT(PS,2);
%     METRIC(S,N)=METRIC(PS,N-1)+real(SYMBOLS(S,1)*Y(N))-INCREMENT(PS,S);
%     SURVIVOR(S,N)=PS;    
%     USED(STATE_CNTR)=S;
%   end
%   PREVIOUS_STATES=USED;
% end
% 
% % NOW THAT WE HAVE MADE THE RUN-IN THE REST OF THE METRICS ARE
% % CALCULATED IN THE STRAIGHT FORWARD MANNER. OBSERVE THAT ONLY
% % THE RELEVANT STATES ARE CALCULATED, THAT IS REAL FOLLOWS COMPLEX
% % AND VICE VERSA.
% 
% N = Lh+1;
% while N <= STEPS,
%   for S = 1:M,
%     PS = PREVIOUS(S,1);
%     M1 = METRIC(PS,N-1) + real(SYMBOLS(S,1)*Y(N) - INCREMENT(PS,S));
%         
%     PS = PREVIOUS(S,2);
%     M2 = METRIC(PS,N-1) + real(SYMBOLS(S,1)*Y(N) - INCREMENT(PS,S)); 
%        
%     if M1 > M2
%       METRIC(S,N) = M1;
%       SURVIVOR(S,N) = PREVIOUS(S,1);
%       DERTA(S,N) = M1 - M2;%update the value 20070921
%     elseif M1 < M2
%         METRIC(S,N) = M2;
%         SURVIVOR(S,N) = PREVIOUS(S,2);
%         DERTA(S,N) = M2 - M1;
%     else
%         METRIC(S,N) = M2;
%         SURVIVOR(S,N) = PREVIOUS(S,2); 
%     end   
%   end
%   N = N+1;
% end  
% 
% %[S1,S] = max(METRIC(:, STEPS));
% if Lh == 4 & (METRIC(1, STEPS))<(METRIC(2, STEPS))%using STOPS
%     S = 2;
% else
%     S = 1;
% end
% % HAVING FOUND THE FINAL STATE, THE MSK SYMBOL SEQUENCE IS ESTABLISHED
% IEST(STEPS) = SYMBOLS(S,1);
% SOFT_VALUE(STEPS) = DERTA(S,STEPS);%   20070921
% N = STEPS - 1;
% 
% while N > 0,
%     S = SURVIVOR(S,N+1);
%     IEST(N) = SYMBOLS(S,1);
%     SOFT_VALUE(N) = DERTA(S,N);%   20070921
%     N = N-1;
% end
% rx_burst_hard = (1-IEST)/2;   %rx_burst_hard = (1-IEST)/2;20070921
% rx_burst_soft = IEST.*SOFT_VALUE;%   20070921    -1->1,1->0


% %method 3,OSA!
% % KNOWLEDGE OF Lh AND M IS NEEDED FOR THE ALGORITHM TO OPERATE
% 
% [ M , Lh ] = size(SYMBOLS);
% 
% % THE NUMBER OF STEPS IN THE VITERBI
% 
% STEPS = length(Y);
% 
% % INITIALIZE TABLES (THIS YIELDS A SLIGHT SPEEDUP).
% 
% METRIC = zeros(M,STEPS);
% SURVIVOR = zeros(M,STEPS);
% LLR = zeros(1,STEPS);  %  20070926
% 
% % DETERMINE PRECALCULATABLE PART OF METRIC
% 
% INCREMENT = make_increment(SYMBOLS,NEXT,Rhh);
% 
% % THE FIRST THING TO DO IS TO ROLL INTO THE ALGORITHM BY SPREADING OUT 
% % FROM 	THE START STATE TO ALL THE LEGAL STATES. 
% 
% PS = START;
% 
% % NOTE THAT THE START STATE IS REFERRED TO AS STATE TO TIME 0
% % AND THAT IT HAS NO METRIC.
% 
% S=NEXT(START,1);
% METRIC(S,1)=real(SYMBOLS(S,1)*Y(1))-INCREMENT(PS,S);
% SURVIVOR(S,1)=START;
% 
% S=NEXT(START,2);
% METRIC(S,1)=real(SYMBOLS(S,1)*Y(1))-INCREMENT(PS,S);
% SURVIVOR(S,1)=START;
% 
% PREVIOUS_STATES=NEXT(START,:);
% 
% % MARK THE NEXT STATES AS REAL. N.B: COMPLEX INDICATES THE POLARITY
% % OF THE NEXT STATE, E.G. STATE 2 IS REAL.
% 
% for N = 2:Lh  
%   STATE_CNTR=0;
%   
%   for PS = PREVIOUS_STATES,
%     STATE_CNTR=STATE_CNTR+1;
%     S=NEXT(PS,1);
%     METRIC(S,N)=METRIC(PS,N-1)+real(SYMBOLS(S,1)*Y(N))-INCREMENT(PS,S);
%     SURVIVOR(S,N)=PS;
%     USED(STATE_CNTR)=S;
%     
%     STATE_CNTR=STATE_CNTR+1;
%     S=NEXT(PS,2);
%     METRIC(S,N)=METRIC(PS,N-1)+real(SYMBOLS(S,1)*Y(N))-INCREMENT(PS,S);
%     SURVIVOR(S,N)=PS;    
%     USED(STATE_CNTR)=S;
%   end
%   PREVIOUS_STATES=USED;
% end
% 
% % NOW THAT WE HAVE MADE THE RUN-IN THE REST OF THE METRICS ARE
% % CALCULATED IN THE STRAIGHT FORWARD MANNER. OBSERVE THAT ONLY
% % THE RELEVANT STATES ARE CALCULATED, THAT IS REAL FOLLOWS COMPLEX
% % AND VICE VERSA.
% 
% N = Lh+1;
% while N <= STEPS,
%   for S = 1:M,
%     PS = PREVIOUS(S,1);
%     M1 = METRIC(PS,N-1) + real(SYMBOLS(S,1)*Y(N) - INCREMENT(PS,S));
%         
%     PS = PREVIOUS(S,2);
%     M2 = METRIC(PS,N-1) + real(SYMBOLS(S,1)*Y(N) - INCREMENT(PS,S)); 
%        
%     if M1 > M2
%       METRIC(S,N) = M1;
%       SURVIVOR(S,N) = PREVIOUS(S,1);
%     else
%         METRIC(S,N) = M2;
%         SURVIVOR(S,N) = PREVIOUS(S,2); 
%     end   
%   end
%   N = N+1;
% end  
% 
% %[S1,S] = max(METRIC(:, STEPS));
% if Lh == 4 & (METRIC(1, STEPS))<(METRIC(2, STEPS))%using STOPS
%     S = 2;
% else
%     S = 1;
% end
% % HAVING FOUND THE FINAL STATE, THE MSK SYMBOL SEQUENCE IS ESTABLISHED
% IEST(STEPS) = SYMBOLS(S,1);
% N = STEPS - 1;
% 
% while N > 0,
%     S = SURVIVOR(S,N+1);
%     IEST(N) = SYMBOLS(S,1);    
%     N = N-1;
% end
% rx_burst_hard = (1-IEST)/2;   %rx_burst_hard = (1-IEST)/2;20070921
%    
% %generator the log likelihood ratio(LLR) of a symbol Sk. %  20070926
% Y1 = [Y1 zeros(1,Lh)];   %in order to compensate the indices. 
% IEST1 = T_SEQ_MSKMAP([zeros(1,Lh) rx_burst_hard zeros(1,Lh)]);
% IEST1 = -IEST1;          %IEST = [zeros(1,Lh) IEST zeros(1,Lh)];
% for stepk = 1:STEPS
%     YHk = 0;              %define a temporary variable!
%     for stepi = 1:Lh      %0:L-1 <=>1:L
%         HSk = 0;          %define a temporary variable!
%         for stepj = 1:Lh   
%             if stepi ~= stepj
%                 HSk = HSk + Rhh(stepj)*IEST1(Lh+stepk+stepi-stepj);%so the indices must be compensated by Lh.                
%             end
%         end
%         YHk = YHk + abs( Y1(stepk+stepi)-HSk-Rhh(stepi)*IEST1(stepk) )^2;%IEST(stepk) = Al 
%     end
%     LLR(stepk) = Lh*log(2*pi*sigma^2)/(4*sigma^2)*YHk;%Lh*log(2*pi*sigma^2)/(4*sigma^2)*YHk; using sigma. 20070927
% end
% LLR = abs(LLR);            %  keeping LLR to be positive. 20070927
% 
% rx_burst_soft = IEST(1:length(IEST)).*LLR;%  20070926  -1->1,1->0