% Single-carrier equalization averaging, Institute of Telecommunications, University of Stuttgart, Germany

% www.inue.uni-stuttgart.de

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

% date: 05.7.2015

% Written for MATLAB

%

% EXAMPLE:  SC_average(3,2048, 20,15, 10,20, 2,0,0,1,filename)

% PARAMETERS:

%  channel_model                   % Mode of channel

%  num_bit,                       % number of bits per symbol

%  number_of_transmissions,       % number of transmissions over channel

%  order_of_equalizer             % Order of Equalizer

%  L                              % Order of channel

%  number_of_channels,			  % number of Channels to average over

%  deactivate_instable            % keep instable, remove some, remove all

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



% default parameters

switch nargin
    
    case 0
        
        channel_model = 3;
        deactivate_instable =2;
        num_bit = 2048;
        number_of_transmissions =20;
        order_of_equalizer =15;
        L=10;
        number_of_channels =40;
        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(number_of_channels>101 && number_of_channels<1)
    warning('Number of channels to average over needs to be in range 1-100')
end
if(deactivate_instable<0 && deactivate_instable>2)
    warning('Only 3 options, sorry pal')
end
if(channel_model<0 && channel_model<3)
    warning('Only 4 options, sorry pal')
end


% END check input parameters for consistency

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





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

% Definition



Eb=1;                                            %Bit energy

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

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

BER_Channel_ZFLS(length(SNR)) =0;

BER_Channel_MMSE(length(SNR)) =0;

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

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

i(Delay_of_channel) = 1;

stable=0;

%% Decay
% No decay channel
if(channel_model==1)
    if(deactivate_instable==0)
        %Energy already averaged
    end
    if(deactivate_instable==1)
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel;
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    h_Channel_Energy(u,:) = h_Channel;
                    if(isstable(1,h_Channel)<1)
                        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                        h_Channel_Energy(u,:) = h_Channel;
                    end
                end
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
    if(deactivate_instable==2)
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel;
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
end

%% L = 0.1 Energy
if(channel_model==2)
    alpha = - (1-sqrt(10))/L;
    if(deactivate_instable==0)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            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));
    end
    if(deactivate_instable==1)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                end
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
    if(deactivate_instable==2)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
end

% Decay L = 0.01 Energy

if(channel_model==3)
    alpha = - (1-10)/L;
    if(deactivate_instable==0)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            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));
    end
    if(deactivate_instable==1)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                    if(isstable(1,h_Channel)<1)
                        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                        h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                    end
                end
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
    if(deactivate_instable==2)
        k = 0:1:L;
        g_k = 1./(alpha*k+1).^2;
        for u=1:100
            h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
            end
        end
        beta = 1/u*sum(sum(abs(h_Channel_Energy).^2));
    end
end




%% Channels

for x=1:number_of_channels
    
    
    
    %% Given Channel
    
    clear h_GivenChannel;

    % Given  channel model
    
    if(channel_model==0)
        h_Channel(1,1)=1;
        for channel_stable_index=2:(L+1)
            h_Channel(1,channel_stable_index) = (-1)^channel_stable_index*1/(2^(channel_stable_index-1));
        end
        h_Channel = h_Channel./sqrt(h_Channel*h_Channel');
    end
    
    % Gaussian channel model
    if(channel_model==1)
        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
        if(deactivate_instable==1)
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    if(isstable(1,h_Channel)<1)
                        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    end
                end
            end
            h_Channel = 1/sqrt(beta)*h_Channel;
        end
        if(deactivate_instable==2)
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
            end
            h_Channel = 1/sqrt(beta)*h_Channel;
        end
        
    end
    
    % Channel model 2
    if(channel_model==2)
        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
        h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
        if(deactivate_instable==1)
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                end
            end
            h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
        end
        if(deactivate_instable==2)
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
            end
        end
        
    end
    
    % Channel model 3
    if(channel_model==3)
        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
        h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
        if(deactivate_instable==1)
            if(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                if(isstable(1,h_Channel)<1)
                    h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                    h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                    if(isstable(1,h_Channel)<1)
                        h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                        h_Channel_Energy(u,:) = h_Channel.*sqrt(g_k);
                    end
                end
            end
            h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
        end
        if(deactivate_instable==2)
            while(isstable(1,h_Channel)<1)
                h_Channel = normrnd(0,1/sqrt(L+1),[1 L+1]);
                h_Channel = 1/sqrt(beta)*h_Channel.*sqrt(g_k);
            end
        end
        
    end
    
    %% 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
    
    f_ZFLS_GivenChannel = (inv(ctranspose(H_Channel)*H_Channel))*ctranspose(H_Channel)*i; % Faster than H\i
    
    
    
    %GivenChannel 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]); %Generate AWGN, Normally distributed random number 1/2 for 2sigma^2
            
            
            
            
            
            %% Calculate Noise received Signals
            
            
            
            Y_Channel = filter(h_Channel, 1,s)+N;
            
            
            
            
            
            %% Linear Equalizing:
            
            
            
            Y_Channel_ZFLS = filter(f_ZFLS_GivenChannel, 1,Y_Channel);
            
            Y_Channel_MMSE = filter(f_MMSE_GivenChannel(:,count),1,Y_Channel);
            
            
            
            
            
            %% Received Data Hard Decision
            
            
            
            r_Channel_ZFLS = sign(Y_Channel_ZFLS);
            
            r_Channel_MMSE = sign(Y_Channel_MMSE);
            
            %% Calculate Errors
            
            Error_GivenChannel_ZFLS =  sum(abs(r_Channel_ZFLS(Delay_of_channel:length(r_Channel_ZFLS))-s(1:length(s)-Delay_of_channel+1)))/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_ZFLS(count) = BER_Channel_ZFLS(count) + Error_GivenChannel_ZFLS/length(r_Channel_ZFLS(Delay_of_channel:length(r_Channel_ZFLS)));
            BER_Channel_MMSE(count) = BER_Channel_MMSE(count) + Error_GivenChannel_MMSE/length(r_Channel_MMSE(Delay_of_channel:length(r_Channel_MMSE)));
            
        end
        
    end
    
    
    %% BER per SNR
    
    BER_Channel_ZFLS_averaged(x,:) =  BER_Channel_ZFLS/number_of_transmissions;
    
    BER_Channel_MMSE_averaged(x,:) =  BER_Channel_MMSE/number_of_transmissions;
    
    clear BER_Channel_ZFLS;
    
    clear BER_Channel_MMSE;
    
    BER_Channel_ZFLS(length(SNR))=0;
    
    BER_Channel_MMSE(length(SNR))=0;
    
end

BER_ZFLS = sum(BER_Channel_ZFLS_averaged)/number_of_channels;

BER_MMSE = sum(BER_Channel_MMSE_averaged)/number_of_channels;

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

%% Plotting

if (keepPlot)
    
    hold on;
    
else
    
    hold off;
    
end
semilogy(SNRdB,BER_ZFLS,'-','color',blue(3,:),'LineWidth',2);

hold on;

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

hold on;

semilogy(SNRdB,BER_th_rayleigh,'--','color',magenta(3,:));
hold on;
semilogy(SNRdB, BER_th, '--','color',red(3,:));

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('ZFLS Avg','MMSE Avg','Rayleigh theoretical','AWGN 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)

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


function [stable,p,z] = isstable(a,b)
    stable = 0;
    p = roots(b)
    if(abs(p)<=1)
        stable = 1;
    end
end


