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 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/283596/image.png)
0 commentaires
Réponses (1)
ForTex
le 11 Fév 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 commentaires
Voir également
Catégories
En savoir plus sur Wavelet Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!