Is there a way to manually plot wavelet(wcoherence) plot with phase arrows using imagesc if we are given wcoh,wcs, f, and coi?
9 views (last 30 days)
Show older comments
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
0 Comments
Answers (1)
ForTex
on 11 Feb 2021
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
0 Comments
See Also
Categories
Find more on Wavelet Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!