dtwexpofn.m

非常好的数字处理教程

    n = N;
    m = M;
		
    while n > 1
        pathtype = psi(n,m);
				if pathtype == 0
					n
					m
					psi
				end
        delta = sum(policy{pathtype},1);
        deltan = delta(1);
        deltam = delta(2);
        if deltan == 1
            m = m-deltam;
            n = n-1;
            p = [p; n m];
        else
            for i=1:deltan-1
                n = n-1;
                p = [p; n m];
            end
            m = m-deltam;
            n = n-1;
            p = [p; n m];
        end
		end
		p = [p; 1 1];
    
    switch type
      case 'Cepstrum'
        ud.DTW.Cep.Qmax = Qmax;
        ud.DTW.Cep.Qmin = Qmin;
        ud.DTW.Cep.dXY = dXY;
        ud.DTW.Cep.D = D;
        ud.DTW.Cep.psi = psi;
        ud.DTW.Cep.p = p;
        ud.DTW.Cep.N = N;
        ud.DTW.Cep.M = M;
      case 'MFCC'
        ud.DTW.MFCC.Qmax = Qmax;
        ud.DTW.MFCC.Qmin = Qmin;
        ud.DTW.MFCC.dXY = dXY;
        ud.DTW.MFCC.D = D;
        ud.DTW.MFCC.psi = psi;
        ud.DTW.MFCC.p = p;
        ud.DTW.MFCC.N = N;
        ud.DTW.MFCC.M = M;
      case 'LPC'
        ud.DTW.LPC.Qmax = Qmax;
        ud.DTW.LPC.Qmin = Qmin;
        ud.DTW.LPC.dXY = dXY;
        ud.DTW.LPC.D = D;
        ud.DTW.LPC.psi = psi;
        ud.DTW.LPC.p = p;
        ud.DTW.LPC.N = N;
        ud.DTW.LPC.M = M;
    end
    
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function updateDTWPlot(ud)
   
    contents = get(ud.datatype,'String');
    
    switch contents{get(ud.datatype,'Value')}
      case 'Cepstrum'
        Qmax = ud.DTW.Cep.Qmax;
        Qmin = ud.DTW.Cep.Qmin;
        dXY = ud.DTW.Cep.dXY;
        D = ud.DTW.Cep.D;
        psi = ud.DTW.Cep.psi;
        p = ud.DTW.Cep.p;
        N = ud.DTW.Cep.N;
        M = ud.DTW.Cep.M;
      case 'MFCC'
        Qmax = ud.DTW.MFCC.Qmax;
        Qmin = ud.DTW.MFCC.Qmin;
        dXY = ud.DTW.MFCC.dXY;
        D = ud.DTW.MFCC.D;
        psi = ud.DTW.MFCC.psi;
        p = ud.DTW.MFCC.p;
        N = ud.DTW.MFCC.N;
        M = ud.DTW.MFCC.M;
      case 'LPC'
        Qmax = ud.DTW.LPC.Qmax;
        Qmin = ud.DTW.LPC.Qmin;
        dXY = ud.DTW.LPC.dXY;
        D = ud.DTW.LPC.D;
        psi = ud.DTW.LPC.psi;
        p = ud.DTW.LPC.p;
        N = ud.DTW.LPC.N;
        M = ud.DTW.LPC.M;
    end
    
    axes(ud.dtwplot);
    hold on;
    
    if get(ud.grid,'Value')
      x = [1:N];
      y = [1:M];
      for i=1:1:N
          dtwgrid((i-1)*M+1:i*M,1) = i;
          dtwgrid((i-1)*M+1:i*M,2) = y;
      end
    
      tpf = plot(dtwgrid(:,1),dtwgrid(:,2),'k.');
      set(tpf,'MarkerSize',4);
    end
    drawnow;
    if get(ud.contours,'Value')
      % Make plot of the surface
      Ds = D;
      Ds(find(D > 10^8))=0;
      graylevel = [0.5 0.5 0.5];
      %tpc = contour([0:N-1]+0.5,[0:M-1]+0.5,Ds',100,'Color',graylevel);
      tpc = contour([0:N-1]+0.5,[0:M-1]+0.5,Ds',100);
      %tpc = pcolor([0:N-1]+0.5,[0:M-1]+0.5,Ds');
      colormap('hot');
      %shading interp;
    end
    drawnow;
    if get(ud.search,'Value')
      % Find limiting frames on searching along reference
      for n=1:N
          mmax = min(Qmax*(n-1)+1+0.5,ceil(Qmin*(n-N)+M)+0.5);
          mmin = max(floor(Qmin*(n-1)+1)-0.5,Qmax*(n-N)+M-0.5);
          r(n,:) = [mmin mmax];
      end
      tpl = plot(r,'k-');
      set(tpl,'LineWidth',2);
    end
    drawnow;
    % Plot optimal path
    if get(ud.optpath,'Value')
      tpp = plot(p(:,1),p(:,2),'b');
      set(tpp,'LineWidth',2);
    end
    
%     legend([tpc(1),tpf,tpl(1),tpp], ...
%     'Distance Contours','Frame','Search Boundary','Optimal Path',2);
    
    set(ud.dtwplot,'XTickLabel','','XTick',[], ...
        'YTickLabel','','YTick',[]);
    axis([0.5 N 0.5 M]);
		tt = text(N*0.5,M*0.1,['Minimum Average Distance: ' num2str(dXY)], ...
      'BackgroundColor','white');
    set(tt,'FontSize',10,'FontWeight','bold');
    

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
 function audiodata = endpoint(audiodata, window_size, window_skip)

    dBthreshold = 25;
    zcthreshold_low = 28;
    zcthreshold = 30;
    
    shape = 'Hamming';
    y = audiodata.data;
    
    win = get_windowdata(window_size, shape);
    
    numwin = floor(size(y,1)/window_skip);
    % Make sure we have an integer number of windows
    if size(y,1) < ((numwin-1)*window_skip + window_size)
        diff = numwin*window_size + window_size - window_skip - length(y);
        y = [y; zeros(diff,1)];
    end
    
    % Get data
    % Gather log energy and zerocrossings
    maxlogen = -1000;
    for j = 1:numwin,
        start = (j-1)*window_skip+1; 
        stop = start+window_size-1;
        frame = y(start:stop);
        zc(j) = zero_crossings(frame);
        logen(j) = logenergy(frame,win);
        if logen(j) > maxlogen
            maxlogen = logen(j);
        end
    end

    % Normalize logen
    logen = (logen - maxlogen)';
    % Avg number per 10 ms interval
    zc = zc'*window_skip/(2*window_size);
    
    tz = [1:size(zc,1)].*window_skip/audiodata.Fs;
    
    logen_i = find(logen > -dBthreshold);
    if numel(logen_i) > 0
        % Add two frames to compensate for window size
        lower1 = logen_i(1)-1;
        higher1 = logen_i(end)+1;
    else
        lower1 = 1;
        higher1 = numel(tz);
    end

    % Revise with zerocrossings
    zc_i = find(zc(1:lower1) > zcthreshold);

    if ~isempty(zc_i)
        lower2 = zc_i(end);
        % Now find where it drops below zcthreshold_low
        zc_i2 = find(zc(1:lower2) < zcthreshold_low);
        if ~isempty(zc_i2)
            lower2 = zc_i2(end);
        end 
    else
        lower2 = lower1;
    end

    zc_i = find(zc(higher1:end) > zcthreshold);
    if ~isempty(zc_i)
        higher2 = zc_i(1)+higher1;
        % Now find where it drops below zcthreshold_low
        zc_i2 = find(zc(higher2:end) < zcthreshold_low);
        if ~isempty(zc_i2)
            higher2 = zc_i2(1)+higher2;
        end 
    else
        higher2 = higher1;
    end
    
    startpt = ceil(tz(lower2)*audiodata.Fs);
    endpt = floor(tz(higher2)*audiodata.Fs);
    audiodata.data = y(startpt:endpt,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function audiodata = getStatistics(audiodata, window_size, window_skip)
    % Get cepstrum coefficients
    audiodata.cepstrum = cepstrum(audiodata, window_size, window_skip);
    audiodata.MFCC = getmfcc(audiodata, window_size, window_skip);
    audiodata.LPC = getlpc(audiodata, window_size, window_skip);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function Nzeros = zero_crossings(signal)
    Nzeros = 0;
    for i=2:length(signal),
        if (sign(signal(i)) ~= sign(signal(i-1)))
            Nzeros = Nzeros + 1;
        end
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function logenergy_value = logenergy(signal,window)
    logenergy_value = 0;
    N = length(signal);
    signal = window.*signal;
    logenergy_value = 10*log10(1e-10+sum(signal.^2));
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function cepdata = cepstrum(audiodata, window_size, window_skip)
    
    y = audiodata.data;
    shape = 'Hamming';
    
    win = get_windowdata(window_size, shape);
    
    numwin = floor(size(y,1)/window_skip);
    % Make sure we have an integer number of windows
    if size(y,1) < ((numwin-1)*window_skip + window_size)
        diff = numwin*window_size + window_size - window_skip - length(y);
        y = [y; zeros(diff,1)];
    end
    
    N = window_size*4;  % to avoid aliasing
    for i=1:numwin
        startN = (i-1)*window_skip+1;
        s = win.*y(startN:startN+window_size-1);
        S = fft(s,N);
        SC = log(abs(S));% + j*angle(S);
        sc = (ifft(SC));
        % take first 16 coefficients except for first one
        cepdata(i,:) = real(sc(2:17));
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function mfccvals = getmfcc(audiodata, windowsize, windowskip)
    mfccvals = mfcc(audiodata.data, audiodata.Fs, windowsize, windowskip)';
    mfccvals = mfccvals(:,2:end);
		
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
% function lpcdata = getlpc(audiodata, window_size, window_skip)
%     y = audiodata.data;
%     shape = 'Hamming';
%     
%     win = get_windowdata(window_size, shape);
%     
%     numwin = floor(size(y,1)/window_skip);
%     % Make sure we have an integer number of windows
%     if size(y,1) < ((numwin-1)*window_skip + window_size)
%         diff = numwin*window_size + window_size - window_skip - length(y);
%         y = [y; zeros(diff,1)];
%     end
%     
%     p = 16;
%     for i=1:numwin
%       startN = (i-1)*window_skip+1;
%       signal = win.*y(startN:startN+window_size-1);
%       R = xcorr(signal,'coeff');
%       R = R(ceil(end/2):end);
%       R = R(1:p+1);
%       Rt = toeplitz(R(1:end-1));
%       alpha = inv(Rt)*R(2:end);
%       lpcdata(i,:) = alpha;
% 		end

		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
function lpcdata = getlpc(audiodata, window_size, window_skip)
    y = audiodata.data;
    shape = 'Hamming';
    
    win = get_windowdata(window_size, shape);
    
    numwin = floor(size(y,1)/window_skip);
    % Make sure we have an integer number of windows
    if size(y,1) < ((numwin-1)*window_skip + window_size)
        diff = numwin*window_size + window_size - window_skip - length(y);
        y = [y; zeros(diff,1)];
    end
    
    p = 16;
    for i=1:numwin
      startN = (i-1)*window_skip+1;
			signal = win.*y(startN:startN+window_size-1);
			R = xcorr(signal);
			R = R(ceil(end/2):end);
			R = R(1:p+1);
			Rt = toeplitz(R(1:end-1));
			a = [ 1; -inv(Rt)*R(2:end)];
			%a = lpc(signal,p)';

			% For likelihood-ratio distance calculation
			lpcdata.R(i,:) = R';
			%b = xcorr([-1 a'],'coeff');
			b = xcorr(-a);
			b = b(ceil(end/2):end);
			b = [b(1)/2; b(2:p+1)]';
			lpcdata.b(i,:) = b;
			%[temp lpcdata.sigma_p3(i)] = lpc(signal,p);
			%lpcdata.sigma_p(i) = a(2:end)'*Rt*a(2:end);
			lpcdata.sigma_p2(i) = 2*b*R(:);
		end
		
		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dist = distance(x,y)
% DIST = DISTANCE(X,Y)

    if size(x) ~= size(y)
        disp('x,y must be same size!');
        return;
    end
    
    if nargin < 2
        dist('Must supply two matrices of the same size');
        return;
    end
    
    foo = (x-y)';
    dist = foo'*conj(foo);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dist = LRdistance(b,Rp,sigma_p)
		dist = 2*b*Rp'/sigma_p - 1;
		
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%		
function plot_distances(ud)
	
	f = figure('Position',[500 300 700 700],'Units','Normalized');
	width = 0.9;
	h = 0.4;
	top = 0.55;
	Caxes1 = axes('position', [0.07 top width h]);
	hold on;
	Caxes2 = axes('position', [0.07 top-h-0.07 width h]);
	hold on;
	
	for mnum=1:3
			clear dist
		
			switch mnum
				case 1
					D = ud.DTW.Cep.D;
					p = ud.DTW.Cep.p;
					ctest = ud.audiodata_test.cepstrum;
					cref = ud.audiodata_ref.cepstrum;
					sig1 = 'k-';
					legtext = 'CC';
				case 2
					D = ud.DTW.MFCC.D;
					p = ud.DTW.MFCC.p;
					ctest = ud.audiodata_test.MFCC;
					cref = ud.audiodata_ref.MFCC;
					sig1 = 'b-';
					legtext = 'MFCC';
				case 3
					D = ud.DTW.LPC.D;
					p = ud.DTW.LPC.p;
					ctest = ud.audiodata_test.LPC;
					cref = ud.audiodata_ref.LPC;
					sig1 = 'r-';
					legtext = 'LPC';
			end

			for i=1:length(p)
				dist(i) = D(p(i,1),p(i,2));
				if strcmp(legtext,'LPC')
					local_distance(i) = LRdistance(cref.b(p(i,2),:),ctest.R(p(i,1),:),ctest.sigma_p2(p(i,1)));
				else
					local_distance(i) = distance(cref(p(i,2),:),ctest(p(i,1),:));
				end
			end
			
% 			dist(find(dist == Inf)) = 1;
% 			local_distance(find(dist == Inf)) = 1;
			
			axes(Caxes1);
			tpp(mnum) = plot(fliplr(dist)./max(dist),sig1);
			leg(mnum) = {legtext};

			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
			% Now plot local distance as a function of node along optimal path

			axes(Caxes2);
			tpld(mnum) = plot(fliplr(local_distance)./max(local_distance),sig1);
			
	end
	
	axes(Caxes1);	
	grid on;
	ylabel('Normalized Distance');
	set(Caxes1,'XTickLabel',[]);
	legend(tpp,leg,'Location','NorthWest');
	tt = title('Normalized Average Distances');
	set(tt,'FontSize',12,'FontWeight','bold');
	axis tight;
	
	axes(Caxes2);
	grid on;
	xlabel('Node');
	ylabel('Normalized Distance');
	tt = title('Normalized Local Distances');
	set(tt,'FontSize',12,'FontWeight','bold');
	axis tight;

% 	avgdist = (max(local_distance)+min(local_distance))/2;
% 	if (local_distance(end) > avgdist)
% 			text(3,min(local_distance)*1.1,'Local Distance',...
% 					'BackgroundColor',[1 1 1],'FontSize',12,'FontWeight','bold');
% 	else
% 			text(3,max(local_distance)*0.85,'Local Distance',...
% 					'BackgroundColor',[1 1 1],'FontSize',12,'FontWeight','bold');
% 	end
% 	axis([0 length(p)+1 0 max(local_distance)+0.05]);
% 	text(3,ceil(dist(1))*0.85,'Accumulated Distance',...
% 			'BackgroundColor',[1 1 1],'FontSize',12,'FontWeight','bold');
% 	axis([0 length(p)+1 0 ceil(dist(1))]);