From: <be...@us...> - 2012-07-25 15:06:01
|
Revision: 10774 http://octave.svn.sourceforge.net/octave/?rev=10774&view=rev Author: benjf5 Date: 2012-07-25 15:05:47 +0000 (Wed, 25 Jul 2012) Log Message: ----------- lsrealwavelet written; not fully tested, but performs as expected thus far. Modified Paths: -------------- trunk/octave-forge/extra/lssa/lsrealwavelet.m Modified: trunk/octave-forge/extra/lssa/lsrealwavelet.m =================================================================== --- trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-25 04:48:43 UTC (rev 10773) +++ trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-25 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_subdiv-1) -# 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 (T-win_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) && ((t-T(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 (T-win_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. |