[]=Blind_recognition_STBC4_simulation_results(N,nr,Modulation,Mod_State_Nb)

function []=Blind_recognition_STBC4_simulation_results(N,nr,Modulation,Mod_State_Nb)

%HELP Blind_recognition_STBC4_simulation_results
%
%Section VI.B of the paper. Perform monte-carlo simulation to approximate
%the average_probability of correct recognition (APCR). 7 STBCs using 4 antennas are
%considered.
%
%Input:    
%
%        - N (Number of received sample)
%        - nr: Number of receiver antennas
%        - Modulation: 'PSK','PAM' or 'QAM'
%        - Mod_State_Nb: Nomber of state of the modulation
%
%Ouput:  - graphics and matlab files are stored in the
%          simulation_results directory.
%
%exemple: Blind_recognition_STBC4_simulation_results(512,8,'PSK',4)
%
%reference: [1] V. Choqueuse, M. Marazin, L. Collin, K.C. Yao, K and G. Burel
%"Blind Recognition of Linear SpaceTime Block Codes: A Likelihood-Based
%Approach" IEEE Transactions on Signal Processing, Vol 58 (3), p 1290-1299,
%2010

code_name_list={'SM','STBC4','STBC4','STBC4','STBC4','OSTBC4','OSTBC4'};
rate_list={'4','3','2','2','1','3/4','1/2'};
num_code_list=[1 1 1 2 1 1 1];
SNR_list=[-15:5:5];
Nb_Monte_Carlo_trials_per_SNR=1;
nt=4;


fprintf('----- Simulation parameters -----\n');
fprintf('Modulation: %i-%s \n',Mod_State_Nb,Modulation);
fprintf('Receiver: %i samples intercepted by %i antennas\n',N,nr);


% for loop for each SNR value

