Is there a way to manually plot wavelet(wcoherence) plot with phase arrows using imagesc if we are given wcoh,wcs, f, and coi?

I'm trying to manually create wavelet -wcoherence plot. I don't know how to add arrows into the plot. Please help!

a= wcoh{1,1};
a1=a(52:104,:); % I'm extracting a range of frequencies
b=coi{1,1};
x=[1:12120];
f1=f{1,1};
f2=f1(52:104,:);
imagesc(x,f2,a1),hold on
set(gca,'YDir','normal')
fillout(x,b); % a code from exchange 
box on
set(gca,'BoxStyle','full')
set(gca,'LineWidth',1), hold off

Answer 1

[ForTex]

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