%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXAMPLE: mlse_viterbi(100, 5, 5, 3, 1, 0, 0, [1], 0, 0, "mlse")
% Stephan ten Brink, Institute of Telecommunications, University of Stuttgart
% PARAMETERS:

% notes
% - assume start in -1 -1 -1 -1 -1 transmitted before (easier)

% notes:
% - output BER (over 10000) at Es/N0=5dB (real) to compare with other equalizers

% EXAMPLES:
% mlse_viterbi(100, 50, -3, 4, 1, 1, 0, 1, [1], 0, 0, "mlse")     % plot interim survivors

% quick test
% mlse_viterbi(100, 6, 3, 2, 1, 1, 0, 0, [1], 0, 0, "mlse")

% autocorr...
% mlse_viterbi(100000, 5, 3, 4, 1, 1, 1, 0, [1], 0, 0, "mlse")


% webdemo
% - show up to time K
% - start_all_minus1_en:
%   + always starts in state 0
%   + this flag indicates, whether for k<0 we assume "0" in the channel shift register, or -1 (if set to 1, DEFAULT)
% examples


% BEGIN, main function (control)
function mlse_viterbi(
  Nsym,                         % nb of BPSK symbols
  K,                            % plot trellis up to K symbols
  EsN0dB,                       % Es/N0 in dB, note: BPSK over real channel here
  h_index,                      % pick one of 4 pre-selected channel impulse responses
  start_all_minus1_en,          % if set, then assume -1,-1..-1 for k<0 in channel shift register (otherwise 0-entries)
  rnd_data_en,                  % new random sequence every time?

  display_signal,               % 0==trellis
                                % 1==t-domain before channel
                                % 2==t-domain after channel, with noise
                                % 3==auto-corr before channel
                                % 4==auto-corr after channel with noise


  plot_interim_survivors_en,    % plot interim survivors; default is 0 (i.e., no)
  plot_en,                      % enable plots 0, 1 [t-domain, f-domain, ampl/phase diagram, h_rc_n]
  svg_en,                       % if 1, then also output svg-file
  webdemo_en,                   % if 0, then test-mode in Octave/Matlab commandline, no webserver running
  filename                      % output file for webdemo
)

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN check input parameters for consistency

  % only one plot at a time
  if (sum(plot_en)!=1)
    plot_en = [1];
    warning("setting  'plot_en = [1]'");
  end

  if (Nsym>100000)  Nsym=100000; warning("setting  'Nsym = 100000'");  end
  if (Nsym<10)      Nsym=10;     warning("setting  'Nsym = 10'");      end

  if (K>100)  K=100;      warning("setting  'K = 100'");  end
  if (K<1)    K=1;        warning("setting  'K = 1'");    end

  if (K>Nsym) Nsym=K;     warning("setting  'Nsym=K'");   end

  if (h_index>5)    h_index=5;   warning("setting  'h_index=5'");  end
  if (h_index<0)    h_index=0;   warning("setting  'h_index=0'");  end

  % END check input parameters for consistency
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm; from dark to light
  gray =    [[0,0,0]; [32,32,32]; [64,64,64]; [96,96,96]; [128,128,128]; [160,160,160]; [192,192,192]; [224,224,224]; [255,255,255]]/255;
  violett = [[51,0,25]; [102,0,51]; [153,0,76]; [204,0,102]; [255,0,127]; [255,51,153]; [255,102,178]; [255,153,204]; [255,204,229]]/255;
  magenta = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
  lilac   = [[25,0,51]; [51,0,102]; [76,0,153]; [102,0,204]; [127,0,255]; [153,51,255]; [178,102,255]; [204,153,255]; [229,204,255]]/255;
  blue    = [[0,0,51]; [0,0,102]; [0,0,153]; [0,0,204]; [0,0,255]; [51,51,255]; [102,102,255]; [153,153,255]; [204,204,255]]/255;
  seablue = [[0,25,51]; [0,51,102]; [0,76,153]; [0,102,204]; [0,128,255]; [51,153,255]; [102,178,255]; [153,204,255]; [204,229,255]]/255;
  cyan    = [[0,51,51]; [0,102,102]; [0,153,153]; [0,204,204]; [0,255,255]; [51,255,255]; [102,255,255]; [153,255,255]; [204,255,255]]/255;
  seagreen = [[0,51,25]; [0,102,51]; [0,153,76]; [0,204,102]; [0,255,128]; [51,255,153]; [102,255,178]; [153,255,204]; [204,255,229]]/255;
  green   = [[0,51,0]; [0,102,0]; [0,153,0]; [0,204,0]; [0,255,0]; [51,255,51]; [102,255,102]; [153,255,153]; [204,255,204]]/255;
  lightgreen = [[25,51,0]; [51,102,0]; [76,153,0]; [102,204,0]; [128,255,0]; [153,255,51]; [178,255,102]; [204,255,153]; [229,255,204]]/255;
  yellow  = [[51,51,0]; [102,102,0]; [153,153,0]; [204,204,0]; [255,255,0]; [255,255,51]; [255,255,102]; [255,255,153]; [255,255,204]]/255;
  amber   = [[51,25,0]; [102,51,0]; [153,76,0]; [204,102,0]; [255,128,0]; [255,153,51]; [255,178,102]; [255,204,153]; [255,229,204]]/255;
  red     = [[51,0,0]; [102,0,0]; [153,0,0]; [204,0,0]; [255,0,0]; [255,51,51]; [255,102,102]; [255,153,153]; [255,204,204]]/255;




  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, set up graphical output

  % Create an invisible figure.
  if (webdemo_en)
    fig = figure(1); set(fig, "visible", "off");
  else
    fig = figure(1); set(fig, "visible", "on");
  end

  landscape = "-S930,350"; landscape_ratio = 930/350;
  portrait  = "-S640,480"; portrait_ratio = 640/480;

  if (display_signal==0)
    output_format = landscape;
  else
    output_format = portrait;
  end

