%A.1 Main Sequence - main.m
%Cooperative Diversity - Main Sequencetic
% --------------
% Set Parametersnr_of_iterations = 10^3;
SNR = [-10:2.5:15];
use_direct_link = 1;
use_relay = 1;

global statistic;
%statistic = generate_statistic_structure;
global signal;
signal = generate_signal_structure;
signal(1).nr_of_bits = 2^10;
signal.modulation_type = 'QPSK'; % 'BPSK', 'QPSK'
calculate_signal_parameter;

channel = generate_channel_structure;
channel(1).attenuation(1).pattern = 'Rayleigh';% ’no’,’Rayleigh’
channel.attenuation.block_length = 1;
channel(2) = channel(1);
channel(3) = channel(1);
channel(1).attenuation.distance = 1;
channel(2).attenuation.distance = 1;
channel(3).attenuation.distance = 1;

rx = generate_rx_structure;
rx(1).combining_type = 'ESNRC'; %'ERC','FRC','SNRC','ESNRC','MRC'
rx(1).sd_weight = 3;

global relay;
relay = generate_relay_structure;
relay(1).mode = 'AAF'; %'AAF', 'DAF'
relay.magic_genie = 0;
relay(1).rx(1) = rx(1); % same beahaviour

% ----------------
% Start Simulation
BER = zeros(size(SNR));

