wolffd@0: function varargout = mirmfcc(orig,varargin)
wolffd@0: %   c = mirmfcc(a) finds the Mel frequency cepstral coefficients (ceps),
wolffd@0: %       a numerical description of the spectrum envelope.
wolffd@0: %
wolffd@0: %   Requires the Auditory Toolbox.
wolffd@0: %
wolffd@0: %   Optional arguments:
wolffd@0: %       c = mirmfcc(...,'Rank',N) computes the coefficients of rank(s) N
wolffd@0: %           (default: N = 1:13).
wolffd@0: %   If a is a frame decomposition, the temporal evolution of the MFCC,
wolffd@0: %       along the successive frames, is returned. In this case, a second
wolffd@0: %       option is available:
wolffd@0: %       mirmfcc(...,'Delta',d) performs d temporal differentiations of
wolffd@0: %           the coefficients, also called delta-MFCC (for d = 1) or
wolffd@0: %           delta-delta-MFCC (for d = 2).
wolffd@0: %       mirmfcc(...,'Delta') corresponds to mirmfcc(...,'Delta',1)
wolffd@0: %   Optional arguments related to the delta computation:
wolffd@0: %       mirmfcc(...,'Radius',r) specifies, for each frame, the number of
wolffd@0: %           successive and previous neighbouring frames taken into
wolffd@0: %           consideration for the least-square approximation.
wolffd@0: %           Usually 1 or 2.
wolffd@0: %           Default value: 2.
wolffd@0: 
wolffd@0:         nbbands.key = 'Bands';
wolffd@0:         nbbands.type = 'Integer';
wolffd@0:         nbbands.default = 40;
wolffd@0:     option.nbbands = nbbands;
wolffd@0: 
wolffd@0:         rank.key = 'Rank';
wolffd@0:         rank.type = 'Integer';
wolffd@0:         rank.default = 1:13;
wolffd@0:     option.rank = rank; 
wolffd@0:         
wolffd@0:         delta.key = 'Delta';
wolffd@0:         delta.type = 'Integer';
wolffd@0:         delta.default = 0;
wolffd@0:         delta.keydefault = 1;
wolffd@0:     option.delta = delta;
wolffd@0:         
wolffd@0:         radius.key = 'Radius';
wolffd@0:         radius.type = 'Integer';
wolffd@0:         radius.default = 2;
wolffd@0:     option.radius = radius;
wolffd@0:         
wolffd@0: specif.option = option;
wolffd@0:      
wolffd@0: varargout = mirfunction(@mirmfcc,orig,varargin,nargout,specif,@init,@main);
wolffd@0: 
wolffd@0: 
wolffd@0: function [x type] = init(x,option)
wolffd@0: if isamir(x,'miraudio') || isamir(x,'mirspectrum')
wolffd@0:     x = mirspectrum(x,'Mel','log','Bands',option.nbbands);
wolffd@0: end
wolffd@0: type = 'mirmfcc';
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 isa(orig,'mirmfcc')
wolffd@0:     c = orig;
wolffd@0:     if option.rank
wolffd@0:         magn = get(c,'Data');
wolffd@0:         rank = get(c,'Rank');
wolffd@0:         for h = 1:length(magn)
wolffd@0:             for k = 1:length(magn{h})
wolffd@0:                 m = magn{h}{k};
wolffd@0:                 r = rank{h}{k};
wolffd@0:                 r1 = r(:,1,1);
wolffd@0:                 range = find(ismember(r1,option.rank));
wolffd@0:                 magn{h}{k} = m(range,:,:);
wolffd@0:                 rank{h}{k} = r(range,:,:);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         c = set(c,'Data',magn,'Rank',rank);
wolffd@0:     end
wolffd@0:     c = modif(c,option);
wolffd@0: else
wolffd@0:     c.delta = 0;
wolffd@0:     %disp('Computing Mel frequency cepstral coefficients...');
wolffd@0:     e = get(orig,'Magnitude');
wolffd@0: 
wolffd@0:     % The following is largely based on the source code from Auditory Toolbox 
wolffd@0:     % (A part that I could not call directly from MIRtoolbox)
wolffd@0:     
wolffd@0:     % (Malcolm Slaney, August 1993, (c) 1998 Interval Research Corporation)
wolffd@0:     
wolffd@0:     try
wolffd@0:         MakeERBFilters(1,1,1); % Just to be sure that the Auditory Toolbox is installed
wolffd@0:     catch
wolffd@0:         error(['ERROR IN MIRFILTERBANK: Auditory Toolbox needs to be installed.']);
wolffd@0:     end  
wolffd@0:     
wolffd@0:     dc = cell(1,length(e));
wolffd@0:     rk = cell(1,length(e));
wolffd@0:     for h = 1:length(e)
wolffd@0:         dc{h} = cell(1,length(e{h}));
wolffd@0:         rk{h} = cell(1,length(e{h}));
wolffd@0:         for i = 1:length(e{h})
wolffd@0:             ei = e{h}{i};
wolffd@0:             totalFilters = size(ei,3); %Number of mel bands.
wolffd@0: 
wolffd@0:             % Figure out Discrete Cosine Transform.  We want a matrix
wolffd@0:             % dct(i,j) which is totalFilters x cepstralCoefficients in size.
wolffd@0:             % The i,j component is given by 
wolffd@0:             %                cos( i * (j+0.5)/totalFilters pi )
wolffd@0:             % where we have assumed that i and j start at 0.
wolffd@0:             mfccDCTMatrix = 1/sqrt(totalFilters/2)*...
wolffd@0:                             cos(option.rank' * ...
wolffd@0:                                 (2*(0:(totalFilters-1))+1) * ...
wolffd@0:                                  pi/2/totalFilters);
wolffd@0:             rank0 = find(option.rank == 0);
wolffd@0:             mfccDCTMatrix(rank0,:) = mfccDCTMatrix(rank0,:) * sqrt(2)/2;
wolffd@0:             ceps = zeros(size(mfccDCTMatrix,1),size(ei,2));
wolffd@0:             for j = 1:size(ei,2)
wolffd@0:                 ceps(:,j) = mfccDCTMatrix * permute(ei(1,j,:),[3 1 2]);
wolffd@0:             end
wolffd@0:             dc{h}{i} = ceps;
wolffd@0:             rk{h}{i} = repmat(option.rank(:),[1 size(ceps,2) size(ceps,3)]);
wolffd@0:         end
wolffd@0:     end
wolffd@0:     c = class(c,'mirmfcc',mirdata(orig));
wolffd@0:     c = purgedata(c);
wolffd@0:     c = set(c,'Title','MFCC','Abs','coefficient ranks','Ord','magnitude',...
wolffd@0:               'Data',dc,'Rank',rk);
wolffd@0:     c = modif(c,option);
wolffd@0: end
wolffd@0: c = {c orig};
wolffd@0: 
wolffd@0: 
wolffd@0: function c = modif(c,option)
wolffd@0: d = get(c,'Data');
wolffd@0: fp = get(c,'FramePos');
wolffd@0: t = get(c,'Title');
wolffd@0: if option.delta
wolffd@0:     M = option.radius;
wolffd@0:     for k = 1:option.delta
wolffd@0:         for h = 1:length(d)
wolffd@0:             for i = 1:length(d{h})
wolffd@0:                 nc = size(d{h}{i},2)-2*M;
wolffd@0:                 di = zeros(size(d{h}{i},1),nc);
wolffd@0:                 for j = 1:M
wolffd@0:                     di = di + j * (d{h}{i}(:,M+j+(1:nc)) ...
wolffd@0:                                  - d{h}{i}(:,M-j+(1:nc)));
wolffd@0:                 end
wolffd@0:                 di = di / 2 / sum((1:M).^2); % MULTIPLY BY 2 INSTEAD OF SQUARE FOR NORMALIZATION ?
wolffd@0:                 d{h}{i} = di;
wolffd@0:                 fp{h}{i} = fp{h}{i}(:,M+1:end-M);
wolffd@0:             end
wolffd@0:         end
wolffd@0:         t = ['Delta-',t];
wolffd@0:     end
wolffd@0: end
wolffd@0: c = set(c,'Data',d,'FramePos',fp,'Delta',get(c,'Delta')+option.delta,...
wolffd@0:           'Title',t);