Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/MIRToolbox/mirflux.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 = mirflux(orig,varargin)
wolffd@0 2 % f = mirflux(x) measures distance between successive frames.
wolffd@0 3 % First argument:
wolffd@0 4 % If x is a spectrum, this corresponds to spectral flux.
wolffd@0 5 % But the flux of any other data can be computed as well.
wolffd@0 6 % If x is an audio file or audio signal, the spectral flux is
wolffd@0 7 % computed by default.
wolffd@0 8 % Optional arguments:
wolffd@0 9 % f = mirflux(x,'Frame',...) specifies the frame parameters, if x is
wolffd@0 10 % not already decomposed into frames. Default values are frame
wolffd@0 11 % length of .2 seconds and hop factor of 1.3.
wolffd@0 12 % f = mirflux(x,'Dist',d) specifies the distance between
wolffd@0 13 % successive frames: (IF 'COMPLEX': DISTANCE = 'CITY' ALWAYS)
wolffd@0 14 % d = 'Euclidian': Euclidian distance (Default)
wolffd@0 15 % d = 'City': City-block distance
wolffd@0 16 % d = 'Cosine': Cosine distance (or normalized correlation)
wolffd@0 17 % f = mirflux(...,'Inc'): Only positive difference between frames are
wolffd@0 18 % summed, in order to focus on increase of energy solely.
wolffd@0 19 % f = mirflux(...,'Complex'), for spectral flux, combines use of
wolffd@0 20 % energy and phase information (Bello et al, 2004).
wolffd@0 21 % f = mirflux(...,'Halfwave'): performs a half-wave rectification on
wolffd@0 22 % the result.
wolffd@0 23 % f = mirflux(...,'Median',l,C): removes small spurious peaks by
wolffd@0 24 % subtracting to the result its median filtering. The median
wolffd@0 25 % filter computes the point-wise median inside a window of length
wolffd@0 26 % l (in seconds), that includes a same number of previous and
wolffd@0 27 % next samples. C is a scaling factor whose purpose is to
wolffd@0 28 % artificially rise the curve slightly above the steady state of
wolffd@0 29 % the signal. If no parameters are given, the default values are:
wolffd@0 30 % l = 0.2 s. and C = 1.3
wolffd@0 31 % f = mirflux(...,'Median',l,C,'Halfwave'): The scaled median
wolffd@0 32 % filtering is designed to be succeeded by the half-wave
wolffd@0 33 % rectification process in order to select peaks above the
wolffd@0 34 % dynamic threshold calculated with the help of the median
wolffd@0 35 % filter. The resulting signal is called "detection function"
wolffd@0 36 % (Alonso, David, Richard, 2004). To ensure accurate detection,
wolffd@0 37 % the length l of the median filter must be longer than the
wolffd@0 38 % average width of the peaks of the detection function.
wolffd@0 39 %
wolffd@0 40 % (Bello et al, 2004) Juan P. Bello, Chris Duxbury, Mike Davies, and Mark
wolffd@0 41 % Sandler, On the Use of Phase and Energy for Musical Onset Detection in
wolffd@0 42 % the Complex Domain, IEEE SIGNAL PROCESSING LETTERS, VOL. 11, NO. 6,
wolffd@0 43 % JUNE 2004
wolffd@0 44
wolffd@0 45 dist.key = 'Dist';
wolffd@0 46 dist.type = 'String';
wolffd@0 47 dist.choice = {'Euclidian','City','Cosine'}; %PDIST???? (euclidean)
wolffd@0 48 dist.default = 'Euclidian';
wolffd@0 49 option.dist = dist;
wolffd@0 50
wolffd@0 51 inc.key = 'Inc';
wolffd@0 52 inc.type = 'Boolean';
wolffd@0 53 inc.default = 0;
wolffd@0 54 option.inc = inc;
wolffd@0 55
wolffd@0 56 complex.key = 'Complex';
wolffd@0 57 complex.type = 'Boolean';
wolffd@0 58 complex.default = 0;
wolffd@0 59 option.complex = complex;
wolffd@0 60
wolffd@0 61 h.key = 'Halfwave';
wolffd@0 62 h.type = 'Boolean';
wolffd@0 63 h.default = 0;
wolffd@0 64 h.when = 'After';
wolffd@0 65 option.h = h;
wolffd@0 66
wolffd@0 67 median.key = 'Median';
wolffd@0 68 median.type = 'Integer';
wolffd@0 69 median.number = 2;
wolffd@0 70 median.default = [0 0];
wolffd@0 71 median.keydefault = [.2 1.3];
wolffd@0 72 median.when = 'After';
wolffd@0 73 option.median = median;
wolffd@0 74
wolffd@0 75 frame.key = 'Frame';
wolffd@0 76 frame.type = 'Integer';
wolffd@0 77 frame.number = 2;
wolffd@0 78 frame.default = [.05 .5];
wolffd@0 79 option.frame = frame;
wolffd@0 80
wolffd@0 81 specif.option = option;
wolffd@0 82
wolffd@0 83 varargout = mirfunction(@mirflux,orig,varargin,nargout,specif,@init,@main);
wolffd@0 84
wolffd@0 85
wolffd@0 86 function [x type] = init(x,option)
wolffd@0 87 if isamir(x,'miraudio')
wolffd@0 88 if isframed(x)
wolffd@0 89 x = mirspectrum(x);
wolffd@0 90 else
wolffd@0 91 x = mirspectrum(x,'Frame',option.frame.length.val,option.frame.length.unit,...
wolffd@0 92 option.frame.hop.val,option.frame.hop.unit);
wolffd@0 93 end
wolffd@0 94 end
wolffd@0 95 if isa(x,'mirdesign')
wolffd@0 96 x = set(x,'Overlap',1);
wolffd@0 97 end
wolffd@0 98 type = 'mirscalar';
wolffd@0 99
wolffd@0 100
wolffd@0 101 function f = main(s,option,postoption)
wolffd@0 102 if iscell(s)
wolffd@0 103 s = s{1};
wolffd@0 104 end
wolffd@0 105 t = get(s,'Title');
wolffd@0 106 if isa(s,'mirscalar') && ...
wolffd@0 107 (strcmp(t,'Harmonic Change Detection Function') || ...
wolffd@0 108 (length(t)>4 && strcmp(t(end-3:end),'flux')) || ...
wolffd@0 109 (length(t)>5 && strcmp(t(end-4:end-1),'flux')))
wolffd@0 110 if not(isempty(postoption))
wolffd@0 111 f = modif(s,postoption);
wolffd@0 112 end
wolffd@0 113 else
wolffd@0 114 if isa(s,'mirspectrum')
wolffd@0 115 t = 'Spectral';
wolffd@0 116 end
wolffd@0 117 m = get(s,'Data');
wolffd@0 118 if option.complex
wolffd@0 119 if not(isa(s,'mirspectrum'))
wolffd@0 120 error('ERROR IN MIRFLUX: ''Complex'' option only defined for spectral flux.');
wolffd@0 121 end
wolffd@0 122 ph = get(s,'Phase');
wolffd@0 123 end
wolffd@0 124 param.complex = option.complex;
wolffd@0 125 param.inc = option.inc;
wolffd@0 126 fp = get(s,'FramePos');
wolffd@0 127 if strcmp(t,'Tonal centroid')
wolffd@0 128 t = 'Harmonic Change Detection Function';
wolffd@0 129 else
wolffd@0 130 t = [t,' flux'];
wolffd@0 131 end
wolffd@0 132 disp(['Computing ' t '...'])
wolffd@0 133 ff = cell(1,length(m));
wolffd@0 134 newsr = cell(1,length(m));
wolffd@0 135 dist = str2func(option.dist);
wolffd@0 136 %[tmp s] = gettmp(s); %get(s,'Tmp');
wolffd@0 137 for h = 1:length(m)
wolffd@0 138 ff{h} = cell(1,length(m{h}));
wolffd@0 139 if not(iscell(m{h}))
wolffd@0 140 m{h} = {m{h}};
wolffd@0 141 end
wolffd@0 142 for i = 1:length(m{h})
wolffd@0 143 mi = m{h}{i};
wolffd@0 144 if size(mi,3) > 1 && size(mi,1) == 1
wolffd@0 145 mi = reshape(mi,size(mi,2),size(mi,3))';
wolffd@0 146 end
wolffd@0 147 if option.complex
wolffd@0 148 phi = ph{h}{i};
wolffd@0 149 end
wolffd@0 150 fpi = fp{h}{i};
wolffd@0 151 %if 0 %not(isempty(tmp)) && isstruct(tmp) && isfield(tmp,'mi')
wolffd@0 152 % mi = [tmp.mi mi];
wolffd@0 153 % fpi = [tmp.fpi fpi];
wolffd@0 154 %end
wolffd@0 155 nc = size(mi,2);
wolffd@0 156 np = size(mi,3);
wolffd@0 157 if nc == 1
wolffd@0 158 warning('WARNING IN MIRFLUX: Flux can only be computed on signal decomposed into frames.');
wolffd@0 159 ff{h}{i} = NaN(1,1,np);
wolffd@0 160 else
wolffd@0 161 if option.complex
wolffd@0 162 fl = zeros(1,nc-2,np);
wolffd@0 163 for k = 1:np
wolffd@0 164 d = diff(phi(:,:,k),2,2);
wolffd@0 165 d = d/(2*pi) - round(d/(2*pi));
wolffd@0 166 g = sqrt(mi(:,3:end,k).^2 + mi(:,2:(end-1),k).^2 ...
wolffd@0 167 - 2.*mi(:,3:end,k)...
wolffd@0 168 .*mi(:,2:(end-1),k)...
wolffd@0 169 .*cos(d));
wolffd@0 170 fl(1,:,k) = sum(g);
wolffd@0 171 end
wolffd@0 172 fp{h}{i} = fpi(:,3:end);
wolffd@0 173 else
wolffd@0 174 fl = zeros(1,nc-1,np);
wolffd@0 175 for k = 1:np
wolffd@0 176 for j = 1:nc-1
wolffd@0 177 fl(1,j,k) = dist(mi(:,j,k),mi(:,j+1,k),option.inc);
wolffd@0 178 end
wolffd@0 179 end
wolffd@0 180 fp{h}{i} = fpi(:,2:end);
wolffd@0 181 end
wolffd@0 182 ff{h}{i} = fl;
wolffd@0 183 %tmp.mi = mi(:,end,:);
wolffd@0 184 %tmp.fpi = fpi(:,end,:);
wolffd@0 185 end
wolffd@0 186 end
wolffd@0 187 %tmp = [];
wolffd@0 188 if size(fpi,2)<2
wolffd@0 189 newsr{h} = 0;
wolffd@0 190 else
wolffd@0 191 newsr{h} = 1/(fpi(1,2)-fpi(1,1));
wolffd@0 192 end
wolffd@0 193 end
wolffd@0 194 f = mirscalar(s,'Data',ff,'FramePos',fp,'Sampling',newsr,...
wolffd@0 195 'Title',t,'Parameter',param); %,'Tmp',tmp);
wolffd@0 196 %f = settmp(f,tmp);
wolffd@0 197 if not(isempty(postoption))
wolffd@0 198 f = modif(f,postoption);
wolffd@0 199 end
wolffd@0 200 end
wolffd@0 201
wolffd@0 202
wolffd@0 203 function f = modif(f,option)
wolffd@0 204 fl = get(f,'Data');
wolffd@0 205 r = get(f,'Sampling');
wolffd@0 206 for h = 1:length(fl)
wolffd@0 207 for i = 1:length(fl{h})
wolffd@0 208 fli = fl{h}{i};
wolffd@0 209 nc = size(fli,2);
wolffd@0 210 np = size(fli,3);
wolffd@0 211 if option.median(1)
wolffd@0 212 ffi = zeros(1,nc,np);
wolffd@0 213 lr = round(option.median(1)*r{i});
wolffd@0 214 for k = 1:np
wolffd@0 215 for j = 1:nc
wolffd@0 216 ffi(:,j,k) = fli(:,j,k) - ...
wolffd@0 217 option.median(2) * median(fli(:,max(1,j-lr):min(nc-1,j+lr),k));
wolffd@0 218 end
wolffd@0 219 end
wolffd@0 220 fli = ffi;
wolffd@0 221 end
wolffd@0 222 if option.h
wolffd@0 223 fli = hwr(fli);
wolffd@0 224 end
wolffd@0 225 fl{h}{i} = fli;
wolffd@0 226 end
wolffd@0 227 end
wolffd@0 228 f = set(f,'Data',fl);
wolffd@0 229
wolffd@0 230
wolffd@0 231 function y = Euclidian(mi,mj,inc)
wolffd@0 232 if inc
wolffd@0 233 y = sqrt(sum(max(0,(mj-mi)).^2));
wolffd@0 234 else
wolffd@0 235 y = sqrt(sum((mj-mi).^2));
wolffd@0 236 end
wolffd@0 237
wolffd@0 238
wolffd@0 239 function y = City(mi,mj,inc)
wolffd@0 240 if inc
wolffd@0 241 y = sum(max(0,(mj-mi)));
wolffd@0 242 else
wolffd@0 243 y = sum(abs(mj-mi));
wolffd@0 244 end
wolffd@0 245
wolffd@0 246
wolffd@0 247 function d = Cosine(r,s,inc)
wolffd@0 248 nr = sqrt(r'*r);
wolffd@0 249 ns = sqrt(s'*s);
wolffd@0 250 if or(nr == 0, ns == 0);
wolffd@0 251 d = 1;
wolffd@0 252 else
wolffd@0 253 d = 1 - r'*s/nr/ns;
wolffd@0 254 end