% S:Qingle-carrier, Institute of Telecommunications, University of Stuttgart, Germany

% www:Q.inue.uni-stuttgart.de

% author: Maximilian Arnold, maximilian.wolfgang.arnold@gmail.com

% date: 27.6.2015

% update: 27.6.2015

%

% Written for MATLAB

%

% EXAMPLE:  SC_noaverage(0,2048, 100, 4, 2, 1,0,0,0,'OFDM_runN')

% PARAMETERS:

% channel_model                   % Modus of Channel

%  num_bit,                       % number of bits per symbol

%  number_of_transmissions,       % number of transmissions over channel

%  order_of_equalizer             % Order of Equalizer

%  L,    		                  % number of Taps of Channel

%  Delay_of_channel,			  % Delay of channel

%  keepPlot,                      % if 1, then the results of the previous plot are shown in the following plot as well.

%  svg_en,                        % if 1, then also output svg-file

%  webdemo_en,                    % if 0, then test-mode in Octave/Matlab commandline, no need for webserver running

%  filename                       % output file for webdemo





% BEGIN, main function (control)



function SC_noaverage(channel_model,num_bit, number_of_transmissions,order_of_equalizer, L,keepPlot,svg_en,webdemo_en,filename)

 

% default parameters

switch nargin

    case 0

        channel_model = 1;

        num_bit = 2048;

        number_of_transmissions =100;

        order_of_equalizer =4;

        L=2;

        hideLegend = 1;

        keepPlot = 0;

        svg_en = 0;

        webdemo_en = 0;

        filename = 'SC_noaverage';

    otherwise

end



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

% 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';

portrait  = '-S640,480';

output_format = portrait; %portrait;landscape;



set(0, 'defaultlinelinewidth', 1);



if (webdemo_en)

    set(0, 'defaultaxesfontsize', 8);

    set(0, 'defaultaxesfontname', 'Helvetica');

    set(0, 'defaulttextfontsize', 8);

    set(0, 'defaulttextfontname', 'Helvetica');

    const_lw = 2;

elseif (svg_en)

    set(0, 'defaultaxesfontsize', 6);

    set(0, 'defaultaxesfontname', 'Helvetica');

    set(0, 'defaulttextfontsize', 6);

    set(0, 'defaulttextfontname', 'Helvetica');

    const_lw = 1;

else

    set(0, 'defaultaxesfontsize', 10);

    set(0, 'defaultaxesfontname', 'Helvetica');

    set(0, 'defaulttextfontsize', 8);

    set(0, 'defaulttextfontname', 'Helvetica');

    const_lw = 2;

end



% RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm

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;





% END, set up graphical output

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





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

% BEGIN check input parameters for consistency

if(num_bit>4097)

    warning('Number of bits is too high')

end

if(num_bit<=0)

    warning('Number of bits needs to positive and bigger than zero')

end

if(number_of_transmissions>101 && number_of_transmissions<5)

    warning('Number of transmission needs to be in range 5-100')

end

if(order_of_equalizer<0 && order_of_equalizer>101)

     warning('Order of equalizer needs to be in range 0-100')

end
if(L<1 && L>50)

    warning('Order of channel needs to be in range 1-50')

end

if(L<1 && L>50)

    warning('Order of channel needs to be in range 1-50')

end


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



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

%Definitions

Eb=1;                                            %Bit energy

SNRdB=0:1:15.5;                                  %Signal to Noise Ratio (in dB)

SNR=10.^(SNRdB/10);                              %Signal to Noise Ratio



% Prealocating for speed

BER_NoChannel(length(SNR))=0;                    %Initialize BER_NoChannel

BER_GivenChannel_ZFLS(length(SNR))=0;            %Initialize BER_GivenChannel_ZFLS

BER_GivenChannel_MMSE(length(SNR))=0;            %Initialize BER_GivenChannel_MMSE

BER_GivenChannel_WithoutEQ(length(SNR))=0;       %Initialize BER_GivenChannel_WithoutEQ



Number_of_States = 2^L;                          %Number_Of_States


% Viterbi Variables:

BER_GivenChannel_Viterbi(length(SNR))=0;

StateBits(L,Number_of_States)=0;                 %Bits Im Register für die States(States 0 =1 , 1 = 2 , ..)

