%%   Short Theoretical Background for the Function:%%   Program for baseband-simulation, using two orthogonal training sequences. %   Adds same training_sequence to all blocks. Simulates the transmitter for all blocks %   first then it simulates the receiver for all blocks.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     %%% Function part of simulation for Space-Time%%% coding project, group Grey-2001.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   Author: Stef and Fred%   Date: 2001-03-19%   Version: 1.0%   Revision (Name & Date):%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%close all, clear%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Initializations%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%N               = 5000;                   % Number of total datablock_size      = N/100;                   % block_size must be multiple of N.training_length1 = 50;                   % training_length + N must be an even number!training_length2 = training_length1;
fs              = 16000;                 % Sampling-frequencyT               = 0.0005;                  % Symboltimepulsetype       = 4;							% 1: rect 2:root-rais-cos 3:hamming 4:raised-cosinemodel		    = 1;							% 0: no Alamouti (1*1) 1: Alamouti (2*1)
rec             = 2;                            % 2: two receivers else one receiver     Eb              = 1;                            % Signalenergisigma2          = 1;gamma0dB        = 10;max_delay       = 880;                          % max length of delay vectorwindow_length   = 100;                          % must be equal to block_size + training_length1nr_plots        = 2;sync_factor     = 3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Running Simulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%w_length_pshaped = length(pulseshape(zeros(1,window_length),fs,pulsetype,T));a1err_mat=zeros(length(gamma0dB),N/block_size);a2err_mat=zeros(length(gamma0dB),N/block_size);t1err_mat=zeros(length(gamma0dB),N/block_size);t2err_mat=zeros(length(gamma0dB),N/block_size);compared_error = [];for(model=1)EbN0error = [];a1err_vec=[];a2err_vec=[];t1err_vec=[];t2err_vec=[];if model == 0     training_sequence1 = random_training(training_length1);else    [training_sequence1, training_sequence2]=make_orthotrain(training_length1);endfor(gamma=1:length(gamma0dB))    gamma0dB(gamma)    gamma0 = 10^(gamma0dB(gamma)/10);    N0 = 2*sigma2*Eb/gamma0;    error = [];     data = random_data(N);      % Generate random source-data        trans1 = [];    trans2 = [];    data_block_vec = [];    indata_est_vec = [];    for(i=0:(length(data)/block_size)-1)    i    data_block = get_datablock(data,block_size,i);    data_block_vec = [data_block_vec data_block];    symbols = bpsk(data_block,model);        [s_antenna1, s_antenna2] = alamouti(symbols,model);         % set = 1, alamouti coding is performed!        if model == 1        %[s_antenna1,s_antenna2,training_sequence1,training_sequence2]=add_training2(s_antenna1,s_antenna2,...        %    training_length1,training_length2,0);               % type = 0, random training-sequence.        %outdata_block = [training_sequence1 training_sequence2 data_block];        [s_antenna1,s_antenna2]=add_orthotrain_mult(s_antenna1,s_antenna2,training_sequence1,training_sequence2); % train1-train2 same length!                outdata_block = [xor(training_sequence1,training_sequence2) data_block]; % xor is just put to get the correct length of vector.            elseif model == 0        s_antenna1=add_training_mult(s_antenna1,training_sequence1);  % type = 0, random training-sequence.        outdata_block = [training_sequence1 data_block];        training_sequence2 = 0;    end % model == 1        %%%    %%% Antenna 1    %%%    
    [s1,t1]=pulseshape(s_antenna1,fs,pulsetype,T);          % Apply pulseshaping        %%%    %%% Antenna 2    %%%    
    [s2,t2]=pulseshape(s_antenna2,fs,pulsetype,T);          % Apply pulseshaping        trans1 = [trans1 s1];    trans2 = [trans2 s2];       end %block    [start_sequence,delayvec_length,sync] = add_start_seq(max_delay,window_length,pulsetype,fs,T,sync_factor);    trans1s = [start_sequence trans1];    trans2s = [zeros(1,length(start_sequence)) trans2];        %%%    %%% Simulating with simple channel 1 & 2 (3&4).    %%%    
    % noise only once !    %s1_thru_channel=s1;    [s1_thru_channel,alpha1,theta1]=rayleigh1(trans1s,sigma2);    %%%    %%%     %%%    [s1_thru_channel]=awgnbase(s1_thru_channel,Eb/N0); % Input argument is SNR in decibel.    %s2_thru_channel=s2;    [s2_thru_channel,alpha2,theta2]=rayleigh2(trans2s,sigma2);    %[s2_thru_channel]=awgnbase(s2_thru_channel,EbN0); % Input argument is SNR in decibel.     if rec == 2         [s1r2_thru_channel,alpha3,theta3]=rayleigh2(trans1s,1);        %[spb1_thru_channel]=channel(spb1,1,1);        %s1r2_thru_channel = s1;        [s1r2_thru_channel]=awgnbase(s1r2_thru_channel,Eb/N0); % Input argument is SNR in decibel.            %s2r2_thru_channel=s2;
        [s2r2_thru_channel,alpha4,theta4]=rayleigh1(trans2s,1);        %[s2r2_thru_channel]=awgnbase(s2r2_thru_channel,EbN0); % Input argument is SNR in decibel.    end % rec == 2                   %%%    %%% receiving in receive-antenna 1.    %%%        if model == 0        receive1 = s1_thru_channel;    else
        receive1 = s1_thru_channel + s2_thru_channel;    end % model == 0        ref_burst = ones(1,window_length);    [ref,t1]=pulseshape(ref_burst,fs,pulsetype,T);    treshold = mean(ref.^2)/6          %channel_mean = calculate_mean(window_length, receive1);           %receive1 = receive1 - channel_mean;         %Trigger         [n,b,w_SNR_vec] = trigger(receive1,sync,w_length_pshaped,treshold,fs,T);        delayvec_length    n    b_facit = delayvec_length + 4*w_length_pshaped + length(sync)    b        axis equal ;    subplot(nr_plots,1,1)           s = 0:w_length_pshaped:w_length_pshaped*(length(w_SNR_vec)-1);    stairs(s,w_SNR_vec)    subplot(nr_plots,1,2)    plot(abs(receive1))        receive1 = receive1(b:end);    if length(receive1) < length(trans1)        receive1 = [receive1 zeros(1,length(trans1)-length(receive1))];    elseif length(receive1) > length(trans1)        receive1 = receive1(1:length(trans1));    end                 if rec ==2        receive2 = s1r2_thru_channel + s2r2_thru_channel;                        [n2,b2,w_SNR_vec2] = trigger(receive2,sync,w_length_pshaped,treshold,fs,T);        delayvec_length;    n2;    b_facit = delayvec_length + 4*w_length_pshaped + length(sync);    b2        axis equal ;    subplot(nr_plots,1,1)           s2 = 0:w_length_pshaped:w_length_pshaped*(length(w_SNR_vec)-1);    stairs(s2,w_SNR_vec2)    subplot(nr_plots,1,2)    plot(abs(receive2))        receive2 = receive2(b2:end);    if length(receive2) < length(trans2)        receive2 = [receive2 zeros(1,length(trans2)-length(receive2))];    elseif length(receive2) > length(trans2)        receive2 = receive2(1:length(trans2));    end     end %if rec==2        % receiver1    for i = (1:w_length_pshaped:length(receive1)-w_length_pshaped+8)    quad = imag(receive1(i:i+w_length_pshaped-1));    inphase = real(receive1(i:i+w_length_pshaped-1));       % figure(2), subplot(2,2,1), plot(1:length(quad),quad);   % title('Quad after down-conv.');   % subplot(2,2,2), plot(1:length(inphase),inphase);   % title('Inphase after down-conv.');        [mf_quad_block]=matched_filter(fs, T, pulsetype, quad);    % Quadrature-part    [mf_inphase_block]=matched_filter(fs, T, pulsetype, inphase);    % Inphase-part        %figure(2)    %subplot(2,1,1), plot(1:2*training_length1,mf_quad_block(1:2*training_length1));    %title('Quad after MF');    %subplot(2,1,2), plot(1:2*training_length1,mf_inphase_block(1:2*training_length1));    %title('Inphase after MF');        %%%    %%% Performing synchronization