comparison toolboxes/MIRtoolbox1.3.2/MIRToolbox/@mirautocor/mirautocor.m @ 0:e9a9cd732c1e tip

first hg version after svn

+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);