vqexpofn.m

非常好的数字处理教程

%
function vqexpofn(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.vectorplot);
      set(handles.info,'Enable','inactive');
			set(handles.show_errors,'Visible','off');
		case 'resize'
      movegui(f,'onscreen');
    case 'load'
					set(handles.show_codebook,'Value',0);
          handles.audiodata = load_audiodata;
          if ~isfield(handles.audiodata, 'filenamepath')
              return;
          end
          if (size(handles.audiodata.data,2) > 1)
              handles.audiodata.data = to_mono(handles.audiodata.data);
          end
          handles.audiodata.data = normalize(handles.audiodata.data);
          if handles.audiodata.Fs > 8000
              handles.audiodata.data = ...
                  resample(handles.audiodata.data, 8000, ...
                  handles.audiodata.Fs);
              handles.audiodata.Fs = 8000;
					end
					% Pre-emphasize
					handles.audiodata.data = filter([1 -.9], 1, handles.audiodata.data);
          windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata.Fs/1000;
          windowskip = str2num(get(handles.windowskip,'String'))* ...
              handles.audiodata.Fs/1000;
          handles.audiodata = endpoint(handles.audiodata, windowsize, windowskip);
          handles.audiodata = getStatistics(handles.audiodata, ...
              windowsize, windowskip);
          cla(handles.vectorplot);
          updateVectorPlot(handles);
          updateInfo(handles);
          drawnow;
          handles = createCodebooks(handles);
					set(handles.show_errors,'Visible','on');
    case 'readinput'

    case 'play'
        if isfield(handles,'audiodata')
            audiodata = handles.audiodata;
            button = handles.play;
            play_audiodata(audiodata, button);
        end
    case 'updatePlot'
      if isfield(handles, 'audiodata')
          cla(handles.vectorplot);
          drawnow;
          updateVectorPlot(handles);
          updateInfo(handles);
      end
    case 'analyze'
			if isfield(handles, 'audiodata')
				set(handles.show_errors,'Visible','off');
				cla(handles.vectorplot);
				windowsize = str2num(get(handles.windowsize,'String'))* ...
						handles.audiodata.Fs/1000;
				windowskip = str2num(get(handles.windowskip,'String'))* ...
						handles.audiodata.Fs/1000;
				handles.audiodata = getStatistics(handles.audiodata, ...
						windowsize, windowskip);
				updateVectorPlot(handles);
				updateInfo(handles);
				drawnow;
				handles = createCodebooks(handles);
				set(handles.show_errors,'Visible','on');
			end
		case 'show_errors'
			f = figure('Position',[500 300 700 700],'Units','Normalized');
			width = 0.9;
			h = 0.89;
			top = 0.065;

			Caxes1 = axes('position', [0.07 top width h]);
			hold on;
				
			for mnum=1:3
				switch mnum
					case 1
						y = handles.audiodata.cepstrum.codebook.distortion;
						sig1 = 'ko--';
						legtext = 'CC';
						l = 2.^[0:size(y,2)-1];
					case 2
						y = handles.audiodata.MFCC.codebook.distortion;
						sig1 = 'bo--';
						legtext = 'MFCC';
						l = 2.^[0:size(y,2)-1];
					case 3
						y = handles.audiodata.LPC.codebook.distortion;
						sig1 = 'ro--';
						legtext = 'LPC';
						l = 2.^[0:size(y,2)-1];
				end
				
				p(mnum) = semilogx(l,y./max(y),sig1);
			end
			
			ylabel('Normalized Total Average Distortion');
			xlabel('Codebook Length');
			set(Caxes1,'XScale','log','XLim',[l(1) l(end)]);
			set(Caxes1,'XTick',l,'XTickLabel',l,'XMinorTick','off');
			axis([l(1) l(end) 0 1.05]);
			legend(p,'CC','MFCC','LPC','Location','NorthWest');
			tt = title('Distortion Values of Codebook Sizes');
			set(tt,'FontSize',14);
			grid on;
			set(Caxes1,'XMinorGrid','off');

    case 'analyze'
      if isfield(handles, 'audiodata_ref') & isfield(handles, 'audiodata')
        
        windowsize = str2num(get(handles.windowsize,'String'))* ...
              handles.audiodata.Fs/1000;
        windowskip = str2num(get(handles.windowskip,'String'))* ...
            handles.audiodata.Fs/1000;
          
        % test
        handles.audiodata = getStatistics(handles.audiodata, ...
            windowsize, windowskip);
        updateInfo(handles,'test');
        
        handles.audiodata_ref = getStatistics(handles.audiodata_ref, ...
              windowsize, windowskip);
        updateInfo(handles,'ref');
        
        cla(handles.dtwplot);
        set(handles.datatype,'Value',1);
        handles = findPath(handles,'Cepstrum');
        updateDTWPlot(handles);
        drawnow;
        handles = findPath(handles,'LPC');
      end
    case {'fourier', 'alias', 'feature', 'firexpo', 'iirexpo', 'formantexpo', 'lpcexpo'}
        if isfield(handles,'audiodata'),
            audiodata = handles.audiodata;
          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 updateInfo(ud)
    box = ud.info;
    data = ud.audiodata;

    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(floor(size(data.data,1)/windowskip))]};
    set(box,'String', info);
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function updateVectorPlot(ud)
   
    axes(ud.vectorplot);
    hold on;
    contents = get(ud.datatype,'String');
    cbcontents = get(ud.codebooksize,'String');
    distortion = 0;
		colors=['b' 'g' 'r' 'c' 'm' 'y'];
		
    switch contents{get(ud.datatype,'Value')}
      case 'Cepstrum'
        data = ud.audiodata.cepstrum.cc;
        codebooksize =  log2(str2num(cbcontents{get(ud.codebooksize,'Value')}))+1;
        if isfield(ud.audiodata.cepstrum,'codebook');
          codebook = ud.audiodata.cepstrum.codebook.codewords{codebooksize};
					distortion = ud.audiodata.cepstrum.codebook.distortion(codebooksize);
        end
      case 'MFCC'
        data = ud.audiodata.MFCC.cc;
        codebooksize =  log2(str2num(cbcontents{get(ud.codebooksize,'Value')}))+1;
        if isfield(ud.audiodata.MFCC,'codebook');
          codebook = ud.audiodata.MFCC.codebook.codewords{codebooksize};
					distortion = ud.audiodata.MFCC.codebook.distortion(codebooksize);
        end
      case 'LPC'
        data = ud.audiodata.LPC.a;
        codebooksize =  log2(str2num(cbcontents{get(ud.codebooksize,'Value')}))+1;
        if isfield(ud.audiodata.LPC,'codebook');
          codebook = ud.audiodata.LPC.codebook.codewords{codebooksize};
					distortion = ud.audiodata.LPC.codebook.distortion(codebooksize);
				end
		end
    l = [0:size(data,2)-1];
		
		switch contents{get(ud.datatype,'Value')}
			case 'LPC'
				if get(ud.original_vectors,'Value') 
					if get(ud.show_codebook,'Value')
						for i=1:size(data,1)
							[ha,f]=freqz(1,data(i,:),512,ud.audiodata.Fs);
							plot(f,20*log10(abs(ha)./max(abs(ha))),'k');
						end
					else
						for i=1:size(data,1)
							[ha,f]=freqz(1,data(i,:),512,ud.audiodata.Fs);
							plot(f,20*log10(abs(ha)./max(abs(ha))),colors(mod(i,length(colors))+1));
						end
					end
				end
				
				if get(ud.show_codebook,'Value')
					for i=1:size(codebook,1)
						[a,b] = parcor_to_lpc(codebook(i,:));
						[ha,f]=freqz(1,a,512,ud.audiodata.Fs);
						p = plot(f,20*log10(abs(ha)./max(abs(ha))), ...
							colors(mod(i,length(colors))+1),'LineWidth',2);
					end
					tt = text(0.01*ud.audiodata.Fs,-47, ...
						['Distortion: ' num2str(distortion)], 'BackgroundColor','white');
					set(tt,'FontSize',12,'FontWeight','bold');
				end
				
				ylabel('Magnitude (dB)');
				axis([0 ud.audiodata.Fs/2 -50 4.8]);
				grid on;

			otherwise
				if get(ud.original_vectors,'Value')
					if get(ud.show_codebook,'Value')
						for i=1:size(data,1)
							plot(l,data(i,:),'k');
						end
					else
						for i=1:size(data,1)
							plot(l,data(i,:),colors(mod(i,length(colors))+1));
						end
					end
				end

				if get(ud.show_codebook,'Value')
					for i=1:size(codebook,1)
						p = plot(l,codebook(i,:),colors(mod(i,length(colors))+1), ...
							'LineWidth',2);
					end
					tt = text(0.5*size(data,2),0.9*max(max(data)), ...
						['Distortion: ' num2str(distortion)], 'BackgroundColor','white');
					set(tt,'FontSize',12,'FontWeight','bold');
				end

				ylabel('Magnitude');
				axis([0 size(data,2)-1 min(min(data)) max(max(data))]);
				grid on;
		end
		
		switch contents{get(ud.datatype,'Value')}
			case 'Cepstrum'
				xlabel('Cepstral Coefficient');
			case 'MFCC'
				xlabel('Mel-Frequency Cepstral Coefficient');
			case 'LPC'
				xlabel('Frequency (Hz)');
		end

     