for index_SNR=1:length(SNR_list)

    SNR=SNR_list(index_SNR);
    fprintf('----- SNR=%i dB (%i sur %i) -----\n',SNR,index_SNR,length(SNR_list));
   
    % For each SNR, counts the number of correct detection
    % Nb_correct_detection_Optimal=0;
    Nb_correct_detection_SOS=0;
    Nb_correct_detection_CP=0;

    % Probability_correct_recognition_Optimal(index_SNR,:)=zeros(1,length(num_code_list));
    Probability_correct_recognition_SOS(index_SNR,:)=zeros(1,length(num_code_list));
    Probability_correct_recognition_CP(index_SNR,:)=zeros(1,length(num_code_list));


    % Monte Carlo trials
    for index_trial=1:Nb_Monte_Carlo_trials_per_SNR
      
        for index_code=1:length(num_code_list)

            code_name=code_name_list(index_code);
            rate=rate_list(index_code);
            num_code=num_code_list(index_code);
            
            %% MIMO Communication Simulation
            
            symbol_rate=str2num(cell2mat(rate));
            
            data=randint(1,N*symbol_rate,Mod_State_Nb);
            
            %modulation
            switch Modulation
                case 'PSK'
                    symbols=pskmod(data,Mod_State_Nb); 
                    alphabet_modulation=pskmod([0:Mod_State_Nb-1],Mod_State_Nb);
                case 'PAM'
                    symbols=pammod(data,Mod_State_Nb);
                    alphabet_modulation=pammod([0:Mod_State_Nb-1],Mod_State_Nb);
                case 'QAM'
                    symbols=qammod(data,Mod_State_Nb);
                    alphabet_modulation=qammod([0:Mod_State_Nb-1],Mod_State_Nb);
            end

            %normalisation $E[|s|^2]=2$.
            variance_signal=mean(abs(alphabet_modulation).^2);
            symbols=symbols/sqrt(variance_signal);    
            coded_signals=space_time_coding(symbols,cell2mat(code_name),cell2mat(rate),num_code);
            coded_signals=coded_signals(:,1:N);   %%make sure that the number of symbols transmitted is equal to N 
            
            total_power_nt=trace((coded_signals*coded_signals')/(size(coded_signals,2)));
            
            % channel matrix
            H=sqrt(0.5)*(randn(nr,nt)+i*randn(nr,nt));
            Y_without_noise=H*coded_signals;

            %Additive noise
            noise_variance=total_power_nt/(10^(SNR/10));
            B=(sqrt(noise_variance/2))*(randn(nr,N)+i*randn(nr,N));
            Y=Y_without_noise+B;
       
            %% Optimal Classifier (too computational-expensive)
            % [code_name_d,Rate_d,Num_code_d]=Optimal_Classifier(Y,H,Modulation,Mod_State_Nb,noise_variance,code_name_list,rate_list,num_code_list,0);
            %             
            % result=(strcmp(code_name,code_name_d)+strcmp(rate,Rate_d)+(num_code==Num_code_d))==3;
            % Nb_correct_detection_Optimal=Nb_correct_detection_Optimal+result;
            % Probability_correct_recognition_Optimal(index_SNR,index_code)=Probability_correct_recognition_Optimal(index_SNR,index_code)+result;
            
            %% SOS-STBC Classifier            
                        
            [code_name_d,Rate_d,Num_code_d]=SOS_STBC_Classifier(Y,H,noise_variance,code_name_list,rate_list,num_code_list,0);
            
            result=(strcmp(code_name,code_name_d)+strcmp(rate,Rate_d)+(num_code==Num_code_d))==3;
            Nb_correct_detection_SOS=Nb_correct_detection_SOS+result;     
            Probability_correct_recognition_SOS(index_SNR,index_code)=Probability_correct_recognition_SOS(index_SNR,index_code)+result;

            %% CP Classifier    
                
            [code_name_d,Rate_d,Num_code_d]=CP_Classifier(Y,code_name_list,rate_list,num_code_list,0);
            
            result=(strcmp(code_name,code_name_d)+strcmp(rate,Rate_d)+(num_code==Num_code_d))==3;
            Nb_correct_detection_CP=Nb_correct_detection_CP+result;     
            Probability_correct_recognition_CP(index_SNR,index_code)=Probability_correct_recognition_CP(index_SNR,index_code)+result;
         
        end          
    end        

    % Compute Average Probability of Correct recognition (AVPCR)

    % APCR_Optimal(index_SNR)=Nb_correct_detection_Optimal/(Nb_Monte_Carlo_trials_per_SNR*length(num_code_list));
    % fprintf('Average Probability of correct detection: Optimal %f\n',APCR_Optimal(index_SNR));
       
    APCR_SOS(index_SNR)=Nb_correct_detection_SOS/(Nb_Monte_Carlo_trials_per_SNR*length(num_code_list));
    fprintf('Average Probability of correct detection: Optimal %f\n',APCR_SOS(index_SNR));
    
    APCR_CP(index_SNR)=Nb_correct_detection_CP/(Nb_Monte_Carlo_trials_per_SNR*length(num_code_list));
    fprintf('Average Probability of correct detection: Optimal %f\n',APCR_CP(index_SNR));
    
    % Compute Probabilty of correct recognition for each code
    % Probability_correct_recognition_Optimal(index_SNR,:)=Probability_correct_recognition_Optimal(index_SNR,:)/Nb_Monte_Carlo_trials_per_SNR;
    Probability_correct_recognition_SOS(index_SNR,:)=Probability_correct_recognition_SOS(index_SNR,:)/Nb_Monte_Carlo_trials_per_SNR;
    Probability_correct_recognition_CP(index_SNR,:)=Probability_correct_recognition_CP(index_SNR,:)/Nb_Monte_Carlo_trials_per_SNR;
    
    
end    


cd simulation_results
              
%% Display Average Probability of Correct recognition (AVPCR)
h_APCR=figure;
hold on;
plot(SNR_list,APCR_SOS,'o-');
plot(SNR_list,APCR_CP,'*-');
hold off;
grid;
xlabel('SNR (dB)');
legend('SOS-STBC','(blind) CP');
ylabel('Average Probability of Correct Detection');
% save figure
instruction=sprintf('Simulation_result_STBC4_APCR_N_%i_nr_%i.fig',N,nr);
saveas(h_APCR,instruction);

%% Display Probability of Correct recognition for each code
for index_code=1:length(num_code_list)
    
    symbol_rate=cell2mat(rate_list(index_code));

    h=figure;
    hold on;
    plot(SNR_list,Probability_correct_recognition_SOS(:,index_code),'o-');
    plot(SNR_list,Probability_correct_recognition_CP(:,index_code),'*-');
    hold off;
    grid;
    xlabel('SNR (dB)');
    legend('SOS-STBC','(blind) CP');
    ylabel(sprintf('Probability of Correct Detection for %s-rate %s (%i)',symbol_rate(find(symbol_rate ~= '/')),code_name_list{index_code},num_code_list(index_code)));
    ylim([0 1.01]);
    
    % save figure
    instruction=sprintf('Simulation_result_STBC4_PCR_%s_rate_%s_%i_N_%i_nr_%i.fig',symbol_rate(find(symbol_rate ~= '/')),code_name_list{index_code},num_code_list(index_code),N,nr);
    saveas(h,instruction);
end

fprintf('OK\n\n');