wolffd@0: function c = mirclassify(a,da,t,dt,varargin)
wolffd@0: %   c = mirclassify(test,features_test,train,features_train) classifies the
wolffd@0: %       audio sequence(s) contained in the audio object test, along the
wolffd@0: %       analytic feature(s) features_test, following the supervised
wolffd@0: %       learning of a training set defined by the audio object train and
wolffd@0: %       the corresponding analytic feature(s) features_train.
wolffd@0: %           * The analytic feature(s) features_test should *not* be frame 
wolffd@0: %               decomposed. Frame-decomposed data should first be 
wolffd@0: %               summarized, using for instance mirmean or mirstd.
wolffd@0: %           * Multiple analytic features have to be grouped into one array 
wolffd@0: %               of cells.
wolffd@0: %       You can also integrate your own arrays of numbers computed outside
wolffd@0: %           MIRtoolbox as part of the features. These arrays should be 
wolffd@0: %           given as matrices where each successive column is the analysis 
wolffd@0: %           of each successive file.
wolffd@0: %   Example:
wolffd@0: %       mirclassify(test, mfcc(test), train, mfcc(train))
wolffd@0: %       mirclassify(test, {mfcc(test), centroid(test)}, ...
wolffd@0: %                train, {mfcc(train), centroid(train)})
wolffd@0: %   Optional argument:
wolffd@0: %       mirclassify(...,'Nearest') uses the minimum distance strategy.
wolffd@0: %           (by default)
wolffd@0: %       mirclassify(...,'Nearest',k) uses the k-nearest-neighbour strategy.
wolffd@0: %           Default value: k = 1, corresponding to the minimum distance
wolffd@0: %               strategy.
wolffd@0: %       mirclassify(...,'GMM',ng) uses a gaussian mixture model. Each class is
wolffd@0: %           modeled by at most ng gaussians.
wolffd@0: %           Default value: ng = 1.
wolffd@0: %           Additionnally, the type of mixture model can be specified,
wolffd@0: %           using the set of value proposed in the gmm function: i.e.,
wolffd@0: %           'spherical','diag','full' (default value) and 'ppca'.
wolffd@0: %               (cf. help gmm)
wolffd@0: %           Requires the Netlab toolbox.
wolffd@0: 
wolffd@0: lab = get(t,'Label');
wolffd@0: c.labtraining = lab;
wolffd@0: rlab = get(a,'Label');
wolffd@0: c.labtest = rlab;
wolffd@0: [k,ncentres,covartype,kmiter,emiter,d,norml,mahl] = scanargin(varargin);
wolffd@0: disp('Classifying...')
wolffd@0: if not(iscell(dt))
wolffd@0:     dt = {dt};
wolffd@0: end
wolffd@0: lvt = length(get(t,'Data'));
wolffd@0: vt = [];
wolffd@0: for i = 1:length(dt)
wolffd@0:     if isnumeric(dt{i})
wolffd@0:         d = cell(1,size(dt{i},2));
wolffd@0:         for j = 1:size(dt{i},2)
wolffd@0:             d{j} = dt{i}(:,j);
wolffd@0:         end
wolffd@0:     else
wolffd@0:         d = get(dt{i},'Data');
wolffd@0:     end
wolffd@0:     vt = integrate(vt,d,lvt,norml);
wolffd@0:     if isa(dt{i},'scalar')
wolffd@0:         m = mode(dt{i});
wolffd@0:         if not(isempty(m))
wolffd@0:             vt = integrate(vt,m,lvt,norml);
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: c.training = vt;
wolffd@0: dim = size(vt,1);
wolffd@0: if not(iscell(da))
wolffd@0:     da = {da};
wolffd@0: end
wolffd@0: lva = length(get(a,'Data'));
wolffd@0: va = [];
wolffd@0: for i = 1:length(da)
wolffd@0:     if isnumeric(da{i})
wolffd@0:         d = cell(1,size(da{i},2));
wolffd@0:         for j = 1:size(da{i},2)
wolffd@0:             d{j} = da{i}(:,j);
wolffd@0:         end
wolffd@0:     else
wolffd@0:         d = get(da{i},'Data');
wolffd@0:     end
wolffd@0:     va = integrate(va,d,lva,norml);
wolffd@0:     if isa(da{i},'scalar')
wolffd@0:         m = mode(da{i});
wolffd@0:         if not(isempty(m))
wolffd@0:             va = integrate(va,m,lva,norml);
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: c.test = va;
wolffd@0: c.nbobs = lvt;
wolffd@0: totva = [vt va];
wolffd@0: mahl = cov(totva');
wolffd@0: if k                % k-Nearest Neighbour
wolffd@0:     c.nbparam = lvt;
wolffd@0:     for l = 1:lva
wolffd@0:         [sv,idx] = sort(distance(va(:,l),vt,d,mahl));
wolffd@0:         labs = cell(0); % Class labels
wolffd@0:         founds = [];    % Number of found elements in each class
wolffd@0:         for i = idx(1:k)
wolffd@0:             labi = lab{i};
wolffd@0:             found = 0;
wolffd@0:             for j = 1:length(labs)
wolffd@0:                 if isequal(labi,labs{j})
wolffd@0:                     found = j;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             if found
wolffd@0:                 founds(found) = founds(found)+1;
wolffd@0:             else
wolffd@0:                 labs{end+1} = labi;
wolffd@0:                 founds(end+1) = 1;
wolffd@0:             end
wolffd@0:         end
wolffd@0:         [b ib] = max(founds);
wolffd@0:         c.classes{l} = labs{ib};
wolffd@0:     end
wolffd@0: elseif ncentres     % Gaussian Mixture Model
wolffd@0:     labs = cell(0);    % Class labels
wolffd@0:     founds = cell(0);  % Elements associated to each label.
wolffd@0:     for i = 1:lvt
wolffd@0:         labi = lab{i};
wolffd@0:         found = 0;
wolffd@0:         for j = 1:length(labs)
wolffd@0:             if isequal(labi,labs{j})
wolffd@0:                 founds{j}(end+1) = i;
wolffd@0:                 found = 1;
wolffd@0:             end
wolffd@0:         end
wolffd@0:         if not(found)
wolffd@0:             labs{end+1} = labi;
wolffd@0:             founds{end+1} = i;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     options      = zeros(1, 18);
wolffd@0:     options(2:3) = 1e-4;
wolffd@0:     options(4)   = 1e-6;
wolffd@0:     options(16)  = 1e-8;
wolffd@0:     options(17)  = 0.1;
wolffd@0:     options(1)   = 0; %Prints out error values, -1 else
wolffd@0:     c.nbparam = 0;
wolffd@0:     OK = 0;
wolffd@0:     while not(OK)
wolffd@0:         OK = 1;
wolffd@0:         for i = 1:length(labs)
wolffd@0:             options(14)  = kmiter;
wolffd@0:             try
wolffd@0:                 mix{i} = gmm(dim,ncentres,covartype);
wolffd@0:             catch
wolffd@0:                 error('ERROR IN CLASSIFY: Netlab toolbox not installed.');
wolffd@0:             end
wolffd@0:             mix{i} = netlabgmminit(mix{i},vt(:,founds{i})',options);
wolffd@0:             options(5)   = 1;
wolffd@0:             options(14)  = emiter;
wolffd@0:             try
wolffd@0:                 mix{i} = gmmem(mix{i},vt(:,founds{i})',options);
wolffd@0:                 c.nbparam = c.nbparam + ...
wolffd@0:                     length(mix{i}.centres(:)) + length(mix{i}.covars(:));
wolffd@0:             catch
wolffd@0:                 err = lasterr;
wolffd@0:                 warning('WARNING IN CLASSIFY: Problem when calling GMMEM:');
wolffd@0:                 disp(err);
wolffd@0:                 disp('Let us try again...');
wolffd@0:                 OK = 0;
wolffd@0:             end
wolffd@0:         end    
wolffd@0:     end
wolffd@0:     pr = zeros(lva,length(labs));
wolffd@0:     for i = 1:length(labs)
wolffd@0:         prior = length(founds{i})/lvt;
wolffd@0:         pr(:,i) = prior * gmmprob(mix{i},va');
wolffd@0:         %c.post{i} = gmmpost(mix{i},va');
wolffd@0:     end
wolffd@0:     [mm ib] = max(pr');
wolffd@0:     for i = 1:lva
wolffd@0:         c.classes{i} = labs{ib(i)};
wolffd@0:     end
wolffd@0: end
wolffd@0: if isempty(rlab)
wolffd@0:     c.correct = NaN;
wolffd@0: else
wolffd@0:     correct = 0;
wolffd@0:     for i = 1:lva
wolffd@0:         if isequal(c.classes{i},rlab{i})
wolffd@0:             correct = correct + 1;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     c.correct = correct / lva;
wolffd@0: end        
wolffd@0: c = class(c,'mirclassify');
wolffd@0: 
wolffd@0: 
wolffd@0: function vt = integrate(vt,v,lvt,norml)
wolffd@0: vtl = [];
wolffd@0: for l = 1:lvt
wolffd@0:     vl = v{l};
wolffd@0:     if iscell(vl)
wolffd@0:         vl = vl{1};
wolffd@0:     end
wolffd@0:     if iscell(vl)
wolffd@0:         vl = vl{1};
wolffd@0:     end
wolffd@0:     if size(vl,2) > 1
wolffd@0:         mirerror('MIRCLASSIFY','The analytic features guiding the classification should not be frame-decomposed.');
wolffd@0:     end
wolffd@0:     vtl(:,l) = vl;
wolffd@0: end
wolffd@0: if norml
wolffd@0:     dnom = repmat(std(vtl,0,2),[1 size(vtl,2)]);
wolffd@0:     dnom = dnom + (dnom == 0);  % In order to avoid division by 0
wolffd@0:     vtl = (vtl - repmat(mean(vtl,2),[1 size(vtl,2)])) ./ dnom;
wolffd@0: end
wolffd@0: vt(end+1:end+size(vtl,1),:) = vtl;
wolffd@0: 
wolffd@0: 
wolffd@0: function [k,ncentres,covartype,kmiter,emiter,d,norml,mahl] = scanargin(v)
wolffd@0: k = 1;
wolffd@0: d = 0;
wolffd@0: i = 1;
wolffd@0: ncentres = 0;
wolffd@0: covartype = 'full';
wolffd@0: kmiter = 10;
wolffd@0: emiter = 100;
wolffd@0: norml = 1;
wolffd@0: mahl = 1;
wolffd@0: while i <= length(v)
wolffd@0:     arg = v{i};
wolffd@0:     if ischar(arg) && strcmpi(arg,'Nearest')
wolffd@0:         k = 1;
wolffd@0:         if length(v)>i && isnumeric(v{i+1})
wolffd@0:             i = i+1;
wolffd@0:             k = v{i};
wolffd@0:         end
wolffd@0:     elseif ischar(arg) && strcmpi(arg,'GMM')
wolffd@0:         k = 0;
wolffd@0:         ncentres = 1;
wolffd@0:         if length(v)>i
wolffd@0:             if isnumeric(v{i+1})
wolffd@0:                 i = i+1;
wolffd@0:                 ncentres = v{i};
wolffd@0:                 if length(v)>i && ischar(v{i+1})
wolffd@0:                     i = i+1;
wolffd@0:                     covartype = v{i};
wolffd@0:                 end
wolffd@0:             elseif ischar(v{i+1})
wolffd@0:                 i = i+1;
wolffd@0:                 covartype = v{i};
wolffd@0:                 if length(v)>i && isnumeric(v{i+1})
wolffd@0:                     i = i+1;
wolffd@0:                     ncentres = v{i};
wolffd@0:                 end
wolffd@0:             end                
wolffd@0:         end
wolffd@0:     elseif isnumeric(arg)
wolffd@0:         k = v{i};
wolffd@0:     else
wolffd@0:         error('ERROR IN MIRCLASSIFY: Syntax error. See help mirclassify.');
wolffd@0:     end    
wolffd@0:     i = i+1;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function y = distance(a,t,d,mahl)
wolffd@0: 
wolffd@0: for i = 1:size(t,2)
wolffd@0:     if det(mahl) > 0  % more generally, uses cond
wolffd@0:         lham = inv(mahl);
wolffd@0:     else
wolffd@0:         lham = pinv(mahl);
wolffd@0:     end
wolffd@0:     y(i) = sqrt((a - t(:,i))'*lham*(a - t(:,i)));        
wolffd@0: end
wolffd@0: %y = sqrt(sum(repmat(a,[1,size(t,2)])-t,1).^2);