mimo2x1.m

alamouti結合ofdm系統下 做的仿真系統。 提供大家討論使用 不懂的再與我聯絡吧

% Transmission - Works well. 27 Mar 03, 21:40.
% Kim Thanh Tung 
clear;

Fs = 44100; %Hz 
Fc = 11025; %Hz
L = 4; 

nr_block = 40;
block_len = 800

half_filter_len = 7;
raised_cos_pulse = root_raised_cosine(L, .22, half_filter_len);
delay = half_filter_len;

guard_bits = [1 1];

training_len = length(training_sequence);
[training_sequence_left training_sequence_right] = training_sequence;

out_count = 1;  % sampling counter (sending to antenna)
modulus = 0;   % mod 4....
count_err = 0;

left = zeros(1, (block_len/4+training_len+4)*L*nr_block); 
right = zeros(1, (block_len/4+training_len+4)*L*nr_block); 


%%% Read data file
%fid = fopen('transmit.dat', 'r');
data = zeros(block_len, nr_block);
for k = 1:nr_block
    data(:, k) = random_data(block_len)';
end

%data = zeros(block_len, nr_block);

for block_count = 1:nr_block
% for each block of data
    b = data(:, block_count);
% 16-QAM bit mapper
    [Ich, Qch] = qam16(b);

    
%%% Alamouti Encoder
    [Ich_left, Qch_left, Ich_right, Qch_right] = alamouti_encoder(Ich, Qch);
    

%%%%%%%%%%%%%%%%%%%% LEFT channel %%%%%%%%%%%%%%%%%%%%
    
% Upsample L 
% Adding training data    
    Ich_left = [guard_bits training_sequence_left Ich_left guard_bits];
    Qch_left = [guard_bits zeros(1, training_len) Qch_left guard_bits];

    Ich_left = upsample(Ich_left, L);

    temp = conv(raised_cos_pulse, Ich_left);
    Ich_left = temp(delay:delay+length(Ich_left)-1);

    Qch_left = upsample(Qch_left, L);
    temp = conv(raised_cos_pulse, Qch_left);
    Qch_left = temp(delay:delay+length(Qch_left)-1);

%%%%%%%%%%%%%%%%%%%% RIGHT channel %%%%%%%%%%%%%%%%%%%%    
    
    Ich_right = [guard_bits training_sequence_right Ich_right guard_bits];
    Qch_right = [guard_bits zeros(1, training_len) Qch_right guard_bits];

    Ich_right = upsample(Ich_right, L);

    temp = conv(raised_cos_pulse, Ich_right);
    Ich_right = temp(delay:delay+length(Ich_right)-1);

    Qch_right = upsample(Qch_right, L);
    temp = conv(raised_cos_pulse, Qch_right);
    Qch_right = temp(delay:delay+length(Qch_right)-1);
   
%%%%%%%% Shift to 11.025kHz

    for k = 1:length(Ich_left);
    % cos  = 1 0 -1 0 * I
    % -sin = 0 1 0 -1 * Q
        switch(modulus)
            case 0
                left(out_count) = +Ich_left(k);
                right(out_count) = +Ich_right(k);                
            case 1
                left(out_count) = -Qch_left(k);            
                right(out_count) = -Qch_right(k);                                
            case 2
                left(out_count) = -Ich_left(k);            
                right(out_count) = -Ich_right(k);                                
            otherwise
                left(out_count) = +Qch_left(k);
                right(out_count) = +Qch_right(k);                                
        end
        modulus = mod(modulus + 1,4);
        out_count = out_count + 1;
    end

end % for block_count

% Pilot Tones:
pilot = zeros(1,12);
pilot(1:4:end) = 6;
pilot(3:4:end) = -6;
left = [pilot left];
right = [zeros(size(pilot)) right];

%%% DSP
%left = (2^15-1)/6*left;
%right = (2^15-1)/6*right;


%%%%%%%%%%%%%%%%
% Delay - Channel
inSeq = [zeros(1, 105) left+right]; 
inSeq = inSeq + .1*rand(1, length(inSeq));

%%%%%%%%%%%%%%%%
%% Receiver
plot(inSeq);


startPos = 1;
% 
for k = 1:500
     if (inSeq(k) > 3) & (inSeq(k+4) > 3) & (inSeq(k+8) > 3)
         startPos = k+12;
         break;
     end
 end
 
% startPos = 118;

%     % The algorithm should be:
%     % - After pilot tone found, search for sync. in a certain time.
%     % - If sync. is not found, search for pilot tones again
% % cos: 1 0 -1 0  I-channel
% % -sin: 0 1 0 -1  Q-channel
len = length(inSeq);
Ich = zeros(1, len);
Qch = zeros(1, len);