%  if (svg_en) || (webdemo_en)
    output_aspect_ratio = landscape_ratio;
%  else
%    output_aspect_ratio = 4/3;
%  end


  set(0, "defaultlinelinewidth", 2);


  markersize_states = 5;

  if (webdemo_en)
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 4;
    fontsize_title = 6;
    fontsize_rx_seq = 7;
    fontsize_bm = 6;
    fontsize_pm_big =8;
    fontsize_pm_small =5;
    fontsize_pm_crossout =10;
  elseif (svg_en)
    set(0, "defaultaxesfontsize", 6);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 6);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 1;
    const_msize = 3;
    fontsize_title = 6;
    fontsize_rx_seq = 8;
    fontsize_bm = 8;
    fontsize_pm_big =10;
    fontsize_pm_small =6;
    fontsize_pm_crossout =12;
  else
    set(0, "defaultaxesfontsize", 10);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");
    const_lw = 2;
    const_msize = 2;
    fontsize_title = 6;
    fontsize_rx_seq = 8;
    fontsize_bm = 8;
    fontsize_pm_big =12;
    fontsize_pm_small =8;
    fontsize_pm_crossout =12;
  end

  % END, set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




  % reset random variable generator
  if (rnd_data_en==0)
    rnd_val_gen_seed = (5:205); % arbitrary vector of length <=625 to init. rnd val generator to always same sequence
    rand("state", rnd_val_gen_seed);
  end


  sigma_n= sqrt(1/ (10^(EsN0dB/10)));


  h0 = [1, 0.5];
  h1 = [1, -3/4, -5/8];      % from script, 1/H(z) instable (ZF, MMSE)
  h2 = [1, 1/2, -1/4];       % from script, 1/H(z) stable   (ZF, MMSE)
  h3 = [1, 1, -0.5];         % from script                  (MLSE, Viterbi)
  h4 = [1, -0.5, 0.75,0.2];
  h5 = [1, 0.5, -0.2, 0.4, -0.8, 0.4, 0.3, -0.2, 0.2];

  if     (h_index==0)   h=h0;
  elseif (h_index==1)   h=h1;
  elseif (h_index==2)   h=h2;
  elseif (h_index==3)   h=h3;
  elseif (h_index==4)   h=h4;
  elseif (h_index==5)   h=h5;
  else                  h=h0;
  end

  L = length(h)-1;      % order of channel FIR filter; memory; L+1 taps (coefficients)

  LUT_metric_increment = compute_LUT_metric_increment(h);


  h_gain = sqrt(sum(h.^2));
  sigma_n = sigma_n*(h_gain);

  % transmitted signal
  s = get_random_BPSK_vector(Nsym);

  % noise on channel
  if (rnd_data_en==0)
    randn("state", rnd_val_gen_seed);
  end

  noise = randn(1,Nsym+L)*sigma_n;


  % convolution with channel
  if (start_all_minus1_en)
    q = conv(h,[-ones(1,L), s])(L+1:end) + noise;
  else
    q = conv(h,s) + noise;
  endif

  if (display_signal==0)
    nb_states = 2^L;

    [bm, pm, pm01, tb_bit] = compute_branch_path_metrics(q, nb_states, Nsym, LUT_metric_increment, h, start_all_minus1_en);

    y=determine_symbol_sequence(pm, tb_bit, nb_states, Nsym);

    BER = sum(abs(s(1:Nsym) .-  y(1:Nsym)))/2.0 / Nsym;
  end




  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, generate plots
  % note: always generate plots from main function

  clf
  hold on

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if (plot_en(1) && (display_signal==0))

    plot_rx_seq_en = 1;
    plot_bm_en = 1;
    plot_pm_en = 1;
    plot_transitions_en = 1;
    plot_states_en = 1;


    if (K<=5)
      xmin=-0.3; xmax=5+0.3;
    else
      xmin=-0.5; xmax=K+0.5;
      plot_rx_seq_en = 0;
      plot_bm_en = 0;
      plot_pm_en = 0;
      markersize_states = 3;
    end

    if (K>50)
      plot_transitions_en = 0;
    end


