Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirsimatrix/mirsimatrix.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function varargout = mirsimatrix(orig,varargin)
wolffd@0 2 % m = mirsimatrix(x) computes the similarity matrix resulting from the
wolffd@0 3 % mutual comparison between each possible frame analysis in x.
wolffd@0 4 % By default, x is the spectrum of the frame decomposition.
wolffd@0 5 % But it can be any other frame analysis.
wolffd@0 6 % Optional argument:
wolffd@0 7 % mirsimatrix(...,'Distance',f) specifies the name of a distance
wolffd@0 8 % function, from those proposed in the Statistics Toolbox
wolffd@0 9 % (help pdist).
wolffd@0 10 % default value: f = 'cosine'
wolffd@0 11 % mirsimatrix(...,'Dissimilarity') return the dissimilarity matrix,
wolffd@0 12 % which is the intermediary result before the computation of the
wolffd@0 13 % actual similarity matrix. It shows the distance between each
wolffd@0 14 % possible frame analysis in x.
wolffd@0 15 % mirsimatrix(...,'Similarity',f) indicates the function f specifying
wolffd@0 16 % the way the distance values in the dissimilarity matrix are
wolffd@0 17 % transformed into similarity values.
wolffd@0 18 % Possible values:
wolffd@0 19 % f = 'oneminus' (default value)
wolffd@0 20 % corresponding to f(x) = 1-x
wolffd@0 21 % f = 'exponential'
wolffd@0 22 % corresponding to f(x)= exp(-x)
wolffd@0 23 % mirsimatrix(...,'Width',w) or more simply dissimatrix(...,w)
wolffd@0 24 % specifies the size of the diagonal bandwidth, in samples,
wolffd@0 25 % outside which the dissimilarity will not be computed.
wolffd@0 26 % if w = inf (default value), all the matrix will be computed.
wolffd@0 27 % mirsimatrix(...,'Horizontal') rotates the matrix 45 degrees in order to
wolffd@0 28 % make the first diagonal horizontal, and to restrict on the
wolffd@0 29 % diagonal bandwidth only.
wolffd@0 30 % mirsimatrix(...,'TimeLag') transforms the (non-rotated) matrix into
wolffd@0 31 % a time-lag matrix, making the first diagonal horizontal as well
wolffd@0 32 % (corresponding to the zero-lag line).
wolffd@0 33 %
wolffd@0 34 % mirsimatrix(M,r) creates a mirsimatrix similarity matrix based on
wolffd@0 35 % the Matlab square matrix M, of frame rate r (in Hz.)
wolffd@0 36 % By default r = 20 Hz.
wolffd@0 37 % mirsimatrix(M,r,'Dissimilarity') creates instead a mirsimatrix
wolffd@0 38 % dissimilarity matrix.
wolffd@0 39 %
wolffd@0 40 % Foote, J. & Cooper, M. (2003). Media Segmentation using Self-Similarity
wolffd@0 41 % Decomposition,. In Proc. SPIE Storage and Retrieval for Multimedia
wolffd@0 42 % Databases, Vol. 5021, pp. 167-75.
wolffd@0 43
wolffd@0 44 % mirsimatrix(...,'Filter',10) filter along the diagonal of the matrix,
wolffd@0 45 % using a uniform moving average filter of size f. The result is
wolffd@0 46 % represented in a time-lag matrix.
wolffd@0 47
wolffd@0 48
wolffd@0 49 distance.key = 'Distance';
wolffd@0 50 distance.type = 'String';
wolffd@0 51 distance.default = 'cosine';
wolffd@0 52 option.distance = distance;
wolffd@0 53
wolffd@0 54 simf.key = 'Similarity';
wolffd@0 55 simf.type = 'String';
wolffd@0 56 simf.default = 'oneminus';
wolffd@0 57 simf.keydefault = 'Similarity';
wolffd@0 58 simf.when = 'After';
wolffd@0 59 option.simf = simf;
wolffd@0 60
wolffd@0 61 dissim.key = 'Dissimilarity';
wolffd@0 62 dissim.type = 'Boolean';
wolffd@0 63 dissim.default = 0;
wolffd@0 64 option.dissim = dissim;
wolffd@0 65
wolffd@0 66 K.key = 'Width';
wolffd@0 67 K.type = 'Integer';
wolffd@0 68 K.default = Inf;
wolffd@0 69 option.K = K;
wolffd@0 70
wolffd@0 71 filt.key = 'Filter';
wolffd@0 72 filt.type = 'Integer';
wolffd@0 73 filt.default = 0;
wolffd@0 74 filt.when = 'After';
wolffd@0 75 option.filt = filt;
wolffd@0 76
wolffd@0 77 view.type = 'String';
wolffd@0 78 view.default = 'Standard';
wolffd@0 79 view.choice = {'Standard','Horizontal','TimeLag'};
wolffd@0 80 view.when = 'After';
wolffd@0 81 option.view = view;
wolffd@0 82
wolffd@0 83 rate.type = 'Integer';
wolffd@0 84 rate.position = 2;
wolffd@0 85 rate.default = 20;
wolffd@0 86 option.rate = rate;
wolffd@0 87
wolffd@0 88 specif.option = option;
wolffd@0 89 specif.nochunk = 1;
wolffd@0 90 varargout = mirfunction(@mirsimatrix,orig,varargin,nargout,specif,@init,@main);
wolffd@0 91
wolffd@0 92
wolffd@0 93 function [x type] = init(x,option)
wolffd@0 94 if isnumeric(x)
wolffd@0 95 m.diagwidth = Inf;
wolffd@0 96 m.view = 's';
wolffd@0 97 m.similarity = NaN;
wolffd@0 98 m.graph = {};
wolffd@0 99 m.branch = {};
wolffd@0 100 m = class(m,'mirsimatrix',mirdata);
wolffd@0 101 m = set(m,'Title','Dissimilarity matrix');
wolffd@0 102 fp = repmat(((1:size(x,1))-.5)/option.rate,[2,1]);
wolffd@0 103 x = set(m,'Data',{x},'Pos',[],...
wolffd@0 104 'FramePos',{{fp}},'Name',{inputname(1)});
wolffd@0 105 else
wolffd@0 106 if not(isamir(x,'mirsimatrix'))
wolffd@0 107 if (isamir(x,'miraudio'))
wolffd@0 108 if isframed(x)
wolffd@0 109 x = mirspectrum(x);
wolffd@0 110 else
wolffd@0 111 x = mirspectrum(x,'Frame',0.05,1);
wolffd@0 112 end
wolffd@0 113 end
wolffd@0 114 end
wolffd@0 115 end
wolffd@0 116 type = 'mirsimatrix';
wolffd@0 117
wolffd@0 118
wolffd@0 119 function m = main(orig,option,postoption)
wolffd@0 120 if iscell(orig)
wolffd@0 121 orig = orig{1};
wolffd@0 122 end
wolffd@0 123 if isa(orig,'mirsimatrix')
wolffd@0 124 d = get(orig,'Data');
wolffd@0 125 for k = 1:length(d)
wolffd@0 126 if iscell(option) && isfield(option,'K') && option.K < orig.diagwidth
wolffd@0 127 nl = size(d{k},1);
wolffd@0 128 if strcmp(orig.view,'h')
wolffd@0 129 dl = (nl - option.K)/2;
wolffd@0 130 dk = d{k}(ceil(dl):nl-floor(dl),:);
wolffd@0 131 else
wolffd@0 132 [spA,spd] = spdiags(d{k},-ceil(option.K/2):ceil(option.K/2));
wolffd@0 133 dk = full(spdiags(spA,spd,nl,size(d{k},2)));
wolffd@0 134 end
wolffd@0 135 d{k} = dk;
wolffd@0 136 orig.diagwidth = option.K;
wolffd@0 137 end
wolffd@0 138 end
wolffd@0 139 m = set(orig,'Data',d);
wolffd@0 140 else
wolffd@0 141 v = get(orig,'Data');
wolffd@0 142 d = cell(1,length(v));
wolffd@0 143 for k = 1:length(v)
wolffd@0 144 vk = v{k};
wolffd@0 145 if mirwaitbar
wolffd@0 146 handle = waitbar(0,'Computing dissimilarity matrix...');
wolffd@0 147 else
wolffd@0 148 handle = 0;
wolffd@0 149 end
wolffd@0 150 if 0 %iscell(vk)
wolffd@0 151 try
wolffd@0 152 vk = cell2mat(vk);
wolffd@0 153 end
wolffd@0 154 end
wolffd@0 155 if 0 %iscell(vk) %&& length(vk) > 1 %%% ATTENTION KL!!<<<<<<<<<<<<
wolffd@0 156 l = length(vk);
wolffd@0 157 dk = NaN(l,l);
wolffd@0 158 hK = floor(option.K/2);
wolffd@0 159 for i = 1:l
wolffd@0 160 if handle
wolffd@0 161 waitbar(i/l,handle);
wolffd@0 162 end
wolffd@0 163 for j = max(1,i-hK):min(l,i+hK)
wolffd@0 164 dk(i,j) = KL(vk{i},vk{j});
wolffd@0 165 end
wolffd@0 166 end
wolffd@0 167 else
wolffd@0 168 if not(iscell(vk))
wolffd@0 169 vk = {vk};
wolffd@0 170 end
wolffd@0 171 for z = 1:length(vk)
wolffd@0 172 vz = vk{z};
wolffd@0 173 ll = size(vz,1);
wolffd@0 174 l = size(vz,2);
wolffd@0 175 nc = size(vz,3);
wolffd@0 176 if ll==1 && nc>1
wolffd@0 177 vz = squeeze(vz)';
wolffd@0 178 ll = nc;
wolffd@0 179 nc = 1;
wolffd@0 180 end
wolffd@0 181 nd = size(vz,4);
wolffd@0 182 if not(isempty(postoption)) && ...
wolffd@0 183 strcmpi(postoption.view,'TimeLag')
wolffd@0 184 lK = ceil(option.K/2);
wolffd@0 185 if isinf(lK)
wolffd@0 186 lK = l;
wolffd@0 187 end
wolffd@0 188 dk{z} = NaN(lK,l,nc);
wolffd@0 189 else
wolffd@0 190 dk{z} = NaN(l,l,nc);
wolffd@0 191 end
wolffd@0 192 for g = 1:nc
wolffd@0 193 if nd == 1
wolffd@0 194 vv = vz;
wolffd@0 195 else
wolffd@0 196 vv = zeros(ll*nd,l);
wolffd@0 197 for h = 1:nd
wolffd@0 198 if iscell(vz)
wolffd@0 199 for i = 1:ll
wolffd@0 200 for j = 1:l
wolffd@0 201 vj = vz{i,j,g,h};
wolffd@0 202 if isempty(vj)
wolffd@0 203 vv((h-1)*ll+i,j) = NaN;
wolffd@0 204 else
wolffd@0 205 vv((h-1)*ll+i,j) = vj;
wolffd@0 206 end
wolffd@0 207 end
wolffd@0 208 end
wolffd@0 209 else
wolffd@0 210 tic
wolffd@0 211 vv((h-1)*ll+1:h*ll,:) = vz(:,:,g,h);
wolffd@0 212 toc
wolffd@0 213 end
wolffd@0 214 end
wolffd@0 215 end
wolffd@0 216 if isinf(option.K) && not(strcmpi(postoption.view,'TimeLag'))
wolffd@0 217 try
wolffd@0 218 manually = 0;
wolffd@0 219 dk{z}(:,:,g) = squareform(pdist(vv',option.distance));
wolffd@0 220 if option.K < Inf
wolffd@0 221 [spA,spd] = spdiags(dk{z},...
wolffd@0 222 -ceil(option.K/2):ceil(option.K/2));
wolffd@0 223 dk{z}(:,:,g) = full(spdiags(spA,spd,size(dk,1),size(dk{z},2)));
wolffd@0 224 end
wolffd@0 225 catch
wolffd@0 226 %err = lasterror;
wolffd@0 227 %warning(err.message)
wolffd@0 228 %disp('Statistics Toolbox does not seem to be
wolffd@0 229 %installed. Recompute the distance matrix manually.');
wolffd@0 230 manually = 1;
wolffd@0 231 end
wolffd@0 232 else
wolffd@0 233 manually = 1;
wolffd@0 234 end
wolffd@0 235 if manually
wolffd@0 236 disf = str2func(option.distance);
wolffd@0 237 if strcmpi(option.distance,'cosine')
wolffd@0 238 for i = 1:l
wolffd@0 239 vv(:,i) = vv(:,i)/norm(vv(:,i));
wolffd@0 240 end
wolffd@0 241 end
wolffd@0 242 if not(isempty(postoption)) && ...
wolffd@0 243 strcmpi(postoption.view,'TimeLag')
wolffd@0 244 lK = ceil(option.K/2);
wolffd@0 245 for i = 1:l
wolffd@0 246 if mirwaitbar && (mod(i,100) == 1 || i == l)
wolffd@0 247 waitbar(i/l,handle);
wolffd@0 248 end
wolffd@0 249 ij = min(i+lK-1,l);
wolffd@0 250 dkij = disf(vv(:,i),vv(:,i:ij));
wolffd@0 251 for j = 1:ij-i+1
wolffd@0 252 dk{z}(j,i+j-1,g) = dkij(j);
wolffd@0 253 end
wolffd@0 254 end
wolffd@0 255 else
wolffd@0 256 for i = 1:l
wolffd@0 257 if mirwaitbar && (mod(i,100) == 1 || i == l)
wolffd@0 258 waitbar(i/l,handle);
wolffd@0 259 end
wolffd@0 260 j = min(i+option.K-1,l);
wolffd@0 261 dkij = disf(vv(:,i),vv(:,i:j));
wolffd@0 262 dk{z}(i,i:j,g) = dkij;
wolffd@0 263 dk{z}(i:j,i,g) = dkij';
wolffd@0 264 end
wolffd@0 265 end
wolffd@0 266 end
wolffd@0 267 end
wolffd@0 268 end
wolffd@0 269 end
wolffd@0 270 d{k} = dk;
wolffd@0 271 if handle
wolffd@0 272 delete(handle)
wolffd@0 273 end
wolffd@0 274 end
wolffd@0 275 m.diagwidth = option.K;
wolffd@0 276 if not(isempty(postoption)) && strcmpi(postoption.view,'TimeLag')
wolffd@0 277 m.view = 'l';
wolffd@0 278 else
wolffd@0 279 m.view = 's';
wolffd@0 280 end
wolffd@0 281 m.similarity = 0;
wolffd@0 282 m.graph = {};
wolffd@0 283 m.branch = {};
wolffd@0 284 m = class(m,'mirsimatrix',mirdata(orig));
wolffd@0 285 m = purgedata(m);
wolffd@0 286 m = set(m,'Title','Dissimilarity matrix');
wolffd@0 287 m = set(m,'Data',d,'Pos',[]);
wolffd@0 288 end
wolffd@0 289 if not(isempty(postoption))
wolffd@0 290 if strcmpi(m.view,'s')
wolffd@0 291 if strcmpi(postoption.view,'Horizontal')
wolffd@0 292 for k = 1:length(d)
wolffd@0 293 for z = 1:length(d{k})
wolffd@0 294 d{k}{z} = rotatesim(d{k}{z},m.diagwidth);
wolffd@0 295 if option.K < m.diagwidth
wolffd@0 296 W = size(d{k}{z},1);
wolffd@0 297 hW = ceil(W/2);
wolffd@0 298 hK = ceil(option.K/2);
wolffd@0 299 d{k}{z} = d{k}{z}(hW-hK:hW+hK,:);
wolffd@0 300 m.diagwidth = option.K;
wolffd@0 301 end
wolffd@0 302 end
wolffd@0 303 end
wolffd@0 304 m = set(m,'Data',d);
wolffd@0 305 m.view = 'h';
wolffd@0 306 elseif strcmpi(postoption.view,'TimeLag') || postoption.filt
wolffd@0 307 W = ceil(m.diagwidth/2);
wolffd@0 308 for k = 1:length(d)
wolffd@0 309 for z = 1:length(d{k})
wolffd@0 310 if isinf(W)
wolffd@0 311 dz = NaN(size(d{k}{z}));
wolffd@0 312 else
wolffd@0 313 dz = NaN(W,size(d{k}{z},2));
wolffd@0 314 end
wolffd@0 315 for l = 1:size(dz,1)
wolffd@0 316 dz(l,l:end) = diag(d{k}{z},l-1)';
wolffd@0 317 end
wolffd@0 318 if option.K < m.diagwidth
wolffd@0 319 dz = dz(1:ceil(option.K/2),:);
wolffd@0 320 m.diagwidth = option.K;
wolffd@0 321 end
wolffd@0 322 d{k}{z}= dz;
wolffd@0 323 end
wolffd@0 324 end
wolffd@0 325 m = set(m,'Data',d);
wolffd@0 326 m.view = 'l';
wolffd@0 327 end
wolffd@0 328 end
wolffd@0 329 if ischar(postoption.simf)
wolffd@0 330 if strcmpi(postoption.simf,'Similarity')
wolffd@0 331 if not(isequal(m.similarity,NaN))
wolffd@0 332 option.dissim = 0;
wolffd@0 333 end
wolffd@0 334 postoption.simf = 'oneminus';
wolffd@0 335 end
wolffd@0 336 if isequal(m.similarity,0) && isstruct(option) ...
wolffd@0 337 && isfield(option,'dissim') && not(option.dissim)
wolffd@0 338 simf = str2func(postoption.simf);
wolffd@0 339 for k = 1:length(d)
wolffd@0 340 for z = 1:length(d{k})
wolffd@0 341 d{k}{z} = simf(d{k}{z});
wolffd@0 342 end
wolffd@0 343 end
wolffd@0 344 m.similarity = postoption.simf;
wolffd@0 345 m = set(m,'Title','Similarity matrix','Data',d);
wolffd@0 346 elseif length(m.similarity) == 1 && isnan(m.similarity) ...
wolffd@0 347 && option.dissim
wolffd@0 348 m.similarity = 0;
wolffd@0 349 end
wolffd@0 350 end
wolffd@0 351 if postoption.filt
wolffd@0 352 fp = get(m,'FramePos');
wolffd@0 353 for k = 1:length(d)
wolffd@0 354 for z = 1:length(d{k})
wolffd@0 355 dz = filter(ones(postoption.filt,1),1,d{k}{z});
wolffd@0 356 d{k}{z} = dz(postoption.filt:end,1:end-postoption.filt+1);
wolffd@0 357 fp{k}{z} = [fp{k}{z}(1,1:end-postoption.filt+1);...
wolffd@0 358 fp{k}{z}(1,postoption.filt:end)];
wolffd@0 359 end
wolffd@0 360 end
wolffd@0 361 m = set(m,'Data',d,'FramePos',fp);
wolffd@0 362 end
wolffd@0 363 end
wolffd@0 364
wolffd@0 365
wolffd@0 366 function S = rotatesim(d,K)
wolffd@0 367 if length(d) == 1;
wolffd@0 368 S = d;
wolffd@0 369 else
wolffd@0 370 K = min(K,size(d,1)*2);
wolffd@0 371 lK = floor(K/2);
wolffd@0 372 S = NaN(K,size(d,2),size(d,3));
wolffd@0 373 for k = 1:size(d,3)
wolffd@0 374 for j = -lK:lK
wolffd@0 375 S(lK+1+j,:,k) = [NaN(1,floor(abs(j)/2)) diag(d(:,:,k),j)' ...
wolffd@0 376 NaN(1,ceil(abs(j)/2))];
wolffd@0 377 end
wolffd@0 378 end
wolffd@0 379 end
wolffd@0 380
wolffd@0 381 function d = cosine(r,s)
wolffd@0 382 d = 1-r'*s;
wolffd@0 383 %nr = sqrt(r'*r);
wolffd@0 384 %ns = sqrt(s'*s);
wolffd@0 385 %if or(nr == 0, ns == 0);
wolffd@0 386 % d = 1;
wolffd@0 387 %else
wolffd@0 388 % d = 1 - r'*s/nr/ns;
wolffd@0 389 %end
wolffd@0 390
wolffd@0 391
wolffd@0 392 function d = KL(x,y)
wolffd@0 393 % Kullback-Leibler distance
wolffd@0 394 if size(x,4)>1
wolffd@0 395 x(end+1:2*end,:,:,1) = x(:,:,:,2);
wolffd@0 396 x(:,:,:,2) = [];
wolffd@0 397 end
wolffd@0 398 if size(y,4)>1
wolffd@0 399 y(end+1:2*end,:,:,1) = y(:,:,:,2);
wolffd@0 400 y(:,:,:,2) = [];
wolffd@0 401 end
wolffd@0 402 m1 = mean(x);
wolffd@0 403 m2 = mean(y);
wolffd@0 404 S1 = cov(x);
wolffd@0 405 S2 = cov(y);
wolffd@0 406 d = (trace(S1/S2)+trace(S2/S1)+(m1-m2)'*inv(S1+S2)*(m1-m2))/2 - size(S1,1);
wolffd@0 407
wolffd@0 408
wolffd@0 409 function s = exponential(d)
wolffd@0 410 s = exp(-d);
wolffd@0 411
wolffd@0 412
wolffd@0 413 function s = oneminus(d)
wolffd@0 414 s = 1-d;