clear; close all; clc;

[a, fs] = audioread('music.wav');
segmentation = 5;


% ===================================================================
% Main Algorithm of Gabor Transform
% ===================================================================

x = a(:,1);
tau = 0 : 1/44100 : 3.2;
dt = 0.01;
df = 1;
t = 0 : dt : max(tau);
f = 20 : df : 2500;
sigma = 200;

dt = t(2)-t(1);
df = f(2)-f(1);
dtau = tau(2)-tau(1);
S = dt/dtau;

C = length(t);
F = length(f);
T = length(tau);

N = 1 / (df*dtau);
B = 1.9143 / (sigma^(1/2));
Q = round( B/dtau );

n0 = tau(1) / dtau;
c0 = t(1) / dt;
m0 = f(1) / df;

X = zeros(C,F);
x1 = zeros(1,N);
window = (sigma^(1/4)) * exp( -sigma*pi*( (Q-( 0:N-1 ))*dtau ).^2 );

for n = c0 :( c0+C-1 )
    window_const = dtau * exp( (-1j*2*pi*( n*S-Q ) ) .* ( m0:(m0+F-1) )./N);
    for q = 0 : N-1
        if ( (q<=(2*Q) ) && ( n*S-Q+q>=0 ) && ( n*S-Q+q+1<=T ) )
            x1(q+1) = x( n*S-Q+q+1 ); 
        else
            x1(q+1) = 0;
        end
    end    
    X1 = fft( (window.*x1), N);   
    X( (n-c0)+1, :) = window_const .* X1( m0+1 : ( (m0+F-1)+1 ) );
end

y = X'; 


% ===================================================================
% Convolution to g[n,m] and find Instantaneous Frequency
% ===================================================================

f0 = 20;
L = (2*(f0/dt));
for n_ = 1:5 % n_ = n+3
    for m_ = 1: (2*L+1)  % m_ = L+(L+1) 
        g(n_,m_) = exp(-(((m_-(L+1))*df)^2))/f0^2;
    end
