close all; clear all;
load mir-st010.mat
load mir-st500.mat
%% Process a given song
songCount=length(ds);
pitchRange=[29, 83];	% Obtain from goGtCheck.m
stdScale=1:100;
stdNoteCount=length(stdScale);
transProb=zeros(stdNoteCount);
observed=cell(stdNoteCount);
for i=1:stdNoteCount
	observed{i}=[];
end
for songId=1:songCount
	fprintf('%d/%d:\n', songId, songCount);
	% Read a song
	note=ds(songId).noteGt;
	pv.pitch=ds(songId).fea.vocal_pitch;
	pv.time=ds(songId).fea.time;
	% Note to PV
	pv2=note2pv(note, 8000/256);
	% Compute transition probability
	pitch2=pv2.pitch;
	pitch2(pitch2==0)=1;	% For easy indexing
	pvCount=length(pitch2);
	for i=1:pvCount-1
		pos1=pitch2(i);
		pos2=pitch2(i+1);
		transProb(pos1, pos2)=transProb(pos1, pos2)+1;
	end
	% Compute state probablility
	for i=1:stdNoteCount
		index=find(pitch2==i);
		observed{i}=[observed{i}; pv.pitch(index)];
	end
end
%% Plot transProb
temp=transProb; temp(temp==0)=nan;
id=[1, pitchRange(1):pitchRange(2)];
temp=temp(id, id);		% Keep only the used pitch
figure; imagesc(temp); axis image; colorbar
%% Plot the length of observations for each state
observationCount=zeros(1, stdNoteCount);
for i=1:stdNoteCount
	observationCount(i)=length(observed{i});
end
figure; bar(observationCount);
xlabel('Semitones'); ylabel('Pitch count');
%%
figure; plot(observed{60}, '.-');
fprintf('Sil percentage=%g%%\n', 100*sum(observed{60}==0)/length(observed{60}));