mimo_ofdm_stbc.m

实现简单的MIMO-OFDM系统性能仿真

clear all;
close all;

%--------------- system parameters----------------
nt = 2;
nr = 2;
fft_len = 64; % FFt size
cp_len = 16; % cp length
L = 5; % channel order
framelen = 50; % frame size
F = 1/sqrt(fft_len)*exp(-j*2*pi/fft_len*[0:fft_len-1]'*[0:fft_len-1]); % FFT matrix

%--------------- simulation parameters-----------
SNR_loop = [0:2:30];
trial_loop = [50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50];

%--------------- Monte Carlo --------------------
E = [];
for i = 1:length(SNR_loop)
SNR = SNR_loop(i);
NPW = 10^(-SNR/10);
trial = trial_loop(i);
error = 0;    
    for i = 1:trial %--------for loop for a specific SNR value
        %---------------- geneate source symbol----------------------
        S_Usr1 = 1/sqrt(2)*( sign(randn(fft_len, framelen)) + j*sign(randn(fft_len, framelen)) ); % frequency-domain symbols of TX 1
        %S_Usr2 = 1/sqrt(2)*( sign(randn(fft_len, framelen)) + j*sign(randn(fft_len, framelen)) ); % frequency-domain symbols of TX 2
        
        S1(:,1:2:framelen) = S_Usr1(:,1:2:framelen);  
        S1(:,2:2:framelen) =-conj(S_Usr1(:,2:2:framelen)); 
        S2(:,1:2:framelen) = S_Usr1(:,2:2:framelen); 
        S2(:,2:2:framelen) = conj(S_Usr1(:,1:2:framelen)); 
        
        S_t1 = F'*S1; % time-domain symbols of TX 1
        S_t2 = F'*S2; % time-domain symbols of TX 2
        S_t_cp1 = [ S_t1( end-cp_len+1 : end ,:); S_t1 ]; % attach cp
        S_t_cp2 = [ S_t2( end-cp_len+1 : end ,:); S_t2 ]; 
        s_tx1 = sqrt( fft_len/(fft_len+cp_len) )*reshape( S_t_cp1, 1, framelen*(fft_len + cp_len) );   % parallel to serial
        s_tx2 = sqrt( fft_len/(fft_len+cp_len) )*reshape( S_t_cp2, 1, framelen*(fft_len + cp_len) );    
        
        %---------------- generate channel tap------------------------
        h11 = 1/sqrt(2*nt*(L+1))*crandn(1,L+1);
        h12 = 1/sqrt(2*nt*(L+1))*crandn(1,L+1);
        h21 = 1/sqrt(2*nt*(L+1))*crandn(1,L+1);
        h22 = 1/sqrt(2*nt*(L+1))*crandn(1,L+1);
        
        %------------ data transmit through multipath channel ---------
        s_rx_1 = 0;
        s_rx_2 = 0;
        for l = 1:L+1
            s_rx_1 = s_rx_1 + h11(l)*[zeros(1,l-1) s_tx1 zeros(1,L-l+1)] + h12(l)*[zeros(1,l-1) s_tx2 zeros(1,L-l+1)];
            s_rx_2 = s_rx_2 + h21(l)*[zeros(1,l-1) s_tx1 zeros(1,L-l+1)] + h22(l)*[zeros(1,l-1) s_tx2 zeros(1,L-l+1)];
        end
        
        %------------- add channel noise ----------------------------------
        n1 = (NPW/2)*crandn( 1, length(s_rx_1) );
        n2 = (NPW/2)*crandn( 1, length(s_rx_2) );
        s_rx_1 = s_rx_1 + n1;
        s_rx_2 = s_rx_2 + n2;
        
        %--------------- receiver processing --------------------------
        S_r1 = [];
        S_r2 = [];
        for k = 1:framelen
            S_r1 = [ S_r1 s_rx_1( (fft_len + cp_len)*(k-1)+1:(fft_len + cp_len)*k ).'  ]; % sreial-to-parallel
            S_r_cp1 = S_r1( cp_len + 1:end,: ); % cp removal
            S_r_f1 = F*S_r_cp1; % receivd frequency-domain OFDM symbol
            
            S_r2 = [ S_r2 s_rx_2( (fft_len + cp_len)*(k-1)+1:(fft_len + cp_len)*k ).'  ]; % sreial-to-parallel
            S_r_cp2 = S_r2( cp_len + 1:end,: );
            S_r_f2 = F*S_r_cp2;
        end
        
        %-------------- channel frequency impulse response-------------
        h_f11 = sqrt(fft_len)*F*[h11 zeros(1,fft_len-(L+1))].';
        h_f12 = sqrt(fft_len)*F*[h12 zeros(1,fft_len-(L+1))].';
        h_f21 = sqrt(fft_len)*F*[h21 zeros(1,fft_len-(L+1))].'; 
        h_f22 = sqrt(fft_len)*F*[h22 zeros(1,fft_len-(L+1))].';
         
        %-------------- Per-tone MMSE VBLAST -------------------------------
        S1_est = [];
        %S2_est = [];
        for p = 1:fft_len
            H = [h_f11(p) h_f12(p); h_f21(p) h_f22(p)]; % MIMO channel freuqncy channel
            r_p = [S_r_f1(p, :); S_r_f2(p, :)];
            
            S_est1(p,1:2:framelen) = S_r_f1(p,1:2:framelen).*conj(h_f11(p))+conj(S_r_f1(p,2:2:framelen)).*h_f12(p)+S_r_f2(p,1:2:framelen).*conj(h_f21(p))+conj(S_r_f2(p,2:2:framelen)).*h_f22(p);
            S_est1(p,2:2:framelen) = S_r_f1(p,1:2:framelen).*conj(h_f12(p))-conj(S_r_f1(p,2:2:framelen)).*h_f11(p)+S_r_f2(p,1:2:framelen).*conj(h_f22(p))-conj(S_r_f2(p,2:2:framelen)).*h_f21(p);
    
                    
            
            
            S1_est = [S1_est; S_est1(p, : )]; % estimated source symbols of TX 1
            %S2_est = [S2_est; S_est(2, :)]; % estimated source symbols of TX 2
        end
        
        %-------------- error counting --------------------------------
        bit_true = reshape(S_Usr1 , 1, framelen*fft_len);
        real_true = sign(real(bit_true));
        imag_true = sign(imag(bit_true));
        
        bit_est = reshape( S1_est , 1, framelen*fft_len);
        real_est = sign(real(bit_est));
        imag_est = sign(imag(bit_est));
        
        error = error + sum( abs(real_true - real_est)/2 ) + sum( abs(imag_true - imag_est)/2 );
    end
    E = [E error/(2*trial*fft_len*framelen)]
    
end

% plot(SNR_loop, E)