diff --git a/GaussMix/ClusterNormalize.m b/GaussMix/ClusterNormalize.m
new file mode 100755
index 0000000000000000000000000000000000000000..d6165f9b5d0b3d12eeb0728505d8d72eddf799d0
--- /dev/null
+++ b/GaussMix/ClusterNormalize.m
@@ -0,0 +1,11 @@
+function mixture1=ClusterNormalize(mixture)
+
+ s = sum([mixture.cluster(:).pb]);
+ for k=1:mixture.K
+ mixture.cluster(k).pb = mixture.cluster(k).pb/s;
+ mixture.cluster(k).invR = inv(mixture.cluster(k).R);
+ mixture.cluster(k).const = -(mixture.M*log(2*pi) +log(det(mixture.cluster(k).R)))/2;
+ end
+ mixture1=mixture;
+
+end
diff --git a/GaussMix/CreerData.m b/GaussMix/CreerData.m
new file mode 100755
index 0000000000000000000000000000000000000000..1498d7929a0a00b193c87878ec03c2335a4c5364
--- /dev/null
+++ b/GaussMix/CreerData.m
@@ -0,0 +1,30 @@
+% close all;
+% clear all;
+% clc;
+
+x = zeros(1,100);
+for i=1:100
+ u = rand(1,1);
+ if (u < 0.3)
+ x(i) = 7+randn(1,1);
+ else if (u<0.8)
+ x(i) = randn(1,1);
+ else
+ x(i) = 4+0.5*randn(1,1);
+ end;
+end;
+end;
+
+I= hist(x,100);
+figure;
+plot(I);
+x=x'
+
+[mtrs,omtr] = GaussianMixture(x, 4, 0, 0);
+disp(sprintf('\toptimal order K*: %d', omtr.K));
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
\ No newline at end of file
diff --git a/GaussMix/EMIterate.m b/GaussMix/EMIterate.m
new file mode 100755
index 0000000000000000000000000000000000000000..f2d95f8e8a3ea991543f8251c793b8dab3d773f1
--- /dev/null
+++ b/GaussMix/EMIterate.m
@@ -0,0 +1,30 @@
+function mixture1 = EMIterate(mixture, pixels)
+% mixture1 = EMIterate(mixture, pixels)
+% perform the EM algorithm with a preassigned fixed order K
+%
+% mixture: a mixture structure pre-initialized
+% mixture.K - order of the mixture
+% mixture.M - dimension of observations
+% mixture.cluster - an array of K cluster structure
+% pixels: a N x M observation matrix, each row is an observation vector
+%
+% mixture1: converged result
+% mixture1.rissanen: MDL(K)
+%
+
+[N M] = size(pixels);
+Lc = 1+M+0.5*M*(M+1);
+epsilon = 0.01*Lc*log(N*M);
+[mixture, llnew] = EStep(mixture, pixels);
+while true
+ llold = llnew;
+ mixture = MStep(mixture, pixels);
+ [mixture, llnew] = EStep(mixture, pixels);
+ if (llnew-llold)<=epsilon
+ break;
+ end
+end
+% mixture = rmfield(mixture, 'pnk');
+mixture.rissanen = -llnew+0.5*(mixture.K*Lc-1)*log(N*M);
+mixture.loglikelihood = llnew;
+mixture1 = mixture;
diff --git a/GaussMix/EStep.m b/GaussMix/EStep.m
new file mode 100755
index 0000000000000000000000000000000000000000..d6a5149e87aa3f869321de02a5c9c75027a3a8b6
--- /dev/null
+++ b/GaussMix/EStep.m
@@ -0,0 +1,29 @@
+function [mixture1 , likelihood] = EStep(mixture, pixels)
+% EStep perform the E-step of the EM algorithm
+% 1) calculate pnk = Prob(Xn=k|Yn=yn, theta)
+% 2) calculate likelihood = log ( prob(Y=y|theta) )
+
+ [N M] = size(pixels);
+ K=mixture.K;
+ pnk=zeros(N,K);
+ for k=1:K
+ Y1=pixels-ones(N,1)*mixture.cluster(k).mu';
+ Y2=-0.5*Y1*mixture.cluster(k).invR;
+ pnk(:,k) = dot(Y1,Y2,2)+mixture.cluster(k).const;
+ end
+ llmax=max(pnk,[],2);
+ pnk =exp( pnk-llmax*ones(1,K) );
+ pnk = pnk.*(ones(N,1)*[mixture.cluster(:).pb]);
+ ss = sum(pnk,2);
+ likelihood = sum(log(ss)+llmax);
+ pnk = pnk./(ss*ones(1,K));
+ mixture1 = mixture;
+ mixture1.pnk = pnk;
+
+
+% function ll = loglike(yn, cluster)
+% Y = yn-cluster.mu;
+% ll = -0.5*Y'*cluster.invR*Y+cluster.const;
+% end
+%
+end
\ No newline at end of file
diff --git a/GaussMix/GMClassLikelihood.m b/GaussMix/GMClassLikelihood.m
new file mode 100755
index 0000000000000000000000000000000000000000..08b36bd42e99499ca1c0caa93a509ffc7eec603b
--- /dev/null
+++ b/GaussMix/GMClassLikelihood.m
@@ -0,0 +1,26 @@
+function ll = GMClassLikelihood(mixture, Y)
+% ll = GMClassLikelihood(class, Y)
+% GMClassLikelihood calculate the log-likelihood of data vectors assuming
+% they are generated by the given Gaussian Mixture
+%
+% mixture: a structure representing a Gaussian mixture
+% Y: a NxM matrix, each row is an observation vector of dimension M
+%
+% ll: a Nx1 matrix with the n-th entry returning the log-likelihood of the
+% n-th observation
+%
+
+ [N, M] = size(Y);
+ K = mixture.K;
+ pnk=zeros(N,K);
+ for k=1:K
+ Y1=Y-ones(N,1)*mixture.cluster(k).mu';
+ Y2=-0.5*Y1*mixture.cluster(k).invR;
+ pnk(:,k) = dot(Y1,Y2,2)+mixture.cluster(k).const;
+ end
+ llmax=max(pnk,[],2);
+ pnk =exp( pnk-llmax*ones(1,K) );
+ pnk = pnk.*(ones(N,1)*[mixture.cluster(:).pb]);
+ ss = sum(pnk,2);
+ ll = log(ss)+llmax;
+end
diff --git a/GaussMix/GaussianMixture.m b/GaussMix/GaussianMixture.m
new file mode 100755
index 0000000000000000000000000000000000000000..394429f0de1ad741eb994add9711b4844a1d9916
--- /dev/null
+++ b/GaussMix/GaussianMixture.m
@@ -0,0 +1,71 @@
+function [mixture, optmixture] = GaussianMixture(pixels, initK, finalK, verbose)
+% [mixture, optmixture] = GaussianMixture(pixels, initK, finalK)
+% perform the EM algorithm to estimate the order, and
+% parameters of a Gaussian Mixture model for a given set of
+% observation.
+%
+% pixels: a N x M matrix of observation vectors with each row being an
+% M-dimensional observation vector, totally N observations
+% initK: the initial number of clusters to start with and will be reduced
+% to find the optimal order or the desired order based on MDL
+% finalK: the desired final number of clusters for the model.
+% Estimate the optimal order if finalK == 0.
+% verbose: true/false, return clustering information if true
+%
+% mixture: an array of mixture structures with each containing the
+% converged Gaussian mixture at a given order
+% mixture(l).K: order of the mixture
+% mixture(l).M: dimension of observation vectors
+% mixture(l).rissanen: converaged MDL(K)
+% mixture(l).loglikelihood: ln( Prob{Y=y|K, theta*} )
+% mixture(l).cluster: an array of cluster structures with each
+% containing the converged cluster parameters
+% mixture(l).cluster(k).pb: pi(k)=Prob(Xn=k|K, theta*)
+% mixture(l).cluster(k).mu: mu(k), mean vector of the k-th cluster
+% mixture(l).cluster(k).R: R(k), covariance matrix of the k-th cluser
+% optmixture: one of the element in the mixture array.
+% If finalK > 0, optmixture = mixture(1) and is the mixture with order finalK.
+% If finalK == 0, optmixture is the one in mixture with minimum MDL
+%
+
+if ~isnumeric(initK) || ~all(size(initK)==[1,1]) || initK<=0 || mod(initK,1)~=0
+ error('GaussianMixture: initial number of clusters initK must be a positive integer');
+end
+if ~isnumeric(finalK) || ~all(size(finalK)==[1,1]) || finalK<0 || mod(finalK,1)~=0
+ error('GaussianMixture: final number of clusters finalK must be a positive integer or zero');
+end
+if finalK > initK
+ error('GaussianMixture: finalK cannot be greater than initK');
+end
+if ~isa(pixels,'double')
+ pixels = double(pixels);
+end
+mtr = initMixtureRange(pixels, initK); % RANDOMIZED INIT
+mtr = EMIterate(mtr, pixels);
+if verbose
+ disp(sprintf('K: %d\t rissanen: %f', mtr.K, mtr.rissanen));
+end
+mixture(mtr.K-max(1,finalK)+1) = mtr;
+while mtr.K > max(1, finalK)
+ mtr = MDLReduceOrder(mtr, verbose);
+ mtr = EMIterate(mtr, pixels);
+ if verbose
+ disp(sprintf('K: %d\t rissanen: %f', mtr.K, mtr.rissanen));
+ end
+ mixture(mtr.K-max(1,finalK)+1) = mtr;
+end
+
+if finalK>0
+ optmixture = mixture(1);
+else
+ minriss = mixture(length(mixture)).rissanen; optl=length(mixture);
+ for l=length(mixture)-1:-1:1
+ if mixture(l).rissanen < minriss
+ minriss = mixture(l).rissanen;
+ optl = l;
+ end
+ end
+ optmixture = mixture(optl);
+end
+
+
diff --git a/GaussMix/MDLReduceOrder.m b/GaussMix/MDLReduceOrder.m
new file mode 100755
index 0000000000000000000000000000000000000000..8702faa6fd4ab24768b8d29b5e979ed3c9dcd84c
--- /dev/null
+++ b/GaussMix/MDLReduceOrder.m
@@ -0,0 +1,45 @@
+function mixture1 = MDLReduceOrder(mixture, verbose)
+% MDLReduceOrder reduce the order of the mixture by 1 by combining the
+% two closest distance clusters
+
+ K=mixture.K;
+ for k1=1:K
+ for k2=k1+1:K
+ dist = distance(mixture.cluster(k1), mixture.cluster(k2));
+ if (k1==1 && k2==2) || (dist < mindist)
+ mink1=k1; mink2=k2;
+ mindist = dist;
+ end
+ end
+ end
+ if verbose
+ disp(['combining cluster: ' int2str(mink1) ' and ' int2str(mink2)]);
+ end
+ cluster = mixture.cluster;
+ cluster(mink1) = addCluster(cluster(mink1), cluster(mink2));
+ cluster(mink2:(K-1)) = cluster((mink2+1):K);
+ cluster = cluster(1:(K-1));
+ mixture1 = mixture;
+ mixture1.cluster = cluster;
+ mixture1.K = K-1;
+ mixture1 = ClusterNormalize(mixture1);
+
+ function d = distance(cluster1, cluster2)
+ cluster3 = addCluster(cluster1, cluster2);
+ d = cluster1.N*cluster1.const + cluster2.N*cluster2.const-cluster3.N*cluster3.const;
+ end
+
+ function cluster3 = addCluster(cluster1, cluster2)
+ wt1 = cluster1.N/(cluster1.N+cluster2.N);
+ wt2 = 1-wt1;
+ M = length(cluster1.mu);
+ cluster3 = cluster1;
+ cluster3.mu = wt1*cluster1.mu+wt2*cluster2.mu;
+ cluster3.R = wt1*(cluster1.R+(cluster3.mu-cluster1.mu)*(cluster3.mu-cluster1.mu)') ...
+ + wt2*(cluster2.R+(cluster3.mu-cluster2.mu)*(cluster3.mu-cluster2.mu)');
+ cluster3.invR = inv(cluster3.R);
+ cluster3.pb = cluster1.pb+cluster2.pb;
+ cluster3.N = cluster1.N+cluster2.N;
+ cluster3.const = -(M*log(2*pi) +log(det(cluster3.R)))/2;
+ end
+end
\ No newline at end of file
diff --git a/GaussMix/MStep.m b/GaussMix/MStep.m
new file mode 100755
index 0000000000000000000000000000000000000000..b14cc96140a88b7a5d424aedb34daf02ec2c3369
--- /dev/null
+++ b/GaussMix/MStep.m
@@ -0,0 +1,24 @@
+function mixture1 = MStep(mixture, pixels)
+% MStep perform the M-step of the EM algorithm
+% from the pnk calculated in the E-step, update parameters of each cluster
+%
+
+ for k=1:mixture.K
+ mixture.cluster(k).N = sum(mixture.pnk(:,k));
+ mixture.cluster(k).pb = mixture.cluster(k).N;
+ mixture.cluster(k).mu = (pixels' * mixture.pnk(:,k))/mixture.cluster(k).N;
+ for r=1:mixture.M
+ for s=r:mixture.M
+ mixture.cluster(k).R(r,s) = ((pixels(:,r)-mixture.cluster(k).mu(r))' ...
+ * ((pixels(:,s)-mixture.cluster(k).mu(s)).*mixture.pnk(:,k))) ...
+ /mixture.cluster(k).N;
+ if r~=s
+ mixture.cluster(k).R(s,r) = mixture.cluster(k).R(r,s);
+ end
+ end
+ end
+ mixture.cluster(k).R = mixture.cluster(k).R+mixture.Rmin*eye(mixture.M);
+ end
+ mixture1=ClusterNormalize(mixture);
+end
+
diff --git a/GaussMix/RStep.m b/GaussMix/RStep.m
new file mode 100755
index 0000000000000000000000000000000000000000..81680acb804bb5abd4d70e7749df10c6bc30fd7c
--- /dev/null
+++ b/GaussMix/RStep.m
@@ -0,0 +1,22 @@
+function mixture1 = RStep(mixture, pixels)
+% MStep perform the M-step of the EM algorithm
+% from the pnk calculated in the E-step, update parameters of each cluster
+%
+ mixture.M=size(pixels,2);
+ for k=1:mixture.K
+ mixture.cluster(k).mu = (pixels' * mixture.pnk(:,k))/mixture.cluster(k).N;
+ for r=1:mixture.M
+ for s=r:mixture.M
+ mixture.cluster(k).R(r,s) = ((pixels(:,r)-mixture.cluster(k).mu(r))' ...
+ * ((pixels(:,s)-mixture.cluster(k).mu(s)).*mixture.pnk(:,k))) ...
+ /mixture.cluster(k).N;
+ if r~=s
+ mixture.cluster(k).R(s,r) = mixture.cluster(k).R(r,s);
+ end
+ end
+ end
+ mixture.cluster(k).R = mixture.cluster(k).R+mixture.Rmin*eye(mixture.M);
+ end
+ mixture1=ClusterNormalize(mixture);
+end
+
diff --git a/GaussMix/SplitClasses.m b/GaussMix/SplitClasses.m
new file mode 100755
index 0000000000000000000000000000000000000000..5d3aabfdf9d7be313bd22e513e5f93199190bf91
--- /dev/null
+++ b/GaussMix/SplitClasses.m
@@ -0,0 +1,21 @@
+function classes = SplitClasses(mixture)
+% classes = SplitClasses(mixture)
+% SplitClasses splits the Gaussian mixture with K subclasses into K
+% Gaussian mixtures, each of order 1 containing each of the subclasses
+%
+% mixture: a structure representing a Gaussian mixture of order
+% mixture.K
+%
+% classes: an array of structures, with each representing a Gaussian
+% mixture of order 1 consisting of one of the original subclasses
+%
+
+K=mixture.K;
+for k=1:K
+ classes(k) = mixture;
+ classes(k).K=1;
+ classes(k).cluster = mixture.cluster(k);
+ classes(k).Rmin = NaN;
+ classes(k).rissanen = NaN;
+ classes(k).loglikelihood = NaN;
+end
diff --git a/GaussMix/TestGMMFCC.mat b/GaussMix/TestGMMFCC.mat
new file mode 100755
index 0000000000000000000000000000000000000000..c86e6cacdb3c7bcd11d0fce55faa989de4fddd3b
Binary files /dev/null and b/GaussMix/TestGMMFCC.mat differ
diff --git a/GaussMix/example1/CreerData.asv b/GaussMix/example1/CreerData.asv
new file mode 100755
index 0000000000000000000000000000000000000000..4de36da90c6abe6af65d8bd72bc0b113805301fe
--- /dev/null
+++ b/GaussMix/example1/CreerData.asv
@@ -0,0 +1,28 @@
+
+
+x = zeros(1,100);
+for i=1:2000
+ u = rand(1,1);
+ if (u < 0.3)
+ x(i) = 7+randn(1,1);
+ else if (u<0.8)
+ x(i) = randn(1,1);
+ else
+ x(i) = 4+0.5*randn(1,1);
+ end;
+end;
+end;
+
+% I= hist(x,100);
+% figure;
+% plot(I);
+x=x'
+
+[mtrs,omtr] = GaussianMixture(x, 20, 0, 0);
+disp(sprintf('\toptimal order K*: %d', omtr.K));
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
\ No newline at end of file
diff --git a/GaussMix/example1/CreerData.m b/GaussMix/example1/CreerData.m
new file mode 100755
index 0000000000000000000000000000000000000000..4de36da90c6abe6af65d8bd72bc0b113805301fe
--- /dev/null
+++ b/GaussMix/example1/CreerData.m
@@ -0,0 +1,28 @@
+
+
+x = zeros(1,100);
+for i=1:2000
+ u = rand(1,1);
+ if (u < 0.3)
+ x(i) = 7+randn(1,1);
+ else if (u<0.8)
+ x(i) = randn(1,1);
+ else
+ x(i) = 4+0.5*randn(1,1);
+ end;
+end;
+end;
+
+% I= hist(x,100);
+% figure;
+% plot(I);
+x=x'
+
+[mtrs,omtr] = GaussianMixture(x, 20, 0, 0);
+disp(sprintf('\toptimal order K*: %d', omtr.K));
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
\ No newline at end of file
diff --git a/GaussMix/example1/GMM.m b/GaussMix/example1/GMM.m
new file mode 100755
index 0000000000000000000000000000000000000000..4fcd1dbcd716f12f4cd4b37d9b665dd4de4e5ded
--- /dev/null
+++ b/GaussMix/example1/GMM.m
@@ -0,0 +1,28 @@
+function x = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3)
+
+pi1 = Pi(1);
+pi2 = Pi(2);
+pi3 = Pi(3);
+
+[V,D] = eig(R1);
+A1 = V*sqrt(D);
+
+[V,D] = eig(R2);
+A2 = V*sqrt(D);
+
+[V,D] = eig(R3);
+A3 = V*sqrt(D);
+
+x1 = A1*randn(2,N) + mu1*ones(1,N);
+
+x2 = A2*randn(2,N) + mu2*ones(1,N);
+
+x3 = A3*randn(2,N) + mu3*ones(1,N);
+
+SwitchVar = ones(2,1)*random('Uniform',0,1,1,N);
+SwitchVar1 = SwitchVar<pi1;
+SwitchVar2 = (SwitchVar>=pi1)&(SwitchVar<(pi1+pi2));
+SwitchVar3 = SwitchVar>=(pi1+pi2);
+
+x = SwitchVar1.*x1 + SwitchVar2.*x2 + SwitchVar3.*x3;
+
diff --git a/GaussMix/example1/data b/GaussMix/example1/data
new file mode 100755
index 0000000000000000000000000000000000000000..763362b59501ec090388560980b507b65430b842
--- /dev/null
+++ b/GaussMix/example1/data
@@ -0,0 +1,500 @@
+ -2.4994293e+000 -4.1082971e+000
+ -2.9671425e+000 -1.3371671e+000
+ 3.5226619e+000 2.5098193e+000
+ 3.0722951e+000 7.7605486e-001
+ -4.5001505e-001 -9.9138714e-001
+ 1.1487160e+000 1.5684598e+000
+ 1.8978815e+000 1.8575634e+000
+ 2.3897087e+000 3.2318664e+000
+ 3.1647261e+000 1.2402088e+000
+ 1.0476926e+000 1.1514808e+000
+ 4.1264519e+000 9.4973383e-001
+ 4.0108082e-001 2.4573201e+000
+ 1.5835736e+000 2.0818306e+000
+ -2.4758251e+000 -2.3757753e+000
+ -2.0908266e+000 -1.6581249e+000
+ 2.4069064e+000 3.8835300e+000
+ 8.3616266e-001 9.5145138e-001
+ 4.7415358e+000 2.8080942e+000
+ -3.5726626e+000 -1.4829159e+000
+ 1.3554735e+000 2.3591238e+000
+ 6.9199495e+000 3.4346874e+000
+ 6.0331966e+000 2.4398066e+000
+ 1.5409646e+000 2.1959455e+000
+ 3.1841053e+000 2.9210719e+000
+ -3.5681563e+000 -1.5778763e+000
+ 5.2829597e+000 1.0719071e+000
+ 2.1520081e+000 3.0866278e+000
+ 2.3234343e+000 1.9964561e+000
+ -2.4439550e+000 -2.6168582e+000
+ -1.6585907e+000 -3.0834572e+000
+ -1.4685333e+000 -3.3949080e-001
+ 1.8131550e+000 2.0644860e+000
+ 1.9260310e-002 -1.9636122e+000
+ 6.0210103e+000 2.1427658e+000
+ 5.2773175e+000 2.8105439e+000
+ 2.5903782e+000 1.8598159e+000
+ 2.7660793e+000 1.9606567e+000
+ 2.7735901e+000 3.7655498e+000
+ -6.8620414e-001 -3.1134468e+000
+ -2.6259834e+000 -1.3972865e+000
+ 3.4318861e+000 1.9934388e+000
+ -1.4158958e+000 -1.8638821e+000
+ -3.8950601e+000 -2.5381157e+000
+ -2.9327580e+000 -9.9066019e-001
+ 2.4690993e+000 3.5762467e+000
+ 2.3623143e+000 9.4261649e-001
+ -2.8539027e+000 -2.1259995e+000
+ 4.8183419e+000 1.7034950e+000
+ 3.1139780e+000 3.6117760e-001
+ 1.3829018e+000 2.7794457e-001
+ 5.9234311e+000 2.8948173e+000
+ -1.8658257e+000 -3.1983343e+000
+ 5.1705721e+000 1.6518858e+000
+ 4.9920924e+000 1.7147543e+000
+ 8.7775965e-001 5.2625589e-001
+ 1.5750626e+000 1.1918216e+000
+ 1.7992611e+000 4.8551802e-001
+ -3.2489881e+000 -2.5787498e+000
+ -1.6014959e+000 -4.6100429e+000
+ -1.8736466e+000 -1.5236121e+000
+ 3.2909016e+000 2.2200528e+000
+ 3.1761381e+000 1.6104832e+000
+ -2.4872971e+000 -1.9009485e+000
+ 2.1772792e+000 1.3870138e+000
+ 1.7380968e+000 6.0048440e-001
+ -3.9767301e+000 -4.6853484e-001
+ -2.1566702e+000 -1.1963494e+000
+ -2.2995271e+000 -3.1264653e-001
+ 4.8126105e-001 1.5430571e+000
+ 2.4882431e+000 -6.3760172e-001
+ 2.0665097e+000 2.5139615e+000
+ 4.3522767e+000 2.1681717e+000
+ 1.5236531e+000 1.7700979e+000
+ 4.1187868e-001 1.5684020e+000
+ 1.0653372e+000 1.6758234e+000
+ 3.3871896e+000 2.9157224e-001
+ -1.9947732e+000 -1.6310953e+000
+ -1.3005316e+000 -4.6750525e+000
+ 6.0101907e+000 2.6049915e+000
+ -1.3409899e+000 -1.7805455e+000
+ -3.2479310e-001 -1.6668795e+000
+ 2.5782965e+000 3.6865760e+000
+ -1.1058942e+000 -1.7315941e+000
+ 5.8395639e+000 8.6635976e-001
+ -2.6851602e+000 1.1284543e-002
+ 5.4012564e+000 3.7084575e+000
+ 6.4876500e+000 6.5784379e-001
+ 6.5756377e+000 2.7666247e+000
+ 1.7034262e+000 1.5074734e+000
+ 2.2215234e-001 2.4789311e+000
+ -4.3923475e-001 -1.1053623e+000
+ 2.6451169e+000 1.0697458e+000
+ -3.6583706e+000 -3.4110735e+000
+ 4.3954146e+000 2.7009695e+000
+ 1.8400307e+000 3.9912944e+000
+ 1.0370402e+000 7.2896396e-001
+ -2.7745089e+000 -1.7416299e+000
+ 6.8834332e+000 3.7827272e+000
+ 5.2886275e+000 1.3362601e+000
+ -2.9592467e+000 -2.9498553e+000
+ 2.1205415e+000 1.1262433e+000
+ 6.6310496e+000 3.1048442e+000
+ -2.7165406e+000 -9.4433044e-001
+ -2.8074161e+000 -2.0057037e+000
+ 4.4858394e+000 2.2253963e+000
+ -1.8775469e+000 -2.3935334e+000
+ -3.8600436e+000 -1.8707910e+000
+ -2.2594070e+000 -9.6137113e-001
+ 6.3681611e+000 2.8750472e+000
+ 2.1929836e+000 2.1516427e+000
+ -2.7938808e+000 -1.4933715e+000
+ 5.2219233e+000 2.0053517e+000
+ -2.2673249e+000 -2.2135980e+000
+ -3.0688212e+000 -4.4197471e-001
+ 2.6787741e+000 4.4606319e-001
+ 2.5000255e+000 2.1269632e+000
+ 3.4701781e+000 1.2293687e+000
+ 2.8626905e+000 1.2354852e+000
+ 1.8049405e+000 2.6801678e+000
+ 2.0198741e+000 3.4723794e+000
+ 5.0153195e+000 2.0015422e+000
+ -2.6667171e+000 -1.2963055e+000
+ -2.2747819e+000 -1.2941752e+000
+ -3.4010743e+000 -3.1685622e+000
+ 1.8608573e+000 2.6585750e+000
+ -3.4723835e+000 -1.4888637e+000
+ 1.3271764e+000 1.6488526e+000
+ 3.3211619e+000 2.4470109e+000
+ -3.0254025e+000 -3.1373110e+000
+ 3.9269207e+000 1.4379519e+000
+ 1.4228050e+000 2.5789187e+000
+ -1.6472955e+000 -1.0035893e+000
+ 6.2462816e+000 2.5207719e+000
+ 3.6151010e+000 5.0345990e-001
+ 2.9414300e+000 1.6689092e+000
+ -4.0542399e+000 -1.4068429e+000
+ 6.6200627e+000 2.0605262e+000
+ -3.2322167e+000 -1.4475055e+000
+ 2.3062268e+000 3.0255570e+000
+ -1.5044377e+000 -1.8478036e+000
+ -3.6595067e+000 -1.1634359e+000
+ -1.4868748e+000 -3.2757052e+000
+ 2.7377191e+000 2.2244277e+000
+ -1.6737283e+000 -1.9988376e+000
+ -1.1719197e+000 -1.6477494e+000
+ 2.4196225e+000 2.0280467e+000
+ 1.0420402e+000 1.5786039e+000
+ 5.3741767e-001 3.0017498e-001
+ 7.6512181e-001 2.3017858e+000
+ 2.6192367e+000 1.0041025e+000
+ 5.9375633e+000 2.1165149e+000
+ 3.9911894e+000 3.2389784e+000
+ -6.1544639e-001 9.1196521e-001
+ 6.6706209e-001 1.2591643e+000
+ -1.5470539e+000 -1.4867496e+000
+ 1.2959077e+000 2.2990982e+000
+ 5.3865407e+000 1.9372642e+000
+ 2.2732782e+000 2.7126818e+000
+ 2.5961858e+000 2.7104783e+000
+ 1.2016437e+000 1.6352796e+000
+ -6.3132877e-001 -1.0130518e+000
+ 2.4821847e+000 1.6711225e+000
+ 6.4111731e+000 2.8796088e+000
+ 2.2600424e+000 1.0462248e+000
+ 3.1806037e+000 3.0130462e+000
+ -2.8974380e+000 -2.3120396e+000
+ 1.6195341e+000 2.0335954e+000
+ -2.8721153e+000 -2.6440364e+000
+ -1.3574849e+000 -1.6247100e+000
+ 6.2521103e+000 1.0760627e+000
+ -2.6826359e+000 -1.8705329e+000
+ -1.8305062e+000 -2.3678353e+000
+ -2.2017280e+000 -1.8803017e+000
+ 3.1030747e+000 4.4929823e+000
+ -2.5405323e+000 -2.7039526e+000
+ 9.4416429e-001 1.5192220e+000
+ 4.4038476e+000 1.5023491e+000
+ -1.9942109e+000 -1.6365715e+000
+ 6.1075054e+000 2.2942330e+000
+ -3.6092914e+000 -1.0057776e+000
+ -3.2368954e+000 -2.0128304e+000
+ 8.8089209e-001 5.9928526e-001
+ -2.0880661e+000 -2.2373511e+000
+ 4.3965292e+000 2.0705245e+000
+ -2.1261464e+000 -3.2506312e+000
+ 6.7842371e+000 2.9638172e+000
+ -2.3015041e+000 -2.3005755e+000
+ -1.7925546e+000 -1.1928439e+000
+ -2.1050996e+000 -8.2983502e-001
+ -8.2905500e-001 -3.2606595e+000
+ 5.0774670e+000 2.2600273e+000
+ 1.2339099e+000 2.3321118e+000
+ 4.3562376e+000 3.5868654e+000
+ 2.3699913e+000 8.9875699e-001
+ -1.4283935e+000 -1.3723794e+000
+ -1.8670244e+000 -1.6451963e+000
+ 3.2719773e+000 3.4186424e+000
+ -3.8567427e+000 -8.7553847e-001
+ -2.9588459e+000 -2.9418173e+000
+ -8.0457590e-002 -2.6596755e+000
+ -2.1739172e+000 -1.8488517e+000
+ 4.8115445e+000 1.9223987e+000
+ -2.5861401e+000 -2.4801926e+000
+ 6.0791681e+000 2.1444257e+000
+ 6.0810213e+000 1.8565889e+000
+ -3.7111499e+000 -2.8075019e+000
+ -2.1317940e+000 -1.7202581e+000
+ -2.1871322e+000 -1.9296521e+000
+ 6.8244703e+000 3.7830800e+000
+ 2.6280031e+000 -8.4739518e-002
+ 5.7508871e+000 2.3361165e+000
+ 2.2859335e+000 7.8499700e-001
+ 1.0961758e+000 1.4309390e+000
+ -2.2225730e+000 -2.4154014e+000
+ 5.1923645e+000 2.8447880e+000
+ 5.9455299e+000 2.0144543e+000
+ 2.0454603e+000 1.0746073e+000
+ 5.4967947e+000 2.1467572e+000
+ -9.0737347e-001 -2.6879145e+000
+ -2.2394370e+000 -1.7345858e+000
+ 5.0307396e+000 1.7456615e+000
+ 5.9399277e+000 2.1752479e+000
+ 1.0030725e+000 2.0221404e+000
+ -6.8558120e-001 -1.5381992e+000
+ 2.1085190e+000 1.6959598e+000
+ 2.3769739e+000 3.1757300e+000
+ 2.3758073e+000 2.1102179e+000
+ 4.5385560e+000 1.8018332e+000
+ 3.7608027e+000 1.1473318e+000
+ 3.9843891e-002 -2.8071851e+000
+ -1.0669756e+000 -2.8841307e+000
+ 3.3513297e+000 1.0460583e+000
+ -2.1429929e+000 -1.5439410e+000
+ -1.5504817e+000 -4.0493829e+000
+ 1.5335221e+000 2.7389541e+000
+ -2.1727474e+000 -1.9386951e+000
+ 1.7802660e+000 2.3705238e+000
+ -8.3922759e-001 -2.0204560e+000
+ 2.3979499e+000 -8.9853418e-002
+ -1.3662227e+000 -2.1153751e+000
+ -3.0690249e+000 -1.4784637e+000
+ 4.3653690e+000 2.7688269e+000
+ 2.8148501e+000 3.8227296e+000
+ 2.2785650e-001 1.8687671e+000
+ -1.1600826e+000 -9.5060212e-001
+ -2.4981155e+000 -3.2322247e+000
+ 3.6279117e+000 3.6755124e+000
+ -1.6643124e+000 -1.9617389e+000
+ 4.5297138e+000 1.3667072e+000
+ 1.5016769e+000 1.5956524e+000
+ -1.1310457e+000 -2.3742534e+000
+ 2.1353357e+000 3.3368753e-001
+ -5.1364193e+000 -3.4352652e+000
+ -1.5802164e+000 -3.4760005e+000
+ -1.1484795e+000 -2.0672066e+000
+ 2.9825913e+000 2.5196367e+000
+ -1.8204054e+000 -2.3251292e+000
+ -1.3828687e+000 -1.3749753e+000
+ 5.1681024e+000 1.5788096e+000
+ -2.6090401e+000 -1.6338400e+000
+ 3.0478763e+000 2.6484737e+000
+ 3.0778957e+000 2.3646524e+000
+ 2.4812899e+000 1.7806173e+000
+ 5.0136691e+000 7.4088631e-001
+ 3.7797550e+000 3.1917534e+000
+ 1.0416715e+000 2.3175632e+000
+ 1.2976376e+000 2.1682829e+000
+ -4.0733788e+000 -1.4998157e+000
+ -2.6353104e+000 -8.6600700e-001
+ -1.5903014e+000 -3.6891517e+000
+ 9.4381385e-001 3.2866631e+000
+ 1.2031462e+000 -4.5662401e-001
+ 6.8770249e+000 2.6427022e+000
+ 2.0938181e+000 2.5196653e+000
+ 3.7776215e-001 -1.3551820e+000
+ 2.5355526e+000 3.3705178e-001
+ 2.9347044e+000 1.8361782e+000
+ 1.1525353e+000 1.1335376e+000
+ 5.5398348e+000 2.2035110e+000
+ 5.9360527e+000 3.1458894e+000
+ 6.3601136e+000 2.0167516e+000
+ 1.5406977e+000 2.7266677e+000
+ 2.1369754e+000 2.8760243e+000
+ 1.5676395e+000 3.1381758e+000
+ 1.0694124e+000 2.9641214e+000
+ 5.9695699e+000 2.7932636e+000
+ -1.7399229e+000 -2.6515495e+000
+ -2.2964794e+000 -1.6852163e+000
+ 5.7277030e+000 2.2633989e+000
+ 1.2276520e+000 1.6188346e+000
+ 5.3536175e+000 2.6043457e+000
+ 2.4355557e+000 2.7635529e+000
+ -3.4575774e+000 -2.6766887e+000
+ -1.8857457e+000 -2.4695457e+000
+ 3.4663421e-001 -1.4268079e+000
+ 7.3955070e-001 1.3103992e+000
+ -2.7262870e+000 -2.2668589e+000
+ -1.7271116e+000 -1.5457363e+000
+ -2.4571638e+000 -2.8338694e+000
+ 5.4956320e+000 2.3725309e+000
+ -1.7864949e+000 -1.7087300e+000
+ 5.5460691e+000 2.4587266e+000
+ -2.1120357e+000 -1.9923365e+000
+ 4.9024754e+000 1.4228082e+000
+ 1.5072825e+000 -6.6505503e-002
+ -2.2308727e+000 -1.7066171e+000
+ 2.5097672e+000 2.6576413e+000
+ -1.7711854e+000 -1.6559841e+000
+ -3.7926587e+000 -2.6101264e+000
+ 2.7419614e+000 2.6831726e+000
+ -5.0561496e+000 -1.0016540e+000
+ 3.2150079e+000 1.0317253e+000
+ 7.5558296e+000 1.0081801e+000
+ 4.2719990e+000 2.7305714e+000
+ 3.2046111e+000 3.8978119e-002
+ 5.4592655e+000 2.0958170e+000
+ 1.8891579e+000 2.3156003e+000
+ 3.0530168e+000 2.0994657e+000
+ 2.0625960e+000 2.8976251e+000
+ 3.0631707e+000 2.4921259e+000
+ -2.5854220e+000 -2.0930551e+000
+ 5.1121862e-001 -2.3170825e+000
+ 5.4756129e+000 6.3495339e-001
+ -9.6832621e-001 -3.2581462e+000
+ -1.6020614e+000 -3.1597720e-001
+ -1.6160363e+000 -2.1889253e+000
+ -2.2922629e+000 -3.4736102e+000
+ 4.7355707e+000 1.1807202e+000
+ -1.9163866e+000 -3.1835974e+000
+ 9.5348377e-001 2.7303333e+000
+ -1.3069116e+000 -6.6321671e-001
+ 3.1595416e+000 1.3010167e+000
+ 2.5409489e+000 2.2062242e+000
+ -2.9646523e+000 -4.1420861e+000
+ -1.9482757e+000 -2.7774977e+000
+ -3.7298518e+000 -2.9219899e+000
+ 1.9492594e+000 2.8750611e+000
+ 1.6148143e+000 1.1183454e-001
+ 1.6865065e+000 2.1356231e-001
+ -3.3230196e+000 -1.0076272e+000
+ -1.0175537e+000 -1.3770098e+000
+ 6.1428186e+000 3.2230843e+000
+ 5.7424349e+000 1.6640713e+000
+ 2.8938332e+000 1.7212667e+000
+ 2.0252725e+000 2.4927555e+000
+ -3.1194636e+000 -2.1353075e+000
+ -2.9663900e+000 -5.6372352e-001
+ 3.5768415e+000 1.9624878e+000
+ -2.3426612e+000 -1.0691368e+000
+ 2.5904333e+000 9.6885642e-001
+ 5.2274555e+000 1.8349025e+000
+ -5.7419221e-001 -1.1911130e+000
+ -2.0096914e+000 -1.8773094e+000
+ 2.7161994e-001 2.8803033e+000
+ 1.0038912e+000 6.6843775e-001
+ 1.3279680e+000 6.8749193e-001
+ 4.5529519e+000 2.3219142e+000
+ -1.2486905e+000 -2.6455937e+000
+ 4.1765762e+000 1.9889733e+000
+ 5.2985362e+000 1.8318100e+000
+ 6.6957657e-001 -2.6358386e+000
+ -2.7606705e+000 -1.0360419e+000
+ 6.7091388e+000 1.4413358e+000
+ 3.6246178e+000 9.9935693e-001
+ -2.7923586e+000 -8.7301821e-001
+ -2.0477387e+000 -2.4379086e+000
+ -3.4545172e+000 -1.9973320e+000
+ 1.5851379e+000 2.6718295e+000
+ -3.1462677e+000 -3.2435903e+000
+ 1.3037452e+000 2.1717478e+000
+ -2.8348318e+000 -1.6370994e+000
+ -2.3588869e+000 -1.2982025e+000
+ -2.8350940e+000 -2.5717469e+000
+ -1.1726033e+000 -3.0330671e+000
+ 1.8079635e+000 8.0323989e-001
+ -1.0502815e+000 -1.2803310e+000
+ 6.6860937e-001 4.4105544e+000
+ 3.1390514e+000 2.1018770e+000
+ 1.7004917e+000 1.8667772e-001
+ 9.4181489e-001 1.9043648e+000
+ 2.8311446e+000 3.2550841e+000
+ -1.2120929e+000 -2.6976946e+000
+ -6.2490996e-001 -2.6295920e+000
+ -1.9089542e+000 -2.5796865e+000
+ 2.8829971e+000 8.1569343e-001
+ 6.3056091e+000 3.4716864e+000
+ 5.1860128e+000 3.1216642e+000
+ 1.7995995e+000 1.3115364e+000
+ -1.6427606e+000 -1.8998332e+000
+ -3.0541267e+000 -1.7727272e+000
+ 5.3583700e+000 2.5195458e+000
+ 6.4533947e+000 9.7055473e-001
+ 2.7491020e+000 1.4902388e+000
+ 1.8759329e+000 1.5161278e+000
+ -2.1471341e+000 -2.4697487e+000
+ 2.9430094e+000 4.6616711e-001
+ 1.6883411e+000 9.5909480e-001
+ 1.0176063e+000 1.8932934e+000
+ 4.9954489e+000 2.3508231e+000
+ -1.7075862e+000 -2.3608536e+000
+ -1.2525656e+000 -1.6090932e+000
+ 6.4682087e+000 1.8036302e+000
+ -2.5343885e+000 -1.3739124e+000
+ 4.2535920e+000 1.7999511e+000
+ 5.9235653e+000 2.2666082e+000
+ -1.3119696e+000 -6.6952007e-001
+ -2.9715463e+000 -2.1128425e+000
+ 2.0459604e+000 1.6613328e+000
+ 6.7529465e-001 1.2194714e+000
+ -4.6953973e-001 -2.7773072e+000
+ 6.4871660e-001 1.4508362e+000
+ 3.0543470e+000 9.8357805e-002
+ 2.5553322e+000 2.8424386e+000
+ 1.1047984e+000 1.2305443e+000
+ 9.2519628e-001 1.3359501e+000
+ -1.3577254e+000 -1.5453333e+000
+ 5.3610518e+000 8.2950077e-001
+ 4.8887129e+000 2.0105734e+000
+ 5.4611479e+000 1.5400073e+000
+ -2.0199150e+000 -2.8027347e+000
+ 4.6353981e+000 4.1226069e+000
+ -1.7243908e+000 -4.3272905e-001
+ 5.2500609e+000 2.6140689e+000
+ 5.1561998e+000 5.4736171e-001
+ 1.8033738e+000 4.9629305e-001
+ -1.5068152e+000 -1.8372543e+000
+ 2.8026475e+000 5.3078854e-001
+ 5.7924020e-001 4.6615051e-001
+ 1.1281560e+000 2.2568055e+000
+ 4.2737184e-001 2.9928316e+000
+ -3.0972141e+000 -1.6276929e+000
+ 5.1282222e+000 2.2819255e+000
+ 2.7021225e+000 3.2288175e+000
+ -2.6313851e+000 -4.7662051e-001
+ -2.9961837e+000 -2.8550736e+000
+ -2.0922997e+000 -1.3346492e+000
+ -1.4670433e+000 -1.6524732e+000
+ 4.3603392e+000 1.1703490e+000
+ 7.7870829e+000 2.0911527e+000
+ -2.0887821e+000 -3.4455226e+000
+ 5.0453742e-001 2.3007790e-001
+ -2.4996947e+000 -2.0166070e+000
+ 2.6488793e+000 2.1306912e+000
+ -2.0636229e+000 -2.8127678e+000
+ -1.5921498e-001 -3.2021723e+000
+ -5.4887003e-001 -2.8496663e+000
+ 2.0638411e+000 3.0905333e+000
+ -1.8121313e+000 -2.3734007e+000
+ 6.2095394e+000 2.8602269e+000
+ 3.4156325e+000 2.4850765e+000
+ 1.6259395e+000 5.3956246e-001
+ 5.4054894e+000 1.2002490e+000
+ -1.4681865e+000 -3.7393502e+000
+ -1.6975912e+000 -2.9794286e+000
+ 5.9329195e+000 3.5411223e+000
+ -2.0786082e+000 -2.3214766e+000
+ 5.3655241e+000 2.0659076e+000
+ 1.6391686e+000 1.6924856e+000
+ 1.7408993e+000 1.8874708e+000
+ 1.3613607e+000 1.1065809e+000
+ -7.7581839e-001 -4.2582739e-001
+ -2.2932760e+000 -6.6236212e-001
+ 6.2569401e+000 1.8894429e+000
+ -8.6559016e-001 -1.8765065e+000
+ -1.6558140e+000 -1.9951723e+000
+ 6.6799991e+000 2.2719981e+000
+ 1.7305565e+000 2.9160533e+000
+ -2.6050169e+000 -1.0229478e+000
+ 5.1558983e+000 1.8501076e+000
+ -1.9716878e+000 -2.5266854e+000
+ 1.7903652e+000 3.3824354e+000
+ -3.0046893e+000 -2.3820382e+000
+ 4.3701908e+000 4.8607177e+000
+ 6.6515634e+000 2.1910670e+000
+ -2.8303812e+000 -1.4199024e+000
+ 4.3756999e+000 3.5064352e+000
+ 6.6400526e+000 1.4397338e+000
+ -4.0009855e+000 -2.5753713e+000
+ 2.8726427e+000 5.0006628e-001
+ 3.1868307e+000 1.7739456e+000
+ -1.7014721e+000 -4.1307275e+000
+ 1.1083038e+000 3.1349628e+000
+ 3.0441527e+000 2.6647854e+000
+ -2.4206111e+000 -1.9813182e+000
+ 2.1519674e+000 1.2542581e+000
+ 2.8775567e+000 1.9762391e+000
+ -2.5814464e+000 -2.4029658e+000
+ 5.5764723e+000 2.1029966e+000
+ 1.4964719e+000 1.9815920e+000
+ 2.1989404e+000 1.9957872e+000
+ 6.1761430e+000 1.1944073e+000
+ 1.9783625e+000 3.3727155e+000
+ -3.7122971e+000 -4.0566829e+000
+ -1.8479871e+000 -2.5371069e+000
+ 7.6459608e+000 1.3928904e+000
+ -3.1317781e+000 -7.3167208e-001
+ -1.0045745e+000 -1.8587473e+000
+ -1.8508324e+000 -2.0473195e+000
+ 3.9335433e+000 1.9983886e+000
+ -1.7877225e-001 1.7128225e+000
diff --git a/GaussMix/example1/mkdata.m b/GaussMix/example1/mkdata.m
new file mode 100755
index 0000000000000000000000000000000000000000..eb40c08edc38cb3cbe68f46763c32babfef1b4d3
--- /dev/null
+++ b/GaussMix/example1/mkdata.m
@@ -0,0 +1,28 @@
+clear all;
+close all;
+clc;
+
+N=500;
+
+R1 = [ 1, 0.1; 0.1, 1];
+mu1 = [2, 2]';
+
+R2 = [ 1, -0.1; -0.1, 1];
+mu2 = [-2, -2]';
+
+R3 = [ 1, 0.2; 0.2, 0.5];
+mu3 = [5.5, 2]';
+
+Pi = [0.4, 0.4, 0.2];
+
+% Generate Data from Gaussian mixture model
+x = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3);
+
+plot(x(1,:),x(2,:),'o');
+title('Scatter Plot of Multimodal Data')
+xlabel('first component')
+ylabel('second component')
+
+x = x';
+save data x /ascii
+
diff --git a/GaussMix/example1/rundemo.m b/GaussMix/example1/rundemo.m
new file mode 100755
index 0000000000000000000000000000000000000000..77bb2320de87c51061b57e0429ccea62826248c1
--- /dev/null
+++ b/GaussMix/example1/rundemo.m
@@ -0,0 +1,48 @@
+clear all;
+
+if exist('GaussianMixture')~=2
+ pathtool;
+ error('the directory containing the Cluster program must be added to the search path');
+end
+
+disp('generating data...');
+mkdata;
+clear all;
+pixels = load('data');
+input('press <Enter> to continue...');
+disp(' ');
+
+% [mtrs, omtr] = GaussianMixture(pixels, initK, finalK, verbose)
+% - pixels is a NxM matrix, containing N training vectors, each of M-dimensional
+% - start with initK=20 initial clusters
+% - finalK=0 means estimate the optimal order
+% - verbose=true displays clustering information
+% - mtrs is an array of structures, each containing the cluster parameters of the
+% mixture of a particular order
+% - omtr is a structure containing the cluster parameters of the mixture with
+% the estimated optimal order
+disp('estimating optimal order and clustering data...');
+[mtrs,omtr] = GaussianMixture(pixels, 20, 0, true);
+disp(sprintf('\toptimal order K*: %d', omtr.K));
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
+input('press <Enter> to continue...');
+disp(' ');
+
+% - finalK=5 means to cluster assuming mixture order of 5
+% - omtr contains the cluster parameters of the mixture of order 5
+disp('clustering data assuming optimal order of 5...');
+[mtrs,omtr] = GaussianMixture(pixels, 20, 5, true);
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
+input('press <Enter> to continue...');
+
+
diff --git a/GaussMix/example2/GMM.m b/GaussMix/example2/GMM.m
new file mode 100755
index 0000000000000000000000000000000000000000..4fcd1dbcd716f12f4cd4b37d9b665dd4de4e5ded
--- /dev/null
+++ b/GaussMix/example2/GMM.m
@@ -0,0 +1,28 @@
+function x = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3)
+
+pi1 = Pi(1);
+pi2 = Pi(2);
+pi3 = Pi(3);
+
+[V,D] = eig(R1);
+A1 = V*sqrt(D);
+
+[V,D] = eig(R2);
+A2 = V*sqrt(D);
+
+[V,D] = eig(R3);
+A3 = V*sqrt(D);
+
+x1 = A1*randn(2,N) + mu1*ones(1,N);
+
+x2 = A2*randn(2,N) + mu2*ones(1,N);
+
+x3 = A3*randn(2,N) + mu3*ones(1,N);
+
+SwitchVar = ones(2,1)*random('Uniform',0,1,1,N);
+SwitchVar1 = SwitchVar<pi1;
+SwitchVar2 = (SwitchVar>=pi1)&(SwitchVar<(pi1+pi2));
+SwitchVar3 = SwitchVar>=(pi1+pi2);
+
+x = SwitchVar1.*x1 + SwitchVar2.*x2 + SwitchVar3.*x3;
+
diff --git a/GaussMix/example2/TestingData b/GaussMix/example2/TestingData
new file mode 100755
index 0000000000000000000000000000000000000000..caf2dfe08cd65d0cb96999590e4d808416d61aa0
--- /dev/null
+++ b/GaussMix/example2/TestingData
@@ -0,0 +1,50 @@
+ 2.1971862e+000 1.1993495e+000
+ -8.3043360e-001 -2.7987701e+000
+ 5.8665530e+000 3.1898402e+000
+ -3.4132031e+000 -7.6621723e-001
+ 2.8730733e+000 5.0462144e-001
+ 1.6707453e+000 1.2947682e+000
+ -1.2146352e+000 -3.1241962e+000
+ -5.1492086e+000 -1.8269926e+000
+ 2.6036572e+000 2.3810209e+000
+ -1.2640568e-001 -1.6500479e+000
+ -2.7104621e+000 -1.5467451e+000
+ 1.3685586e+000 2.3863738e+000
+ -2.5655313e+000 -8.3511554e-001
+ -2.0774062e+000 3.9271006e-001
+ 2.3308628e+000 3.1675733e+000
+ 4.5698037e+000 1.7263580e+000
+ -2.3594068e+000 -1.6137032e+000
+ 2.2523826e+000 2.3207585e+000
+ -8.6639573e-002 -2.1695638e+000
+ 2.5863777e+000 1.4053144e+000
+ -1.8543777e+000 -1.5331829e+000
+ 5.1418028e+000 2.7684652e+000
+ 5.7766524e+000 1.6944210e+000
+ 6.3652491e+000 1.6645575e+000
+ 5.7823120e+000 7.9581296e-001
+ -2.7636660e+000 1.6935928e+000
+ -4.3570958e-001 4.4371331e+000
+ 1.7101286e+000 -2.5487073e+000
+ 2.8951547e+000 -3.1986172e+000
+ -2.2143221e+000 1.6318553e+000
+ -4.0074078e+000 2.5685027e+000
+ 1.7145779e+000 -2.5571522e+000
+ 1.7178859e+000 -2.2144649e+000
+ -6.4713146e+000 2.6336545e+000
+ -7.9566172e-001 1.1551709e+000
+ 1.5224425e+000 -8.2340341e-001
+ 4.1321043e+000 -1.2519017e+000
+ -5.0331120e+000 3.5826263e+000
+ 4.6795003e+000 -2.0810188e+000
+ 2.1414433e+000 -2.6374452e+000
+ 1.2357672e+000 -2.5872921e+000
+ -5.4049599e+000 1.3408878e+000
+ -2.4818577e+000 2.5167067e+000
+ -3.3851656e+000 3.3410145e+000
+ -1.0861205e+000 2.4655977e+000
+ 2.7245175e+000 -2.4150836e+000
+ 3.2744279e+000 -3.1687109e+000
+ -5.6289618e+000 1.2555843e+000
+ 1.4646357e+000 -2.8137472e+000
+ -1.6799700e+000 1.7721904e+000
diff --git a/GaussMix/example2/TrainingData1 b/GaussMix/example2/TrainingData1
new file mode 100755
index 0000000000000000000000000000000000000000..52da52cee185f4e3c0e1dabff434702402072d17
--- /dev/null
+++ b/GaussMix/example2/TrainingData1
@@ -0,0 +1,500 @@
+ -8.2691760e-001 -1.2050459e+000
+ -3.9239203e+000 -2.0161882e+000
+ 2.4879670e-001 1.3844788e+000
+ 4.4622537e-003 -2.8884597e+000
+ 5.5956608e+000 1.7152329e+000
+ 9.7420878e-001 7.2623676e-001
+ 4.1758622e+000 2.3983047e+000
+ 3.8345761e+000 1.6565882e+000
+ 5.4856455e+000 9.0859513e-001
+ 1.4821777e+000 2.2860795e+000
+ 2.8648435e+000 2.6921048e+000
+ -3.4392392e+000 -2.2906714e+000
+ 7.9879281e-001 2.5323459e+000
+ -2.9776505e+000 -2.5009997e+000
+ 3.4136143e+000 3.7066177e+000
+ 6.1058968e+000 2.0012507e+000
+ -3.0344391e+000 -1.2340273e+000
+ 4.9556701e+000 1.6499631e+000
+ 7.3941501e+000 3.1326202e+000
+ -3.0191433e+000 -1.2086618e+000
+ -1.5729548e+000 -2.9380380e+000
+ 4.3745711e+000 1.5469558e+000
+ 5.7839054e+000 3.0041349e+000
+ 5.6721776e+000 1.3579906e+000
+ -1.6443709e+000 -2.5504814e+000
+ 9.9005870e-001 4.0772085e-001
+ -3.0709658e+000 3.6791992e-001
+ -2.0692758e+000 -3.9017494e+000
+ -3.6390477e+000 -2.6026315e+000
+ -2.1984618e+000 -3.4126216e+000
+ 1.4358937e+000 3.1557529e+000
+ -1.8685433e+000 -1.7996294e+000
+ -3.0881358e+000 -3.3513457e+000
+ -1.8835385e+000 -3.1859952e+000
+ -2.2010907e+000 -1.2703709e+000
+ 2.6059290e+000 2.5827685e+000
+ -1.9235420e+000 -3.3851239e+000
+ 1.7448817e+000 6.2499439e-001
+ -1.5948131e+000 -2.7398277e+000
+ 3.0825595e+000 2.9173824e+000
+ 3.7672093e+000 5.9095598e-001
+ -2.6432716e+000 -2.7013228e+000
+ -9.4938042e-001 -2.5131612e+000
+ 4.9853493e+000 2.6886787e+000
+ 1.7663306e+000 1.6747798e+000
+ 3.1396720e+000 2.5144983e+000
+ 5.2599451e+000 1.8719227e+000
+ -6.7534215e-001 -4.1068096e+000
+ -1.7491434e+000 -4.0765221e+000
+ 6.5244821e+000 2.2127293e+000
+ 1.8687486e+000 2.5331976e+000
+ 1.1100627e+000 1.8897138e+000
+ 5.3385926e+000 1.7769978e+000
+ 7.6690500e-001 2.1148008e+000
+ -2.1609282e-001 -2.3028118e+000
+ -1.7495128e+000 -1.6880956e+000
+ -1.2276791e+000 -2.8854234e+000
+ 4.3904839e+000 2.4438104e+000
+ 2.4442360e+000 2.2399624e+000
+ -3.9414221e+000 -3.3299493e+000
+ 5.0365202e+000 1.8322486e+000
+ 2.9415843e+000 9.3112285e-001
+ -2.4768028e+000 -7.9160859e-001
+ -3.2500562e+000 -2.1201864e+000
+ 1.3530053e+000 1.4906673e+000
+ 8.7075787e-001 1.1299585e+000
+ -3.4904035e-001 -2.1185310e+000
+ -2.0825560e+000 -2.4134370e+000
+ 6.5150615e+000 1.4266579e+000
+ -1.9383740e+000 -2.7892562e+000
+ -1.3072571e+000 -2.6673535e+000
+ 2.5443428e+000 3.3878468e+000
+ -3.2151182e+000 -1.8245487e+000
+ -2.1380190e+000 -1.8304378e+000
+ 7.0886822e-001 3.2732594e+000
+ 2.0869857e+000 3.3499188e+000
+ 4.8576381e+000 1.1451490e+000
+ 1.5097927e+000 3.3564114e+000
+ 6.2099705e+000 3.2888367e+000
+ -2.2835024e+000 -3.5069877e+000
+ 5.7801804e+000 5.4443514e-001
+ -1.1897903e+000 -2.4067846e+000
+ 1.8243693e+000 1.2322049e+000
+ 1.9498543e+000 1.4427639e+000
+ -3.0974587e-001 -2.3149875e+000
+ 1.3240267e+000 1.9027747e+000
+ -2.4457520e+000 -2.3857379e+000
+ 9.5619265e-001 9.6334668e-001
+ 4.7209609e+000 2.6739515e+000
+ 4.5226896e-001 3.3999852e+000
+ 2.1990562e+000 2.9844166e+000
+ 5.2299872e+000 1.3525643e+000
+ 5.4114134e+000 2.1733263e+000
+ 1.7869545e+000 3.3292522e+000
+ 1.9789918e+000 3.3629438e+000
+ -5.6124903e-001 2.0971276e-001
+ -3.0349487e+000 -3.0877115e+000
+ 1.9440506e+000 4.9572413e-001
+ 6.7731324e+000 3.4389370e+000
+ -2.6261429e+000 -2.2111932e+000
+ 2.2820509e+000 1.9479899e+000
+ 6.3317721e-001 8.0238616e-001
+ 4.8274070e+000 2.3250519e+000
+ -3.1261923e+000 -1.2801039e+000
+ 4.7468637e+000 2.7772636e+000
+ 5.4841382e+000 1.3103176e+000
+ -1.4529454e+000 -3.4302434e+000
+ -1.4688218e+000 -3.5411790e+000
+ 7.4973540e+000 2.7672322e+000
+ 2.2035000e+000 1.6310828e+000
+ 3.7065309e+000 1.7022017e+000
+ 2.8548444e+000 5.4874939e-001
+ -2.4485611e+000 -1.7206309e+000
+ -1.2053791e+000 -2.1786090e+000
+ 2.0085849e+000 8.4256934e-001
+ -3.8484637e+000 -2.1086352e+000
+ 1.4963432e+000 9.5426817e-001
+ -2.3155129e+000 -6.8639134e-001
+ -5.9169757e-001 -2.6401730e+000
+ -1.4560087e+000 -1.1915188e+000
+ 1.3222089e+000 2.1378826e+000
+ 6.1329397e+000 1.6444192e+000
+ -2.5218939e+000 -1.1597346e+000
+ -3.7632821e+000 -1.3149767e+000
+ -1.3571426e+000 -2.2831657e+000
+ -2.5682677e+000 -2.3125672e+000
+ -2.0947102e+000 -1.8106456e+000
+ 1.4273530e+000 1.3800586e+000
+ -1.1018982e-001 2.5708174e+000
+ -1.9090598e+000 -6.4605216e-001
+ 6.2341419e+000 1.9653485e+000
+ -3.1083165e+000 -3.2911502e+000
+ -1.9259640e+000 -3.6348208e+000
+ 1.7567169e+000 3.5362587e+000
+ 1.8785097e+000 3.6366313e+000
+ -8.7864182e-001 -6.6712643e-001
+ 6.6830071e+000 2.9572910e+000
+ -3.9366542e-001 -1.7748178e+000
+ -1.4997982e+000 -1.6473558e+000
+ 1.4210170e+000 -4.8005848e-001
+ -1.4921494e+000 -8.8918532e-001
+ 5.9607234e-001 1.3074072e-001
+ -1.6147713e+000 -2.6982766e+000
+ 3.4993155e+000 1.8435193e+000
+ 2.5659185e+000 2.2664854e+000
+ -2.7905414e+000 -1.4483953e+000
+ 1.8071201e+000 2.0070071e+000
+ -1.9030077e+000 -1.0091412e+000
+ 9.4530755e-002 9.0722602e-002
+ 6.5851610e+000 2.8508732e+000
+ 1.0121456e+000 3.4038983e+000
+ -4.0070898e+000 -3.0364394e+000
+ 5.5475805e+000 2.0493909e+000
+ 2.4508220e+000 1.0582038e+000
+ -5.3008887e+000 -3.1708097e+000
+ 1.2845523e+000 1.6571738e+000
+ -2.0673754e+000 -2.1732206e+000
+ -1.6679258e+000 -3.7392126e+000
+ 7.2763012e+000 1.7068612e+000
+ 2.8098209e+000 1.5267310e+000
+ 3.3099675e+000 2.5383703e+000
+ -3.1242335e+000 -1.4331335e+000
+ -1.1408029e-001 -1.3527264e+000
+ 2.9880139e+000 1.8625962e+000
+ -2.3994658e+000 -2.6319580e+000
+ -2.2627003e+000 -1.8269851e+000
+ -1.4338143e+000 -2.2929704e+000
+ -1.7645332e-001 -1.9316447e+000
+ 1.9783909e+000 1.3288269e+000
+ 5.4831015e+000 1.4833859e+000
+ 4.0732148e+000 2.5689306e+000
+ 5.5704042e+000 1.6344585e+000
+ 5.4951593e+000 2.3375716e+000
+ 8.1247362e+000 2.8796003e+000
+ -2.3996555e+000 -1.2632269e+000
+ -1.7058610e+000 -1.4606382e+000
+ -1.0673503e+000 -9.5361569e-001
+ 2.4130297e+000 2.1037873e+000
+ 1.0900289e+000 7.9144859e-001
+ 2.7739222e+000 2.2919625e+000
+ 5.6101899e+000 1.8194125e+000
+ 2.7866755e+000 2.4078463e+000
+ 2.5631004e+000 1.2855027e+000
+ 7.7465426e-001 1.3768056e+000
+ -1.7590903e-002 -1.9506899e+000
+ -2.0284008e-001 3.6814302e+000
+ 5.1207441e+000 1.7705272e+000
+ 3.2497475e+000 3.1144415e+000
+ -2.9353943e+000 -3.0437857e+000
+ 5.2272459e+000 2.5673588e+000
+ 1.0988246e+000 1.9263787e+000
+ 2.2658184e+000 3.2118223e+000
+ 1.3895628e+000 1.4343233e+000
+ -1.9936968e+000 -2.9391481e+000
+ 2.4306576e+000 2.5053218e+000
+ -2.6668859e-001 -2.2462924e+000
+ -2.5741085e+000 -2.3077924e+000
+ -3.8521339e+000 -2.3539076e+000
+ -2.7932098e+000 -3.0416005e+000
+ -1.3355889e+000 -1.5445132e+000
+ -9.0139079e-001 -3.4063362e+000
+ 4.6648165e+000 1.7666716e+000
+ 5.8715743e+000 1.8932575e+000
+ 4.9850853e+000 1.3002173e+000
+ 1.4800089e+000 3.6784384e+000
+ 5.4911578e+000 1.9613693e+000
+ 4.4971385e+000 2.2754214e+000
+ 3.0038844e+000 2.6594978e+000
+ 1.3818516e+000 3.1156151e+000
+ 1.7945146e+000 2.1624320e+000
+ -1.3070513e+000 -3.6830996e+000
+ 2.1222048e+000 1.4249542e+000
+ 7.6309758e+000 2.8408909e+000
+ 2.0851244e+000 3.0998481e+000
+ -3.3141936e+000 -1.7103208e+000
+ 2.5886441e+000 2.1135229e+000
+ -2.6638383e+000 -1.3847331e+000
+ -5.2930656e-001 -1.6597116e+000
+ 2.1111343e+000 3.6143535e+000
+ -1.7389801e+000 -4.4606746e+000
+ 2.9482771e+000 1.1911044e+000
+ 1.6976269e+000 3.7333735e+000
+ -2.8446389e-001 -4.3106962e+000
+ 6.7957701e+000 1.6033057e+000
+ -2.0820390e+000 -1.0229105e+000
+ 2.3648434e+000 2.7733358e+000
+ -4.4280365e-001 -1.9246605e+000
+ 2.8269315e+000 1.5943273e+000
+ -2.5467514e+000 -2.3270069e+000
+ 3.7885050e+000 2.6697983e+000
+ -3.5331401e+000 -2.2700009e+000
+ -2.6410138e+000 -2.0310298e+000
+ -2.6593576e+000 -1.9881452e+000
+ 5.0137369e+000 1.3025287e+000
+ -3.1177807e+000 -1.2290185e+000
+ 4.1973888e+000 2.6696191e+000
+ 5.5752286e+000 1.3825153e+000
+ 1.5650922e+000 3.1686117e+000
+ -1.7386775e+000 -3.3908867e+000
+ -2.0224344e+000 -2.1496544e+000
+ 7.8053880e+000 2.6864705e+000
+ 5.8818521e-001 1.9534034e+000
+ 5.3325392e+000 1.4371639e+000
+ -2.3489978e+000 -1.3456940e+000
+ 3.0221691e+000 1.7989975e+000
+ 1.6883578e+000 -5.4709710e-001
+ 8.2627229e-001 1.6146088e+000
+ 7.0639815e+000 1.8590896e+000
+ 3.4099588e+000 2.5455145e+000
+ -1.7709315e+000 -2.9817047e+000
+ 5.7504659e+000 1.9082913e+000
+ -1.0944338e+000 -2.1118448e+000
+ -2.1429571e+000 -3.2043568e+000
+ 8.1435317e-001 1.8027974e+000
+ 2.2813733e-001 1.7991854e+000
+ 1.4022263e+000 3.0138376e+000
+ -2.8648045e+000 -7.8212366e-001
+ -2.1381743e+000 -6.4541059e-001
+ -2.6768040e+000 -6.1500850e-001
+ 1.4524378e+000 2.0915126e+000
+ 3.4141867e+000 2.5366305e+000
+ -3.9624627e+000 -3.9184399e-001
+ -1.5733103e+000 -1.6202287e+000
+ 2.0938662e+000 3.0005550e+000
+ -1.3340284e+000 -1.2626483e+000
+ 5.6034156e+000 1.9454240e+000
+ -1.9020521e+000 -2.7223490e+000
+ -2.3173239e+000 -3.8286721e+000
+ -2.3477349e+000 -3.4712393e+000
+ -1.9006307e+000 -2.6166183e+000
+ 5.5608490e+000 2.1209451e+000
+ 3.9706314e+000 1.4831739e+000
+ -3.6213542e+000 -3.6099191e+000
+ -1.8486038e+000 -3.4576120e+000
+ 6.1982287e+000 2.9441109e+000
+ -1.6789138e+000 -2.9075877e+000
+ -2.2720625e+000 -1.0520459e+000
+ -2.0336044e+000 -2.8188498e+000
+ 5.3168921e+000 8.0860877e-001
+ 1.5890300e+000 2.4358238e+000
+ -1.7266437e-001 2.3396882e+000
+ -2.7189453e+000 -1.9024024e+000
+ 2.7455172e+000 2.5009448e+000
+ 3.1539106e+000 1.7995102e+000
+ -1.9781743e+000 -1.7236511e+000
+ 5.8029030e-001 4.0235214e-001
+ 2.7993007e+000 4.1489746e+000
+ 3.6373453e+000 2.8709177e+000
+ 5.1695180e+000 2.6675799e+000
+ 1.7518959e+000 2.9405589e+000
+ -2.7749800e+000 -3.7380497e+000
+ -2.1925918e+000 -5.5872177e-001
+ 4.9680205e-002 1.0009258e+000
+ 1.5088386e+000 -6.3183873e-001
+ 5.8313662e+000 1.7344139e+000
+ 6.0940067e+000 4.4160856e+000
+ 7.4119343e+000 1.8170559e+000
+ 1.6728295e+000 9.9283878e-001
+ -1.3498552e+000 -2.3857744e+000
+ -1.6930807e+000 -3.2048284e+000
+ 3.1677043e-001 3.6582128e-001
+ 4.1179164e+000 2.6827506e+000
+ 2.7387767e+000 1.3233874e+000
+ 2.1455813e+000 3.1223854e+000
+ 6.2808013e+000 9.6439219e-001
+ 3.3286875e+000 1.9634466e+000
+ -2.8052367e+000 -1.6994091e+000
+ 4.6426313e+000 2.1390181e+000
+ 9.2322499e-001 1.1521913e+000
+ 3.2984268e+000 3.3150415e+000
+ -2.4755271e+000 -1.4209597e+000
+ 2.9549956e+000 2.2586466e+000
+ 1.8416754e+000 1.4605806e+000
+ 1.6953848e+000 1.5910868e+000
+ -2.4832723e+000 -1.7951173e+000
+ 2.7937346e+000 2.8970104e+000
+ -1.4807577e+000 -1.1895410e+000
+ 1.4021233e+000 5.3249333e-001
+ -3.4959091e+000 -9.6476898e-001
+ -1.9768578e+000 -1.0747765e+000
+ -1.8113868e+000 -1.2608870e+000
+ 1.6163478e+000 1.8018900e+000
+ 2.8051045e+000 3.0924390e+000
+ 9.0410544e-001 2.6425777e+000
+ -1.9088397e+000 -1.1575430e+000
+ 1.6523039e+000 6.8678379e-001
+ 2.4121994e+000 1.4706600e+000
+ -2.3872486e+000 -3.7786637e-001
+ 1.9225292e+000 1.8030803e+000
+ 8.2935625e-001 1.7847496e+000
+ -1.4686285e+000 -2.0790260e+000
+ 1.8872933e+000 3.1501926e+000
+ 5.3252071e-001 2.9203923e+000
+ 5.4058398e+000 2.5232075e+000
+ -2.5586116e+000 -1.5706654e+000
+ 2.8718177e+000 2.4386103e+000
+ 5.6064262e+000 1.0461471e+000
+ 1.9552967e+000 3.4415908e-001
+ 2.6367513e+000 1.9379351e+000
+ 4.6343913e+000 1.3367959e+000
+ -3.0768948e+000 -2.1924586e+000
+ -2.0295573e+000 -1.0985035e+000
+ -2.2767839e+000 -1.8612956e+000
+ 3.3866500e+000 1.4294696e+000
+ 3.6638890e-001 3.6043911e+000
+ -3.3547148e+000 -2.2427830e+000
+ 8.4451343e-001 2.2120448e+000
+ -2.2733553e+000 -2.1275993e+000
+ 3.7077744e+000 1.4532053e+000
+ 2.0750831e+000 1.7255075e+000
+ -3.9541444e-001 -5.2310894e-001
+ 1.7078826e+000 1.4087915e+000
+ 2.0020426e+000 4.3416068e+000
+ 1.7627536e+000 -5.7888970e-002
+ 1.9400134e+000 1.9151391e+000
+ -2.6556402e+000 -2.0695271e+000
+ 3.3815487e+000 3.2055909e+000
+ 5.1807595e+000 1.9997884e+000
+ -1.2597354e+000 1.3378316e-001
+ -3.0422685e+000 -3.0964077e+000
+ -7.7716060e-001 -3.2480660e+000
+ 3.7515311e+000 1.4632459e+000
+ 2.6944685e+000 9.7869790e-001
+ -2.1145158e+000 -1.6380541e+000
+ 4.7811731e+000 1.6536735e+000
+ 2.8538678e+000 9.6729593e-001
+ -1.9622444e+000 -1.4903151e+000
+ 3.3843012e+000 2.0393416e+000
+ 1.1045424e+000 3.8625096e+000
+ 5.3249839e+000 2.7923426e+000
+ 3.8722647e+000 1.2604236e+000
+ -1.9538865e+000 -2.8704042e+000
+ 6.4992627e+000 2.7720824e+000
+ -9.2363681e-001 -1.3351875e+000
+ 9.9129619e-001 1.7340198e+000
+ 3.2274900e+000 7.9237024e-001
+ -2.4717007e+000 -1.6165435e+000
+ -2.7237882e+000 -2.2389348e+000
+ -3.2644551e+000 -2.1864121e+000
+ -3.6688457e+000 -3.8340328e+000
+ 3.7240952e+000 2.8045789e+000
+ -2.5660161e+000 -1.1099793e+000
+ -1.4290780e+000 -1.6888942e+000
+ 5.8244186e+000 2.0928257e+000
+ 9.3302711e-001 1.6501715e+000
+ 1.9475450e+000 4.9729855e-001
+ 7.4455072e-001 1.9752648e+000
+ -1.3642994e+000 -2.8881724e+000
+ -3.4342014e+000 -3.9112118e-001
+ 3.1646656e+000 9.5026059e-001
+ 5.3127705e+000 2.5641602e+000
+ -2.5063746e+000 -1.8929695e+000
+ 4.6639911e+000 1.9161795e+000
+ 8.0302263e-001 1.5481499e+000
+ -1.9070068e+000 -1.9313386e+000
+ 2.5571912e+000 1.6621060e+000
+ 5.1409968e+000 1.4730887e+000
+ 3.3476570e+000 1.9673857e+000
+ 1.8358165e+000 1.0285612e+000
+ -1.1315108e+000 -2.9811669e+000
+ -8.6992124e-001 -3.0647003e+000
+ -3.5423969e+000 -2.9739972e+000
+ 5.5608368e+000 1.6360458e+000
+ -9.6189622e-001 -3.4805733e+000
+ -8.1712876e-001 -1.9558336e+000
+ -1.4681838e+000 -3.6384896e+000
+ 1.3660560e+000 2.3828135e+000
+ 2.0696215e+000 1.6372076e+000
+ 4.8528723e+000 1.7102151e+000
+ -2.1629972e+000 -3.5021330e-001
+ -2.9858755e+000 -2.7885458e+000
+ 2.9964730e+000 2.1350066e+000
+ 6.5141146e+000 2.3729766e+000
+ 4.4473609e-001 -2.2042648e+000
+ 6.3946126e+000 1.7683332e+000
+ -1.7777494e+000 1.5873189e-001
+ 4.8749357e+000 1.5012762e+000
+ 2.7333896e+000 4.1793878e+000
+ 2.3038458e+000 1.3971225e+000
+ 5.5531269e+000 1.7464027e+000
+ -2.0140448e+000 -2.6381175e-001
+ 5.6840220e+000 1.7853064e+000
+ -2.5686222e+000 -1.8168543e+000
+ -1.8798701e+000 -3.0767154e+000
+ -2.3492370e+000 -1.6402518e+000
+ 3.2913429e+000 2.3068625e+000
+ 1.2216923e+000 1.6979741e+000
+ 1.6164298e+000 2.5992161e+000
+ 2.3924495e+000 1.7298108e+000
+ -6.3507905e-001 -2.5262922e+000
+ -2.2537553e+000 -2.7913213e+000
+ -1.7101551e+000 -2.9582734e+000
+ 6.8338858e+000 2.2970488e+000
+ 3.7612112e+000 3.5229264e+000
+ 2.6055554e+000 2.6437519e+000
+ -2.6646117e+000 -1.0063774e+000
+ -2.7959630e+000 -1.9289798e+000
+ 3.4492294e+000 1.8299931e+000
+ 1.4055738e+000 2.1684868e+000
+ 1.9859642e+000 2.8128380e-001
+ 6.4602200e+000 2.1482633e+000
+ -3.4891465e+000 -6.3899783e-001
+ -1.2673496e+000 -1.0419070e+000
+ -3.3559617e+000 -1.8620568e+000
+ -1.8518284e+000 -2.6489226e+000
+ 3.1123639e+000 1.4127378e+000
+ -2.2887663e+000 -3.3066939e+000
+ 1.2962181e+000 1.8610075e+000
+ 6.3971671e+000 2.1942553e+000
+ 4.6402583e+000 2.0502507e+000
+ 2.7617788e+000 4.5934726e+000
+ 4.9511852e+000 6.3315327e-001
+ -1.1279799e+000 -2.3304796e+000
+ 2.1915568e+000 2.1577962e+000
+ -2.0310133e+000 -1.8120516e+000
+ 4.4891877e+000 1.7510403e+000
+ 3.5811662e-001 -1.9870999e+000
+ 6.0941246e+000 2.4164995e+000
+ -2.6487619e+000 -2.3043661e+000
+ 1.3207411e+000 1.1351708e+000
+ 2.5918641e+000 2.1906745e+000
+ 5.2927967e+000 2.2018438e+000
+ 5.3067742e+000 2.6189213e+000
+ 2.2627728e+000 1.0902338e+000
+ 2.7551893e-001 2.4462372e+000
+ 8.1053939e-001 3.8392179e+000
+ 1.4403479e+000 2.6155152e+000
+ 1.9611798e+000 1.6502774e+000
+ -2.1055303e+000 -2.8985010e+000
+ -2.7755140e+000 -4.7437548e-001
+ 5.1770484e+000 1.6804222e+000
+ -3.0334386e+000 -2.7373318e+000
+ 6.0637433e-001 2.7778633e+000
+ -1.7452359e+000 -2.1321238e+000
+ 4.8518966e+000 2.5198720e+000
+ -1.8299290e+000 -1.8400312e+000
+ 2.4149744e+000 3.0206535e+000
+ 6.3465231e+000 1.4053450e+000
+ -2.8943339e+000 -2.4860720e+000
+ 5.7884153e+000 1.8736963e+000
+ 2.0305510e+000 6.4606536e-001
+ -3.7321743e+000 -2.3785652e+000
+ 1.8434431e+000 1.5505547e+000
+ -3.1976895e+000 -1.8064030e+000
+ 4.5573185e+000 1.7248030e+000
+ -1.8609879e+000 -3.8910721e-001
+ 2.3178734e+000 1.1252684e+000
+ -1.6480719e+000 -1.8864007e+000
+ -2.9393917e+000 -2.9873542e+000
+ 1.2534204e+000 7.0288814e-001
+ 9.9984186e-001 1.1632261e+000
+ -1.6411692e+000 -3.0248960e+000
+ 3.0504916e+000 1.3565138e+000
+ 2.7923730e+000 1.4833777e+000
+ 4.4248742e+000 2.6004514e+000
+ -1.1830057e+000 -2.7075284e+000
+ 5.7230205e+000 2.9874484e+000
+ 6.9695587e+000 3.3848155e+000
+ 2.6580319e+000 2.1575699e+000
+ 3.7965996e+000 2.0124439e+000
diff --git a/GaussMix/example2/TrainingData2 b/GaussMix/example2/TrainingData2
new file mode 100755
index 0000000000000000000000000000000000000000..00d1c239d1ccda052f661a781d9ebf0da74c5d03
--- /dev/null
+++ b/GaussMix/example2/TrainingData2
@@ -0,0 +1,500 @@
+ -2.1676544e+000 7.2724912e-001
+ -2.2021094e+000 1.6599544e+000
+ -4.1996061e+000 1.9044286e+000
+ 1.0107508e+000 4.0330519e+000
+ -3.8789592e+000 1.4577770e+000
+ -2.3304724e+000 3.6440477e-001
+ -2.2782990e+000 1.0236861e+000
+ -5.4486511e+000 1.4020687e+000
+ -5.4208090e+000 1.4469357e+000
+ 2.4725909e+000 1.0549263e-001
+ -2.3289652e+000 1.7148352e+000
+ -1.7209313e+000 1.0914628e+000
+ -5.6965405e+000 3.4105391e+000
+ 4.2054913e+000 -3.1484994e+000
+ -1.6021021e+000 2.4054353e+000
+ 1.6635118e+000 -2.1337700e+000
+ -4.4200434e-001 2.7968322e+000
+ -9.9996573e-001 4.1078072e+000
+ -1.5412461e+000 2.3529769e+000
+ -5.6475509e+000 1.1828212e+000
+ 5.8742792e+000 -2.1919377e-001
+ -3.1054012e+000 4.2197632e-001
+ 2.2958504e+000 -3.8211168e+000
+ -3.7882433e+000 2.1567046e+000
+ -1.5128058e+000 3.9813699e+000
+ 1.2090458e+000 -1.6553278e+000
+ -2.8474535e+000 1.7660527e+000
+ -2.8726977e+000 5.3859478e-001
+ 2.4244914e+000 -1.0172181e+000
+ 1.0430026e+000 -1.9432317e+000
+ -5.3100208e+000 1.5015664e+000
+ -5.1761773e+000 2.4991713e+000
+ -6.1415075e+000 1.5881890e+000
+ -1.6258894e+000 6.0114982e-001
+ 1.1925656e+000 -2.0203504e+000
+ -2.2087392e+000 -8.9080716e-002
+ 2.0241082e+000 -3.3732095e+000
+ -1.6613388e+000 1.9323536e+000
+ -1.8295480e+000 5.9332688e-001
+ 1.3598777e+000 -1.7556601e-004
+ 7.7880968e-001 -1.6939482e+000
+ -2.7644346e+000 2.1947783e+000
+ 1.2324791e+000 -1.5593517e+000
+ 1.7475843e+000 -3.1296889e+000
+ -1.8393914e+000 1.9064585e+000
+ 1.0439233e+000 -2.0355201e+000
+ -2.9460185e+000 2.5039615e+000
+ 2.6409956e+000 -2.0646046e-001
+ -2.7541621e+000 2.5313755e+000
+ 1.4621537e+000 -2.2849190e+000
+ -6.7819693e+000 1.9361884e+000
+ 1.9339727e+000 -2.4119725e+000
+ -2.6210369e+000 2.0986504e+000
+ -2.3312682e+000 8.6016611e-001
+ 1.4256515e+000 -5.0501904e-001
+ 1.3568230e+000 -1.7702061e+000
+ 2.5878266e+000 -1.9930017e+000
+ -5.8069391e+000 1.9239306e+000
+ -1.5805892e+000 -6.0959397e-001
+ -1.2781726e+000 3.6903877e+000
+ -6.7263108e+000 1.9160322e+000
+ -2.7345208e+000 -3.3408696e-001
+ 3.4692311e+000 -2.7828553e+000
+ 3.0710302e+000 -2.3614453e+000
+ -3.0984970e+000 3.5655461e+000
+ 1.2961392e+000 -2.7175179e+000
+ -3.3204093e+000 7.5518525e-001
+ -6.9055978e-002 7.1634074e-001
+ 1.3432192e+000 -2.2941179e+000
+ -5.3809156e+000 1.8100490e+000
+ -6.1975239e+000 1.2365340e+000
+ -6.5514826e+000 1.9671028e+000
+ 2.1327706e+000 -3.0154789e+000
+ -1.7464415e+000 2.4624324e+000
+ -4.5960539e+000 2.5926307e+000
+ 2.4040046e+000 -2.3422233e+000
+ 1.8241043e+000 -1.2482035e+000
+ -2.7925874e+000 1.7824508e+000
+ 3.7409566e-001 -1.4669049e+000
+ 2.5401382e+000 -3.9814063e-002
+ 1.2359872e+000 -3.0280581e+000
+ -3.8537901e-001 7.2259292e-001
+ 2.8331064e+000 -3.0314877e+000
+ 1.8714528e+000 -1.1256756e+000
+ -5.1038604e+000 1.4229145e+000
+ 3.0681665e+000 -2.7055532e+000
+ -5.5917015e+000 7.4500076e-001
+ 9.1510063e-001 -1.6171342e+000
+ -1.4837529e+000 2.0464535e+000
+ 1.7907894e+000 -2.3393889e+000
+ -4.8038677e+000 1.0875485e+000
+ -6.2281033e+000 1.6463275e+000
+ 1.4578973e+000 -2.2452449e+000
+ -1.7761689e+000 3.4368541e+000
+ 7.5470585e-001 -4.0815469e+000
+ -1.3616265e-003 -5.4460165e-001
+ -1.6598466e+000 2.1128249e+000
+ -8.3257214e-001 -1.1974065e-001
+ -8.2578555e-001 4.0132822e+000
+ -2.6529080e+000 1.2171562e+000
+ 1.8521896e+000 -2.8500123e+000
+ 2.2409172e+000 -4.2718929e+000
+ -1.6682557e-001 -2.0739330e+000
+ 2.8708619e+000 -1.8848118e+000
+ -1.5234618e+000 2.6116790e+000
+ 2.0446424e+000 -5.4153692e-001
+ -6.9581139e-001 3.2948647e+000
+ 2.9341036e+000 -9.2567876e-001
+ -5.1896941e-001 9.2282710e-001
+ -8.7899110e-002 -3.2418894e+000
+ -2.3935435e+000 1.8313160e+000
+ -8.2906145e-001 4.2326294e-001
+ -4.9075976e+000 2.2638321e+000
+ 3.5472163e+000 -2.3683960e+000
+ 2.8912217e+000 -2.1787764e+000
+ -5.6779281e+000 2.5815444e+000
+ 3.0238135e-001 -1.6874784e+000
+ 1.8526717e-001 2.0676069e+000
+ -4.5011558e-001 -1.4671461e-001
+ -1.8981916e+000 2.8982182e+000
+ -3.2223280e+000 2.6312330e+000
+ -6.9014583e+000 2.5055681e+000
+ -2.7001484e+000 2.5064444e+000
+ -4.7623341e+000 2.7741123e+000
+ 1.5956862e+000 -1.8823092e+000
+ 3.0481034e+000 -1.4585313e+000
+ 1.7359365e+000 -1.5389986e+000
+ 1.6950405e+000 -2.0109643e+000
+ -1.8051536e+000 -4.6449935e-001
+ -2.0658223e+000 3.1575390e+000
+ 2.5205639e+000 -2.8708585e+000
+ -6.6502194e+000 1.4026320e+000
+ -4.7845097e+000 2.7395966e+000
+ 8.6258419e-001 -1.5602537e+000
+ -1.8525669e+000 2.1920748e+000
+ 1.6256108e+000 -1.7521262e+000
+ -2.9188315e+000 5.2052763e-001
+ -2.7311700e+000 2.5319829e+000
+ 2.2445654e+000 -2.5114241e+000
+ -1.6239417e+000 2.7683064e+000
+ 7.7968579e-001 -7.9589156e-001
+ 1.6195688e+000 -1.9512063e+000
+ -4.2775232e+000 2.0647848e+000
+ -6.5157689e+000 2.6507175e+000
+ -2.0470925e+000 1.2058557e+000
+ -6.5666890e+000 1.4707612e+000
+ -7.8903656e-001 1.3534165e+000
+ 2.4583356e+000 -3.1693256e+000
+ 2.7504045e+000 -4.0060268e-001
+ -8.5103678e-001 1.4596248e+000
+ 9.1992445e-001 -1.6161139e+000
+ 1.4566911e+000 -3.0546541e+000
+ -1.8986316e+000 1.4420503e+000
+ 2.5338245e+000 -3.0192621e+000
+ 2.8167800e+000 -2.3131932e+000
+ 2.4577494e+000 -1.4165095e+000
+ -3.5032227e+000 1.6840617e+000
+ -1.7605767e+000 -9.8985808e-001
+ -3.5852859e+000 3.0309145e+000
+ 1.7267245e+000 -1.4508199e+000
+ -5.9770031e-002 2.0433837e+000
+ 1.9591338e+000 -1.2454845e+000
+ -1.6621864e+000 2.4359607e+000
+ -2.8969952e+000 1.7847007e-001
+ -2.4075774e+000 1.3060376e+000
+ -2.9544135e+000 1.4190357e+000
+ 2.3932471e+000 -4.4411209e+000
+ 1.4336152e+000 -2.4483331e+000
+ 2.0275232e+000 -1.4394640e+000
+ 7.6676449e-001 -3.1327761e+000
+ -7.5474719e-001 1.1723951e+000
+ 2.1233294e+000 -8.7659978e-001
+ -3.1540385e+000 2.2729508e+000
+ 2.5413090e-001 -1.4152179e+000
+ -6.5945952e+000 2.1106674e+000
+ 6.2425861e-002 3.5318482e-001
+ 4.0751814e+000 -2.6548648e+000
+ -4.5854439e+000 1.8006036e+000
+ -2.7453582e+000 2.9244188e+000
+ 2.1666379e+000 -4.7403909e-001
+ -4.6679148e+000 1.7756871e+000
+ 1.0888032e+000 -3.9636496e+000
+ -5.1348687e+000 3.8794186e+000
+ -5.6681376e+000 2.0832871e+000
+ -5.9954573e+000 1.4158661e+000
+ 1.8533864e+000 -9.3465573e-001
+ 1.5781670e+000 -1.5059691e-001
+ 3.0004589e+000 -4.3959419e+000
+ -9.3114583e-001 2.2872923e+000
+ -2.2507734e+000 3.3390525e+000
+ 9.9088145e-001 -3.2409501e+000
+ 1.0514559e+000 -3.7135556e+000
+ 1.0123068e+000 -1.3890368e+000
+ -4.5620657e+000 3.9066963e+000
+ 2.3769804e+000 -1.5833669e+000
+ -6.1506551e+000 1.6961295e+000
+ 6.1387758e-001 -2.0395316e+000
+ -1.5749111e+000 1.2484395e+000
+ -1.9532372e+000 1.9053486e+000
+ -2.5928959e+000 2.8441799e+000
+ -3.2790305e+000 4.5361662e-002
+ -5.5124509e+000 1.5167045e+000
+ -1.5329760e+000 2.4863803e+000
+ -3.0122685e+000 2.3984111e+000
+ 3.0659268e-001 1.7256360e+000
+ 1.5809877e+000 -3.3411063e+000
+ -2.6259081e+000 3.4351030e+000
+ 1.8474106e+000 -8.3237545e-001
+ 2.4576598e+000 -2.3751503e+000
+ 1.1207859e+000 -2.3325954e+000
+ 1.9587322e+000 -1.1760887e+000
+ -2.2110303e+000 1.9814177e+000
+ -2.9509701e+000 2.6387979e+000
+ -2.9017903e+000 1.6663783e+000
+ 2.6338990e+000 -1.4892379e+000
+ 3.2998999e+000 -3.0112377e+000
+ 1.3395184e+000 -3.3287122e+000
+ -1.9391178e+000 2.8683140e+000
+ -4.9874024e+000 1.1529442e+000
+ -1.6317410e-001 2.3471695e+000
+ 2.5727615e+000 -1.5583052e+000
+ 2.0148416e+000 -3.0736133e+000
+ 2.8343926e+000 -1.9513714e+000
+ -2.6189279e+000 2.5228178e+000
+ 4.4991759e-001 -3.1892390e+000
+ -8.5823286e-001 3.2558474e+000
+ -2.9911366e+000 2.5642154e+000
+ 2.6688936e+000 -1.0339082e+000
+ -4.2991335e+000 1.3584499e+000
+ -3.3770308e+000 7.4409642e-001
+ 1.3019788e+000 -3.7369064e+000
+ 3.1000690e+000 -2.1990308e+000
+ -3.0335165e+000 6.3044597e-001
+ -9.9872160e-001 4.3625315e-001
+ 2.5017067e+000 -2.9251836e+000
+ 3.1699718e+000 -1.6991676e+000
+ -1.6628983e+000 1.2102858e+000
+ -5.8505178e+000 1.5983912e+000
+ 1.0037421e+000 -2.7624714e+000
+ -2.0844665e+000 1.8522690e+000
+ 2.1709219e+000 -7.3666672e-001
+ -2.2796164e+000 1.5626564e+000
+ 2.2811151e+000 -2.1342902e+000
+ -1.6752886e+000 1.2028144e+000
+ -2.6096121e+000 2.5319809e+000
+ 3.5814789e-001 -2.8110915e+000
+ -5.0372794e+000 1.5623347e+000
+ 2.0298215e+000 -8.3147923e-001
+ 2.2782645e+000 -2.8877736e+000
+ -3.4881764e+000 2.9199784e+000
+ 2.2039885e+000 -2.4547136e+000
+ 1.0968223e+000 -3.3404408e+000
+ 2.7036679e+000 -3.0465614e+000
+ -5.8617857e+000 1.3820514e+000
+ 1.9194568e+000 -3.2085153e+000
+ -1.5254846e+000 2.4405091e+000
+ -3.4345758e+000 1.8141590e+000
+ -3.7279401e+000 2.2604515e+000
+ 2.1290572e+000 -1.8710352e+000
+ -4.0723798e+000 1.9418600e+000
+ -6.3040560e-001 3.8419941e+000
+ -2.1586730e+000 3.0083083e+000
+ -9.9527102e-001 2.9667980e+000
+ 2.6223061e+000 -2.5934817e+000
+ -4.9397590e+000 1.9780623e+000
+ -6.2369039e+000 1.8003624e+000
+ -6.4290560e+000 2.7712555e+000
+ 1.4482711e+000 -8.1344993e-001
+ -6.6184573e+000 1.4599140e+000
+ -3.6045096e+000 2.8595287e+000
+ 1.7681283e+000 -1.7266080e+000
+ -2.3335568e+000 1.0827217e+000
+ 1.2867179e+000 -3.2625336e+000
+ -9.2455813e-001 4.0209850e+000
+ 1.9933454e+000 -1.5704483e+000
+ 2.5639669e+000 -3.4105423e+000
+ 1.3887063e+000 -4.4313006e+000
+ -4.3328679e+000 4.2522646e-001
+ 1.8178065e+000 -3.1978440e+000
+ -2.7840824e+000 1.9298148e+000
+ 3.3194759e+000 -9.0262417e-001
+ -9.5278334e-001 3.5754862e+000
+ 1.7242768e+000 -8.0856886e-001
+ -1.9450839e+000 2.1997442e+000
+ 1.7104614e+000 -2.4824556e+000
+ -4.6889698e+000 1.5358873e+000
+ 2.8441577e+000 -5.6668600e-001
+ 1.8027410e+000 -3.3694643e+000
+ 1.3346451e+000 -2.9255453e+000
+ 1.0397731e+000 -1.6834729e+000
+ -1.6467661e+000 1.9999782e+000
+ -2.1819960e+000 -4.0019552e-001
+ -4.8583724e+000 1.5845082e+000
+ 2.0910299e+000 -2.3308899e+000
+ -3.3523693e+000 2.3699674e+000
+ -4.9197002e+000 2.0139531e+000
+ -2.0932192e+000 2.0995407e+000
+ -2.4426226e+000 1.7733960e+000
+ 1.4766185e-003 3.5787343e+000
+ 3.4502809e+000 -3.8711454e+000
+ 2.2134327e+000 -2.0565704e+000
+ -5.8048356e+000 1.3897030e+000
+ -2.3163499e+000 1.5032838e+000
+ -2.0940587e+000 1.5496989e+000
+ 1.3480146e+000 -2.5609421e+000
+ -4.6903045e+000 2.6514904e+000
+ 2.3446133e+000 -8.0380776e-001
+ 1.5650832e+000 -2.7804089e+000
+ 1.4471924e+000 -6.7079514e-001
+ -4.4759503e+000 3.1643315e+000
+ 2.0148164e+000 -1.1357063e+000
+ 1.8485598e+000 -1.8512816e+000
+ 2.8042447e+000 -2.1250771e+000
+ -4.7107714e+000 2.6126934e+000
+ 1.7992479e+000 -2.4952918e+000
+ -6.0357597e+000 1.6716722e+000
+ 9.5642012e-001 -2.3503747e+000
+ 1.9487231e+000 -2.0808151e+000
+ -3.1815369e+000 1.5074330e+000
+ -7.8541774e+000 2.1418548e+000
+ 2.3223498e+000 -3.5492377e+000
+ -2.9974539e+000 1.8832369e+000
+ 8.3873575e-001 -2.3505409e+000
+ -3.5175899e+000 1.2320267e+000
+ -3.0189285e+000 3.8605051e+000
+ 2.1131138e+000 -1.3086544e+000
+ -2.3261940e+000 3.0496692e+000
+ 3.5421193e+000 -2.7721793e+000
+ -2.2891964e+000 3.1861181e+000
+ 2.4734674e+000 -2.3477099e+000
+ -1.4128343e+000 3.1602779e+000
+ 3.3021818e+000 -2.8195258e+000
+ 1.9035311e+000 -1.6749456e+000
+ 1.3082625e+000 -2.8544585e+000
+ -5.2673158e+000 2.0492032e+000
+ -2.7188138e+000 1.2383612e+000
+ 1.0356042e+000 -1.1404494e+000
+ -5.0585118e+000 1.6950885e+000
+ -2.8744535e+000 1.5186952e+000
+ 2.1183767e+000 -1.7233601e+000
+ -3.0712938e+000 1.4818908e+000
+ 2.7749948e+000 -1.6587319e+000
+ -2.4547641e+000 1.7712970e+000
+ -2.2402051e+000 1.4802455e+000
+ -1.9626847e+000 2.1868115e+000
+ -3.3213157e+000 3.0109268e+000
+ 1.4640134e+000 -1.8598314e+000
+ -2.2062488e+000 2.6040769e+000
+ 2.2587109e-001 3.0775939e+000
+ 2.5109662e+000 -1.2267337e+000
+ -3.7456291e+000 1.7845522e+000
+ -1.3096232e+000 2.4817701e+000
+ -4.6888863e+000 2.2578949e+000
+ -1.4978460e+000 2.1700829e+000
+ 1.1977385e+000 -1.8277932e+000
+ -5.6885959e+000 1.1021093e+000
+ 2.1043876e+000 -1.6881736e+000
+ -4.2236324e+000 2.6493849e+000
+ -1.9853225e+000 9.8061757e-001
+ -5.1317505e+000 3.5436603e+000
+ 2.9033653e+000 -1.7000101e+000
+ -4.0332649e+000 1.5698715e+000
+ -7.3930051e+000 9.2185487e-001
+ 2.7476843e+000 -2.8768288e+000
+ -3.4410611e+000 2.0218007e+000
+ -1.7253245e+000 2.8429494e+000
+ -3.3392850e+000 9.8176753e-001
+ 2.3263611e+000 -5.3257268e-001
+ 1.9466126e+000 -2.7863122e+000
+ 3.0263195e+000 -2.2216725e+000
+ -5.9470897e+000 2.6146257e+000
+ -6.9074646e-001 9.0320140e-001
+ -2.9540331e-002 2.2258842e+000
+ -6.2826929e+000 1.1343384e+000
+ -5.2157719e-001 1.2647725e+000
+ -2.8635868e+000 9.1171395e-001
+ 4.7247713e+000 -3.3273389e+000
+ -1.3634923e+000 2.2502324e+000
+ -2.5099342e+000 2.2419773e+000
+ -3.0223070e+000 1.8040082e+000
+ -2.3506676e+000 2.2463678e-001
+ -2.1162973e+000 9.1791454e-001
+ -3.0101287e+000 1.7357305e+000
+ 2.2583022e+000 -3.5637595e+000
+ 2.3424337e+000 -3.2268255e+000
+ 9.2488474e-001 -2.6333742e+000
+ -5.4612692e+000 1.4600298e+000
+ -5.9999225e-001 1.5649417e+000
+ 2.3854943e+000 -5.8940445e-002
+ 2.9160289e+000 -1.9549155e+000
+ 1.2465840e+000 -1.9124123e+000
+ 1.5187157e+000 -1.5136322e+000
+ -5.4419571e+000 1.4361358e+000
+ 2.4932512e+000 -1.1788646e+000
+ -7.6094411e-001 1.0276082e+000
+ -3.3041943e+000 1.7704488e+000
+ -1.7301287e+000 3.3272934e+000
+ -6.3235921e+000 1.7287941e+000
+ -2.2843920e+000 3.6025497e+000
+ 1.9099843e+000 -2.7857718e+000
+ -5.1441693e+000 1.7010004e+000
+ -1.6369830e+000 2.6302618e+000
+ 3.3623845e+000 -2.0901974e+000
+ -5.8351682e+000 1.5365242e+000
+ 2.5045822e+000 -2.4401132e+000
+ 2.9372926e+000 -1.0307887e+000
+ -3.1187666e+000 3.2178220e+000
+ 1.1317001e+000 -1.4572731e+000
+ -2.5215315e+000 1.9947905e+000
+ 2.4809648e+000 -2.6411335e+000
+ 1.3324044e+000 -1.8003361e+000
+ 2.7578116e+000 -3.3578071e-001
+ -1.2156542e+000 3.7800405e+000
+ -3.8921009e+000 1.8846481e+000
+ 9.1211295e-001 -9.6082375e-001
+ -1.0396915e+000 6.2798394e-001
+ 2.1355459e+000 -1.8149112e+000
+ -2.4507887e+000 3.0725035e+000
+ -1.5673305e+000 1.4073507e+000
+ -3.3281293e+000 2.6199991e+000
+ -2.6075175e+000 3.3492980e+000
+ -5.2666275e+000 3.2610538e-001
+ -1.7419816e+000 2.9121948e+000
+ -4.9458458e+000 1.6884673e+000
+ -2.0687261e+000 2.1135684e+000
+ -2.7928258e+000 1.1815234e+000
+ -1.5698826e+000 2.8220756e+000
+ -1.8970752e+000 2.4549007e+000
+ 2.5585756e+000 -3.1280968e+000
+ 1.9873646e+000 -3.8968824e+000
+ -3.1643912e+000 1.8497893e+000
+ 2.5809581e+000 -2.0357697e+000
+ 1.0422829e+000 -1.0119894e+000
+ 4.5268782e-001 3.1636245e+000
+ -2.1484651e+000 1.9616253e+000
+ -1.6585188e+000 3.1767369e+000
+ -9.8745907e-001 1.2290280e+000
+ -4.7119564e+000 1.9489852e+000
+ 1.7637686e+000 -2.3271390e+000
+ -3.3857605e+000 1.5368706e+000
+ 1.5096165e+000 -1.3438428e+000
+ -1.4632883e+000 2.3760604e+000
+ 3.1113228e+000 -8.6507602e-001
+ 1.2817012e+000 -3.6538279e-001
+ -6.5891322e+000 1.3498170e+000
+ -5.6568279e+000 9.0660653e-001
+ 2.6187611e+000 -1.3304174e+000
+ 1.4892796e+000 -1.8202912e+000
+ -2.4259056e+000 1.4058776e+000
+ -1.1296872e+000 8.7096637e-001
+ -1.7561140e+000 1.1498979e+000
+ 1.4567467e+000 -2.6633505e+000
+ 1.2692843e+000 2.8744692e-001
+ -1.6029110e+000 2.3193387e+000
+ -1.2622046e+000 5.0325290e-001
+ -6.2990757e+000 2.3433287e+000
+ 3.2330669e+000 -2.9176328e+000
+ -2.1309259e+000 1.1293462e+000
+ -6.7048757e+000 9.5979888e-001
+ 1.8684839e+000 -3.4170178e+000
+ -2.4939771e+000 1.6319685e+000
+ -5.2412834e+000 2.4549733e+000
+ -5.9657743e+000 1.0863658e+000
+ 3.2579653e+000 -1.8180766e+000
+ 3.5736034e+000 -2.0098089e+000
+ 9.9946088e-001 -1.4850167e+000
+ 2.6688424e+000 -1.8165807e+000
+ -2.4736052e+000 3.7510348e+000
+ 3.1200308e+000 -2.2481138e+000
+ 2.5828428e+000 -2.2578799e+000
+ -2.2862417e+000 1.7999704e+000
+ -4.3055708e+000 1.9845527e+000
+ -1.8802485e+000 3.1818966e+000
+ -3.3446473e+000 3.3598045e+000
+ 1.0468711e+000 -2.1830001e+000
+ 1.3060193e+000 -3.6338176e+000
+ 2.5571198e+000 -3.2733204e+000
+ 3.3239060e+000 -2.3111518e+000
+ 2.5806153e+000 -2.3612028e+000
+ -1.1414758e+000 4.3062739e+000
+ -2.8611504e+000 7.8207517e-001
+ 8.1025286e-001 -4.5564056e+000
+ -1.0536141e+000 1.8137619e+000
+ 1.0038846e+000 -2.5692194e+000
+ -2.0162120e+000 2.5699164e+000
+ 1.5257350e+000 -1.9221790e-001
+ -1.5598715e+000 1.0872521e+000
+ 1.0908929e+000 -1.7086302e+000
+ -6.5388936e+000 8.4605395e-001
+ -3.9671925e+000 2.9819935e+000
+ 2.0225894e+000 -7.0042646e-001
+ -8.4756067e-001 1.5323072e+000
+ 3.4309161e+000 -2.9699768e+000
+ 1.0088562e+000 -5.0425607e+000
+ -6.0803525e+000 1.6557803e+000
+ -1.9215240e+000 2.8981590e+000
+ -2.3888767e-001 4.3991276e+000
+ 3.2583366e+000 -2.5158776e+000
+ -6.3186218e+000 2.4757644e+000
+ 3.8861088e+000 -3.3307862e+000
diff --git a/GaussMix/example2/mkdata.m b/GaussMix/example2/mkdata.m
new file mode 100755
index 0000000000000000000000000000000000000000..a5c1bd94f406a31ee7c627b32d0daf8e4a94e6d3
--- /dev/null
+++ b/GaussMix/example2/mkdata.m
@@ -0,0 +1,74 @@
+
+close all
+clear all
+
+N=500;
+
+R1 = [ 1, 0.1; 0.1, 1];
+mu1 = [2, 2]';
+
+R2 = [ 1, -0.1; -0.1, 1];
+mu2 = [-2, -2]';
+
+R3 = [ 1, 0.2; 0.2, 0.5];
+mu3 = [5.5, 2]';
+
+Pi = [0.4, 0.4, 0.2];
+
+% Generate Data from Gaussian mixture model
+x1 = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3); % training
+y1 = GMM(N/20,Pi,mu1,mu2,mu3,R1,R2,R3); % testing
+
+axis([-8,10,-7,7])
+plot(x1(1,:),x1(2,:),'bo');
+title('Scatter Plot of Multimodal Traning Data for Class 0')
+xlabel('first component')
+ylabel('second component')
+hold on
+
+
+R1 = [ 1, 0.1; 0.1, 1];
+mu1 = [-2, 2]';
+
+R2 = [ 1, -0.1; -0.1, 1];
+mu2 = [2, -2]';
+
+R3 = [ 1, 0.2; 0.2, 0.5];
+mu3 = [-5.5, 2]';
+
+Pi = [0.4, 0.4, 0.2];
+
+% Generate Data from Gaussian mixture model
+x2 = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3); % training
+y2 = GMM(N/20,Pi,mu1,mu2,mu3,R1,R2,R3); % testing
+
+plot(x2(1,:),x2(2,:),'rx');
+title('Scatter Plot of Multimodal Training Data for Class 1')
+xlabel('first component')
+ylabel('second component')
+
+
+% Concatenate testing data from two distributions
+y = [y1,y2];
+
+figure
+plot(y(1,:),y(2,:),'rx');
+title('Scatter Plot of Multimodal Testing Data')
+xlabel('first component')
+ylabel('second component')
+
+
+
+% Save data to files
+x1 = x1';
+save TrainingData1 x1 /ascii
+
+x2 = x2';
+save TrainingData2 x2 /ascii
+
+y = y';
+save TestingData y /ascii
+
+
+
+
diff --git a/GaussMix/example2/rundemo.m b/GaussMix/example2/rundemo.m
new file mode 100755
index 0000000000000000000000000000000000000000..ea6cd1388ab8573051b08e5d73c0662f18e4bae6
--- /dev/null
+++ b/GaussMix/example2/rundemo.m
@@ -0,0 +1,65 @@
+clear all;
+
+if exist('GaussianMixture')~=2
+ pathtool;
+ error('the directory containing the Cluster program must be added to the search path');
+end
+
+disp('generating data...');
+mkdata;
+clear all;
+traindata1 = load('TrainingData1');
+traindata2 = load('TrainingData2');
+testdata = load('TestingData');
+input('press <Enter> to continue...');
+disp(' ');
+
+% [mtrs, omtr] = GaussianMixture(pixels, initK, finalK, verbose)
+% - pixels is a NxM matrix, containing N training vectors, each of M-dimensional
+% - start with initK=20 initial clusters
+% - finalK=0 means estimate the optimal order
+% - verbose=true displays clustering information
+% - mtrs is an array of structures, each containing the cluster parameters of the
+% mixture of a particular order
+% - omtr is a structure containing the cluster parameters of the mixture with
+% the estimated optimal order
+disp('clustering class 1...');
+[mtrs,class1] = GaussianMixture(traindata1, 20, 0, false);
+disp(sprintf('\toptimal order K*: %d', class1.K));
+for i=1:class1.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', class1.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(class1.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(class1.cluster(i).R,6)]);
+end
+input('press <Enter> to continue...');
+disp(' ');
+
+disp('clustering class 2...');
+[mtrs,class2] = GaussianMixture(traindata2, 20, 0, false);
+disp(sprintf('\toptimal order K*: %d', class2.K));
+for i=1:class2.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', class2.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(class2.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(class2.cluster(i).R,6)]);
+end
+input('press <Enter> to continue...');
+disp(' ');
+
+disp('performing maximum likelihood classification...');
+disp('for each test vector, the following calculates the log-likelihood given each of the two classes, and classify');
+disp('the first half of the samples are generated from class 1, the remaining half from class 2');
+disp(' ');
+likelihood=zeros(size(testdata,1), 2);
+likelihood(:,1) = GMClassLikelihood(class1, testdata);
+likelihood(:,2) = GMClassLikelihood(class2, testdata);
+class=ones(size(testdata,1),1);
+class(find(likelihood(:,1)<=likelihood(:,2)))=2;
+for n=1:size(testdata,1)
+ disp([mat2str(testdata(n,:),4), sprintf('\tlikelihood: '), mat2str(likelihood(n,:),4), sprintf('\tclass: %d', class(n))]);
+end
+
+
+
+
diff --git a/GaussMix/example3/GMM.m b/GaussMix/example3/GMM.m
new file mode 100755
index 0000000000000000000000000000000000000000..4fcd1dbcd716f12f4cd4b37d9b665dd4de4e5ded
--- /dev/null
+++ b/GaussMix/example3/GMM.m
@@ -0,0 +1,28 @@
+function x = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3)
+
+pi1 = Pi(1);
+pi2 = Pi(2);
+pi3 = Pi(3);
+
+[V,D] = eig(R1);
+A1 = V*sqrt(D);
+
+[V,D] = eig(R2);
+A2 = V*sqrt(D);
+
+[V,D] = eig(R3);
+A3 = V*sqrt(D);
+
+x1 = A1*randn(2,N) + mu1*ones(1,N);
+
+x2 = A2*randn(2,N) + mu2*ones(1,N);
+
+x3 = A3*randn(2,N) + mu3*ones(1,N);
+
+SwitchVar = ones(2,1)*random('Uniform',0,1,1,N);
+SwitchVar1 = SwitchVar<pi1;
+SwitchVar2 = (SwitchVar>=pi1)&(SwitchVar<(pi1+pi2));
+SwitchVar3 = SwitchVar>=(pi1+pi2);
+
+x = SwitchVar1.*x1 + SwitchVar2.*x2 + SwitchVar3.*x3;
+
diff --git a/GaussMix/example3/mkdata.m b/GaussMix/example3/mkdata.m
new file mode 100755
index 0000000000000000000000000000000000000000..5631d7bac80ab37fb5399356d1ea010a06a29af5
--- /dev/null
+++ b/GaussMix/example3/mkdata.m
@@ -0,0 +1,24 @@
+N=500;
+
+R1 = [ 1, 0.1; 0.1, 1];
+mu1 = [2, 2]';
+
+R2 = [ 1, -0.1; -0.1, 1];
+mu2 = [-2, -2]';
+
+R3 = [ 1, 0.2; 0.2, 0.5];
+mu3 = [5.5, 2]';
+
+Pi = [0.4, 0.4, 0.2];
+
+% Generate Data from Gaussian mixture model
+x = GMM(N,Pi,mu1,mu2,mu3,R1,R2,R3);
+
+plot(x(1,:),x(2,:),'o');
+title('Scatter Plot of Multimodal Data')
+xlabel('first component')
+ylabel('second component')
+
+x = x';
+save data x /ascii
+
diff --git a/GaussMix/example3/rundemo.m b/GaussMix/example3/rundemo.m
new file mode 100755
index 0000000000000000000000000000000000000000..c3cd55e2628a25f1196470480a3013e935de731a
--- /dev/null
+++ b/GaussMix/example3/rundemo.m
@@ -0,0 +1,58 @@
+clear all;
+
+if exist('GaussianMixture')~=2
+ pathtool;
+ error('the directory containing the Cluster program must be added to the search path');
+end
+
+disp('generating data...');
+mkdata;
+clear all;
+pixels = load('data');
+input('press <Enter> to continue...');
+disp(' ');
+
+% [mtrs, omtr] = GaussianMixture(pixels, initK, finalK, verbose)
+% - pixels is a NxM matrix, containing N training vectors, each of M-dimensional
+% - start with initK=20 initial clusters
+% - finalK=0 means estimate the optimal order
+% - verbose=true displays clustering information
+% - mtrs is an array of structures, each containing the cluster parameters of the
+% mixture of a particular order
+% - omtr is a structure containing the cluster parameters of the mixture with
+% the estimated optimal order
+disp('estimating optimal order and clustering data...');
+[mtrs,omtr] = GaussianMixture(pixels, 20, 0, true);
+disp(sprintf('\toptimal order K*: %d', omtr.K));
+for i=1:omtr.K
+ disp(sprintf('\tCluster %d:', i));
+ disp(sprintf('\t\tpi: %f', omtr.cluster(i).pb));
+ disp([sprintf('\t\tmean: '), mat2str(omtr.cluster(i).mu',6)]);
+ disp([sprintf('\t\tcovar: '), mat2str(omtr.cluster(i).R,6)]);
+end
+input('press <Enter> to continue...');
+disp(' ');
+
+% with omtr containing the optimal clustering with order K*, the following
+% split omtr into K* classes, each containing one of the subclusters of
+% omtr
+disp('split the optimal clustering into classes each containing a subcluster...');
+mtrs = SplitClasses(omtr);
+input('press <Enter> to continue...');
+disp(' ');
+
+disp('performing maximum likelihood classification...');
+disp('for each test vector, the following calculates the log-likelihood given each of the classes, and classify');
+disp(' ');
+likelihood=zeros(size(pixels,1), length(mtrs));
+for k=1:length(mtrs)
+ likelihood(:,k) = GMClassLikelihood(mtrs(k), pixels);
+end
+class=ones(size(pixels,1),1);
+for k=1:length(mtrs)
+ class(find(likelihood(:,k)==max(likelihood,[],2)))=k;
+end
+for n=1:size(pixels,1)
+ str = [mat2str(pixels(n,:),4), sprintf('\tlikelihood: '), mat2str(likelihood(n,:),4), sprintf('\tclass: %d', class(n))];
+ disp(str);
+end
diff --git a/GaussMix/initMixture.m b/GaussMix/initMixture.m
new file mode 100755
index 0000000000000000000000000000000000000000..7506efa1e36c17f2b44e441af4b56ce34bead909
--- /dev/null
+++ b/GaussMix/initMixture.m
@@ -0,0 +1,23 @@
+function mixture = initMixture(pixels, K)
+
+[N M] = size(pixels);
+mixture.K=K;
+mixture.M=M;
+R = (N-1)*cov(pixels)/N;
+mixture.Rmin = mean(diag(R))/1e5;
+cluster(1).N=0;
+cluster(1).pb = 1/K;
+cluster(1).mu = pixels(1,:)';
+cluster(1).R=R+mixture.Rmin*eye(mixture.M);
+if K>1
+ period = (N-1)/(K-1);
+ for k=2:K
+ cluster(k).N=0;
+ cluster(k).pb = 1/K;
+ cluster(k).mu = pixels(floor((k-1)*period+1),:)';
+ cluster(k).R=R+mixture.Rmin*eye(mixture.M);
+ end
+end
+mixture.cluster = cluster;
+mixture = ClusterNormalize(mixture);
+
diff --git a/GaussMix/initMixtureRange.m b/GaussMix/initMixtureRange.m
new file mode 100755
index 0000000000000000000000000000000000000000..1389f550ada03dd2683d144b9789cb8040126099
--- /dev/null
+++ b/GaussMix/initMixtureRange.m
@@ -0,0 +1,28 @@
+function mixture = initMixtureRange(pixels, K)
+
+[N, M] = size(pixels);
+mixture.K=K;
+mixture.M=M;
+R = (N-1)*cov(pixels)/N;
+mixture.Rmin = mean(diag(R))/1e5;
+
+[~,imax]=max(mean(abs(pixels)));
+[~,Ipixels]=sort(pixels(:,imax));
+Nslice=floor(N/K);
+res=mod(N,Nslice);
+istart=floor(res/2);
+Nrand=ceil(Nslice/10);
+
+for k=1:K
+ Ik=Ipixels(istart+(k-1)*Nslice+1:istart+k*Nslice);
+ Irand=randi(Nslice,Nrand,1); %Nrand samples among Nslice in slice Ik
+ cluster(k).N=0;
+ cluster(k).pb = 1/K;
+ cluster(k).mu = mean(pixels(Ik(Irand),:))';
+ cluster(k).R=R+mixture.Rmin*eye(mixture.M);
+end
+
+
+mixture.cluster = cluster;
+mixture = ClusterNormalize(mixture);
+
diff --git a/LFPSpikeDemixing.m b/LFPSpikeDemixing.m
new file mode 100644
index 0000000000000000000000000000000000000000..8caebace8136c22f4504322db4a8d043c23062d2
--- /dev/null
+++ b/LFPSpikeDemixing.m
@@ -0,0 +1,76 @@
+function [paramVBDS,Spike,w]=LFPSpikeDemixing(xbr,xb,Fs,ideltasr,gB,wfeat,paramVBDS,indsupp,verbose)
+
+T=length(xb);
+Tr=length(xbr);
+gBr(1:Tr/2+1,1)=gB(round(((1:Tr/2+1)/(Tr/2+1)*(T/2+1))));
+gBr(Tr/2+2:Tr,1)=flipud(gBr(2:Tr/2));
+ns=paramVBDS.ds*Fs/1000;
+wname='sym6';
+L=6;
+K_init=paramVBDS.K;
+Lc = 1+paramVBDS.nbfeat+paramVBDS.nbfeat*(paramVBDS.nbfeat-1)/2;
+
+
+Klow=max(K_init-paramVBDS.rangeK,1);
+Khigh=K_init+paramVBDS.rangeK;
+
+risstrial=NaN(Khigh,paramVBDS.nbiter);
+
+
+% classif/despiking iterations
+for trial=1:paramVBDS.nbiter
+ for K=Klow:Khigh
+
+ [~, optmixture] = GaussianMixture(wfeat, K+10, K, false);
+
+ % VBDS
+ optmixture=RStep(optmixture,paramVBDS.Cs);
+ [wr,param] = VBClassifDespike7w(xbr',optmixture,ideltasr,paramVBDS,gBr,Fs,wname,L);
+
+ risstrial(K,trial)=(-param.ll+0.5*(K*Lc-1)*log(length(ideltasr)*paramVBDS.nbfeat));
+
+ if verbose
+ disp(['fin K=',num2str(K),', trial=',num2str(trial), '- riss: ', num2str(risstrial(K,trial)), '- nbr: ', num2str(sum(param.r))])
+ end
+
+
+ paramK{K,trial}.param=param;
+
+ end
+end
+
+
+% Best result selection (min rissanen) and shaping
+[~,imin]=min(risstrial(:));
+tbest=ceil(imin/Khigh);
+Kbest=rem(imin-1,Khigh)+1;
+param=paramK{Kbest,tbest}.param;
+
+% final w building
+w=ifft(fft(xb'-param.mb/2^(L/2))./fftshift((param.gam*gB*T+param.varb)));
+w=param.gam*ifft(fft(w).*fftshift(T*gB));
+indSpike=setdiff((1:T),indsupp);
+w(indSpike)=wr;
+w=ifft(fft(w-param.mb/2^(L/2))./fftshift((param.gam*gB*T+param.varb)));
+w=param.gam*ifft(fft(w).*fftshift(T*gB));
+
+% posterior spike temporal waveforms
+csr=zeros(size(paramVBDS.c));
+csr(paramVBDS.tabcs)=param.wSpike;
+
+yr=waverec(csr,paramVBDS.l,wname);
+idx = bsxfun(@plus,ideltasr',-ns/2+1:ns/2)';
+Spike=yr(idx)';
+
+% mean and std of spike clusters
+S0=param.piapost'*Spike./repmat(sum(param.piapost)',1,ns);
+varskt=zeros(Kbest,1);
+for k=1:Kbest
+ varskt(k)=sum(param.piapost(:,k)'*(Spike-S0(k,:)).^2)/(ns*sum(param.piapost(:,k)));
+end
+
+% output structure
+paramVBDS.K=Kbest;
+paramVBDS.S0=S0;
+paramVBDS.varskt=varskt;
+paramVBDS.param=param;
diff --git a/VBClassifDespike7w.m b/VBClassifDespike7w.m
new file mode 100755
index 0000000000000000000000000000000000000000..8fb8156162541f48a53c972c43b881b7ee8fac7b
--- /dev/null
+++ b/VBClassifDespike7w.m
@@ -0,0 +1,159 @@
+function [w,param] = VBClassifDespike7w(y,minit,ideltas,paramFeat,g,Sr,wname,level)
+%UNTITLED Summary of this function goes here
+% Detailed explanation goes here
+% Passage en ondelettes pour la classif
+
+nbiter=5; %usually enough
+T=length(y);
+K=minit.K;
+Ns=length(ideltas);
+ns=size(paramFeat.tabcs,2);
+gam=1;
+r=zeros(Ns,1);
+covsmem=zeros(nbiter,1);
+
+
+% y in wav
+[cy,l] = wavedec(y,level,wname);
+
+% w initialization: filtered version of y
+[b,a]=butter(1,50/Sr*2);
+w=filtfilt(b,a,y);
+
+% hyperparameter initialization
+pik=ones(K,1)/K;
+mb=mean(y);
+cmb=[2^(level/2)*mb*ones(1,l(1)) zeros(1,sum(l(2:level+1)))]';
+cw=wavedec(w,level,wname);
+varb=var(y)/1000;
+
+
+%wav spike init
+Cs=cy(paramFeat.tabcs)-cw(paramFeat.tabcs)-cmb(paramFeat.tabcs);
+
+%variable initialization
+piapost=minit.pnk;
+Sw0=(piapost'*Cs)./sum(piapost)';
+eyeK=eye(ns);
+eyeK=repmat(eyeK(:)',Ns,1);
+eyeK=reshape(eyeK,Ns,ns,ns);
+invSw=zeros(K,ns,ns);
+CovSw=zeros(K,ns,ns);
+for k=1:K
+ invSw(k,:,:)=minit.cluster(k).invR;%*1e1; % may help for more tightened spike classes
+ CovSw(k,:,:)=minit.cluster(k).R;%*1e-1;
+end
+reg=max(max(max(abs(CovSw))));
+wSpike=zeros(Ns,ns);
+covSwapost=zeros(Ns,ns,ns);
+
+% Some structures of ones to facilitate the posterior calculus
+circone=zeros(ns*ns,ns);
+tocirc=[ones(ns,1);zeros((ns-1)*ns,1)];
+tocirc2=[ones(ns,1);zeros((K-1)*ns,1)];
+for i=1:ns
+ circone(:,i)=circshift(tocirc',[0 (i-1)*ns])';
+end
+for i=1:K
+ circoneK2(i,:)=circshift(tocirc2',[0 (i-1)*ns])';
+end
+circoneK=repmat(circone,1,K);
+indtrace=(1:ns:ns*ns)+(0:ns-1);
+
+it=0;
+while (it<nbiter)
+ it=it+1;
+ indr=find(r==1);
+ indk=find(r==0);
+ %% spike posterior:
+ % data
+ Cs=cy(paramFeat.tabcs)-cw(paramFeat.tabcs)-cmb(paramFeat.tabcs);
+ invSwapost=1/varb*eyeK+reshape(piapost*invSw(:,:),Ns,ns,ns);
+ % post
+ covSW0=reshape((reshape(permute(invSw,[1 3 2]),K,ns*ns)*(circoneK.*reshape(repmat(Sw0,1,ns*ns)',ns*ns,ns*K))).*circoneK2,1,K^2*ns);
+ covSW0(covSW0==0)=[];
+ covSW0=piapost*reshape(covSW0',ns,K)';
+
+ for n=1:Ns
+ covSwapost(n,:,:)=eye(ns)/squeeze(invSwapost(n,:,:));
+ wSpike(n,:)=squeeze(covSwapost(n,:,:))*(1/varb*Cs(n,:)+covSW0(n,:))';
+ end
+
+
+ %% posterior on labels Z
+ wSpikez=wSpike(:,paramFeat.plage);
+ CovSwz=CovSw(:,paramFeat.plage,paramFeat.plage);
+ SW0z=Sw0(:,paramFeat.plage);
+
+
+ piapost=zeros(Ns,K);
+ for k=1:K
+ Covsk=squeeze(CovSwz(k,:,:));
+ piapost(:,k)=1/sqrt(det(Covsk))*exp(-1/2*diag((wSpikez-SW0z(k,:))/Covsk*(wSpikez-SW0z(k,:))'));
+ end
+
+ piapostnn(indk,:)=repmat(pik',Ns-sum(r),1).*piapost(indk,:); % piapost non normalized for ll
+ piapost(indk,:)=piapostnn(indk,:)./repmat(sum(piapostnn(indk,:),2),1,K);
+ piapost(indr,:)=0;
+
+ r(find(isnan(sum(piapost,2))))=1;
+ indr=find(r==1);
+ indk=find(r==0);
+ piapost(indr,:)=0;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ %% w update: fft filtering
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ cw=cy-cmb;
+ cw(paramFeat.tabcs)=cw(paramFeat.tabcs)-wSpike;
+ w=waverec(cw,l,wname);
+ w=ifft(fft(w)./fftshift((gam*g*T+varb)));
+ w=gam*ifft(fft(w).*fftshift(T*g));
+ cw=wavedec(w,level,wname);
+
+
+ %% hyperparameters
+ %wav spikes
+ for k=1:K
+ if (sum(piapost(indk,k))>1&&pik(k)>0.001)
+ Sw0(k,:)=1/sum(piapost(indk,k))*(piapost(indk,k))'*wSpike(indk,:);
+ CovSw(k,:,:)=(repmat(piapost(indk,k),1,ns).*(wSpike(indk,:)-Sw0(k,:)))'*(wSpike(indk,:)-Sw0(k,:))/sum(piapost(indk,k));
+ invSw(k,:,:)=eye(ns)/(squeeze(CovSw(k,:,:))+reg*1e-5*eye(ns,ns));
+ else
+ CovSw(k,:,:)=0;
+ end
+ end
+ covsmem(it)=sum(sum(sum((CovSw.^2))));
+
+ % noise (wav)
+ cb=cy-cw;
+ cb(paramFeat.tabcs)=cb(paramFeat.tabcs)-wSpike;
+ mb=mean(cb(1:l(1)));
+ cmb=[mb*ones(1,l(1)) zeros(1,sum(l(2:level+1)))]';
+ eigbm=sum(1./(1/gam*ones(T,1)./(g*T)+1/varb*ones(T,1)));
+ varb=var(cb(1:l(4))-cmb(1:l(4)))+eigbm/T+sum(sum(covSwapost(:,indtrace)))/T;
+
+ % gamma (time)
+ eigb=sum(1./(gam*g*T+varb));
+ gam=varb*gam*eigb/T+w'*ifft(fft(w)./fftshift((gam*g*T)))/T;
+
+ %% les proportions
+ pik=sum(piapost(indk,:))'/(Ns-sum(r));
+end
+
+
+param.it=it;
+param.varb=varb;
+param.mb=mb;
+param.vars=covsmem;
+param.r=r;
+param.Sw0=Sw0;
+param.wSpike=wSpike;
+param.gam=gam;
+param.piapost=piapost;
+param.pik=pik;
+param.ll=sum(log(sum(piapostnn(indk,:),2)));
+
+
+
+
diff --git a/VBDS.m b/VBDS.m
new file mode 100644
index 0000000000000000000000000000000000000000..b42b0ba74f15a453c799a86327a0ec83307c603f
--- /dev/null
+++ b/VBDS.m
@@ -0,0 +1,47 @@
+function [paramVBDS,Spike,w]=VBDS(xb,Fs,paramVBDS,display,verbose)
+
+
+
+%%
+%Band pass filtering to highlight the spikes and thresholding
+[xb_spkband,spkband,spk,paramVBDS]=spikeExtractor(xb,Fs,paramVBDS,display,verbose);
+
+
+%% Spike forms registering and data reduction
+% leave intact data outside spike intervals
+[xbr,ideltasr,indsupp]=spikeRegistering(xb,xb_spkband,Fs,spkband,spk,paramVBDS.ideltas,verbose);
+
+
+if isempty(paramVBDS.ideltas)
+ disp('!! NO (OR VERY FEW) SPIKES DETECTED: PROCESS ABORTED !!')
+else
+%% wav feat extraction (with outliers rejection)
+[wfeat,paramVBDS,ideltasr]=featureExtraction(xbr,Fs,ideltasr,paramVBDS,verbose);
+
+
+%% Find a good value for g (LFP power Spectrum)
+% Zanos, 2011 according to the method in the appendix
+gB = fitLFPpowerSpectrum(xb',.01,2*10e3,Fs,display);
+gB=fftshift(gB);
+
+%% Initial guess on the number of classes (using GMM-EM)
+
+if paramVBDS.mode
+
+KGMM=zeros(10,1);
+for t=1:10
+ [~,mixture] = GaussianMixture(wfeat, 20, 0, false);
+ KGMM(t)=mixture.K;
+end
+
+paramVBDS.K=ceil(mean(KGMM));
+
+end
+
+
+%%
+% Classif and despiking (with iterations if mode=1)
+%
+[paramVBDS,Spike,w]=LFPSpikeDemixing(xbr,xb,Fs,ideltasr,gB,wfeat,paramVBDS,indsupp,verbose);
+
+end
\ No newline at end of file
diff --git a/example.mat b/example.mat
new file mode 100644
index 0000000000000000000000000000000000000000..0e15c40bd97efaa6fb2d55299ee1889b7494f839
Binary files /dev/null and b/example.mat differ
diff --git a/featureExtraction.m b/featureExtraction.m
new file mode 100644
index 0000000000000000000000000000000000000000..c2b92d78925c912dab3013f2e7b4720c00ec31b4
--- /dev/null
+++ b/featureExtraction.m
@@ -0,0 +1,107 @@
+function [wfeat,paramVBDS,ideltasr]=featureExtraction(xbr,Fs,ideltasr,paramVBDS,verb)
+%% Wavelet feature extraction (within [100 3000] frequency band - approximation excluded
+%% Rejection procedure in feature space based on an "ad hoc" nearest neighbors strategy
+%% !! Matlab Wavelet toolbox required !!
+%
+% input:matrix
+% xbr: raw data (pre-processed by 'spikeRegistering.m' function)
+% Fs: Frequency sammpling
+% ideltasr: maximum time samples of detected spikes in xbr
+% paramVBDS: wavelet feature extraction parameters
+% verb: verbose mode
+%
+% output:
+% wfeat: raw data with spike waveforms registered
+% paramVBDS: wavelet feature extraction parameters
+
+
+wname='sym6';
+[hL,hD,rL,rD] = wfilters(wname);
+dwtmode('per');
+
+L=6;
+[c,l]=wavedec(xbr,L,wname);
+cs=zeros(size(c));
+
+Ns=length(ideltasr);
+nsw=(-1.5*Fs/1000:3*Fs/1000);
+indcs=indwavsupp(ideltasr(1)+min(nsw),ideltasr(1)+max(nsw)-1,l,wname,12);
+tabcs=zeros(Ns,length(indcs));
+tabcs(1,:)=indcs;
+
+
+for i=2:Ns
+ tabcs(i,:)=indwavsupp(ideltasr(i)+min(nsw),ideltasr(i)+max(nsw)-1,l,wname,12);
+ if verb
+ if rem(i,1000)==0
+ disp([num2str(i) '/' num2str(Ns) ' spikes processed']);
+ end
+ end
+end
+Cs=reshape(c(tabcs(:)),Ns,size(tabcs,2));
+
+
+
+% selection of highest variance coeffs
+indselec=indwavsupp(ideltasr(1)-Fs/2000,ideltasr(1)+Fs/2000-1,l,wname,16);
+indselec=mod(find(indselec==indcs'),length(indcs));
+vcs=var(Cs(:,indselec));
+[~,Is]=sort(vcs,'desc');
+plage=indselec(Is(1:paramVBDS.nbfeat));
+wfeat=Cs(:,plage);
+
+
+
+
+% élimination outlier - ppv
+mind=zeros(Ns,1);
+for ii=1:Ns
+ dist=sort(sqrt(sum(((wfeat-wfeat(ii,:)).^2)')));
+ mind(ii)=dist(ceil(Ns/100));
+end
+th=mean(mind)+2.64*std(mind);
+Ir=find(mind>th);
+
+wfeat(Ir,:)=[];
+Cs(Ir,:)=[];
+tabcs(Ir,:)=[];
+
+paramVBDS.c=c;
+paramVBDS.l=l;
+paramVBDS.Cs=Cs;
+paramVBDS.tabcs=tabcs;
+paramVBDS.plage=plage;
+ideltasr(Ir)=[];
+
+% ideltas update
+paramVBDS.ideltas(Ir)=[];
+
+end
+
+
+
+function ind = indwavsupp(start,stop,l,wname,extent)
+% finds the wavelet coeff indices supporting the [start-stop] time interval (in
+% samples) - provide a fixed number of coeff for a given interval length
+% l = auxiliary vector (output of wavedec)
+% wname = wavelet name
+% extent: factor restricting side coefficients
+
+[h,~]=wfilters(wname);
+Nf=length(h);
+L=length(l)-2;
+interval=stop-start+1;
+ind=[];
+
+for ll=1:L
+ cinf=max(round((start-Nf/extent)/2),1);
+ csup=min(round(cinf+interval/2^ll+Nf/12),l(end-ll));
+ ind=[ind, sum(l(1:L-ll+1))+cinf:sum(l(1:L-ll+1))+csup];
+ start=cinf;
+end
+%ind=[ind cinf:csup]; % add approx level for features
+ind(ind>sum(l(1:L)))=[]; % remove highest details level for features
+
+end
+
+
diff --git a/fitLFPpowerSpectrum.m b/fitLFPpowerSpectrum.m
new file mode 100755
index 0000000000000000000000000000000000000000..80330c36f5532a3b68461fd64f7619db06a7ef86
--- /dev/null
+++ b/fitLFPpowerSpectrum.m
@@ -0,0 +1,75 @@
+function [theg,x] = fitLFPpowerSpectrum(y,Fl,Fh,Fs,display)
+ %[g,params] = fitLFPpowerSpectrum(y,Fl,Fh,Fs)
+ %
+ %Fits the PSD of an LFP y with a function:
+ %psd = a*log(1+exp(b*(logfreq-c)) + d
+ %
+ %The function is fit only in the range Fl to Fh. It goes to d at low
+ %frequencies and to d+a*b*(logfreq-c)) (log-linear) at high frequencies
+ %Fs is the sampling frequency
+ %
+ % g can then be fed to despikeLFP
+ %
+ % History:
+ %
+ % 07/07/2010: Updated first line of documentation and fixed second
+ % argument
+ % 28/06/2010: v1.0
+ %
+ % Created by Patrick Mineault
+ %
+ % See also despikeLFP
+
+
+ ffty = fft(y);
+ allfreqs = (1:length(y)/2)'/length(y)*Fs;
+ freqs = allfreqs(allfreqs > Fl & allfreqs < Fh);
+ rg = 1+round(freqs*length(y)/Fs);
+ logpsdy = log(ffty(rg).*conj(ffty(rg)));
+
+ g = @(fs,x) x(1)*log( 1 + exp(x(2)*(log(fs)-x(3)))) + x(4);
+
+ %Guestimate the optimal parameters
+ smoothpsdy = ksr(log(freqs),logpsdy);
+ [gm,maxidx] = max(diff(smoothpsdy.f));
+ c0 = smoothpsdy.x(maxidx+1);
+ d0 = max(smoothpsdy.f);
+ b0 = 1;
+ a0 = (min(smoothpsdy.f)-d0)/(log(Fh)-c0);
+
+ %Initialize guesses for parameters
+ x0 = [a0;b0;c0;d0];
+
+ %The inverse density of points in a spot is ~ diff(log(freqs)), meaning that
+ %higher frequencies disproportionally affect the fit on a log scale.
+ %Correct for this.
+ idensity = sqrt(diff(log(freqs)));
+ idensity = [idensity(1);idensity];
+
+ objectivefun = @(x) idensity.*(g(freqs,x) - logpsdy);
+ %Feed that into lsqnonlin
+ x = lsqnonlin(objectivefun,x0);
+
+ allrg = 1+round(allfreqs*length(y)/Fs);
+ logpsdy = log(ffty(allrg).*conj(ffty(allrg)));
+ smoothpsdy = ksr(log(allfreqs),logpsdy);
+ dbfactor = 10/log(10);
+ smallfreqs = exp(log(min(allfreqs)):.05:log(max(allfreqs)));
+
+ themax = max(dbfactor*smoothpsdy.f);
+
+ if display
+ figure;
+ semilogx(allfreqs,dbfactor*logpsdy-themax,exp(smoothpsdy.x),dbfactor*smoothpsdy.f-themax,smallfreqs,dbfactor*g(smallfreqs,x)-themax);
+ box off;
+ xlabel('Frequency (Hz)');
+ ylabel('PSD (dB relative to max)');
+ legend('Empirical PSD of wideband signal','Smoothed PSD of wideband signal','Estimated PSD of LFP');
+ legend(gca,'Location','SouthWest');
+ legend boxoff;
+ xlim([min(allfreqs)/3,max(allfreqs)*3]);
+ end
+
+ gx = exp(g(allfreqs,x));
+ theg = [gx(1);gx;gx(end-1:-1:1)];
+end
diff --git a/ksr.m b/ksr.m
new file mode 100755
index 0000000000000000000000000000000000000000..9ace0844ee8d6cd606285a4fd6a564d4e3418a7c
--- /dev/null
+++ b/ksr.m
@@ -0,0 +1,104 @@
+function r=ksr(x,y,h,Norx)
+% KSR Kernel smoothing regression
+%
+% r=ksr(x,y) returns the Gaussian kernel regression in structure r such that
+% r.f(r.x) = y(x) + e
+% The bandwidth and number of samples are also stored in r.h and r.n
+% respectively.
+%
+% r=ksr(x,y,h) performs the regression using the specified bandwidth, h.
+%
+% r=ksr(x,y,h,n) calculates the regression in n points (default n=100).
+%
+% Without output, ksr(x,y) or ksr(x,y,h) will display the regression plot.
+%
+% Algorithm
+% The kernel regression is a non-parametric approach to estimate the
+% conditional expectation of a random variable:
+%
+% E(Y|X) = f(X)
+%
+% where f is a non-parametric function. Based on the kernel density
+% estimation, this code implements the Nadaraya-Watson kernel regression
+% using the Gaussian kernel as follows:
+%
+% f(x) = sum(kerf((x-X)/h).*Y)/sum(kerf((x-X)/h))
+%
+% See also gkde, ksdensity
+
+% Example 1: smooth curve with noise
+%{
+x = 1:100;
+y = sin(x/10)+(x/50).^2;
+yn = y + 0.2*randn(1,100);
+r=ksr(x,yn);
+plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--','linewidth',2)
+legend('true','data','regression','location','northwest');
+title('Gaussian kernel regression')
+%}
+% Example 2: with missing data
+%{
+x = sort(rand(1,100)*99)+1;
+y = sin(x/10)+(x/50).^2;
+y(round(rand(1,20)*100)) = NaN;
+yn = y + 0.2*randn(1,100);
+r=ksr(x,yn);
+plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--','linewidth',2)
+legend('true','data','regression','location','northwest');
+title('Gaussian kernel regression with 20% missing data')
+%}
+
+% By Yi Cao at Cranfield University on 12 March 2008.
+%
+
+% Check input and output
+error(nargchk(2,4,nargin));
+error(nargoutchk(0,1,nargout));
+if numel(x)~=numel(y)
+ error('x and y are in different sizes.');
+end
+
+x=x(:);
+y=y(:);
+% clean missing or invalid data points
+inv=(x~=x)|(y~=y);
+x(inv)=[];
+y(inv)=[];
+
+% Default parameters
+if nargin<4
+ Norx=100;
+end
+r.n=length(x);
+if nargin<3 || isnan(h)
+ % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
+ hx=median(abs(x-median(x)))/0.6745*(4/3/r.n)^0.2;
+ hy=median(abs(y-median(y)))/0.6745*(4/3/r.n)^0.2;
+ h=sqrt(hy*hx);
+elseif ~isscalar(h)
+ error('h must be a scalar.')
+end
+r.h=h;
+
+% Gaussian kernel function
+kerf=@(z)exp(-z.*z/2)/sqrt(2*pi);
+
+if numel(Norx) == 1
+ r.x=linspace(min(x),max(x),Norx);
+else
+ r.x = Norx(:)';
+end
+N = length(r.x);
+r.f=zeros(1,N);
+for k=1:N
+ z=kerf((r.x(k)-x)/h);
+ r.f(k)=sum(z.*y)/sum(z);
+end
+
+% Plot
+if ~nargout
+ plot(r.x,r.f,'r',x,y,'bo')
+ ylabel('f(x)')
+ xlabel('x')
+ title('Kernel Smoothing Regression');
+end
diff --git a/main.m b/main.m
new file mode 100644
index 0000000000000000000000000000000000000000..6825569fb8814f83d9e95634ab08cbb7910b32d2
--- /dev/null
+++ b/main.m
@@ -0,0 +1,126 @@
+clearvars;
+close all;
+
+%% Some parameters %%
+
+notch=true; % line noise filtering if true
+Fline=50; % frequency of the line noise
+nFline=3; % number of line noise harmonics to be removed
+paramVBDS.spike_band=[300 3000]; % set the frequency band for spikes detection
+paramVBDS.coef=4; %coef for robust threholding (spike detection)
+paramVBDS.coefH=5; % coef for outlier rejection (rejection thr = coefH*detection thr)
+paramVBDS.ds=3; % set spike duration in ms (3ms or above)
+paramVBDS.nbfeat=6; %number of wav features for classif
+paramVBDS.mode=1; % 0: number of class given by user, 1: automatic mode (with nbiter iterations)
+paramVBDS.K=2; % considered if mode=0
+paramVBDS.nbiter=10;
+paramVBDS.rangeK=1;
+
+display=true; % display some figures during processing
+verbose=true; % verbose mode
+
+%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%
+
+%% DATA LOADING AND PRE-PROCESSING
+
+% Two data examples: 20sec snippets of simulation from Le Cam et. al. 2023 and
+% simulated data from Camunas et. al. 2014. (De-)comment one or the other.
+% !!Warning!!: second example has 10 minutes duration: much much longer to
+% process
+
+%%% SIMULATED DATA Le Cam et. al.
+% DATA_PATH='./';
+% DATA_NAME='example.mat';
+% load([DATA_PATH DATA_NAME],'xb','Fs');
+%%%
+
+%%% SIMULATED DATA CAMUNAS (EXAMPLE 4)
+paramVBDS.coef=4.5; % better parameters for these simulated data
+paramVBDS.nbfeat=12;
+load('simulation_4.mat');
+xb=data;
+Fs=24000;
+%%%
+
+[l,c]=size(xb);
+
+if (l~=1 && c~=1)
+ warning('data should be a single channel');
+end
+
+
+
+% Resampling (12kHz)
+xb=resample(xb,12e3,Fs);
+Fs=12e3;
+
+
+% Notch filter
+if notch
+ for i=1:nFline
+ xb=notchfilter(xb,Fs,i*Fline);
+ end
+end
+clear i
+
+T=2*floor(length(xb)/2);;
+xb=xb(1:T); %make the number of samples even
+
+[l,c]=size(xb);
+if l>c
+ xb=xb';
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%
+
+%% VBDS
+[paramVBDS,Spike,w]=VBDS(xb,Fs,paramVBDS,display,verbose);
+
+
+%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%
+
+%% RESULTS DISPLAY
+
+if display
+
+[~,classif]=max(paramVBDS.param.piapost');
+K=paramVBDS.K;
+ns=paramVBDS.ds*Fs/1000;
+
+miny=min(min(Spike));
+maxy=max(max(Spike));
+
+% Clusters
+figure,
+title('Spike clusters')
+rootK=floor(sqrt(K))+1;
+rootK2=ceil(K/rootK);
+
+for k=1:K
+ subplot(rootK,rootK2,k);
+ piapostk=paramVBDS.param.piapost(:,k).*(1-paramVBDS.param.r);
+ Ik=find(piapostk>0.5);
+ plot((-ns/2+1:ns/2)/Fs*1000,Spike(Ik,:)');
+ hold on,
+ plot((-ns/2+1:ns/2)/Fs*1000,paramVBDS.S0(k,:),'r--','Linewidth',2);
+ plot((-ns/2+1:ns/2)/Fs*1000,paramVBDS.S0(k,:)+2.2*sqrt(paramVBDS.varskt(k)),'k--','Linewidth',1.5);
+ plot((-ns/2+1:ns/2)/Fs*1000,paramVBDS.S0(k,:)-2.2*sqrt(paramVBDS.varskt(k)),'k--','Linewidth',1.5);
+ ylim(1.05*[miny maxy]) %ylim(1.1*[-maxy maxy]);
+ xlim([-1 1]);
+ title(['#' num2str(length(Ik))],'FontSize',16)
+ xlabel('Time (ms)','FontSize',16)
+ ylabel('Voltage','FontSize',16)
+end
+
+% Original vs despikes signal (with spike instant as *)
+figure, plot((1:T)/Fs,xb,'Color',[1 1 1]/2); hold on, plot((1:T)/Fs,w,'r');
+plot(paramVBDS.ideltas/Fs,xb(paramVBDS.ideltas),'*g')
+legend({'raw data', 'despiked data',''});
+title('Despiking results with detected spikes as *')
+end
+
+
diff --git a/notchfilter.m b/notchfilter.m
new file mode 100644
index 0000000000000000000000000000000000000000..fc4a8fb1c50bfddc95a5e9072f78a2be75b1c374
--- /dev/null
+++ b/notchfilter.m
@@ -0,0 +1,61 @@
+function [filt] = notchfilter(dat,Fs,Fl,N)
+
+% NOTCHFILTER line noise reduction filter for EEG/MEG data
+%
+% [filt] = notchfilter(dat, Fsample, Fline)
+%
+% where
+% dat data matrix (Nchans X Ntime)
+% Fsample sampling frequency in Hz
+% Fline line noise frequency (would normally be 50Hz)
+% N optional filter order, default is 4
+%
+% if Fline is specified as 50, a band of 48-52 is filtered out
+% if Fline is specified as [low high], that band is filtered out
+
+% original (c) 2003, Pascal Fries
+% modifications (c) 2003, Robert Oostenveld
+%
+% This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
+% for the documentation and details.
+%
+% FieldTrip is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% FieldTrip is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with FieldTrip. If not, see <http://www.gnu.org/licenses/>.
+%
+% $Id: notchfilter.m 7123 2012-12-06 21:21:38Z roboos $
+
+if nargin<4
+ % set the default filter order
+ N = 1;
+end
+
+Nchans = size(dat,1);
+Nsamples = size(dat,2);
+
+if Nchans>Nsamples
+ dat=dat';
+end
+
+Nchans = size(dat,1);
+Nsamples = size(dat,2);
+
+% use a digital FIR filter
+Fn = Fs/2; % Nyquist frequency
+if length(Fl)==1
+ % default use a notch-width of 2Hz in both directions
+ % otherwise use the specified band
+ Fl = [Fl-0.5 Fl+0.5];
+end
+[B, A] = butter(N, [min(Fl)/Fn max(Fl)/Fn], 'stop');
+filt = filtfilt(B, A, dat')';
+
diff --git a/simulation_4.mat b/simulation_4.mat
new file mode 100755
index 0000000000000000000000000000000000000000..2a6de1ec63231aba8013b44133f2bc527039fe1b
Binary files /dev/null and b/simulation_4.mat differ
diff --git a/spikeExtractor.m b/spikeExtractor.m
new file mode 100644
index 0000000000000000000000000000000000000000..2a7b114d148edef3257dc45ad654e381bb397290
--- /dev/null
+++ b/spikeExtractor.m
@@ -0,0 +1,103 @@
+function [xb_spkband,spkband,spk,paramVBDS]=spikeExtractor(xb,Fs,paramVBDS,display,verb)
+%% Extract spike waveforms from data
+%% Reject overlapping events and high amplitude artifacts
+%
+% input:
+% xb: raw data
+% Fs: Frequency sammpling
+% paramVBDS: extraction parameters
+% display: plot of detected events and spikes if true
+%
+% output:
+% xb_spkband: filtered data in frequency band spike_band
+% spkband: matrix of band-pass filtered spike waveforms
+% spk: matrix of raw (unfiltered) spike waveforms
+% paramVBDS: updated extraction parameters (with detected max spike sample vector ideltas)
+
+[l,c]=size(xb);
+T=max(l,c);
+
+% band pass filtering
+[b,a]=butter(2,2*paramVBDS.spike_band/Fs);
+xb_spkband=filtfilt(b,a,xb);
+
+% low and high threshold
+thrQ=paramVBDS.coef*median(abs(xb_spkband)/0.6745);
+thrQH=paramVBDS.coefH*thrQ;
+ideltas=find(abs(xb_spkband)>thrQ);
+
+%spike duration in samples (ds ms)
+ns=max(Fs/1e3*paramVBDS.ds,2^5);
+
+
+% eliminate neighbouring sample
+[mini,imini]=min(diff(ideltas));
+while (mini<ns)
+ [~,isuppr]=min(abs(xb_spkband(ideltas(imini:imini+1))));
+ ideltas(imini+isuppr-1)=[];
+ [mini,imini]=min(diff(ideltas));
+end
+
+% to avoid border effect
+ideltas(ideltas<=ns)=[];
+ideltas(ideltas>=T-ns)=[];
+
+
+% eliminate positive or negative spikes if very low (arbitray set to <50)
+if sum(xb_spkband(ideltas)<0)<50
+ ideltas(xb_spkband(ideltas)<0)=[];
+end
+if sum(xb_spkband(ideltas)>0)<50
+ ideltas(xb_spkband(ideltas)>0)=[];
+end
+
+if display
+figure, plot(xb_spkband); hold,
+plot(ideltas,xb_spkband(ideltas),'r*');
+title('Filtered data with detected spikes as *')
+figure, plot(xb); hold,
+plot(ideltas,xb(ideltas),'r*');
+title('Raw data with detected spikes as *')
+end
+
+%%
+%extract spike waveforms
+
+Ns=length(ideltas);
+deltas=zeros(1,max([l c]));
+deltas(ideltas)=1;
+deltafen=conv(deltas,ones(ns,1));
+deltafen=deltafen(ns/2:end-ns/2);
+spkband=deltafen.*xb_spkband;
+spk=deltafen.*xb;
+
+iNZ=find(spkband~=0);
+spkband=spkband(iNZ);
+spk=spk(iNZ);
+spkband=reshape(spkband,ns,Ns)';
+spk=reshape(spk,ns,Ns)';
+
+% outliers filtering
+mspkband=max(abs(spkband'));
+spkband(mspkband>thrQH,:)=[];
+spk(mspkband>thrQH,:)=[];
+ideltas(mspkband>thrQH)=[];
+Ns=length(ideltas);
+deltas=zeros(max([l c]),1);
+deltas(ideltas)=1;
+
+% paramVBDS update with ideltas
+paramVBDS.ideltas=ideltas;
+
+if display
+figure, plot(spkband')
+title('Superimposed extracted spike waveforms (from filtered data)')
+end
+
+if verb
+ disp('####################')
+ disp('END OF SPIKE EXTRACTION')
+ disp(['Number of detected events : ',num2str(Ns)]);
+ disp('####################')
+ disp(' ');
+end
diff --git a/spikeRegistering.m b/spikeRegistering.m
new file mode 100644
index 0000000000000000000000000000000000000000..7c83a9cecac7dfbfaccc50ee3d4c796f75be0524
--- /dev/null
+++ b/spikeRegistering.m
@@ -0,0 +1,74 @@
+function [xbr,ideltasr,indsupp]=spikeRegistering(xb,xb_spkband,Fs,spkband,spk,ideltas,verb)
+%% spike registering (try to correct jitter due to sampling effect) and
+%% data reconstruction with the registered spikes (aligned with dwt - level 5)
+%
+% input:
+% xb: raw data
+% xb_spkband: filtered data (spike frequency band)
+% Fs: Frequency sammpling
+% spkband: matrix of band-pass filtered spike waveforms
+% spk: matrix of raw (unfiltered) spike waveforms
+% ideltas: maximum time samples of detected spikes in xb
+% Ns: number of detected events (length of ideltas)
+% ds: spike duration in ms
+% verb: verbose mode
+%
+% output:
+% xbr: raw data with spike waveforms registered
+% xb_spkbandr: band-pass filtered data with spike waveforms registered
+% idletasr: instant of registered spikes in xbr/xb_spkbandr
+
+
+
+
+
+% spike registering by upsampling (factor 10)
+
+Ns=length(ideltas);
+T=length(xb);
+vecun=ones(1,20);
+ns=max(Fs/1e3*3,2^5);
+ups=10;
+for ii=1:Ns
+ y=interp(spkband(ii,:),ups);
+ [~,im]=max(abs(y(ns*ups/2-2*ups:ns*ups/2+2*ups)));
+ difm=2*ups-im+1;
+ yr=[y(1)*vecun(1:difm) y(max(1,-difm+1):min(ns*ups,ns*ups-difm)) y(end)*vecun(1:-difm)];
+ spkband(ii,:)=yr(ups:ups:end);
+
+ y=interp(spk(ii,:),ups);
+ yr=[y(1)*vecun(1:difm) y(max(1,-difm+1):min(ns*ups,ns*ups-difm)) y(end)*vecun(1:-difm)];
+ spk(ii,:)=yr(ups:ups:end);
+end
+
+% data rebuilding, preparing for wavelet feature extraction...
+xb_spkbandr=xb_spkband;
+xbr=xb;
+ideltass=[0 ideltas];
+ideltasr=[0 ideltas];
+lsupp=zeros(Ns,1);
+indsupp=[];
+tm=zeros(Ns,1);
+for i=2:Ns+1
+ rdec=rem(ideltasr(i),2^5);
+ lsupp(i-1)=rdec;
+ tm(i-1)=floor((ideltass(i)+ideltass(i-1))/2-rdec/2);
+ ideltasr(i:end)=ideltasr(i:end)-rdec;
+ indsupp=[indsupp tm(i-1)+1:tm(i-1)+rdec];
+ xbr(ideltass(i)-ns/2+1:ideltass(i)+ns/2)=spk(i-1,:);
+ xb_spkbandr(ideltass(i)-ns/2+1:ideltass(i)+ns/2)=spkband(i-1,:);
+end
+if mod(length(indsupp),2)==1
+ indsupp=[indsupp T];
+end
+xb_spkbandr(indsupp)=[];
+xbr(indsupp)=[];
+ideltasr=ideltasr(2:end);
+
+if verb
+ disp(' ')
+ disp('####################')
+ disp('DATA READY FOR FEATURE EXTRACTION')
+ disp('####################')
+ disp(' ')
+end