vqexpofn.m

非常好的数字处理教程

					end
							
				else
					clear tempbook;
					tempbook = codebook;

					% Step 2: split centroid vectors into two codewords each
					for j=1:2^(k-1)
							codebook(2*j-1,:) = tempbook(j,:).*(1-e);
							codebook(2*j,:) = tempbook(j,:).*(1+e);
					end

					% Step 3: map each vector to these codewords
					matches = zeros(1,numvec);

					testortion(1:numvec) = 1.e15;
					for l=1:2^(k-1)
						for j=1:numvec
							if mnum == 3
								[a,b] = parcor_to_lpc(codebook(l,:));
								d = 2*b*y(j,:)' - 1;
							else
								d = sum((y(j,:)-codebook(l,:)).^2);
							end
							if d < testortion(j)
								testortion(j) = d;
								matches(j) = l;
							end
						end
					end

					dist(2) = 1.e15;
					count = 1;
					while (abs(diff(dist)) > distortion_threshold & count < MAX_ITERATIONS)
						% Step 4: Now find centroids of these matches and match
						dist(2) = 0;
						for j=1:2^k
								ind{j} = find(matches == j);

								if ~isempty(ind{j})
									centroid = zeros(1,numcc);
									for m=1:numcc
										 centroid(m) = sum(y(ind{j},m))/numel(ind{j});
									end

									% Find new distortion for set of quantized vectors
									if mnum == 3
										% Convert centroid R to parcor
										[a,b] = parcor_to_lpc(codebook(j,:));
										for m=1:numel(ind{j})
											dist(2) = dist(2) + 2*b*y(ind{j}(m),:)' - 1;
										end
										codebook(j,:) = R_to_parcor(centroid);
									else
										for m=1:numel(ind{j})
												dist(2) = dist(2) + ...
														sum((y(ind{j}(m),:)-codebook(j,:)).^2);
										end
										codebook(j,:) = centroid;
									end

								end

						end
						dist(2) = dist(2)/numvec;

						% map each cc vector to these codewords
						matches = zeros(1,numvec);

						testortion(1:numvec) = 1.e15;
						for l=1:2^(k-1)
							for j=1:numvec
								if mnum == 3
									[a,b] = parcor_to_lpc(codebook(l,:));
									d = 2*b*y(j,:)' - 1;
								else
									d = sum((y(j,:)-codebook(l,:)).^2);
								end
								if d < testortion(j)
									testortion(j) = d;
									matches(j) = l;
								end
							end
						end

						dist(1) = sum(testortion)/numvec;

						count = count + 1;
					end

					% Find final distortion
					dist(2) = 0;
					for j=1:2^k
							ind{j} = find(matches == j);

							if ~isempty(ind{j})
								% Find new distortion for set of quantized vectors
								if mnum == 3
									[a,b] = parcor_to_lpc(codebook(j,:));
									for m=1:numel(ind{j})
											dist(2) = dist(2) + 2*b*y(ind{j}(m),:)' - 1;
									end
								else
									for m=1:numel(ind{j})
											dist(2) = dist(2) + sum((y(ind{j}(m),:)-codebook(j,:)).^2);
									end
								end
							end

					end
					dist(2) = dist(2)/numvec;
					
					tempbook = [];
					% Get rid of unused codewords
					for j=1:2^k
						ind{j} = find(matches == j);
						if ~isempty(ind{j})
								tempbook = [tempbook; codebook(j,:)];
						end
					end

				end
				
				if count >= MAX_ITERATIONS
					disp('Maximum iterations exceeded!');
				end
				% Place codebook into array
				disp(['Distortion = ' num2str(dist(2))]);
				contents{k} = num2str(2^(k-1));
				if num_codewords > numvec
					disp(['Codebook size exceeds number of vectors (' num2str(numvec) ')!']);
				end

        switch mnum
          case 1
            ud.audiodata.cepstrum.codebook.codewords{k} = tempbook;
            ud.audiodata.cepstrum.codebook.distortion(k) = dist(2);
          case 2
            ud.audiodata.MFCC.codebook.codewords{k} = tempbook;
            ud.audiodata.MFCC.codebook.distortion(k) = dist(2);
					case 3
						ud.audiodata.LPC.codebook.codewords{k} = tempbook;
            ud.audiodata.LPC.codebook.distortion(k) = dist(2);
				end
				%tempbook
        k = k+1;
				waitbar((mnum-1)/3,h,['Creating ' mssg]);
    end

  end
  set(ud.codebooksize,'String',contents);
	waitbar(1,h);
	close(h);


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

    dBthreshold = 15;
    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
    if higher2 > numel(tz)
			higher2 = numel(tz);
		end
		if lower2 == 0
			lower2 = 1;
		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)
    audiodata.cepstrum.cc = cepstrum(audiodata, window_size, window_skip);
    audiodata.MFCC.cc = 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);
			R = R(ceil(end/2):end);
			R = R(1:p+1);
			Rt = toeplitz(R(1:end-1));
			lpcdata.a(i,:) = [ -1; inv(Rt)*R(2:end)]';
			%a = lpc(signal,p)';

			% For likelihood-ratio distance calculation
			lpcdata.R(i,:) = R';
			b = xcorr(lpcdata.a(i,:));
			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 k = R_to_parcor(R)
		p = 16;
		E=zeros(1,p);
    k=zeros(1,p);
    alpha=zeros(p,p);
    
    E(1)=R(1);
    ind=1;
    k(ind)=R(ind+1)/E(ind);
    alpha(ind,ind)=k(ind);
    E(ind+1)=(1-k(ind).^2)*E(ind);
    for ind=2:p
        k(ind)=(R(ind+1)-sum(alpha(1:ind-1,ind-1)'.*R(ind:-1:2)))/E(ind);
        alpha(ind,ind)=k(ind);
        for jnd=1:ind-1
            alpha(jnd,ind)=alpha(jnd,ind-1)-k(ind)*alpha(ind-jnd,ind-1);
        end
        E(ind+1)=(1-k(ind).^2)*E(ind);
		end

    %G=sqrt(E(p+1));

function [a,b] = parcor_to_lpc(k)
	p = 16;
	alpha(1:p,1:p)=0;

	for i=1:p
			alpha(i,i)=k(i);
			if (i > 1)
					for j=1:i-1
							alpha(j,i)=alpha(j,i-1)-k(i)*alpha(i-j,i-1);
					end
			end
	end
	a = [-1; alpha(1:p,p)];
	b = xcorr(a);
	b = b(ceil(end/2):end)';
	b = [b(1)/2 b(2:p+1)];
	
		
function b = R_to_b(R)
	p = 16;
	Rt = toeplitz(R(1:end-1));
	a = [ -1; inv(Rt)*R(2:end)']';

	b = xcorr(a);
	b = b(ceil(end/2):end);
	b = [b(1)/2 b(2:end)];