《MATLAB语音信号分析与合成（第二版）》语音工具包

function [imf,ort,nbit] = emd_online(x,t,stop,nbpresift,tst,tst2) % EMD_ONLINE (On Line Empirical Mode Decomposition) computes on-line an EMD % % stopping criterion for sifting : % at each point : mean amplitude < threshold*envelope amplitude % if mean amplitude > max(envelope amplitude)/tolerance % & % at each point : mean amplitude < threshold2*enveloppe amplitude % if mean amplitude > max(envelope amplitude)/tolerance2 % %[imf,ort,nbits] = EMD_ONLINE(x,t,stop,nbpresift,tst,tst2) % inputs: % - x : analyzed signal % - t (optional) : sampling times (default : 1:length(x)) % - stop (optional) : threshold, and threshold2 (optional) % tolerance, and tolerance2 (both optional) % for sifting stopping criterion % default : [0.05,0.5,20,100] % - nbpresift (optional) : number of sifting by pieces iterations (default 4) % - tst (optional) : if equals to 1 shows sifting steps % - tst2 (optional) : if equals to 1 shows sifting by pieces steps % % outputs: % - imf : intrinsic mode functions (last line = residual) % - ort : index of orthogonality % - nbits : number of iterations for each mode % % calls: % - extr (finds extrema and zero-crossings) % - io : computes the index of orthogonality % % G. Rilling, July 2002 DEFSTOP = [0.05,0.5,20,100];% default parameters for sifting stop NBPRESIFT = 4;%number of sifting iterations per block if(nargin==1) t = 1:length(x); stop = DEFSTOP; tst = 0; tst2 = 0; end if(nargin==2) stop = DEFSTOP; tst = 0; tst2 = 0; end if (nargin==3) tst=0; tst2 = 0; end if (nargin==4) tst=0; tst2 = 0; end if (nargin==5) tst2 = 0; end if nargin > 3 NBPRESIFT = nbpresift; end S = size(x); if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2) error('x must have only one row or one column') end if S(1) > 1 x = x'; end S = size(t); if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2) error('t must have only one row or one column') end if S(1) > 1 t = t'; end if (length(t)~=length(x)) error('x and t must have the same length') end S = size(stop); if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 4) | (S(2) > 4) | (length(S) > 2) error('stop must have only one row or one column of max four elements') end if S(1) > 1 stop = stop'; S = size(stop); end if S(2) < 4 stop(4) = DEFSTOP(4); end if S(2) < 3 stop(3) = DEFSTOP(3); end if S(2) < 2 stop(2) = DEFSTOP(2); end if S(2) == 1 stop=[stop, DEFSTOP(2)]; end sd = stop(1); sd2 = stop(2); tol = stop(3); tol2 = stop(4); if tst figure figures(1) = gcf; figure figures(3) = gcf; end if tst2 figure figures(2) = gcf; end MAXITERATIONS=10000; LARGMIN = 5; NBSYM = 2;% maximum number of symmetrized points for interpolations LARGTRANS = 10; LARGTRANSPS = 5; PAS = 20; STEP = 5;% maximal number of iterations on a mode LX = length(x); % for display sdt(LX) = 0; sdt = sdt+sd; sd2t(LX) = 0; sd2t = sd2t+sd2; % number of minima and maxima on the considered zone lm = 0; lM = 0; % number of minima and maxima right of the considered zone, % after "stop" or "stopps" lmr = 0; lMr = 0; % same, but left before "start" lml = 0; lMl = 0; % total number of extrema, left and right nem = 0; nemr = 0; neml = 0; k = 1; nbit = 0; % number of modes, and number of modes on which block siftings are completed nbmodes = 1; nbmodes_psdone = 0; start = 1; % end of the constant part of the window stop = min(PAS+1,LX); % start and end of the considered zone stopr = 1; startl = 1; % end of available data on the considered zone fin = 1; % start of the considered zone for block sifting limpsl(1,1:NBPRESIFT) = 1; % start and end of segment to which sifting is applied startps = 1; stopps(1,NBPRESIFT) = 0; stopps = stopps + 1; % end of available data for block sifting finps(1,1:NBPRESIFT) = 1; finps(1,1) = 10*PAS; % tests if all data are available for an iteration of block sifting lafinps(1,NBPRESIFT) = 0; % allows to interrupt a mode extraction to process to the next mode % interrupts also if not enough data available suspps(1,NBPRESIFT) = 0; % tests for the termination of one iteration of block sifting stoptestps(1,NBPRESIFT) = 0; indmin = []; indmax = []; % tests if all data are available for on-line sifting lafin = 0; % tests if a mode is entirely extracted stoptest = 0; % allows to interrupt a mode extraction to process to the next mode % interrupts also if not enough data available susp = 0; % tells if the considered zone has to be moved forward for having enough extrema needextr = 1; % idem for block sifting needextrps(1,1:NBPRESIFT) = 1; % tells how many iterations of block sifting have been initiated nbstartedpresift = 1; % modes concerned by block sifting mps = x; % mode concerned by on-line sifting m(LX) = 0; trig = 0; if tst | tst2 disp('appuyer sur une touche pour commencer') pause end while sum(stoptest) < nbmodes % global loop for k = 1:nbmodes nsteps = 0; waittest = 0; if k == 1 & trig suspps(1,1) = 0; trig = 0; end while sum(stoptestps(k,:)) < NBPRESIFT & ~waittest & sum(suspps(k,:)) < nbstartedpresift(k) % boucle de presifting for i = 1:nbstartedpresift(k) if needextrps(k,i) == 1 [indmintmp,indmaxtmp] = extr(mps(k,max([(limpsl(k,i)-1),1]):finps(k,i),i)); nb = sum(indmintmp > stopps(k,i))+sum(indmaxtmp > stopps(k,i)); stoprps(k,i) = finps(k,i); if nb < 8*LARGTRANSPS & finps(k,i) < LX suspps(k,i) = 1; if k == 1 & i == 1 finps(1,1) = min(LX,finps(1,1) + 10*PAS); [indmintmp,indmaxtmp] = extr(mps(k,max([(limpsl(k,i)-1),1]):finps(k,i),i)); nb = sum(indmintmp > stopps(k,i))+sum(indmaxtmp > stopps(k,i)); stoprps(k,i) = finps(k,i); trig = 1; end else lmt = length(indmintmp); lMt = length(indmaxtmp); if lmt > 0 indminps(k,1:lmt,i) = indmintmp + max([(limpsl(k,i)-1),1])-1; end if lMt > 0 indmaxps(k,1:lMt,i) = indmaxtmp + max([(limpsl(k,i)-1),1])-1; end if lmt < size(indminps,2) indminps(k,length(indmintmp)+1:end,i) = 0; end if lMt < size(indmaxps,2) indmaxps(k,length(indmaxtmp)+1:end,i) = 0; end needextrps(k,i) = 0; end if stoprps(k,i) >= LX lafinps(k,i) = 1; needextrps(k,i) = 0; end end if ~suspps(k,i) curindminps = indminps(k,find(indminps(k,:,i) >= limpsl(k,i)),i); curindmaxps = indmaxps(k,find(indmaxps(k,:,i) >= limpsl(k,i)),i); nemps = length(curindminps) + length(curindmaxps); end % loop of block (pre)sifting while (~needextrps(k,i) | lafinps(k,i)) & ~stoptestps(k,i) & ~waittest & ~suspps(k,i) if nemps < 3 & lafinps(k,i) stoptestps(k,:) = 1; stoptest(k) = 1; m(k,:) = mps(k,:,i); if i > 1 m(k+1,:) = x - sum(m(1:k,:)); end break end if limpsl(k,i) == 1 startps(k,i) = 1; else startps(k,i) = stopps(k,i); end if lafinps(k,i) stopps(k,i) = LX; stoptestps(k,i) = 1; else stopps(k,i) = min(curindminps(max([1,end - LARGTRANSPS+1])),curindmaxps(max([1,end - LARGTRANSPS+1]))); % si ~lafinps(k,i) end if startps(k,i) == stopps(k,i) pause needextrps(k,i) = 1; break end lmr = sum(curindminps > stopps(k,i)); lMr = sum(curindmaxps > stopps(k,i)); nemrps(k,i) = lmr + lMr; if nemrps(k,i) < 8*LARGTRANSPS needextrps(k,i) = 1; end if limpsl(k,i) == 1 margeml = 0; margeMl = 0; if curindmaxps(1) < curindminps(1) if mps(k,1,i) > mps(k,curindminps(1),i) lmax = fliplr(curindmaxps(2:min(end,NBSYM+1))); lmin = fliplr(curindminps(1:min(end,NBSYM))); lsym = curind