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);