Daniel@0: function varargout = mirenvelope(orig,varargin) Daniel@0: % e = mirenvelope(x) extracts the envelope of x, showing the global shape Daniel@0: % of the waveform. Daniel@0: % mirenvelope(...,m) specifies envelope extraction method. Daniel@0: % Possible values: Daniel@0: % m = 'Filter' uses a low-pass filtering. (Default strategy) Daniel@0: % m = 'Spectro' uses a spectrogram. Daniel@0: % Daniel@0: % Options related to the 'Filter' method: Daniel@0: % mirenvelope(...,'Hilbert'): performs a preliminary Hilbert Daniel@0: % transform. Daniel@0: % mirenvelope(...,'PreDecim',N) downsamples by a factor N>1, where Daniel@0: % N is an integer, before the low-pass filtering (Klapuri, 1999). Daniel@0: % Default value: N = 1. Daniel@0: % mirenvelope(...,'Filtertype',f) specifies the filter type. Daniel@0: % Possible values are: Daniel@0: % f = 'IIR': filter with one autoregressive coefficient Daniel@0: % (default) Daniel@0: % f = 'HalfHann': half-Hanning (raised cosine) filter Daniel@0: % (Scheirer, 1998) Daniel@0: % Option related to the 'IIR' option: Daniel@0: % mirenvelope(...,'Tau',t): time constant of low-pass filter in Daniel@0: % seconds. Daniel@0: % Default value: t = 0.02 s. Daniel@0: % mirenvelope(...,'PostDecim',N) downsamples by a factor N>1, where Daniel@0: % N is an integer, after the low-pass filtering. Daniel@0: % Default value: N = 16 if 'PreDecim' is not used, else N = 1. Daniel@0: % mirenvelope(...,'Trim'): trims the initial ascending phase of the Daniel@0: % curves related to the transitory state. Daniel@0: % Daniel@0: % Options related to the 'Spectro' method: Daniel@0: % mirenvelope(...,b) specifies whether the frequency range is further Daniel@0: % decomposed into bands. Possible values: Daniel@0: % b = 'Freq': no band decomposition (default value) Daniel@0: % b = 'Mel': Mel-band decomposition Daniel@0: % b = 'Bark': Bark-band decomposition Daniel@0: % b = 'Cents': decompositions into cents Daniel@0: % mirenvelope(...,'Frame',...) specifies the frame configuration. Daniel@0: % Default value: length: .1 s, hop factor: 10 %. Daniel@0: % mirenvelope(...,'UpSample',N) upsamples by a factor N>1, where Daniel@0: % N is an integer. Daniel@0: % Default value if 'UpSample' called: N = 2 Daniel@0: % mirenvelope(...,'Complex') toggles on the 'Complex' method for the Daniel@0: % spectral flux computation. Daniel@0: % Daniel@0: % Other available for all methods: Daniel@0: % mirenvelope(...,'Sampling',r): resamples to rate r (in Hz). Daniel@0: % 'Down' and 'Sampling' options cannot therefore be combined. Daniel@0: % mirenvelope(...,'Halfwave'): performs a half-wave rectification. Daniel@0: % mirenvelope(...,'Center'): centers the extracted envelope. Daniel@0: % mirenvelope(...,'HalfwaveCenter'): performs a half-wave Daniel@0: % rectification on the centered envelope. Daniel@0: % mirenvelope(...,'Log'): computes the common logarithm (base 10) of Daniel@0: % the envelope. Daniel@0: % mirenvelope(...,'Mu',mu): computes the logarithm of the Daniel@0: % envelope, before the eventual differentiation, using a mu-law Daniel@0: % compression (Klapuri, 2006). Daniel@0: % Default value for mu: 100 Daniel@0: % mirenvelope(...,'Log'): computes the logarithm of the envelope. Daniel@0: % mirenvelope(...,'Power'): computes the power (square) of the Daniel@0: % envelope. Daniel@0: % mirenvelope(...,'Diff'): computes the differentation of the Daniel@0: % envelope, i.e., the differences between successive samples. Daniel@0: % mirenvelope(...,'HalfwaveDiff'): performs a half-wave Daniel@0: % rectification on the differentiated envelope. Daniel@0: % mirenvelope(...,'Normal'): normalizes the values of the envelope by Daniel@0: % fixing the maximum value to 1. Daniel@0: % mirenvelope(...,'Lambda',l): sums the half-wave rectified envelope Daniel@0: % with the non-differentiated envelope, using the respective Daniel@0: % weight 01 % Already in bands (channels in 3d dim) Daniel@0: d{h}{i} = permute(sum(d{h}{i}),[2 1 3]); Daniel@0: ph{h}{i} = permute(ph{h}{i},[2 1 3]); Daniel@0: else % Simple spectrogram, frequency range sent to 3d dim Daniel@0: d{h}{i} = permute(d{h}{i},[2 3 1]); Daniel@0: ph{h}{i} = permute(ph{h}{i},[2 3 1]); Daniel@0: end Daniel@0: p{h}{i} = mean(fp{h}{i})'; Daniel@0: if not(sr{h}) && size(fp{h}{i},2)>1 Daniel@0: sr{h} = 1/(fp{h}{i}(1,2)-fp{h}{i}(1,1)); Daniel@0: end Daniel@0: end Daniel@0: if not(sr{h}) Daniel@0: warning('WARNING IN MIRENVELOPE: The frame decomposition did not succeed. Either the input is of too short duration, or the chunk size is too low.'); Daniel@0: end Daniel@0: ch{h} = (1:size(d{h}{1},3))'; Daniel@0: end Daniel@0: e.downsampl = 0; Daniel@0: e.hwr = 0; Daniel@0: e.diff = 0; Daniel@0: %e.todelete = {}; Daniel@0: e.method = 'Spectro'; Daniel@0: e.phase = ph; Daniel@0: e = class(e,'mirenvelope',mirtemporal(orig)); Daniel@0: e = set(e,'Title','Envelope','Data',d,'Pos',p,... Daniel@0: 'Sampling',sr,'Channels',ch); Daniel@0: postoption.trim = 0; Daniel@0: postoption.ds = 0; Daniel@0: e = post(e,postoption); Daniel@0: else Daniel@0: if isnan(option.zp) Daniel@0: if strcmpi(option.filter,'IIR') Daniel@0: option.zp = 1; Daniel@0: else Daniel@0: option.zp = 0; Daniel@0: end Daniel@0: end Daniel@0: if option.zp == 1 Daniel@0: option.decim = 0; Daniel@0: end Daniel@0: e.downsampl = 1; Daniel@0: e.hwr = 0; Daniel@0: e.diff = 0; Daniel@0: e.method = option.filter; Daniel@0: e.phase = {}; Daniel@0: e = class(e,'mirenvelope',mirtemporal(orig)); Daniel@0: e = purgedata(e); Daniel@0: e = set(e,'Title','Envelope'); Daniel@0: sig = get(e,'Data'); Daniel@0: x = get(e,'Pos'); Daniel@0: sr = get(e,'Sampling'); Daniel@0: disp('Extracting envelope...') Daniel@0: d = cell(1,length(sig)); Daniel@0: [state e] = gettmp(orig,e); Daniel@0: for k = 1:length(sig) Daniel@0: if option.decim Daniel@0: sr{k} = sr{k}/option.decim; Daniel@0: end Daniel@0: if strcmpi(option.filter,'IIR') Daniel@0: a2 = exp(-1/(option.tau*sr{k})); % filter coefficient Daniel@0: a = [1 -a2]; Daniel@0: b = 1-a2; Daniel@0: elseif strcmpi(option.filter,'HalfHann') Daniel@0: a = 1; Daniel@0: b = hann(sr{k}*.4); Daniel@0: b = b(ceil(length(b)/2):end); Daniel@0: end Daniel@0: d{k} = cell(1,length(sig{k})); Daniel@0: for i = 1:length(sig{k}) Daniel@0: sigi = sig{k}{i}; Daniel@0: if option.zp == 2 Daniel@0: sigi = flipdim(sigi,1); Daniel@0: end Daniel@0: if option.hilb Daniel@0: try Daniel@0: for h = 1:size(sigi,2) Daniel@0: for j = 1:size(sigi,3) Daniel@0: sigi(:,h,j) = hilbert(sigi(:,h,j)); Daniel@0: end Daniel@0: end Daniel@0: catch Daniel@0: disp('Signal Processing Toolbox does not seem to be installed. No Hilbert transform.'); Daniel@0: end Daniel@0: end Daniel@0: sigi = abs(sigi); Daniel@0: Daniel@0: if option.decim Daniel@0: dsigi = zeros(ceil(size(sigi,1)/option.decim),... Daniel@0: size(sigi,2),size(sigi,3)); Daniel@0: for f = 1:size(sigi,2) Daniel@0: for c = 1:size(sigi,3) Daniel@0: dsigi(:,f,c) = decimate(sigi(:,f,c),option.decim); Daniel@0: end Daniel@0: end Daniel@0: sigi = dsigi; Daniel@0: clear dsigi Daniel@0: x{k}{i} = x{k}{i}(1:option.decim:end,:,:); Daniel@0: end Daniel@0: Daniel@0: % tmp = filtfilt(1-a,[1 -a],sigi); % zero-phase IIR filter for smoothing the envelope Daniel@0: Daniel@0: % Manual filtfilt Daniel@0: emptystate = isempty(state); Daniel@0: tmp = zeros(size(sigi)); Daniel@0: for c = 1:size(sigi,3) Daniel@0: if emptystate Daniel@0: [tmp(:,:,c) state(:,c,1)] = filter(b,a,sigi(:,:,c)); Daniel@0: else Daniel@0: [tmp(:,:,c) state(:,c,1)] = filter(b,a,sigi(:,:,c),... Daniel@0: state(:,c,1)); Daniel@0: end Daniel@0: end Daniel@0: Daniel@0: tmp = max(tmp,0); % For security reason... Daniel@0: if option.zp == 2 Daniel@0: tmp = flipdim(tmp,1); Daniel@0: end Daniel@0: d{k}{i} = tmp; Daniel@0: %td{k} = round(option.tau*sr{k}*1.5); Daniel@0: end Daniel@0: end Daniel@0: e = set(e,'Data',d,'Pos',x,'Sampling',sr); %,'ToDelete',td Daniel@0: e = settmp(e,state); Daniel@0: if not(option.zp == 2) Daniel@0: e = post(e,postoption); Daniel@0: end Daniel@0: end Daniel@0: if isfield(option,'presel') && ischar(option.presel) && ... Daniel@0: strcmpi(option.presel,'Klapuri06') Daniel@0: e = mirsum(e,'Adjacent',10); Daniel@0: end Daniel@0: Daniel@0: Daniel@0: function e = post(e,postoption) Daniel@0: if isempty(postoption) Daniel@0: return Daniel@0: end Daniel@0: if isfield(postoption,'lambda') && not(postoption.lambda) Daniel@0: postoption.lambda = 1; Daniel@0: end Daniel@0: d = get(e,'Data'); Daniel@0: tp = get(e,'Time'); Daniel@0: sr = get(e,'Sampling'); Daniel@0: ds = get(e,'DownSampling'); Daniel@0: ph = get(e,'Phase'); Daniel@0: for k = 1:length(d) Daniel@0: if isfield(postoption,'sampling') Daniel@0: if postoption.sampling Daniel@0: newsr = postoption.sampling; Daniel@0: elseif isfield(postoption,'ds') && postoption.ds>1 Daniel@0: newsr = sr{k}/postoption.ds; Daniel@0: else Daniel@0: newsr = sr{k}; Daniel@0: end Daniel@0: end Daniel@0: if isfield(postoption,'up') && postoption.up Daniel@0: [z,p,gain] = butter(6,10/newsr/postoption.up*2,'low'); Daniel@0: [sos,g] = zp2sos(z,p,gain); Daniel@0: Hd = dfilt.df2tsos(sos,g); Daniel@0: end Daniel@0: for i = 1:length(d{k}) Daniel@0: if isfield(postoption,'sampling') && postoption.sampling Daniel@0: if and(sr{k}, not(sr{k} == postoption.sampling)) Daniel@0: dk = d{k}{i}; Daniel@0: for j = 1:size(dk,3) Daniel@0: if not(sr{k} == round(sr{k})) Daniel@0: mirerror('mirenvelope','The ''Sampling'' postoption cannot be used after using the ''Down'' postoption.'); Daniel@0: end Daniel@0: rk(:,:,j) = resample(dk(:,:,j),postoption.sampling,sr{k}); Daniel@0: end Daniel@0: d{k}{i} = rk; Daniel@0: tp{k}{i} = repmat((0:size(d{k}{i},1)-1)',... Daniel@0: [1 1 size(tp{k}{i},3)])... Daniel@0: /postoption.sampling + tp{k}{i}(1,:,:); Daniel@0: if not(iscell(ds)) Daniel@0: ds = cell(length(d)); Daniel@0: end Daniel@0: ds{k} = round(sr{k}/postoption.sampling); Daniel@0: end Daniel@0: elseif isfield(postoption,'ds') && postoption.ds>1 Daniel@0: if not(postoption.ds == round(postoption.ds)) Daniel@0: mirerror('mirenvelope','The ''Down'' sampling rate should be an integer.'); Daniel@0: end Daniel@0: ds = postoption.ds; Daniel@0: tp{k}{i} = tp{k}{i}(1:ds:end,:,:); % Downsampling... Daniel@0: d{k}{i} = d{k}{i}(1:ds:end,:,:); Daniel@0: end Daniel@0: if isfield(postoption,'sampling') Daniel@0: if not(strcmpi(e.method,'Spectro')) && postoption.trim Daniel@0: tdk = round(newsr*.1); Daniel@0: d{k}{i}(1:tdk,:,:) = repmat(d{k}{i}(tdk,:,:),[tdk,1,1]); Daniel@0: d{k}{i}(end-tdk+1:end,:,:) = repmat(d{k}{i}(end-tdk,:,:),[tdk,1,1]); Daniel@0: end Daniel@0: if postoption.log Daniel@0: d{k}{i} = log10(d{k}{i}); Daniel@0: end Daniel@0: if postoption.mu Daniel@0: dki = max(0,d{k}{i}); Daniel@0: mu = postoption.mu; Daniel@0: dki = log(1+mu*dki)/log(1+mu); Daniel@0: dki(~isfinite(d{k}{i})) = NaN; Daniel@0: d{k}{i} = dki; Daniel@0: end Daniel@0: if postoption.power Daniel@0: d{k}{i} = d{k}{i}.^2; Daniel@0: end Daniel@0: if postoption.up Daniel@0: dki = zeros(size(d{k}{i},1).*postoption.up,... Daniel@0: size(d{k}{i},2),size(d{k}{i},3)); Daniel@0: dki(1:postoption.up:end,:,:) = d{k}{i}; Daniel@0: dki = filter(Hd,[dki;... Daniel@0: zeros(6,size(d{k}{i},2),size(d{k}{i},3))]); Daniel@0: d{k}{i} = dki(1+ceil(6/2):end-floor(6/2),:,:); Daniel@0: tki = zeros(size(tp{k}{i},1).*postoption.up,... Daniel@0: size(tp{k}{i},2),... Daniel@0: size(tp{k}{i},3)); Daniel@0: dt = repmat((tp{k}{i}(2)-tp{k}{i}(1))... Daniel@0: /postoption.up,... Daniel@0: [size(tp{k}{i},1),1,1]); Daniel@0: for j = 1:postoption.up Daniel@0: tki(j:postoption.up:end,:,:) = tp{k}{i}+dt*(j-1); Daniel@0: end Daniel@0: tp{k}{i} = tki; Daniel@0: newsr = sr{k}*postoption.up; Daniel@0: end Daniel@0: if (postoption.diffhwr || postoption.diff) && ... Daniel@0: not(get(e,'Diff')) Daniel@0: tp{k}{i} = tp{k}{i}(1:end-1,:,:); Daniel@0: order = max(postoption.diffhwr,postoption.diff); Daniel@0: if postoption.complex Daniel@0: dph = diff(ph{k}{i},2); Daniel@0: dph = dph/(2*pi);% - round(dph/(2*pi)); Daniel@0: ddki = sqrt(d{k}{i}(3:end,:,:).^2 + d{k}{i}(2:end-1,:,:).^2 ... Daniel@0: - 2.*d{k}{i}(3:end,:,:)... Daniel@0: .*d{k}{i}(2:end-1,:,:)... Daniel@0: .*cos(dph)); Daniel@0: d{k}{i} = d{k}{i}(2:end,:,:); Daniel@0: tp{k}{i} = tp{k}{i}(2:end,:,:); Daniel@0: elseif order == 1 Daniel@0: ddki = diff(d{k}{i},1,1); Daniel@0: else Daniel@0: b = firls(order,[0 0.9],[0 0.9*pi],'differentiator'); Daniel@0: ddki = filter(b,1,... Daniel@0: [repmat(d{k}{i}(1,:,:),[order,1,1]);... Daniel@0: d{k}{i};... Daniel@0: repmat(d{k}{i}(end,:,:),[order,1,1])]); Daniel@0: ddki = ddki(order+1:end-order-1,:,:); Daniel@0: end Daniel@0: if postoption.diffhwr Daniel@0: ddki = hwr(ddki); Daniel@0: end Daniel@0: d{k}{i} = (1-postoption.lambda)*d{k}{i}(1:end-1,:,:)... Daniel@0: + postoption.lambda*sr{k}/10*ddki; Daniel@0: end Daniel@0: if postoption.aver Daniel@0: y = filter(ones(1,postoption.aver),1,... Daniel@0: [d{k}{i};zeros(postoption.aver,... Daniel@0: size(d{k}{i},2),... Daniel@0: size(d{k}{i},3))]); Daniel@0: d{k}{i} = y(1+ceil(postoption.aver/2):... Daniel@0: end-floor(postoption.aver/2),:,:); Daniel@0: end Daniel@0: if postoption.gauss Daniel@0: sigma = postoption.gauss; Daniel@0: gauss = 1/sigma/2/pi... Daniel@0: *exp(- (-4*sigma:4*sigma).^2 /2/sigma^2); Daniel@0: y = filter(gauss,1,[d{k}{i};zeros(4*sigma,1,size(d{k}{i},3))]); Daniel@0: y = y(4*sigma:end,:,:); Daniel@0: d{k}{i} = y(1:size(d{k}{i},1),:,:); Daniel@0: end Daniel@0: if postoption.chwr Daniel@0: d{k}{i} = center(d{k}{i}); Daniel@0: d{k}{i} = hwr(d{k}{i}); Daniel@0: end Daniel@0: if postoption.hwr Daniel@0: d{k}{i} = hwr(d{k}{i}); Daniel@0: end Daniel@0: if postoption.c Daniel@0: d{k}{i} = center(d{k}{i}); Daniel@0: end Daniel@0: if postoption.norm Daniel@0: d{k}{i} = d{k}{i}./repmat(max(abs(d{k}{i})),... Daniel@0: [size(d{k}{i},1),1,1]); Daniel@0: end Daniel@0: end Daniel@0: end Daniel@0: if isfield(postoption,'sampling') Daniel@0: sr{k} = newsr; Daniel@0: end Daniel@0: end Daniel@0: if isfield(postoption,'ds') && postoption.ds>1 Daniel@0: e = set(e,'DownSampling',postoption.ds,'Sampling',sr); Daniel@0: elseif isfield(postoption,'sampling') && postoption.sampling Daniel@0: e = set(e,'DownSampling',ds,'Sampling',sr); Daniel@0: elseif isfield(postoption,'up') && postoption.up Daniel@0: e = set(e,'Sampling',sr); Daniel@0: end Daniel@0: if isfield(postoption,'sampling') Daniel@0: if postoption.hwr Daniel@0: e = set(e,'Halfwave',1); Daniel@0: end Daniel@0: if postoption.diff Daniel@0: e = set(e,'Diff',1,'Halfwave',0,'Title','Differentiated envelope'); Daniel@0: end Daniel@0: if postoption.diffhwr Daniel@0: e = set(e,'Diff',1,'Halfwave',1,'Centered',0); Daniel@0: end Daniel@0: if postoption.c Daniel@0: e = set(e,'Centered',1); Daniel@0: end Daniel@0: if postoption.chwr Daniel@0: e = set(e,'Halfwave',1,'Centered',1); Daniel@0: end Daniel@0: end Daniel@0: e = set(e,'Data',d,'Time',tp); Daniel@0: %if isfield(postoption,'frame') && isstruct(postoption.frame) Daniel@0: % e = mirframenow(e,postoption); Daniel@0: %end