Daniel@0: function varargout = mircentroid(x,varargin)
Daniel@0: %   c = mircentroid(x) calculates the centroid (or center of gravity) of x.
Daniel@0: %   x can be either:
Daniel@0: %       - a spectrum (spectral centroid),
Daniel@0: %       - an envelope (temporal centroid)
Daniel@0: %       - a histogram,
Daniel@0: %       - or any data. Only the positive ordinates of the data are taken
Daniel@0: %           into consideration.
Daniel@0: %   c = mircentroid(x,'Peaks') calculates the centroid of the peaks only.
Daniel@0: 
Daniel@0: % Beauchamp 1982 version?
Daniel@0: 
Daniel@0:         peaks.key = 'Peaks';
Daniel@0:         peaks.type = 'String';
Daniel@0:         peaks.choice = {0,'NoInterpol','Interpol'};
Daniel@0:         peaks.default = 0;
Daniel@0:         peaks.keydefault = 'NoInterpol';
Daniel@0:     option.peaks = peaks;
Daniel@0:     
Daniel@0: specif.option = option;
Daniel@0: 
Daniel@0: varargout = mirfunction(@mircentroid,x,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: if not(isamir(x,'mirdata')) || isamir(x,'miraudio')
Daniel@0:     x = mirspectrum(x);
Daniel@0: end
Daniel@0: type = 'mirscalar';
Daniel@0: 
Daniel@0: 
Daniel@0: function c = main(x,option,postoption)
Daniel@0: if iscell(x)
Daniel@0:     x = x{1};
Daniel@0: end
Daniel@0: if option.peaks
Daniel@0:     if strcmpi(option.peaks,'Interpol')
Daniel@0:         pt = get(x,'PeakPrecisePos');
Daniel@0:         pv = get(x,'PeakPreciseVal');
Daniel@0:     else
Daniel@0:         pt = get(x,'PeakPos');
Daniel@0:         pv = get(x,'PeakVal');
Daniel@0:     end
Daniel@0:     cx = cell(1,length(pt));
Daniel@0:     for h = 1:length(pt)
Daniel@0:         cx{h} = cell(1,length(pt{h}));
Daniel@0:         for i = 1:length(pt{h})
Daniel@0:             pti = pt{h}{i};
Daniel@0:             pvi = pv{h}{i};
Daniel@0:             %if isempty(pti)
Daniel@0: 
Daniel@0:             nfr = size(pti,2);
Daniel@0:             nbd = size(pti,3);
Daniel@0:             ci = zeros(1,nfr,nbd);
Daniel@0:             for j = 1:nfr
Daniel@0:                 for k = 1:nbd
Daniel@0:                     ptk = pti{1,j,k};
Daniel@0:                     pvk = pvi{1,j,k};
Daniel@0:                     sk = sum(pvk);
Daniel@0:                     ci(1,j,k) = sum(ptk.*pvk) ./ sk;
Daniel@0:                 end
Daniel@0:             end
Daniel@0:             cx{h}{i} = ci;
Daniel@0:         end
Daniel@0:     end
Daniel@0: else
Daniel@0:     cx = peaksegments(@centroid,get(x,'Data'),get(x,'Pos'));
Daniel@0: end
Daniel@0: if isa(x,'mirspectrum')
Daniel@0:     t = 'Spectral centroid';
Daniel@0: elseif isa(x,'mirenvelope')
Daniel@0:     t = 'Temporal centroid';
Daniel@0: else
Daniel@0:     t = ['centroid of ',get(x,'Title')];
Daniel@0: end
Daniel@0: c = mirscalar(x,'Data',cx,'Title',t);
Daniel@0: 
Daniel@0: 
Daniel@0: function c = centroid(d,p)
Daniel@0: c = (p'*d) ./ sum(d);