Daniel@0: function varargout = mirpulseclarity(orig,varargin)
Daniel@0: %   r = mirpulseclarity(x) estimates the rhythmic clarity, indicating the
Daniel@0: %       strength of the beats estimated by the mirtempo function.
Daniel@0: %   Optional arguments:
Daniel@0: %       mirpulseclarity(...,s): specifies a strategy for pulse clarity
Daniel@0: %           estimation.
Daniel@0: %           Possible values: 'MaxAutocor' (default), 'MinAutocor',
Daniel@0: %               'KurtosisAutocor', MeanPeaksAutocor', 'EntropyAutocor', 
Daniel@0: %               'InterfAutocor', 'TempoAutocor', 'ExtremEnvelop', 
Daniel@0: %               'Attack', 'Articulation'
Daniel@0: %       mirpulseclarity(...,'Frame',l,h): orders a frame decomposition of
Daniel@0: %           the audio input of window length l (in seconds) and hop factor
Daniel@0: %           h, expressed relatively to the window length.
Daniel@0: %           Default values: l = 5 seconds and h = .1
Daniel@0: %       Onset detection strategies: 'Envelope' (default), 'DiffEnvelope', 
Daniel@0: %           'SpectralFlux', 'Pitch'.
Daniel@0: %       Options related to the autocorrelation computation can be specified
Daniel@0: %           as well: 'Min', 'Max', 'Resonance', 'Enhanced'
Daniel@0: %       Options related to the tempo estimation can be specified here
Daniel@0: %           as well: 'Sum', 'Total', 'Contrast'.
Daniel@0: %       cf. User's Manual for more details.
Daniel@0: %   [r,a] = mirpulseclarity(x) also returns the beat autocorrelation.
Daniel@0:         
Daniel@0:         model.key = 'Model';
Daniel@0:         model.type = 'Integer';
Daniel@0:         model.default = 0;
Daniel@0:     option.model = model;
Daniel@0: 
Daniel@0:         stratg.type = 'String';
Daniel@0:         stratg.choice = {'MaxAutocor','MinAutocor','MeanPeaksAutocor',...
Daniel@0:                          'KurtosisAutocor','EntropyAutocor',...
Daniel@0:                          'InterfAutocor','TempoAutocor','ExtremEnvelop',...
Daniel@0:                          'Attack','Articulation'};    ...,'AttackDiff'
Daniel@0:         stratg.default = 'MaxAutocor';
Daniel@0:     option.stratg = stratg;
Daniel@0: 
Daniel@0:         frame.key = 'Frame';
Daniel@0:         frame.type = 'Integer';
Daniel@0:         frame.number = 2;
Daniel@0:         frame.keydefault = [5 .1];
Daniel@0:         frame.default = [0 0];
Daniel@0:     option.frame = frame;
Daniel@0:         
Daniel@0: %% options related to mironsets:  
Daniel@0: 
Daniel@0:         fea.type = 'String';
Daniel@0:         fea.choice = {'Envelope','DiffEnvelope','SpectralFlux','Pitch'};
Daniel@0:         fea.default = 'Envelope';
Daniel@0:     option.fea = fea;
Daniel@0:     
Daniel@0:     
Daniel@0:     %% options related to 'Envelope':
Daniel@0:     
Daniel@0:             envmeth.key = 'Method';
Daniel@0:             envmeth.type = 'String';
Daniel@0:             envmeth.choice = {'Filter','Spectro'};
Daniel@0:             envmeth.default = 'Spectro';
Daniel@0:         option.envmeth = envmeth;
Daniel@0: 
Daniel@0:             %% options related to 'Filter':
Daniel@0: 
Daniel@0:                 ftype.key = 'FilterType';
Daniel@0:                 ftype.type = 'String';
Daniel@0:                 ftype.choice = {'IIR','HalfHann'};
Daniel@0:                 ftype.default = 'IIR';
Daniel@0:             option.ftype = ftype;
Daniel@0: 
Daniel@0:                 fb.key = 'Filterbank';
Daniel@0:                 fb.type = 'Integer';
Daniel@0:                 fb.default = 20;
Daniel@0:             option.fb = fb;
Daniel@0: 
Daniel@0:                 fbtype.key = 'FilterbankType';
Daniel@0:                 fbtype.type = 'String';
Daniel@0:                 fbtype.choice = {'Gammatone','Scheirer','Klapuri'};
Daniel@0:                 fbtype.default = 'Scheirer';
Daniel@0:             option.fbtype = fbtype;
Daniel@0: 
Daniel@0:             %% options related to 'Spectro':
Daniel@0: 
Daniel@0:                 band.type = 'String';
Daniel@0:                 band.choice = {'Freq','Mel','Bark','Cents'};
Daniel@0:                 band.default = 'Freq';
Daniel@0:             option.band = band;
Daniel@0: 
Daniel@0: 
Daniel@0:             diffhwr.key = 'HalfwaveDiff';
Daniel@0:             diffhwr.type = 'Integer';
Daniel@0:             diffhwr.default = 0;
Daniel@0:             diffhwr.keydefault = 1;
Daniel@0:         option.diffhwr = diffhwr;
Daniel@0: 
Daniel@0:             lambda.key = 'Lambda';
Daniel@0:             lambda.type = 'Integer';
Daniel@0:             lambda.default = 1;
Daniel@0:         option.lambda = lambda;
Daniel@0: 
Daniel@0:             aver.key = 'Smooth';
Daniel@0:             aver.type = 'Integer';
Daniel@0:             aver.default = 0;
Daniel@0:             aver.keydefault = 30;
Daniel@0:         option.aver = aver;
Daniel@0: 
Daniel@0:             oplog.key = 'Log';
Daniel@0:             oplog.type = 'Boolean';
Daniel@0:             oplog.default = 0;
Daniel@0:         option.log = oplog;
Daniel@0:         
Daniel@0:             mu.key = 'Mu';
Daniel@0:             mu.type = 'Boolean';
Daniel@0:             mu.default = 1;
Daniel@0:         option.mu = mu;
Daniel@0: 
Daniel@0:     %% options related to 'SpectralFlux'
Daniel@0:     
Daniel@0:             inc.key = 'Inc';
Daniel@0:             inc.type = 'Boolean';
Daniel@0:             inc.default = 1;
Daniel@0:         option.inc = inc;
Daniel@0: 
Daniel@0:             median.key = 'Median';
Daniel@0:             median.type = 'Integer';
Daniel@0:             median.number = 2;
Daniel@0:             median.default = [0 0]; % Not same default as in mirtempo
Daniel@0:         option.median = median;
Daniel@0: 
Daniel@0:             hw.key = 'Halfwave';
Daniel@0:             hw.type = 'Boolean';
Daniel@0:             hw.default = 0; %NaN; %0; % Not same default as in mirtempo
Daniel@0:         option.hw = hw;    
Daniel@0:         
Daniel@0:     
Daniel@0: %% options related to mirattackslope
Daniel@0:         slope.type = 'String';
Daniel@0:         slope.choice = {'Diff','Gauss'};
Daniel@0:         slope.default = 'Diff';
Daniel@0:     option.slope = slope;
Daniel@0:     
Daniel@0: %% options related to mirautocor:    
Daniel@0:     
Daniel@0:         enh.key = 'Enhanced';
Daniel@0:         enh.type = 'Integers';
Daniel@0:         enh.default = [];
Daniel@0:         enh.keydefault = 2:10;
Daniel@0:     option.enh = enh;
Daniel@0:     
Daniel@0:         r.key = 'Resonance';
Daniel@0:         r.type = 'String';
Daniel@0:         r.choice = {'ToiviainenSnyder','vonNoorden',0,'off','no'};
Daniel@0:         r.default = 'ToiviainenSnyder';
Daniel@0:     option.r = r;
Daniel@0:         
Daniel@0:         mi.key = 'Min';
Daniel@0:         mi.type = 'Integer';
Daniel@0:         mi.default = 40;
Daniel@0:     option.mi = mi;
Daniel@0:         
Daniel@0:         ma.key = 'Max';
Daniel@0:         ma.type = 'Integer';
Daniel@0:         ma.default = 200;
Daniel@0:     option.ma = ma;    
Daniel@0: 
Daniel@0: %% options related to mirtempo:
Daniel@0: 
Daniel@0:         sum.key = 'Sum';
Daniel@0:         sum.type = 'String';
Daniel@0:         sum.choice = {'Before','After','Adjacent'};
Daniel@0:         sum.default = 'Before';
Daniel@0:     option.sum = sum;
Daniel@0:     
Daniel@0:         m.key = 'Total';
Daniel@0:         m.type = 'Integer';
Daniel@0:         m.default = 1;
Daniel@0:     option.m = m;
Daniel@0:         
Daniel@0:         thr.key = 'Contrast';
Daniel@0:         thr.type = 'Integer';
Daniel@0:         thr.default = 0.01; % Not same default as in mirtempo
Daniel@0:     option.thr = thr;
Daniel@0: 
Daniel@0: specif.option = option;
Daniel@0: 
Daniel@0: varargout = mirfunction(@mirpulseclarity,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: %% Initialisation
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: %if isframed(x)
Daniel@0: %    warning('WARNING IN MIRPULSECLARITY: The input should not be already decomposed into frames.');
Daniel@0: %    disp(['Suggestion: Use the ''Frame'' option instead.'])
Daniel@0: %end
Daniel@0: if iscell(x)
Daniel@0:     x = x{1};
Daniel@0: end
Daniel@0: if isamir(x,'mirautocor')
Daniel@0:     type = {'mirscalar','mirautocor'};
Daniel@0: elseif length(option.model) > 1
Daniel@0:     a = x;
Daniel@0:     type = {'mirscalar'};
Daniel@0:     for m = 1:length(option.model)
Daniel@0:         if option.frame.length.val
Daniel@0:             y = mirpulseclarity(a,'Model',option.model(m),...
Daniel@0:                                   'Frame',option.frame.length.val,...
Daniel@0:                                           option.frame.length.unit,...
Daniel@0:                                           option.frame.hop.val,...
Daniel@0:                                           option.frame.hop.unit);
Daniel@0:         else
Daniel@0:             y = mirpulseclarity(a,'Model',option.model(m));
Daniel@0:         end
Daniel@0:         if m == 1
Daniel@0:             x = y;
Daniel@0:         else
Daniel@0:             x = x + y;
Daniel@0:         end
Daniel@0:     end
Daniel@0: else
Daniel@0:     if option.model
Daniel@0:         switch option.model
Daniel@0:             case 1
Daniel@0:             case 2
Daniel@0:                 option.envmeth = 'Filter';
Daniel@0:                 option.fbtype = 'Gammatone';
Daniel@0:                 option.mu = 0;
Daniel@0:                 option.r = 0;
Daniel@0:                 option.lambda = .8;
Daniel@0:                 option.sum = 'After';
Daniel@0:         end
Daniel@0:     end
Daniel@0:     if length(option.stratg)>7 && strcmpi(option.stratg(end-6:end),'Autocor')
Daniel@0:         if (strcmpi(option.stratg,'MaxAutocor') || ...
Daniel@0:             strcmpi(option.stratg,'MinAutocor') || ...
Daniel@0:             strcmpi(option.stratg,'EntropyAutocor'))
Daniel@0:             option.m = 0;
Daniel@0:         end
Daniel@0:         if strcmpi(option.stratg,'MinAutocor')
Daniel@0:             option.enh = 0;
Daniel@0:         end
Daniel@0:         if option.frame.length.val
Daniel@0:             [t,x] = mirtempo(x,option.fea,'Method',option.envmeth,...
Daniel@0:                        option.band,...
Daniel@0:                        'Sum',option.sum,'Enhanced',option.enh,...
Daniel@0:                        'Resonance',option.r,'Smooth',option.aver,...
Daniel@0:                        'HalfwaveDiff',option.diffhwr,...
Daniel@0:                        'Lambda',option.lambda,...
Daniel@0:                        'Frame',option.frame.length.val,...
Daniel@0:                                option.frame.length.unit,...
Daniel@0:                                option.frame.hop.val,...
Daniel@0:                                option.frame.hop.unit,...
Daniel@0:                        'FilterbankType',option.fbtype,...
Daniel@0:                        'FilterType',option.ftype,...
Daniel@0:                        'Filterbank',option.fb,'Mu',option.mu,...
Daniel@0:                        'Log',option.log,...
Daniel@0:                        'Inc',option.inc,'Halfwave',option.hw,...
Daniel@0:                        'Median',option.median(1),option.median(2),...
Daniel@0:                        'Min',option.mi,'Max',option.ma,...
Daniel@0:                        'Total',option.m,'Contrast',option.thr);
Daniel@0:         else
Daniel@0:             [t,x] = mirtempo(x,option.fea,'Method',option.envmeth,...
Daniel@0:                        option.band,...
Daniel@0:                        'Sum',option.sum,'Enhanced',option.enh,...
Daniel@0:                        'Resonance',option.r,'Smooth',option.aver,...
Daniel@0:                        'HalfwaveDiff',option.diffhwr,...
Daniel@0:                        'Lambda',option.lambda,...
Daniel@0:                        'FilterbankType',option.fbtype,...
Daniel@0:                        'FilterType',option.ftype,...
Daniel@0:                        'Filterbank',option.fb,'Mu',option.mu,...
Daniel@0:                        'Log',option.log,...
Daniel@0:                        'Inc',option.inc,'Halfwave',option.hw,...
Daniel@0:                        'Median',option.median(1),option.median(2),...
Daniel@0:                        'Min',option.mi,'Max',option.ma,...
Daniel@0:                        'Total',option.m,'Contrast',option.thr);
Daniel@0:         end
Daniel@0:         type = {'mirscalar','mirautocor'};
Daniel@0:     elseif strcmpi(option.stratg,'ExtremEnvelop')
Daniel@0:         x = mironsets(x,'Filterbank',option.fb);
Daniel@0:         type = {'mirscalar','mirenvelope'};
Daniel@0:     elseif strcmpi(option.stratg,'Attack')
Daniel@0:         x = mirattackslope(x,option.slope);
Daniel@0:         type = {'mirscalar','mirenvelope'};
Daniel@0: %    elseif strcmpi(option.stratg,'AttackDiff')
Daniel@0: %        type = {'mirscalar','mirenvelope'};
Daniel@0:     elseif strcmpi(option.stratg,'Articulation')
Daniel@0:         x = mirlowenergy(x,'ASR');
Daniel@0:         type = {'mirscalar','mirscalar'};
Daniel@0:     else
Daniel@0:         type = {'mirscalar','miraudio'};
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: %% Main function
Daniel@0: 
Daniel@0: function o = main(a,option,postoption)
Daniel@0: if option.model == 2
Daniel@0:     option.stratg = 'InterfAutocor';
Daniel@0: end
Daniel@0: if isa(a,'mirscalar') && not(strcmpi(option.stratg,'Attack')) % not very nice test... to improve.
Daniel@0:     o = {a};
Daniel@0:     return
Daniel@0: end
Daniel@0: if option.m == 1 && ...
Daniel@0:         (strcmpi(option.stratg,'InterfAutocor') || ...
Daniel@0:          strcmpi(option.stratg,'MeanPeaksAutocor'))
Daniel@0:     option.m = Inf;
Daniel@0: end
Daniel@0: if iscell(a)
Daniel@0:     a = a{1};
Daniel@0: end
Daniel@0: if strcmpi(option.stratg,'MaxAutocor')
Daniel@0:     d = get(a,'Data');
Daniel@0:     rc = mircompute(@max,d);
Daniel@0: elseif strcmpi(option.stratg,'MinAutocor')
Daniel@0:     d = get(a,'Data');
Daniel@0:     rc = mircompute(@minusmin,d);
Daniel@0: elseif strcmpi(option.stratg,'MeanPeaksAutocor')
Daniel@0:     m = get(a,'PeakVal');
Daniel@0:     rc = mircompute(@meanpeaks,m);    
Daniel@0: elseif strcmpi(option.stratg,'KurtosisAutocor')
Daniel@0:     a = mirpeaks(a,'Extract','Total',option.m,'NoBegin','NoEnd');
Daniel@0:     k = mirkurtosis(a);
Daniel@0:     %d = get(k,'Data');
Daniel@0:     %rc = mircompute(@meanpeaks,d);
Daniel@0:     rc = mirmean(k);
Daniel@0: elseif strcmpi(option.stratg,'EntropyAutocor')
Daniel@0:     rc = mirentropy(a);
Daniel@0: elseif strcmpi(option.stratg,'InterfAutocor')
Daniel@0:     a = mirpeaks(a,'Total',option.m,'NoBegin','NoEnd');
Daniel@0:     m = get(a,'PeakVal');
Daniel@0:     p = get(a,'PeakPosUnit');
Daniel@0:     rc = mircompute(@interf,m,p);
Daniel@0: elseif strcmpi(option.stratg,'TempoAutocor')
Daniel@0:     a = mirpeaks(a,'Total',1,'NoBegin','NoEnd');
Daniel@0:     p = get(a,'PeakPosUnit');
Daniel@0:     rc = mircompute(@tempo,p);
Daniel@0: elseif strcmpi(option.stratg,'ExtremEnvelop')
Daniel@0:     a = mirenvelope(a,'Normal');
Daniel@0:     p = mirpeaks(a,'Order','Abscissa');
Daniel@0:     p = get(p,'PeakPreciseVal');
Daniel@0:     n = mirpeaks(a,'Valleys','Order','Abscissa');
Daniel@0:     n = get(n,'PeakPreciseVal');
Daniel@0:     rc = mircompute(@shape,p,n);
Daniel@0: elseif strcmpi(option.stratg,'Attack')
Daniel@0:     rc = mirmean(a);
Daniel@0: %elseif strcmpi(option.stratg,'AttackDiff')
Daniel@0: %    a = mirpeaks(a);
Daniel@0: %    m = get(a,'PeakVal');
Daniel@0: %    rc = mircompute(@meanpeaks,m);    
Daniel@0: elseif strcmpi(option.stratg,'Articulation')
Daniel@0:     rc = a;
Daniel@0: end
Daniel@0: 
Daniel@0: if iscell(rc)
Daniel@0:     pc = mirscalar(a,'Data',rc,'Title','Pulse clarity');
Daniel@0: else
Daniel@0:     pc = set(rc,'Title',['Pulse clarity (',get(rc,'Title'),')']);
Daniel@0: end
Daniel@0: 
Daniel@0: if option.model
Daniel@0:     switch option.model
Daniel@0:         case 1
Daniel@0:             alpha = 0;
Daniel@0:             beta = 2.2015;
Daniel@0:             lambda = .1;
Daniel@0:         case 2
Daniel@0:             alpha = 0;
Daniel@0:             beta = 3.5982;
Daniel@0:             lambda = 1.87;
Daniel@0:     end
Daniel@0:     if not(lambda == 0)
Daniel@0:         pc = (pc+alpha)^lambda * beta;
Daniel@0:     else
Daniel@0:         pc = log(pc+alpha) * beta;
Daniel@0:     end
Daniel@0:     title = ['Pulse clarity (Model ',num2str(option.model),')'];
Daniel@0:     pc = set(pc,'Title',title);
Daniel@0: end
Daniel@0: 
Daniel@0: o = {pc a};
Daniel@0: 
Daniel@0: 
Daniel@0: %% Routines
Daniel@0: 
Daniel@0: function r  = shape(p,n)
Daniel@0: p = p{1};
Daniel@0: n = n{1};
Daniel@0: if length(p)>length(n)
Daniel@0:     d = sum(p(1:end-1) - n) + sum(p(2:end) - n);
Daniel@0:     r  = d/(2*length(n));
Daniel@0: elseif length(p)<length(n)
Daniel@0:     d = sum(p - n(1:end-1)) + sum(p - n(2:end));
Daniel@0:     r  = d/(2*length(p));
Daniel@0: else
Daniel@0:     d = sum(p(2:end) - n(1:end-1)) + sum(p(1:end-1) - n(2:end));
Daniel@0:     r  = d/(2*(length(p)-1));
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function rc = minusmin(ac)
Daniel@0: rc = -min(ac);
Daniel@0: 
Daniel@0: 
Daniel@0: function rc = meanpeaks(ac)
Daniel@0: rc = zeros(1,length(ac));
Daniel@0: for j = 1:length(ac)
Daniel@0:     if isempty(ac{j})
Daniel@0:         rc(j) = NaN;
Daniel@0:     else
Daniel@0:         rc(j) = mean(ac{j});
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function rc = interf(mk,pk)
Daniel@0: rc = zeros(size(mk));
Daniel@0: for j = 1:size(mk,3)
Daniel@0:     for i = 1:size(mk,2)
Daniel@0:         pij = pk{1,i,j};
Daniel@0:         mij = mk{1,i,j};
Daniel@0:         if isempty(pij)
Daniel@0:             rc(1,i,j) = 0;
Daniel@0:         else
Daniel@0:             high = max(pij(2:end),pij(1));
Daniel@0:             low = min(pij(2:end),pij(1));
Daniel@0:             quo = rem(high,low)./low;
Daniel@0:             nomult = quo>.15 & quo<.85;
Daniel@0:             fij = mij(2:end)/mij(1) .*nomult;
Daniel@0:             fij(fij<0) = 0;
Daniel@0:             rc(1,i,j) = exp(-sum(fij)/4); % Pulsations that are not in integer ratio
Daniel@0:                                           % with dominant pulse decrease clarity
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function rc = tempo(pk)
Daniel@0: rc = zeros(size(pk));
Daniel@0: for j = 1:size(pk,3)
Daniel@0:     for i = 1:size(pk,2)
Daniel@0:         pij = pk{1,i,j};
Daniel@0:         if isempty(pij)
Daniel@0:             rc(1,i,j) = 0;
Daniel@0:         else
Daniel@0:             rc(1,i,j) = exp(-pij(1)/4)/exp(-.33/4); % Fast dominant pulse
Daniel@0:                                                     % increases clarity
Daniel@0:         end
Daniel@0:     end
Daniel@0: end