function [a,b,c] = corrplot( x,y,wavelet,scales,fs )
% For determining the correlation between two signals at a range of
% periodicities
% Author: Tom D Pering - University of Sheffield
% You are free to use and alter with acknowledgement

% Input Details
% Both signals must be the same sampling rate and the same length
% The signals can be scalars or vectors
% "x" is signal 1
% "y" is signal 2
% "wavelet" is the mother wavelet, example input: 'morl' 
% You can change the mother wavelet, 'gaus8' also works well
% The in-built Matlab Wavelet Toolbox allows visualisation of different
% mother wavelets
% "scales" indicates the maximum range of steps used in the continuous
% wavelet and is a time component (e.g. seconds, hours etc)
% transform (the units are in sampling rate) and the 'window' and 'nfft used in psd welch analysis
% An error will show if the maximum "scales" value set is above the Nyquist
% Criterion e.g. [Nyquist, H., 2002. Certain topics in telegraph transmission 
% theory (Reprinted from Transactions of the A. I. E. E., February, 
% pg 617-644, 1928). Proceedings IEEE. 90 (2), 280-305,
% doi:10.1109/5.989875]
% fs is the sampling rate of the signal in hertz

% Output Details
% "a" is a matrix of correlation coefficients as generated by the
% "corrplot" code
% "b" is the diagonal of output "a" or the "1:1 best-fit line"
% "c" is a matrix of cross-correlation coefficients over the range of lags
% determined by the code
% A number of plots are also auto-generated, see accompanying paper for
% full details

if scales>((length(x)/2))
    error('Scales above Nyquist limit')
end

% Wavelet Transform
cwt1=cwt(x/max(x),1:scales,wavelet);
cwt2=cwt(y/max(y),1:scales,wavelet);

% Shift the data
cwt1=ctranspose(cwt1);
cwt2=ctranspose(cwt2);

% Correlate the data
a=corr(cwt1,cwt2,'type','Spearman');

% Extract the "1:1 best-fit" line
b=diag(a);

% Extract max and min correlation location
[max_corr,loc_max_corr]=max(b)
[min_corr,loc_min_corr]=min(b)
[M1,N1]=ind2sub(size(b),loc_max_corr); % Need M1 for location of largest correlation
[M2,N2]=ind2sub(size(b),loc_min_corr); % Need M2 for location of smallest 

% Individual coefficients at max and min location
wave_coeff1_max=cwt1(:,M1); % Individual normalised coefficeints
wave_coeff1_min=cwt1(:,M2);
wave_coeff2_max=cwt2(:,M1);
wave_coeff2_min=cwt2(:,M2);

% Power spectral densities
[b1,freq1]=pwelch(x/max(x),scales,0,scales,fs);
[b2,freq1]=pwelch(y/max(y),scales,0,scales,fs);

% Xcorr lag plot 
cwt1=ctranspose(cwt1);
cwt2=ctranspose(cwt2);
for ls=1:scales;
    s1=cwt1(ls,:);
    s2=cwt2(ls,:);
    maxlags=scales/2;
    lag_corr=xcorr(s1,s2,maxlags,'coeff');
    c(ls,:)=horzcat(lag_corr);
end
c=ctranspose(c);

% Plot the data
figure1=figure;

% Periodicity Correlation Plot
axes1 = axes('Parent',figure1,'Position',[0.053 0.110 0.459 0.826],...
    'Layer','top');
xlim(axes1,[0.5 scales]);
ylim(axes1,[0.5 scales]);
box(axes1,'on');
hold(axes1,'all');
title('Periodicity Correlation Plot');
imagesc(a,'Parent',axes1,'CDataMapping','scaled');
xlabel('Scales ''x''');
ylabel('Scales ''y''');

% Welch plot - x data
axes2 = axes('Parent',figure1,...
    'Position',[0.645 0.593 0.331 0.341]);
box(axes2,'on');
hold(axes2,'all');
title('PSD Signal ''x''');
xlabel('Frequency (Hz)');
plot(freq1,b1,'Parent',axes2,...
    'DisplayName','b1');

% Welch plot - y data
axes3 = axes('Parent',figure1,...
    'Position',[0.645 0.112 0.334 0.341]);
