Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

To resolve issues starting MATLAB on Mac OS X 10.10 (Yosemite) visit: http://www.mathworks.com/matlabcentral/answers/159016

How Can I change the colorbar value.

Asked by Besma on 27 Jul 2013
This is my Matlab code of wavelet coherence:

function varargout=wtc(x,y,varargin) %% Wavelet coherence % % USAGE: [Rsq,period,scale,coi,sig95]=wtc(x,y,[,settings]) % % % Settings: Pad: pad the time series with zeros? % . Dj: Octaves per scale (default: '1/12') % . S0: Minimum scale % . J1: Total number of scales % . Mother: Mother wavelet (default 'morlet') % . MaxScale: An easier way of specifying J1 % . MakeFigure: Make a figure or simply return the output. % . BlackandWhite: Create black and white figures % . AR1: the ar1 coefficients of the series % . (default='auto' using a naive ar1 estimator. See ar1nv.m) % . MonteCarloCount: Number of surrogate data sets in the significance calculation. (default=300) % . ArrowDensity (default: [30 30]) % . ArrowSize (default: 1) % . ArrowHeadSize (default: 1) % % Settings can also be specified using abbreviations. e.g. ms=MaxScale. % For detailed help on some parameters type help wavelet. % % Example: % t=1:200; % wtc(sin(t),sin(t.*cos(t*.01)),'ms',16) % % Please acknowledge the use of this software in any publications: % "Crosswavelet and wavelet coherence software were provided by % A. Grinsted." % % (C) Aslak Grinsted 2002-2004 % % http://www.pol.ac.uk/home/research/waveletcoherence/

% ------------------------------------------------------------------------- % Copyright (C) 2002-2004, Aslak Grinsted % This software may be used, copied, or redistributed as long as it is not % sold and this copyright notice is reproduced on each copy made. This % routine is provided as is without any express or implied warranties % whatsoever.

% ------validate and reformat timeseries. [x,dt]=formatts(x); [y,dty]=formatts(y); if dt~=dty error('timestep must be equal between time series') end t=(max(x(1,1),y(1,1)):dt:min(x(end,1),y(end,1)))'; %common time period if length(t)<4 error('The two time series must overlap.') end

n=length(t);

%----------default arguments for the wavelet transform----------- Args=struct('Pad',1,... % pad the time series with zeroes (recommended) 'Dj',1/12, ... % this will do 12 sub-octaves per octave 'S0',2*dt,... % this says start at a scale of 2 years 'J1',[],... 'Mother','Morlet', ... 'MaxScale',[],... %a more simple way to specify J1 'MakeFigure',(nargout==0),... 'MonteCarloCount',10000,... 'BlackandWhite',1,... 'AR1',[0 0],... 'ArrowDensity',[30 30],... 'ArrowSize',0,... 'ArrowHeadSize',0); Args=parseArgs(varargin,Args,{'BlackandWhite'}); if isempty(Args.J1) if isempty(Args.MaxScale) Args.MaxScale=(n*.17)*2*dt; %auto maxscale end Args.J1=round(log2(Args.MaxScale/Args.S0)/Args.Dj); end

ad=mean(Args.ArrowDensity); Args.ArrowSize=Args.ArrowSize*30*.03/ad; Args.ArrowHeadSize=Args.ArrowHeadSize*Args.ArrowSize*220;

if ~strcmpi(Args.Mother,'morlet') warning('Smoothing operator is designed for morlet wavelet.') end

if strcmpi(Args.AR1,'auto') Args.AR1=[ar1nv(x(:,2)) ar1nv(y(:,2))]; if any(isnan(Args.AR1)) error('Automatic AR1 estimation failed. Specify it manually (use arcov or arburg).') end end

nx=size(x,1); sigmax=std(x(:,2));

ny=size(y,1); sigmay=std(y(:,2));

%-----------:::::::::::::--------- ANALYZE ----------::::::::::::------------

