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