MATLAB:相位谱代码

admin

[mw_shl_code=applescript,true]%% 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',[1 1 1])
%tlim=[min(d1(1,1),d2(1,1)) max(d1(end,1),d2(end,1))];
%subplot(1,1,1);
wt(d1)
[power,period]=wt(d1);
[m,n]=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',[1 1 1])
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',[1 1 1]);
%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',[200 200 500 480]);
set(gca,'Position', [0.2 0.2 0.7 0.7]);
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.[/mw_shl_code]