modulus = 0;
count = 1;
for k = startPos:len  
     switch(modulus)
         case 0
             Ich(count) = inSeq(k);
             %Qch(k) = 0
         case 1
             %Ich(k) = 0
             Qch(count) = -inSeq(k);            
         case 2
             Ich(count) = -inSeq(k);            
             %Qch(k) = 0
         otherwise
             %Ich(k) = 0
             Qch(count) = +inSeq(k);                        
     end            
     modulus = mod(modulus+1, 4);
     count = count + 1;
end

% figure;
% plot(Ich(1:100));
 
% % Low pass filter
lpf = low_pass_filter(Fs);

Ich = conv(lpf, Ich);
Qch = conv(lpf, Qch);

% figure
% psd(Ich, 1024, Fs);
 
% % Matched filter
Ich = conv(raised_cos_pulse, Ich);        
Qch = conv(raised_cos_pulse, Qch);        

% figure;
% psd(Ich, 1024, Fs);
 
%%%%%% Synchronization &&& Detection
block_start_pos = 1;
for block_count = 1:nr_block
    sync_I_from_left = zeros(1, 60);
    sync_I_from_right = zeros(1, 60);    
    sync_Q_from_left = zeros(1, 60);
    sync_Q_from_right = zeros(1, 60);    
    
    sync_I = zeros(1, 60);    
    
    pos_I = 1;
    for k = block_start_pos:block_start_pos + 59
        sync_I_from_left(pos_I) = dot(training_sequence_left, Ich(k:L:k+(training_len-1)*L));
        sync_Q_from_left(pos_I) = dot(training_sequence_left, Qch(k:L:k+(training_len-1)*L));        
        
        sync_I_from_right(pos_I) = dot(training_sequence_right, Ich(k:L:k+(training_len-1)*L));        
        sync_Q_from_right(pos_I) = dot(training_sequence_right, Qch(k:L:k+(training_len-1)*L));                
        
        
        % May need to perform on Q as well
        pos_I = pos_I + 1;
    end
    
    sync_I  = sync_I_from_left.*sync_I_from_left + sync_I_from_right.*sync_I_from_right;
    sync_Q  = sync_Q_from_left.*sync_Q_from_left + sync_Q_from_right.*sync_Q_from_right;
    
    [Imax, idx] = max(sync_I + sync_Q);
    T_samp = block_start_pos + idx - 1;
%    [Imax_right, idx_right] = max(sync_I_from_right);    
    
%      if Imax_left > Imax_right
%          T_samp = block_start_pos + idx_left - 1;
%      else
%          T_samp = block_start_pos + idx_right - 1;        
%      end
%    
%    T_samp = T_samp + 1
      
%    [Qmax_left, idx] = max(sync_Q_from_left);       
%    [Qmax_right, idx_right] = max(sync_Q_from_right);        
    
    if block_count == 1
        figure
        plot(1:60, sync_I+sync_Q)

    end

    
    %%% Downsampling
    bits_buff_I = zeros(1, block_len/4 + training_len);
    bits_buff_Q = zeros(1, block_len/4 + training_len);
    next_pos = T_samp;
    k = 1;
    while k <= training_len + block_len/4
        bits_buff_I(k)  = Ich(next_pos);
        bits_buff_Q(k)  = Qch(next_pos);    
        next_pos = next_pos + L;
        k = k+1;
    end

    %% update block_start_pos
    block_start_pos  = next_pos + 1; %% skip the two guard bits

  %% Channel Estimation
    h_I_from_left = dot(training_sequence_left, bits_buff_I(1:training_len))/training_len/9;
    h_Q_from_left = dot(training_sequence_left, bits_buff_Q(1:training_len))/training_len/9;

    h_I_from_right = dot(training_sequence_right, bits_buff_I(1:training_len))/training_len/9;
    h_Q_from_right = dot(training_sequence_right, bits_buff_Q(1:training_len))/training_len/9;
    
    % % Combination
    [shifted_bits_buff_I, shifted_bits_buff_Q] = alamouti_decoder21(bits_buff_I(training_len+1:end), bits_buff_Q(training_len+1:end), h_I_from_left, h_Q_from_left, h_I_from_right, h_Q_from_right);

    chGain = h_I_from_left^2 + h_Q_from_left^2 + h_I_from_right^2 + h_Q_from_right^2;
    % % Detector
% % Detect should take the channel gain as a parameters, otherwise we wont
% % get the best threshold
    bhat = detect(shifted_bits_buff_I(1:end), shifted_bits_buff_Q(1:end), chGain);
% 
% % Count errors
%   err_count = 0;
    for k = 1:block_len
        if (data(k, block_count) ~= bhat(k)) 
            count_err = count_err + 1;
  %         k
        end
    end
    %if block_count == 1
    %    count_err
    %end

end %% for block_count
 
% 
% 
count_err/nr_block/block_len