www.pudn.com > CooperativeDiversity.rar > CooperativeDiversity.m


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