[log_LF]=compute_LF(received_signal,H,Modulation,Mod_State_Nb,noise_variance,code_name,rate,num_code)

function [log_LF]=compute_LF(received_signal,H,Modulation,Mod_State_Nb,noise_variance,code_name,rate,num_code)

%HELP compute_LF
%Compute the exact LF for a STBC C and a given modulation. 
%
%Input:     - received_signal
%           - H (Matrice de canal)
%           - Modulation ('PSK','PAM','QAM')
%           - Mod_State_Nb: 2^(n)
%           - noise-variance
%           - code_name ('SM','Alamouti','OSTBC3')
%           - rate ('1/2','3/4','1')
%           - num_code (1 ou 2)
%
%Output:    - LF: Likelihood Function
%
%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


%% Load the code parameters (nt,l,n) and the coding matrices A
N=size(received_signal,2);
nr=size(H,1);
[nt,l]=size(space_time_coding(0,code_name,rate,num_code,1));
n=l*str2num(rate);

%load the matrice A
A_real=zeros(l,nt,n);
A_imag=zeros(l,nt,n);
A=zeros(l,nt,2*n);
for indice=1:n
    real_part=zeros(1,n);
    real_part(indice)=1;
    imag_part=zeros(1,n);
    imag_part(indice)=i; 
    [block_code_real]=space_time_coding(real_part,code_name,rate,num_code);  
    [block_code_imag]=-i*space_time_coding(imag_part,code_name,rate,num_code);  
    for matrice_index=1:l
        A_real(matrice_index,:,indice)=block_code_real(:,matrice_index);
        A_imag(matrice_index,:,indice)=block_code_imag(:,matrice_index);
        A(matrice_index,:,[indice indice+n])=[block_code_real(:,matrice_index) i*block_code_imag(:,matrice_index)];
    end
end  

%% Construct all the possible combinations S of n symbols of the given
%% modulation

% Construct the Modulation alphabet
switch Modulation
    case 'PSK'
        alphabet_modulation=pskmod([0:Mod_State_Nb-1],Mod_State_Nb);
    case 'PAM'
        alphabet_modulation=pammod([0:Mod_State_Nb-1],Mod_State_Nb);
    case 'QAM'
        alphabet_modulation=qammod([0:Mod_State_Nb-1],Mod_State_Nb);
end        
%respect the conditon E[|s|^2]=In
variance_signal=mean(abs(alphabet_modulation).^2);
alphabet_modulation=alphabet_modulation/sqrt(variance_signal);

%Construct all the possible combination of n symbols
S=PermsRep(alphabet_modulation,n);
Nb_S=size(S,1);

%% Precompute the term (I\otimes H_bar)Mc to reduce computation time
Mc=[];
for column_number=1:l
        column_index=1+2*(column_number-1)*nt:2*column_number*nt;
        A_real_temp2(:,:)=real(A(column_number,:,:));
        A_imag_temp2(:,:)=imag(A(column_number,:,:));
        Mc=[Mc;A_real_temp2;A_imag_temp2];
end
H_bar=[real(H) -imag(H);imag(H) real(H)];
I_H_MC=kron(eye(l),H_bar)*Mc;


%% Compute the log likelihood
log_LF=0;
    
Nb=N/l;
log_LF=-((n*Nb)*log(Mod_State_Nb)+nr*N*log(pi*noise_variance)); %compared to equation (18); i've developped the logarithm

%compute the first sum with a for loop
for block_number=0:Nb-1
       
        index=[(block_number*l):((block_number+1)*l-1)]+1; %index of the columns for the vth block
        Yv=received_signal(:,index);
        yv_tilde=[real(Yv);imag(Yv)];
        yv_tilde=yv_tilde(:);
        
        %compute the second sum with a for loop
        inner_term=0;
        for combinaison_element=1:Nb_S
            s=S(combinaison_element,:);
            s_tilde=[real(s(:)); imag(s(:))];
            inner_term=inner_term+exp(-1*(sum(sum(abs(yv_tilde-I_H_MC*s_tilde).^2)))/noise_variance);   %see equation (18)
        end 
        log_LF=log_LF+log(inner_term);                                                                  %see equation (18)
end    
    
%Remark: If C is an OSTBC, the LF can be easily simplified
%(I_H_MC*I_H_MC=alpha*I)