dtwexpofn.m

非常好的数字处理教程

%
function dtwexpofn(action, datastruct)

    if nargin < 1
        action='init';
    end

  name = mfilename;
  figname = [name(1:end-2) '_fig'];
  f=findobj('Tag',figname);
  handles = get(f,'UserData');

  switch action
    case 'help'
      display_help(figname);
    case 'init'
      setdefaults;
      movegui(f,'center');
      reset(handles.timeplot_ref);
      reset(handles.timeplot_test);
      reset(handles.dtwplot);
    case 'resize'
      movegui(f,'onscreen');
    case 'load'
      switch datastruct
        case 'test'
          handles.audiodata_test = load_audiodata;
          if ~isfield(handles.audiodata_test, 'filenamepath')
              return;
          end
          if (size(handles.audiodata_test.data,2) > 1)
              handles.audiodata_test.data = to_mono(handles.audiodata_test.data);
          end
          handles.audiodata_test.data = normalize(handles.audiodata_test.data);
          if handles.audiodata_test.Fs > 8000
              handles.audiodata_test.data = ...
                  resample(handles.audiodata_test.data, 8000, ...
                  handles.audiodata_test.Fs);
              handles.audiodata_test.Fs = 8000;
					end
					% Pre-emphasize
					%handles.audiodata_test.data = filter([1 -.97], 1, handles.audiodata_test.data);
          windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata_test.Fs/1000;
          windowskip = str2num(get(handles.windowskip,'String'))* ...
              handles.audiodata_test.Fs/1000;
          % Endpoint sound
          handles.audiodata_test = endpoint(handles.audiodata_test, windowsize, windowskip);
    
          handles.audiodata_test = getStatistics(handles.audiodata_test, ...
              windowsize, windowskip);
          makeTimePlot(handles,'test');
          updateInfo(handles,'test');
              
        case 'ref'
          handles.audiodata_ref = load_audiodata;
          if ~isfield(handles.audiodata_ref, 'filenamepath')
              return;
          end
          if (size(handles.audiodata_ref.data,2) > 1)
              handles.audiodata_ref.data = to_mono(handles.audiodata_ref.data);
          end
          handles.audiodata_ref.data = normalize(handles.audiodata_ref.data);
          if handles.audiodata_ref.Fs > 8000
              handles.audiodata_ref.data = ...
                  resample(handles.audiodata_ref.data, 8000, ...
                  handles.audiodata_ref.Fs);
              handles.audiodata_ref.Fs = 8000;
					end
					% Pre-emphasize
					%handles.audiodata_ref.data = filter([1 -.97], 1, handles.audiodata_ref.data);
          windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata_ref.Fs/1000;
          windowskip = str2num(get(handles.windowskip,'String'))* ...
              handles.audiodata_ref.Fs/1000;
          % Endpoint sound
          handles.audiodata_ref = endpoint(handles.audiodata_ref, windowsize, windowskip);
          handles.audiodata_ref = getStatistics(handles.audiodata_ref, ...
              windowsize, windowskip);
          makeTimePlot(handles,'ref');
          updateInfo(handles,'ref');
      end
      drawnow;
      
      if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata_test')
          cla(handles.dtwplot);
          handles = findPath(handles,'Cepstrum');
					handles = findPath(handles,'MFCC');
          handles = findPath(handles,'LPC');
          updateDTWPlot(handles);
          drawnow;
          %handles.updateAnalysis = false;
      end

    case 'readinput'

    case 'play'
      switch datastruct
          case 'test'
              if isfield(handles,'audiodata_test')
                  audiodata = handles.audiodata_test;
                  button = handles.play_test;
                  play_audiodata(audiodata, button);
              end
          case 'ref'
              if isfield(handles,'audiodata_ref')
                  audiodata = handles.audiodata_ref;
                  button = handles.play_ref;
                  play_audiodata(audiodata, button);
              end
      end
