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