clear; close all;
N_type = 5; N_slice = 1; N_data = 100;
N_clip = N_type*N_slice*N_data;

oct_diff = [0 200 400 800 1600 3200 6400 11025];  % Freq Normalized interval
oct_alpha = 0.02;  % tunable

N_freq_intvl = (length(oct_diff)-1);
feature = zeros(N_clip, N_freq_intvl*2);
style = zeros(N_clip, 1);

%

% Data Generating
data_type = ["classical", "country", "jazz", "hiphop", "metal"];
%data_diff = [0 5 10 15 20 25 30]/30;    % Time Normalized interval
data_diff = [0 30]/30;    % Time Normalized interval
for i = 1:1:N_type
    display("Now processing "+data_type(i));
    for j = 0:1:N_data-1
        num = sprintf('%02.0f', j);
        path = "Data/genres_original/"+data_type(i)+"/"+data_type(i)+".000"+num+".wav";
        % display(path);
        [y,Fs] = audioread(path);
        N_y = length(y);
        for k = 1:1:N_slice
            y_slice = y(floor(N_y*data_diff(k))+1:1:floor(N_y*data_diff(k+1))); % Slice
            Y = fft(y_slice); % FFT
            N_Y = length(Y);
            % Octave-Scale Filter
            oct_diff_norm = oct_diff/(Fs/2);  % Normalized interval
            sc = zeros(1,N_freq_intvl); valley = zeros(1,N_freq_intvl); peak = zeros(1,N_freq_intvl);
            for m = 1 : length(oct_diff_norm)-1
                [peak(m), sc(m), valley(m)] = oct_feature(oct_diff_norm(m), oct_diff_norm(m+1), Y(floor(N_Y/2):N_Y), oct_alpha);
            end
            feature((i-1)*N_data*N_slice+j*N_slice+(k-1)+1,:) = [sc, valley];
            style((i-1)*N_data*N_slice+j*N_slice+(k-1)+1,1) = i;
        end
    end
end

save('feature.mat','feature');
save('style.mat','style');

%

load('feature.mat','feature');
load('style.mat','style');

%

writematrix(feature,'feature.csv');
writematrix(style,'style.csv');

%

% clear; close all;
% load('feature.mat','feature');
% load('style.mat','style');
% X = feature(1:500,:);    
% [N,p] = size(X);
% G = cov(X,1);
% [U,L,A] = svd(G);
% lambda=diag(L)
% Mx = repmat(mean(X),N,1);
% C = (X-Mx)*A;
% min = 10000;
% for m=1:1:p
%     m
%     Xhat = C(:,1:m)*A(:,1:m)'+Mx;
%     E = X(:) - Xhat(:);
%     MSE1 = sum(E.^2)/N 
%     % mean square error
%     MSE2 = sum(lambda(m+1:p))  
%     % mean square error
%     MSE = (MSE1+MSE2)/2;
%     if (MSE<min)
%         min = MSE
%         Xhat_opt = Xhat;
%     end
% end

%

for i = 1:1:N_type
    feature_klt((N_slice*N_data*(i-1)+1):(N_slice*N_data*i),:) = ...
        klt_audio(feature((N_slice*N_data*(i-1)+1):(N_slice*N_data*i),:));
end

%

alpha_spilt = 0.2;   % training & testing data split factor
num_train = N_clip*(1-alpha_spilt);
num_test = N_clip*alpha_spilt;
% x means feature; y means label
x_train = zeros(num_train, N_freq_intvl*2);
x_test = zeros(num_test, N_freq_intvl*2);
y_train = zeros(num_train, 1);
y_test = zeros(num_test, 1);
% Random and class-uniform distribution & 
for i = 1:1:N_type
    indices=randperm(N_slice*N_data);
    x_tmp = feature_klt((N_slice*N_data*(i-1)+1):(N_slice*N_data*i),:);
    y_tmp = style((N_slice*N_data*(i-1)+1):(N_slice*N_data*i),:);
    x_tmp = x_tmp(indices,:);
    y_tmp = y_tmp(indices,:);
    x_train((N_slice*N_data*(1-alpha_spilt)*(i-1)+1):(N_slice*N_data*(1-alpha_spilt)*i),:) = ...
        x_tmp(1:N_slice*N_data*(1-alpha_spilt),:);
    x_test((N_slice*N_data*alpha_spilt*(i-1)+1):(N_slice*N_data*alpha_spilt*i),:) = ...
        x_tmp(N_slice*N_data*(1-alpha_spilt)+1:N_slice*N_data,:);
    y_train((N_slice*N_data*(1-alpha_spilt)*(i-1)+1):(N_slice*N_data*(1-alpha_spilt)*i),:) = ...
        y_tmp(1:N_slice*N_data*(1-alpha_spilt),:);
    y_test((N_slice*N_data*alpha_spilt*(i-1)+1):(N_slice*N_data*alpha_spilt*i),:) = ...
        y_tmp(N_slice*N_data*(1-alpha_spilt)+1:N_slice*N_data,:);
end

%

% Data
% x_train, y_train => training feature/label
% x_test, y_test   => testing  feature/label

% train
[y_pred_train,model,llh] = mixGaussEm(feature_klt',style');
num_correct = sum(style==y_pred_train')
acc_train = num_correct/(N_type*N_slice*N_data)*100

% train
[y_pred_train,model,llh] = mixGaussEm(x_train',y_train');
num_correct = sum(y_train==y_pred_train')
acc_train = num_correct/(N_type*N_slice*N_data*(1-alpha_spilt))*100    % training performance

% predict
y_pred_test = mixGaussPred(x_test',model);
num_correct = sum(y_test==y_pred_test')
acc_test = num_correct/(N_type*N_slice*N_data*alpha_spilt)*100