 %BER OFDM over different Subcarrier, Institute of Telecommunications, University of Stuttgart, Germany

% www.inue.uni-stuttgart.de

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

% date: 27.6.2015

% update: 27.6.2015

%

% Written for MATLAB

%

% EXAMPLE:  OFDM_L( 0,2048, 100, 5, 10, 1024,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

%  number_of_subcarrier           % number of subcarrier

%  number_of_channels,		  % number of Channels to average over

%  Lmin,		          % lowest number of Taps of Channel

%  Lmax,			  % highest number of Taps 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 OFDM_L(channel_model,num_bit, number_of_transmissions,number_of_subcarrier, number_of_channels, Lmin, Lmax,keepPlot,svg_en,webdemo_en,filename)



% default parameters

switch nargin

    case 0
        channel_model =0;
        num_bit = 2048;
        number_of_transmissions =100;
        Number_of_Data_Subcariers = 64;
        number_of_subcarrier= Number_of_Data_Subcariers; 
        number_of_channels =10;
        Lmin=4;
        Lmax = 10;
        hideLegend = 1;
        keepPlot = 0;
        svg_en = 0;
        webdemo_en = 0;
        filename = 'OFDM_runN';
    otherwise
        
end
CP=0;
Noisy_Estimate =0;
         Number_of_Data_Subcariers = number_of_subcarrier; 

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

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

clr  ={red(2,:), blue(2,:), green(2,:), magenta(4,:), cyan(4,:),  violett(4,:), seagreen(6,:), yellow(4,:) amber(4,:) ,lightgreen(4,:), seagreen(4,:)};




% 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_subcarrier>2048)
   warning('Number of subcarriers is too high')
end
if(number_of_subcarrier>num_bit-1)
    num_bit = 2*number_of_subcarrier;
end

if(number_of_transmissions>101 && number_of_transmissions<5)

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

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 (Lmin<1)
   warning('Sorry L is not allowed to be smaller than 1')
end
if (Lmax>20) 
   warning('Sorry L is not allowed to be over 20, if needed download source code and use your own matlab')
end
if (Lmax<Lmin);

    warning('Lmax needs to be bigger than Lmin');
end

% END check input parameters for consistency

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





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

%Bit energy

Eb=1;

%Signal to Noise Ratio (in dB)

SNRdB=0:5:35.5;

%Signal to Noise Ratio

SNR=10.^(SNRdB/10);

% Num Bits auf 2^Potenz bringen:

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





% Prealocating for speed

BER_GivenChannel_OFDM(1000,length(SNR))=0;

BER_GivenChannel_OFDM_L((Lmax-Lmin)+1,length(SNR))=0;

%Decay 

%% OFDM

% Loop over all Subcarriers

for x=1:(Lmax-Lmin)+1
    %Decay 
    clear h_Channel_Energy;
    if(channel_model==1)
    alpha = - (1-sqrt(10))/(Lmin+x-1);
    k = 0:1:(Lmin+x-1);
    g_k = 1./(alpha*k+1).^2;
    for u=1:1000
        h_GivenChannel = normrnd(0,1/sqrt(Lmin+x),[1 Lmin+x]);
        h_Channel_Energy(u,:) = h_GivenChannel.*sqrt(g_k);
    end
    beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
end

