rx.asv

mimo2x2天线选择系统的全系统matlab程序,先前的是dsp程序.

 function BER = rx(detect_type,select_type)
% Transmission
% Maxime Maury
% 05-04-25

% If argument is omitted, the channel is simulated
if (nargin < 1)
    detect_type =3;
end
if (nargin < 2)
    select_type = 1;
end

close all;
disp('-- Start RX --')

% ---------------------------------------------------
% Global parameters
% ---------------------------------------------------

disp('Loading...');

% Real or simulated channel?
load('SimType');

% Load global parameters
load('tx_param');

% Time bewteen 2 switches within a frame (before upsampling)
switch_len = guard_len + training_s_len + training_len;

% Delta Symbols set for exhaustive searching
[Delta_S_Set, Set_len] = D_Symbol_Set;

% Initial Antenna selection from {1,2,3,4} but only 4 different choices:
% that are: 1&3, 1&4, 2&3, 2&4
A = 1;
B = 3;

% Type of different detectors, parameters for Detector.m 
ML = 1;
JMMSE = 2;
ZF = 3;

% Type of different antenna selection criteria methods 
MBER = 1;
MMI = 2;
MNP = 3; %Minimum Noise Power
MMNP = 4;

% (De) activate antenna selection
antenna_select = 1;

% H channel matrices Initial
H = zeros(4,nr_frames*2+2);
H([A B],1:2) = eye(2); %Initial
H([3-A 7-B],1:2) = eye(2); 

% H channel matrices Initial (simulated channel)
if (real_ch==0)
    H_s = zeros(2,frames_len*2);
end

% Memory for the channel estimation
lambda_H = 0.95; 

% Load the channel applied 
if (real_ch==0)
    load('ChannelMatrix')
end

% Read data file (binary data sent via the TX)
fid = fopen('transmit.dat', 'r');
total_data_len = data_len_bits*2*nr_frames;
data_sent = fread(fid, total_data_len, 'int16');
data_sent = data_sent';
fclose(fid);

% Read output of the channel
load('ChannelOutput');

% data_rec is the data received in bits
data_rec = zeros(1,data_len_bits*2*nr_frames);

% data_rec_split is the data received in bits on Channel 1 and on Channel 2
data_rec_split = zeros(2,data_len_bits*nr_frames);

% M_sync is the analysis window length for synchronization
M_sync = 250;
    
% Out_match is the output of the match filter
Out_match = zeros(4,Lf+M_sync);

% Out_data data read in one frame
Out_data = zeros(4,data_len); % Data in symbols

% For plotting the constellation
r1_data_I_block = [];
r1_data_Q_block = [];

% ---------------------------------------------------
% Noise estimation
% ---------------------------------------------------

disp('Nois estimation...');

% Estimate the noise variance from the first values
sigma_sqr_noise = (var(Y_R(A,1:init_len)) + var(Y_R(B,1:init_len)) )/2;
sigma_sqr_noise_bb = 4 * sigma_sqr_noise;
N0 = sigma_sqr_noise_bb/2/L;
Var_symbol =  (var(Y_R(A,2000:2000+init_len)) + var(Y_R(B,2000:2000+init_len)) )/2;
Var_symbol = Var_symbol - sigma_sqr_noise;
Es_bb = Var_symbol * norm(pulse_shape)^2/L;
SNRbb = 10 * log10(Es_bb/sigma_sqr_noise_bb);
SNRVar = 10 * log10(Var_symbol/sigma_sqr_noise);

disp(['Estimated bb noise variance: ', num2str(sigma_sqr_noise_bb)]);
disp(['Var(Signal)/Var(Noise):      ', num2str(SNRVar)]);

% For plots
Nc_sync = ceil(sqrt(nr_sync_frames));
Nr_sync = ceil(nr_sync_frames/Nc_sync);
Nc_est = ceil(sqrt(nr_frames*2));
Nr_est = ceil(nr_frames*2/Nc_est);
Nc = 3;
Nr = 2;

% Length of the signal
signal_len = size(Y_R,2);

% ---------------------------------------------------
% Start of the signal detector
% ---------------------------------------------------

disp('Detection of the signal...');

% Plot the beginning of the signal
figure(1);
subplot(Nr,Nc,1)
plot(Y_R(A,1:init_len + time_end + 500));

% Determine the beginning of the signal
sig_start = 0;

% Magnitude of the signal at the beginning
signal_mag = sigma_sqr_noise_bb;

% Magnitude threshold of signal detected
mag_T = 1.65*sqrt(Var_symbol); 
if SNRVar > 10
    mag_T = 1.45*sqrt(Var_symbol);
