'update'
[epclust.git] / temp / wave.m
CommitLineData
b7cd987d
BA
1t = 1:128
2
3n=length(t);
4
5%----------default arguments for the wavelet transform-----------
6Args=struct('Pad',1,... % pad the time series with zeroes (recommended)
7 'Dj',1/12, ... % this will do 12 sub-octaves per octave
8 'S0',2*dt,... % this says start at a scale of 2 years
9 'J1',[],...
10 'Mother','Morlet', ...
11 'MaxScale',[],... %a more simple way to specify J1
12
13J1 = noctave * nvoice - 1
14
15nx=size(x,1);
16ny=size(y,1);
17
18[X,period,scale] = wavelet(x,Args.Dj,Args.S0,Args.J1);
19[Y,period,scale] = wavelet(y,Args.Dj,Args.S0,Args.J1);
20
21
22
23
24
25## pad with 0 ?
26## compare Rwave and biwavelets
27x = runif(128) #runif(100)
28w1 = wt(cbind(1:128,x), pad = TRUE, dt = NULL, do.sig=FALSE, dj=1/12, J1=51)
29w2 = Rwave::cwt(x, 13, 4, w0=2*pi, twoD=TRUE, plot=FALSE)
30
31
32
33
34
35
36
37
38%Smooth X and Y before truncating! (minimize coi)
39sinv=1./(scale');
40sX=smoothwavelet(sinv(:,ones(1,nx)).*(abs(X).^2),dt,period,Args.Dj,scale);
41sY=smoothwavelet(sinv(:,ones(1,ny)).*(abs(Y).^2),dt,period,Args.Dj,scale);
42
43% -------- Cross wavelet -------
44Wxy=X.*conj(Y);
45
46% ----------------------- Wavelet coherence ---------------------------------
47sWxy=smoothwavelet(sinv(:,ones(1,n)).*Wxy,dt,period,Args.Dj,scale);
48
49%%%%%% SMOOTHWAVELET
50n=size(wave,2);
51
52%swave=zeros(size(wave));
53twave=zeros(size(wave));
54
55%zero-pad to power of 2... Speeds up fft calcs if n is large
56npad=2.^ceil(log2(n));
57
58k = 1:fix(npad/2);
59k = k.*((2.*pi)/npad);
60k = [0., k, -k(fix((npad-1)/2):-1:1)];
61
62k2=k.^2;
63snorm=scale./dt;
64for ii=1:size(wave,1)
65 F=exp(-.5*(snorm(ii)^2)*k2); %Thanks to Bing Si for finding a bug here.
66 smooth=ifft(F.*fft(wave(ii,:),npad));
67 twave(ii,:)=smooth(1:n);
68end
69
70if isreal(wave)
71 twave=real(twave); %-------hack-----------
72end
73
74%scale smoothing (boxcar with width of .6)
75
76%
77% TODO: optimize. Because this is done many times in the monte carlo run.
78%
79
80
81dj0=0.6;
82dj0steps=dj0/(dj*2);
83% for ii=1:size(twave,1)
84% number=0;
85% for l=1:size(twave,1);
86% if ((abs(ii-l)+.5)<=dj0steps)
87% number=number+1;
88% swave(ii,:)=swave(ii,:)+twave(l,:);
89% elseif ((abs(ii-l)+.5)<=(dj0steps+1))
90% fraction=mod(dj0steps,1);
91% number=number+fraction;
92% swave(ii,:)=swave(ii,:)+twave(l,:)*fraction;
93% end
94% end
95% swave(ii,:)=swave(ii,:)/number;
96% end
97
98kernel=[mod(dj0steps,1); ones(2 * round(dj0steps)-1,1); ...
99 mod(dj0steps,1)]./(2*round(dj0steps)-1+2*mod(dj0steps,1));
100swave=conv2(twave,kernel,'same'); %thanks for optimization by Uwe Graichen
101
102
103
104
105
106
107
108
109Rsq=abs(sWxy).^2./(sX.*sY);
110return Rsq
111
112%varargout={Rsq,period,scale,coi,wtcsig,t};
113%varargout=varargout(1:nargout);