if(channel_model==2)
    alpha = - (1-10)/(Lmin+x-1);
    k = 0:1:Lmin+x-1;
    g_k = 1./(alpha*k+1).^2;
    for u=1:1000
        h_GivenChannel = normrnd(0,1/sqrt(Lmin+x),[1 Lmin+x]);
        h_Channel_Energy(u,:) = h_GivenChannel.*sqrt(g_k);
    end
    beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end

    

    % Number of symbols

    OFDM_Symbols = num_bit./number_of_subcarrier;

    

    % Loop over

    for i=1:number_of_channels

        

        %% Given Channel

        clear h_GivenChannel;

       % Gaussian channel model

        if(channel_model==0)

        h_GivenChannel = normrnd(0,1/sqrt(Lmin+x),[1 Lmin+x]);


        end

        % decay 0.1 model

         if(channel_model==1)
         h_GivenChannel = normrnd(0,1/sqrt(Lmin+x),[1 Lmin+x]);

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

         end

       % decay 0.01  model

        if(channel_model==2)
        h_GivenChannel = normrnd(0,1/sqrt(Lmin+x),[1 Lmin+x]);

        h_GivenChannel = 1/sqrt(beta)*h_GivenChannel.*sqrt(g_k);
        end




        

        %% Calculating BER

        

        % Sending Packets over Channel

        

        for transmission_counter=1:number_of_transmissions

            for count=1:length(SNR)
                %% Modulate the Signal and Preconditions
                %Generate binary data source
                data=randi([0,1],1,num_bit);
                
                %BPSK modulation
                s=2*data-1;
                
                %Generate Noise
                No=Eb/SNR(count);
                N_OFDM = normrnd(0,sqrt(No/2),[1,num_bit+OFDM_Symbols* floor(Number_of_Data_Subcariers*CP)])+1i*normrnd(0,sqrt(No/2),[1,num_bit+OFDM_Symbols* floor(Number_of_Data_Subcariers*CP)]);
                               
                if(Noisy_Estimate==0)
                    % Remove Noise for Pilot to keep Rayleigh Curve correct
                    N_OFDM(floor(Number_of_Data_Subcariers*(1+CP))+1:2*(floor(Number_of_Data_Subcariers*(1+CP)))) =0;
                end
                
                %% ODFM Transmitter
                %Reshape into symbols
                data_s = reshape(s,Number_of_Data_Subcariers,OFDM_Symbols);
                
                % iFFT
                if (Number_of_Data_Subcariers==1)
                    ofdm_normalized_ifft = sqrt(num_bit)*ifft(data_s);
                else
                    ofdm_normalized_ifft = sqrt(Number_of_Data_Subcariers)*ifft(data_s);
                end
                % Add CP
                ofdm_normalized_ifft = vertcat(ofdm_normalized_ifft(end-floor(CP*Number_of_Data_Subcariers)+1:end,:), ofdm_normalized_ifft);
                % Concatenating multiple symbols to form a long vector
                ofdm_normalized_ifft = reshape(ofdm_normalized_ifft,OFDM_Symbols*(floor(Number_of_Data_Subcariers*(1+CP))),1);
                
                %% OFDM
                Y_Given_Channel_OFDM = filter(h_GivenChannel.', 1,ofdm_normalized_ifft)+N_OFDM';
                
                %% OFDM Receiver
                % Receiver reshape into Symbols
                r_GivenChannel_moded_cycled_OFDM = reshape(Y_Given_Channel_OFDM,floor(Number_of_Data_Subcariers*(1+CP)),OFDM_Symbols); % formatting the received vector into symbols
                % Remove CP
                r_GivenChannel_moded_cycled_OFDM(1:CP*Number_of_Data_Subcariers,:)=[];
                % FFT
                r_GivenChannel_moded_OFDM = fft(r_GivenChannel_moded_cycled_OFDM)*1/sqrt(Number_of_Data_Subcariers);
                
                % Reshape intto long Vector
                Channel_Equaliser = (r_GivenChannel_moded_OFDM(:,2)./transpose(s(Number_of_Data_Subcariers+1:2*Number_of_Data_Subcariers)))*ones(1,OFDM_Symbols);
                Channel_Equaliser = reshape(Channel_Equaliser,Number_of_Data_Subcariers*OFDM_Symbols,1);
                
                r_GivenChannel_moded_OFDM = reshape(r_GivenChannel_moded_OFDM ,Number_of_Data_Subcariers*OFDM_Symbols,1);
                % Single Tap Equalizer
                r_GivenChannel_moded_OFDM = r_GivenChannel_moded_OFDM./Channel_Equaliser;
                
                %% Received Data Hard Decision
                r_GivenChannel_OFDM = sign(real(r_GivenChannel_moded_OFDM.'));
                
                %% Calculate Errors
                Error_GivenChannel_OFDM = sum(abs(r_GivenChannel_OFDM(2*(num_bit/OFDM_Symbols+1):num_bit)-s(2*(num_bit/OFDM_Symbols+1):num_bit)))/2;
                
                %% Calcuate BER for each Transmission
                BER_GivenChannel_OFDM(i,count) = BER_GivenChannel_OFDM(i,count) + Error_GivenChannel_OFDM/length((num_bit/OFDM_Symbols+1):num_bit);
            end
        end

        

        %% BER per SNR

        BER_GivenChannel_OFDM(i,:) =  BER_GivenChannel_OFDM(i,:)/number_of_transmissions;

    end

    BER_GivenChannel_OFDM_L(x,:) = sum(BER_GivenChannel_OFDM)/number_of_channels;

end



%% Plotting

if (keepPlot)

    hold on;

else

    hold off;

end



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



for i=1:(Lmax-Lmin)+1

    liste{Lmax-Lmin+2-i}=['Channel order L=' int2str(i-1+Lmin)]; 

    semilogy(SNRdB,BER_GivenChannel_OFDM_L(Lmax-Lmin+2-i,:),'-o','color', clr{i});

    hold on;

end

liste{Lmax-Lmin+2}=['Theory Rayleigh'];

semilogy(SNRdB,BER_th_rayleigh,'-','color', gray(1,:));

hold on;



grid on;



% Limits and Labels

xlim([0 max(SNRdB)]);

ylim([1e-6 1]);

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

ylabel('BER');

legend(liste{:},'Location', 'southwest');

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

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









