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

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

%HELP: subspace_channel_estimation_STBC 
%
%Blind Subspace channel identification for MIMO-STBC systems. The method
%uses the knowledge of the Space time code to extract the channel. The
%method doesn't work for non-identifiable STBCs (see publication [1] for
%details). For identifiable STBCs, there is one 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,...)
%Ouput: - channel
%       - separated data
%
%Reference: 
%
%[1] Ammar, N.; Ding, Z. "Blind Channel Identifiability for Generic
%Linear Space-Time Block Codes" IEEE Transactions on Signal Processing
%Volume 55, Issue 1, Jan. 2007 Page(s):202 - 217
%
%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);
[nb_receiver,N]	= size(received_signal); 
Nb_bloc=N/code_length;

%extract space time codign information

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(:,indice)=Ak(:);
    C(:,indice+nb_symbols_block)=Bk(:);
end

%% rewrite the receiving sample
R=[];
for indice=1:Nb_bloc
   received_bloc=received_signal(:,(indice-1)*code_length+1:indice*code_length);
   R(:,indice)=received_bloc(:);
end    
size(R);

%% extract a basis NL for the left null space of R without noise
[U,D,V]=svd(R);
NL=U(:,2*nb_symbols_block+1:end);

%% Construct the Phi matrix
Big_E=[];
for indice=1:code_length
   %contruct Ek
   Ek=[zeros(nb_receiver*(indice-1),nb_receiver);...
       eye(nb_receiver);...
       zeros(nb_receiver*(code_length-indice),nb_receiver)];
   Big_E=[Big_E;kron(eye(nb_emitters),Ek)];    
end   

%% perform channel estimation
phi=kron(C.',NL')*Big_E;
[U2,D2,V2]=svd(phi);
h_est=V2(:,nb_receiver*nb_emitters);
channel=reshape(h_est,nb_receiver,nb_emitters);

%Display the result
if verbose==1
   transmitted_signal=pinv(channel)*received_signal;
   for indice=1:nb_emitters
        scatterplot(transmitted_signal(indice,:));
        title(sprintf('Extracted signal nb %d',indice));
        grid;
   end 
end