gusucode.com > 声音的处理有：LPC，FFT，共振峰，频谱源码程序 > siganlandsystemusingMatlab/SSUM/speech/dtw/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
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))]);