%% Project of EEL6586 --- LPC

clear;
close all;
trainData = load('trainTemplate.mat');
trainNames = fieldnames(trainData);


N = 10;
windowSize = 200;
overlap = 100;
hamWindow = 0.54 - 0.46*cos(2*pi*(0:windowSize-1)/windowSize);
Alpc = zeros(6,30,500);

for s = 1 : length(trainNames)%% s is the number of sentences
    name = trainNames{s};
    data = getfield(trainData,name);
    
    %% noise reduction
    for m = 1 : size(data,1) %% m is number of training data for each sentence
        %% pre-emphasis
        preData(m,:) = filter([1 -.97], 1, data(m,:));
        drvData(m,:) = conv2(preData(m,:), [0.333 0.666 0 -0.666 -0.333], 'same');
        winData(m,:) = wiener2(drvData(m,:));
        
        %% use end-point to delete silence
        Qs = filter(hamWindow,1, winData(m,:).* winData(m,:)*100);
        [dfoo, rang] = find(Qs >= 0.005);
        
        if rang(301) - rang(1) == 300
            rang1 = rang(1);
        else
            mm = 1;
            while rang(mm+1) - rang(mm) ==1 & rang(mm+300) - rang(mm) >300
                mm=mm+1;
            end
            rang1 = rang(mm+300);
        end
        
        if rang(length(rang)) - rang(length(rang)-300) == 300;
            rang2 = rang(length(rang));
        else
            mm = 0;
            while rang(length(rang)-mm) - rang(length(rang)-mm-1) ==1 & rang(length(rang)-mm) - rang(length(rang)-mm-300) > 300 
                mm = mm + 1;
            end
            rang2 = rang(length(rang)-mm-1);
        end
        
        trunkData = winData(m,rang1 : rang2);
        
        %% add windows to each training data and calculate LPC for each window
        LPCtemp = [];
        for n = 1 : (windowSize-overlap) : floor(length(trunkData) / windowSize-1) * windowSize
            windowed = trunkData(n : n + windowSize - 1);
            windowed = windowed .* hamWindow;
            LPCtemp((n-1) / (windowSize-overlap) + 1,:) = real(lpc(windowed, N));
        end
        LPCtemp(:,1) = []; 
        
        %% calculate autocorrelation of each training data's 4th LPC coeff.
        Alpc(s,m,:) = [LPCtemp(:,4)', zeros(1,500-size(LPCtemp,1))];
        
        LPCtrain = squeeze(Alpc(s,m,:)); %% LPCtrain is current data to be trained for model
        R(s,m,:) = xcorr(LPCtrain,LPCtrain);
        %aLPC(s,m,:) = R(s,m,500:509);
        aLPC(s,m,:) = sum(R(s,m,:));
    end
end

%% Bayes
a = length(trainNames);
for s = 1 : a       %% tatol a sentenses
    A = squeeze(aLPC(s,:,:))';
    miu(s,:) = mean(A);
    b = size(A,1);  %% there are b speakers saying each sentense
    sigma = zeros(size(A,2));
    
    for n = 1 : b
        temp = (A(n,:) - miu(s,:))' * (A(n,:) - miu(s,:));
        sigma = sigma + temp;
    end
    
    sigma = sigma / (b-1);
    
    Sigma(s,:,:) = sigma;
    
end


%% Fisher-linear
%% w{m,n} is the projection line for class m and n
for m = 1 : a-1
    for n = m+1 : a
        Sw{m,n} = squeeze(Sigma(m,:,:)) + squeeze(Sigma(n,:,:));
        w{m,n} = inv(Sw{m,n}) * (miu(m,:) - miu(n,:))';
        w{m,n} = w{m,n} / sqrt((w{m,n})' * w{m,n});
        w{n,m} = w{m,n};
    end
end

%% training data projection on w
%% proj{m,n} is projection of class m on w{m,n}
%% proj{n,m} is projection of class n on w{m,n}
for s = 1 : a
    A = squeeze(aLPC(s,:,:));
    for n = 1 : a
        if s ~= n
            for k = 1 : b
                position(k) =  w{s,n}' * A(k,:)'; 
            end
            proj{s,n} = position;
        end
    end
end

%% W0 for each classifier
for m = 1: a-1
    for n = m+1 : a
        MN = sort(proj{m,n});
        NM = sort(proj{n,m});
        if MN(1) <= NM(b) & MN(1) > NM(1)
            MNmix = sort([MN, NM]);
            errMN = []; errNM = [];
            tail = find(MNmix == NM(b));
            head = find(MNmix == MN(1));
            t = tail - head +1;
            errMN(1) = 1 ; errNM = t;
            for i = 1 : t
                if ~isempty(find(MNmix(head + i - 1) == MN)) %% "*" found
                    errMN(i+1) = errMN(i) + 1;
                    errNM(i+1) = errNM(i) - 1;
                end
                if ~isempty(find(MNmix(head + i - 1) == NM)) %% "o" found
                    errMN(i+1) = errMN(i);
                    errNM(i+1) = errNM(i) -1;
                end
            end
            [foo1, foo2] = min(errMN + errNM);
            w0{m,n} = MNmix(head+foo2-1)/2 + MNmix(head+foo2-2)/2;
            w0{n,m} = w0{m,n};
            
        elseif NM(1) <= MN(b) & NM(1) > MN(1)
            NMmix = sort([MN, NM]);
            errNM = []; errMN = [];
            tail = find(NMmix == MN(b));
            head = find(NMmix == NM(1));
            t = tail - head +1;
            errNM(1) = 1 ; errMN = t;
            for i = 1 : t
                if ~isempty(find(NMmix(head + i - 1) == NM)) %% "*" found
                    errNM(i+1) = errNM(i) + 1;
                    errMN(i+1) = errMN(i) - 1;
                end
                if ~isempty(find(NMmix(head + i - 1) == MN)) %% "o" found
                    errNM(i+1) = errNM(i);
                    errMN(i+1) = errMN(i) -1;
                end
            end
            [foo1, foo2] = min(errNM + errMN);
            w0{m,n} = NMmix(head+foo2-1)/2 + NMmix(head+foo2-2)/2;
            w0{n,m} = w0{m,n};
            
        elseif MN(1) > NM(b)
            w0{m,n} = MN(1)/2 + NM(b)/2;
            w0{n,m} = w0{m,n};
            
        elseif NM(1) > MN(b)
            w0{m,n} = NM(1)/2 + MN(b)/2;
            w0{n,m} = w0{m,n};
            
        end
    end
end



%%##################################################################################################

%% test recorded data
testData = load('testTemplate.mat');
testNames = fieldnames(testData);

for s = 1 : length(testNames)%% s is the number of sentences
    name = testNames{s};
    data = getfield(testData,name);
    
    %% noise reduction
    for m = 1 : size(data,1) %% m is number of training data for each sentence
        %% pre-emphasis
        preData(m,:) = filter([1 -.97], 1, data(m,:));
        drvData(m,:) = conv2(preData(m,:), [0.333 0.666 0 -0.666 -0.333], 'same');
        winData(m,:) = wiener2(drvData(m,:));
        
        %% use end-point to delete silence
        Qs = filter(hamWindow,1, winData(m,:).* winData(m,:)*100);
        [dfoo, rang] = find(Qs >= 0.005);
        
        if rang(301) - rang(1) == 300
            rang1 = rang(1);
        else
            mm = 1;
            while rang(mm+1) - rang(mm) ==1 & rang(mm+300) - rang(mm) >300
                mm=mm+1;
            end
            rang1 = rang(mm+300);
        end
        
        if rang(length(rang)) - rang(length(rang)-300) == 300;
            rang2 = rang(length(rang));
        else
            mm = 0;
            while rang(length(rang)-mm) - rang(length(rang)-mm-1) ==1 & rang(length(rang)-mm) - rang(length(rang)-mm-300) > 300 
                mm = mm + 1;
            end
            rang2 = rang(length(rang)-mm-1);
        end
        
        trunkData = winData(m,rang1 : rang2);
        
        %% LPC coeff
        
        temp2 = real(lpc(trunkData,N));
        training = temp2(2:N+1);
        
        flag = 0;
        
        k = 1;
        while k <= a
            
            n = 1;
            while n <= a
                
                if k ~= n
                    testScore = training * w{k,n};
                    if k<n & testScore>=w0{k,n}
                        claim(s,m) = k;
                        
                    elseif k>n & testScore<=w0{k,n}
                        claim(s,m) = k;
                        
                    elseif k<n & testScore<=w0{k,n}
                        claim(s,m) = n; k = n-1;
                        flag = 0.5;  %% normal jump
                        n = a + 2;
                        
                    elseif k>n & testScore>=w0{k,n} 
                        if flag == 0.5              
                            claim(s,m) = 0;         %% reject temperaly
                            n = a+2; %k = a+2;       %% set value to get next k
                            flag = 1;               %% abnormal jump
                        else                           
                            claim(s,m) = n-1;
                            %flag = 1;
                        end
                    end
                end
                n = n + 1;
                
            end
            
            if flag ~= 1 & n == a+1        %% normal finish
                claim(s,m) = k;            %% jump out to get next k
                k = a+1; n = a+1;
            elseif flag ==1 & n~=a+1       %% if currently rejected
                flag = 0;                   %% clear flag
                %k = k+1;                   %% and go to the next entity
            end
            
            k = k+1;                       %% go to netx entity   
        end
        
        fprintf('result = %d\n\n', claim(s,m));
    end
end


%############################################################################################


%% test training data
for s = 1 : a 
    A = squeeze(aLPC(s,:,:));
    
    for k = 1 : b
        training = A(k,:);
        
        flag = 0;
        
        m = 1;
        while m <= a
            
            n = 1;
            while n <= a
                
                if m ~= n
                    testScore = training * w{m,n};
                    if m<n & testScore>=w0{m,n}
                        claim(s,k) = m;
                        
                    elseif m>n & testScore<=w0{m,n}
                        claim(s,k) = m;
                        
                    elseif m<n & testScore<=w0{m,n}
                        claim(s,k) = n; m = n-1;
                        flag = 0.5;  %% normal jump
                        n = a + 2;
                        
                    elseif m>n & testScore>=w0{m,n} 
                        if flag == 0.5              
                            claim(s,k) = 0;         %% reject temperaly
                            n = a+2; %m = a+2;       %% set value to get next k
                            flag = 1;               %% abnormal jump
                        else                           
                            claim(s,k) = n-1;
                            %flag = 1;
                        end
                    end
                end
                n = n + 1;
                
            end
            
            if flag ~= 1 & n == a+1        %% normal finish
                claim(s,k) = m;            %% jump out to get next k
                m = a+1; n = a+1;
            elseif flag ==1 & n~=a+1       %% if currently rejected
                flag = 0;                   %% clear flag
                %m = m+1;                   %% and go to the next entity
            end
            
            m = m+1;                       %% go to netx entity   
        end
        
    end
end