SNR_direct = estimate_channel_SNR(channel(1), ...
                                signal.modulation_type, relay.mode);
                        SNR_via = estimate_channel_SNR([channel(2),...
                                channel(3)], signal.modulation_type, relay.mode);
                        if (signal.modulation_type == 'QPSK')
                            SNR_via = [SNR_via, SNR_via];
                            SNR_direct = [SNR_direct, SNR_direct];
                        end
                        
                        switch rx.combining_type
                            case 'SNRC'
                                bit_sequence_ratio = (sum([SNR_direct; SNR_via] .* ...
                                    values2analyse,1)>=0)*2-1;
                                bit_sequence_inf = (mean(values2analyse,1)>=0)*2-1;
                                SNR_equal_inf = ((SNR_via == inf) &...
                                    (SNR_direct == inf));
                                bit_sequence = SNR_equal_inf .* bit_sequence_inf +...
                                    not(SNR_equal_inf) .* bit_sequence_ratio;
                                
                            case 'ESNRC'
                                % .1 < SNR_direct/SNR_via < 10 : the to channels are
                                % weighted equally otherwise only the
                                % channel with the
                                % higher SNRR is used.
                                use_direct = (SNR_direct == inf) & (SNR_via ~= inf)...
                                    | ((SNR_direct ./ SNR_via) > 10);
                                use_via = (SNR_via == inf) & (SNR_direct ~= inf) | ...
                                    ((SNR_via ./ SNR_direct) > 10);
                                use_equal_ratio = not(use_direct + use_via);
                                
                                bit_sequence_equal_ratio =...
                                    (mean(values2analyse,1)>=0)*2-1;
                                bit_sequence_direct = (values2analyse(1,:)>=0)*2-1;
                                bit_sequence_via = (values2analyse(2,:)>=0)*2-1;
                                bit_sequence = ...
                                    use_equal_ratio .* bit_sequence_equal_ratio + ...
                                use_direct .* bit_sequence_direct + ...
                                use_via .* bit_sequence_via;
                        end
                    end
                    
                otherwise
                    error(['Combining-type unknown: ',rx.combining_type])
            end
    end
    
    symbol_sequence = bit2symbol(bit_sequence);
    
    %A.5 Relay - prepare relay2send.m
    function relay = prepare_relay2send(relay,channel);
    % Relay: Prepare received signal to make it ready to send
    global signal;
    switch relay.mode
        case 'AAF'
            % Amplify and Forward
            % Normalise signal power to the power of the original signal
            
            xi = abs(signal.symbol_sequence(1))^2;
            relay.amplification = sqrt(xi ./ (xi .*...
                channel(1).attenuation.h_mag .^ 2 + 2 .*...
                channel(1).noise.sigma .^ 2));
            relay.signal2send = ...
                relay.rx.received_signal .* relay.amplification;
            
        case 'DAF'
            % Decode and Forward
            relay.rx = rx_correct_phaseshift(relay.rx,...
                channel.attenuation.phi);
            relay.symbol_sequence = rx_combine(relay.rx, channel, 0);
            relay.signal2send = relay.symbol_sequence;
            
        otherwise
            error(['Unknown relay-mode: ', relay.mode])
    end
    
    %A.6 Structures
    %A.6.1 Signal - generate signal structure.m
    %function [signal_structure] = generate_signal_structure();
    % Creates the structure for all signal parameters
    
    signal_structure = struct(...
        'nr_of_bits',{},...% nr of bits to transfer
    'nr_of_symbols',{},...% nr of symbols to transfer
    'bits_per_symbol',{},...% BPSK (1 bit/symbol)
    ...% QPSK (2 bits/symbol)
    'modulation_type',{},...% ’BPSK’, ’QPSK’
    'bit_sequence',{},...% bit sequence of the signal
    'symbol_sequence',{},...% symbol sequence of the signal
    'received_bit_sequence',{});% bit sequence after transmission

