MATLAB:相位谱代码
%% Example% This example illustrates how simple it is to do
% continuous wavelet transform (CWT), Cross wavelet transform (XWT)
% and Wavelet Coherence (WTC) plots of your own data.
%
% The time series we will be analyzing are the winter
% Arctic Oscillation index (AO) and
% the maximum sea ice extent in the Baltic (BMI).
%
% http://www.pol.ac.uk/home/research/waveletcoherence/
%% Load the data
% First we load the two time series into the matrices d1 and d2.
clear;
clc;
seriesname={'ET' 'Ta' 'ga'};
% d1=load('jao.txt');
% d2=load('jbaltic.txt');
d1=load('.\2009\ET.txt');
d2=load('.\2009\vpd.txt');
%d3(:,1)=d1(:,1);
%d3(:,2)=d1(:,2)-d2(:,2);
%% Change the pdf.
% The time series of Baltic Sea ice extent is highly bi-modal and we
% therefore transform the timeseries into a series of percentiles. The
% transformed series probably reacts 'more linearly' to climate.
d1(:,2)=boxpdf(d1(:,2));
d2(:,2)=boxpdf(d2(:,2));
%% Continuous wavelet transform (CWT)
% The CWT expands the time series into time
% frequency space.
figure('color',)
%tlim=;
%subplot(1,1,1);
wt(d1)
=wt(d1);
=size(power);
power_mean=[];
for i=1:m
power_mean(i)=mean(power(i,:));
end
% title(seriesname{1});
% set(gca,'xlim',tlim/48,'fontweight','bold','fontsize',14);
% set(gca,'XTick',8:20:148)
% set(gca,'XTicklabel',120:20:260)
% xlabel('Day of year','fontweight','bold','fontsize',16)
% ylabel('Period (d)','fontweight','bold','fontsize',16)
% subplot(2,1,2);
% wt(d2);
% title(seriesname{2});
% set(gca,'xlim',tlim/48,'fontweight','bold','fontsize',11);
%% Cross wavelet transform (XWT)
% The XWT finds regions in time frequency space where
% the time series show high common power.
figure('color',)
xwt(d1,d2)
title(['XWT: ' seriesname{1} '-' seriesname{2} ] )
%% Wavelet coherence (WTC)
% The WTC finds regions in time frequency space where the two
% time series co-vary (but does not necessarily have high power).
figure('color',);
%wtc(d1,d2)
% subplot(2,1,1);
% wtc(d1,d2,'mcc',100); %MCC:MonteCarloCount
% title(['WTC: ' seriesname{1} '-' seriesname{2} ] )
% set(gca,'fontweight','bold','fontsize',11)
% xlabel('Days after sowing (d)','fontweight','bold','fontsize',12)
% subplot(2,1,2);
wtc(d1,d2,'mcc',50); %MCC:MonteCarloCount
% title(['WTC: ' seriesname{1} '-' seriesname{3} ] )
set(gcf,'Position',);
set(gca,'Position', );
set(gca,'fontweight','bold','fontsize',11)
xlabel('DAS(d)','fontweight','bold','fontsize',12)
%imwrite(hh,'aa.tif');
% imsave
% print -dmeta aa
%% Copyright notice
% 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.
页:
[1]