close all;clear all;
binNumber = 60;     % Number of bins in the histogram
roomNumber = 40;    % Number of rooms in the training data
recNumber = 30;     % Number of recording for each room in the training data  
fs = 16000;         % Sampling frequency
epsilon = 128;      % Scaling factor for the kernel
epsilonTest = 1024;   % Sacling factor for the test kernel

% Parameters
RTtrain = linspace(0.2,0.8,40);
lambdaTrain = ((-6)*log(10))./(RTtrain);
binRange = linspace(-10,10,binNumber);
covariance = zeros(binNumber,binNumber,roomNumber);
h = zeros(binNumber,roomNumber,recNumber);
% Loading the Training data
load('trainData.mat');
figure
hold on

for r =1:roomNumber
    % Initializing the empirical mean to zero for each
    hrbar = zeros(binNumber,1);
    fprintf('Room no.%d\n',r);
    for i =1:recNumber
        temp = trainData(:,r,i);    % ith recording for rth room
        [s,f,t,p] = spectrogram(temp,512,384);  % STFT 
        lambda = zeros(length(f),1);            % Initializing the lambdas to zero    
        energyEnvelope = log(abs(s).*abs(s));   % Calculating the energy envelope
        for k = 1:length(f)
            P = polyfit(t/fs,energyEnvelope(k,:),1);    % LS fit to the lambda
            lambda(k) = P(1);                   % Storing the different lambda
        end
        lambdaMatrix(r,i,:) = lambda;           %     
        h(:,r,i) = histc(lambda,binRange);      % Plotting the histogram     
        hrbar = hrbar + h(:,r,i);               % Accumulating the hist_feature for hbar
    end
    hrbar = hrbar/length(recNumber);
    plot(hrbar);
    % Calculating Co-variance matrix for each room
    for i=1:recNumber
        covariance(:,:,r) = covariance(:,:,r) + (h(:,r,i) - hrbar)*((h(:,r,i)-hrbar)');
    end
    covariance(:,:,r)  = covariance(:,:,r)/recNumber;
end
% Starting the kernal Matrix Calculation 
kernalDim = roomNumber * recNumber;
kernalMatrix = zeros(kernalDim);
distanceMatrix = zeros(kernalDim);
covarianceInv = zeros(binNumber,binNumber,roomNumber);
% Pre-calculating the inv of covariance matrix
for i=1:roomNumber
    covarianceInv(:,:,i) = pinv(covariance(:,:,i));
end

fprintf('Calculating the Kernel Matrix\n');
for r = 0:(roomNumber-1)
    for rho = 0:(roomNumber-1)
        for i = 1:recNumber
            for j = 1:recNumber
                n = r*recNumber + i;
                m = rho*recNumber + j;
                distance = (mahaDistance(h(:,r+1,i),h(:,rho+1,j),covarianceInv(:,:,r+1),covarianceInv(:,:,rho+1)));
                distanceMatrix(n,m) = distance;
                kernalMatrix(n,m) = exp(distance/((-1)*epsilon));
            end
        end
    end
    fprintf('room %d\n',r);
end
% Extracting the sum of rows into
diagElement = sum(kernalMatrix,2);
D = diag(((diagElement).^(-0.5)));  %Calculating the D matrix

Wbar = D*kernalMatrix*D;            % Calculating the W matrix
[eigenVector, eigenValue] = eig(Wbar);  % Eigendecomposition
% Plotting the values of eigenvector v/s reverberation time
x = 1:30:1200;
figure
plot(RTtrain,eigenVector(x,2),'o');
xlabel('Reverberation Time');ylabel('Eigenvector');
fprintf('Eigendecomposition completed\n');
% Test code starts
fprintf('Starting the testing procedure\n');
load('testData.mat')        % Loading the test data    
numTestCase = size(testData,2); % 
for u=1:numTestCase
    temp = testData(:,u);       % Loading the uth test case    
    [s,f,t,p] = spectrogram(temp,256,192);  % STFT
    lambda = zeros(length(f),1);            
    energyEnvelope = log(abs(s).*abs(s));   % Calculating energy envelope
    for i = 1:length(f)
        P = polyfit(t/fs,energyEnvelope(i,:),1);  % LS fit to calculate the lambda  
        lambda(i) = P(1);
    end
    hTest(:,u) = histc(lambda,binRange);     % Calculating the histogram   
end
% Calculating non-symmetric affinity matrix 
A = zeros(numTestCase,recNumber*roomNumber);
for u=1:numTestCase
    for r=0:(roomNumber-1)
        for i=1:recNumber
            n = r*recNumber + i;
            A(u,n) = exp((mahaDistance(hTest(:,u), h(:,r+1,i), covarianceInv(:,:,r+1), covarianceInv(:,:,r+1)))/(2*epsilonTest));
        end
    end
end
% Calculating the Da diagnoal matrix containing sum of rows of A
DtestVec = zeros(numTestCase,1);
for n=1:kernalDim
    DtestVec = DtestVec + A(:,n);
end
% Calculating the w matrix
DA = diag(DtestVec);
tempMat = (inv(DA))*A;
w = sum(tempMat,1);
W = diag(w);
% Calculatin the Abar
Abar = tempMat * (inv(W));
% Calculatin the new affinity matrix
Wu = Abar*Abar';
% Performing the EVD on the affnity matrix
[eigVecWu, eigValWu] = eig(Wu);
% Generating the lambdaR vector
lambdaTrainVec = zeros(kernalDim,1);
count = 1;
for i=1:kernalDim
    lambdaTrainVec(i) = lambdaTrain(count);
    if mod(i,30)==0
        count = count+1;
    end
end
% Estimatin the lambda unseen room
lambdaEstimated = (Abar * lambdaTrainVec)/(sqrt(eigValWu(2,2)));
rtEstimated = -6*log(10)./lambdaEstimated;
% Plotting the results
plot(rt,rtEstimated*10^4,'o');
xlabel('True Reverberation Time');ylabel('Estimated Reverberation Time');