[channel]=SOS_channel_estimation_OSTBC(received_signal,code_name,rate,num_code)

function [channel]=SOS_channel_estimation_OSTBC(received_signal,code_name,rate,num_code)

%HELP: channel estimation OSTBC
%Determine the channel matrix for a space time coded communication using Orthogonal STBC. 
%Channel estimation is performed through the eigenvalue of a modified covariance matrix[1,2].
%The method doesn't work for non-identifiable OSTBCs including Alamouti code (see publication [1,2]
%for details). For identifiable OSTBCs, there is one eal scalar ambiguity after
%channel estimation. 
%
%Input: - received signal
%       - code_name: The code name ('Alamouti','OSTBC3',...)
%       - code_rate: The code rate ('1','1/2',...)
%       - num_code: The code number (1,2,...)
%Output: - channel
%       - separated data
%
%Reference: 

% [1] S. Shahbazpanahi, A.B. Gershman, and J.H. Manton, "Closed
% form blind mimo channel estimation for othogonal space-time
% codes," IEEE Transactions on Signal Processing, vol. 53, no.
% 12, pp. 4506{4517, 2005.
% [2] J. Va and I. Santamara, "On the blind identifiability of 
% orthogonal space-time block codes from second order statistics",
% IEEE Transactions on Information Theory, vol.54, no.2, pp.709-722,
% February 2008 
%
%Programmed by V. Choqueuse ([e-mail])

verbose=1;


%% extract space time block coding information
Rendement=str2num(rate);
[nb_emitters,code_length]=size(space_time_coding(0,code_name,rate,num_code,1));
nb_symbols_block=code_length*str2num(rate);


for indice=1:nb_symbols_block
    real_part=zeros(1,nb_symbols_block);
    real_part(indice)=1;
    imag_part=zeros(1,nb_symbols_block);
    imag_part(indice)=i; 
    Ak=space_time_coding(real_part,code_name,rate,num_code).';  
    Bk=space_time_coding(imag_part,code_name,rate,num_code).';   
    C_tilde(indice,:,:)=[real(Ak) -imag(Ak);imag(Ak) real(Ak)];
    C_tilde(indice+nb_symbols_block,:,:)=[real(Bk) -imag(Bk);imag(Bk) real(Bk)];
    
end

%% Construct the matrix Z
[nb_receiver,N]	= size(received_signal); 
Nb_bloc=N/code_length;
Z=[];
for indice=1:Nb_bloc;
    
   Y_bloc=received_signal(:,(indice-1)*code_length+1:indice*code_length).';
   Z_temp=zeros(2*nb_symbols_block, 2*nb_emitters*nb_receiver);
   
   for indice_ligne=1:2*nb_symbols_block 
       C_tilde_temp(:,:)=C_tilde(indice_ligne,:,:);
       for indice_colonne=1:nb_receiver 
           y_tilde=[real(Y_bloc(:,indice_colonne)); imag(Y_bloc(:,indice_colonne))];
           Z_temp(indice_ligne,(indice_colonne-1)*2*nb_emitters+1:indice_colonne*2*nb_emitters)=...
                       y_tilde'*C_tilde_temp;
       end         
   end
   Z=[Z;Z_temp];
end    

%% perform channel estimation
[U,V]=eig(Z'*Z);
[V,V_sort]=sort(diag(V),'descend');
h_est=U(:,V_sort(1));
channel=reshape(h_est,2*nb_emitters,nb_receiver);
channel=channel(1:nb_emitters,:)+j*channel(nb_emitters+1:end,:);
channel=channel.';

separated_data=Z*h_est;
separated_data=reshape(separated_data,2*nb_symbols_block,Nb_bloc);
separated_data=[eye(nb_symbols_block) i*eye(nb_symbols_block)]*separated_data;   

%Display the result
if verbose==1
   for indice=1:nb_symbols_block
        scatterplot(separated_data(indice,:));
        title(sprintf('Extracted signal nb %d',indice));
        grid;
   end 
end