function ud = createCodebooks(ud)
  
  e = str2num(get(ud.epsilon,'String'));
  distortion_threshold = str2num(get(ud.error_thresh,'String'));
  %distortion_threshold = 200;
	
  % First for cepstrum
  h = waitbar(0,'Creating Codebooks');
	
  for mnum=1:3
    clear y;
    switch mnum
      case 1
        y = ud.audiodata.cepstrum.cc;
        disp('Creating codebooks for cepstrum vectors.');
				mssg = 'Cepstral Codebook';
      case 2
        y = ud.audiodata.MFCC.cc;
        disp('Creating codebooks for mel-frequency cepstrum vectors.');
				mssg = 'MF Cepstral Codebook';
			case 3
				R = ud.audiodata.LPC.R;
				%a = ud.audiodata.LPC.a;
				%b = ud.audiodata.LPC.b;
				sigma_p2 = ud.audiodata.LPC.sigma_p2;
				% Normalize autocorrelations
				for i=1:size(R,1)
					y(i,:) = R(i,:)./sigma_p2(i);
				end
        disp('Creating codebooks for LPC vectors.');
				mssg = 'LPC Codebook';
    end
    
    [numvec, numcc] = size(y);
    
    % Step 1: determine centroid vector
    centroid = zeros(1,numcc);
    for j=1:numcc
        centroid(j) = sum(y(:,j))/numvec;
    end
    codebook = centroid; % 1 element codebook
    dist = [0 0];
		
		% Find initial distortion
		if mnum == 3 % LPC
			% Convert R-centroid to b for distance
			b = R_to_b(centroid);
			for j=1:numvec
					dist(1) = dist(1) + 2*b*y(j,:)' - 1;
			end
			codebook = R_to_parcor(centroid);
		else
			for j=1:numvec
					dist(1) = dist(1) + sum((y(j,:)-centroid).^2);
			end
		end
    dist(1) = dist(1)/numvec;
    disp(['Initial distortion = ' num2str(dist(1))]);
    
    % for each codebook size, using the previous designed codebook
    MAX_ITERATIONS = 100;
    k = 1;

    while (2^(k-1) < 256)
        num_codewords = 2^(k-1);      % Must be 2^i, where i E Z

        disp(['Determining ' num2str(num_codewords) ' codewords...']);
				
				if num_codewords >= numvec
					dist(2) = 0;
					switch mnum 
						case 3
							for i=1:numvec
								tempbook(i,:) = R_to_parcor(y(i,:));
							end
						otherwise
							tempbook = y;