#29: first version of B-EM and SB-EM

Algorithms/Algo_Mode/Blind/EM/B_EM_Nonlinear_MIMO.m

@@ -0,0 +1,176 @@
+function Err = B_EM_Nonlinear_MIMO(Op, SNR_i, Output_type)
+
+%% Blind Expectation-Maximization for Non-linear MIMO communications
+%
+%% Input:
+    % + 1. Ns: number of sample data
+    % + 2. Nt: number of transmit antennas
+    % + 3. Nr: number of receive antennas
+    % + 4. M: length of the channel
+    % + 5. Ch_type: type of the channel (real, complex,
+    % user' input
+    % + 6. Mod_type: type of modulation (All)
+    % + 7. Threshold: threshold of EM algorithm
+    % + 8. Np: number of pilot symbols
+    % + 9. N_iter_max: number of maximum iterations EM 
+    % + 10. SNR_i: signal noise ratio
+    % + 11. Output_type: SER / BER / MSE Channel
+%
+%% Output:
+    % + 1. Err: SER / BER / MSE Channel
+%
+%% Algorithm:
+    % Step 1: Initialize variables
+    % Step 2: Return 
+%
+% Ref: Ouahbi Rekik, Karim Abed-Meraim, Mohamed Nait-Meziane,
+% Anissa Mokraoui, and Nguyen Linh Trung, "Maximum likelihood 
+% based identification for nonlinear multichannel communications 
+% systems," Signal Processing, vol. 189, p. 108297, Dec. 2021.
+%
+%% Require R2006A
+
+% Author: Ouahbi Rekik, L2TI, UR 3043, Université Sorbonne Paris Nord, France
+
+% Adapted for InSI by Do Hai Son, 1-Aug-2023
+% InSI: A MatLab Toolbox for Informed System Identification in 
+% Wireless Communications
+% https://avitech-vnu.github.io/InSI
+% Project: NAFOSTED 01/2019/TN on Informed System Identification
+% PI: Nguyen Linh Trung, Vietnam National University, Hanoi, Vietnam
+% Co-PI: Karim Abed-Meraim, Université d’Orléans, France
+
+
+% Initialize variables
+Ns        = Op{1};     % number of sig sequences
+Nt        = Op{2};     % number of transmit antennas
+Nr        = Op{3};     % number of receive antennas
+M         = Op{4};     % length of the channel
+Ch_type   = Op{5};     % complex
+Mod_type  = Op{6};     
+threshold = Op{7};     % threshold of EM algorithm
+Np        = Op{8};     % number of pilot symbols 
+N_iter_max= Op{9};     % number of maximum iterations EM 
+Np_init   = Np;
+
+
+% Generate input signal
+modulation = {'Bin', 'QPSK', 'QAM4', 'QAM16', 'QAM64', 'QAM128', 'QAM256'};
+
+% Generate channel
+if Ch_type == 2
+    channels_mat = Generate_channel(Nr, Nt * 2*(M+1) - 1, Ch_type) * sqrt(2);
+else
+    channels_mat = Generate_channel(Nr, Nt * 2*(M+1) - 1, Ch_type);
+end
+channels_vec = reshape(channels_mat, [], 1); % vectorize channel matrix
+
+% Generate signals
+[sig_src, data] = eval(strcat(modulation{Mod_type}, '(1000)'));
+alphabet = unique(sig_src);
+P = length(alphabet);
+average_energy_per_symb = norm(alphabet)^2/P; % or: 2/3*(P-1);
+alphabet = alphabet/sqrt(average_energy_per_symb); % normalization
+
+%% transitions, predecessors, successors
+trans1 = SB_EM_cal_trans(alphabet, Nt*(M+1));
+trans2 = trans1.^2;                % quadratic model
+trans  = [trans1; trans2];
+nbr_etats    = P^(Nt*M); %number of states
+successors   = zeros(nbr_etats,P^Nt);
+predecessors = zeros(nbr_etats,P^Nt);
+tableau = zeros(P^(Nt*(M+1)),2); % tableau(:,1)-->predecesseurs (i) ,tableau(:,2)-->successeurs (j)
+
+for i=1:nbr_etats
+    predecessors(i,:) = SB_EM_PredecessorsMIMO( i ,P, M, Nt );
+    [itransPred,~]    = SB_EM_FindTransitionsPredMIMO(predecessors(i,:),i,M, P, Nt );
+    tableau(itransPred,2) = i; 
+
+    successors(i,:)   = SB_EM_SuccessorsMIMO( i ,P, M, Nt );
+    [itransSuccess,~] = SB_EM_FindTransitionsSuccessMIMO(i,successors(i,:),M, P, Nt );
+    tableau(itransSuccess,1) = i;
+end
+
+%% Transmitted signal
+S = [];
+D = [];
+
+for i = 1:Nt
+    [sig_src, data] = eval(strcat(modulation{Mod_type}, '(Ns + M)'));
+    S = [S; sig_src.'];
+    D = [D; data.'];
+end
+S_decim = D;
+S       = S/sqrt(average_energy_per_symb);
+
+Stop1 = toeplitz(S(1,M+1:-1:1),S(1,M+1:Ns+M));
+D_Stop1 = toeplitz(S_decim(1,M+1:-1:1),S_decim(1,M+1:Ns+M));
+
+if Nt==1
+    S_mat = [Stop1; Stop1.^2];
+
+    D_mat = [D_Stop1; D_Stop1.^2];
+else    
+    Stop2 = toeplitz(S(2,M+1:-1:1),S(2,M+1:Ns+M));
+    S_mat = [Stop1; Stop2; Stop1.^2; Stop2.^2];
+
+    D_Stop2 = toeplitz(S_decim(2,M+1:-1:1), S_decim(2,M+1:Ns+M));
+    D_mat = [D_Stop1; D_Stop2; D_Stop1.^2; D_Stop2.^2];
+end
+
+%% Received signal
+sig_cap_sans_bruit = channels_mat*S_mat;
+
+% ... pilots
+Sp = [];
+Dp = [];
+
+for i = 1:Nt
+    [sig_src, data] = eval(strcat(modulation{Mod_type}, '(Ns + M)'));
+    Sp = [Sp; sig_src.'];
+    Dp = [Dp; data.'];
+end
+Sp_decim = Dp;
+Sp       = Sp/sqrt(average_energy_per_symb);
+
+Stop1p = toeplitz(Sp(1,M+1:-1:1),Sp(1,M+1:Np+M));
+if Nt==1
+    Sp_mat = [Stop1p; Stop1p.^2];   
+else 
+    Stop2p = toeplitz(Sp(2,M+1:-1:1), Sp(2,M+1:Np+M));
+    Sp_mat = [Stop1p; Stop2p; Stop1p.^2; Stop2p.^2];
+end
+sigp_cap_sans_bruit = channels_mat*Sp_mat;
+
+% --- received signal strength ----------------------------------
+Pr = norm(sig_cap_sans_bruit,'fro').^2/numel(sig_cap_sans_bruit);
+
+sigmav2 = 10^((10*log10(Pr)-SNR_i)/10); % noise variance for received signal
+
+noise = sqrt(sigmav2)*(randn(Nr,Ns) + randn(Nr,Ns)*1i)/sqrt(2);
+
+sig_cap = sig_cap_sans_bruit + noise;
+sigcap_pilot = sigp_cap_sans_bruit + noise(:,1:Np); %noisep;
+
+%----------------- pilot-based initialization -------------------
+canauxp  = sig_cap(:,1:Np_init)*pinv(S_mat(:,1:Np_init));
+
+%-------------------------- B EM -------------------------------
+[probaB, ~, channels_mat_EM_B, ~] = ...
+    B_EM_func(sig_cap.', trans, canauxp, sigmav2, predecessors, ...
+    successors, tableau, N_iter_max, threshold, Nt);
+
+%% Vectorize estimated channels
+channels_vec_EM_B = reshape(channels_mat_EM_B, [],1); % h_{EM-SB} semi-blind
+
+%% Data detection
+est_src = SB_EM_min_symbol_errorMIMO(probaB, alphabet, M, Nt);
+
+% Compute SER / BER / MSE channel
+if Output_type ~= 4
+    Err = ER_func(D_mat(1,:).', est_src, Mod_type, Output_type);
+else
+    Err = ER_func(channels_vec.', channels_vec_EM_B, Mod_type, Output_type);
+end
+
+end

Algorithms/Algo_Mode/Blind/EM/Params/B_EM_Nonlinear_MIMO_params.m

@@ -0,0 +1,55 @@
+classdef B_EM_Nonlinear_MIMO_params
+    %Params Summary of this class goes here
+    %   Detailed explanation goes here
+
+    properties
+        % Parameters
+        num_params = 9
+        params = {'No. samples', 'No. transmitters', 'No. receivers', 'Channel order', 'Channel type', 'Modulation', 'Threshold', 'No. pilots', 'No. iter max'}
+        notations = {'Ns', 'Nt', 'Nr', 'M', 'ChType', 'Mod', 'threshold', 'Np', 'N_iter_max'}
+        tooltips = {}
+        % Type of the UIControl: edit_text   = 1
+        %                        popup_menu  = 2
+        %                        button      = 3
+        params_type = [1, 1, 1, 1, 2, 2, 1, 1, 1]
+        values = {100, 2, 4, 2, {'Real', 'Complex', 'Input'}, {'Binary', 'QPSK', '4-QAM', '16-QAM', '64-QAM', '128-QAM', '256-QAM'}, 0.001, 10, 20}
+        default_values = {100, 2, 4, 2, 2, 3, 0.001, 10, 20}
+
+        % Default SNR and Monte
+        default_Monte = 10
+        default_SNR = '-10:5:20'
+
+        % Output
+        % Type of the outputs: SER Sig = 1
+        %                      BER Sig = 2
+        %                      MSE Sig = 3
+        %                      MSE Cha = 4
+        outputs = [1, 2, 4]
+        default_output = 1
+
+        % Figure
+        sys_model = 'Default.png'
+        title     = {'B-EM Non-linear MIMO'}
+        xlabel    = {'SNR (dB)', 'SNR (dB)', 'SNR (dB)', 'SNR (dB)'}
+        ylabel    = {'SER', 'BER', 'MSE Signal (dB)', 'MSE Channel (dB)'}
+        trigger   = false 
+        position
+        linewidth = 1
+        color     = 'k'
+
+        % Triggers/Flags
+        has_inter     = [true, true, true, false, false, false, false, false, false]
+        rect = {}
+        rect_position = {[150 480 100 100], [260 1130 290 100], [1470 1130 290 100], 0, 0, 0, 0, 0, 0}
+        rect_linewidth = {2, 2, 2, 2, 2, 2, 2, 2, 2}
+        rect_color     = {'b', 'g', 'r', 'b', 'r', 'g', 'g', 'b', 'b'}      
+
+        % Reference website
+        web_url = 'https://www.sciencedirect.com/science/article/abs/pii/S0165168421003340'
+    end
+
+    methods (Access = private)
+
+    end
+
+end

Algorithms/Algo_Mode/Blind/EM/Params/B_EM_func.m

@@ -0,0 +1,59 @@
+% Calcul des parametres (H,sigma) par l'algorithme EM maximisant la
+% vraisemblance des signaux observes.
+%
+% SYNTAXE [proba,estsigma1,estcan1]=EM(sig_cap,trans,estcan,estsigma,predess,success,nb_iter,seuil)
+%
+% sig_cap : tableau des signaux recus par les capteurs
+% trans : tableau ((M+1) X P^(M+1)) des transitions
+% estcan : estimee initiale des coefficients des canaux
+% estsigma : estimmee initiale de la puissance du bruit
+% predess : tableau des predecesseurs d'etats
+% success : tableau des successeurs d'etats
+% nb_iter : nbr maximal d'iterations autorise
+% seuil : seuil d'arret des iterations.
+%
+% sortie:
+% proba : tableau des densites de probabilites L(Y,Xt=En) a la fin des iterations.
+% estsigma1 : reestimee de sigma a la fin des iterations
+% estcan1 : reestimee de can a la fin des iterations
+%
+function [proba,estsigma1,estcan1,compt] = B_EM_func(sig_cap,trans,estcan,estsigma,predess,success,tableau,nb_iter,seuil,Nt)
+
+old_vrais = -1;
+test      = seuil+1;
+compt     = 0;
+M         = length(estcan(1,:))/(2*Nt)-1;
+Ncap      = size(estcan,1);
+P         = (length(predess(1,:)))^(1/Nt);
+T         = length(sig_cap(:,1));
+
+while(test > seuil) && (compt<nb_iter)
+    compt=compt+1;
+    %%%%%%% Algorithme
+
+    %tic
+    tab = SB_EM_cal_tab(sig_cap,trans,estcan,estsigma);
+
+    [alpha,beta,Rho] = SB_EM_alpha_beta2(tableau,tab,M,T,Nt,P);
+
+    proba = SB_EM_cal_prob2(alpha,beta,tab,tableau,P,M,T,Nt);
+    %t=toc
+
+    new_vrais = sum(proba(:,1))/prod(Rho);
+    if (old_vrais>new_vrais)
+        error('The likelihood is diminishing!!!')
+    end
+
+    old_vrais = new_vrais;
+    [estsigma1,estcan1] = SB_EM_cal_sigmaH_B(proba,sig_cap,trans,Nt);
+
+    vec_new  = reshape(estcan1.',Ncap*Nt*2*(M+1),1);
+    vec_old  = reshape(estcan.',Ncap*Nt*2*(M+1),1);
+
+    test = norm(vec_new-vec_old)/norm(vec_old); %+abs(estsigma1-estsigma);
+
+    estcan   = estcan1;
+    estsigma1 = estsigma;
+
+end
+end

Algorithms/Algo_Mode/Blind/SS/Params/B_SS_params.m

@@ -5,7 +5,7 @@
     properties
         % Parameters
         num_params = 6
-        params = {'No. bits', 'No. channels', 'Channel order', 'Channel type', 'Modulation', 'Window length'}
+        params = {'No. samples', 'No. channels', 'Channel order', 'Channel type', 'Modulation', 'Window length'}
         notations = {'N', 'Nr', 'ChL', 'ChType', 'Mod', 'L'}
         tooltips = {}
         % Type of the UIControl: edit_text   = 1

Algorithms/Algo_Mode/Non-blind/ZF/Params/NB_ZF_params.m

@@ -5,7 +5,7 @@
     properties
         % Parameters
         num_params = 6
-        params = {'No. bits', 'No. transmitters', 'No. receivers', 'Channel order', 'Channel type', 'Modulation'}
+        params = {'No. samples', 'No. transmitters', 'No. receivers', 'Channel order', 'Channel type', 'Modulation'}
         notations = {'N', 'Nt', 'Nr', 'ChL', 'ChType', 'Mod'}
         tooltips = {}
         % Type of the UIControl: edit_text   = 1

Algorithms/Algo_Mode/Semi-blind/EM/Params/SB_EM_Dec2Alphabet.m

@@ -0,0 +1,15 @@
+function result = SB_EM_Dec2Alphabet( x, alphabet, M )
+% this function converts the decimal "x" into its equivalent representation
+% with elements of alphabet
+
+P   = length(alphabet);
+mat = (dec2base(x-1,P,M))';
+
+result = zeros(M,1);
+
+ for i = 1:M
+     ind = mat(i,1);
+     result(i,1) = alphabet(1,str2double(ind)+1); 
+  end
+end
+

Algorithms/Algo_Mode/Semi-blind/EM/Params/SB_EM_EstimeCov.m

@@ -0,0 +1,32 @@
+%
+% script EstimeCov.m: estimation of the covariance matrix
+% spatio-temporal 0 temporal lag
+% Syntaxe R = EstimeCov(sig_cap,T,L,N)
+% Entrees
+%         sig_cap = number of signal samples
+%         T = number of signal samples
+%         L = no of receivers
+%         N = size snapshots of each sensor ( correlation window
+% Exits
+%         R = spatio-temporal covariance matrix (dimension LN x LN)
+%
+function R = SB_EM_EstimeCov(sig_cap,T,L,N)
+
+%
+%... verification of entries
+%
+if (size(sig_cap) ~= [T L])
+   error('Size of signals is incompatible with parameters');
+end
+%
+%... Calculation of covariances
+%
+R = zeros(L*N);
+for t=1:T-N+1
+ x=sig_cap(t+N-1:-1:t,:);
+ y=x(:);
+ R=R+y*y';
+end
+R=(1/(T-N+1))*R;
+
+end

Algorithms/Algo_Mode/Semi-blind/EM/Params/SB_EM_FindTransitionsPred.m

@@ -0,0 +1,30 @@
+function [Trans,y] = SB_EM_FindTransitionsPred(s1,s2,M, P )
+
+% find the transitions vector : x(t)=s1 and x(t+1) = s2
+% s1(t) ---> s1(t+1)=s2  s1 contains the predeccessors of s2
+
+l=length(s1);
+x2 =(dec2base(s2-1,P,M))';
+for i=1:l
+x1(:,i) =(dec2base(s1(i)-1,P,M))';
+end
+
+x=[];
+for i=1:l
+x=[x [x2(1);x1(:,i)]];
+end
+%x = [x2(1);x1];
+for i=1:l
+for j=1:M+1
+   y(j,i)= str2double(x(j,i));
+end
+end
+y=y';
+Trans=zeros(1,l);
+for i=1:l
+ B = sum(y(i,:).*(10.^(numel(y(i,:))-1:-1:0)));   
+ xx = num2str(B);
+ Trans(1,i) = base2dec(xx,P)+1;  
+end
+
+end