for iSNR = 1:size(SNR,2)
    channel(1).noise(1).SNR = SNR(iSNR);
    channel(2).noise(1).SNR = SNR(iSNR);
    channel(3).noise(1).SNR = SNR(iSNR);
    
    disp(['progress: ',int2str(iSNR),'/',int2str(size(SNR,2))])
    
    for it = 1:nr_of_iterations;
        % --------------
        % Reset receiver
        rx = rx_reset(rx);
        relay.rx = rx_reset(relay.rx);
        % -----------
        % Direct link
        if (use_direct_link == 1)
            [channel(1), rx] = add_channel_effect(channel(1), rx,...
                signal.symbol_sequence);
            rx = rx_correct_phaseshift(rx, channel(1).attenuation.phi);
        end
        
        % ---------
        % Multi-hop
        if (use_relay == 1)
            % Sender to relay
            [channel(2), relay.rx] = add_channel_effect(channel(2),...
                relay.rx, signal.symbol_sequence);
                relay = prepare_relay2send(relay,channel(2));
                % Relay to destination
                [channel(3), rx] = add_channel_effect(channel(3), rx,...
                    relay.signal2send);
                
                switch relay.mode
                    % Correct phaseshift
                    case 'AAF'
                        rx = rx_correct_phaseshift(rx,...
                            channel(3).attenuation.phi + channel(2).attenuation.phi);
                    case 'DAF'
                        rx = rx_correct_phaseshift(rx,channel(3).attenuation.phi);
                end
            end
            
            % Receiver
            [received_symbol, signal.received_bit_sequence] = ...
                rx_combine(rx, channel, use_relay);
            
            BER(iSNR) = BER(iSNR) + sum(not(...
                signal.received_bit_sequence == signal.bit_sequence));
            
            if (BER(iSNR) > 10000)
                % Stop iterate
                break;
            end
        end % Iteration
        
        if (BER(iSNR)<100)
            warning(['Result might not be precise when SNR equal ',...
                    num2str(SNR(iSNR))])
        end
        BER(iSNR) = BER(iSNR) ./ it ./ signal.nr_of_bits;
    end
    % ------------------------------------
    % Present the result of the simulation
    txt_distance = [' - distance: ',...
            num2str(channel(1).attenuation.distance), ':',...
            num2str(channel(2).attenuation.distance), ':',...
            num2str(channel(3).attenuation.distance)];
    txt_distance='';
    if (use_relay == 1)
        if (relay.magic_genie == 1)
            txt_genie = ' - Magic Genie';
            elsetxt_genie = '';
        end
        
        txt_combining = [' - combining: ', rx(1).combining_type];
        switch rx(1).combining_type
            case 'FRC'
                txt_combining = [txt_combining, ' ',...
                        num2str(rx(1).sd_weight),':1'];
        end
        add2statistic(SNR,BER,[signal.modulation_type, ' - ',...
                relay.mode, txt_combining, txt_distance, txt_genie])
    else
        switch channel(1).attenuation.pattern
            case 'no'
                txt_fading = ' - no fading';
            otherwise
                txt_fading = ' - Rayleigh fading';
                enda
                dd2statistic(SNR,BER,[signal.modulation_type,txt_fading])
        end
        
        % % -----------------
        % % Graphs to compare
        SNR_linear = 10.^(SNR/10);
        % add2statistic(SNR,ber(SNR_linear,’BPSK’, ’Rayleigh’),’BPSK - single link transmission’)
        % add2statistic(SNR,ber_2_senders(SNR_linear, ’QPSK’),’QPSK - 2 senders’)
        show_statistic;
        toc
        
        %A.2 Initialise
        %A.2.1 Signal Parameter - calculate signal parameter.m
        %function calculate_signal_parameter
        % Calculates some additional signal parameters
        global signal;
        % Bits per symbol
        switch signal.modulation_type
            case 'BPSK'
                signal.bits_per_symbol = 1;
                
            case 'QPSK'
                signal.bits_per_symbol = 2;
                if (signal.nr_of_bits/2 ~= ceil(signal.nr_of_bits/2))
                    error(['Using QPSK, number of bits must be a multiple of 2'])
                end
            otherwise
                error(['Modulation-type unknown: ', signal.modulation_type])
        end
        
        % Number of symbols to transfer
        signal.nr_of_symbols = signal.nr_of_bits/signal.bits_per_symbol;
        
        % Bit sequence (random sequence of -1 and 1)
        signal.bit_sequence = floor(rand(1,signal.nr_of_bits)*2)*2-1;
        
        % Symbol sequence
        signal.symbol_sequence = bit2symbol(signal.bit_sequence);
        
        %A.2.2 Reset Receiver - rx reset.m
        function [rx] = rx_reset(rx);
        % Reset the receiver
        
        rx.signal2analyse = [];
        
        %A.3 Channel - add_channel_effect.m
        function [channel, rx] = add_channel_effect(channel,rx,signal_sequence)
        % Add noise fading and path loss
        
        global signal;
        %---------------------
        % Fading and path loss
        
        channel.attenuation.d = 1 / (channel.attenuation.distance ^ 2);
        
        % Path loss is constant for the whole transmission
        
        switch channel.attenuation.pattern
            case 'no'
                % No fading at all (only path loss)
                channel.attenuation.phi = zeros(size(signal_sequence));
                channel.attenuation.h = ones(size(signal_sequence)) * ...
                    channel.attenuation.d;
                channel.attenuation.h_mag = channel.attenuation.h;
                
            case 'Rayleigh'
                % Rayleigh fading and path loss
                nr_of_blocks = ceil(size(signal_sequence,2) /...
                    channel.attenuation.block_length);
                h_block = (randn(nr_of_blocks,1) + j * randn(nr_of_blocks...
                    ,1)) * channel.attenuation.d;
                
                h = reshape((h_block * ...
                    ones(1, channel.attenuation.block_length)), 1,...
                    channel.attenuation.block_length * nr_of_blocks);
                channel.attenuation.h = h(1:(size(signal_sequence,2)));
                
                [channel.attenuation.phi, channel.attenuation.h_mag] =...
                    cart2pol(real(channel.attenuation.h),...
                    imag(channel.attenuation.h));
                
                channel.attenuation.phi = -channel.attenuation.phi;
                
            otherwise
                error(['Fading-pattern unknown: ',...
                        channel.attenuation.pattern])
        end
        
        % ------------
        % Noise (AVGN)
        
        S = mean(abs(signal_sequence).^2);
        SNR_linear = 10^(channel.noise.SNR/10);
        
        %SNR = a^2/(2*sigma^2)
        channel.noise.sigma = sqrt(S / (2 * SNR_linear));
        noise_vector = (randn(size(signal_sequence)) +...
            j * randn(size(signal_sequence))) * channel.noise.sigma;
        % Add fading, path loss and noise to the signal
        rx.received_signal = signal_sequence .* channel.attenuation.h...
            + noise_vector;
        
        %A.4 Receiver
        %A.4.1 Correct Phase Shift - rx correct phaseshift.m
        function [rx] = rx_correct_phaseshift(rx, phi);
        % Correct phaseshift of the received signal
        
        switch rx.combining_type
            case 'MRC'
                % No phaseshift correction in MRC mode.
                % Phaseshift will be corrected when the received signal are
                % combinedrx.signal2analyse = [rx.signal2analyse; rx.received_signal];
                
            otherwise
                % Assuming that perfect phaseshift estimation possible
                rx.signal2analyse = [rx.signal2analyse;...
                        rx.received_signal .* exp(j * (phi))];
        end
        
        %A.4.2 Combine Received Signals - rx combine.m
        function [symbol_sequence, bit_sequence] = rx_combine(rx, channel, use_relay);
        % Combine all received signals
        global signal;
        global relay;
        values2analyse = rx.signal2analyse;
        if (use_relay == 1) & (relay.magic_genie == 1)
            switch relay.mode
                case 'DAF'
                    values2analyse(2,:) = (signal.symbol_sequence ==...
                        relay.symbol_sequence) .* values2analyse(2,:);
                otherwise
                    error(['Magic Genie works only with "DAF"'])
            end
        end
        switch rx.combining_type
            case 'MRC'
                switch relay.mode
                    case 'DAF'
                        if (use_relay == 0)
                            h = conj(channel(1).attenuation.h);
                        else
                            h = conj([channel(1).attenuation.h; 
                                channel(3).attenuation.h]);
                        end
                        bit_sequence = (mean(symbol2bit(h .*...
                            values2analyse),1)>=0)*2-1;
                        
                    otherwise
                        error('Maximum ratio combining works only with DAF')
                end
                
            case {'ERC', 'FRC', 'SNRC', 'ESNRC'}
                % The received values are already in phase
                values2analyse = symbol2bit(values2analyse);
                
                switch rx.combining_type
                    case 'ERC'
                        % Equal Ratio Combining
                        bit_sequence = (mean(values2analyse,1)>=0)*2-1;
                    case 'FRC'
                        % Fixed Ratio Combining
                        if (use_relay == 0)
                            bit_sequence = (mean(values2analyse,1)>=0)*2-1;
                        else
                            bit_sequence = (mean([rx.sd_weight;1] *...
                                ones(1,size(values2analyse,2)) .*...
                                values2analyse,1)>=0)*2-1;
                        end
                        
                    case {'SNRC', 'ESNRC'}
                        % Ratio depending on the SNR
                        if (use_relay == 0)
                            bit_sequence = (mean(values2analyse,1)>=0)*2-1;
                        else