Music information retrieve for music classification

classfication_method.R

## Classification starts now
library(MASS)
library(e1071)
library(rda)
####################################
 qda.model=function(traindata){
 	qda.result=qda(Y~.,data=traindata)
 	return(qda.result)
 }
 qda.pred=function(model,dataSets.X){
 	predresult=predict(model,dataSets.X)
   # predresult=lapply(dataSets.X,function(X) predict(model,X)$class) 
	 # errorInfo(predresult,"qda")
	 return(predresult)
 }
####################################
 rda.model=function(traindata){
  rda.result=rda(Y~.,data=traindata)
 	return(rda.result)
 }
 rda.pred=function(model,dataSets.X){
   # predresult=lapply(dataSets.X,function(X) predict(model,X)$class) 
	 # errorInfo(predresult,"rda")
 	predresult=predict(model,dataSets.X)
	 return(predresult)
 }
####################################
 baye.model=function(traindata){
	 baye.result=naiveBayes(Y~.,data=traindata)
  return(baye.result)	
 }
 baye.pred=function(model,dataSets.X){
   # predresult=lapply(dataSets.X,function(X) predict(model,X)) 
 	predresult=predict(model,dataSets.X)
	 # errorInfo(predresult,"baye")
	 return(predresult)
 }

#####################################
 confusionMatrix=function(test_label,test_pred){

 }

classification.R

setwd("~/Developer/Courses/CS598PS_project/");
source('classfication_method.R');


featureArray=c(2:73)
n_pca=min(10,length(featureArray))

get_classification_model=function(featureArray=c(2:73),n_pca=10,cv=FALSE){
	cpop=read.csv("~/Developer/Courses/CS598PS_project/FeatureFiles/cpop.txt",header=F)
	opera=read.csv("~/Developer/Courses/CS598PS_project/FeatureFiles/opera.txt",header=F)

	cpop[,1]=1;
	opera[,1]=0;
	
	names(cpop)[1]="Y";
	names(opera)[1]="Y";

	num_cpop=dim(cpop)[1];
	num_opera=dim(opera)[1];

	print(num_cpop)
	print(num_opera)

	training_index_cpop=sample(c(1:num_cpop),0.5*num_cpop,replace=FALSE)
	training_index_opera=sample(c(1:num_opera),0.5*num_opera,replace=FALSE)

	testing_index_cpop=c(1:num_cpop)[-training_index_cpop];
	testing_index_opera=c(1:num_opera)[-training_index_opera];

	cpop_training_data=cpop[training_index_cpop,];
	cpop_testing_data=cpop[testing_index_cpop,];

	opera_training_data=opera[training_index_opera,];
	opera_testing_data=opera[testing_index_opera,];
	
	if(cv==TRUE){
		trainData_raw=rbind(cpop_training_data,opera_training_data);
		training_label=trainData_raw[,1]
		trainData_raw=trainData_raw[,featureArray]
		

		testData_raw=rbind(cpop_testing_data,opera_testing_data);
		testing_label=testData_raw[,1]
		testData_raw=testData_raw[,featureArray]
		

		pca_model=prcomp(trainData_raw,center=TRUE,scale=TRUE);
		trainData=as.data.frame(predict(pca_model,trainData_raw)[,c(1:n_pca)])
		testData=as.data.frame(predict(pca_model,testData_raw)[,c(1:n_pca)])
		trainData$Y=training_label

	###
		model_qda=qda.model(trainData)
		qda_result=qda.pred(model_qda,testData)

		cm=table(qda_result$class,testing_label)
		# print(1-sum(abs(as.numeric(as.character(qda_result$class))-testing_label))/length(testing_label));
		return(list(cm,1-sum(abs(as.numeric(as.character(qda_result$class))-testing_label))/length(testing_label)))
	}else{
		trainData_raw=rbind(cpop,opera);
		training_label=trainData_raw[,1]
		trainData_raw=trainData_raw[,featureArray];

		library(ggplot2)
		plot_data=cbind(as.factor(training_label),trainData_raw)
		names(plot_data)[1]="label";
		p=ggplot(data=plot_data)
		p=p+geom_point(aes(V2,V3,color=label))+xlab("MFCC_1")+ylab("MFCC_2")
		# print(p)
		pdf("mfcc.pdf",height=5,width=5)
		print(p)
		dev.off()

		pca_model=prcomp(trainData_raw,center=TRUE,scale=TRUE);
		trainData=as.data.frame(predict(pca_model,trainData_raw)[,c(1:n_pca)])
		trainData$Y=training_label
		model_qda=qda.model(trainData)
		return(list(model_qda,pca_model));
	}
}