for States=1:Number_of_States

    StateBits(:,States)= de2bi(States-1,L);

    Shiftbits(States,:) = de2bi(States-1,L);

    Path_TraceBack(States,:) = [bi2de([Shiftbits(States,L:-1:2) 0]), bi2de([Shiftbits(States,L:-1:2) 1])];

end

StateBits=2*StateBits-1;                         %Bits for Output Y

Initial_RegisterState = 0;                       %Initial RegisterState

Path_Bits(:,num_bit) = 0;

next_path_metric(2) = 0;

Path_Metric=1000*ones(Number_of_States,num_bit);  %Initialize PathMetric





%% Channels

% Given Channel

if (channel_model==0)

h_GivenChannel(1,1)=1;

for channel_stable_index=2:(L+1)

    h_GivenChannel(1,channel_stable_index) = (-1)^channel_stable_index*1/(2^(channel_stable_index-1));
    i = zeros(L+order_of_equalizer+1,1);
    Delay_of_channel =1;
    i(Delay_of_channel,1) = 1;

end

end

if (channel_model==1)

    h_GivenChannel(1,:) = normrnd(0,1/(L+1),[1,L+1]);
    [value Delay_of_channel] = max(abs(h_GivenChannel));
    i = zeros(L+order_of_equalizer+1,1);
    i(Delay_of_channel,1) = 1;
end


%% Normalize Channels

%GivenChannel

