wolffd@0: function varargout = mirsimatrix(orig,varargin)
wolffd@0: %   m = mirsimatrix(x) computes the similarity matrix resulting from the 
wolffd@0: %       mutual comparison between each possible frame analysis in x.
wolffd@0: %           By default, x is the spectrum of the frame decomposition.
wolffd@0: %           But it can be any other frame analysis.
wolffd@0: %   Optional argument:
wolffd@0: %       mirsimatrix(...,'Distance',f) specifies the name of a distance
wolffd@0: %           function, from those proposed in the Statistics Toolbox
wolffd@0: %               (help pdist).
wolffd@0: %           default value: f = 'cosine'
wolffd@0: %       mirsimatrix(...,'Dissimilarity') return the dissimilarity matrix,
wolffd@0: %           which is the intermediary result before the computation of the 
wolffd@0: %           actual similarity matrix. It shows the distance between each
wolffd@0: %           possible frame analysis in x.
wolffd@0: %       mirsimatrix(...,'Similarity',f) indicates the function f specifying
wolffd@0: %           the way the distance values in the dissimilarity matrix are
wolffd@0: %           transformed into similarity values.
wolffd@0: %           Possible values:
wolffd@0: %               f = 'oneminus' (default value) 
wolffd@0: %                   corresponding to f(x) = 1-x
wolffd@0: %               f = 'exponential'
wolffd@0: %                   corresponding to f(x)= exp(-x)
wolffd@0: %       mirsimatrix(...,'Width',w) or more simply dissimatrix(...,w) 
wolffd@0: %           specifies the size of the diagonal bandwidth, in samples,
wolffd@0: %           outside which the dissimilarity will not be computed.
wolffd@0: %           if w = inf (default value), all the matrix will be computed.
wolffd@0: %       mirsimatrix(...,'Horizontal') rotates the matrix 45 degrees in order to
wolffd@0: %           make the first diagonal horizontal, and to restrict on the
wolffd@0: %           diagonal bandwidth only.
wolffd@0: %       mirsimatrix(...,'TimeLag') transforms the (non-rotated) matrix into
wolffd@0: %           a time-lag matrix, making the first diagonal horizontal as well
wolffd@0: %           (corresponding to the zero-lag line).
wolffd@0: %
wolffd@0: %       mirsimatrix(M,r) creates a mirsimatrix similarity matrix based on
wolffd@0: %           the Matlab square matrix M, of frame rate r (in Hz.)
wolffd@0: %               By default r = 20 Hz.
wolffd@0: %           mirsimatrix(M,r,'Dissimilarity') creates instead a mirsimatrix
wolffd@0: %               dissimilarity matrix.
wolffd@0: %
wolffd@0: %   Foote, J. & Cooper, M. (2003). Media Segmentation using Self-Similarity
wolffd@0: % Decomposition,. In Proc. SPIE Storage and Retrieval for Multimedia
wolffd@0: % Databases, Vol. 5021, pp. 167-75.
wolffd@0: 
wolffd@0: %       mirsimatrix(...,'Filter',10) filter along the diagonal of the matrix,
wolffd@0: %           using a uniform moving average filter of size f. The result is
wolffd@0: %           represented in a time-lag matrix.
wolffd@0: 
wolffd@0: 
wolffd@0:         distance.key = 'Distance';
wolffd@0:         distance.type = 'String';
wolffd@0:         distance.default = 'cosine';
wolffd@0:     option.distance = distance;
wolffd@0: 
wolffd@0:         simf.key = 'Similarity';
wolffd@0:         simf.type = 'String';
wolffd@0:         simf.default = 'oneminus';
wolffd@0:         simf.keydefault = 'Similarity';        
wolffd@0:         simf.when = 'After';
wolffd@0:     option.simf = simf;
wolffd@0: 
wolffd@0:         dissim.key = 'Dissimilarity';
wolffd@0:         dissim.type = 'Boolean';
wolffd@0:         dissim.default = 0;
wolffd@0:     option.dissim = dissim;
wolffd@0:     
wolffd@0:         K.key = 'Width';
wolffd@0:         K.type = 'Integer';
wolffd@0:         K.default = Inf;
wolffd@0:     option.K = K;
wolffd@0:     
wolffd@0:         filt.key = 'Filter';
wolffd@0:         filt.type = 'Integer';
wolffd@0:         filt.default = 0;
wolffd@0:         filt.when = 'After';
wolffd@0:     option.filt = filt;
wolffd@0: 
wolffd@0:         view.type = 'String';
wolffd@0:         view.default = 'Standard';
wolffd@0:         view.choice = {'Standard','Horizontal','TimeLag'};
wolffd@0:         view.when = 'After';
wolffd@0:     option.view = view;
wolffd@0:     
wolffd@0:         rate.type = 'Integer';
wolffd@0:         rate.position = 2;
wolffd@0:         rate.default = 20;
wolffd@0:     option.rate = rate;
wolffd@0:     
wolffd@0: specif.option = option;
wolffd@0: specif.nochunk = 1;
wolffd@0: varargout = mirfunction(@mirsimatrix,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: if isnumeric(x)
wolffd@0:     m.diagwidth = Inf;
wolffd@0:     m.view = 's';
wolffd@0:     m.similarity = NaN;
wolffd@0:     m.graph = {};
wolffd@0:     m.branch = {};
wolffd@0:     m = class(m,'mirsimatrix',mirdata);
wolffd@0:     m = set(m,'Title','Dissimilarity matrix');
wolffd@0:     fp = repmat(((1:size(x,1))-.5)/option.rate,[2,1]);
wolffd@0:     x = set(m,'Data',{x},'Pos',[],...
wolffd@0:               'FramePos',{{fp}},'Name',{inputname(1)});
wolffd@0: else
wolffd@0:     if not(isamir(x,'mirsimatrix'))
wolffd@0:         if (isamir(x,'miraudio'))
wolffd@0:             if isframed(x)
wolffd@0:                 x = mirspectrum(x);
wolffd@0:             else
wolffd@0:                 x = mirspectrum(x,'Frame',0.05,1);
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: type = 'mirsimatrix';
wolffd@0: 
wolffd@0: 
wolffd@0: function m = main(orig,option,postoption)
wolffd@0: if iscell(orig)
wolffd@0:     orig = orig{1};
wolffd@0: end
wolffd@0: if isa(orig,'mirsimatrix')
wolffd@0:     d = get(orig,'Data');
wolffd@0:     for k = 1:length(d)
wolffd@0:         if iscell(option) && isfield(option,'K') && option.K < orig.diagwidth
wolffd@0:             nl = size(d{k},1);
wolffd@0:             if strcmp(orig.view,'h')
wolffd@0:                 dl = (nl - option.K)/2;
wolffd@0:                 dk = d{k}(ceil(dl):nl-floor(dl),:);
wolffd@0:             else
wolffd@0:                 [spA,spd] = spdiags(d{k},-ceil(option.K/2):ceil(option.K/2));
wolffd@0:                 dk = full(spdiags(spA,spd,nl,size(d{k},2)));
wolffd@0:             end
wolffd@0:             d{k} = dk;
wolffd@0:             orig.diagwidth = option.K;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     m = set(orig,'Data',d);
wolffd@0: else
wolffd@0:     v = get(orig,'Data');
wolffd@0:     d = cell(1,length(v));
wolffd@0:     for k = 1:length(v)
wolffd@0:         vk = v{k};
wolffd@0:         if mirwaitbar
wolffd@0:             handle = waitbar(0,'Computing dissimilarity matrix...');
wolffd@0:         else
wolffd@0:             handle = 0;
wolffd@0:         end
wolffd@0:         if 0 %iscell(vk)
wolffd@0:             try
wolffd@0:                 vk = cell2mat(vk);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         if 0 %iscell(vk) %&& length(vk) > 1 %%% ATTENTION KL!!<<<<<<<<<<<<
wolffd@0:             l = length(vk);
wolffd@0:             dk = NaN(l,l);
wolffd@0:             hK = floor(option.K/2);
wolffd@0:             for i = 1:l
wolffd@0:                 if handle
wolffd@0:                     waitbar(i/l,handle);
wolffd@0:                 end
wolffd@0:                 for j = max(1,i-hK):min(l,i+hK)
wolffd@0:                     dk(i,j) = KL(vk{i},vk{j});
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         else
wolffd@0:             if not(iscell(vk))
wolffd@0:                 vk = {vk};
wolffd@0:             end    
wolffd@0:             for z = 1:length(vk)
wolffd@0:                 vz = vk{z};
wolffd@0:                 ll = size(vz,1);
wolffd@0:                 l = size(vz,2);
wolffd@0:                 nc = size(vz,3);
wolffd@0:                 if ll==1 && nc>1
wolffd@0:                     vz = squeeze(vz)';
wolffd@0:                     ll = nc;
wolffd@0:                     nc = 1;
wolffd@0:                 end
wolffd@0:                 nd = size(vz,4);
wolffd@0:                 if not(isempty(postoption)) && ...
wolffd@0:                     strcmpi(postoption.view,'TimeLag')
wolffd@0:                     lK = ceil(option.K/2);
wolffd@0:                     if isinf(lK)
wolffd@0:                         lK = l;
wolffd@0:                     end
wolffd@0:                     dk{z} = NaN(lK,l,nc);
wolffd@0:                 else
wolffd@0:                     dk{z} = NaN(l,l,nc);
wolffd@0:                 end
wolffd@0:                 for g = 1:nc
wolffd@0:                     if nd == 1
wolffd@0:                         vv = vz;
wolffd@0:                     else
wolffd@0:                         vv = zeros(ll*nd,l);
wolffd@0:                         for h = 1:nd
wolffd@0:                             if iscell(vz)
wolffd@0:                                 for i = 1:ll
wolffd@0:                                     for j = 1:l
wolffd@0:                                         vj = vz{i,j,g,h};
wolffd@0:                                         if isempty(vj)
wolffd@0:                                             vv((h-1)*ll+i,j) = NaN;
wolffd@0:                                         else
wolffd@0:                                             vv((h-1)*ll+i,j) = vj;
wolffd@0:                                         end
wolffd@0:                                     end
wolffd@0:                                 end
wolffd@0:                             else
wolffd@0:                                 tic
wolffd@0:                                 vv((h-1)*ll+1:h*ll,:) = vz(:,:,g,h);
wolffd@0:                                 toc
wolffd@0:                             end
wolffd@0:                         end
wolffd@0:                     end
wolffd@0:                     if isinf(option.K) && not(strcmpi(postoption.view,'TimeLag'))
wolffd@0:                         try
wolffd@0:                             manually = 0;
wolffd@0:                             dk{z}(:,:,g) = squareform(pdist(vv',option.distance));
wolffd@0:                             if option.K < Inf
wolffd@0:                                 [spA,spd] = spdiags(dk{z},...
wolffd@0:                                     -ceil(option.K/2):ceil(option.K/2));
wolffd@0:                                 dk{z}(:,:,g) = full(spdiags(spA,spd,size(dk,1),size(dk{z},2)));
wolffd@0:                             end
wolffd@0:                         catch
wolffd@0:                             %err = lasterror;
wolffd@0:                             %warning(err.message)
wolffd@0:                             %disp('Statistics Toolbox does not seem to be
wolffd@0:                             %installed. Recompute the distance matrix manually.');
wolffd@0:                             manually = 1;
wolffd@0:                         end
wolffd@0:                     else
wolffd@0:                         manually = 1;
wolffd@0:                     end
wolffd@0:                     if manually
wolffd@0:                         disf = str2func(option.distance);
wolffd@0:                         if strcmpi(option.distance,'cosine')
wolffd@0:                             for i = 1:l
wolffd@0:                                 vv(:,i) = vv(:,i)/norm(vv(:,i));
wolffd@0:                             end
wolffd@0:                         end
wolffd@0:                         if not(isempty(postoption)) && ...
wolffd@0:                                 strcmpi(postoption.view,'TimeLag')
wolffd@0:                             lK = ceil(option.K/2);
wolffd@0:                             for i = 1:l
wolffd@0:                                 if mirwaitbar && (mod(i,100) == 1 || i == l)
wolffd@0:                                     waitbar(i/l,handle);
wolffd@0:                                 end
wolffd@0:                                 ij = min(i+lK-1,l);
wolffd@0:                                 dkij = disf(vv(:,i),vv(:,i:ij));
wolffd@0:                                 for j = 1:ij-i+1
wolffd@0:                                     dk{z}(j,i+j-1,g) = dkij(j);
wolffd@0:                                 end
wolffd@0:                             end
wolffd@0:                         else
wolffd@0:                             for i = 1:l
wolffd@0:                                 if mirwaitbar && (mod(i,100) == 1 || i == l)
wolffd@0:                                     waitbar(i/l,handle);
wolffd@0:                                 end
wolffd@0:                                 j = min(i+option.K-1,l);
wolffd@0:                                 dkij = disf(vv(:,i),vv(:,i:j));
wolffd@0:                                 dk{z}(i,i:j,g) = dkij;
wolffd@0:                                 dk{z}(i:j,i,g) = dkij';
wolffd@0:                             end
wolffd@0:                         end
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:         d{k} = dk;
wolffd@0:         if handle
wolffd@0:             delete(handle)
wolffd@0:         end
wolffd@0:     end
wolffd@0:     m.diagwidth = option.K;
wolffd@0:     if not(isempty(postoption)) && strcmpi(postoption.view,'TimeLag')
wolffd@0:         m.view = 'l';
wolffd@0:     else
wolffd@0:         m.view = 's';
wolffd@0:     end
wolffd@0:     m.similarity = 0;
wolffd@0:     m.graph = {};
wolffd@0:     m.branch = {};
wolffd@0:     m = class(m,'mirsimatrix',mirdata(orig));
wolffd@0:     m = purgedata(m);
wolffd@0:     m = set(m,'Title','Dissimilarity matrix');
wolffd@0:     m = set(m,'Data',d,'Pos',[]);
wolffd@0: end
wolffd@0: if not(isempty(postoption))
wolffd@0:     if strcmpi(m.view,'s')
wolffd@0:         if strcmpi(postoption.view,'Horizontal')
wolffd@0:             for k = 1:length(d)
wolffd@0:                 for z = 1:length(d{k})
wolffd@0:                     d{k}{z} = rotatesim(d{k}{z},m.diagwidth);
wolffd@0:                     if option.K < m.diagwidth
wolffd@0:                         W = size(d{k}{z},1);
wolffd@0:                         hW = ceil(W/2);
wolffd@0:                         hK = ceil(option.K/2);
wolffd@0:                         d{k}{z} = d{k}{z}(hW-hK:hW+hK,:);
wolffd@0:                         m.diagwidth = option.K;
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m = set(m,'Data',d);
wolffd@0:             m.view = 'h';
wolffd@0:         elseif strcmpi(postoption.view,'TimeLag') || postoption.filt
wolffd@0:             W = ceil(m.diagwidth/2);
wolffd@0:             for k = 1:length(d)
wolffd@0:                 for z = 1:length(d{k})
wolffd@0:                     if isinf(W)
wolffd@0:                         dz = NaN(size(d{k}{z}));
wolffd@0:                     else
wolffd@0:                         dz = NaN(W,size(d{k}{z},2));
wolffd@0:                     end
wolffd@0:                     for l = 1:size(dz,1)
wolffd@0:                         dz(l,l:end) = diag(d{k}{z},l-1)';
wolffd@0:                     end
wolffd@0:                     if option.K < m.diagwidth
wolffd@0:                         dz = dz(1:ceil(option.K/2),:);
wolffd@0:                         m.diagwidth = option.K;
wolffd@0:                     end
wolffd@0:                     d{k}{z}= dz;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m = set(m,'Data',d);
wolffd@0:             m.view = 'l';
wolffd@0:         end
wolffd@0:     end
wolffd@0:     if ischar(postoption.simf)
wolffd@0:         if strcmpi(postoption.simf,'Similarity')
wolffd@0:             if not(isequal(m.similarity,NaN))
wolffd@0:                 option.dissim = 0;
wolffd@0:             end
wolffd@0:             postoption.simf = 'oneminus';
wolffd@0:         end
wolffd@0:         if isequal(m.similarity,0) && isstruct(option) ...
wolffd@0:                 && isfield(option,'dissim') && not(option.dissim)
wolffd@0:             simf = str2func(postoption.simf);
wolffd@0:             for k = 1:length(d)
wolffd@0:                 for z = 1:length(d{k})
wolffd@0:                      d{k}{z} = simf(d{k}{z});
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             m.similarity = postoption.simf;
wolffd@0:             m = set(m,'Title','Similarity matrix','Data',d);
wolffd@0:         elseif length(m.similarity) == 1 && isnan(m.similarity) ...
wolffd@0:                 && option.dissim
wolffd@0:             m.similarity = 0;
wolffd@0:         end
wolffd@0:     end
wolffd@0:     if postoption.filt
wolffd@0:         fp = get(m,'FramePos');
wolffd@0:         for k = 1:length(d)
wolffd@0:             for z = 1:length(d{k})
wolffd@0:                 dz = filter(ones(postoption.filt,1),1,d{k}{z});
wolffd@0:                 d{k}{z} = dz(postoption.filt:end,1:end-postoption.filt+1);
wolffd@0:                 fp{k}{z} = [fp{k}{z}(1,1:end-postoption.filt+1);...
wolffd@0:                             fp{k}{z}(1,postoption.filt:end)];
wolffd@0:             end
wolffd@0:         end
wolffd@0:         m = set(m,'Data',d,'FramePos',fp);
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function S = rotatesim(d,K)
wolffd@0: if length(d) == 1;
wolffd@0:     S = d;
wolffd@0: else
wolffd@0:     K = min(K,size(d,1)*2);
wolffd@0:     lK = floor(K/2);
wolffd@0:     S = NaN(K,size(d,2),size(d,3));
wolffd@0:     for k = 1:size(d,3)
wolffd@0:         for j = -lK:lK
wolffd@0:             S(lK+1+j,:,k) = [NaN(1,floor(abs(j)/2)) diag(d(:,:,k),j)' ...
wolffd@0:                                                     NaN(1,ceil(abs(j)/2))];
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: 
wolffd@0: function d = cosine(r,s)
wolffd@0: d = 1-r'*s;
wolffd@0: %nr = sqrt(r'*r);
wolffd@0: %ns = sqrt(s'*s);
wolffd@0: %if or(nr == 0, ns == 0);
wolffd@0: %    d = 1;
wolffd@0: %else
wolffd@0: %    d = 1 - r'*s/nr/ns;
wolffd@0: %end
wolffd@0: 
wolffd@0: 
wolffd@0: function d = KL(x,y)
wolffd@0: % Kullback-Leibler distance
wolffd@0: if size(x,4)>1
wolffd@0:     x(end+1:2*end,:,:,1) = x(:,:,:,2);
wolffd@0:     x(:,:,:,2) = [];
wolffd@0: end
wolffd@0: if size(y,4)>1
wolffd@0:     y(end+1:2*end,:,:,1) = y(:,:,:,2);
wolffd@0:     y(:,:,:,2) = [];
wolffd@0: end
wolffd@0: m1 = mean(x);
wolffd@0: m2 = mean(y);
wolffd@0: S1 = cov(x);
wolffd@0: S2 = cov(y);
wolffd@0: d = (trace(S1/S2)+trace(S2/S1)+(m1-m2)'*inv(S1+S2)*(m1-m2))/2 - size(S1,1);
wolffd@0:     
wolffd@0: 
wolffd@0: function s = exponential(d)
wolffd@0:     s = exp(-d);
wolffd@0:     
wolffd@0:     
wolffd@0: function s = oneminus(d)
wolffd@0:     s = 1-d;