Daniel@0: function varargout = mirroughness(x,varargin)
Daniel@0: %   r = mirroughness(x) calculates the roughness, or sensory dissonance,
Daniel@0: %           due to beating phenomenon between close frequency peaks. 
Daniel@0: %       The frequency components are supposed to remain sufficiently 
Daniel@0: %           constant throughout each frame of each audio file.
Daniel@0: %   r = mirroughness(...,'Contrast',c) specifies the contrast parameter
Daniel@0: %       used for peak picking (cf. mirpeaks).
Daniel@0: %       Default value: c = .01
Daniel@0: %   [r,s] = mirroughness(x) also displays the spectrum and its peaks, used
Daniel@0: %           for the computation of roughness.
Daniel@0: %   Optional arguments:
Daniel@0: %       Method used:
Daniel@0: %           mirroughness(...,'Sethares') (default): based on the summation
Daniel@0: %               of roughness between all pairs of sines (obtained through
Daniel@0: %               spectral peak-picking).
Daniel@0: %           mirroughness(...,'Vassilakis'): variant of 'Sethares' model
Daniel@0: %               with a more complex weighting (Vassilakis, 2001, Eq. 6.23).
Daniel@0: 
Daniel@0:         meth.type = 'String';
Daniel@0:         meth.choice = {'Sethares','Vassilakis'};
Daniel@0:         meth.default = 'Sethares';
Daniel@0:     option.meth = meth;
Daniel@0:     
Daniel@0:         cthr.key = 'Contrast';
Daniel@0:         cthr.type = 'Integer';
Daniel@0:         cthr.default = .01;
Daniel@0:     option.cthr = cthr;
Daniel@0:     
Daniel@0:         frame.key = 'Frame';
Daniel@0:         frame.type = 'Integer';
Daniel@0:         frame.number = 2;
Daniel@0:         frame.default = [.05 .5];
Daniel@0:     option.frame = frame;
Daniel@0:     
Daniel@0: specif.option = option;
Daniel@0: specif.defaultframelength = .05;
Daniel@0: specif.defaultframehop = .5;
Daniel@0: 
Daniel@0: 
Daniel@0: varargout = mirfunction(@mirroughness,x,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: if isamir(x,'miraudio') && not(isframed(x))
Daniel@0:     x = mirframenow(x,option);
Daniel@0: end
Daniel@0: x = mirspectrum(x);
Daniel@0: if not(haspeaks(x))
Daniel@0:     x = mirpeaks(x,'Contrast',option.cthr);
Daniel@0: end
Daniel@0: type = 'mirscalar';
Daniel@0: 
Daniel@0: 
Daniel@0: function r = main(p,option,postoption)
Daniel@0: if iscell(p)
Daniel@0:     p = p{1};
Daniel@0: end
Daniel@0: if strcmpi(option.meth,'Sethares') || strcmpi(option.meth,'Vassilakis')
Daniel@0:     pf = get(p,'PeakPosUnit');
Daniel@0:     pv = get(p,'PeakVal');
Daniel@0:     rg = cell(1,length(pf));
Daniel@0:     for h = 1:length(pf)
Daniel@0:         rg{h} = cell(1,length(pf{h}));
Daniel@0:         for i = 1:length(pf{h})
Daniel@0:             pfi = pf{h}{i};
Daniel@0:             pvi = pv{h}{i};
Daniel@0:             rg{h}{i} = zeros(1,length(pfi));
Daniel@0:             for k = 1:size(pfi,3)
Daniel@0:                 for j = 1:size(pfi,2)
Daniel@0:                     pfj = pfi{1,j,k};
Daniel@0:                     pvj = pvi{1,j,k};
Daniel@0:                     f1 = repmat(pfj,[1 length(pfj)]);
Daniel@0:                     f2 = repmat(pfj',[length(pfj) 1]);
Daniel@0:                     v1 = repmat(pvj,[1 length(pvj)]);
Daniel@0:                     v2 = repmat(pvj',[length(pvj) 1]);
Daniel@0:                     rj = plomp(f1,f2);
Daniel@0:                     if strcmpi(option.meth,'Sethares')
Daniel@0:                         rj = v1.*v2.*rj;
Daniel@0:                     elseif strcmpi(option.meth,'Vassilakis')
Daniel@0:                         rj = (v1.*v2).^.1.*.5.*(2*v2./(v1+v2)).^3.11.*rj;
Daniel@0:                     end
Daniel@0:                     rg{h}{i}(1,j,k) = sum(sum(rj));
Daniel@0:                 end
Daniel@0:             end
Daniel@0:         end
Daniel@0:     end
Daniel@0: else
Daniel@0: end    
Daniel@0: r = mirscalar(p,'Data',rg,'Title','Roughness');
Daniel@0: r = {r,p};
Daniel@0: 
Daniel@0: 
Daniel@0: function pd = plomp(f1, f2)
Daniel@0: % returns the dissonance of two pure tones at frequencies f1 & f2 Hz
Daniel@0: % according to the Plomp-Levelt curve (see Sethares)
Daniel@0:     b1 = 3.51;
Daniel@0:     b2 = 5.75;
Daniel@0:     xstar = .24;
Daniel@0:     s1 = .0207;
Daniel@0:     s2 = 18.96;
Daniel@0:     s = tril(xstar ./ (s1 * min(f1,f2) + s2 ));
Daniel@0:     pd = exp(-b1*s.*abs(f2-f1)) - exp(-b2*s.*abs(f2-f1));
Daniel@0: return