function varargout = PCA(PCA_sel,M,t,NumTrig)

res = size(PCA_sel,1);

eg = 4;
sel = 5;    %Number of principal components for reconstruction

psi = mean(PCA_sel,2);  %psi = average signal
figure(1)
plot(t,psi,'r-')              %psi = displays average face
title('Average of all signals')
xlabel('Time (ms)'); ylabel('Voltage (mV');
hold on

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

% Enter a name for each set of averaged signals in the lines below, 
% eg. 'psi_Cannabinoid'. Comment out all the other lines by typing "%" in
% front. Then click 'Edit' , 'Find and Replace' and change all the names
% accordingly (eg. replace 'basal' with 'Cannabinoid') throughout PCA.m by
% clicking 'Replace All'

psi_basal = mean(PCA_sel(:,1:NumTrig),2);  %average basal signal
psi_drug1 = mean(PCA_sel(:,NumTrig+1:2*NumTrig),2);  %average drug1 signal
%psi_drug2 = mean(PCA_sel(:,2*NumTrig+1:3*NumTrig),2);  %average drug2 signal
% psi_something-interesting1 = mean(PCA_sel(:,3*NumTrig+1:4*NumTrig),2);  %average drug signal
% psi_something-interesting2 = mean(PCA_sel(:,4*NumTrig+1:5*NumTrig),2);  %average drug signal


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


plot(t,psi_basal,'k:') 
hold on
plot(t,psi_drug1,'b--') 

legend ('Average of all responses', 'Average of first set of responses', 'Average of second set of responses')


psi2 = repmat(psi,1,M);    %stacks psi (avg signal) 
A = PCA_sel - psi2;    %A, the difference between each 
                              %signal and the average signal, with phi_i arranged in columns 
    
L = A'*A;        %L is equivalent to the covariance matrix
[V,K]=eig(L);   %Calculates the eigenvalues and eigenvectors of L

prenormU = A*V;        %Calculates the set of eigensignals, arranged in M columns


%Normalise the eigensignals
for jj = 1:size(prenormU,2)
U(:,jj) = prenormU(:,jj)./sqrt(sum(prenormU(:,jj)).^2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%selects eigensignals with highest eigenvalues
eigvals = diag(K);      %selects all eigenvalues    
sorted = sort(eigvals,1,'descend'); %sorts eigenvalues into descending order
lim = sorted(sel);                  %selects the values of the smallest eigenenvalue to not be discarded
trunc_eig = eigvals.*(eigvals>=lim);  %truncates array of eigenvalues to discard unwanted small eigenvalues
selU_1 = zeros(res,1);              %creates a column vector of zeros, length = no. of eigenvalues kept

for ss = 1:M
    if trunc_eig(ss) ~= 0
    selU_1 = [selU_1 U(:,ss)]; %selects eigensignals corresponding with highest 'sel' eigenvalues
    end
end
selU = fliplr(selU_1(:,2:sel+1));   %selected eigensignals in columns
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure(2)
for zz = 1:sel
subplot(sel,1,zz)             %Displays eigensignals
plot(t,selU(:,zz))
title(strcat('PC', num2str(zz), '   Percentage variance =  ',num2str(sorted(zz)/sum(sorted)*100), ' %'))
end
xlabel(strcat('Total variance in first', num2str(sel), 'components', '= ',num2str(sum(sorted(1:zz)/sum(sorted)*100)), ' %'))


omegas = selU'*A;   %caluclates a set of weights (projection of eigensignals onto each mean-adjusted signal)
% x = omegas(1);
% y = omegas(2);
PCX = 1;
PCY = 2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plotting the data in Principal Component space

princ_comps_basal = omegas(PCX:PCY-PCX:PCY,1:NumTrig);  %Selects the 3 principal components for the first set of data entered
PC1_basal = princ_comps_basal(1,:);
PC2_basal = princ_comps_basal(2,:);

figure(3)
plot(PC1_basal,PC2_basal,'k.')  %Plots PC1 vs. PC2 vs. PC3 for the first set of data entered
hold on

princ_comps_drug1 = omegas(PCX:PCY-PCX:PCY,NumTrig+1:2*NumTrig); %Selects the 3 principal components for the second set of data entered
PC1_drug1 = princ_comps_drug1(1,:);
PC2_drug1 = princ_comps_drug1(2,:);

plot(PC1_drug1,PC2_drug1, 'bx')


grid on
xlabel(strcat('PC',num2str(PCX))); ylabel(strcat('PC',num2str(PCY)))
title('Principal Component comparison')
legend ('First set of responses', 'Second set of responses')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S T A T S
SdeviationPCs_basal = std(PC1_basal+PC2_basal)
SdeviationPCs_drug1 = std(PC1_drug1+PC2_drug1)

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


 reconsignal = selU*omegas;   %reconstructs signals from eigensignals
% 
 reconstructed =(reconsignal(:,eg)+psi); %add psi to bring back any d.c. offset

 figure(4)
 plot(reconstructed)
 title(strcat('Reconstructed example response', num2str(eg)))


figure(5)


reconstructed_basal = mean(reconsignal(:,1:NumTrig),2)+psi;  %finds the PCA reconstruction for basal mean

reconstructed_drug1 = mean(reconsignal(:,NumTrig+1:2*NumTrig),2)+psi;   %finds the PCA reconstruction for drug1 mean
plot(t,reconstructed_basal, 'k:')
hold on
plot(t,reconstructed_drug1, 'b--')
xlabel('Time (s)'); ylabel('Voltage (mV');
title(strcat('Comparison of means reconstructed using ', num2str(sel), ' principal components'))
legend ('Reconstructed first set of responses', 'Reconstructed second set of responses')