//READING 
  stacksize(2e7);
  temp1=wavread('e:\momo\bss\s1.wav');
  s1=temp1(5000:11000);
  temp2=wavread('e:\momo\bss\s1.wav');
  s2=temp2(5000:11000);
  
  sources=[s1';s2'];

  [N,signal_length]=size(sources);         //caculating the number of sources to be mixed
  xset("window",1);
  xname("Source sounds");
  for n=1:N,               // Ploting the beginning of the N source signals  
    subplot(N,1,n); plot(sources(n,:));
    if n==1, xtitle('The beginning of the speech signals'); end;
  end 

//MIXING
  
//generating mixing matrics,remember it should not be singular
  rand("uniform");
  flag=1;
  epsilon=1e-1;
  while flag==1
    MIX_MATRIX=rand(N,N);
    if abs(det(MIX_MATRIX))<=epsilon then, flag=1; else flag=0; end
  end  
  mprintf(['\n***The source signals are linearily mixed together using the matrix:***\n']);
  A=MIX_MATRIX;
  //A=[.3321719  .5015342 ; .5935095  .4368588]
  //A=[.6723950  .3911574 ; .2017173  .8300317]
  A=[.3873779     .9488184 
     .9222899     .3435337] 

//Mixed signals  
  mixing=A*sources;    
  xset("window",2);
  xname("Mixed sounds");
  for n=1:N,               // Plotting the mixed signals  
    subplot(N,1,n); plot(mixing(n,:));
    if n==1, xtitle('The mixed signals'); end;
  end

//***********************************Prewhiten!*******************************************
  for i=1:N,
    mixing(i,:)=mixing(i,:)-mean(mixing(i,:));
  end
  cov=(mixing*mixing')/signal_length;
  [D,E]=bdiag(cov);
  mixing=E*inv(sqrt(D))*E'*mixing;

//************************************Begin!*********************************************** 

//initialize  
  ite=50;           //ite is the maximum Number of Iterations
  Er=[];            //Recording the perfomance index during the learning process
  W=eye(N,N);
  L=0.00004;        //learning step
  I=eye(N,N);                          

//training
  for k=1:ite,                                 
    Y=W*mixing;
    W=W+L*(I-2*tanh(Y)*Y')*W;
//************************************************************************************
    p=W*E*inv(sqrt(D))*E'*A;
    e=0;
    for i=1:N
       e=e+sum(abs(p(i,:)))/max(abs(p(i,:)))-1;
       e=e+sum(abs(p(:,i)))/max(abs(p(:,i)))-1;
    end
    Er=[Er,e];
//*************************** Er Caculating end!! **************************************                                       
  end

  Y=W*mixing;
  mprintf(['\n***The perfomance matrix p:***\n']);
  p=W*E*inv(sqrt(D))*E'*A
  xset("window",3);
  xname("Unmixed sounds");
  for n=1:N,               // Ploting the Unmixed speech signals 
    subplot(N,1,n); plot(Y(n,:));
    if n==1, xtitle('The unmixed speech signals'); end;
  end 
//*********************************** plot Er *********************************************
  xset("window",4);
  xname("Error measure Er");
  plot(Er);
  xtitle("Error measure Er during the learning process");