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