%    xmin=-0.3; xmax=5+0.3;
    ymin=-0.4; ymax=nb_states-0.5 + 0.3;

    set(gca,'YTick',(0:1:nb_states-1)); % 'fontsize',12
    xtickstep=ceil(K/20);
    set(gca,'XTick',(0:xtickstep:K));

    xlim([xmin, xmax]);
    ylim([ymin, ymax]);
    xlabel("discrete-time k");
    ylabel("state number S");

    box off;
%    set(gca,'Box','off');   % does not work in linux...
    set(gca,'LineWidth',-1); % does not work in linux...
%    set(gca,'color',gray(9,:)); % does this work?

%    set(gca,'XColor',lightgreen(4,:))
%    set(gca,'YColor',amber(4,:))



    %%%%%%%%%%%%%%%%%%%%%%%%%%% plot title

    Hstr=["H(z)=",num2str(h(1))];
    for k=1:L
      if (h(k+1)<0)
        Hstr=[Hstr,num2str(h(k+1)),"z^{",num2str(-k),"}"];
      else
        Hstr=[Hstr,"+",num2str(h(k+1)),"z^{",num2str(-k),"}"];
      end
    end

    title([Hstr,",  E_s/N_0=",num2str(EsN0dB),"dB,  BER=", num2str(BER)],'fontsize', fontsize_title);

    % compute ratio for correct plotting of rotated text
    rot_ratio =  (xmax-xmin)/(ymax-ymin) * 1/output_aspect_ratio;


    %%%%%%%%%%%%%%%%%%%%%%%%%%% plot noisy received sequence
    if (plot_rx_seq_en)
      ypos=ymax-0.3;
      text(0, ypos, '(x_k|z_k)', 'fontsize', fontsize_rx_seq, 'color', green(3,:), 'horizontalalignment', 'center');

      for k=1:K
        xpos=k-0.5;
        text(xpos, ypos, ["(",num2str(s(k)),"|",num2str(q(k),'%5.2f'),")"], 'fontsize', fontsize_rx_seq, 'color', green(3,:), 'horizontalalignment', 'center');
      end
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%% plot surviving path

    if (plot_interim_survivors_en)
      kk_start=2;
    else
      kk_start=K;
    end