[X,period,scale,coix] = wavelet(x(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother); [Y,period,scale,coiy] = wavelet(y(:,2),dt,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);

%Smooth X and Y before truncating! (minimize coi) sinv=1./(scale');

sX=smoothwavelet(sinv(:,ones(1,nx)).*(abs(X).^2),dt,period,Args.Dj,scale); sY=smoothwavelet(sinv(:,ones(1,ny)).*(abs(Y).^2),dt,period,Args.Dj,scale);

% truncate X,Y to common time interval (this is first done here so that the coi is minimized) dte=dt*.01; %to cricumvent round off errors with fractional timesteps idx=find((x(:,1)>=(t(1)-dte))&(x(:,1)<=(t(end)+dte))); X=X(:,idx); sX=sX(:,idx); coix=coix(idx);

idx=find((y(:,1)>=(t(1))-dte)&(y(:,1)<=(t(end)+dte))); Y=Y(:,idx); sY=sY(:,idx); coiy=coiy(idx);

coi=min(coix,coiy);

% -------- Cross wavelet ------- Wxy=X.*conj(Y);

% ----------------------- Wavelet coherence --------------------------------- sWxy=smoothwavelet(sinv(:,ones(1,n)).*Wxy,dt,period,Args.Dj,scale); Rsq=abs(sWxy).^2./(sX.*sY);

if (nargout>0)||(Args.MakeFigure) wtcsig=wtcsignif(Args.MonteCarloCount,Args.AR1,dt,length(t)*2,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother,.6); wtcsig=(wtcsig(:,2))*(ones(1,n)); wtcsig=Rsq./wtcsig; end

if Args.MakeFigure

    Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
    if Args.BlackandWhite
      levels = [0 .5 .6 .7 .8 .9 1];
   %     levels = [0 .7 .75 .8 .85 .9 .95 1];
   %     levels = [-1 -.7 -.5 .0 .5 .9 .95 1];
        [cout,H]=safecontourf(t,log2(period),Rsq,levels);
        colorbarf(cout,H)
        cmap=[0 1;.5 .9;.6 .8;.7 .7;.8 .6;.9 .5;1 .4];
        cmap=interp1(cmap(:,1),cmap(:,2),(0:.1:1)');
%        cmap=[0 1;.7 .9;.75 .8;.8 .7;.85 .6;.9 .5;.95 .4;1 .3];
%        cmap=interp1(cmap(:,1),cmap(:,2),(0:.05:1)');
        cmap=cmap(:,[1 1 1]);
        colormap(cmap)
        set(gca,'YLim',log2([min(period),max(period)]), ...
            'YDir','reverse', 'layer','top', ...
            'YTick',log2(Yticks(:)), ...
            'YTickLabel',num2str(Yticks'), ...
            'layer','top')
        ylabel('Frequency (in years)')
        xlabel('Time')
        hold on
        %phase plot
        aWxy=angle(Wxy);
        lowRidx=find(Rsq<.5); %remove phase indication where Rsq is low
        aaa=aWxy;
        aaa(lowRidx)=NaN;
        %[xx,yy]=meshgrid(t(1:5:end),log2(period));
        phs_dt=round(length(t)/30); tidx=max(floor(phs_dt/2),1):phs_dt:length(t);
        phs_dp=round(length(period)/30); pidx=max(floor(phs_dp/2),1):phs_dp:length(period);
%        phaseplot(t(tidx),log2(period(pidx)),aaa(pidx,tidx),Args.ArrowSize,Args.ArrowHeadSize);
        if ~all(isnan(wtcsig))
            [c,h] = contour(t,log2(period),wtcsig,[1 1],'k');
            set(h,'linewidth',3)                                    %
        end
        %suptitle([sTitle ' coherence']);
%        plot(t,log2(coi),'k','linewidth',2,'LineStyle','--')         %
        hold off
    else
        H=imagesc(t,log2(period),Rsq);
        %[c,H]=safecontourf(t,log2(period),Rsq,[0:.05:1]);
        %set(H,'linestyle','none')
        set(gca,'clim',[0 1])
        HCB=safecolorbar;
        set(gca,'YLim',log2([min(period),max(period)]), ...
            'YDir','reverse', 'layer','top', ...
            'YTick',log2(Yticks(:)), ...
            'YTickLabel',num2str(Yticks'), ...
            'layer','top')
        ylabel('')
        hold on
        %phase plot
        aWxy=angle(Wxy);
        lowRidx=find(Rsq<.5); %remove phase indication where Rsq is low
        aaa=aWxy;
        aaa(lowRidx)=NaN;
        %[xx,yy]=meshgrid(t(1:5:end),log2(period));
        phs_dt=round(length(t)/30); tidx=max(floor(phs_dt/2),1):phs_dt:length(t);
        phs_dp=round(length(period)/30); pidx=max(floor(phs_dp/2),1):phs_dp:length(period);
        phaseplot(t(tidx),log2(period(pidx)),aaa(pidx,tidx),Args.ArrowSize,Args.ArrowHeadSize);
        if ~all(isnan(wtcsig))
            [c,h] = contour(t,log2(period),wtcsig,[1 1],'k');
            set(h,'linewidth',2)
        end
        %suptitle([sTitle ' coherence']);
        tt=[t([1 1])-dt*.5;t;t([end end])+dt*.5];
        hcoi=fill(tt,log2([period([end 1]) coi period([1 end])]),'w');
        set(hcoi,'alphadatamapping','direct','facealpha',.5)
        hold off
    end
end

varargout={Rsq,period,scale,coi,wtcsig}; varargout=varargout(1:nargout);

function [cout,H]=safecontourf(varargin) vv=sscanf(version,'%i.'); if (version('-release')<14)|(vv(1)<7) [cout,H]=contourf(varargin{:}); else [cout,H]=contourf('v6',varargin{:}); end

function hcb=safecolorbar(varargin) vv=sscanf(version,'%i.');

if (version('-release')<14)|(vv(1)<7) hcb=colorbar(varargin{:}); else hcb=colorbar('v6',varargin{:}); end

Besma

Products

0 Answers

Contact us