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