%     case 'findPaths'
%       if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata_test')
%           cla(handles.dtwplot);
%           handles = findPath(handles);
%           updateDTWPlot(handles);
%       end
    case 'updateDTWPlot'
      if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata_test')
          cla(handles.dtwplot);
          updateDTWPlot(handles);
			end
    case 'analyze'
      if isfield(handles, 'audiodata_ref') 
        windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata_ref.Fs/1000;
        windowskip = str2num(get(handles.windowskip,'String'))* ...
            handles.audiodata_ref.Fs/1000;
				handles.audiodata_ref = getStatistics(handles.audiodata_ref, ...
              windowsize, windowskip);
        updateInfo(handles,'ref');
				
			elseif isfield(handles, 'audiodata_test')
				windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata_test.Fs/1000;
        windowskip = str2num(get(handles.windowskip,'String'))* ...
            handles.audiodata_test.Fs/1000;
				handles.audiodata_test = getStatistics(handles.audiodata_test, ...
            windowsize, windowskip);
        updateInfo(handles,'test');
				
			else
				return;
			end
			
			if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata_test')
        cla(handles.dtwplot);
        set(handles.datatype,'Value',1);
        handles = findPath(handles,'Cepstrum');
        updateDTWPlot(handles);
        drawnow;
        handles = findPath(handles,'LPC');
			end
		case 'policy'
			if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata_test')
        cla(handles.dtwplot);
        set(handles.datatype,'Value',1);
        handles = findPath(handles,'Cepstrum');
        updateDTWPlot(handles);
        drawnow;
				handles = findPath(handles,'MFCC');
        handles = findPath(handles,'LPC');
			end
		case 'show_distances'
			plot_distances(handles);

    case {'fourier', 'alias', 'feature', 'firexpo', 'iirexpo', 'formantexpo', 'lpcexpo'}
        if isfield(handles,'audiodata_test'),
            audiodata = handles.audiodata_test;
          switch action
            case 'fourier'
                fourierexpogui(audiodata);
            case 'alias'
                aliasexpogui(audiodata);
            case 'feature'
                featurexpogui(audiodata);
            case 'firexpo'
                firexpogui(audiodata);
            case 'iirexpo'
                iirexpogui(audiodata);
            case 'formantexpo'
                formantexpogui(audiodata);
            case 'lpcexpo'
              lpcexpogui(audiodata);
          end
         end
    case 'print'
      print_figure(f);
    case 'close'
      close_figure(f,figname(1:end-4));
      return;
end
set(f,'UserData',handles);


% --------------------------------------------------------------------