get_feature_alongTime=function(pca_model,classification_model,fileName,featureArray=c(2:73),n_pca=10,saveName="tmp.pdf"){
	dataSet=read.csv(fileName,header=F);
	dataSet=as.data.frame(dataSet)
	dataSet_pca=as.data.frame(predict(pca_model,dataSet[,featureArray]))[1:n_pca];
	prediction_result=predict(classification_model,dataSet_pca)
	prob=prediction_result$posterior[,1]
	# c(1:length(prob)*10-5)
	# plot(prob)
	plot_data=c()
	time=c()
	for(i in 1:length(prob) ){
		plot_data=c(plot_data,prob[i],prob[i])
		time=c(time,(i-1)*10,i*10)
	}
	time_hs=strptime(time,"%s")
	time_hs= format(time_hs, format = "%M:%S")

	plot_data=data.frame(data=plot_data,time=time_hs,second=time)

	# print(plot_data)
	library(ggplot2)
	p=ggplot(data=plot_data)
	p=p+geom_line(aes(second,data))+scale_x_discrete(breaks=plot_data$second,labels=as.character(plot_data$time))
	p=p+xlab("Time")+ylab("Probability to Chinese Opera")+theme(axis.text.x  = element_text(angle=45,vjust=0.5))+ylim(c(0,1))
	pdf(saveName,width=8,height=4)
	print(p)
	dev.off()
	# print(p)
	return(plot_data)
}

accuracy=c()

cm=matrix(c(0,0,0,0),ncol=2)
for(i in 1:100){
	model_validate_result=get_classification_model(featureArray=featureArray,cv=TRUE,n_pca=n_pca)
	accuracy=c(accuracy,model_validate_result[[2]])
	cm=cm+model_validate_result[[1]]
}

print(sum(accuracy/100));
print(cm/100)

model=get_classification_model(featureArray=featureArray,cv=FALSE,n_pca=n_pca);

prob_to_opera_mix1=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/1_mix.txt",featureArray=featureArray,n_pca=n_pca,"1_mix.pdf");
prob_to_opera_mix2=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/2_mix.txt",featureArray=featureArray,n_pca=n_pca,"2_mix.pdf");
prob_to_opera_mix3=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/3_mix.txt",featureArray=featureArray,n_pca=n_pca,"3_mix.pdf");

prob_to_opera_pop1=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/pop1.txt",featureArray=featureArray,n_pca=n_pca,"pop1.pdf");
prob_to_opera_pop2=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/pop2.txt",featureArray=featureArray,n_pca=n_pca,"pop2.pdf");
prob_to_opera_pop3=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/pop3.txt",featureArray=featureArray,n_pca=n_pca,"pop3.pdf");

prob_to_opera_opera1=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/opera1.txt",featureArray=featureArray,n_pca=n_pca,"opera1.pdf");
prob_to_opera_opera2=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/opera2.txt",featureArray=featureArray,n_pca=n_pca,"opera2.pdf");
prob_to_opera_opera3=get_feature_alongTime(model[[2]],model[[1]],"./FeatureFiles/opera3.txt",featureArray=featureArray,n_pca=n_pca,"opera3.pdf");




# model_rda=rda.model(trainData)
# rda_result=rda.pred(model_qda,testData)

# model_bay=baye.model(trainData)
# bay_result=baye.pred(model_qda,testData)


## logistical regression
# logreg.result=glm(Y~.,data=traindata,family=binomial,maxit=200) 
# log_result=predict(log_result,testData)

expandFeatureStruct.m

function [structData] = expandFeatureStruct(featureStruct)
	SFname=fieldnames(featureStruct);
	structData=zeros(1,featureStruct.num);
	start=1;
	for i=2:numel(SFname)
		featureData=featureStruct.(SFname{i});
		featureSize=size(featureData);
		structData(start:start+featureSize(2)-1)=featureData;
		start=start+featureSize(2);
	end;

extractSingleSong.m

function[] = extractSingleSong(sourceFile,targetFile,filepath)

