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