end
% Pourcent of the last value kept
lambda = .95; 

values = [];
while ( signal_mag < mag_T )

    sig_start = sig_start + 1;
    
    % If we are at the end the detection has failed
    if (sig_start>signal_len)
        error('Signal detection has failed')
        break
    else
        signal_mag = lambda*signal_mag + (1-lambda) * abs(Y_R(A,sig_start));
        values = [values signal_mag];
    end
end

disp(['The signal starts at instant:',num2str(init_len+1)]);
disp(['Estimation:                  ',num2str(sig_start)]);

% Excepted time instant when the 1st frame starts
t_frame1 = sig_start + time_end + 1;

% The signal detection algorithm will not detect the very begining of
% the sinusoid. Instead, a delay will be genereated. Then, we need to
% compensate for this delay and stop the DPPL in advance:
estimated_delay = 100;
DPLL_end = time_end - estimated_delay;

% Plot the estimated times
plot(sig_start,0,'pg');
plot(DPLL_end + sig_start, 0 ,'sy');
plot(t_frame1,0,'*k');
legend('Signal','Detected start','End of DPLL','Frame')

% Plot the signal detection output
subplot(Nr,Nc,2);
plot(values)
hold on;
plot(mag_T*ones(1,length(values)),'r');
title('Signal detection');
grid on;

% ---------------------------------------------------
% DPLL
% ---------------------------------------------------

subplot(Nr,Nc,3);
hold on;
[deltaA, thetaA] = lock_phase(2*pi*(Fc/Fs), 0, 0, Y_R(A,sig_start:DPLL_end+sig_start));
[deltaB, thetaB] = lock_phase(2*pi*(Fc/Fs), 0, 0, Y_R(B,sig_start:DPLL_end+sig_start));
title('Frequency offset estimation: Output of the DPLL');
grid on;
drawnow

% %test the zero-crossing method
% Fc_1 = frequency_estimation(Y_R(A,sig_start:DPLL_end+sig_start),Fs);
% Fc_2 = frequency_estimation(Y_R(B,sig_start:DPLL_end+sig_start),Fs);
% Fc_zero = (Fc_1 + Fc_2)/2

freq_offset_hat = (deltaA + deltaB)/(2*2*pi);
disp(['Estimated frequency:         ', num2str(Fc+freq_offset_hat*Fs)]);
Fc_hat = Fc + freq_offset_hat*Fs;

for fr=1:nr_frames
    
    n_fr_sync = floor((fr-2)/refresh_fr) + 1;
    n_fr_reg =  (fr-1) - n_fr_sync;
    
    disp(' ');
    disp(['Frame ',num2str(fr),' Using antennas ', num2str(A),num2str(B)]);
    
    % Take a bunch of values selected by a window
    % All the following times are used for windowing. They are not expected
    % to be accurate
    t_start =   t_frame1 + n_fr_reg * Lf + n_fr_sync * (Lf+training_s_len*L) ; % Position in the frame  (in symbols)
    t_switch =  t_start + switch_len*L + M_sync;
    
    % Switch to the complementary set of the best set of antenna
    % Worst set selection (I=3-A, J=7-B)
    
    % Demodulation: Downconversion and Pulse-shaping
    
    if (antenna_select)
        I = 3 - A;
        J = 7 - B;
    else
        I = A;
        J = B;        
    end
    
    % ---------------------------------------------------
    % Demodulation and pulse-shaping
    % ---------------------------------------------------

    [rI_I, rI_Q, rJ_I, rJ_Q] = demodulate(Y_R(:,t_start:t_switch),I,J,Fc_hat, Fs, t_start, pulse_shape);        

    % ---------------------------------------------------
    % Synchronize
    % ---------------------------------------------------
    
    if (mod((fr-1),refresh_fr) == 0)
        
        % First synchronization 
%         figure(2);
%         subplot(Nr_sync,Nc_sync, floor((fr-1)/refresh_fr) + 1)
        t_first_train = synchronize(    rI_I, rI_Q, ...
            rJ_I, rJ_Q, ...
            tr_sync1_I, tr_sync1_Q , ...
            M_sync, L , 0  );
        
%         title(['Synchronization - Frame ', num2str(fr)]);
%         drawnow
        
        % Where the first sychronization training sequence stops 
        t_first_train_e  = t_first_train + training_s_len * L - 1;     
    else
        t_first_train_e  = t_first_train - 1;  
    end
          
    
    % Switch at t_first_train_e to the best set! (A, B)
    
    % ---------------------------------------------------
    % Channel estimation
    % ---------------------------------------------------
 
    % Start of the channel estimation training sequence