end
Xs = conv2(abs(y),g','same'); 

% Create Threshold matrix of instantaneous frequency
threshold_value = max(max(abs(Xs))) * 0.01;
threshold = zeros(F,C);
threshold = threshold + threshold_value;

% Create Instantaneous Frequency zero matrix for later storing data  
InstantaneousFreq = zeros( size(Xs) );

for m = 1+1 : F-1 % to avoid reaching invalid index 
    for n = 1 : C  
        % find max frequency to become instantaneous frequency
        if ( Xs(m,n)>Xs(m-1,n) ) && ( Xs(m,n)>Xs(m+1,n) ) && ( abs((Xs(m,n))) > threshold(m,n) )
            InstantaneousFreq(m,n) = 1;
        else 
            InstantaneousFreq(m,n) = 0;
        end
    end
end


% ===================================================================
% Connect [m,n] to create clear Instantaneous Frequency line
% ===================================================================

% create Clear Frequency data matrix for later operation
ClearFrequency = []; 

for i = 1 : F
    % extract every frequency line row by row 
    tmp = InstantaneousFreq(i,:); 
    
    % identify the index of clear frequency line
    index = find( abs( real(tmp) ) > 0.0001 );
    
    for j = 1 : C
        if ismember(j, index)
            % let the clear frequency line element to 1
            tmp(j) = 1;
        else
            % remove the unclear frequemcy line
            tmp(j) = 0;
        end
    end
    
    % label the different region of clear frequemcy line
    tmp = bwlabel(tmp, 8);
    % find how many number of clear frequency line region
    tmp_max = max(tmp, [], 'all');
    
    for k = 1 : tmp_max
        if size( find( tmp==k ) ) < segmentation % Segmentation Variable 
            tmp( find( tmp==k ) ) = 0; 
        else
            tmp( find( tmp==k ) ) = 1;
        end
    end
    
    % store data in Clear Frequecny matrix
    ClearFrequency = [ ClearFrequency ; tmp ]; 
end

% show the plot of siganl after Gabor Transform 
figure; 
image(t,f,abs(y)/max(max(abs(y)))*400)			
colormap(gray(256))								
set(gca,'Ydir','normal')						
set(gca,'Fontsize',12)							
xlabel('Time (Sec)','Fontsize',12)				
ylabel('Frequency (Hz)','Fontsize',12)	

% show the plot of Instantaneous Frequency 
figure; 
image(t,f,abs(InstantaneousFreq)/max(max(abs(InstantaneousFreq)))*400)			
colormap(gray(256))								
set(gca,'Ydir','normal')						
set(gca,'Fontsize',12)							
xlabel('Time (Sec)','Fontsize',12)				
ylabel('Frequency (Hz)','Fontsize',12)	

% show the plot of Clear Instantaneous Frequency 
figure; 
image(t,f,abs(ClearFrequency)/max(max(abs(ClearFrequency)))*400)			
colormap(gray(256))								
set(gca,'Ydir','normal')						
set(gca,'Fontsize',12)							
xlabel('Time (Sec)','Fontsize',12)				
ylabel('Frequency (Hz)','Fontsize',12)	


% ===================================================================
% Find the MAX time sample number of all instantaneous frequency
% ===================================================================

% create FreqTimeSample data matrix for later operation
FreqTimeSample = []; 

% label the different region of clear frequemcy line
ClearFrequency_label = bwlabel( ClearFrequency, 8 );

% find how many number of clear frequency line region
label_max = max(ClearFrequency_label, [], 'all');

for i = 1 : label_max
    % find the number of time sample of every instantaneous frequency
    FreqTimeSample_tmp = size( find(ClearFrequency_label == i) , 1);
    % store data in FreqTimeSample matrix
    FreqTimeSample = [ FreqTimeSample FreqTimeSample_tmp ];
end

% find the MAX time sample number of instantaneous frequency 
MAX_FreqTimeSample = max(FreqTimeSample);


% ===================================================================
% Find Real Frequency Value and its Time data of all frequency
% ===================================================================

% create Real Frequency Value data matrix for later operation
RealFreqValue = [];

% TimeRange will lose 1 data after sampling (round down) 
MaxSampleNum = MAX_FreqTimeSample - 1;

% initial setting of Time Data for all time sample of frequency  
TimeData = zeros(label_max, MaxSampleNum);

% initial setting of Frequency Value Data for all time sample of frequency 
% indexes of FreqValueData are correspond to its TimeData
FreqValueData = zeros(label_max , MaxSampleNum);

[ row, column ] = size(ClearFrequency);

for i = 1 : label_max 
    
    [r, c] = find( ClearFrequency_label == i );
    % calculate the Real Frequency Value of Instaneous frequency
    FreqValue = mean(r) + 20; % correction coefficient ylim [20 1001] 
    
    % store data in Real Frequency Value matrix
    RealFreqValue = [ RealFreqValue FreqValue ];
    
    % calculate the Real Time Value of every sample and find its min and max
    c_min = min(c) * (1.6/column);
    c_max = max(c) * (1.6/column);
    
    % TimeRange store the data of all time sample for all Freq region
    TimeRange = [ c_min : dt : c_max ];
    
    for k = 1 : size( TimeRange,2 ) % how many time sample in this region
        TimeData(i,k) = TimeRange(k);
        FreqValueData(i,k) = FreqValue;
    end
end


% ===================================================================
% Find frequency energy ratio of all primary frequency region 
% ===================================================================

E1_set = [];
E2_set = [];
E3_set = [];

EnergyRatioE2E1_set = [];
EnergyRatioE3E1_set = [];

E_onset_set = [];
E_offset_set = [];
E_all_set =[];

EnergyRatioOn_set = [];
EnergyRatioOff_set = [];


for FrequencyRegion = 1 : label_max
    
    y_value = (abs(y)).^2; 
    [r, c] = find( ClearFrequency_label == FrequencyRegion );
    
    fmin1 = 0.5 * RealFreqValue(FrequencyRegion);
    fmax1 = 1.5 * RealFreqValue(FrequencyRegion);
    
    min_row1 = round(fmin1) - 20;
    max_row1 = round(fmax1) - 20;
    
    if min_row1 < 1
        min_row1 = 1;
    elseif max_row1 > F
        max_row1 = F;
    end
    
    E1 = sum( y_value( min_row1 : max_row1 , min(c) : max(c) ) , 'all');
    E1_set = [ E1_set E1 ]; 
    
    fmin2 = 1.5 * RealFreqValue(FrequencyRegion);
    fmax2 = 2.5 * RealFreqValue(FrequencyRegion);
    
    min_row2 = round(fmin2) - 20;
    max_row2 = round(fmax2) - 20;
    
    if min_row2 < 1
        min_row2 = 1;
    elseif max_row2 > F
        max_row2 = F;
    end
    
    E2 = sum( y_value( min_row2 : max_row2 , min(c) : max(c) ) , 'all');
    E2_set = [ E2_set E2 ];
    
    fmin3 = 2.5 * RealFreqValue(FrequencyRegion);
    fmax3 = 3.5 * RealFreqValue(FrequencyRegion);
    
    min_row3 = round(fmin3) - 20;
    max_row3 = round(fmax3) - 20;
    
    if min_row3 < 1
        min_row3 = 1;
    elseif max_row3 > F
        max_row3 = F;
    end

​    E3 = sum( y_value( min_row3 : max_row3 , min(c) : max(c) ) , 'all');
​    E3_set = [ E3_set E3 ];
​    
​    EnergyRatioE2E1 = E2 / E1;
​    EnergyRatioE2E1_set = [ EnergyRatioE2E1_set EnergyRatioE2E1 ];
​    
​    EnergyRatioE3E1 = E3 / E1;
​    EnergyRatioE3E1_set = [ EnergyRatioE3E1_set EnergyRatioE3E1 ];
​    
​    begin_all = min(c); 
​    finish_all = max(c);
​    begin_onset = min(c);
​    finish_onset = min(c) + round(( (max(c)-min(c)) * 0.1 ));
​    begin_offset = min(c) + round(( (max(c)-min(c)) * 0.9 ));
​    finish_offset = max(c);
​        
​    LowLimit = round( ( RealFreqValue(FrequencyRegion) -20 ) - 5 );  
​    HighLimit = round( ( RealFreqValue(FrequencyRegion) -20 ) + 5 );
​    
​    if ( LowLimit < 1 )
​        LowLimit = 1;
​    elseif (HighLimit > F )
​        HighLimit = F;
​    end
​    
​    E_onset = sum( y_value( LowLimit:HighLimit, begin_onset:finish_onset ) , 'all');
​    E_onset_set = [ E_onset_set E_onset ];
​    
​    E_offset = sum( y_value( LowLimit:HighLimit, begin_offset:finish_offset ) , 'all');
​    E_offset_set = [ E_offset_set E_offset ];
​      
​    E_all = sum( y_value( LowLimit:HighLimit, begin_all:finish_all ) , 'all');
​    E_all_set =[ E_all_set E_all ];
​     
​    EnergyRatioOn = E_onset / E_all;
​    EnergyRatioOn_set = [ EnergyRatioOn_set EnergyRatioOn ];
​    
​    EnergyRatioOff = E_offset / E_all;
​    EnergyRatioOff_set = [ EnergyRatioOff_set EnergyRatioOff ];
​    
end


% ===================================================================
% Find error between Real Frequency & its second-order approximation
% ===================================================================

% create error data matrix of each frequency region for later operation
error_set = []; 

for FrequencyRegion = 1 : label_max
    
    % create Frequency column vector data of every sample point
    B = [];

​    for i = 1 : MaxSampleNum
​        
​        % extract the TimeData of every sample point row by row 
​        n_FreqTimeSample = TimeData( FrequencyRegion, : );
​        % find B = [ m m*n m*n^2 ]' column vector for approximation
​        B_tmp = RealFreqValue( FrequencyRegion ) * n_FreqTimeSample(i) .^ (0:2);
​        
​        % the element of B_tmp that not exist is time = 0
​        if ( B_tmp(2) == 0 ) && ( B_tmp(3) == 0) 
​            B_tmp(1) = 0;
​        end
​        % change the element of not exist B_tmp become zero column vector 
​        
​        B = [ B B_tmp'];
​    
​    end 
​        
​    % create approximation coefficient for every sample point
​    S = [];

​    for i = 1 : MaxSampleNum 
​        
​        % extract the Frequemcy data of every region column by column
​        B_column = B(:,i);
​        
​        % create A matrix for second-order approximation
​        for n = 1:3 
​            for m = 1:3
​                for k = 1:3
​                    A(m,k) = ( n_FreqTimeSample(i) )^( m+k-2 );
​                end
​            end
​        end
​        
​        % lower precision for lower computation
​        A = round(A,1);
​        % inverse A matrix for calculate approximation coefficient
​        S_tmp = pinv(A) * B_column;           
​        S = [S S_tmp];
​        
​    end

​    % create approzimation B of each frequency region for later operation
​    B_approximation_set= [];
​    
​    for i = 1 : MaxSampleNum 
​        
​        % extract the Coefficient data of every region column by column
​        S_column = S(:,i);
​        
​        % create A matrix for second-order approximation 
​        for n = 1:3 
​            for m = 1:3
​                for k = 1:3
​                    A(m,k) = ( n_FreqTimeSample(i) )^( m+k-2 );
​                end
​            end
​        end
​        
​        % approximate the frequency of each time samplefor each region 
​        B_approximation_matrix = A*S_column;
​        B_approximation_set = [ B_approximation_set B_approximation_matrix ];
​        
​    end
​    
​    for i = 1 : MaxSampleNum
​        
​        % some sample point does not exist so its B_approximation will be 0
​        % convert 0 to eps to avoid later operation problem (0/0=InF) 
​        if ( B_approximation_set(1,i) == 0 )
​            B_approximation_set(1,i) = eps;
​        end
​        
​        % some frequency point does not exist so its Freq data will be 0
​        % convert 0 to eps to avoid later operation problem (0/0=InF) 
​        if ( FreqValueData(FrequencyRegion,i) == 0 )
​            FreqValueData(FrequencyRegion,i) = eps;
​        end
​        
​    end
​    
​    % extract the data of approximation frequency we need in the first row
​    B_approximation = B_approximation_set(1,:);
​    % extract the data of real frequency we need row by row
​    FreqValueData_row = FreqValueData(FrequencyRegion,:);
​    
​    % calculate error data row by row for each region 
​    error = sum(abs((B_approximation-FreqValueData_row) ./ FreqValueData_row)) / (MaxSampleNum);
​    error_percentage = round( error*100, 2 );
​    
​    error_set = [ error_set error_percentage ];
end

% find total error for all frequency region 
total_error = sum(error_set)/size(error_set,2);

AllData = [ RealFreqValue ; EnergyRatioE2E1_set ; EnergyRatioE3E1_set; ...
            EnergyRatioOn_set ; EnergyRatioOff_set ; error_set ]

HighFreqE2E1 = [];  
HighFreqE3E1 = [];
HighFreqOn = [];  
HighFreqOff = [];  
HighFreqError = [];

MiddleFreqE2E1 = [];
MiddleFreqE3E1 = [];
MiddleFreqOn = [];
MiddleFreqOff = [];
MiddleFreqError = [];

LowFreqE2E1 = [];  
LowFreqE3E1 = [];
LowFreqOn = [];  
LowFreqOff = [];
LowFreqError = [];

for i = 1 : label_max
    
    if AllData(1,i) > 1000       
        HighFreqE2E1 = [ HighFreqE2E1  AllData(2,i) ];  
        HighFreqE3E1 = [ HighFreqE3E1  AllData(3,i) ];  
        HighFreqOn = [ HighFreqOn  AllData(4,i) ] ; 
        HighFreqOff = [ HighFreqOff  AllData(5,i) ];  
        HighFreqError = [ HighFreqError  AllData(6,i) ]; 
        
    elseif AllData(1,i) < 500
        LowFreqE2E1 = [ LowFreqE2E1  AllData(2,i) ];  
        LowFreqE3E1 = [ LowFreqE3E1  AllData(3,i) ];  
        LowFreqOn = [ LowFreqOn  AllData(4,i) ] ; 
        LowFreqOff = [ LowFreqOff  AllData(5,i) ];  
        LowFreqError = [ LowFreqError  AllData(6,i) ]; 
        
    else 
        MiddleFreqE2E1 = [ MiddleFreqE2E1  AllData(2,i) ];  
        MiddleFreqE3E1 = [ MiddleFreqE3E1  AllData(3,i) ];  
        MiddleFreqOn = [ MiddleFreqOn  AllData(4,i) ] ; 
        MiddleFreqOff = [ MiddleFreqOff  AllData(5,i) ];  
        MiddleFreqError = [ MiddleFreqError  AllData(6,i) ]; 
    end
end


MeanValue1 = 0;
HighFreqE2E1_sort_set = [];
for i = 1 : size( HighFreqE2E1, 2 )
    HighFreqE2E1_sort = sort( HighFreqE2E1 );
    if HighFreqE2E1_sort(:,i) < 2 
       MeanValue1 = ( MeanValue1*(i-1) + HighFreqE2E1_sort(:,i) )/i;
    else
       HighFreqE2E1_sort(:,i) = MeanValue1;
    end
    HighFreqE2E1_sort_set = [ HighFreqE2E1_sort_set  HighFreqE2E1_sort(i) ];
end


MeanValue2 = 0;
HighFreqE3E1_sort_set = [];
for i = 1 : size( HighFreqE3E1, 2 )
    HighFreqE3E1_sort = sort( HighFreqE3E1 );
    if HighFreqE3E1_sort(:,i) < 2 
       MeanValue2 = ( MeanValue2*(i-1) + HighFreqE3E1_sort(:,i) )/i;
    else
       HighFreqE3E1_sort(:,i) = MeanValue1;
    end
    HighFreqE3E1_sort_set = [ HighFreqE3E1_sort_set  HighFreqE3E1_sort(i) ];
end


MeanValue3 = 0;
HighFreqOn_sort_set = [];
for i = 1 : size( HighFreqOn, 2 )
    HighFreqOn_sort = sort( HighFreqOn );
    if HighFreqOn_sort(:,i) < 2 
       MeanValue3 = ( MeanValue3*(i-1) + HighFreqOn_sort(:,i) )/i;
    else
       HighFreqOn_sort(:,i) = MeanValue3;
    end
    HighFreqOn_sort_set = [ HighFreqOn_sort_set  HighFreqOn_sort(i) ];
end


MeanValue4 = 0;
HighFreqOff_sort_set = [];
for i = 1 : size( HighFreqOff, 2 )
    HighFreqOff_sort = sort( HighFreqOff );
    if HighFreqOff_sort(:,i) < 2 
       MeanValue4 = ( MeanValue4*(i-1) + HighFreqOff_sort(:,i) )/i;
    else
       HighFreqOff_sort(:,i) = MeanValue4;
    end
    HighFreqOff_sort_set = [ HighFreqOff_sort_set  HighFreqOff_sort(i) ];
end


MeanValue5 = 0;
HighFreqError_sort_set = [];
for i = 1 : size( HighFreqError, 2 )
    HighFreqError_sort = sort( HighFreqError );
    if HighFreqError_sort(i) < 10 
       MeanValue5 = ( MeanValue5*(i-1) + HighFreqError_sort(i) )/i;
    else
       HighFreqError_sort(i) = MeanValue5;
    end
    HighFreqError_sort_set = [ HighFreqError_sort_set HighFreqError_sort(i) ];
end


MeanValue6 = 0;
MiddleFreqE2E1_sort_set = [];
for i = 1 : size( MiddleFreqE2E1, 2 )
    MiddleFreqE2E1_sort = sort( MiddleFreqE2E1 );
    if MiddleFreqE2E1_sort(:,i) < 2 
       MeanValue6 = ( MeanValue6*(i-1) + MiddleFreqE2E1_sort(:,i) )/i;
    else
       MiddleFreqE2E1_sort(:,i) = MeanValue1;
    end
    MiddleFreqE2E1_sort_set = [ MiddleFreqE2E1_sort_set  MiddleFreqE2E1_sort(i) ];
end


MeanValue7 = 0;
MiddleFreqE3E1_sort_set = [];
for i = 1 : size( MiddleFreqE3E1, 2 )
    MiddleFreqE3E1_sort = sort( MiddleFreqE3E1 );
    if MiddleFreqE3E1_sort(:,i) < 2 
       MeanValue7 = ( MeanValue7*(i-1) + MiddleFreqE3E1_sort(:,i) )/i;
    else
       MiddleFreqE3E1_sort(:,i) = MeanValue7;
    end
    MiddleFreqE3E1_sort_set = [ MiddleFreqE3E1_sort_set  MiddleFreqE3E1_sort(i) ];
end


MeanValue8 = 0;
MiddleFreqOn_sort_set = [];
for i = 1 : size( MiddleFreqOn, 2 )
    MiddleFreqOn_sort = sort( MiddleFreqOn );
    if MiddleFreqOn_sort(:,i) < 2 
       MeanValue8 = ( MeanValue8*(i-1) + MiddleFreqOn_sort(:,i) )/i;
    else
       MiddleFreqOn_sort(:,i) = MeanValue8;
    end
    MiddleFreqOn_sort_set = [ MiddleFreqOn_sort_set  MiddleFreqOn_sort(i) ];
end


MeanValue9 = 0;
MiddleFreqOff_sort_set = [];
for i = 1 : size( MiddleFreqOff, 2 )
    MiddleFreqOff_sort = sort( MiddleFreqOff );
    if MiddleFreqOff_sort(:,i) < 2 
       MeanValue9 = ( MeanValue9*(i-1) + MiddleFreqOff_sort(:,i) )/i;
    else
       MiddleFreqOff_sort(:,i) = MeanValue9;
    end
    MiddleFreqOff_sort_set = [ MiddleFreqOff_sort_set  MiddleFreqOff_sort(i) ];
end


MeanValue10 = 0;
MiddleFreqError_sort_set = [];
for i = 1 : size( MiddleFreqError, 2 )
    MiddleFreqError_sort = sort( MiddleFreqError );
    if MiddleFreqError_sort(i) < 10 
       MeanValue10 = ( MeanValue10*(i-1) + MiddleFreqError_sort(i) )/i;
    else
       MiddleFreqError_sort(i) = MeanValue10;
    end
    MiddleFreqError_sort_set = [ MiddleFreqError_sort_set  MiddleFreqError_sort(i) ];
end


MeanValue11 = 0;
LowFreqE2E1_sort_set = [];
for i = 1 : size( LowFreqE2E1, 2 )
    LowFreqE2E1_sort = sort( LowFreqE2E1 );
    if LowFreqE2E1_sort(:,i) < 2 
       MeanValue11 = ( MeanValue11*(i-1) + LowFreqE2E1_sort(:,i) )/i;
    else
       LowFreqE2E1_sort(:,i) = MeanValue11;
    end
    LowFreqE2E1_sort_set = [ LowFreqE2E1_sort_set LowFreqE2E1_sort(i) ];
end


MeanValue12 = 0;
LowFreqE3E1_sort_set = [];
for i = 1 : size( LowFreqE3E1, 2 )
    LowFreqE3E1_sort = sort( LowFreqE3E1 );
    if LowFreqE3E1_sort(:,i) < 2 
       MeanValue12 = ( MeanValue12*(i-1) + LowFreqE3E1_sort(:,i) )/i;
    else
       LowFreqE3E1_sort(:,i) = MeanValue12;
    end
    LowFreqE3E1_sort_set = [ LowFreqE3E1_sort_set LowFreqE3E1_sort(i) ];
end


MeanValue13 = 0;
LowFreqOn_sort_set = [];
for i = 1 : size( LowFreqOn, 2 )
    LowFreqOn_sort = sort( LowFreqOn );
    if LowFreqOn_sort(:,i) < 2 
       MeanValue13 = ( MeanValue13*(i-1) + LowFreqOn_sort(:,i) )/i;
    else
       lowFreqOn_sort(:,i) = MeanValue13;
    end
    LowFreqOn_sort_set = [ LowFreqOn_sort_set LowFreqOn_sort(i) ];
end


MeanValue14 = 0;
LowFreqOff_sort_set = [];
for i = 1 : size( LowFreqOff, 2 )
    LowFreqOff_sort = sort( LowFreqOff );
    
    if LowFreqOff_sort(:,i) < 2 
       MeanValue14 = ( MeanValue14*(i-1) + LowFreqOff_sort(:,i) )/i;
    else
       LowFreqOff_sort(:,i) = MeanValue14;
    end
    LowFreqOff_sort_set = [ LowFreqOff_sort_set  LowFreqOff_sort(i) ];
end


MeanValue15 = 0;
LowFreqError_sort_set = [];
for i = 1 : size( LowFreqError, 2 )
    LowFreqError_sort = sort( LowFreqError );
    if LowFreqError_sort(i) < 10 
       MeanValue15 = ( MeanValue15*(i-1) + LowFreqError_sort(i) )/i;
    else
       LowFreqError_sort(i) = MeanValue15;
    end
    LowFreqError_sort_set = [ LowFreqError_sort_set  LowFreqError_sort(i) ];
end


HighFreqE2E1_mean = mean( HighFreqE2E1_sort_set )
HighFreqE3E1_mean = mean( HighFreqE3E1_sort_set )
HighFreqOn_mean = mean( HighFreqOn_sort_set )
HighFreqOff_mean = mean( HighFreqOff_sort_set )
HighFreqError_mean = mean( HighFreqError_sort_set )

MiddleFreqE2E1_mean = mean( MiddleFreqE2E1_sort_set )
MiddleFreqE3E1_mean = mean( MiddleFreqE3E1_sort_set )
MiddleFreqOn_mean = mean( MiddleFreqOn_sort_set )
MiddleFreqOff_mean = mean( MiddleFreqOff_sort_set ) 
MiddleFreqError_mean = mean( MiddleFreqError_sort_set )

LowFreqE2E1_mean = mean( LowFreqE2E1_sort_set )
LowFreqE3E1_mean = mean( LowFreqE3E1_sort_set )
LowFreqOn_mean = mean( LowFreqOn_sort_set )
LowFreqOff_mean = mean( LowFreqOff_sort_set ) 
LowFreqError_mean = mean( LowFreqError_sort_set )

SignalFeature = [ HighFreqE2E1_mean, HighFreqE3E1_mean, ... 
                  HighFreqOn_mean, HighFreqOff_mean, HighFreqError_mean, ...
                  MiddleFreqE2E1_mean, MiddleFreqE3E1_mean, ... 
                  MiddleFreqOn_mean, MiddleFreqOff_mean, MiddleFreqError_mean, ...
                  LowFreqE2E1_mean, LowFreqE3E1_mean, ... 
                  LowFreqOn_mean, LowFreqOff_mean, LowFreqError_mean ]


SignalFeature_kmeans = lbgVQ( SignalFeature, 4);

load('library_cello1.mat')
d1_set = [];
for i = 1 : size(library_cello1,3)    
    d1(i) = distanceEu( SignalFeature, library_cello1(:,:,i) );
    d1_set = [ d1_set d1(i) ];
end
d1_data = sum( d1_set )/size(library_cello1,3);


load('library_cello2.mat')
d2_set = [];
for i = 1 : size(library_cello2,3)    
    d2(i) = distanceEu( SignalFeature, library_cello2(:,:,i) );
    d2_set = [ d2_set d2(i) ];
end
d2_data = sum( d2_set )/size(library_cello2,3);


load('library_guitar1.mat')
d3_set = [];
for i = 1 : size(library_guitar1,3)
    d3(i) = distanceEu( SignalFeature, library_guitar1(:,:,i) );
    d3_set = [ d3_set d3(i) ];
end
d3_data = sum( d3_set )/size(library_guitar1,3);


load('library_guitar2.mat')
d4_set = [];
for i = 1 : size(library_guitar2,3)
    d4(i) = distanceEu( SignalFeature, library_guitar2(:,:,i) );
    d4_set = [ d4_set d4(i) ];
end
d4_data = sum( d4_set )/size(library_guitar2,3);


load('library_guitar3.mat')
d5_set = [];
for i = 1 : size(library_guitar3,3)
    d5(i) = distanceEu( SignalFeature, library_guitar3(:,:,i) );
    d5_set = [ d5_set d5(i) ];
end
d5_data = sum( d5_set )/size(library_guitar3,3);


load('library_piano1.mat')
d6_set = [];
for i = 1 : size(library_piano1,3)
    d6(i) = distanceEu( SignalFeature, library_piano1(:,:,i) );
    d6_set = [ d6_set d6(i) ];
end
d6_data = sum( d6_set )/size(library_piano1,3);


load('library_piano2.mat')
d7_set = [];
for i = 1 : size(library_piano2,3)
    d7(i) = distanceEu( SignalFeature, library_piano2(:,:,i) );
    d7_set = [ d7_set d7(i) ];
end
d7_data = sum( d7_set )/size(library_piano2,3);


load('library_piano3.mat') % electro piano 
d8_set = [];
for i = 1 : size(library_piano3,3)
    d8(i) = distanceEu( SignalFeature, library_piano3(:,:,i) );
    d8_set = [ d8_set d8(i) ];
end
d8_data = sum( d8_set )/size(library_piano3,3);


% Recoginition
RecogData = [ d1_data d2_data d3_data d4_data d5_data d6_data d7_data d8_data ]

if ( min( RecogData ) == d1_data ) || ( min( RecogData ) == d2_data )
    fprintf('The signal is Cello. ')

elseif ( min( RecogData ) == d3_data ) || ( min( RecogData ) == d4_data ) || ( min( RecogData ) == d5_data )
    fprintf('The signal is Guitar. ')

elseif ( min( RecogData ) == d6_data ) || ( min( RecogData ) == d7_data ) || ( min( RecogData ) == d8_data )
    fprintf('The signal is Piano. ')
end