wolffd@0: function d = mirdist(x,y,dist)
wolffd@0: %   d = mirdist(x,y) evaluates the distance between x and y.
wolffd@0: %   x is the feature values corresponding to one audio file, and y is the
wolffd@0: %       values (for the same feature) corrdponding to one (or several)
wolffd@0: %       audio files.
wolffd@0: %   If x and y are not decomposed into frames,
wolffd@0: %       d = mirdist(x,y,f) specifies distance function.
wolffd@0: %           Default value: f = 'Cosine'
wolffd@0: %   If x and y are composed of clustered frames (using mircluster), the
wolffd@0: %       cluster signatures are compared using Earth Mover Distance.
wolffd@0: %       (Logan, Salomon, 2001)
wolffd@0: %   If x and y contains peaks, the vectors representing the peak
wolffd@0: %       distributions are compared using Euclidean distance.
wolffd@0: %       (used with mirnovelty in Jacobson, 2006)
wolffd@0: %
wolffd@0: % The Earth Mover Distance is based on the implementation by Yossi Rubner,
wolffd@0: % wrapped for Matlab by Elias Pampalk.
wolffd@0: 
wolffd@0: if not(isa(x,'mirdata'))
wolffd@0:     x = miraudio(x);
wolffd@0: end
wolffd@0: if not(isa(y,'mirdata'))
wolffd@0:     y = miraudio(y);
wolffd@0: end
wolffd@0: 
wolffd@0: clx = get(x,'Clusters');
wolffd@0: if isempty(clx{1})
wolffd@0:     px = get(x,'PeakPos');
wolffd@0:     if not(iscell(px)) || isempty(px{1}) || ...
wolffd@0:             not(iscell(px{1})) || isempty(px{1}{1}) || not(iscell(px{1}{1}))
wolffd@0:         if nargin < 3
wolffd@0:             dist = 'Cosine';
wolffd@0:         end
wolffd@0: 
wolffd@0:         d = get(x,'Data');
wolffd@0:         dd = d{1}{1};
wolffd@0:         if iscell(dd)
wolffd@0:             dd = dd{1};
wolffd@0:         end
wolffd@0:         if size(dd,2)>1
wolffd@0:             if size(dd,1)>1
wolffd@0:                 error('ERROR IN MIRDIST: If the input is decomposed into frames, they should first be clustered.'); 
wolffd@0:             else
wolffd@0:                 dd = dd';
wolffd@0:             end
wolffd@0:         end
wolffd@0: 
wolffd@0:         e = get(y,'Data');
wolffd@0:         dt = cell(1,length(e));
wolffd@0:         for h = 1:length(e)
wolffd@0:             ee = e{h}{1};
wolffd@0:             if iscell(ee)
wolffd@0:                 ee = ee{1};
wolffd@0:             end
wolffd@0:             if size(ee,2)>1
wolffd@0:                 if size(ee,1)>1
wolffd@0:                     error('ERROR IN MIRDIST: If the input is decomposed into frames, they should first be clustered.'); 
wolffd@0:                 else
wolffd@0:                     ee = ee';
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             if isempty(ee)
wolffd@0:                 if isempty(dd)
wolffd@0:                     dt{h}{1} = 0;
wolffd@0:                 else
wolffd@0:                     dt{h}{1} = Inf;
wolffd@0:                 end
wolffd@0:             else
wolffd@0:                 if length(dd)<length(ee)
wolffd@0:                     dd(length(ee)) = 0;
wolffd@0:                     %ee = ee(1:length(d));
wolffd@0:                 elseif length(ee)<length(dd)
wolffd@0:                     ee(length(dd)) = 0;
wolffd@0:                     %dd = dd(1:length(ee));
wolffd@0:                 end
wolffd@0:                 if length(dd) == 1
wolffd@0:                     dt{h}{1} = abs(dd-ee);
wolffd@0:                 elseif norm(dd) && norm(ee)
wolffd@0:                     dt{h}{1} = pdist([dd(:)';ee(:)'],dist);
wolffd@0:                 else
wolffd@0:                     dt{h}{1} = NaN;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     else
wolffd@0:         % Euclidean distance between vectors to compare data with peaks
wolffd@0:         % (used with mirnovelty in Jacobson, 2006).
wolffd@0:         sig = pi/4;
wolffd@0:         dx = get(x,'Data');
wolffd@0:         nx = length(px{1}{1}{1});
wolffd@0:         cx = mean(px{1}{1}{1}/length(dx{1}{1}));
wolffd@0:         dy = get(y,'Data');
wolffd@0:         py = get(y,'PeakPos');
wolffd@0:         dt = cell(1,length(py));
wolffd@0:         for h = 1:length(py)
wolffd@0:             ny = length(py{h}{1}{1});
wolffd@0:             cy = mean(py{h}{1}{1}/length(dy{h}{1}));
wolffd@0:             dt{h}{1} = sqrt((nx*cos(sig*cx)-ny*cos(sig*cy))^2 ...
wolffd@0:                            +(nx*sin(sig*cx)-ny*sin(sig*cy))^2);
wolffd@0:         end
wolffd@0:     end
wolffd@0: else
wolffd@0:     % Earth Mover's Distance to compare clustered data.
wolffd@0:     cly = get(y,'Clusters');
wolffd@0:     dt = cell(1,length(cly));
wolffd@0:     for h = 1:length(cly)
wolffd@0:         cost = zeros(length(clx{1}.weight),length(cly{h}.weight));
wolffd@0:         for i = 1:length(clx{1}.weight)
wolffd@0:             for j = 1:length(cly{h}.weight)
wolffd@0:                 covx = clx{1}.covar(:,i);
wolffd@0:                 covy = cly{h}.covar(:,j);
wolffd@0:                 mux = clx{1}.centr(:,i);
wolffd@0:                 muy = cly{h}.centr(:,j);
wolffd@0:                 cost(i,j) = sum(covx./covy + covy./covx + ...
wolffd@0:                         (mux-muy).^2.*(1./covx + 1./covy) - 2);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         dt{h}{1} = emd_wrapper(cost,clx{1}.weight,cly{h}.weight);
wolffd@0:     end
wolffd@0: end
wolffd@0: d = mirscalar(y,'Data',dt,'Title',[get(y,'Title'),' Distance'],...
wolffd@0:                 'Name',get(x,'Name'),'Name2',get(y,'Name'));