wolffd@0: function varargout = mirchromagram(orig,varargin)
wolffd@0: %   c = mirchromagram(x) computes the chromagram, or distribution of energy 
wolffd@0: %       along pitches, of the audio signal x.
wolffd@0: %       (x can be the name of an audio file as well, or a spectrum, ...)
wolffd@0: %   Optional argument:
wolffd@0: %       c = mirchromagram(...,'Tuning',t): specifies the central frequency
wolffd@0: %           (in Hz.) associated to  chroma C.
wolffd@0: %               Default value, t = 261.6256 Hz
wolffd@0: %       c = mirchromagram(...,'Wrap',w): specifies whether the chromagram is
wolffd@0: %           wrapped or not.
wolffd@0: %           w = 1: groups all the pitches belonging to same pitch classes
wolffd@0: %               (default value)
wolffd@0: %           w = 0: pitches are considered as absolute values.
wolffd@0: %       c = mirchromagram(...,'Frame',l,h) orders a frame decomposition of window
wolffd@0: %           length l (in seconds) and hop factor h, expressed relatively to
wolffd@0: %           the window length. For instance h = 1 indicates no overlap.
wolffd@0: %           Default values: l = .2 seconds and h = .05
wolffd@0: %       c = mirchromagram(...,'Center'): centers the result.
wolffd@0: %       c = mirchromagram(...,'Normal',n): performs a n-norm of the
wolffd@0: %           resulting chromagram. Toggled off if n = 0
wolffd@0: %               Default value: n = Inf (corresponding to a normalization by
wolffd@0: %                   the maximum value).
wolffd@0: %       c = mirchromagram(...,'Pitch',p): specifies how to label chromas in
wolffd@0: %           the figures.
wolffd@0: %               p = 1: chromas are labeled using pitch names (default)
wolffd@0: %                   alternative syntax: chromagram(...,'Pitch')
wolffd@0: %               p = 0: chromas are labeled using MIDI pitch numbers
wolffd@0: %       c = mirchromagram(...,'Triangle'): weight the contribution of each
wolffd@0: %           frequency with respect to the distance with the actual
wolffd@0: %           frequency of the corresponding chroma.
wolffd@0: %       c = mirchromagram(...,'Weight',o): specifies the relative radius of
wolffd@0: %           the weighting window, with respect to the distance between
wolffd@0: %           frequencies of successive chromas.
wolffd@0: %           o = 1: each window begins at the centers of the previous one.
wolffd@0: %           o = .5: each window begins at the end of the previous one.
wolffd@0: %               (default value)
wolffd@0: %       mirchromagram(...,'Min',mi) indicates the lowest frequency taken into
wolffd@0: %           consideration in the spectrum computation, expressed in Hz.
wolffd@0: %           Default value: 100 Hz. (Gomez, 2006)
wolffd@0: %       mirchromagram(...,'Max',ma) indicates the highest frequency taken into
wolffd@0: %           consideration in the spectrum computation, expressed in Hz.
wolffd@0: %           This upper limit is further shifted to a highest value until
wolffd@0: %           the frequency range covers an exact multiple of octaves.
wolffd@0: %           Default value: 5000 Hz. (Gomez, 2006)
wolffd@0: %       mirchromagram(...,'Res',r) indicates the resolution of the
wolffd@0: %           chromagram in number of bins per octave.
wolffd@0: %               Default value, r = 12.
wolffd@0: %
wolffd@0: % Gómez, E. (2006). Tonal description of music audio signal. Phd thesis, 
wolffd@0: %   Universitat Pompeu Fabra, Barcelona .
wolffd@0: 
wolffd@0:         cen.key = 'Center';
wolffd@0:         cen.type = 'Boolean';
wolffd@0:         cen.default = 0;
wolffd@0:     option.cen = cen;
wolffd@0:     
wolffd@0:         nor.key = {'Normal','Norm'};
wolffd@0:         nor.type = 'Integer';
wolffd@0:         nor.default = Inf;
wolffd@0:     option.nor = nor;
wolffd@0:     
wolffd@0:         wth.key = 'Weight';
wolffd@0:         wth.type = 'Integer';
wolffd@0:         wth.default = .5;
wolffd@0:     option.wth = wth;
wolffd@0:     
wolffd@0:         tri.key = 'Triangle';
wolffd@0:         tri.type = 'Boolean';
wolffd@0:         tri.default = 0;
wolffd@0:     option.tri = tri;
wolffd@0:     
wolffd@0:         wrp.key = 'Wrap';
wolffd@0:         wrp.type = 'Boolean';
wolffd@0:         wrp.default = 1;
wolffd@0:     option.wrp = wrp;
wolffd@0:     
wolffd@0:         plabel.key = 'Pitch';
wolffd@0:         plabel.type = 'Boolean';
wolffd@0:         plabel.default = 1;
wolffd@0:     option.plabel = plabel;
wolffd@0:     
wolffd@0:         thr.key = {'Threshold','dB'};
wolffd@0:         thr.type = 'Integer';
wolffd@0:         thr.default = 20;
wolffd@0:     option.thr = thr;
wolffd@0:     
wolffd@0:         min.key = 'Min';
wolffd@0:         min.type = 'Integer';
wolffd@0:         min.default = 100;
wolffd@0:     option.min = min;
wolffd@0:     
wolffd@0:         max.key = 'Max';
wolffd@0:         max.type = 'Integer';
wolffd@0:         max.default = 5000;
wolffd@0:     option.max = max;
wolffd@0: 
wolffd@0:         res.key = 'Res';
wolffd@0:         res.type = 'Integer';
wolffd@0:         res.default = 12;
wolffd@0:     option.res = res;
wolffd@0: 
wolffd@0:         origin.key = 'Tuning';
wolffd@0:         origin.type = 'Integer';
wolffd@0:         origin.default = 261.6256;
wolffd@0:     option.origin = origin;
wolffd@0:    
wolffd@0: specif.option = option;
wolffd@0: specif.defaultframelength = .2;
wolffd@0: specif.defaultframehop = .05;
wolffd@0: 
wolffd@0: varargout = mirfunction(@mirchromagram,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: if isamir(x,'mirtemporal') || isamir(x,'mirspectrum')
wolffd@0:     freqmin = option.min;
wolffd@0:     freqmax = freqmin*2;
wolffd@0:     while freqmax < option.max
wolffd@0:         freqmax = freqmax*2;
wolffd@0:     end
wolffd@0:     %freqres = freqmin*(2.^(1/option.res)-1);
wolffd@0:         % Minimal frequency resolution should correspond to frequency range
wolffd@0:         %   between the first two bins of the chromagram 
wolffd@0:         
wolffd@0:     x = mirspectrum(x,'dB',option.thr,'Min',freqmin,'Max',freqmax,...
wolffd@0:                       'NormalInput','MinRes',option.res,'OctaveRatio',.85);
wolffd@0:                   %freqres*.5,...
wolffd@0:                   %    'WarningRes',freqres);
wolffd@0: end
wolffd@0: type = 'mirchromagram';
wolffd@0: 
wolffd@0: 
wolffd@0: function c = main(orig,option,postoption)
wolffd@0: if iscell(orig)
wolffd@0:     orig = orig{1};
wolffd@0: end
wolffd@0: if option.res == 12
wolffd@0:     chromascale = {'C','C#','D','D#','E','F','F#','G','G#','A','A#','B'};
wolffd@0: else
wolffd@0:     chromascale = 1:option.res;
wolffd@0:     option.plabel = 0;
wolffd@0: end
wolffd@0: if isa(orig,'mirchromagram')
wolffd@0:     c = modif(orig,option,chromascale);
wolffd@0: else
wolffd@0:     c.plabel = 1;
wolffd@0:     c.wrap = 0;
wolffd@0:     c.chromaclass = {};
wolffd@0:     c.chromafreq = {};
wolffd@0:     c.register = {};
wolffd@0:     c = class(c,'mirchromagram',mirdata(orig));
wolffd@0:     c = purgedata(c);
wolffd@0:     c = set(c,'Title','Chromagram','Ord','magnitude','Interpolable',0);
wolffd@0:     if option.wrp
wolffd@0:         c = set(c,'Abs','chroma class');
wolffd@0:     else
wolffd@0:         c = set(c,'Abs','chroma');
wolffd@0:     end
wolffd@0:     m = get(orig,'Magnitude');
wolffd@0:     f = get(orig,'Frequency');
wolffd@0:     %disp('Computing chromagram...')
wolffd@0:     fs = get(orig,'Sampling');
wolffd@0:     n = cell(1,length(m));  % The final structured list of magnitudes.
wolffd@0:     cc = cell(1,length(m));  % The final structured list of chroma classes.
wolffd@0:     o = cell(1,length(m));  % The final structured list of octave registers.
wolffd@0:     p = cell(1,length(m));  % The final structured list of chromas.
wolffd@0:     cf = cell(1,length(m));  % The final structured list of central frequencies related to chromas.
wolffd@0:     for i = 1:length(m)
wolffd@0:         mi = m{i};
wolffd@0:         fi = f{i};
wolffd@0:         if not(iscell(mi))
wolffd@0:             mi = {mi};
wolffd@0:             fi = {fi};
wolffd@0:         end
wolffd@0:         ni = cell(1,length(mi));    % The list of magnitudes.
wolffd@0:         ci = cell(1,length(mi));    % The list of chroma classes.
wolffd@0:         oi = cell(1,length(mi));    % The list of octave registers.
wolffd@0:         pi = cell(1,length(mi));    % The list of absolute chromas.
wolffd@0:         cfi = cell(1,length(mi));    % The central frequency of each chroma.
wolffd@0:         for j = 1:length(mi)
wolffd@0:             mj = mi{j};
wolffd@0:             fj = fi{j};
wolffd@0:                         
wolffd@0:             % Let's remove the frequencies exceeding the last whole octave.
wolffd@0:             minfj = min(min(min(fj)));
wolffd@0:             maxfj = max(max(max(fj)));
wolffd@0:             maxfj = minfj*2^(floor(log2(maxfj/minfj)));
wolffd@0:             fz = find(fj(:,1,1,1) > maxfj);
wolffd@0:             mj(fz,:,:,:) = [];      
wolffd@0:             fj(fz,:,:,:) = [];
wolffd@0:             
wolffd@0:             [s1 s2 s3] = size(mj);
wolffd@0:             
wolffd@0:             cj = freq2chro(fj,option.res,option.origin);
wolffd@0:             if not(ismember(min(cj)+1,cj))
wolffd@0:                 warning('WARNING IN MIRCHROMAGRAM: Frequency resolution of the spectrum is too low.');    
wolffd@0:                 display('The conversion of low frequencies into chromas may be incorrect.');
wolffd@0:             end
wolffd@0:             ccj = min(min(min(cj))):max(max(max(cj)));
wolffd@0:             sc = length(ccj);   % The size of range of absolute chromas.
wolffd@0:             mat = zeros(s1,sc);
wolffd@0:             fc = chro2freq(ccj,option.res,option.origin);   % The absolute chromas in Hz.
wolffd@0:             fl = chro2freq(ccj-1,option.res,option.origin); % Each previous chromas in Hz.
wolffd@0:             fr = chro2freq(ccj+1,option.res,option.origin); % Each related next chromas in Hz.
wolffd@0:             for k = 1:sc
wolffd@0:                 rad = find(and(fj(:,1) > fc(k)-option.wth*(fc(k)-fl(k)),...
wolffd@0:                                fj(:,1) < fc(k)-option.wth*(fc(k)-fr(k))));
wolffd@0:                 if option.tri
wolffd@0:                     dist = fc(k) - fj(:,1,1,1);
wolffd@0:                     rad1 = dist/(fc(k) - fl(k))/option.wth;
wolffd@0:                     rad2 = dist/(fc(k) - fr(k))/option.wth;
wolffd@0:                     ndist = max(rad1,rad2);
wolffd@0:                     mat(:,k) = max(min(1-ndist,1),0)/length(rad);
wolffd@0:                 else
wolffd@0:                     mat(rad,k) = ones(length(rad),1)/length(rad);
wolffd@0:                 end
wolffd@0:                 if k ==1 || k == sc
wolffd@0:                     mat(:,k) = mat(:,k)/2;
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             nj = zeros(sc,s2,s3);
wolffd@0:             for k = 1:s2
wolffd@0:                 for l = 1:s3
wolffd@0:                     nj(:,k,l) = (mj(:,k,l)'*mat)';
wolffd@0:                 end
wolffd@0:             end
wolffd@0:             cj = mod(ccj',option.res);
wolffd@0:             oi{j} = floor(ccj/option.res)+4;
wolffd@0:             if option.plabel
wolffd@0:                 pj = strcat(chromascale(cj+1)',num2str(oi{j}'));
wolffd@0:             else
wolffd@0:                 pj = ccj'+60;
wolffd@0:             end
wolffd@0:             ci{j} = repmat(cj,[1,s2,s3]);
wolffd@0:             pi{j} = repmat(pj,[1,s2,s3]);
wolffd@0:             ni{j} = nj;
wolffd@0:             cfi{j} = fc;
wolffd@0:         end
wolffd@0:         n{i} = ni;
wolffd@0:         cc{i} = ci;
wolffd@0:         o{i} = oi;
wolffd@0:         p{i} = pi;
wolffd@0:         cf{i} = cfi;
wolffd@0:     end
wolffd@0:     c = set(c,'Magnitude',n,'Chroma',p,'ChromaClass',cc,...
wolffd@0:               'ChromaFreq',cf,'Register',o);
wolffd@0:     c = modif(c,option,chromascale);
wolffd@0:     c = {c orig};
wolffd@0: end
wolffd@0: 
wolffd@0:    
wolffd@0: function c = modif(c,option,chromascale)
wolffd@0: if option.plabel
wolffd@0:     c = set(c,'PitchLabel',1);
wolffd@0: end            
wolffd@0: if option.cen || option.nor || option.wrp
wolffd@0:     n = get(c,'Magnitude');
wolffd@0:     p = get(c,'Chroma');
wolffd@0:     cl = get(c,'ChromaClass');
wolffd@0:     fp = get(c,'FramePos');
wolffd@0:     n2 = cell(1,length(n));
wolffd@0:     p2 = cell(1,length(n));
wolffd@0:     wrp = option.wrp && not(get(c,'Wrap'));
wolffd@0:     for i = 1:length(n)
wolffd@0:         ni = n{i};
wolffd@0:         pi = p{i};
wolffd@0:         cli = cl{i};
wolffd@0:         if not(iscell(ni))
wolffd@0:             ni = {ni};
wolffd@0:             pi = {pi};
wolffd@0:             cli = {cli};
wolffd@0:         end
wolffd@0:         if wrp
wolffd@0:             c = set(c,'Wrap',option.wrp);
wolffd@0:         end
wolffd@0:         n2i = cell(1,length(ni));
wolffd@0:         p2i = cell(1,length(ni));
wolffd@0:         for j = 1:length(ni)
wolffd@0:             nj = ni{j};
wolffd@0:             pj = pi{j};
wolffd@0:             clj = cli{j};
wolffd@0:             if wrp
wolffd@0:                 n2j = zeros(option.res,size(nj,2),size(nj,3));
wolffd@0:                 for k = 1:size(pj,1)
wolffd@0:                     n2j(clj(k)+1,:,:) = n2j(clj(k)+1,:,:) + nj(k,:,:);  % squared sum (parameter)
wolffd@0:                 end
wolffd@0:                 p2i{j} = chromascale';
wolffd@0:             else
wolffd@0:                 n2j = nj;
wolffd@0:                 p2i{j} = pi{j};
wolffd@0:             end
wolffd@0:             if option.cen
wolffd@0:                 n2j = n2j - repmat(mean(n2j),[size(n2j,1),1,1]);
wolffd@0:             end
wolffd@0:             if option.nor
wolffd@0:                 n2j = n2j ./ repmat(vectnorm(n2j,option.nor) + ...
wolffd@0:                     repmat(1e-6,[1,size(n2j,2),size(n2j,3)] )...
wolffd@0:                     ,[size(n2j,1),1,1]);
wolffd@0:             end
wolffd@0:             n2i{j} = n2j;
wolffd@0:         end
wolffd@0:         n2{i} = n2i;
wolffd@0:         p2{i} = p2i;
wolffd@0:     end
wolffd@0:     c = set(c,'Magnitude',n2,'Chroma',p2,'FramePos',fp);
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function c = freq2chro(f,res,origin)
wolffd@0: c = round(res*log2(f/origin));
wolffd@0: 
wolffd@0: 
wolffd@0: function f = chro2freq(c,res,origin)
wolffd@0: f = 2.^(c/res)*origin;
wolffd@0: 
wolffd@0: 
wolffd@0: function y = vectnorm(x,p)
wolffd@0: if isinf(p)
wolffd@0:     y = max(x);
wolffd@0: else
wolffd@0:     y = sum(abs(x).^p).^(1/p);
wolffd@0: end