h_GivenChannel = h_GivenChannel./sqrt(h_GivenChannel*h_GivenChannel');





%% Calculation Coeffiencents of Filters

%GivenChannel ZFLS

H_GivenChannel = zeros(order_of_equalizer+L+1,order_of_equalizer+1);

stable_matrix_j=1;

for stable_matrix_i=1:(order_of_equalizer+1)

    k=1;

    for stable_matrix_j=stable_matrix_j:(stable_matrix_j+L)

        H_GivenChannel(stable_matrix_j,stable_matrix_i) = h_GivenChannel(1,k);

        k=k+1;

    end

    stable_matrix_j = stable_matrix_j+1-L;

end

f_ZFLS_GivenChannel = (inv(ctranspose(H_GivenChannel)*H_GivenChannel))*ctranspose(H_GivenChannel)*i; % Faster than H\i



%GivenChannel MMSE

for count=1:length(SNR)

    f_MMSE_GivenChannel(:,count) = (inv(ctranspose(H_GivenChannel)*H_GivenChannel+(1/(SNR(count)))*eye(order_of_equalizer+1, order_of_equalizer+1)))*ctranspose(H_GivenChannel)*i;

end





%% Normalize filter Coeficients:

%MMSE

for count=1:length(SNR)

    f_MMSE_GivenChannel(:,count) = f_MMSE_GivenChannel(:,count)./sqrt(f_MMSE_GivenChannel(:,count)'*f_MMSE_GivenChannel(:,count));

end



%% Calculating Y_past for all States



Y_Past_GivenChannel_Viterbi = StateBits'*h_GivenChannel(2:L+1)';

Y_Current_Given_Channel_Viterbi(:,1) =  Y_Past_GivenChannel_Viterbi-1*h_GivenChannel(1);% 1 ist Input -1/ 2 ist Input 1

Y_Current_Given_Channel_Viterbi(:,2) =  Y_Past_GivenChannel_Viterbi+1*h_GivenChannel(1);% 2 ist Input 1



%% Calculating BER



% Sending Packets over Channel

for transmission_counter=1:number_of_transmissions

    for count=1:length(SNR)

        %% Modulate the Signal and Preconditions

        

        No=Eb/SNR(count);                    %Calculate noise power from SNR to Symbol for SC

        data=randi([0,1],1,num_bit);         %Generate binary data source

        s=2*data-1;                          %BPSK modulation

        N=normrnd(0,sqrt(No/2),[1,num_bit]); %Generate AWGN, Normally distributed random number 1/2 for 2sigma^2



 

        %% Calculate Noise received Signals

        

        Y_NoChannel =s+N;

        Y_GivenChannel_WithoutEq = filter(h_GivenChannel, 1,s) + N;

        Y_GivenChannel = filter(h_GivenChannel, 1,s)+N;



        

        %% Linear Equalizing:

        

        Y_GivenChannel_ZFLS = filter(f_ZFLS_GivenChannel, 1,Y_GivenChannel);

        Y_GivenChannel_MMSE = filter(f_MMSE_GivenChannel(:,count),1,Y_GivenChannel);



        

        %% Viterbi Equalizing:

        % Pilot Bits /Aktueller Zustand

        Path_Metric(bi2de(data(2:-1:1))+1,L+1)=0;

        

        % Alle PathMetriken rechnen

        for bit=1+L:num_bit

            for state=1:Number_of_States

                next_path_metric(1) = Path_Metric(Path_TraceBack(state,1)+1,bit) + (Y_GivenChannel(1,bit)-Y_Current_Given_Channel_Viterbi(Path_TraceBack(state,1)+1,Shiftbits(state,1)+1))^2;

                next_path_metric(2) = Path_Metric(Path_TraceBack(state,2)+1,bit) + (Y_GivenChannel(1,bit)-Y_Current_Given_Channel_Viterbi(Path_TraceBack(state,2)+1,Shiftbits(state,1)+1))^2;

                if(next_path_metric(1) < next_path_metric(2))

                    Path_Metric(state,bit+1) =  next_path_metric(1);

                else

                    Path_Metric(state,bit+1) =  next_path_metric(2);

                end

            end

        end

        

        

        %Traceback

        [~,statecounter(1,num_bit+1)] = min(Path_Metric(:,num_bit+1));

        for k=num_bit:-1:1+L

            [~, possible_pre_state] = min(Path_Metric(:,k));

            for states_runner=1:Number_of_States

                if((mod(2*(possible_pre_state-1),2^L)+0)==statecounter(1,k+1)-1)

                    Path_Bits(1,k)=-1;

                    statecounter(1,k) = possible_pre_state;

                    break;

                else if((mod(2*(possible_pre_state-1),2^L)+1)==statecounter(1,k+1)-1)

                        Path_Bits(1,k)=1;

                        statecounter(1,k) = possible_pre_state;

                        break;

                    else

                        Path_Metric(possible_pre_state,k)=1e9;

                        [~, possible_pre_state] = min(Path_Metric(:,k));

                    end

                end

            end

        end

        

        

        %% Received Data Hard Decision

        

        r_NoChannel = sign(Y_NoChannel);

        r_GivenChannel_WithoutEq = sign(Y_GivenChannel_WithoutEq);

        r_GivenChannel_ZFLS = sign(Y_GivenChannel_ZFLS);

        r_GivenChannel_MMSE = sign(Y_GivenChannel_MMSE);

        r_GivenChannel_Viterbi = Path_Bits;



        %% Calculate Errors

        

        Error_NoChannel = sum(abs(r_NoChannel-s))/2;

        Error_GivenChannel_without_eq = sum(abs(r_GivenChannel_WithoutEq-s))/2;

         Error_GivenChannel_ZFLS = sum(abs(r_GivenChannel_ZFLS(Delay_of_channel:length(r_GivenChannel_ZFLS))-s(1:length(s)-Delay_of_channel+1)))/2;

         Error_GivenChannel_MMSE = sum(abs(r_GivenChannel_MMSE(Delay_of_channel:length(r_GivenChannel_MMSE))-s(1:length(s)-Delay_of_channel+1)))/2;

        Error_GivenChannel_Viterbi = sum(abs(r_GivenChannel_Viterbi(L+1:num_bit)-s(L+1:num_bit)))/2;



        %% Calcuate BER for each Transmission

        

        BER_NoChannel(count) = BER_NoChannel(count) + Error_NoChannel/num_bit;

        BER_GivenChannel_WithoutEQ(count) = BER_GivenChannel_WithoutEQ(count) + Error_GivenChannel_without_eq/num_bit;

        BER_GivenChannel_ZFLS(count) = BER_GivenChannel_ZFLS(count) + Error_GivenChannel_ZFLS/length(r_GivenChannel_ZFLS(Delay_of_channel:length(r_GivenChannel_ZFLS)));
        
        BER_GivenChannel_MMSE(count) = BER_GivenChannel_MMSE(count) + Error_GivenChannel_MMSE/length(r_GivenChannel_MMSE(Delay_of_channel:length(r_GivenChannel_MMSE)));
            
        BER_GivenChannel_Viterbi(count) = BER_GivenChannel_Viterbi(count) + Error_GivenChannel_Viterbi/(length(r_GivenChannel_Viterbi(L+1:num_bit)));



    end

end



%% BER per SNR

BER_th=0.5*erfc(sqrt(SNR));

BER_NoChannel =  BER_NoChannel/number_of_transmissions;

BER_GivenChannel_WithoutEQ =  BER_GivenChannel_WithoutEQ/number_of_transmissions;

BER_GivenChannel_ZFLS =  BER_GivenChannel_ZFLS/number_of_transmissions;

BER_GivenChannel_MMSE =  BER_GivenChannel_MMSE/number_of_transmissions;

BER_GivenChannel_Viterbi =  BER_GivenChannel_Viterbi/number_of_transmissions;





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

%% Plotting

if (keepPlot)

    hold on;

else

    hold off;

end





semilogy(SNRdB,BER_GivenChannel_WithoutEQ,'-d','color', green(4,:));

hold on;

semilogy(SNRdB,BER_GivenChannel_ZFLS,'-+','color',violett(4,:));

hold on;

semilogy(SNRdB,BER_GivenChannel_MMSE,'-*','color',blue(4,:));

hold on;

semilogy(SNRdB,BER_GivenChannel_Viterbi,'r-x','color',red(4,:));

hold on;


if(channel_model==0)

semilogy(SNRdB,BER_NoChannel,'-o','color',magenta(4,:));

hold on;

end
semilogy(SNRdB,BER_th,'--','color',cyan(4,:));

hold on;


grid on;



% Limits and Labels

xlim([0 max(SNRdB)]);

ylim([1e-6 1]);

xlabel('SNR=E_s/N_0 [dB]');

ylabel('BER');

if(channel_model==0)

legend('GivenChannel without EQ','ZFLS GivenChannel','MMSE GivenChannel','Viterbi GivenChannel','No Channel','AWGN BER theoretical','Location','southwest')

else

legend('GivenChannel without EQ','ZFLS RandomChannel','MMSE RandomChannel','Viterbi GivenChannel','AWGN BER theoretical','Location','southwest')

end

% Create an image file. This image is displayed on the webpage.

if (webdemo_en)

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

end

if (svg_en)

    print(['fig_',filename], '-dsvg',output_format);

end



end



% END, main program (control)

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


%% Function Definitions:
function b = de2bi(d, n, p, f)

  if (nargin == 1)
    p = 2;
    n = floor ( log (max (max (d), 1)) ./ log (p) ) + 1;
    f = 'right-msb';
  elseif (nargin == 2)
    p = 2;
    if (isstr(n))
      f = n;
      n = floor ( log (max (max (d), 1)) ./ log (p) ) + 1;
    else
      f = 'right-msb';
    endif     
  elseif (nargin == 3)
    if (isstr(n))
      f = n;
      n = p
      p = 2;
    elseif (isstr(p))
      f = p;
      p = 2;
    else
      f = 'right-msb';
    endif     
  elseif (nargin == 3)
    if (isstr(n))
      tmp = f
      f = n;
      n = p
      p = tmp;
    elseif (isstr(p))
      tmp = f;
      f = p;
      p = tmp;
    endif     
  else
    error ("usage: b = de2bi (d [, n [, p]])");
  endif

  d = d(:);
  if ( any (d < 0) || any (d != floor (d)) )
    error ("de2bi: only handles non-negative integers");
  endif

  power = ones (length (d), 1) * (p .^ [0 : n-1] );
  d = d * ones (1, n);
  b = floor (rem (d, p*power) ./ power);

  if (strcmp(f,'left-msb'))
    b = b(:,length(b):-1:1);
  elseif (!strcmp(f,'right-msb'))
    error("de2bi: unrecognized flag");
  endif
  
endfunction;

function d = bi2de (b, p, f)

  switch (nargin)
    case 1,
      p = 2;
      f = 'right-msb';
     case 2,
      if (isstr(p))
        f = p;
        p = 2;
      else
        f = 'right-msb';
      endif
     case 3,
      if (isstr(p))
        tmp = f;
        f = p;
        p = tmp;
      endif
   otherwise
      error ("usage: d = bi2de (b, [p])");
  endswitch

  if ( any (b (:) < 0) || any (b (:) > p - 1) )
    error ("bi2de: d must only contain value in [0, p-1]");
  endif

  if (strcmp(f,'left-msb'))
    b = b(:,size(b,2):-1:1);
  elseif (!strcmp(f,'right-msb'))
    error("bi2de: unrecognized flag");
  endif

  if (length (b) == 0)
    d = [];
  else
    d = b * ( p .^ [ 0 : (columns(b)-1) ]' );
  endif

endfunction;






