[Octave-cvsupdate] SF.net SVN: octave:[10777] trunk/octave-forge/extra/lssa

[<be...@us...>]

Revision: 10777
          http://octave.svn.sourceforge.net/octave/?rev=10777&view=rev
Author:   benjf5
Date:     2012-07-26 23:07:34 +0000 (Thu, 26 Jul 2012)
Log Message:
-----------
Added lscomplexwavelet, still working on it.

Modified Paths:
--------------
    trunk/octave-forge/extra/lssa/lsrealwavelet.m

Added Paths:
-----------
    trunk/octave-forge/extra/lssa/lscomplexwavelet.m

Added: trunk/octave-forge/extra/lssa/lscomplexwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/lscomplexwavelet.m	                        (rev 0)
+++ trunk/octave-forge/extra/lssa/lscomplexwavelet.m	2012-07-26 23:07:34 UTC (rev 10777)
@@ -0,0 +1,66 @@
+## Copyright (C) 2012 Ben Lewis <be...@gm...>
+##
+## This code is part of GNU Octave.
+##
+## This software is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 2 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+
+
+function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, minimum_window_count, sigma)
+
+## This is a transform based entirely on the simplified complex-valued transform
+## in the Mathias paper, page 10. My problem with the code as it stands is, it
+## doesn't have a good way of determining the window size. Sigma is currently up
+## to the user, and sigma determines the window width (but that might be best.)
+##
+## Currently the code does not apply a time-shift, which needs to be fixed so
+## that it will work correctly over given frequencies.
+##
+## Additionally, each octave up adds one to the number of frequencies.
+
+transform = cell(ncoeff*noctave,1);
+for octave_iter = 1:noctave
+  current_octave = maxfreq * 2 ^ ( - octave_iter );
+  current_window_number = minimum_window_count + noctave - octave_iter;
+  window_step = ( tmax - tmin ) / current_window_number;
+  for coeff_iter = 1: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.
+    o = current_frequency = maxfreq * 2 ^ ( - octave_iter*coeff_iter / ncoeff );
+    current_radius = 1 / ( current_octave * sigma );
+    
+    transform{iter} = zeros(1,current_window_number);
+    win_t = window_min + ( window_step / 2);
+    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.
+      zeta = sum( cubicwgt ((windowed_t - win_t) ./ current_radius) .* exp( - i *
+									      o .* windowed_t ) .* X ) / sum( cubicwgt ((windowed_t - win_t ) ./
+															current_radius) .* exp ( -i * o .* windowed_t ));
+      transform{iter}(iter_window) = zeta;
+      window_min += window_step ;
+    ## I remain hesitant about this value, since it is entirely possible necessary precision will be lost. Should I try to reduce that?
+    endfor
+  endfor
+  
+endfor
+
+endfunction
+

Modified: trunk/octave-forge/extra/lssa/lsrealwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/lsrealwavelet.m	2012-07-26 13:15:46 UTC (rev 10776)
+++ trunk/octave-forge/extra/lssa/lsrealwavelet.m	2012-07-26 23:07:34 UTC (rev 10777)
@@ -10,7 +10,6 @@
 
 function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number )
 
-len_tx = length (T)
 omegamult = 2 ^ ( - 1/ ncoeff )
 omegamult_inv = 1 / omegamult
 
@@ -63,7 +62,7 @@
     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;
 endfor
   
 

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.