%%   Short Theoretical Background for the Function:%%   Program for baseband-simulation, using two orthogonal training sequences.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     %%% 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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%block_size      = 30;                  Nr_of_blocks    = 1;N               = Nr_of_blocks*block_size;                   % Number of total datatraining_length1 = 20;                   % 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             = 1;                            % 2: two receivers else one receiver     Eb              = 1;                            % Signalenergisigma2          = 1;gamma0dB        = 20;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Running Simulation%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
compared_error = [];for (rec=1)    for(model=0)         if (model~=0 | rec~=2)            EbN0error = []                        for(gamma=1:length(gamma0dB))                    gamma0dB(gamma)                gamma0 = 10^(gamma0dB(gamma)/10);                N0 = 2*sigma2*Eb/gamma0;                error = zeros(1,Nr_of_blocks);                    data = random_data(N);      % Generate random source-data                   for(i=0:Nr_of_blocks-1)                                               data_block = get_datablock(data,block_size,i);                    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_orthotrain(s_antenna1,s_antenna2,training_length1); % 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,training_sequence1]=add_training(s_antenna1,training_length1,0);  % 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                               %%%                    %%% Simulating with simple channel 1 & 2 (3&4).                    %%%    
                    % noise only once !                                        [s1_thru_channel,alpha1,theta1]=rayleigh(s1,sigma2);                    [s1_thru_channel]=awgnbase(s1_thru_channel,Eb/N0); % Input argument is SNR in decibel.                                        [s2_thru_channel,alpha2,theta2]=rayleigh(s2,sigma2);                        if rec == 2                                 [s1r2_thru_channel,alpha3,theta3]=rayleigh(s1,sigma2);                        [s1r2_thru_channel]=awgnbase(s1r2_thru_channel,Eb/N0); % Input argument is SNR in decibel.                                                [s2r2_thru_channel,alpha4,theta4]=rayleigh(s2,sigma2);                    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                             quad = imag(receive1);                    inphase = real(receive1);                            [mf_quad_block]=matched_filter(fs, T, pulsetype, quad);    % Quadrature-part                    [mf_inphase_block]=matched_filter(fs, T, pulsetype, inphase);    % Inphase-part                                              %%%                    %%% Performing synchronization                    %%%                        [corr_maxi,sample_nri,sum_of_r_times_c1i,sum_of_r_times_d2i]=orthosynchronization(mf_inphase_block, fs, T, training_sequence1, training_sequence2, model);	                [corr_maxq,sample_nrq,sum_of_r_times_c1q,sum_of_r_times_d2q]=orthosynchronization(mf_quad_block, fs, T, training_sequence1, training_sequence2, model);                                %if corr_maxq > corr_maxi       % choose the correlation with the biggest max value                    %    sample_nr = sample_nrq;                    %else                    %    sample_nr = sample_nri;                    %end                        sample_nr = 81;                     block_length=length(outdata_block);                    [inphase_symbols]=down_sampler(mf_inphase_block,sample_nr,fs,T,block_length);                    [quad_symbols]=down_sampler(mf_quad_block,sample_nr,fs,T,block_length);     
	                %%%                    %%% Estimating Channel 1 and 2.                    %%%
                          %[alpha1est,theta1est] = channel_estimator_(quad_symbols(1:training_length1), ...                     %   inphase_symbols(1:training_length1), training_sequence1, model, 1);                    if model == 1
                          %  [alpha2est,theta2est]=channel_estimator_(quad_symbols(1:training_length1),...                       %     inphase_symbols(1:training_length1),training_sequence2,model,1);                            elseif model == 0                        alpha2est = 1;                        theta2est = 0;                    end % model==1                           % set true values instead:                    alpha1est = alpha1;                    theta1est = theta1;                    alpha2est = alpha2;                    theta2est = theta2;                                                %a1err=alpha1est-alpha1                    %a2err=alpha2est-alpha2                    %t1err=theta1est-theta1                    %t2err=theta2est-theta2
                                       %a1err_vec=[a1err_vec a1err];                    %a2err_vec=[a2err_vec a2err];                       %t1err_vec=[t1err_vec t1err];                    %t2err_vec=[t2err_vec t2err];                                                                                %%%
                    %%% The combiner for receiver1
                    %%%
                            [combined_signal]=combiner(alpha1est, theta1est, alpha2est, theta2est, ...                        inphase_symbols, quad_symbols, model);

                        inphase_symbols = real(combined_signal);                                        %%%                    %%% if rec ==2 receiver2                    %%%                                        if rec ==2                        if model == 0                            receive2 = s1r2_thru_channel;                        else
                            receive2 = s1r2_thru_channel + s2r2_thru_channel;                        end % model == 0                          quad2 = imag(receive2);                        inphase2 = real(receive2);                               [mf_quad_block2]=matched_filter(fs, T, pulsetype, quad2);    % Quadrature-part                        [mf_inphase_block2]=matched_filter(fs, T, pulsetype, inphase2);    % Inphase-part                                                                   %%%                        %%%     Performing synchronization                        %%%                        %[corr_maxi2,sample_nri2,sum_of_r_times_c1i2,sum_of_r_times_d2i2]=orthosynchronization(mf_inphase_block2, fs, T, training_sequence1, training_sequence2, model);					    %[corr_maxq2,sample_nrq2,sum_of_r_times_c1q2,sum_of_r_times_d2q2]=orthosynchronization(mf_quad_block2, fs, T, training_sequence1, training_sequence2, model);                                                sample_nr2 = 81;                                %if corr_maxq2 > corr_maxi2       % choose the correlation with the biggest max value                        %   sample_nr2 = sample_nrq2;                        %else                        %   sample_nr2 = sample_nri2;                        %end                                                 block_length=length(outdata_block);                        [inphase_symbols2]=down_sampler(mf_inphase_block2,sample_nr2,fs,T,block_length);                        [quad_symbols2]=down_sampler(mf_quad_block2,sample_nr2,fs,T,block_length);                                        
					    %%%                        %%% Estimating Channel 3 and 4.                        %%%                        %[alpha3est,theta3est] = channel_estimator_(quad_symbols2(1:training_length1), ...                         %   inphase_symbols2(1:training_length1), training_sequence1, model, 1);                               if model == 1
                                 %  [alpha4est,theta4est] = channel_estimator_(quad_symbols2(1:training_length1), ...                           %     inphase_symbols2(1:training_length1), training_sequence2, model, 2);                        elseif model == 0                            alpha4 = 1;                            theta4 = 0;                        end % model == 1                                 alpha3est = alpha3;                        theta3est = theta3;                        alpha4est = alpha4;                        theta4est = theta4;                                                               %a3err=alpha3est-alpha3                        %a4err=alpha4est-alpha4                        %t3err=theta3est-theta3                        %t4err=theta4est-theta4
                                                                   [combined_signal2]=combiner(alpha3est, theta3est, alpha4est, theta4est, ...
                            inphase_symbols2, quad_symbols2, model);				                        inphase_symbols2 = real(combined_signal2);                                            %%%                        %%% add the signals from receiver1 and receiver2                         %%%                                            inphase_symbols = inphase_symbols + inphase_symbols2;                    end %receiver2                            %%%
                    %%% Detector
                    %%%
                     
                    if model == 1   % if orthogonal training-sequences.                        [indata_est,training1_est,training2_est]=detector_ortho(inphase_symbols,training_length1,training_length2,model);                    else                        [indata_est,training1_est,training2_est]=detector(inphase_symbols,training_length1,training_length2,model);                    end                                        %nr_of_errors = sum(abs(data_block-indata_est));                    %nr_of_errors_in_rate = nr_of_errors/length(data_block);                                       %error(i+1)   = nr_of_errors_in_rate;                     error(i+1)  =   sum(abs(data_block-indata_est))/length(data_block);                end % block                 error_mean = mean(error);                EbN0error = [EbN0error error_mean];            end % EbN0                       compared_error = [compared_error EbN0error']        end % if !...    end %for(model)end %for(rec)