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 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
ForTex
am 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 Kommentare
Siehe auch
Kategorien
Mehr zu Wavelet Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!