% evtl. include interim survivors
for kk=kk_start:K

    % 1. determine best path metric (smallest)
    pm_min=1000;
    s_min=0;

    for s=0:nb_states-1
      if (pm(s+1,kk+1)<pm_min)
        pm_min=pm(s+1,kk+1);
        s_min=s;
      end
    end

    % 2. trace back surviving path

    s=s_min;     % end-state
    for k=kk:-1:1
      % prev state
      if (tb_bit(s+1, k+1)==0)
        s_prev = floor(s/2);
      else
        s_prev = floor(s/2)+nb_states/2;
      end

      if (kk==K)
        trans_color = yellow(5,:);
      else
        trans_color = yellow(3,:);
      end

      plot([k-1, k], [s_prev,s], "-", 'linewidth',8, "color", trans_color);
      s=s_prev;
    end

end

    %%%%%%%%%%%%%%%%%%%%%%%%%%% plot transitions
    if (plot_transitions_en)
      for k=1:K
        for s=0:nb_states-1
          if (mod(s,2)==0)
            if (plot_interim_survivors_en)  trans_color = blue(9,:);  else  trans_color = blue(4,:); end
          else
            if (plot_interim_survivors_en)  trans_color = red(9,:);  else  trans_color = red(4,:); end
          end
          plot([k-1, k], [floor(s/2),s], "-", 'linewidth',1, "color", trans_color);
          plot([k-1, k], [floor(s/2)+floor(nb_states/2),s], "-", 'linewidth',1, "color", trans_color);
        end
      end
    end


    %%%%%%%%%%%%%%%%%%%%%%%%%%% plot states
    if (plot_states_en)
      for k=0:K
        for s=0:nb_states-1
          if (plot_interim_survivors_en)
            plot(k, s, "o", "markersize", markersize_states, 'linewidth',2, "color", gray(7,:), 'MarkerFaceColor',gray(8,:));
          else
            plot(k, s, "o", "markersize", markersize_states, 'linewidth',2, "color", "k", 'MarkerFaceColor',gray(7,:));
          end
        end
      end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%% annotate with bm
    if (plot_bm_en)

      for k=1:K
        for s=0:nb_states-1

          if (mod(s,2)==0)
            trans_color = 'b';
          else
            trans_color = 'r';
          end

          slope0=(s-floor(s/2))/1  * rot_ratio;
          slope0_deg = atan(slope0)*180/pi;

          xpos = k-0.5;
          ypos = (s-floor(s/2))/2+floor(s/2)-0.2   + 0.05;
          text(xpos, ypos, num2str(bm(s+1,k),'%5.2f'), 'fontsize', fontsize_bm, 'color', trans_color,...
               'horizontalalignment', 'center',...
               'rotation',slope0_deg);

          slope1=(s-(floor(s/2)+nb_states/2))/1  * rot_ratio;
          slope1_deg = atan(slope1)*180/pi;

          xpos = k-0.25;
          ypos = -(s-(floor(s/2)+nb_states/2))*0.25+s+0.1  + 0.02;
          text(xpos, ypos, num2str(bm(s+1+nb_states,k),'%5.2f'), 'fontsize', fontsize_bm, 'color', trans_color,...
               'horizontalalignment', 'center',...
               'rotation',slope1_deg);
        end
      end

    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%% annotate with accumulated pm

    if (plot_pm_en)

      if (L==2)
        cross_out_pos = 1/4;
      elseif (L==3)
        cross_out_pos = 1/3;
      else
        cross_out_pos = 1/3;
      end

      for k=1:K+1
        for s=0:nb_states-1
          if (pm(s+1,k)==pm01(s+1,k))
            pm0_color='m';
            pm0_fontsize=fontsize_pm_big;
            pm1_color=gray(6,:);
            pm1_fontsize=fontsize_pm_small;

            % respective other path crossed out
            xpos_crossout = k-1-cross_out_pos;
            ypos_crossout = -(s-(floor(s/2)+nb_states/2))*cross_out_pos+s;

          else
            pm0_color=gray(6,:);
            pm0_fontsize=fontsize_pm_small;
            pm1_color='m';
            pm1_fontsize=fontsize_pm_big;

            % respective other path crossed out
            xpos_crossout = k-1-cross_out_pos;
            ypos_crossout = (s-floor(s/2))*(1-cross_out_pos)+floor(s/2);

          end

          xpos = k-1;
          ypos = s-0.2;

          text(xpos, ypos, num2str(pm01(s+1,k),'%5.2f'), 'fontsize', pm0_fontsize, 'color', pm0_color,...
               'horizontalalignment', 'center');

          ypos = s+0.2;
          text(xpos, ypos, num2str(pm01(s+1+nb_states,k),'%5.2f'), 'fontsize', pm1_fontsize, 'color', pm1_color,...
               'horizontalalignment', 'center');

          if (k>1)
            text(xpos_crossout, ypos_crossout+0.05, '||', 'fontsize', fontsize_pm_crossout, 'color',  gray(3,:),'horizontalalignment', 'center');
          end

        end
      end

    end

  end

  %
  %%%%%%%%%%%
  %

  %%%%%%%%%%%%%%%%%%%%%%%%%%% plot title

  Hstr=["H(z)=",num2str(h(1))];
  for k=1:L
    if (h(k+1)<0)
      Hstr=[Hstr,num2str(h(k+1)),"z^{",num2str(-k),"}"];
    else
      Hstr=[Hstr,"+",num2str(h(k+1)),"z^{",num2str(-k),"}"];
    end
  end


  % display different signals, TIME-DOMAIN
  if (plot_en(1) && ((display_signal==1) || (display_signal==2)))
    title([Hstr,",  E_s/N_0=",num2str(EsN0dB),"dB"],'fontsize', fontsize_title);