function makeTimePlot(ud,type)
    switch type
        case 'test'
            axes(ud.timeplot_test);
            ud.t = [0:1/ud.audiodata_test.Fs:(length(ud.audiodata_test.data)-1)/...
                ud.audiodata_test.Fs];
            tplot = plot(ud.t,ud.audiodata_test.data);
            maxtime = length(ud.t)/ud.audiodata_test.Fs;
            set(ud.timeplot_test,'XLim',[0 maxtime]);
            set(ud.timeplot_test,'YLim',[-1.0 1.0],'YTickLabel','');
            grid;
            xlabel('time (s)');
        case 'ref'
            axes(ud.timeplot_ref);
            ud.t = [0:1/ud.audiodata_ref.Fs:(length(ud.audiodata_ref.data)-1)/...
                ud.audiodata_ref.Fs];
            tplot = plot(ud.audiodata_ref.data,ud.t);
            maxtime = length(ud.t)/ud.audiodata_ref.Fs;
            set(ud.timeplot_ref,'YLim',[0 maxtime], ...
                'YAxisLocation','left');
            set(ud.timeplot_ref,'XLim',[-1.0 1.0],'XTickLabel','');
            grid;
            ylabel('time (s)');
    end
     

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function updateInfo(ud,type)
  switch type
    case 'test'
      box = ud.testinfo;
      data = ud.audiodata_test;
    case 'ref'
      box = ud.refinfo;
      data = ud.audiodata_ref;
  end
  
  windowskip = str2num(get(ud.windowskip,'String'))*data.Fs/1000;
  filename = data.filenamepath;
  [t,r] = strtok(filename,'/');
  while ~isempty(r)
      [t,r] = strtok(r,'/');
  end
  info = {...
      ['Name: ' t]; ...
      ['Fs: ' num2str(data.Fs)]; ...
      ['Length: ' num2str(size(data.data,1))]; ...
      ['Data Frames: ' num2str(ceil(size(data.data,1)/windowskip))]};
  set(box,'String', info);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ud = findPath(ud,type)
  
    switch type
      case 'Cepstrum'
        ctest = ud.audiodata_test.cepstrum;
        cref = ud.audiodata_ref.cepstrum;
      case 'MFCC'
				ctest = ud.audiodata_test.MFCC;
        cref = ud.audiodata_ref.MFCC;
      case 'LPC'
        ctest = ud.audiodata_test.LPC;
        cref = ud.audiodata_ref.LPC;
		end
    
		if strcmp(type,'LPC')
			N = size(ctest.b,1);
			M = size(cref.b,1);
		else
			N = size(ctest,1);
			M = size(cref,1);
		end
		
    D = zeros(N,M) + Inf;
    psi = zeros(N,M);
		if strcmp(type,'LPC')
			D(1,1) = LRdistance(ctest.b(1,:), cref.R(1,:), cref.sigma_p2(1));
		else
			D(1,1) = distance(ctest(1,:), cref(1,:));
		end
    psi(1,1) = 0;
    % Policies for determining paths [right in ntest dir., up in ntrain
    % dir.]

		policy_contents = get(ud.policymenu,'String');
		switch policy_contents{get(ud.policymenu,'Value')}
			case 'I'
				policy{1} = [1 0];			% Path I
				policy{2} = [1 1];      % Path II
				policy{3} = [0 1];      % Path III
				psi(1,1:M) = 3;
				psi(1:N,1) = 1;
				Qmax = inf;
				Qmin = 0;
				for n=2:N
					if strcmp(type,'LPC')
						D(n,1) = LRdistance(ctest.b(n,:), cref.R(1,:), cref.sigma_p2(1));
					else
						D(n,1) = distance(ctest(n,:), cref(1,:));
					end
				end
				for m=2:M
					if strcmp(type,'LPC')
						D(1,m) = LRdistance(ctest.b(1,:), cref.R(m,:), cref.sigma_p2(m));
					else
						D(1,m) = distance(ctest(1,:), cref(m,:));
					end
				end
				for n=2:N
					%mmax = min(Qmax*(n-1)+1,ceil(Qmin*(n-N)+M));
					%mmin = max(ceil(Qmin*(n-1)+1),Qmax*(n-N)+M);
					mmax = M;
					mmin = 2;
					% [mmin mmax]
					% Find distance along allowed paths
					for m = mmin:mmax
							d = ones(1,length(policy))*Inf;
							if strcmp(type,'LPC')
								local_distance = LRdistance(ctest.b(n,:), cref.R(m,:), cref.sigma_p2(m));
							else
								local_distance = distance(cref(m,:),ctest(n,:));
							end

							d(1) = D(n-1,m) + local_distance;
							d(2) = D(n-1,m-1) + 2*local_distance;
							d(3) = D(n,m-1) + local_distance;
							
							D(n,m) = min(d);

							% Now path of origin for minimum distance path
							psi(n,m) = find(d==min(d));

					end
				end

			case 'II'
				policy{1} = [1 1; 1 0]; % Path I
				policy{2} = [1 1];      % Path II
				policy{3} = [1 1; 0 1];      % Path III
				Qmax = 2;
				Qmin = 0.5;
				
				for n=2:N
					mmax = min(Qmax*(n-1)+1,ceil(Qmin*(n-N)+M));
					mmin = max(ceil(Qmin*(n-1)+1),Qmax*(n-N)+M);
					% [mmin mmax]
					% Find distance along allowed paths
					for m = mmin:mmax
							d = ones(1,length(policy))*Inf;
							if strcmp(type,'LPC')
								local_distance = LRdistance(ctest.b(n,:), cref.R(m,:), cref.sigma_p2(m));
							else
								local_distance = distance(cref(m,:),ctest(n,:));
							end
							
							if n > 2,
								if strcmp(type,'LPC')
									d(1) = D(n-2,m-1) + 0.5*LRdistance(ctest.b(n-1,:), cref.R(m,:), cref.sigma_p2(m)) + 0.5*local_distance;
								else
									d(1) = D(n-2,m-1) + 0.5*distance(cref(m,:),ctest(n-1,:)) + 0.5*local_distance;
								end
							end
							
							d(2) = D(n-1,m-1) + local_distance;
							
							if m > 2
								if strcmp(type,'LPC')
									d(3) = D(n-1,m-2) + 0.5*LRdistance(ctest.b(n,:), cref.R(m-1,:), cref.sigma_p2(m-1)) + 0.5*local_distance;
								else
									d(3) = D(n-1,m-2) + 0.5*distance(cref(m-1,:),ctest(n,:)) + 0.5*local_distance;
								end
							end

							D(n,m) = min(d);

							% Now path of origin for minimum distance path
							psi(n,m) = find(d==min(d));

					end
				end
			case 'III'
				policy{1} = [2 1]; % Path I
				policy{2} = [1 1]; % Path II
				policy{3} = [1 2]; % Path III
				Qmax = 2;
				Qmin = 0.5;
				
				for n=2:N
					mmax = min(Qmax*(n-1)+1,ceil(Qmin*(n-N)+M));
					mmin = max(ceil(Qmin*(n-1)+1),Qmax*(n-N)+M);
					% [mmin mmax]
					% Find distance along allowed paths
					for m = mmin:mmax
							d = ones(1,length(policy))*Inf;
							if strcmp(type,'LPC')
								local_distance = LRdistance(ctest.b(n,:), cref.R(m,:), cref.sigma_p2(m));
							else
								local_distance = distance(cref(m,:),ctest(n,:));
							end
							d(2) = D(n-1,m-1) + 2*local_distance;
							if n > 2,
									d(1) = D(n-2,m-1) + 3*local_distance;
							end

							if m > 2,
									d(3) = D(n-1,m-2) + 3*local_distance;
							end

							D(n,m) = min(d);

							% Now path of origin for minimum distance path
							psi(n,m) = find(d==min(d));

					end
				end
				
			case 'IV'
				policy{1} = [1 1; 1 0]; % Path I
				policy{2} = [1 2; 1 0]; % Path II
				policy{3} = [1 1];      % Path III
				policy{4} = [1 2];      % Path IV
				Qmax = 2;
				Qmin = 0.5;

				for n=2:N
					mmax = min(Qmax*(n-1)+1,ceil(Qmin*(n-N)+M));
					mmin = max(ceil(Qmin*(n-1)+1),Qmax*(n-N)+M);
					% [mmin mmax]
					% Find distance along allowed paths

					for m = mmin:mmax
							d = ones(1,length(policy))*Inf;

							if strcmp(type,'LPC')
								local_distance = LRdistance(ctest.b(n,:), cref.R(m,:), cref.sigma_p2(m));
							else
								local_distance = distance(cref(m,:),ctest(n,:));
							end
							if n > 2,
									if strcmp(type,'LPC')
										pre_local_distance = LRdistance(ctest.b(n-1,:), cref.R(m,:), cref.sigma_p2(m));
									else
										pre_local_distance = distance(cref(m,:), ctest(n-1,:));
									end
									d(1) = D(n-2,m-1) + pre_local_distance + local_distance;
									if m > 2,
											d(2) = D(n-2,m-2) + pre_local_distance + local_distance;
									end
							end

							if m > 1,
									d(3) = D(n-1,m-1) + local_distance;
									if m > 2,
											d(4) = D(n-1,m-2) + local_distance;
									end
							end
							D(n,m) = min(d);

							% Now path of origin for minimum distance path
							psi(n,m) = find(d==min(d));

					end
				end
		
		end

    % Global pattern dissimilarity
    dXY = D(N,M)./N;

    % Now find alignment path backwards
    p(1,:) = [N M];   % Start from end