function mixing = kmeans(data)
%it is help of the function mixing to estimate the mixing matrix from two
%channel mixtures by K-means clustering algorithm
%  Tao Xu and Wenwu Wang
%  CVSSP, Univeristy of Surrey
%  Aug. 2010

%The two mixtures of four speech signals
%data=wavread('dev1_female4_inst_mix.wav');

%transform the two mixtures to TF-domain
f1=data(:,1);
f2=data(:,2);

%f1=f1(1:61440);
%f2=f2(1:61440);

f1=f1(1:160000);
f2=f2(1:160000);

window = hamming(512);
F1 = spectrogram(f1,window,256,511,8000);
F2 = spectrogram(f2,window,256,511,8000);
%F1 = real(F1);
%F2 = real(F2);
F1 = abs(F1);
F2 = abs(F2);
f1=abs(f1);
f2=abs(f2);
figure(1);plot(f1,f2,'.');
%gtext('x_2')
xlabel('x_1'),ylabel('x_2');
%TF_data = zeros(14080,2);
%TF_data(:,1) = reshape(F1,61184,1);
%TF_data(:,2) = reshape(F2,61184,1);

TF_data(:,1) = reshape(F1,159744,1);
TF_data(:,2) = reshape(F2,159744,1);

%display the scatter of mixtures in the TF plane
x1=TF_data(:,1);
x2=TF_data(:,2);

figure(2);plot(x1,x2,'.');
xlabel('x_1'),ylabel('x_2');

%normalize data to a circle
a=sqrt(x1.^2 + x2.^2);
A=[a,a];

data=[x1,x2];
data=data./A;


%set the pimitive value of centres
%%  nomally
%centres=[ mean(data(1:15360,:),1); mean(data(15361:30720,:),1); mean(data(30721:33000,:),1); mean(data(33001:61184,:),1) ];
%%  for dev1
%%centres=[ mean(data(1:20000,:),1); mean(data(10001:30720,:),1); mean(data(20721:40000,:),1); mean(data(30001:61184,:),1) ];
%%  for test
%centres=[ mean(data(1:30000,:),1); mean(data(20001:40720,:),1); mean(data(30721:50000,:),1); mean(data(40001:61184,:),1) ];



%female
centres=[ mean(data(1:140000,:),1); mean(data(40001:80000,:),1); mean(data(80001:150000,:),1); mean(data(120001:159744,:),1) ];
%male
%centres=[ mean(data(1:20000,:),1); mean(data(40001:80000,:),1); mean(data(10001:150000,:),1); mean(data(70001:159744,:),1) ];


%set some variables
ndata = size(data,1);
ncentre = size(centres,1);
assignments = zeros(ndata,1)';


%estimate the centres by K-means
while 1
  
  %call the function dist2(calculate the eula distance between data and centres)
  d2 = dist2(data,centres); 
  
  %find the closest centre
  [minvals, index] = min(d2',[],1);
  
  %the condition of procedure completing
  if index==assignments
   
     break;

  else
  
     assignments = index;

  end

  %recalculate the centres
  for c=1:ncentre

    centres(c,:)=mean(data(find(index==c),:),1);

  end

  centres %print the iteration process of estimating mixing matrix

end

%display the estimate mixing matrix
mixing=centres'