%    ylim([-2,2]);
    xlabel("discrete time index k");
    ylabel("discrete-time sample-value");

    if (display_signal==1)
      y=s;
      ylim([-2, 2]);
    elseif (display_signal==2)
      y=q;
    end

    stem((1:K)-1, y(1:K), 'linewidth',2, "markerfacecolor", blue(8,:), "markersize", 3, "color", blue(3,:));
  end

  % display different signals, AUTO-CORRELATION
  if (plot_en(1) && ((display_signal==3) || (display_signal==4)))
    title([Hstr,",  E_s/N_0=",num2str(EsN0dB),"dB"],'fontsize', fontsize_title);
    xlabel("discrete time lag k");
    ylabel("auto-correlation value");

    if     (display_signal==3) y=s;
    elseif (display_signal==4) y=q;
    end

    nb_samples_display=K;
    acorr = compute_auto_correlation_normalized(y, nb_samples_display);
    stem((-(length(acorr)-1):1:length(acorr)-1), [acorr(end:-1:1), acorr(2:end)], 'linewidth',2, "markerfacecolor", green(8,:), "markersize", 3, "color", green(3,:));
  end


  % END, generate plots
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % BEGIN, output to file

  print(filename, "-dpng",output_format);

  if (svg_en)
    print(["fig_",filename], "-dsvg",output_format);
  end

  % END, output to file
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


end

% END, main program (control)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%














%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN, function definitions

function y=determine_symbol_sequence(pm, tb_bit, nb_states, K)
  % 1. determine best path metric (smallest)
  pm_min=1000;
  s_min=0;

  y=zeros(1,K);

  for s=0:nb_states-1
    if (pm(s+1,K+1)<pm_min)
      pm_min=pm(s+1,K+1);
      s_min=s;
    end
  end

  % 2. trace back surviving path

  s=s_min;     % end-state
  for k=K:-1:1
    % prev state
    if (tb_bit(s+1, k+1)==0)
      s_prev = floor(s/2);
    else
      s_prev = floor(s/2)+nb_states/2;
    end

    symbol = 2*mod(s,2)-1;
    y(k)=symbol;

    s=s_prev;
  end

end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=get_random_BPSK_vector(N)
  M = 2;                                  % nb of ampl. levels
  ampl_vec = (-M+1:2:M-1);                % ampl. levels
  prob_vec = ones(1,M)*1/M;               % all equally probable
  P_av = 1/M*sum(ampl_vec.^2);
  y=discrete_rnd(ampl_vec/sqrt(P_av), prob_vec, 1, N);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LUT is organized as [bit transition 0 or 1 (input to FIR), state of FIR tapped delay line]
