Mercurial > hg > camir-aes2014
diff toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirautocor/mirautocor.m @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirautocor/mirautocor.m Tue Feb 10 15:05:51 2015 +0000 @@ -0,0 +1,604 @@ +function varargout = mirautocor(orig,varargin) +% a = mirautocor(x) computes the autocorrelation function related to x. +% Optional parameters: +% mirautocor(...,'Min',mi) indicates the lowest delay taken into +% consideration. The unit can be precised: +% mirautocor(...,'Min',mi,'s') (default unit) +% mirautocor(...,'Min',mi,'Hz') +% Default value: 0 s. +% mirautocor(...,'Max',ma) indicates the highest delay taken into +% consideration. The unit can be specified as for 'Min'. +% Default value: +% if x is a signal, the highest delay is 0.05 s +% (corresponding to a minimum frequency of 20 Hz). +% if x is an envelope, the highest delay is 2 s. +% mirautocor(...,'Resonance',r) multiplies the autocorrelation function +% with a resonance curve: +% Possible values: +% 'Toiviainen' from (Toiviainen & Snyder, 2003) +% 'vanNoorden' from (van Noorden & Moelants, 2001) +% mirautocor(...,'Center',c) assigns the center value of the +% resonance curve, in seconds. +% Works mainly with 'Toiviainen' option. +% Default value: c = 0.5 +% mirautocor(...,'Enhanced',a) reduces the effect of subharmonics. +% The original autocorrelation function is half-wave rectified, +% time-scaled by the factor a (which can be a factor list as +% well), and substracted from the original clipped function. +% (Tolonen & Karjalainen) +% If the 'Enhanced' option is not followed by any value, +% default value is a = 2:10 +% mirautocor(...,'Halfwave') performs a half-wave rectification on the +% result. +% mirautocor(...,'Freq') represents the autocorrelation function in the +% frequency domain. +% mirautocor(...,'NormalWindow',w): applies a window to the input +% signal and divides the autocorrelation by the autocorrelation of +% that window (Boersma 1993). +% Possible values: any windowing function proposed in the Signal +% Processing Toolbox (help window) plus 'rectangle' (no +% windowing) +% Default value: w = 'hanning' +% mirautocor(...,'NormalWindow',0): toggles off this normalization +% (which is on by default). +% All the parameters described previously can be applied to an +% autocorrelation function itself, in order to arrange the results +% after the actual computation of the autocorrelation computations. +% For instance: a = mirautocor(a,'Resonance','Enhanced') +% Other optional parameter: +% mirautocor(...,'Compres',k) computes the autocorrelation in the +% frequency domain and includes a magnitude compression of the +% spectral representation. A normal autocorrelation corresponds +% to the value k=2, but values lower than 2 are suggested by +% (Tolonen & Karjalainen, 2000). +% Default value: k = 0.67 +% mirautocor(...,'Normal',n) or simply mirautocor(...,n) specifies +% the normalization strategy. Accepted values are 'biased', +% 'unbiased', 'coeff' (default value) and 'none'. +% See help xcorr for an explanation. + + min.key = 'Min'; + min.type = 'Integer'; + min.unit = {'s','Hz'}; + if isamir(orig,'mirspectrum') + min.defaultunit = 'Hz'; + else + min.defaultunit = 's'; + end + min.default = 0; + min.opposite = 'max'; + option.min = min; + + max.key = 'Max'; + max.type = 'Integer'; + max.unit = {'s','Hz'}; + if isamir(orig,'mirspectrum') + max.defaultunit = 'Hz'; + else + max.defaultunit = 's'; + end + if isamir(orig,'mirenvelope') || isamir(orig,'mirdiffenvelope') + max.default = 2; % for envelopes, longest period: 2 seconds. + elseif isamir(orig,'miraudio') || ischar(orig) % for audio signal,lowest frequency: 20 Hz. + max.default = 1/20; + else + max.default = Inf; + end + max.opposite = 'min'; + option.max = max; + + scaleoptbw.key = 'Normal'; %'Normal' keyword optional + scaleoptbw.key = 'Boolean'; + option.scaleoptbw = scaleoptbw; + scaleopt.type = 'String'; + scaleopt.choice = {'biased','unbiased','coeff','none'}; + scaleopt.default = 'coeff'; + option.scaleopt = scaleopt; + + gener.key = {'Generalized','Compres'}; + gener.type = 'Integer'; + gener.default = 2; + gener.keydefault = .67; + option.gener = gener; + + ni.key = 'NormalInput'; %% Normalize before frame or chunk?? + ni.type = 'Boolean'; + ni.default = 0; + option.ni = ni; + + reso.key = 'Resonance'; + reso.type = 'String'; + reso.choice = {'ToiviainenSnyder','Toiviainen','vanNoorden','no','off',0}; + reso.keydefault = 'Toiviainen'; + reso.when = 'After'; + reso.default = 0; + option.reso = reso; + + resocenter.key = {'Center','Centre'}; + resocenter.type = 'Integer'; + resocenter.when = 'After'; + option.resocenter = resocenter; + + h.key = 'Halfwave'; + h.type = 'Boolean'; + h.when = 'After'; + h.default = 0; + option.h = h; + + e.key = 'Enhanced'; + e.type = 'Integers'; + e.default = []; + e.keydefault = 2:10; + e.when = 'After'; + option.e = e; + + fr.key = 'Freq'; + fr.type = 'Boolean'; + fr.default = 0; + fr.when = 'After'; + option.fr = fr; + + nw.key = 'NormalWindow'; + nw.when = 'Both'; + if isamir(orig,'mirspectrum') + nw.default = 0; + elseif isamir(orig,'mirenvelope') + nw.default = 'rectangular'; + else + nw.default = 'hanning'; + end + option.nw = nw; + + win.key = 'Window'; + win.type = 'String'; + win.default = NaN; + option.win = win; + +specif.option = option; + +specif.defaultframelength = 0.05; +specif.defaultframehop = 0.5; +specif.eachchunk = @eachchunk; +specif.combinechunk = @combinechunk; + +if isamir(orig,'mirscalar') || isamir(orig,'mirenvelope') + specif.nochunk = 1; +end + +varargout = mirfunction(@mirautocor,orig,varargin,nargout,specif,@init,@main); + + +function [x type] = init(x,option) +type = 'mirautocor'; + + +function a = main(orig,option,postoption) +if iscell(orig) + orig = orig{1}; +end +if isa(orig,'mirautocor') + a = orig; + if not(isempty(option)) && ... + (option.min || iscell(option.max) || option.max < Inf) + coeff = get(a,'Coeff'); + delay = get(a,'Delay'); + for h = 1:length(coeff) + if a.freq + mi = 1/option.max; + ma = 1/option.min; + else + mi = option.min; + ma = option.max; + end + for k = 1:length(coeff{h}) + range = find(and(delay{h}{k}(:,1,1) >= mi,... + delay{h}{k}(:,1,1) <= ma)); + coeff{h}{k} = coeff{h}{k}(range,:,:); + delay{h}{k} = delay{h}{k}(range,:,:); + end + end + a = set(a,'Coeff',coeff,'Delay',delay); + end + if not(isempty(postoption)) && not(isequal(postoption,0)) + a = post(a,postoption); + end +elseif ischar(orig) + a = mirautocor(miraudio(orig),option,postoption); +else + if nargin == 0 + orig = []; + end + a.freq = 0; + a.ofspectrum = 0; + a.window = {}; + a.normalwindow = 0; + a = class(a,'mirautocor',mirdata(orig)); + a = purgedata(a); + a = set(a,'Ord','coefficients'); + sig = get(orig,'Data'); + if isa(orig,'mirspectrum') + a = set(a,'Title','Spectrum autocorrelation','OfSpectrum',1,... + 'Abs','frequency (Hz)'); + pos = get(orig,'Pos'); + else + if isa(orig,'mirscalar') + a = set(a,'Title',[get(orig,'Title') ' autocorrelation']); + pos = get(orig,'FramePos'); + for k = 1:length(sig) + for l = 1:length(sig{k}) + sig{k}{l} = sig{k}{l}'; + pos{k}{l} = pos{k}{l}(1,:,:)'; + end + end + else + if isa(orig,'mirenvelope') + a = set(a,'Title','Envelope autocorrelation'); + elseif not(isa(orig,'mirautocor')) + a = set(a,'Title','Waveform autocorrelation'); + end + pos = get(orig,'Pos'); + end + a = set(a,'Abs','lag (s)'); + end + f = get(orig,'Sampling'); + + if isnan(option.win) + if isequal(option.nw,0) || ... + strcmpi(option.nw,'Off') || strcmpi(option.nw,'No') + option.win = 0; + elseif isequal(option.nw,1) || strcmpi(option.nw,'On') || ... + strcmpi(option.nw,'Yes') + option.win = 'hanning'; + else + option.win = postoption.nw; + end + end + + coeff = cell(1,length(sig)); + lags = cell(1,length(sig)); + wind = cell(1,length(sig)); + for k = 1:length(sig) + s = sig{k}; + p = pos{k}; + fk = f{k}; + if iscell(option.max) + mi = option.min{k}; + ma = option.max{k}; + else + mi = option.min; + ma = option.max; + end + coeffk = cell(1,length(s)); + lagsk = cell(1,length(s)); + windk = cell(1,length(s)); + for l = 1:length(s) + sl = s{l}; + sl(isnan(sl)) = 0; + if option.ni + mxsl = repmat(max(sl),[size(sl,1),1,1]); + mnsl = repmat(min(sl),[size(sl,1),1,1]); + sl = (sl-mnsl)./(mxsl-mnsl); + end + pl = p{l}; + pl = pl-repmat(pl(1,:,:),[size(pl,1),1,1]); + ls = size(sl,1); + if mi + misp = find(pl(:,1,1)>=mi); + if isempty(misp) + warning('WARNING IN MIRAUTOCOR: The specified range of delays exceeds the temporal length of the signal.'); + disp('Minimum delay set to zero.') + misp = 1; % misp is the lowest index of the lag range + mi = 0; + else + misp = misp(1); + end + else + misp = 1; + end + if ma + masp = find(pl(:,1,1)>=ma); + if isempty(masp) + masp = Inf; + else + masp = masp(1); + end + else + masp = Inf; + end + masp = min(masp,ceil(ls/2)); + if masp <= misp + if size(sl,2) > 1 + warning('WARNING IN MIRAUTOCOR: Frame length is too small.'); + else + warning('WARNING IN MIRAUTOCOR: The audio sequence is too small.'); + end + display('The autocorrelation is not defined for this range of delays.'); + end + sl = center(sl); + if not(ischar(option.win)) || strcmpi(option.win,'Rectangular') + kw = ones(size(sl)); + else + N = size(sl,1); + winf = str2func(option.win); + try + w = window(winf,N); + catch + if strcmpi(option.win,'hamming') + disp('Signal Processing Toolbox does not seem to be installed. Recompute the hamming window manually.'); + w = 0.54 - 0.46 * cos(2*pi*(0:N-1)'/(N-1)); + else + warning(['WARNING in MIRAUTOCOR: Unknown windowing function ',option.win,' (maybe Signal Processing Toolbox is not installed).']); + disp('No windowing performed.') + w = ones(size(sl,1),1); + end + end + kw = repmat(w,[1,size(sl,2),size(sl,3)]); + sl = sl.* kw; + end + + if strcmpi(option.scaleopt,'coeff') + scaleopt = 'none'; + else + scaleopt = option.scaleopt; + end + c = zeros(masp,size(sl,2),size(sl,3)); + for i = 1:size(sl,2) + for j = 1:size(sl,3) + if option.gener == 2 + cc = xcorr(sl(:,i,j),masp-1,scaleopt); + c(:,i,j) = cc(masp:end); + else + ss = abs(fft(sl(:,i,j))); + ss = ss.^option.gener; + cc = ifft(ss); + ll = (0:masp-1); + c(:,i,j) = cc(ll+1); + end + end + if strcmpi(option.scaleopt,'coeff') && option.gener == 2 + % to be adapted to generalized autocor + c(:,i,:) = c(:,i,:)/xcorr(sum(sl(:,i,:),3),0); + % This is a kind of generalization of the 'coeff' + % normalization for multi-channels signals. In the + % original 'coeff' option, the autocorrelation at zero + % lag is identically 1.0. In this multi-channels + % version, the autocorrelation at zero lag is such that + % the sum over channels becomes identically 1.0. + end + end + coeffk{l} = c(misp:end,:,:); + pl = pl(find(pl(:,1,1) >=mi),:,:); + lagsk{l} = pl(1:min(size(coeffk{l},1),size(pl,1)),:,:); + windk{l} = kw; + end + coeff{k} = coeffk; + lags{k} = lagsk; + wind{k} = windk; + end + a = set(a,'Coeff',coeff,'Delay',lags,'Window',wind); + if not(isempty(postoption)) + a = post(a,postoption); + end +end + + +function a = post(a,option) +debug = 0; +coeff = get(a,'Coeff'); +lags = get(a,'Delay'); +wind = get(a,'Window'); +freq = option.fr && not(get(a,'FreqDomain')); +if isequal(option.e,1) + option.e = 2:10; +end +if max(option.e) > 1 + pa = mirpeaks(a,'NoBegin','NoEnd','Contrast',.01); + va = mirpeaks(a,'Valleys','Contrast',.01); + pv = get(pa,'PeakVal'); + vv = get(va,'PeakVal'); +end +for k = 1:length(coeff) + for l = 1:length(coeff{k}) + c = coeff{k}{l}; % Coefficients of autocorrelation + t = lags{k}{l}; % Delays of autocorrelation + if not(isempty(c)) + if not(isequal(option.nw,0) || strcmpi(option.nw,'No') || ... + strcmpi(option.nw,'Off') || a.normalwindow) % 'NormalWindow' option + xw = zeros(size(c)); + lc = size(c,1); + for j = 1:size(c,3) + for i = 1:size(c,2) + xwij = xcorr(wind{k}{l}(:,i,j),lc,'coeff'); + xw(:,i,j) = xwij(lc+2:end); + end + end + c = c./ xw; + a.normalwindow = 1; + end + if ischar(option.reso) && ... + (strcmpi(option.reso,'ToiviainenSnyder') || ... + strcmpi(option.reso,'Toiviainen') || ... + strcmpi(option.reso,'vanNoorden')) + if isa(a,'mirautocor') && get(a,'FreqDomain') + ll = 1./t; + else + ll = t; + end + if not(option.resocenter) + option.resocenter = .5; + end + if strcmpi(option.reso,'ToiviainenSnyder') || ... + strcmpi(option.reso,'Toiviainen') + w = max(1 - 0.25*(log2(max(ll,1e-12)/option.resocenter)).^2, 0); + elseif strcmpi(option.reso,'vanNoorden') + f0=2.193; b=option.resocenter; + f=1./ll; a1=(f0*f0-f.*f).^2+b*f.^2; a2=f0^4+f.^4; + w=(1./sqrt(a1))-(1./sqrt(a2)); + end + if max(w) == 0 + warning('The resonance curve, not defined for this range of delays, will not be applied.') + else + w = w/max(w); + c = c.* repmat(w,[1,size(c,2),size(c,3)]); + end + end + if option.h + c = hwr(c); + end + if max(option.e) > 1 + if a.freq + freq = 1; + for i = 1:size(c,3) + c(:,:,i) = flipud(c(:,:,i)); + end + t = flipud(1./t); + end + + for g = 1:size(c,2) + for h = 1:size(c,3) + cgh = c(:,g,h); + if length(cgh)>1 + pvk = pv{k}{l}{1,g,h}; + mv = []; + if not(isempty(pvk)) + mp = min(pv{k}{l}{1,g,h}); %Lowest peak + vvv = vv{k}{l}{1,g,h}; %Valleys + mv = vvv(find(vvv<mp,1,'last')); + %Highest valley below the lowest peak + + if not(isempty(mv)) + cgh = cgh-mv; + end + end + cgh2 = cgh; + tgh2 = t(:,g,1); + coef = cgh(2)-cgh(1); % initial slope of the autocor curve + tcoef = tgh2(2)-tgh2(1); + deter = 0; + inter = 0; + + repet = find(not(diff(tgh2))); % Avoid bug if repeated x-values + if repet + warning('WARNING in MIRAUTOCOR: Two successive samples have exactly same temporal position.'); + tgh2(repet+1) = tgh2(repet)+1e-12; + end + + if coef < 0 + % initial descending slope removed + deter = find(diff(cgh2)>0,1)-1; + % number of removed points + if isempty(deter) + deter = 0; + end + cgh2(1:deter) = []; + tgh2(1:deter) = []; + coef = cgh2(2)-cgh2(1); + end + + if coef > 0 + % initial ascending slope prolonged to the left + % until it reaches the x-axis + while cgh2(1) > 0 + coef = coef*1.1; + % the further to the left, ... + % the more ascending is the slope + % (not sure it always works, though...) + inter = inter+1; + % number of added points + cgh2 = [cgh2(1)-coef; cgh2]; + tgh2 = [tgh2(1)-tcoef; tgh2]; + end + cgh2(1) = 0; + end + + for i = option.e % Enhancing procedure + % option.e is the list of scaling factors + % i is the scaling factor + if i + be = find(tgh2 & tgh2/i >= tgh2(1),1); + % starting point of the substraction + % on the X-axis + + if not(isempty(be)) + ic = interp1(tgh2,cgh2,tgh2/i); + % The scaled autocorrelation + ic(1:be-1) = 0; + ic(find(isnan(ic))) = Inf; + % All the NaN values are changed + % into 0 in the resulting curve + ic = max(ic,0); + + if debug + hold off,plot(tgh2,cgh2) + end + + cgh2 = cgh2 - ic; + % The scaled autocorrelation + % is substracted to the initial one + + cgh2 = max(cgh2,0); + % Half-wave rectification + + if debug + hold on,plot(tgh2,ic,'r') + hold on,plot(tgh2,cgh2,'g') + drawnow + figure + end + end + end + end + + % The temporary modifications are + % removed from the final curve + if inter>=deter + c(:,g,h) = cgh2(inter-deter+1:end); + if not(isempty(mv)) + c(:,g,h) = c(:,g,h) + mv; + end + else + c(:,g,h) = [zeros(deter-inter,1);cgh2]; + end + end + end + end + end + if freq + if t(1,1) == 0 + c = c(2:end,:,:); + t = t(2:end,:,:); + end + for i = 1:size(c,3) + c(:,:,i) = flipud(c(:,:,i)); + end + t = flipud(1./t); + end + coeff{k}{l} = c; + lags{k}{l} = t; + end + end +end +a = set(a,'Coeff',coeff,'Delay',lags,'Freq'); +if freq + a = set(a,'FreqDomain',1,'Abs','frequency (Hz)'); +end + + +function [y orig] = eachchunk(orig,option,missing,postchunk) +option.scaleopt = 'none'; +y = mirautocor(orig,option); + + +function y = combinechunk(old,new) +do = get(old,'Data'); +do = do{1}{1}; +dn = get(new,'Data'); +dn = dn{1}{1}; +if abs(size(dn,1)-size(do,1)) <= 2 % Probleme of border fluctuation + mi = min(size(dn,1),size(do,1)); + dn = dn(1:mi,:,:); + do = do(1:mi,:,:); +elseif length(dn) < length(do) + dn(length(do),:,:) = 0; % Zero-padding +end +y = set(old,'ChunkData',do+dn);