Daniel@0: function varargout = mirflux(orig,varargin)
Daniel@0: %   f = mirflux(x) measures distance between successive frames.
Daniel@0: %   First argument:
Daniel@0: %       If x is a spectrum, this corresponds to spectral flux.
Daniel@0: %       But the flux of any other data can be computed as well.
Daniel@0: %       If x is an audio file or audio signal, the spectral flux is
Daniel@0: %           computed by default.
Daniel@0: %   Optional arguments:
Daniel@0: %       f = mirflux(x,'Frame',...) specifies the frame parameters, if x is
Daniel@0: %           not already decomposed into frames. Default values are frame 
Daniel@0: %           length of .2 seconds and hop factor of 1.3.
Daniel@0: %       f = mirflux(x,'Dist',d) specifies the distance between 
Daniel@0: %           successive  frames:                 (IF 'COMPLEX': DISTANCE = 'CITY' ALWAYS)
Daniel@0: %           d = 'Euclidian': Euclidian distance (Default)
Daniel@0: %           d = 'City': City-block distance
Daniel@0: %           d = 'Cosine': Cosine distance (or normalized correlation)
Daniel@0: %       f = mirflux(...,'Inc'): Only positive difference between frames are
Daniel@0: %           summed, in order to focus on increase of energy solely.
Daniel@0: %       f = mirflux(...,'Complex'), for spectral flux, combines use of
Daniel@0: %           energy and phase information (Bello et al, 2004).
Daniel@0: %       f = mirflux(...,'Halfwave'): performs a half-wave  rectification on
Daniel@0: %           the result.
Daniel@0: %       f = mirflux(...,'Median',l,C): removes small spurious peaks by
Daniel@0: %           subtracting to the result its median filtering. The median
Daniel@0: %           filter computes the point-wise median inside a window of length
Daniel@0: %           l (in seconds), that  includes a same number of previous and
Daniel@0: %           next samples. C is a scaling factor whose purpose is to
Daniel@0: %           artificially rise the curve slightly above the steady state of
Daniel@0: %           the signal. If no parameters are given, the default values are:
Daniel@0: %           l = 0.2 s. and C = 1.3
Daniel@0: %       f = mirflux(...,'Median',l,C,'Halfwave'): The scaled median
Daniel@0: %           filtering is designed to be succeeded by the  half-wave 
Daniel@0: %           rectification process in order to select peaks above the
Daniel@0: %           dynamic threshold calculated with the help of the median
Daniel@0: %           filter. The resulting signal is called "detection function"
Daniel@0: %           (Alonso, David, Richard, 2004). To ensure accurate detection,
Daniel@0: %           the length l of the median filter must be longer than the
Daniel@0: %           average width of the peaks of the detection function.
Daniel@0: %
Daniel@0: %   (Bello et al, 2004) Juan P. Bello, Chris Duxbury, Mike Davies, and Mark
Daniel@0: %   Sandler, On the Use of Phase and Energy for Musical Onset Detection in
Daniel@0: %   the Complex Domain, IEEE SIGNAL PROCESSING LETTERS, VOL. 11, NO. 6,
Daniel@0: %   JUNE 2004
Daniel@0: 
Daniel@0:         dist.key = 'Dist';
Daniel@0:         dist.type = 'String';
Daniel@0:         dist.choice = {'Euclidian','City','Cosine'};  %PDIST???? (euclidean)
Daniel@0:         dist.default = 'Euclidian';
Daniel@0:     option.dist = dist;
Daniel@0:     
Daniel@0:         inc.key = 'Inc';
Daniel@0:         inc.type = 'Boolean';
Daniel@0:         inc.default = 0;
Daniel@0:     option.inc = inc;
Daniel@0:     
Daniel@0:         complex.key = 'Complex';
Daniel@0:         complex.type = 'Boolean';
Daniel@0:         complex.default = 0;
Daniel@0:     option.complex = complex;
Daniel@0:     
Daniel@0:         h.key = 'Halfwave';
Daniel@0:         h.type = 'Boolean';
Daniel@0:         h.default = 0;
Daniel@0:         h.when = 'After';
Daniel@0:     option.h = h;
Daniel@0:     
Daniel@0:         median.key = 'Median';
Daniel@0:         median.type = 'Integer';
Daniel@0:         median.number = 2;
Daniel@0:         median.default = [0 0];
Daniel@0:         median.keydefault = [.2 1.3];
Daniel@0:         median.when = 'After';
Daniel@0:     option.median = median;
Daniel@0:     
Daniel@0:         frame.key = 'Frame';
Daniel@0:         frame.type = 'Integer';
Daniel@0:         frame.number = 2;
Daniel@0:         frame.default = [.05 .5];
Daniel@0:     option.frame = frame;
Daniel@0:     
Daniel@0: specif.option = option;
Daniel@0:      
Daniel@0: varargout = mirfunction(@mirflux,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: if isamir(x,'miraudio') 
Daniel@0:     if isframed(x)
Daniel@0:         x = mirspectrum(x);
Daniel@0:     else
Daniel@0:         x = mirspectrum(x,'Frame',option.frame.length.val,option.frame.length.unit,...
Daniel@0:                                   option.frame.hop.val,option.frame.hop.unit);
Daniel@0:     end
Daniel@0: end
Daniel@0: if isa(x,'mirdesign')
Daniel@0:     x = set(x,'Overlap',1);
Daniel@0: end
Daniel@0: type = 'mirscalar';
Daniel@0: 
Daniel@0: 
Daniel@0: function f = main(s,option,postoption)
Daniel@0: if iscell(s)
Daniel@0:     s = s{1};
Daniel@0: end
Daniel@0: t = get(s,'Title');
Daniel@0: if isa(s,'mirscalar') && ...
Daniel@0:         (strcmp(t,'Harmonic Change Detection Function') || ...
Daniel@0:          (length(t)>4 && strcmp(t(end-3:end),'flux')) || ...
Daniel@0:          (length(t)>5 && strcmp(t(end-4:end-1),'flux')))
Daniel@0:     if not(isempty(postoption))
Daniel@0:         f = modif(s,postoption);
Daniel@0:     end
Daniel@0: else
Daniel@0:     if isa(s,'mirspectrum')
Daniel@0:         t = 'Spectral';
Daniel@0:     end
Daniel@0:     m = get(s,'Data'); 
Daniel@0:     if option.complex
Daniel@0:         if not(isa(s,'mirspectrum'))
Daniel@0:             error('ERROR IN MIRFLUX: ''Complex'' option only defined for spectral flux.');
Daniel@0:         end
Daniel@0:         ph = get(s,'Phase');
Daniel@0:     end
Daniel@0:     param.complex = option.complex;
Daniel@0:     param.inc = option.inc;
Daniel@0:     fp = get(s,'FramePos');
Daniel@0:     if strcmp(t,'Tonal centroid')
Daniel@0:         t = 'Harmonic Change Detection Function';
Daniel@0:     else
Daniel@0:         t = [t,' flux'];
Daniel@0:     end
Daniel@0:     disp(['Computing ' t '...'])
Daniel@0:     ff = cell(1,length(m));
Daniel@0:     newsr = cell(1,length(m));
Daniel@0:     dist = str2func(option.dist);
Daniel@0:     %[tmp s] = gettmp(s); %get(s,'Tmp');
Daniel@0:     for h = 1:length(m)
Daniel@0:         ff{h} = cell(1,length(m{h}));
Daniel@0:         if not(iscell(m{h}))
Daniel@0:             m{h} = {m{h}};
Daniel@0:         end
Daniel@0:         for i = 1:length(m{h})
Daniel@0:             mi = m{h}{i};
Daniel@0:             if size(mi,3) > 1 && size(mi,1) == 1
Daniel@0:                 mi = reshape(mi,size(mi,2),size(mi,3))';
Daniel@0:             end
Daniel@0:             if option.complex
Daniel@0:                 phi = ph{h}{i};
Daniel@0:             end
Daniel@0:             fpi = fp{h}{i};
Daniel@0:             %if 0 %not(isempty(tmp)) && isstruct(tmp) && isfield(tmp,'mi')
Daniel@0:             %    mi = [tmp.mi mi];
Daniel@0:             %    fpi = [tmp.fpi fpi];
Daniel@0:             %end
Daniel@0:             nc = size(mi,2);
Daniel@0:             np = size(mi,3);
Daniel@0:             if nc == 1
Daniel@0:                 warning('WARNING IN MIRFLUX: Flux can only be computed on signal decomposed into frames.');
Daniel@0:                 ff{h}{i} = NaN(1,1,np);
Daniel@0:             else
Daniel@0:                 if option.complex
Daniel@0:                     fl = zeros(1,nc-2,np);
Daniel@0:                     for k = 1:np
Daniel@0:                         d = diff(phi(:,:,k),2,2);
Daniel@0:                         d = d/(2*pi) - round(d/(2*pi));
Daniel@0:                         g = sqrt(mi(:,3:end,k).^2 + mi(:,2:(end-1),k).^2 ...
Daniel@0:                                                   - 2.*mi(:,3:end,k)...
Daniel@0:                                                      .*mi(:,2:(end-1),k)...
Daniel@0:                                                      .*cos(d));
Daniel@0:                         fl(1,:,k) = sum(g);
Daniel@0:                     end
Daniel@0:                     fp{h}{i} = fpi(:,3:end);
Daniel@0:                 else
Daniel@0:                     fl = zeros(1,nc-1,np);
Daniel@0:                     for k = 1:np
Daniel@0:                         for j = 1:nc-1
Daniel@0:                             fl(1,j,k) = dist(mi(:,j,k),mi(:,j+1,k),option.inc);
Daniel@0:                         end
Daniel@0:                     end
Daniel@0:                     fp{h}{i} = fpi(:,2:end);
Daniel@0:                 end
Daniel@0:                 ff{h}{i} = fl;
Daniel@0:                 %tmp.mi = mi(:,end,:);
Daniel@0:                 %tmp.fpi = fpi(:,end,:);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         %tmp = [];
Daniel@0:         if size(fpi,2)<2
Daniel@0:             newsr{h} = 0;
Daniel@0:         else
Daniel@0:             newsr{h} = 1/(fpi(1,2)-fpi(1,1));
Daniel@0:         end
Daniel@0:     end
Daniel@0:     f = mirscalar(s,'Data',ff,'FramePos',fp,'Sampling',newsr,...
Daniel@0:                         'Title',t,'Parameter',param); %,'Tmp',tmp);
Daniel@0:     %f = settmp(f,tmp);
Daniel@0:     if not(isempty(postoption))
Daniel@0:         f = modif(f,postoption);
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function f = modif(f,option)
Daniel@0: fl = get(f,'Data'); 
Daniel@0: r = get(f,'Sampling');
Daniel@0: for h = 1:length(fl)
Daniel@0:     for i = 1:length(fl{h})
Daniel@0:         fli = fl{h}{i};
Daniel@0:         nc = size(fli,2);
Daniel@0:         np = size(fli,3);
Daniel@0:         if option.median(1)
Daniel@0:             ffi = zeros(1,nc,np);
Daniel@0:             lr = round(option.median(1)*r{i});
Daniel@0:             for k = 1:np
Daniel@0:                 for j = 1:nc
Daniel@0:                     ffi(:,j,k) = fli(:,j,k) - ...
Daniel@0:                         option.median(2) * median(fli(:,max(1,j-lr):min(nc-1,j+lr),k));
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             fli = ffi;
Daniel@0:         end
Daniel@0:         if option.h
Daniel@0:             fli = hwr(fli);
Daniel@0:         end
Daniel@0:         fl{h}{i} = fli;
Daniel@0:     end
Daniel@0: end
Daniel@0: f = set(f,'Data',fl);
Daniel@0: 
Daniel@0: 
Daniel@0: function y = Euclidian(mi,mj,inc)
Daniel@0: if inc
Daniel@0:     y = sqrt(sum(max(0,(mj-mi)).^2));
Daniel@0: else
Daniel@0:     y = sqrt(sum((mj-mi).^2));
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function y = City(mi,mj,inc)
Daniel@0: if inc
Daniel@0:     y = sum(max(0,(mj-mi)));
Daniel@0: else
Daniel@0:     y = sum(abs(mj-mi));
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function d = Cosine(r,s,inc)
Daniel@0: nr = sqrt(r'*r);
Daniel@0: ns = sqrt(s'*s);
Daniel@0: if or(nr == 0, ns == 0);
Daniel@0:     d = 1;
Daniel@0: else
Daniel@0:     d = 1 - r'*s/nr/ns;
Daniel@0: end