filename=[filepath sourceFile];
fileID = fopen(targetFile,'w');
[y,Fs]=audioread(filename);
% y_music=(y(:,1)-y(:,2))/2;
y_single=(y(:,1)+y(:,2))/2;

for j=10:10:floor(length(y_single)/Fs),
	Song=miraudio(y_single((j-10)*Fs+1:j*Fs),Fs,'Normal');
%     j/10
%     mirplay(Song)
 	Song_feature=featureExtract(Song);
 	data=expandFeatureStruct(Song_feature);
 	fprintf(fileID,',%10.4f',data);
 	fprintf(fileID,'\n');
end;
fclose(fileID);

extractStruct.m

function [ structData ] = extractStruct(structureObject)
	
	SFname=fieldnames(structureObject);
	% structData=zeros(1,numel(SFname)-1);
	structData=zeros(1,2); % Only extract mean and var
    for i=1:2,
		structData(i)=structureObject.(SFname{i});
	end;

featureExtract.m

function [ feature ] = featureExtract(mirAudioObject)
% This function is to extract features from a single mirObject
% 	len=get(mirAudioObject,'Length');
% 	len=len{1,1}{1,1};

% 	fs=get(mirAudioObject,'Sampling');
% 	fs=fs{1,1};
	
% 	segments_index=0:fs*10:len;
% 	segments_index=segments_index(2:end);
	

	%%%%%%%%%%%%%%%%%%%%%%%
	% Segments/Frames/bands/spectrogram/cepstrum
	%%%%%%%%%%%%%%%%%%%%%%%

	% segments=mirsegment()
	segments=mirAudioObject;
	frames=mirframe(segments,'length',2048,'sp');
		
	spec=mirspectrum(frames);
% 	spec_db=mirspectrum(frames,'dB');
	
	bands=mirfilterbank(segments,'Scheirer');  % Divide the data into Scheirer Bands
	envelop=mirenvelope(bands);
% 	onsets_curve=(envelop);  % onset_curve

	ceptrum=mircepstrum(spec,'Freq','Min',0.0012,'Max',0.02);
% 	peaks=mirpeaks(ceptrum,'Track');

	%%%%%%%%%%%%%%%%%%%%%%%
	% Extract features 
	%%%%%%%%%%%%%%%%%%%%%%%
		feature.num=0;
		% Extract whole features
		mfcc_feature=mirmfcc(frames); % MFCC --> 13*2 features
		feature.mfcc=reshape([mean(mirgetdata(mfcc_feature),2) var(mirgetdata(mfcc_feature),0,2)],[1,26]);
		feature.num=feature.num+26;

		lowenergy_feature=mirlowenergy(frames); % lowenergy ---> 1 feature
		feature.lowenergy=mirgetdata(lowenergy_feature);
		feature.num=feature.num+1;

		% Extract mean values of frequencies band feature
		% envelopData=mirgetdata(envelop);
		% envelop_size=size(envelopData);
		% envelop_feature=zeros(2,envelop_size(3)); % envelop size of different bins ---> 2 * bin number features
		
		% for i=1:envelop_size(3),
			% envelop_feature(1,i)=mean(envelopData(:,1,i));
			% envelop_feature(2,i)=var(envelopData(:,1,i));
		% end;

% 		feature.envelop=reshape(envelop_feature,[1,2*envelop_size(3)]);
% 		feature.num=feature.num+2*envelop_size(3);

		% Rhythm
