wolffd@0
|
1 function varargout = mirroughness(x,varargin)
|
wolffd@0
|
2 % r = mirroughness(x) calculates the roughness, or sensory dissonance,
|
wolffd@0
|
3 % due to beating phenomenon between close frequency peaks.
|
wolffd@0
|
4 % The frequency components are supposed to remain sufficiently
|
wolffd@0
|
5 % constant throughout each frame of each audio file.
|
wolffd@0
|
6 % r = mirroughness(...,'Contrast',c) specifies the contrast parameter
|
wolffd@0
|
7 % used for peak picking (cf. mirpeaks).
|
wolffd@0
|
8 % Default value: c = .01
|
wolffd@0
|
9 % [r,s] = mirroughness(x) also displays the spectrum and its peaks, used
|
wolffd@0
|
10 % for the computation of roughness.
|
wolffd@0
|
11 % Optional arguments:
|
wolffd@0
|
12 % Method used:
|
wolffd@0
|
13 % mirroughness(...,'Sethares') (default): based on the summation
|
wolffd@0
|
14 % of roughness between all pairs of sines (obtained through
|
wolffd@0
|
15 % spectral peak-picking).
|
wolffd@0
|
16 % mirroughness(...,'Vassilakis'): variant of 'Sethares' model
|
wolffd@0
|
17 % with a more complex weighting (Vassilakis, 2001, Eq. 6.23).
|
wolffd@0
|
18
|
wolffd@0
|
19 meth.type = 'String';
|
wolffd@0
|
20 meth.choice = {'Sethares','Vassilakis'};
|
wolffd@0
|
21 meth.default = 'Sethares';
|
wolffd@0
|
22 option.meth = meth;
|
wolffd@0
|
23
|
wolffd@0
|
24 cthr.key = 'Contrast';
|
wolffd@0
|
25 cthr.type = 'Integer';
|
wolffd@0
|
26 cthr.default = .01;
|
wolffd@0
|
27 option.cthr = cthr;
|
wolffd@0
|
28
|
wolffd@0
|
29 frame.key = 'Frame';
|
wolffd@0
|
30 frame.type = 'Integer';
|
wolffd@0
|
31 frame.number = 2;
|
wolffd@0
|
32 frame.default = [.05 .5];
|
wolffd@0
|
33 option.frame = frame;
|
wolffd@0
|
34
|
wolffd@0
|
35 specif.option = option;
|
wolffd@0
|
36 specif.defaultframelength = .05;
|
wolffd@0
|
37 specif.defaultframehop = .5;
|
wolffd@0
|
38
|
wolffd@0
|
39
|
wolffd@0
|
40 varargout = mirfunction(@mirroughness,x,varargin,nargout,specif,@init,@main);
|
wolffd@0
|
41
|
wolffd@0
|
42
|
wolffd@0
|
43 function [x type] = init(x,option)
|
wolffd@0
|
44 if isamir(x,'miraudio') && not(isframed(x))
|
wolffd@0
|
45 x = mirframenow(x,option);
|
wolffd@0
|
46 end
|
wolffd@0
|
47 x = mirspectrum(x);
|
wolffd@0
|
48 if not(haspeaks(x))
|
wolffd@0
|
49 x = mirpeaks(x,'Contrast',option.cthr);
|
wolffd@0
|
50 end
|
wolffd@0
|
51 type = 'mirscalar';
|
wolffd@0
|
52
|
wolffd@0
|
53
|
wolffd@0
|
54 function r = main(p,option,postoption)
|
wolffd@0
|
55 if iscell(p)
|
wolffd@0
|
56 p = p{1};
|
wolffd@0
|
57 end
|
wolffd@0
|
58 if strcmpi(option.meth,'Sethares') || strcmpi(option.meth,'Vassilakis')
|
wolffd@0
|
59 pf = get(p,'PeakPosUnit');
|
wolffd@0
|
60 pv = get(p,'PeakVal');
|
wolffd@0
|
61 rg = cell(1,length(pf));
|
wolffd@0
|
62 for h = 1:length(pf)
|
wolffd@0
|
63 rg{h} = cell(1,length(pf{h}));
|
wolffd@0
|
64 for i = 1:length(pf{h})
|
wolffd@0
|
65 pfi = pf{h}{i};
|
wolffd@0
|
66 pvi = pv{h}{i};
|
wolffd@0
|
67 rg{h}{i} = zeros(1,length(pfi));
|
wolffd@0
|
68 for k = 1:size(pfi,3)
|
wolffd@0
|
69 for j = 1:size(pfi,2)
|
wolffd@0
|
70 pfj = pfi{1,j,k};
|
wolffd@0
|
71 pvj = pvi{1,j,k};
|
wolffd@0
|
72 f1 = repmat(pfj,[1 length(pfj)]);
|
wolffd@0
|
73 f2 = repmat(pfj',[length(pfj) 1]);
|
wolffd@0
|
74 v1 = repmat(pvj,[1 length(pvj)]);
|
wolffd@0
|
75 v2 = repmat(pvj',[length(pvj) 1]);
|
wolffd@0
|
76 rj = plomp(f1,f2);
|
wolffd@0
|
77 if strcmpi(option.meth,'Sethares')
|
wolffd@0
|
78 rj = v1.*v2.*rj;
|
wolffd@0
|
79 elseif strcmpi(option.meth,'Vassilakis')
|
wolffd@0
|
80 rj = (v1.*v2).^.1.*.5.*(2*v2./(v1+v2)).^3.11.*rj;
|
wolffd@0
|
81 end
|
wolffd@0
|
82 rg{h}{i}(1,j,k) = sum(sum(rj));
|
wolffd@0
|
83 end
|
wolffd@0
|
84 end
|
wolffd@0
|
85 end
|
wolffd@0
|
86 end
|
wolffd@0
|
87 else
|
wolffd@0
|
88 end
|
wolffd@0
|
89 r = mirscalar(p,'Data',rg,'Title','Roughness');
|
wolffd@0
|
90 r = {r,p};
|
wolffd@0
|
91
|
wolffd@0
|
92
|
wolffd@0
|
93 function pd = plomp(f1, f2)
|
wolffd@0
|
94 % returns the dissonance of two pure tones at frequencies f1 & f2 Hz
|
wolffd@0
|
95 % according to the Plomp-Levelt curve (see Sethares)
|
wolffd@0
|
96 b1 = 3.51;
|
wolffd@0
|
97 b2 = 5.75;
|
wolffd@0
|
98 xstar = .24;
|
wolffd@0
|
99 s1 = .0207;
|
wolffd@0
|
100 s2 = 18.96;
|
wolffd@0
|
101 s = tril(xstar ./ (s1 * min(f1,f2) + s2 ));
|
wolffd@0
|
102 pd = exp(-b1*s.*abs(f2-f1)) - exp(-b2*s.*abs(f2-f1));
|
wolffd@0
|
103 return |