% Dmitriy Yakovlev
% speacorr.m - speaker identification based on neural networks

% WINDOWED LPC VERSION

% get the data from the files

labels = ['d';'y';'b';'s';'z';'e';'l']; 

% these are the available labels.. the amount of samples in each group is
% arbitrary and automatically determined ((to be 5))

% d Dmitriy Yakovlev
% y Mary Van Vleet
% b Autumn Petros-Good
% s Leah Anderson
% z Seanna Vine
% e Denis Aleshin
% l Stephen Lee

% if a different problem is desired (e.g. male-female classification), change 
% the labels to two classes and rename the voice sample files to match them

[data, Fs] = getdata(labels);

%multiply data by 10 to get bigger numbers
data = data * 10;

% targets are encoded by index in the labels array (e.g. a network output
% of 3 would indicate that the speaker is Autumn.

% targets is a one-hot encoding

d = [1 0 0 0 0 0 0];
y = [0 1 0 0 0 0 0];
b = [0 0 1 0 0 0 0];
s = [0 0 0 1 0 0 0];
z = [0 0 0 0 1 0 0];
e = [0 0 0 0 0 1 0];
l = [0 0 0 0 0 0 1];

targets = [d' d' d' d' d' y' y' y' y' y' b' b' b' b' b' s' s' s' s' s' z' z' z' z' z' e' e' e' e' e' l' l' l' l' l'];

% break data into windows, and run lpc on windows

% divide up a signal into windows

w = hamming(1600);                          % hamming window
nw = length(w);                            % size of window
nx = size(data,1);                         % length of all columns

lpcCoeffsacc = zeros(35,32);
frames = 0;

pos=1;

% run lpc on each window
while (pos+nw <= nx)                       % while enough signal left
    y = data(pos:(pos+nw-1),:);

    for z=1:35
        y(:,z) = y(:,z).*w;           % make window y
    end
   
    % run LPC on y
    lpcCoeffs = lpc(y, 31); % get 32nd order polynomial estimate of the data

    % lpc coeffs is now a matrix of 32x35 numbers 
    lpcCoeffsacc = (lpcCoeffsacc * frames + lpcCoeffs) / (frames + 1);

    frames = frames + 1;

    pos = pos + nw/2;                 % next window
end

% coeffs are now rows of 32 things, not columns, but newff takes input
% where the data vectors are columns, so flip lpcCoeffs accordingly


% NEURAL NETWORK PART

% phase 1 - load training data

[num_samples, num_input_attributes] = size(lpcCoeffs');
[num_samples, num_output_attributes] = size(targets');

% phase 2 - set up network
net = newff(lpcCoeffs', targets, [32],{'logsig'},'traincgb');
net.trainParam.show = 250; % show every epoch
net.trainParam.lr = 0.001; % learning rate
net.trainParam.epochs = 1000; 
net.trainParam.goal = 1e-6; % MSE to stop at
net.trainParam.mem_reduc = 2; % use less memory
net.trainParam.max_fail = 50000; % allow to fail lots of time... weird.
net.trainParam.minstep = 1e-8;

net.inputweights{1,1}.initFcn = 'rands';
net.biases{1}.initFcn = 'rands';
net = init(net);


% phase 3 - train network on data

[net,tr]=train(net,lpcCoeffs',targets);

% check how it does on the given data
calc_ans = sim(net,lpcCoeffs');

num_correct_samples = 0;

for i = 1:num_samples
    anscol = calc_ans(:,i); % grab the column
    if (max(anscol) == targets(:,i)'*anscol)
        num_correct_samples = num_correct_samples + 1;
    end  
end

calc_ans

fprintf('%d correct samples identified out of %d in training set\n', num_correct_samples, num_samples);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% now wav files are columns in the data matrix.  Compute DFT of all columns
% and get frequencies

% divide up a signal into windows

%w = hamming(160)';                          % hamming window
%nw = length(w);                            % size of window
%nx = size(data,1);                         % length of all columns

%fftdata = zeros(nx/(200), size(data,2));  % holds fft results

%pos=1;
%while (pos+nw <= nx)                       % while enough signal left
%    y = data(pos:pos+nw-1,:).*w;           % make window y

%    % perform DFT on that window (matrix DFT does it by columns)
%    Y = fft(y,nw);
%    Y = 2*abs(Y(1:(nw/2),:));  
    
    % take 4000 amplitudes and make 40 numbers out of them for each column
    
    % need a representative # every 100 Hz, so take the average of every
    % 100 entries of every column
    
%    blockavg = @(x) (sum(x)/size(x))*ones(size(x),1);
    
%    Yblock = blkproc(Y, [100 1], blockavg);
    
    % collapse array to get just the averages
%    for k=1:40
%        Yblock(((k-1)*100+1),:)
        
%    end
    

%    pos = pos + nw/2;                 % next window
%end

%for i=1:size(data,2)
%    y = data(:,i);                % actual signal samples
%    T = 1/Fs;                     % Sample time
%    L = size(y,1);                  % Length of signal
%    t = (0:L-1)*T;                % Time vector
%    
%    Y = fft(y,80);
%    f = Fs/2*linspace(0,1,40);
%    
%    plot(f,2*abs(Y(1:40))) 
%    title('Single-Sided Amplitude Spectrum of y(t)')
%    xlabel('Frequency (Hz)')
%    ylabel('|Y(f)|')
%    
%    pause
%end


