% Comparision of Single-carrier modulati, Institute of Telecommunications, University of Stuttgart, Germany



% www.inue.uni-stuttgart.de



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



% date: 05.07.2015





% Written for MATLAB



%



% EXAMPLE: SCvsMC(2048, 20,15, 10,50,0.01,0,0,1,filename)



% PARAMETERS:



%  num_bit                        % Number of bit per transmission



%  number_of_transmissions        % Number of transmissions



%  order_of_equalizer             % Order of Equalizer



%  L                              % Order of channel



%  number_of_channels             % Number of channels to average over



%  BER_target                     % Target Biterrorratio



%  hideLegend                     % Hide Legend



%  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 SCvsMC(num_bit, number_of_transmissions,order_of_equalizer, L,number_of_channels,BER_target,keepPlot,svg_en,webdemo_en,filename)



% default parameters



switch nargin

    

    case 0

        num_bit = 2048;                  % Number of bit per transmission

        number_of_transmissions =20;     % Number of transmissions

        order_of_equalizer =15;          % Order of Equalizer

        L=10;                            % Order of channel

        number_of_channels =100;         % Number of channels to average over

        BER_target = 0.01;               % Target Biterrorratio

        hideLegend = 1;                  % Hide Legend

        keepPlot = 0;                    % Keep plot

        svg_en = 0;                      % Debug option

        webdemo_en = 0;                  % Debug option

        filename = 'SCvsMC';             % filename

        

    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<2047)
    num_bit = 4096;
end


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(number_of_channels>101 && number_of_channels<1)

     warning('Number of channels to average over needs to be in range 1-100')

end

if(BER_target>0.5 && BER_target<0.01)

    warning('BER needs to be in range 0.01-0.5')

end



% END check input parameters for consistency



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











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



%% Definitions

%SNR

Eb=1;                                            %Bit energy

SNRdB=0:0.5:35.5;                                  %Signal to Noise Ratio (in dB)

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



% number of bits tranmission



num_bit = 2^(ceil(log2(num_bit)));



% BER allocation for speed

BER_Channel_MMSE(length(SNR)) =0;

BER_Channel_OFDM(length(SNR))=0;

BER_MMSE(length(SNR))=0;

BER_ZFLS(length(SNR))=0;

BER_OFDM(length(SNR))=1;



% BER target vectorization

BER_target = ones(1,length(SNR))*BER_target;



% Received vector wish

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

Delay_of_channel = ceil((L+1)/2);

i(Delay_of_channel) = 1;



% Timer initialization

timer_OFDM=0;

timer_MMSE=0;



%% Decay

% Decay L = 0.01 Energy

alpha = - (1-10)/L;

k = 0:1:L;

g_k = 1./(alpha*k+1).^2;

for u=1:1000

    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);

    h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);

end

beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));



%% Subcarrier prechoice

if(BER_target(1)>0.1)

subcarrier_array =[4 8 16 32 64 128 256 512 1024];

end

if(0.05<BER_target(1)<=0.1)

subcarrier_array =[32 64 128 256 512 1024];   

end

if(BER_target(1)<=0.05)

subcarrier_array =[256 512 1024];   

end

index=1;



%% Channels

