Daniel@0: function varargout = mirinharmonicity(orig,varargin)
Daniel@0: %   ih = mirinharmonicity(x) estimates the inharmonicity of x, i.e., the
Daniel@0: %       amount of partials that are not multiples of the fundamental
Daniel@0: %       frequency.
Daniel@0: %       x can be either an audio file, a miraudio or a mirspectrum object.
Daniel@0: %   WARNING: This function presupposes that there is only one fundamental
Daniel@0: %       frequency.
Daniel@0: %   Optional argument:
Daniel@0: %       mirinharmonicity(...,'f0',f) bases the computation of the
Daniel@0: %           inharmonicity on the fundamental frequency indicated by f.
Daniel@0: %           Default value: f = mirpitch(...,'Mono')
Daniel@0: %   [ih,s] = mirinharmonicity(x) also displays the spectrum used for the
Daniel@0: %       computation of the inharmonicity.
Daniel@0: %   [ih,s,p] = mirinharmonicity(x) also displays the result of the
Daniel@0: %       estimation of the fundamental frequency.
Daniel@0: 
Daniel@0:         f0.key = 'f0';
Daniel@0:         f0.default = [];
Daniel@0:     option.f0 = f0;
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 = [.1 .125];
Daniel@0:     option.frame = frame;
Daniel@0:         
Daniel@0: specif.option = option;
Daniel@0:      
Daniel@0: varargout = mirfunction(@mirinharmonicity,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [i type] = init(x,option)
Daniel@0: if isamir(x,'miraudio')
Daniel@0:     if option.frame.length.val
Daniel@0:         s = mirspectrum(x,'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:         s = mirspectrum(x);
Daniel@0:     end
Daniel@0: else
Daniel@0:     s = x;
Daniel@0: end
Daniel@0: if isempty(option.f0)
Daniel@0:     if option.frame.length.val
Daniel@0:         p = mirpitch(x,'Mono','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:         p = mirpitch(x,'Mono');
Daniel@0:     end
Daniel@0: else
Daniel@0:     p = option.f0;
Daniel@0: end
Daniel@0: i = {s,p};
Daniel@0: type = {'mirscalar','mirspectrum','mirscalar'};
Daniel@0: 
Daniel@0: 
Daniel@0: function ih = main(x,option,postoption)
Daniel@0: if isa(x{2},'mirdesign')
Daniel@0:     x = x{1};
Daniel@0: end
Daniel@0: s = x{1};
Daniel@0: p = x{2};
Daniel@0: if iscell(p)
Daniel@0:     p = p{1};
Daniel@0: end
Daniel@0: m = get(s,'Magnitude');
Daniel@0: f = get(s,'Frequency');
Daniel@0: fp1 = get(s,'FramePos');
Daniel@0: if isnumeric(p)
Daniel@0:     pf = {{{p}}};
Daniel@0: else
Daniel@0:     pf = get(p,'Data');
Daniel@0:     fp2 = get(p,'FramePos');
Daniel@0: end
Daniel@0: v = cell(1,length(m));
Daniel@0: for h = 1:length(m)
Daniel@0:     v{h} = cell(1,length(m{h}));
Daniel@0:     for i = 1:length(m{h})
Daniel@0:         mi = m{h}{i};
Daniel@0:         fi = f{h}{i};
Daniel@0:         pfi = pf{h}{i};
Daniel@0:         v{h}{i} = zeros(1,size(mi,2),size(mi,3));
Daniel@0:         if not(size(mi,2) == size(pfi,2))
Daniel@0:             beg = find(fp2{h}{i}(1,:) == fp1{h}{i}(1,1));
Daniel@0:             if isempty(beg) || (beg + size(mi,2)-1 > size(pfi,2))
Daniel@0:                 error('ERROR IN MIRINHARMONICITY: The ''f0'' argument should have the same frame decomposition than the main input.');
Daniel@0:             end
Daniel@0:             pfi = pfi(:,beg:beg+size(mi,2)-1);
Daniel@0:         end
Daniel@0:         for j = 1:size(mi,3)
Daniel@0:             for k = 1:size(mi,2)
Daniel@0:                 mk = mi(:,k,j);
Daniel@0:                 fk = fi(:,k,j);
Daniel@0:                 pfk = pfi(:,k);
Daniel@0:                 if isempty(pfk{1})
Daniel@0:                     v{h}{i}(1,k,j) = NaN;
Daniel@0:                 else
Daniel@0:                     r = fk/pfk{1}(1);
Daniel@0:                     rr = 2*abs(r-round(r));
Daniel@0:                     if isempty(rr)
Daniel@0:                         v{h}{i}(1,k,j) = NaN;
Daniel@0:                     else
Daniel@0:                         v{h}{i}(1,k,j) = sum(rr.*mk) / sum(mk);
Daniel@0:                     end
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0: end
Daniel@0: ih = mirscalar(s,'Data',v,'Title','Inharmonicity');
Daniel@0: if isa(p,'mirdata')
Daniel@0:     ih = {ih s p};
Daniel@0: else
Daniel@0:     ih = {ih s};
Daniel@0: end