dtwplot.m

来自「一个关于数据聚类和模式识别的程序,在生物化学,化学中因该都可以用到.希望对大家有」· M 代码 · 共 167 行

function dtwplot(vec1, vec2, DTWpath, DTWpath2) 
%DTWPLOT Plot the result of DTW of two pitch/MFCC vectors
%	Usage:
%	dtwplot(vec1, vec2, DTWpath)
%
%	For example:
%		vec1=[71 73 75 80 80 80 78 76 75 73 71 71 71 73 75 76 76 68 76 76 75 73 71 70 70 69 68 68 72 74 78 79 80 80 78];
%		vec2=[69 69 73 75 79 80 79 78 76 73 72 71 70 70 69 69 69 71 73 75 76 76 76 76 76 75 73 71 70 70 71 73 75 80 80 80 78];
%		[minDist, dtwPath, dtwTable] = dtw1(vec1, vec2);
%		dtwplot(vec1, vec2, dtwPath);

%	Roger Jang, 20020428, 20070522

if nargin<1; selfdemo; return; end

if size(vec1,1)>1 & size(vec1,2)>1 & size(vec2,1)>1 & size(vec2,2)>1	% Inputs are MFCC matrices for ASR
	frame1num = size(vec1, 2);
	frame2num = size(vec2, 2);

	% Vertical yellow grid line
	for i = 1:frame1num,
		line([i, i], [1, frame2num], 'color', 'y');
	end
	% Horizontal yellow grid line
	for j = 1:frame2num,
		line([1, frame1num], [j, j], 'color', 'y');
	end

	% ====== Solution obtained by DTW
	for i = 1:size(DTWpath,2)-1,
		line([DTWpath(1,i), DTWpath(1,i+1)], [DTWpath(2,i), DTWpath(2,i+1)], 'color', 'r', 'marker', '.');
	end
	if nargin==4,	% Another DTW path via different DTW method
		for i = 1:size(DTWpath2,2)-1,
			line([DTWpath2(1,i), DTWpath2(1,i+1)], [DTWpath2(2,i), DTWpath2(2,i+1)], 'color', 'k', 'marker', '.');
		end
	end
	xlabel(inputname(1));
	ylabel(inputname(2));
	if ~isempty(DTWpath),
		temp = vec1(:,DTWpath(1,:))-vec2(:,DTWpath(2,:));
		dtwDist = sum(sqrt(sum(temp.^2)));
	else
		dtwDist = inf;
	end
	title(['DTW total distance = ', num2str(dtwDist)]);
	box on
	axis tight
	axis image
	return
end

vec1=vec1(:);
vec2=vec2(:);
pos = get(gcf, 'pos');
width = pos(3);
height = pos(4);
% Assume width > height
% Assume the aspect ratio of the signal axis is b:1
b = 4;
k = height/(b+2.5);

vec1Axis = subplot(2,2,4);
plot(1:length(vec1), vec1, 1:length(vec1), vec1, 'r.');
axis tight
ylabel('vec1');
set(vec1Axis, 'unit', 'pixels');
set(vec1Axis, 'position', [width-(b+0.5)*k, height-(b+2)*k, b*k, k]); 

vec2Axis = subplot(2,2,1);
plot(min(vec2)+max(vec2)-vec2, 1:length(vec2), min(vec2)+max(vec2)-vec2, 1:length(vec2), 'r.');
axis tight;
xlabel('vec2');
set(vec2Axis, 'unit', 'pixels');
set(vec2Axis, 'position', [width-(b+2)*k, height-(b+0.5)*k, k, b*k]); 
% The following should be done after axis is resized
xticklabel = get(vec2Axis, 'xticklabel');
set(vec2Axis, 'xticklabel', flipud(xticklabel));

pathAxis = subplot(2,2,2);
box on;
set(pathAxis, 'unit', 'pixels');
position = [width-(b+0.5)*k, height-(b+0.5)*k, b*k, b*k]; 
set(pathAxis, 'position', position); 
% Plot the yellow grid line
for i = 1:length(vec1),
	line([i, i], [1, length(vec2)], 'color', 'y');
end
for i = 1:length(vec2),
	line([1, length(vec1)], [i, i], 'color', 'y');
end
% ====== Solution obtained by DTW
for i = 1:size(DTWpath,2)-1,
	line(DTWpath(1, i:i+1), DTWpath(2, i:i+1), 'color', 'r', 'marker', '.');
end
if ~isempty(DTWpath),
	dtwDist = sum(abs(vec1(DTWpath(1,:))-vec2(DTWpath(2,:))));
else
	dtwDist = inf;
end
title(sprintf('DTW total distance = %g', dtwDist));

if nargin==4,	% Another DTW path via different DTW method
	for i = 1:size(DTWpath2,2)-1,
		line(DTWpath2(1, i:i+1), DTWpath2(2, i:i+1), 'color', 'k', 'marker', '.');
	end
	if ~isempty(DTWpath2),
		dtwDist2 = sum(abs(vec1(DTWpath2(1,:))-vec2(DTWpath2(2,:))));
	else
		dtwDist2 = inf;
	end
	title(sprintf('DTW total distances = %g, %g', dtwDist, dtwDist2));
end

axis equal
axis tight

if length(vec2)>length(vec1),
	% Adjust the axis of vec1
	dtwPlotWidth = position(3)*length(vec1)/length(vec2);
	oldPos = get(vec1Axis, 'pos');
	position = oldPos;
	position(1) = oldPos(1)+(oldPos(3)-dtwPlotWidth)/2;
	position(3) = dtwPlotWidth;
	set(vec1Axis, 'pos', position);
else
	% Adjust the axis of vec2
	dtwPlotHeight = position(3)*length(vec2)/length(vec1);
	oldPos = get(vec2Axis, 'pos');
	position = oldPos;
	position(2) = oldPos(2)+(oldPos(4)-dtwPlotHeight)/2;
	position(4) = dtwPlotHeight;
	set(vec2Axis, 'pos', position);
end

set(findobj(gcf, 'unit', 'pixel'), 'unit', 'normalize');

% ====== self demo ======
function selfdemo

addpath ../audioProcessing -end
file1 = 'waveData/singer1.wav';
file2 = 'waveData/singer2.wav';
y1 = wavread(file1);
y2 = wavread(file2);
mfcc1 = wave2mfccmex(y1);
mfcc2 = wave2mfccmex(y2);
rmpath ../audioProcessing

% ====== Invoke DTW
[minDist, dtwPath1] = dtw1mex(mfcc1, mfcc2);
fprintf('dtw1mex(): total distance = %g, path length = %g\n', minDist, size(dtwPath1,2));
[minDist, dtwPath2] = dtw1meanDistMex(mfcc1, mfcc2);
fprintf('dtw1meanDistMex(): average distance = %g, path length = %g\n', minDist, size(dtwPath2,2));

% ====== Plot the DTW paths
figure;
feval(mfilename, mfcc1, mfcc2, dtwPath1, dtwPath2);	% Plot two paths
title('DTW paths (red: min. total distance, green: min. average distance)');
xlabel(file1); ylabel(file2);

v1=randn(1, 25);
v2=randn(1, 30);
[minDist, dtwPath]=dtw1mex(v1, v2);
figure;
feval(mfilename, v1, v2, dtwPath);