Hey. I got it to work by dissecting the wcoherence function. You need the plotcoherenceperiod, plotPhaseVectors and determinearrowextent functions. Its surprisingly complicated
Below is an example:
[wcohB,wcsB,fB]=wcoherence(normalize(datax),normalize(datay),1/3600,'VoicesPerOctave',32,'PhaseDisplayThreshold',0.7);
N = size(wcohFP,2);
dt=1/3600;
t = 0:dt:N*dt-dt;
FourierFactor = (2*pi)/6;
sigmawav = 1/sqrt(2);
coiScalar = FourierFactor/sigmawav;
coitmp = coiScalar*dt*[1E-5,1:((N+1)/2-1),fliplr((1:(N/2-1))),1E-5];
plotcoherenceperiod(wcohFP,wcsFP,fFP,t,coitmp,nv,mc,'secs')
function plotcoherenceperiod(wcoh,wcs,frequencies,t,coitmp,nv,mc,Units)
period = (1./frequencies);
switch Units
case 'years'
Yticks = 2.^(round(log2(min(period))):round(log2(max(period))));
logYticks = log2(Yticks(:));
YtickLabels = num2str(sprintf('%g\n',Yticks));
case 'days'
Yticks = 2.^(round(log2(min(period))):round(log2(max(period))));
logYticks = log2(Yticks(:));
YtickLabels = num2str(sprintf('%g\n',Yticks));
case 'hrs'
Yticks = 2.^(round(log2(min(period))):round(log2(max(period))));
logYticks = log2(Yticks(:));
YtickLabels = num2str(sprintf('%g\n',Yticks));
case 'mins'
Yticks = 2.^(round(log2(min(period)),1):round(log2(max(period)),1));
logYticks = log2(Yticks(:));
YtickLabels = num2str(sprintf('%g\n',Yticks));
case 'secs'
Yticks = 2.^(round(log2(min(period)),2):round(log2(max(period)),2));
logYticks = log2(Yticks(:));
YtickLabels = num2str(sprintf('%g\n',Yticks));
end
%
AX = newplot;
f = ancestor(AX,'figure');
setappdata(AX,'evstruct',[]);
cla(AX,'reset');
imagesc(t,log2(period),wcoh);
AX.CLim = [0 1];
AX.YLim = log2([min(period), max(period)]);
AX.YTick = logYticks;
AX.YDir = 'normal';
set(AX,'YLim',log2([min(period),max(period)]), ...
'layer','top', ...
'YTick',logYticks, ...
'YTickLabel',YtickLabels, ...
'layer','top')
ylabel([getString(message('Wavelet:wcoherence:Period')) ' (' Units ') ']);
xlabel([getString(message('Wavelet:wcoherence:Time')) ' (' Units ')']);
title(getString(message('Wavelet:wcoherence:CoherenceTitle')));
hold(AX,'on');
hcol = colorbar;
hcol.Label.String = 'Magnitude-Squared Coherence';
plot(AX,t,log2(coitmp),'w--','linewidth',2);
theta = angle(wcs);
theta(wcoh< mc)= NaN;
if all(isnan(theta))
return;
end
% Create mesh grid for phase plot
tspace = ceil(size(theta,2)/40);
pspace = round(2^log2(size(theta,1)/nv/2));
tax = t(1:tspace:size(theta,2));
pax = period(1:pspace:size(theta,1));
plotPhaseVectors(AX,theta,tax,pax,tspace,pspace);
hzoom = zoom(f);
cbzoom = @(~,evd)zoomArrows(evd,theta,tax,pax,tspace,pspace);
cbfig = @(hobject,evd)ResizeFig(hobject,evd,theta,tax,pax,tspace,pspace);
evstruct.sclistener = event.listener(f,'SizeChanged',cbfig);
evstruct.ylimlistener = event.proplistener(AX,AX.findprop('YLim'),...
'PostSet',cbfig);
evstruct.xlimlistener = event.proplistener(AX,AX.findprop('XLim'),...
'PostSet',cbfig);
setappdata(AX,'evstruct',evstruct);
set(hzoom,'ActionPostCallback',cbzoom);
% Set NextPlot property to 'replace'
f.NextPlot = 'replace';
end
function plotPhaseVectors(axhandle,theta,tax,pax,tspace,pspace)
if ~isempty(findobj(axhandle,'type','patch'))
delete(findobj(axhandle, 'type', 'patch'));
end
[tgrid,pgrid]=meshgrid(tax,log2(pax));
theta = theta(1:pspace:size(theta,1),1:tspace:size(theta,2));
idx = find(~any(isnan([tgrid(:) pgrid(:) theta(:)]),2));
tgrid = tgrid(idx);
pgrid = pgrid(idx);
theta = theta(idx);
% Determine extent of phase arrows in plot
[dx,dy] = determinearrowextent(axhandle);
%
% Create the arrow patch object for plotting the phase
arrowpatch = [-1 0 0 1 0 0 -1; 0.1 0.1 0.5 0 -0.5 -0.1 -0.1]';
for ii=numel(tgrid):-1:1
% Multiply each arrow by the rotation matrix for the given theta
rotarrow = arrowpatch*[cos(theta(ii)) sin(theta(ii));-sin(theta(ii)) cos(theta(ii))];
patch(tgrid(ii)+rotarrow(:,1)*dx,pgrid(ii)+rotarrow(:,2)*dy,[0 0 0],...
'edgecolor','none' ,'Parent',axhandle);
end
end
function [dx,dy] = determinearrowextent(axhandle)
% Get the data aspect ratio of the y and x axis
dataaspectratio = get(axhandle,'DataAspectRatio');
axesposition = get(axhandle,'position');
widthheight = axesposition(3:4);
ar = widthheight./dataaspectratio(1:2);
ar(2)=ar(1)/ar(2);
ar(1)=1;
xlim = axhandle.XLim;
dxlim = xlim(2)-xlim(1);
dx=ar(1).*0.02*dxlim;
dy=ar(2).*0.02*dxlim;
end
