Daniel@0: function varargout = mircepstrum(orig,varargin)
Daniel@0: %   s = mircepstrum(x) computes the cepstrum, which indicates
Daniel@0: %       periodicities, and is used for instance for pitch detection.
Daniel@0: %   x can be either a spectrum, an audio signal, or the name of an audio file.
Daniel@0: %   Optional parameter:
Daniel@0: %       mircepstrum(...,'Min',min) specifies the lowest delay taken into
Daniel@0: %           consideration, in seconds.
Daniel@0: %           Default value: 0.0002 s (corresponding to a maximum frequency of 
Daniel@0: %               5 kHz).
Daniel@0: %       mircepstrum(...,'Max',max) specifies the highest delay taken into
Daniel@0: %           consideration, in seconds.
Daniel@0: %           Default value: 0.05 s (corresponding to a minimum frequency of 
Daniel@0: %               20 Hz).
Daniel@0: %       mircepstrum(...,'Freq') represents the cepstrum in the frequency 
Daniel@0: %           domain.
Daniel@0: 
Daniel@0:         mi.key = 'Min';
Daniel@0:         mi.type = 'Integer';
Daniel@0:         mi.default = 0.0002;
Daniel@0:         mi.unit = {'s','Hz'};
Daniel@0:         mi.defaultunit = 's';
Daniel@0:         mi.opposite = 'ma';
Daniel@0:     option.mi = mi;
Daniel@0: 
Daniel@0:         ma.key = 'Max';
Daniel@0:         ma.type = 'Integer';
Daniel@0:         ma.default = .05;
Daniel@0:         ma.unit = {'s','Hz'};
Daniel@0:         ma.defaultunit = 's';
Daniel@0:         ma.opposite = 'mi';
Daniel@0:     option.ma = ma;
Daniel@0:         
Daniel@0:         fr.key = 'Freq';
Daniel@0:         fr.type = 'Boolean';
Daniel@0:         fr.default = 0;
Daniel@0:     option.fr = fr;
Daniel@0: 
Daniel@0:         complex.key = 'Complex';
Daniel@0:         complex.type = 'Boolean';
Daniel@0:         complex.default = 0;
Daniel@0:     option.complex = complex;
Daniel@0: 
Daniel@0: specif.option = option;
Daniel@0: 
Daniel@0: specif.defaultframelength = 0.05;
Daniel@0: specif.defaultframehop = 0.5;
Daniel@0: 
Daniel@0: varargout = mirfunction(@mircepstrum,orig,varargin,nargout,specif,@init,@main);
Daniel@0: 
Daniel@0: 
Daniel@0: function [x type] = init(x,option)
Daniel@0: if not(isamir(x,'mircepstrum'))
Daniel@0:     x = mirspectrum(x);
Daniel@0: end
Daniel@0: type = 'mircepstrum';
Daniel@0: 
Daniel@0: 
Daniel@0: function c = main(orig,option,postoption)
Daniel@0: if iscell(orig)
Daniel@0:     orig = orig{1};
Daniel@0: end
Daniel@0: c.phase = [];
Daniel@0: if isa(orig,'mircepstrum')
Daniel@0:     c.freq = orig.freq;
Daniel@0: else
Daniel@0:     c.freq = 0;
Daniel@0: end
Daniel@0: c = class(c,'mircepstrum',mirdata(orig));
Daniel@0: c = purgedata(c);
Daniel@0: c = set(c,'Title','Cepstrum','Abs','quefrency (s)','Ord','magnitude');
Daniel@0: 
Daniel@0: if isa(orig,'mircepstrum')
Daniel@0:     if option.ma < Inf || option.mi > 0 || get(orig,'FreqDomain')
Daniel@0:         mag = get(orig,'Magnitude');
Daniel@0:         pha = get(orig,'Phase');
Daniel@0:         que = get(orig,'Quefrency');
Daniel@0:         for h = 1:length(mag)
Daniel@0:             for k = 1:length(mag{h})
Daniel@0:                 if get(orig,'FreqDomain')
Daniel@0:                     mag{h}{k} = flipud(mag{h}{k});
Daniel@0:                     que{h}{k} = flipud(1./que{h}{k});
Daniel@0:                     pha{h}{k} = flipud(pha{h}{k});
Daniel@0:                 end
Daniel@0:                 range = find(que{h}{k}(:,1,1) <= option.ma & ...
Daniel@0:                              que{h}{k}(:,1,1) >= option.mi);
Daniel@0:                 mag{h}{k} = mag{h}{k}(range,:,:);
Daniel@0:                 pha{h}{k} = pha{h}{k}(range,:,:);
Daniel@0:                 que{h}{k} = que{h}{k}(range,:,:);
Daniel@0:             end
Daniel@0:         end
Daniel@0:         c = set(c,'Magnitude',mag,'Phase',pha,'Quefrency',que,'FreqDomain',0);
Daniel@0:     end
Daniel@0:     c = modif(c,option);
Daniel@0: elseif isa(orig,'mirspectrum')
Daniel@0:     mag = get(orig,'Magnitude');
Daniel@0:     pha = get(orig,'Phase');
Daniel@0:     f = get(orig,'Sampling');
Daniel@0:     q = cell(1,length(mag));
Daniel@0:     for h = 1:length(mag)
Daniel@0:         len = ceil(option.ma*f{h});
Daniel@0:         start = ceil(option.mi*f{h})+1;
Daniel@0:         q{h} = cell(1,length(mag{h}));
Daniel@0:         for k = 1:length(mag{h})
Daniel@0:             m = mag{h}{k}.*exp(1i*pha{h}{k});
Daniel@0:             m = [m(1:end-1,:) ; conj(flipud(m))];  % Reconstitution of the complete abs(FFT)
Daniel@0:             if not(option.complex)
Daniel@0:                 m = abs(m);
Daniel@0:             end
Daniel@0:             m = log(m);
Daniel@0:             c0=fft(m);
Daniel@0:             q0=repmat((0:(size(c0,1)-1))'/f{k},[1,size(m,2),size(m,3)]);
Daniel@0:             len = min(len,floor(size(c0,1)/2));
Daniel@0:             mag{h}{k} = abs(c0(start:len,:,:));
Daniel@0:             if option.complex
Daniel@0:                 pha{h}{k} = unwrap(angle(c0(start:len,:,:)));    
Daniel@0:             else
Daniel@0:                 pha{h}{k} = nan(size(c0(start:len,:,:)));
Daniel@0:             end
Daniel@0:             q{h}{k} = q0(start:len,:,:);
Daniel@0:         end
Daniel@0:     end
Daniel@0:     c = set(c,'Magnitude',mag,'Phase',pha,'Quefrency',q);
Daniel@0:     c = modif(c,option);
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: function c = modif(c,option)
Daniel@0: mag = get(c,'Magnitude');
Daniel@0: que = get(c,'Quefrency');
Daniel@0: if option.fr && not(get(c,'FreqDomain'))
Daniel@0:     for k = 1:length(mag)
Daniel@0:         for l = 1:length(mag{k})
Daniel@0:             m = mag{k}{l};
Daniel@0:             q = que{k}{l};
Daniel@0:             if not(isempty(m))
Daniel@0:                 if q(1,1) == 0
Daniel@0:                     m = m(2:end,:,:);
Daniel@0:                     q = q(2:end,:,:);
Daniel@0:                 end
Daniel@0:                 m = flipud(m);
Daniel@0:                 q = flipud(1./q);
Daniel@0:             end
Daniel@0:             mag{k}{l} = m;
Daniel@0:             que{k}{l} = q;
Daniel@0:         end
Daniel@0:     end
Daniel@0:     c = set(c,'FreqDomain',1,'Abs','frequency (Hz)');
Daniel@0: end
Daniel@0: c = set(c,'Magnitude',mag,'Quefrency',que,'Freq');