while (BER_OFDM(length(SNR))-0.5*BER_target(length(SNR))>0)

    number_of_subcarriers = subcarrier_array(index);

    % Number of symbols

    OFDM_Symbols = num_bit./number_of_subcarriers;

    

    for x=1:number_of_channels

        %% Given Channel

        clear h_Channel;

        % Create channel

        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);

        h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);

        

        %% Calculation Coeffiencents of Filters

        

        %GivenChannel ZFLS

        H_Channel = 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_Channel(stable_matrix_j,stable_matrix_i) = h_Channel(1,k);

                

                k=k+1;

                

            end

            

            stable_matrix_j = stable_matrix_j+1-L;

            

        end
       
              

        %% Filtercoef MMSE     

        for count=1:length(SNR)      

            f_MMSE_GivenChannel(:,count) = (inv(ctranspose(H_Channel)*H_Channel+(1/(SNR(count)))*eye(order_of_equalizer+1, order_of_equalizer+1)))*ctranspose(H_Channel)*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 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]);

                N_OFDM = normrnd(0,sqrt(No/2),[1,num_bit])+1i*normrnd(0,sqrt(No/2),[1,num_bit]);

                

                

                %% ODFM Transmitter

                

                %Reshape into symbols

                data_s = reshape(s,number_of_subcarriers,OFDM_Symbols);

                

                % First Symbol as Pilot(Only ones)

                data_s(:,1)= abs(data_s(:,1));

                

                % iFFT

                tic

                ofdm_normalized_ifft = sqrt(number_of_subcarriers)*ifft(data_s);

                timer_OFDM = timer_OFDM +toc/(OFDM_Symbols);

                

                % Concatenating multiple symbols to form a long vector

                ofdm_normalized_ifft = reshape(ofdm_normalized_ifft,OFDM_Symbols*number_of_subcarriers,1);

                

                %% Calculate Noise received Signals

                Y_Channel = filter(h_Channel, 1,s)+N;

                Y_Channel_OFDM = filter(h_Channel.', 1,ofdm_normalized_ifft)+N_OFDM';

                test = filter(h_Channel.', 1,ofdm_normalized_ifft);

                

                

                %% OFDM Receiver

                

                % Receiver reshape into Symbols

                r_GivenChannel_moded_cycled_OFDM = reshape(Y_Channel_OFDM,number_of_subcarriers,OFDM_Symbols); % formatting the received vector into symbols

                test = reshape(test,number_of_subcarriers,OFDM_Symbols);

                

                

                % FFT

                tic

                r_GivenChannel_moded_OFDM = fft(r_GivenChannel_moded_cycled_OFDM)*1/sqrt(number_of_subcarriers);

                timer_OFDM = timer_OFDM +toc/(OFDM_Symbols);

                test = fft(test)*1/sqrt(number_of_subcarriers);

                

                % Reshape intto long Vector

                r_GivenChannel_moded_OFDM = reshape(r_GivenChannel_moded_OFDM ,number_of_subcarriers*OFDM_Symbols,1);

                test = test(:,1)*ones(1,OFDM_Symbols);

                test = reshape(test,number_of_subcarriers*OFDM_Symbols,1);

                

                % Single Tap Equalizer

                tic

                r_GivenChannel_moded_OFDM = r_GivenChannel_moded_OFDM.*conj(test);

                timer_OFDM = timer_OFDM +toc/(OFDM_Symbols);

                

                %% Linear Equalizing:            

                 tic

                Y_Channel_MMSE = filter(f_MMSE_GivenChannel(:,count),1,Y_Channel);

                timer_MMSE = timer_MMSE+toc;

                

                %% Received Data Hard Decision

                r_GivenChannel_OFDM = sign(real(r_GivenChannel_moded_OFDM.'));

                 r_Channel_MMSE = sign(Y_Channel_MMSE);

                

                %% Calculate Errors

                Error_GivenChannel_OFDM = sum(abs(r_GivenChannel_OFDM((num_bit/OFDM_Symbols+1):num_bit)-s((num_bit/OFDM_Symbols+1):num_bit)))/2;
            
                Error_GivenChannel_MMSE = sum(abs(r_Channel_MMSE(Delay_of_channel:length(r_Channel_MMSE))-s(1:length(s)-Delay_of_channel+1)))/2;

                

                %% Calcuate BER for each Transmission
             
                BER_Channel_MMSE(count) = BER_Channel_MMSE(count) + Error_GivenChannel_MMSE/length(r_Channel_MMSE(Delay_of_channel:length(r_Channel_MMSE)));

                BER_Channel_OFDM(count) = BER_Channel_OFDM(count) + Error_GivenChannel_OFDM/length((num_bit/OFDM_Symbols+1):num_bit);

            end

        end

        

        

        %% BER per transmission 

        BER_Channel_MMSE_averaged(x,:) =  BER_Channel_MMSE/number_of_transmissions;

        BER_Channel_OFDM_averaged(x,:) =  BER_Channel_OFDM/number_of_transmissions;
      
        clear BER_Channel_MMSE;

        clear BER_Channel_OFDM;

        BER_Channel_MMSE(length(SNR))=0;

        BER_Channel_OFDM(length(SNR))=0;

       

    end

    BER_MMSE = sum(BER_Channel_MMSE_averaged)/number_of_channels;

    BER_OFDM = sum(BER_Channel_OFDM_averaged)/number_of_channels;

    index = index+1;

end



% Theoretic Rayleigh

BER_th_rayleigh=0.5*(1-sqrt((SNR)./(SNR+1)));



% Intersection finding

intersection_MMSE = BER_MMSE-BER_target;

intersection_OFDM = BER_OFDM-BER_target;

ix = find(intersection_MMSE <0);
if(isempty(ix))
   if( BER_MMSE(length(SNR))-BER_target>0.01)
      ix = 2*length(SNR);
   else 
   ix = length(SNR);
   end
end
iy = find(intersection_OFDM <0);

%Flops counting

flops_STE = number_of_subcarriers;

flops_OFDM = 2*number_of_subcarriers*log2(number_of_subcarriers)+flops_STE;

flops_MMSE = number_of_subcarriers^2;



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



%% Plotting

if (keepPlot)

    hold on;

else

    hold off;

end



semilogy(SNRdB,BER_MMSE,'-','color', green(3,:),'LineWidth',3);

hold on;

semilogy(SNRdB,BER_OFDM,'-','color', cyan(3,:),'LineWidth',3);

hold on;
semilogy(SNRdB,BER_th_rayleigh,'--','color',magenta(1,:));
hold on;

semilogy(SNRdB,BER_target,'--','color', red(3,:));

hold on;

semilogy([SNRdB(ix(1)) SNRdB(ix(1))],[1e-6 1],'--','color', red(3,:));

hold on;

semilogy([SNRdB(iy(1)) SNRdB(iy(1))],[1e-6 1],'--','color', red(3,:));

hold on;

%% Add text

if(SNRdB(ix(1))-SNRdB(iy(1))<6)

    text(SNRdB(ix(1))+0.5,0.3,['SNR_{req,SC}=' num2str(SNRdB(ix(1))) 'dB'],'FontSize',8,'EdgeColor','g');

else

    text(SNRdB(ix(1))-0.5,0.3,['SNR_{req,SC}= ' num2str(SNRdB(ix(1))) 'dB'],'FontSize',8,'EdgeColor','g','HorizontalAlignment','right');

end


text(SNRdB(iy(1))-0.5,0.3,['SNR_{req,MC}=' num2str(SNRdB(iy(1))) 'dB'],'FontSize',8,'EdgeColor','g','HorizontalAlignment','right');

text(1,0.2e-4,['\Delta SNR=' num2str(SNRdB(ix(1))-SNRdB(iy(1))) 'dB'],'FontSize',8,'EdgeColor','g');

text(1,0.8e-5, ['NOPS_{MC/SC}='   num2str(round(100*flops_OFDM/flops_MMSE)/100)],'FontSize',8,'EdgeColor','g');

text(1,0.3e-5,['t_{MC}/t_{SC}=' num2str(round(100*timer_OFDM/timer_MMSE)/100)],'FontSize',8,'EdgeColor','g');

grid on;







% Limits and Labels



xlim([0 max(SNRdB)]);



ylim([1e-6 1]);



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



ylabel('BER');



legend('MMSE',['OFDM N_{sub}=' num2str(subcarrier_array(index))],'Rayleigh theoretical','Location','southeast')







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



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














