wolffd@0: function varargout = mirpitch(orig,varargin)
wolffd@0: %   p = mirpitch(x) evaluates the pitch frequencies (in Hz).
wolffd@0: %   Specification of the method(s) for pitch estimation (these methods can
wolffd@0: %       be combined):
wolffd@0: %       mirpitch(...,'Autocor') computes an autocorrelation function
wolffd@0: %           (Default method)
wolffd@0: %           mirpitch(...'Enhanced',a) computes enhanced autocorrelation
wolffd@0: %               (see help mirautocor)
wolffd@0: %              toggled on by default
wolffd@0: %           mirpitch(...,'Compress',k) performs magnitude compression
wolffd@0: %               (see help mirautocor)
wolffd@0: %           mirpitch(...,fb) specifies a type of filterbank.
wolffd@0: %               Possible values:
wolffd@0: %                   fb = 'NoFilterBank': no filterbank decomposition
wolffd@0: %                   fb = '2Channels' (default value)
wolffd@0: %                   fb = 'Gammatone' 
wolffd@0: %       mirpitch(...,'AutocorSpectrum') computes the autocorrelation of
wolffd@0: %           the FFT spectrum
wolffd@0: %       mirpitch(...,'Cepstrum') computes the cepstrum
wolffd@0: %       Alternatively, an autocorrelation or a cepstrum can be directly
wolffd@0: %           given as first argument of the mirpitch function.
wolffd@0: %   Peak picking options:
wolffd@0: %       mirpitch(...,'Total',m) selects the m best pitches.
wolffd@0: %           Default value: m = Inf, no limit is set concerning the number
wolffd@0: %           of pitches to be detected.
wolffd@0: %       mirpitch(...,'Mono') corresponds to morpitch(...,'Total',1)
wolffd@0: %       mirpitch(...,'Min',mi) indicates the lowest frequency taken into
wolffd@0: %           consideration.
wolffd@0: %           Default value: 75 Hz. (Praat)
wolffd@0: %       mirpitch(...,'Max',ma) indicates the highest frequency taken into
wolffd@0: %           consideration. 
wolffd@0: %           Default value: 2400 Hz. Because there seems to be some problems
wolffd@0: %           with higher frequency, due probably to the absence of 
wolffd@0: %           pre-whitening in our implementation of Tolonen and Karjalainen
wolffd@0: %           approach (used by default, cf. below).
wolffd@0: %       mirpitch(...,'Contrast',thr) specifies a threshold value.
wolffd@0: %           (see help peaks)
wolffd@0: %           Default value: thr = .1
wolffd@0: %       mirpitch(...,'Order',o) specifies the ordering for the peak picking.
wolffd@0: %           Default value: o = 'Amplitude'.
wolffd@0: %       Alternatively, the result of a mirpeaks computation can be directly
wolffd@0: %           given as first argument of the mirpitch function.
wolffd@0: %   Post-processing options:
wolffd@0: %       mirpitch(...,'Sum','no') does not sum back the channels at the end 
wolffd@0: %           of the computation. The resulting pitch information remains
wolffd@0: %           therefore decomposed into several channels.
wolffd@0: %       mirpitch(...,'Median') performs a median filtering of the pitch
wolffd@0: %           curve. When several pitches are extracted in each frame, the
wolffd@0: %           pitch curve contains the best peak of each successive frame.
wolffd@0: %       mirpitch(...,'Stable',th,n) remove pitch values when the difference 
wolffd@0: %           (or more precisely absolute logarithmic quotient) with the
wolffd@0: %           n precedent frames exceeds the threshold th. 
wolffd@0: %           if th is not specified, the default value .1 is used
wolffd@0: %           if n is not specified, the default value 3 is used
wolffd@0: %       mirpitch(...'Reso',r) removes peaks whose distance to one or
wolffd@0: %           several higher peaks is lower than a given threshold.
wolffd@0: %           Possible value for the threshold r:
wolffd@0: %               'SemiTone': ratio between the two peak positions equal to
wolffd@0: %                   2^(1/12)
wolffd@0: %       mirpitch(...,'Frame',l,h) orders a frame decomposition of window
wolffd@0: %           length l (in seconds) and hop factor h, expressed relatively to
wolffd@0: %           the window length. For instance h = 1 indicates no overlap.
wolffd@0: %           Default values: l = 46.4 ms and h = 10 ms (Tolonen and
wolffd@0: %           Karjalainen, 2000)
wolffd@0: %   Preset model:
wolffd@0: %       mirpitch(...,'Tolonen') implements (part of) the model proposed in
wolffd@0: %           (Tolonen & Karjalainen, 2000). It is equivalent to
wolffd@0: %           mirpitch(...,'Enhanced',2:10,'Generalized',.67,'2Channels')
wolffd@0: %   [p,a] = mirpitch(...) also displays the result of the method chosen for
wolffd@0: %       pitch estimation, and shows in particular the peaks corresponding
wolffd@0: %       to the pitch values.
wolffd@0: %   p = mirpitch(f,a,<r>) creates a mirpitch object based on the frequencies
wolffd@0: %       specified in f and the related amplitudes specified in a, using a
wolffd@0: %       frame sampling rate of r Hz (set by default to 100 Hz).
wolffd@0: %
wolffd@0: %   T. Tolonen, M. Karjalainen, "A Computationally Efficient Multipitch 
wolffd@0: %       Analysis Model", IEEE TRANSACTIONS ON SPEECH AND AUDIO PROCESSING,
wolffd@0: %       VOL. 8, NO. 6, NOVEMBER 2000
wolffd@0: 
wolffd@0:         ac.key = 'Autocor';
wolffd@0:         ac.type = 'Boolean';
wolffd@0:         ac.default = 0;
wolffd@0:     option.ac = ac;
wolffd@0:     
wolffd@0:             enh.key = 'Enhanced';
wolffd@0:             enh.type = 'Integer';
wolffd@0:             enh.default = 2:10;
wolffd@0:         option.enh = enh;
wolffd@0: 
wolffd@0:             filtertype.type = 'String';
wolffd@0:             filtertype.choice = {'NoFilterBank','2Channels','Gammatone'};
wolffd@0:             filtertype.default = '2Channels';
wolffd@0:         option.filtertype = filtertype;
wolffd@0: 
wolffd@0:             gener.key = {'Generalized','Compress'};
wolffd@0:             gener.type = 'Integer';
wolffd@0:             gener.default = .5;
wolffd@0:         option.gener = gener;
wolffd@0: 
wolffd@0:         as.key = 'AutocorSpectrum';
wolffd@0:         as.type = 'Boolean';
wolffd@0:         as.default = 0;
wolffd@0:     option.as = as;
wolffd@0:     
wolffd@0:         s.key = 'Spectrum';
wolffd@0:         s.type = 'Boolean';
wolffd@0:         s.default = 0;
wolffd@0:     option.s = s;
wolffd@0:         
wolffd@0:         ce.key = 'Cepstrum';
wolffd@0:         ce.type = 'Boolean';
wolffd@0:         ce.default = 0;
wolffd@0:     option.ce = ce;
wolffd@0:         
wolffd@0: %% peak picking options
wolffd@0: 
wolffd@0:         m.key = 'Total';
wolffd@0:         m.type = 'Integer';
wolffd@0:         m.default = Inf;
wolffd@0:     option.m = m;
wolffd@0:     
wolffd@0:         multi.key = 'Multi';
wolffd@0:         multi.type = 'Boolean';
wolffd@0:         multi.default = 0;
wolffd@0:     option.multi = multi;
wolffd@0: 
wolffd@0:         mono.key = 'Mono';
wolffd@0:         mono.type = 'Boolean';
wolffd@0:         mono.default = 0;
wolffd@0:     option.mono = mono;
wolffd@0: 
wolffd@0:         mi.key = 'Min';
wolffd@0:         mi.type = 'Integer';
wolffd@0:         mi.default = 75;
wolffd@0:     option.mi = mi;
wolffd@0:         
wolffd@0:         ma.key = 'Max';
wolffd@0:         ma.type = 'Integer';
wolffd@0:         ma.default = 2400;
wolffd@0:     option.ma = ma;
wolffd@0:         
wolffd@0:         thr.key = 'Contrast';
wolffd@0:         thr.type = 'Integer';
wolffd@0:         thr.default = .1;
wolffd@0:     option.thr = thr;
wolffd@0:     
wolffd@0:         order.key = 'Order';
wolffd@0:         order.type = 'String';
wolffd@0:         order.choice = {'Amplitude','Abscissa'};
wolffd@0:         order.default = 'Amplitude';
wolffd@0:     option.order = order;    
wolffd@0: 
wolffd@0:         reso.key = 'Reso';
wolffd@0:         reso.type = 'String';
wolffd@0:         reso.choice = {0,'SemiTone'};
wolffd@0:         reso.default = 0;
wolffd@0:     option.reso = reso;
wolffd@0:         
wolffd@0:         track.key = 'Track';        % Not used yet
wolffd@0:         track.type = 'Boolean';
wolffd@0:         track.default = 0;
wolffd@0:     option.track = track;
wolffd@0: 
wolffd@0: %% post-processing options
wolffd@0:         
wolffd@0:         stable.key = 'Stable';
wolffd@0:         stable.type = 'Integer';
wolffd@0:         stable.number = 2;
wolffd@0:         stable.default = [Inf 0];
wolffd@0:         stable.keydefault = [.1 3];
wolffd@0:     option.stable = stable;
wolffd@0:     
wolffd@0:         median.key = 'Median';
wolffd@0:         median.type = 'Integer';
wolffd@0:         median.default = 0;
wolffd@0:         median.keydefault = .1;
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 = [0 0];
wolffd@0:         frame.keydefault = [NaN NaN];
wolffd@0:     option.frame = frame;
wolffd@0: 
wolffd@0:         sum.key = 'Sum';
wolffd@0:         sum.type = 'Boolean';
wolffd@0:         sum.default = 1;
wolffd@0:     option.sum = sum;
wolffd@0:     
wolffd@0: %% preset model
wolffd@0: 
wolffd@0:         tolo.key = 'Tolonen';
wolffd@0:         tolo.type = 'Boolean';
wolffd@0:         tolo.default = 0;
wolffd@0:     option.tolo = tolo;
wolffd@0:     
wolffd@0: specif.option = option;
wolffd@0: specif.chunkframebefore = 1;
wolffd@0: 
wolffd@0: if isnumeric(orig)
wolffd@0:     if nargin<3
wolffd@0:         f = 100;
wolffd@0:     else
wolffd@0:         f = varargin{2};
wolffd@0:     end
wolffd@0:     fp = (0:size(orig,1)-1)/f;
wolffd@0:     fp = [fp;fp+1/f];
wolffd@0:     p.amplitude = {{varargin{1}'}};
wolffd@0:     s = mirscalar([],'Data',{{orig'}},'Title','Pitch','Unit','Hz',...
wolffd@0:                      'FramePos',{{fp}},'Sampling',f,'Name',{inputname(1)});
wolffd@0:     p = class(p,'mirpitch',s);
wolffd@0:     varargout = {p};
wolffd@0: else
wolffd@0:     varargout = mirfunction(@mirpitch,orig,varargin,nargout,specif,@init,@main);
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: function [y type] = init(orig,option)
wolffd@0: if option.tolo
wolffd@0:     option.enh = 2:10;
wolffd@0:     option.gener = .67;
wolffd@0:     option.filtertype = '2Channels';
wolffd@0: end
wolffd@0: if not(option.ac) && not(option.as) && not(option.ce) && not(option.s)
wolffd@0:     option.ac = 1;
wolffd@0: end
wolffd@0: if isnan(option.frame.length.val)
wolffd@0:     option.frame.length.val = .0464;
wolffd@0: end
wolffd@0: if isnan(option.frame.hop.val)
wolffd@0:     option.frame.hop.val = .01;
wolffd@0:     option.frame.hop.unit = 's';
wolffd@0: end
wolffd@0: if isamir(orig,'mirscalar') || haspeaks(orig)
wolffd@0:     y = orig;
wolffd@0: else
wolffd@0:     if isamir(orig,'mirautocor')
wolffd@0:         y = mirautocor(orig,'Min',option.mi,'Hz','Max',option.ma,'Hz','Freq');
wolffd@0:     elseif isamir(orig,'mircepstrum')
wolffd@0:         y = orig;
wolffd@0:     elseif isamir(orig,'mirspectrum')
wolffd@0:         if not(option.as) && not(option.ce) && not(option.s)
wolffd@0:             option.ce = 1;
wolffd@0:         end
wolffd@0:         if option.as
wolffd@0:             y = mirautocor(orig,...
wolffd@0:                             'Min',option.mi,'Hz','Max',option.ma,'Hz');
wolffd@0:         end
wolffd@0:         if option.ce
wolffd@0:             ce = mircepstrum(orig,'freq',...
wolffd@0:                             'Min',option.mi,'Hz','Max',option.ma,'Hz');
wolffd@0:             if option.as
wolffd@0:                 y = y*ce;
wolffd@0:             else
wolffd@0:                 y = ce;
wolffd@0:             end
wolffd@0:         end
wolffd@0:     else
wolffd@0:         if option.ac
wolffd@0:             x = orig;
wolffd@0:             if not(strcmpi(option.filtertype,'NoFilterBank'))
wolffd@0:                 x = mirfilterbank(x,option.filtertype);
wolffd@0:             end
wolffd@0:             x = mirframenow(x,option);
wolffd@0:             y = mirautocor(x,'Generalized',option.gener,...
wolffd@0:                                 'Min',option.mi,'Hz','Max',option.ma,'Hz');
wolffd@0:             if option.sum
wolffd@0:                 y = mirsummary(y);
wolffd@0:             end
wolffd@0:             y = mirautocor(y,'Enhanced',option.enh,'Freq');
wolffd@0:         end
wolffd@0:         if option.as || option.ce || option.s
wolffd@0:             x = mirframenow(orig,option);
wolffd@0:             y = mirspectrum(x);
wolffd@0:             if option.as
wolffd@0:                 as = mirautocor(y,...
wolffd@0:                                 'Min',option.mi,'Hz','Max',option.ma,'Hz');
wolffd@0:                 if option.ac
wolffd@0:                     y = y*as;
wolffd@0:                 else
wolffd@0:                     y = as;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             if option.ce
wolffd@0:                 ce = mircepstrum(y,'freq',...
wolffd@0:                                 'Min',option.mi,'Hz','Max',option.ma,'Hz');
wolffd@0:                 if option.ac || option.as
wolffd@0:                     y = y*ce;
wolffd@0:                 else
wolffd@0:                     y = ce;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: type = {'mirpitch',mirtype(y)};
wolffd@0:     
wolffd@0: 
wolffd@0: function o = main(x,option,postoption)
wolffd@0: if option.multi && option.m == 1
wolffd@0:     option.m = Inf;
wolffd@0: end
wolffd@0: if option.mono && option.m == Inf
wolffd@0:     option.m = 1;
wolffd@0: end
wolffd@0: if iscell(x)
wolffd@0:     x = x{1};
wolffd@0: end
wolffd@0: if not(isa(x,'mirpitch'))
wolffd@0:     x = mirpeaks(x,'Total',option.m,'Track',option.track,...
wolffd@0:                    'Contrast',option.thr,'Threshold',.4,...
wolffd@0:                    'Reso',option.reso,'NoBegin','NoEnd',...
wolffd@0:                    'Order',option.order);
wolffd@0: end
wolffd@0: if isa(x,'mirscalar')
wolffd@0:     pf = get(x,'Data');
wolffd@0: else
wolffd@0:     pf = get(x,'PeakPrecisePos');
wolffd@0:     pa = get(x,'PeakPreciseVal');
wolffd@0: end
wolffd@0: fp = get(x,'FramePos');
wolffd@0: if option.stable(1) < Inf
wolffd@0:     for i = 1:length(pf)
wolffd@0:         for j = 1:length(pf{i})
wolffd@0:             for k = 1:size(pf{i}{j},3)
wolffd@0:                 for l = size(pf{i}{j},2):-1:option.stable(2)+1
wolffd@0:                     for m = length(pf{i}{j}{1,l,k}):-1:1
wolffd@0:                         found = 0;
wolffd@0:                         for h = 1:option.stable(2)
wolffd@0:                             for n = 1:length(pf{i}{j}{1,l-h,k})
wolffd@0:                                 if abs(log10(pf{i}{j}{1,l,k}(m) ...
wolffd@0:                                             /pf{i}{j}{1,l-h,k}(n))) ...
wolffd@0:                                        < option.stable(1)
wolffd@0:                                     found = 1;
wolffd@0:                                 end
wolffd@0:                             end
wolffd@0:                         end
wolffd@0:                         if not(found)
wolffd@0:                             pf{i}{j}{1,l,k}(m) = [];
wolffd@0:                         end
wolffd@0:                     end
wolffd@0:                     pf{i}{j}{1,1,k} = zeros(1,0);
wolffd@0:                 end
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if option.median
wolffd@0:     sr = get(x,'Sampling');
wolffd@0:     for i = 1:length(pf)
wolffd@0:         for j = 1:length(pf{i})
wolffd@0:             if size(fp{i}{j},2) > 1
wolffd@0:                 npf = zeros(size(pf{i}{j}));
wolffd@0:                 for k = 1:size(pf{i}{j},3)
wolffd@0:                     for l = 1:size(pf{i}{j},2)
wolffd@0:                         if isempty(pf{i}{j}{1,l,k})
wolffd@0:                             npf(1,l,k) = NaN;
wolffd@0:                         else
wolffd@0:                             npf(1,l,k) = pf{i}{j}{1,l,k}(1);
wolffd@0:                         end
wolffd@0:                     end
wolffd@0:                 end
wolffd@0:                 pf{i}{j} = medfilt1(npf,...
wolffd@0:                      round(option.median/(fp{i}{j}(1,2)-fp{i}{j}(1,1))));
wolffd@0:             end
wolffd@0:         end
wolffd@0:     end
wolffd@0: end
wolffd@0: if isa(x,'mirscalar')
wolffd@0:     p.amplitude = 0;
wolffd@0: else
wolffd@0:     p.amplitude = pa;
wolffd@0: end
wolffd@0: s = mirscalar(x,'Data',pf,'Title','Pitch','Unit','Hz');
wolffd@0: p = class(p,'mirpitch',s);
wolffd@0: o = {p,x};