% 		[tempo_feature ac] = mirtempo(mirAudioObject); % tempo
% 		onsets_feature=(mirAudioObject, ‘Diff', ‘Sum', ‘no’, ‘Filterbank', 5, 'H alfwavediff', 'Detect', 'no')


		% Timbre
		%% Timbre statistics
		% autocor_feature=mirautocor(frames); % autocorrelation
        
		zerocross_feature=mirzerocross(frames); % zerocross ---> stat (6 features)
		feature.zerocross=extractStruct(mirstat(zerocross_feature));
			
		rolloff_feature=mirrolloff(spec); % rolloff ---> stat (6 features)
		feature.rolloff=extractStruct(mirstat(rolloff_feature));

		brightness_feature=mirbrightness(spec); % brightness --->  stat (6 features)
		feature.brightness=extractStruct(mirstat(brightness_feature));
		
		rms_feature=mirrms(frames);  % rms --->  stat (6 features)
		feature.rms=extractStruct(mirstat(rms_feature));
		
		centroid_feature=mircentroid(spec); % spectrogram centroid --->  stat (6 features)
		feature.centroid=extractStruct(mirstat(centroid_feature));
		
		entropy_feature=mirentropy(spec); % entropy --->  stat (6 features)
		feature.entropy=extractStruct(mirstat(entropy_feature));


		feature.num=feature.num+6*2;


		%% Attack features
		% attacktime_feature=mirattacktime(onsets_curve);
		% attackleap_feature=mirattackleap(onsets_curve);
		
		% Extract Pitch
		% pitch=mirpitch(ceptrum,'Min',100,'Max',1000,'Total',2)
		% pitch=mirpitch(ceptrum,'Min',100,'Max',1000,'Total',2)
		
		% peaks=mirpeaks(ceptrum)
		% pitch=mirpitch(spec,'Mono','Max',1000); % fundemental frequency

		pitch_feature=mirpitch(ceptrum,'Mono');

		feature.pitch_stat=extractStruct(mirstat(pitch_feature));
		pitch_hist=hist(mirgetdata(pitch_feature),[200 300 400 500 600 700 800]);
		feature.pitch_hist=pitch_hist/sum(pitch_hist);
		feature.num=feature.num+7+2;
		% Tonality
		frames=mirframe(segments,'length',4096,'sp');
		spec_octave=mirspectrum(frames,'dB',20,'OctaveRatio',.85);
		chromagram= mirchromagram(spec_octave, 'Wrap', 'yes');
		chromagram_feature=mirgetdata(chromagram);   % 12 chromagram strength  ---> 12 features
		feature.chromagram=reshape([mean(chromagram_feature,2) var(chromagram_feature,0,2)],[1,24]);
		feature.num=feature.num+24;

featureExtraction_main.m

path_cpop='~/Developer/Project_IO/Music/cpop/';
path_opera='~/Developer/Project_IO/Music/opera/';


files=dir([path_cpop '*.mp3']);
n=length(files);
fileID = fopen('cpop.txt','w');
for i=1:n,
	filename=[path_cpop files(i).name];
	[y,Fs]=audioread(filename);
	y_music=(y(:,1)-y(:,2))/2;
	y_single=(y(:,1)+y(:,2))/2;

	% segmentIndex=[1:10:floor(length(y_single)/Fs)];
	for j=10:10:floor(length(y_single)/Fs),
		Song=miraudio(y_single((j-10)*Fs+1:j*Fs),Fs,'Normal');
		Song_feature=featureExtract(Song);
		data=expandFeatureStruct(Song_feature);
		fprintf(fileID,',%10.4f',data);
		fprintf(fileID,'\n');
	end;
	% pureMusic=miraudio(y_music(1:10*Fs),Fs);
end;
fclose(fileID);


files=dir([path_opera '*.mp3']);
n=length(files);
fileID = fopen('opera.txt','w');
for i=1:n,
	filename=[path_opera files(i).name];
	[y,Fs]=audioread(filename);
	y_music=(y(:,1)-y(:,2))/2;
	y_single=(y(:,1)+y(:,2))/2;
	
	for j=10:10:floor(length(y_single)/Fs),
		Song=miraudio(y_single((j-10)*Fs+1:j*Fs),Fs,'Normal');
		Song_feature=featureExtract(Song);
		data=expandFeatureStruct(Song_feature);
		fprintf(fileID,',%10.4f',data);
		fprintf(fileID,'\n');
	end;

end;
fclose(fileID);

% extractSingleSong(sourceFile,targetFile,filepath)

path_mix='~/Developer/Project_IO/Music/mix/';
files=dir([path_mix '*.mp3']);
n=length(files);

for i=1:n,
	extractSingleSong(files(i).name,[num2str(i) '_mix.txt'],path_mix);
end;


path_twoversion='~/Developer/Project_IO/Music/twoversion/';
extractSingleSong('opera_version_1.mp3','opera1.txt',path_twoversion)
extractSingleSong('opera_version_2.mp3','opera2.txt',path_twoversion)
extractSingleSong('opera_version_3.mp3','opera3.txt',path_twoversion)
extractSingleSong('pop_version_1.mp3','pop1.txt',path_twoversion)
extractSingleSong('pop_version_2.mp3','pop2.txt',path_twoversion)
extractSingleSong('pop_version_3.mp3','pop3.txt',path_twoversion)