basesim_mult_blocks.m

这个程序主要是实现空时编码的matlab编解码仿真。

%
%   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 data
block_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-frequency
T               = 0.0005;                  % Symboltime

pulsetype       = 4;							% 1: rect 2:root-rais-cos 3:hamming 4:raised-cosine
model		    = 1;							% 0: no Alamouti (1*1) 1: Alamouti (2*1)
rec             = 2;                            % 2: two receivers else one receiver     
Eb              = 1;                            % Signalenergi
sigma2          = 1;
gamma0dB        = 10;

max_delay       = 880;                          % max length of delay vector
window_length   = 100;                          % must be equal to block_size + training_length1
nr_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);
end

for(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