% 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_N(0, 2048, 100, 5, 10, 1024,0,0,0,'OFDM_runN')

% PARAMETERS:

%  channel_model                  % Channel Modus

%  num_bit,                       % number of bits per symbol

%  number_of_transmissions,       % number of transmissions over channel

%  number_of_channels,		  % number of Channels to average over

%  L,		                  % number of Taps of Channel

%  Nmax  		          % Max Subcarrier

%  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 

%  filename                       % output file for webdemo


% BEGIN, main function (control)



function OFDM_N(channel_model, num_bit, number_of_transmissions, number_of_channels, L, Nmax,CP,Noisy_Estimate,keepPlot,svg_en,webdemo_en,filename)



	% default parameters

	switch nargin

	case 0

	channel_model =0;

	num_bit = 2048;

	number_of_transmissions =100;

	number_of_channels =15;

	L=10;

	Nmax = 1024;
    CP=0;
    Noisy_Estimate=0;
	hideLegend = 1;

	keepPlot = 0;

	svg_en = 0;

	webdemo_en = 0;

	filename = 'OFDM_runN';

	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;

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



% END, set up graphical output

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





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

% BEGIN check input parameters for consistency


if(num_bit<=0)

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

if(Nmax>num_bit)
   warning('Number of subcarriers need to be lower than bits')
end


if(number_of_transmissions>101 && number_of_transmissions<5)

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

end

if(Nmax>2048)

     warning('Number of subcarrier is too high')

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(Nmax>num_bit-1)
    num_bit = 2*Nmax;
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));



% Subcarriers ausgewählt
if(channel_model==1)
alpha = - (1-sqrt(10))/L;
    k = 0:1:L;
    g_k = 1./(alpha*k+1).^2;
    for u=1:1000
        h_GivenChannel = normrnd(0,1/sqrt(L+1),[1 L+1]);
        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)/L;
    k = 0:1:L;
    g_k = 1./(alpha*k+1).^2;
    for u=1:1000
        h_GivenChannel = normrnd(0,1/sqrt(L+1),[1 L+1]);
        h_Channel_Energy(u,:) = h_GivenChannel.*sqrt(g_k);
    end
    beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
end

subcarrier_array = [0 1 2 4 6 8 10];

if (ceil(log2(Nmax)>10))
    subcarrier_array = [subcarrier_array ceil(log2(Nmax))];
else
   subcarrier_array(floor(log2(Nmax)):end) =[];
end



% Prealocating for speed
BER_GivenChannel_OFDM(1000,length(SNR))=0;
BER_GivenChannel_OFDM_L(length(subcarrier_array),length(SNR))=0;



%% OFDM

% Loop over all Subcarriers

for x=1:length(subcarrier_array)

    Number_of_Data_Subcariers = 2^subcarrier_array(x);

    

    % Number of symbols

    OFDM_Symbols = num_bit./Number_of_Data_Subcariers;

    

    % 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(L+1),[1 L+1]);
        
        
        end

        % decay 0.1 model

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

         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(L+1),[1 L+1]);

        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:length(subcarrier_array)

    semilogy(SNRdB,BER_GivenChannel_OFDM_L(i,:),'-x','color',clr{i});

    hold on;

end

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


if (Nmax==64)
    legend('N_{sub}_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','Theory Rayleigh','location', 'southwest');

end
if (Nmax==256)
    legend('N_{sub}_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','N_{sub}=256','Theory Rayleigh','location', 'southwest');
end
if (Nmax==512)
    legend('N_{sub}_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','N_{sub}=256','N_{sub}=512','Theory Rayleigh','location', 'southwest');
end
if (Nmax==16)
    legend('N_{sub}_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','Theory Rayleigh','location', 'southwest');
end

if (Nmax==1024)

    legend('N_{sub}_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','N_{sub}=256','N_{sub}=1024','Theory Rayleigh','location', 'southwest');

end

if (2048>=Nmax>1024)

    legend('N_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','N_{sub}=256','N_{sub}=1024','2048','Theory Rayleigh','location', 'southwest');

end

if (4096>=Nmax>2048)

    legend('N_{sub}=1','N_{sub}=2','N_{sub}=4','N_{sub}=16','N_{sub}=64','N_{sub}=256','N_{sub}=1024','2048','Theory Rayleigh','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)

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