%A.6.2 Channel - generate channel structure.m
%function [channel_structure] = generate_channel_structure();
% Creates the structure for all channel parameters
attenuation_structure = generate_attenuation_structure;
noise_structure = generate_noise_structure;
channel_structure = struct(...
    'attenuation', attenuation_structure,... % fading
    'noise', noise_structure);% noise

%function [fading_structure] = generate_attenuation_structure();
% Creates the structure for all fading parameters

fading_structure = struct(...
    'pattern',{},...% ’no’, ’Rayleigh’
    'distance', {},... % distance
    'd', {},...% path loss
    'h',{},...% attenuation incl. phaseshift
    'h_mag',{},...% magnitude
    'phi',{},...% phaseshift
    'block_length',{}); % lenth of the block (bit/block)

%function [noise_structure] = generate_noise_structure();
% Creates the structure for all noise parameters

noise_structure = struct(...
    'SNR',{},... % Signal to Noise Ratio (dB)
    'sigma',{}); % sigma of AVGN

%A.6.3 Receiver - generate rx structure.m
%function [rx_structure] = generate_rx_structure();
% Creates the structure for all receiver (Rx) parameters
rx_structure = struct(...
    'combining_type',{},... % ’ERC’, ’SNRC’, ’ESNRC’, ’MRC’
    'sd_weight',{},...% used for ’FRC’
                   ...% relay link is weighted one
    'received_signal',{},...% signal originally received. after
                          ...phaseshift is undone, saved in
                            ...signal2analyse
    'signal2analyse',{});% one row per incomming signal, which
                           % then are combined to estimate the
                           % bit-sequence
                %A.6.4 Relay - generate relay structure.m
                %function [relay_structure] = generate_relay_structure();
                % Creates the structure for all relay parameters
                rx_structure = generate_rx_structure;
                relay_structure = struct(...
                    'mode',{},...% ’AAF’ (Amplify and Forward)
                    ...              ’DAF’ (Decode and Forward)
                    'magic_genie',{},...% ’Magic Genie
                    'amplification',{},...% used in AAF mode
                    'symbol_sequence',{},...% used in DAF mode
                    'signal2send',{},...% Signal to be send
                    'rx',struct(rx_structure)); % Receiver\
                
                %A.6.5 Statistic - generate statistic structure.m
                %function [statistic_structure] = generate_statistic_structure();
                % Creates the structure for all statistic parameters
                statistic_structure = struct(...
                    'xlabel','SNR [dB]',...% label x-axis
                    'ylabel','Probability of error',... % label y-axis
                    'x',[],...% one graph per row x-axis
                    'y',[],...%y-axis
                    'legend','');
                %legend
                
                %A.7 ConversionsA.7.1 SNR to BER - ber2snr.m
                %function y = snr2ber(x)
                % Calculates the BER of the channel
                global signal;
                switch signal.modulation_type
                    case 'QPSK'
                        y = q(sqrt(x));
                    case 'BPSK'
                        y = q(sqrt(2 * x));
                    otherwise
                        error(['Modulation-type unknown: ', signal.modulation_type])
                end
                %A.7.2 BER to SNR - ber2snr.m
                %function y = ber2snr(x);
                % Calculates the SNR of the channel
                %
                % The SNR of the channel can be estimated/calculated when the
                % BER of the channel is known.
                global signal;
                switch signal.modulation_type
                    case 'QPSK'
                        y = qinv(x) .^ 2;
                    case 'BPSK'
                        y = qinv(x) .^ 2 / 2;
                    otherwise
                        error(['Modulation-type unknown: ', signal.modulation_type])
                end
                %A.7.3 Symbol Sequence to Bit Sequence - symbol2bit.m
                %function [bit_sequence] = symbol2bit(symbol_sequence);
                % Calculates bit_sequence from the symbol_sequence depending
                % on the modulation type
                global signal;
                switch signal.modulation_type
                    case 'BPSK'
                        bit_sequence = symbol_sequence;
                        
                    case 'QPSK'
                        bit_sequence = [real(symbol_sequence), imag(symbol_sequence)];
                        
                    otherwise
                        error(['Modulation-type unknown: ', signal.modulation_type])
                end
                %A.7.4 Bit Sequence to Symbol Sequence - bit2symbol.m
                %function [symbol_sequence] = bit2symbol(bit_sequence);
                % Calculates symbol_sequence from the bit_sequence depending on
                % the modulation type
                
                global signal;
                switch signal.modulation_type
                    case 'BPSK'
                        symbol_sequence = bit_sequence;
                        
                    case 'QPSK'
                        symbol_sequence = bit_sequence(1:signal.nr_of_symbols) + j*...
                            bit_sequence(signal.nr_of_symbols + 1 : signal.nr_of_bits);
                    otherwise
                        error(['Modulation-type unknown: ', signal.modulation_type])
                end
                %A.8 Statistic
                %A.8.1 Add Statistic - add2statistic.m
                %function add2statistic(x,y,leg);
                % Add graph to statistic
                global statistic;
                
                statistic.x = [statistic.x;x];
                statistic.y = [statistic.y;y];
                statistic.legend = strvcat(statistic.legend,leg);
                %A.8.2 Show Statistic - show statistic.m
                %function [handle] = show_statistic(colour_bw, order);
                % Shows the result in a plot
                global statistic;
                if (nargin<1), colour_bw = 0;
                end
                if (nargin<2), order = 1:size(statistic.x,1); end
                if (colour_bw == 1)
                    colours = ['k-o';'k-*';'k-s';'k-+';'k-^';'k-h';'k-v';'k-p'];
                else
                    colours = ['b-o';'r-d';'g-s';'k-v';'m-^';'b-<';'r->';'g-p'];
                end
                legend_ordered = [];
                handle = figure;
                colour = 0;
                for n = order
                    colour = colour + 1;
                    semilogy(statistic.x(n,:),statistic.y(n,:),colours(colour,:));
                    legend_ordered = strvcat(legend_ordered,statistic.legend(n,:));
                    hold on
                end
                grid on;
                legend (legend_ordered,3)
                xlabel (statistic.xlabel)
                ylabel (statistic.ylabel)
                %A.9 Theoretical BER
                %A.9.1 Single Link Channel - ber.m
                %function [y] = ber(snr, modulation_type, fading_type);
                % Calculates the BER(SNR) depending on the modulation-type and
                % the fading-type
                
                switch fading_type
                    case 'Rayleigh'
                        switch modulation_type
                            case 'BPSK'
                                y = (1 - sqrt(snr ./ (1 / 2 + snr))) / 2;
                            case 'QPSK'
                                y = (1 - sqrt(snr ./ (1 + snr))) / 2;
                                
                            otherwise
                                error(['Modulation-type unknown: ', modulation_type])
                        end
                    case 'no'
                        switch modulation_type
                            case 'BPSK'
                                y = q(sqrt(2 * snr));
                            case 'QPSK'
                                y = q(sqrt(snr));
                            otherwise
                                error(['Modulation-type unknown: ', modulation_type])
                        end
                    otherwise
                        error(['Fading-type unknown: ', fading_type])
                end
                %A.9.2 Two Independent Senders - ber 2 senders.m
                %function y = ber_2_senders(SNR_avg, modulation_type);
                % BER(SNR) using two senders. The (average) SNR is assumed to be
                % equal for both channel
                switch modulation_type
                    case 'BPSK'
                        mu = sqrt(SNR_avg ./ (1 / 2 + SNR_avg));
                    case 'QPSK'
                        mu = sqrt(SNR_avg ./ (1 + SNR_avg));
                    otherwise
                        error(['Modulation-type unknown: ', modulation_type])
                end
                
                y = 1 / 4 * (1 - mu) .^ 2 .* (2 + mu);
                %A.10 Math functions
                %A.10.1 Q-function - q.m
                function [y] = q(x);
                % Q-probability function
                y = erfc(x / sqrt(2)) / 2;
                
                %A.10.2 Inverse Q-function - invq.m
                function [y] = qinv(x);
                % Inverse Q-probability function
                y = erfcinv(x*2) *sqrt(2);