ce4f1f0ab7944a839050436b78be2a87cbef1ddd
[epclust.git] / temp / wave.m
1 t = 1:128
2
3 n=length(t);
4
5 %----------default arguments for the wavelet transform-----------
6 Args=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
13 J1 = noctave * nvoice - 1
14
15 nx=size(x,1);
16 ny=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
27 x = runif(128) #runif(100)
28 w1 = wt(cbind(1:128,x), pad = TRUE, dt = NULL, do.sig=FALSE, dj=1/12, J1=51)
29 w2 = 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)
39 sinv=1./(scale');
40 sX=smoothwavelet(sinv(:,ones(1,nx)).*(abs(X).^2),dt,period,Args.Dj,scale);
41 sY=smoothwavelet(sinv(:,ones(1,ny)).*(abs(Y).^2),dt,period,Args.Dj,scale);
42
43 % -------- Cross wavelet -------
44 Wxy=X.*conj(Y);
45
46 % ----------------------- Wavelet coherence ---------------------------------
47 sWxy=smoothwavelet(sinv(:,ones(1,n)).*Wxy,dt,period,Args.Dj,scale);
48
49 %%%%%% SMOOTHWAVELET
50 n=size(wave,2);
51
52 %swave=zeros(size(wave));
53 twave=zeros(size(wave));
54
55 %zero-pad to power of 2... Speeds up fft calcs if n is large
56 npad=2.^ceil(log2(n));
57
58 k = 1:fix(npad/2);
59 k = k.*((2.*pi)/npad);
60 k = [0., k, -k(fix((npad-1)/2):-1:1)];
61
62 k2=k.^2;
63 snorm=scale./dt;
64 for 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);
68 end
69
70 if isreal(wave)
71 twave=real(twave); %-------hack-----------
72 end
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
81 dj0=0.6;
82 dj0steps=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
98 kernel=[mod(dj0steps,1); ones(2 * round(dj0steps)-1,1); ...
99 mod(dj0steps,1)]./(2*round(dj0steps)-1+2*mod(dj0steps,1));
100 swave=conv2(twave,kernel,'same'); %thanks for optimization by Uwe Graichen
101
102
103
104
105
106
107
108
109 Rsq=abs(sWxy).^2./(sX.*sY);
110 return Rsq
111
112 %varargout={Rsq,period,scale,coi,wtcsig,t};
113 %varargout=varargout(1:nargout);