From: <benjf5@us...>  20120725 15:06:01

Revision: 10774 http://octave.svn.sourceforge.net/octave/?rev=10774&view=rev Author: benjf5 Date: 20120725 15:05:47 +0000 (Wed, 25 Jul 2012) Log Message:  lsrealwavelet written; not fully tested, but performs as expected thus far. Modified Paths:  trunk/octaveforge/extra/lssa/lsrealwavelet.m Modified: trunk/octaveforge/extra/lssa/lsrealwavelet.m ===================================================================  trunk/octaveforge/extra/lssa/lsrealwavelet.m 20120725 04:48:43 UTC (rev 10773) +++ trunk/octaveforge/extra/lssa/lsrealwavelet.m 20120725 15:05:47 UTC (rev 10774) @@ 8,18 +8,23 @@ ## @end deffun function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number, sigma=0.1 ) +function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number ) delta_t = (t_max  t_min)/(t_subdiv1) # fails if you're trying to use less than 2 windows. Why use the wavelet tool at that point?  len_tx = length (T) omegamult = 2 ^ (  1/ ncoeff ) omegamult_inv = 1 / omegamult winradius = pi / ( sigma * maxfreq ) current_radius = winradius; +minimum_window_width = ( t_max  t_min ) / minimum_window_number; +minimum_window_radius = minimum_window_width / 2; +sigma = maxfreq * 2 ^ ( noctave ) / minimum_window_radius; +## sigma needs to be such that  t _ k  t  = minimum_window_radius implies +## that sigma * maxfrequency * 2 ^ (  noctave ) * minimum_window_radius = 1 + +## for a specific other frequency, sigma * frequency * window_radius = 1 means +## window_radius = 1 / ( frequency * sigma ) + + o = maxfreq; # zeta _ ( t , omega ) = sum(w(sigma omega (t_k  t )e^(i omega (t_k  t))xi_k) @@ 33,57 +38,34 @@ transform = cell(noctave*ncoeff,1); ## in this, win_t is the centre of the window in question window_min = t_min; ## Although that will vary depending on the window. This is just an implementation for the first window. win_t = window_min + window_radius;  windowed_t = ((abs (Twin_t) < window_radius) .* T); zeta = sum( cubicwgt ((windowed_t  win_t) ./ window_radius) .* exp(  i * o .* windowed_t ) .* X ) / sum ( cubicwgt( windowed_t ./ window_radius ) ); iota = sum( cubicwgt ((windowed_t  win_t) ./ window_radius) .* exp(  i * 2 * o .* windowed_t) .* X ) / sum ( cubicwgt( windowed_t .* 2 * o ) ) ## this will of course involve an additional large memory allocation, at least in the short term, ## but it's faster to do the operation once on the time series, then multiply it by the data series.   +for iter = 1:(noctave*ncoeff) + ## in this, win_t is the centre of the window in question + window_min = t_min; + ## Although that will vary depending on the window. This is just an + ## implementation for the first window. + current_frequency = maxfreq * 2 ^ (  iter / ncoeff ); + current_radius = 1 / ( current_frequency * sigma ); + current_window_number = ceil( ( t_max  t_min ) / current_radius); for iter_freq = 1:( noctave * ncoeff )  t = t_min;  tx_point = 1;  for iter_elem = 1:t_subdiv  if tx_point > len_tx  break  endif  while ( tx_point+1 <= len_tx) && ((tT(tx_point+1))<winradius)  tx_point++;  endwhile  zeta = 0;  iota = 0;  iota0 = 0;  count = 0;  ot = o * ( T(tx_point)  t);  sot = sigma * ot;  while ( tx_point <= len_tx) && ( (sot = sigma * (ot = o*(T(tx_point)t))) < pi )  w = 0.5 * ( 1 + cos (sot));  iota0 += w;  zeta += w * ( cos (ot) + sin (ot) * i ) * X(tx_point);  ot *= 2;  iota += w * ( cos(ot) + sin (ot) * i );  count++;  tx_point++;  endwhile  if count > 0  result = 2 * ( conj (zeta) * iota0 + zeta * conj (iota) ) / ( ( count + iota0 ) ** 2  real(iota) ** 2  imag(iota) ** 2 )  printf("%d\n",result);  transform(iter_freq,iter_elem) = result;  else  transform(iter_freq,iter_elem) = 0;  endif  t += delta_t + transform{iter} = zeros(1,current_window_number); + win_t = window_min + current_radius; + for iter_window = 1:current_window_number + ## the beautiful part of this code is that if parts end up falling outside the + ## vector, it won't matter (although it's wasted computations.) + ## I may add a trap to avoid too many wasted cycles. + windowed_t = ((abs (Twin_t) < current_radius) .* T); + ## this will of course involve an additional large memory allocation, at least in the short term, + ## but it's faster to do the operation once on the time series, then multiply it by the data series. + iota0 = sum ( cubicwgt (windowed_t ./ current_radius ) ); + zeta = sum( cubicwgt ((windowed_t  win_t) ./ current_radius) .* exp(  i * o .* windowed_t ) .* X ) / iota0; + iota = sum( cubicwgt ((windowed_t  win_t) ./ current_radius) .* exp(  i * 2 * o .* windowed_t) .* X ) / sum ( cubicwgt( windowed_t .* 2 * o ) ); + transform{iter}(iter_window) = 2 * ( conj(zeta) * iota0 + zeta * conj(iota) ) / ( ( length ( find (windowed_t)) + iota0 ) ^ 2  real(iota) ^ 2  imag(iota) ^ 2 ); + window_min += 2 * current_radius; + ## I remain hesitant about this value, since it is entirely possible necessary precision will be lost. Should I try to reduce that? endfor  o *= omegamult  winradius *= omegamult_inv  iter_freq + endfor + endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. 