function y = compute_LUT_metric_increment(h)
  L = length(h)-1;
  nb_states = 2^L;

  y=zeros(2, nb_states);

  for xval = 0:nb_states-1
    xvec = 2*bitget(xval,1:L)-1; % vector of -1, 1 (BPSK); LSB of xval is on the leftmost side, ok.
    inp_vec = [-1, xvec];
    y(0 +1, xval+1) = sum(inp_vec .* h);    % inp_vec(1) stands for x_k, inp_vec(1) for x_{k-1} etc...
    inp_vec = [+1, xvec];
    y(1 +1, xval+1) = sum(inp_vec .* h);    % inp_vec(1) stands for x_k, inp_vec(1) for x_{k-1} etc...
  end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% account for 0-entries in tapped delay line at beginning  (not only -1, +1)
function y=get_metric_increment(bit_transition, prev_state, k, h, LUT_metric_increment, start_all_minus1_en)
  L = length(h)-1;

  if ((start_all_minus1_en) || (k-1>=L)) % then easier to understand
    y=LUT_metric_increment(bit_transition +1, prev_state +1);
  else % (k-1<L)  % special case for boundary (0 entries in tapped delay line, not just -1, +1)
    if (k-1==0)
      xvec=[];
    else
      xvec = 2*bitget(prev_state, 1:(k-1))-1;   % vector of -1, 1 (BPSK); LSB of xval is on the leftmost side, ok.
    end

    inp_vec = [2*bit_transition-1, xvec, zeros(1,L-(k-1))];
    y = sum(inp_vec .* h);    % inp_vec(1) stands for x_k, inp_vec(1) for x_{k-1} etc...
  endif

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  [bm, pm, pm01, tb_bit] = compute_branch_path_metrics(q, nb_states, K, LUT_metric_increment, h,start_all_minus1_en)

  pm01 = zeros(2*nb_states, K+1);          % path metrics
  bm = zeros(2*nb_states, K+1);            % branch metrics
  pm = zeros(nb_states, K+1);              % path metrics
  tb_bit = zeros(nb_states, K+1);          % path metrics

  % initialize for discrete time k=0
  for s = 1:nb_states
    if (s==1)
      pm(1,0+1) = 0; % start off from state 0
    else
      pm(s,0+1) = 1000;  % no start from other states
    end

    % just init; keep for plotting
    pm01(s,1)=pm(s,1);
    pm01(s+nb_states,1)=pm(s,1);
  end



  % note: in index +1... since not from 0
  % also: LSB of s is input to FIR channel, i.e. bit transition

  % loop over discrete-time k
  for k = 2:K+1
    for s = 0:nb_states-1

      bit_transition = mod(s,2); % 0, 1

      prev_state0 = floor(s/2);
      prev_state1 = prev_state0+floor(nb_states/2);

      % metric increment; Euclidean distance
      w_hyp0 = (q(k-1) - get_metric_increment(bit_transition, prev_state0, k-1, h, LUT_metric_increment,start_all_minus1_en))^2;
      w_hyp1 = (q(k-1) - get_metric_increment(bit_transition, prev_state1, k-1, h, LUT_metric_increment,start_all_minus1_en))^2;

      % branch metric
      bm(s+1,k-1)           = w_hyp0;
      bm(nb_states+s+1,k-1) = w_hyp1;

      % accumulated path metric
      pm0 = pm(prev_state0+1,k-1)+bm(s+1,k-1);
      pm1 = pm(prev_state1+1,k-1)+bm(nb_states+s+1,k-1);

      % store for plotting
      pm01(s+1,k)           = pm0;
      pm01(nb_states+s+1,k) = pm1;

      % determine surviving path metric at this node and set traceback-bit
      % smaller acc. path metric wins (minimize Euclidean distance)
      if (pm0<pm1)
        pm(s+1,k) = pm0;
        tb_bit(s+1,k) = 0;
      else
        pm(s+1,k) = pm1;
        tb_bit(s+1,k) = 1;
      end

    end % for s
  end % for k

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y_vec = compute_auto_correlation_normalized(x, max_lag)
  len = length(x);
  norm=sum(x.*conj(x));    % so that we get a normalized auto-corr.

  y_vec = zeros(1,max_lag);
  for t_lag = (0:max_lag-1)
    x1 = [ x, zeros(1,t_lag) ];
    x2 = [ zeros(1,t_lag), x ];

    y_vec(1+t_lag)=sum(x1 .* conj(x2))/norm;
  end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% END, function definitions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
