tito_main.m

用MATLAB实现MIMO系统盲辨识

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                         TITO_main.m
% This program gives a demo of the algorithm proposed in the following paper. 
% Two source signals were mixed by a Two-Input-Two-Output (TITO) dynamic channel.
% Frequency domain method was used to estimate the corresponding channel
% response blindly. 
% An Iterative algorithm was adopted here to estimate the 
% correct H1(w) and H2(w) through R and S, by modifying E(w).
%
% Reference: 
% K. Diamantaras, A. Petropulu and B. Chen, "Blind Two-Input-Two-Output 
% FIR Channel Identification Based on Second-Order Statistics,"
% IEEE Transactions on Signal Processing, vol. 48, No. 2, pp. 534-542, 
% Feb. 2000.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


TIME_OF_PROGRAM_START = clock;    %%% 

format compact;
i=sqrt(-1);

%%% Rx=E{x(w)x'(w)}
Rx_cond_selection_percentage=90;  %%% Using the top 90% frequencies based on the condition
                                  %%% number of the matrix Rx(w)

Minimum_Phase_Channel=0;          %%% Whether the TITO system is minimum phase
                                  %%% If it is, we can resonstruct it based on the
                                  %%% amplitude response alone, if it is not, we need
                                  %%% use both amplitude and phase response.
                                  
SNR=30;                           %%% Signal to Noise ratio
RUN_TIMES=5;                      %%% Monte-Carlo running times.
ITER_TIMES=1;                     %%% Iteration times for the EE(w) in the paper.
L = 3;                            %%% Mixing Channel length
n = 2;                            %%% Number of inputs/outputs
N = 4096*2;                       %%% Input/Output data length
NF = 128;                         %%% Size of FFT used in the estimation.
w1=1;                             %%% Frequency dufference used in constructing the Rx1
                                  %%% Rx1(w, w+w1)=E{x(w)x'(w+w1)}

Iter_HH2_abs=0;

MSE_p_h1=zeros(1,RUN_TIMES);      %%% MSE of the estimated h1 before iteration
MSE_p_h2=zeros(1,RUN_TIMES);      %%% MSE of the estimated h2 before iteration

MSE_h1=zeros(ITER_TIMES,RUN_TIMES);%%% MSE of the estimated h2 after iteration
MSE_h2=zeros(ITER_TIMES,RUN_TIMES);%%% MSE of the estimated h2 after iteration

S_his = zeros(RUN_TIMES,NF/2+1);
R_his = zeros(RUN_TIMES,NF/2+1);

window = kaiser(NF,12);           %%% one dimensional window function
winmat = window*window';          %%% two dimensional window function

F = exp(-i*2*pi*(0:NF-1)'*(0:NF-1)/NF); %%% Fourier matrix

h=zeros(L,n,n);
h(:,1,1) = [1, zeros(1,L-1)];
h(:,2,2) = [1, zeros(1,L-1)];
h(:,1,2)=[1.0000    0.345    0.1];
h(:,2,1)=[1.0000   -0.5766    0.6575];

%h(:,1,2)= [1.0000   -0.6420   -0.1948    0.2082   -0.8049];
%h(:,2,1)=[1.0000    1.6038   -0.6740    0.0091    1.8725];

%%%%% Minimum Phase channel
%h(:,1,2) = [1.0000    0.3998   -0.8962    0.0339   -0.0107   -0.3405];
%h(:,2,1) = [1.0000    0.7557    0.4338    0.1477   -0.0431   -0.1299];

h2=h(:,2,1)';
h1=h(:,1,2)';
H2=fft(h2,NF);
H1=fft(h1,NF);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%  The Monte-Carlo test begin   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iloop=1:RUN_TIMES
   
   iloop
   rand('state',sum(100*clock))
   s = zeros(n,N);
   
   g(1,:) = [ -0.5327    0.5594   -0.1026    0.9871   -0.5474    1.3319   -1.8898   -0.3999   -0.3299    0.2608   -0.5087   -0.7627   -0.3299   -0.5442   -1.0783   -0.6630    1.1211    1.9612    1.8482    0.0176    1.2493   -1.0560   -2.2745   -0.9096    2.0539    0.7513    1.1254    0.2137   -0.8023    0.0878];
   g(2,:) = [ -0.0197    0.6498   -0.7554    0.9009    0.2750    1.5394    2.1862    0.2524   -0.2311    1.4686   -0.6528    2.0335    1.8471   -0.0590   -0.0193   -2.3757    0.9451    0.0510   -0.5135   -0.8881    0.5549   -2.0494    1.5916    0.3167   -0.1147   -1.9982    0.1660    1.2596   -0.6778    0.7646];
   
   for (ii=1:n)
      s(ii,:)  = randn(1,N);
   end
   s = s-(mean(s'))'*ones(1,N);
   s=orth(s')';
   
   for (ii=1:n)
      s(ii,:) = filter(g(ii,:),1,s(ii,:));
   end
   
   s = s-(mean(s'))'*ones(1,N);
   s(1,:) = s(1,:)/std(s(1,:));
   s(2,:) = s(2,:)/std(s(2,:));
   
   %
   %Multichannel mixing
   %===================
   x = zeros(size(s));     %Observed mixture
   y = zeros(size(s));     %Observed mixture
   H = randn(NF,n,n);      %Mixing matrix
   %H(frequency, signalIndex, signalIndex)
   HH = randn(n,n,NF);     %HH(signalIndex, signalIndex, frequency)
   for (ii=1:n)
      for (jj=1:n)
         x(ii,:) = x(ii,:) + filter(h(:,ii,jj),1,s(jj,:));
         H(:,ii,jj) = fft(h(:,ii,jj), NF);
      end
   end
   
   noise=randn(1,N);
   x(1,:) = x(1,:)+noise/((10^(SNR/20))*std(noise)/std(x(1,:)));
   noise=randn(1,N);
   x(2,:) = x(2,:)+noise/((10^(SNR/20))*std(noise)/std(x(2,:)));
   
   for (w=1:NF)
      for (ii=1:n)
         for (jj=1:n)
            HH(ii,jj,w) = H(w,ii,jj);
         end
      end
   end
   
   %
   %rx(t) estimation (t = time lag, rx = matrix)
   %=================================================
   rx = zeros(2*NF-1,n,n);
   for (ii=1:n)
      for (jj=1:n)
         for (lag=-NF+1:-1)
            rx(lag+NF,ii,jj) = x(ii,-lag+1:N)*x(jj,1:N+lag)'/N;
         end
         for (lag=0:NF-1)
            rx(lag+NF,ii,jj) = x(ii,1:N-lag)*x(jj,lag+1:N)'/N;
         end
      end
   end
   
   %
   %Rx(w) estimation (w = frequency, Rx = matrix)
   %=================================================
   Rx = zeros(n,n,NF);                     % = Rx(w,w)
   Rx1 = zeros(n,n,NF/2+1);        % = Rx(w,w+1)
   
   for (ii=1:n)
      for (jj=1:n)
         RR = toeplitz(rx(NF:-1:1,ii,jj),rx(NF:2*NF-1,ii,jj));
         RR = winmat.*RR;
         for (w=1:NF)
            if (jj==ii)
               Rx(ii,jj,w) = real(F(w,:)*RR*F(w,:)');
            else
               Rx(ii,jj,w) = F(w,:)*RR*F(w,:)';
            end                  
         end
         for (w=1:NF/2+1)
            Rx1(ii,jj,w) = F(w,:)*RR*F(w+w1,:)';
         end
      end
   end
   
   for (w=1:NF)
      Rx_eig(:,w)=eig(Rx(:,:,w));
      Rx_cond(w,1)=real(max(Rx_eig(:,w))/min(Rx_eig(:,w)));
   end
   
   [Rx_cond_sort Rx_cond_sort_index]=sort(Rx_cond);
   
   Rx_cond_max=Rx_cond(Rx_cond_sort_index(floor(NF*Rx_cond_selection_percentage/100)));
   
   Rx_cond_Index=[];
   for (w=1:NF)
      if Rx_cond(w,1) < Rx_cond_max
         Rx_cond_Index=[Rx_cond_Index w];
      end
   end
   
   %
   %Prewhitening
   %==============
   V = zeros(n,n,NF);              %V(signalIndex, signalIndex, frequency)
   for (w=1:NF)
      V(:,:,w) = inv(sqrtm(Rx(:,:,w)));
   end
   
   for (ii=1:n)
      for (jj=1:n)
         v(:,ii,jj) = fftshift(ifft(V(ii,jj,:), NF));
      end
   end
   
   for (ii=1:n)
      for (jj=1:n)
         y(ii,:) = y(ii,:) + filter(v(:,ii,jj),1,x(jj,:));
      end
   end
   
   ry = zeros(2*NF-1,n,n);
   for (ii=1:n)
      for (jj=1:n)
         for (lag=-NF+1:-1)
            ry(lag+NF,ii,jj) = y(ii,-lag+1:N)*y(jj,1:N+lag)'/N;
         end
         for (lag=0:NF-1)
            ry(lag+NF,ii,jj) = y(ii,1:N-lag)*y(jj,lag+1:N)'/N;
         end
      end                 
   end
   
   for (ii=1:n)
      for (jj=1:n)
         RR = toeplitz(ry(NF:-1:1,ii,jj),ry(NF:2*NF-1,ii,jj));
         RR = winmat.*RR;
         
         for (w=1:NF)
            if (jj==ii)
               Ry(ii,jj,w) = real(F(w,:)*RR*F(w,:)');
            else
               Ry(ii,jj,w) = F(w,:)*RR*F(w,:)';
            end
         end
      end
   end
   
   % Compute Ry(w,w+1), Ry(w,w+2) (w = frequency, Ry = matrix)
   %============================================================
   Ry1 = zeros(n,n,NF/2+1);
   for (w=1:NF/2+1)
      Ry1(:,:,w) = V(:,:,w)*Rx1(:,:,w)*V(:,:,w+w1)';
   end
   
   Cy1 = zeros(n,n,NF/2+1);
   for (w=1:NF/2+1)
      Cy1(:,:,w) = Ry1(:,:,w)*Ry1(:,:,w)';
   end
   
   %