box(axes3,'on');
hold(axes3,'all');
title('PSD Signal ''y''');
xlabel('Frequency (Hz)');
plot(freq1,b2,'Parent',axes3,...
    'DisplayName','b2 vs freq');

% Colourbar
colorbar('peer',axes1,...
    [0.522 0.108 0.026 0.829]);

% Plot 3d
figure2=figure;
axes1 = axes('Parent',figure2,...
    'Position',[0.094 0.108 0.775 0.815]);
view(axes1,[-37.5 30]);
grid(axes1,'on');
hold(axes1,'all');
surf(a,'Parent',axes1,'LineStyle','none');
xlabel('Scales ''x''');
ylabel('Scales ''y''');
zlabel('Correlation Coefficient');
colorbar('peer',axes1,...
    [0.923 0.112 0.026 0.815]);

% Plot Best-Fit Line and Coefficients;

figure3 = figure('PaperSize',[20.98 29.68]);
axes1 = axes('Parent',figure3,'Position',[0.098 0.709 0.775 0.216],...
    'FontSize',12,...
    'FontName','Times New Roman'); box('on'); hold('all');
    %1
plot(b,'Parent',axes1,'LineWidth',2,'DisplayName','b','Color',[0 0 0]);
xlabel('Scales ''a''','FontSize',12,'FontName','Times New Roman');
ylabel('Correlation - Best Fit Line','FontSize',12,...
    'FontName','Times New Roman');
axes2 = axes('Parent',figure3,'Position',[0.09833 0.4069 0.775 0.2157],...
    'FontSize',12,...
    'FontName','Times New Roman');
box('on'); hold('all');
xlabel('Time [units]','FontSize',12,...
    'FontName','Times New Roman');
ylabel('Wavelet Coefficients','FontSize',12,'FontName','Times New Roman');
    %2
plot1 = plot([wave_coeff1_max wave_coeff2_max],'Parent',axes2,'LineWidth',2);
set(plot1(1),'DisplayName','Max ''cwt1''','Color',[0 0 0]);
set(plot1(2),'LineStyle','--','Color',[0.502 0.502 0.502],...
    'DisplayName','Max ''cwt2''');
axes3 = axes('Parent',figure3,'Position',[0.09917 0.11 0.775 0.2157],...
    'FontSize',12,...
    'FontName','Times New Roman');
box('on'); hold('all');
    %3
plot2 = plot([wave_coeff1_min wave_coeff2_min],'Parent',axes3,'LineWidth',2);
set(plot2(1),'DisplayName','Min ''cwt1''','Color',[0 0 0]);
set(plot2(2),'LineStyle','--','Color',[0.502 0.502 0.502],...
    'DisplayName','Min ''cwt2''');
xlabel('Time [units]','FontSize',12,...
    'FontName','Times New Roman');
ylabel('Wavelet Coefficients','FontSize',12,'FontName','Times New Roman');

legend1 = legend(axes2,'show');
set(legend1,'Position',[0.8865 0.5572 0.1092 0.06648]);
legend2 = legend(axes3,'show');
set(legend2,'Position',[0.889 0.2585 0.1075 0.06648]);

% Plot xcorr lag plot
figure4=figure;
axes1 = axes('Parent',figure4,'YDir','normal',...
    'Position',[0.131 0.11 0.677 0.815],...
    'Layer','top',...
    'FontSize',12,...
    'FontName','Times New Roman');
Xticks=(1:1:scales); Yticks=(-scales/2:1:scales/2);
xlim(axes1,[0.5 scales]); ylim(axes1,[(-scales/2) (scales/2)]);
box(axes1,'on');
hold(axes1,'all');
title('Xcorr Lag Plot');
image(Xticks,Yticks,c,'Parent',axes1,'CDataMapping','scaled'); 
xlabel('Scales','FontSize',12,'FontName','Times New Roman');
ylabel('Lags','FontSize',12,'FontName','Times New Roman');
colorbar('peer',axes1,...
    [0.860 0.108 0.0378 0.816],...
    'FontSize',12,...
    'FontName','Times New Roman');
end


