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

SourceForge Headquarters 225 Broadway Suite 1600 San Diego, CA 92101 +1 (858) 422-6466

Revision: 10485
          http://octave.svn.sourceforge.net/octave/?rev=10485&view=rev
Author:   benjf5
Date:     2012-05-22 02:02:32 +0000 (Tue, 22 May 2012)
Log Message:
-----------
Added cubicwgt and the wavelet coefficient functions.

Added Paths:
-----------
    trunk/octave-forge/extra/lssa/cubicwgt.m
    trunk/octave-forge/extra/lssa/nucorrcoeff.m
    trunk/octave-forge/extra/lssa/nuwaveletcoeff.m

Added: trunk/octave-forge/extra/lssa/cubicwgt.m
===================================================================

--- trunk/octave-forge/extra/lssa/cubicwgt.m	                        (rev 0)
+++ trunk/octave-forge/extra/lssa/cubicwgt.m	2012-05-22 02:02:32 UTC (rev 10485)
@@ -0,0 +1,22 @@
+## Copyright (C) 2012 Benjamin Lewis <be...@gm...>
+##
+## This program 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 3 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/>.
+
+## This function implements the windowing function on page 10 of the doc.
+## if t is in [-1,1] then the windowed term is a = 1 + ( |t|^2 * ( 2|t| - 3 )
+## else the windowed term is 0.
+function a = cubicwgt(t) ## where t is the set of time values
+  a = abs(t);
+  a = ifelse( ( a < 1 ), 1 + ( ( a .^ 2 ) .* ( 2 .* a - 3 ) ), a = 0);
+endfunction

Added: trunk/octave-forge/extra/lssa/nucorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/nucorrcoeff.m	                        (rev 0)
+++ trunk/octave-forge/extra/lssa/nucorrcoeff.m	2012-05-22 02:02:32 UTC (rev 10485)
@@ -0,0 +1,38 @@
+## Copyright (C) 2012 Benjamin Lewis  <be...@gm...>
+## 
+## This program 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 3 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/>.
+
+## nucorrcoeff, computes a coefficient of the wavelet correlation of two time series
+
+function coeff = nucorrcoeff(X1, Y1, X2, Y2, t, o, wgt = cubicwgt, wgtrad = 1)
+  so = 0.05 * o;
+  ## The first solution that comes to mind is admittedly slightly ugly and has a data footprint of O(2n)
+  ## but it is vectorised. 
+  mask = ( abs( X1 - t ) * so ) < wgtrad;
+  mask = mask .* [ 1 : length(mask) ];
+  rx1 = X1(mask); ## I've kept the variable names from the R function here
+  ry1 = Y1(mask); ## Needs to have a noisy error if length(Y1) != length(X1) -- add this!
+  mask = ( abs( X2 - t ) * so ) < wgtrad;
+  mask = mask .* [ 1 : length(mask) ];
+  rx2 = X2(mask);
+  ry2 = Y2(mask);
+  ## I've used the same mask for all of these as it's an otherwise unimportant variable ... can this leak memory?
+  length(rx1) ##printing this length is probably used as a warning if 0 is returned; I inculded it
+  ## in particular to maintain an exact duplicate of the R function.
+  s = sum( wgt( ( rx1 - t ) .* so ) ) * sum( wgt( ( rx2 - t ) .* so );
+  coeff = ifelse( s != 0 , ( sum( wgt( ( rx1 - t ) .* so ) .* exp( i .* o .* rx1 ) .* ry1 )
+		 * sum( wgt( ( rx2 - t ) .* so ) .* exp( i .* o .* rx2 ) .* ry2 ) ) / s, 0 );
+
+
+endfunction
\ No newline at end of file

Added: trunk/octave-forge/extra/lssa/nuwaveletcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/nuwaveletcoeff.m	                        (rev 0)
+++ trunk/octave-forge/extra/lssa/nuwaveletcoeff.m	2012-05-22 02:02:32 UTC (rev 10485)
@@ -0,0 +1,27 @@
+## Copyright (C) 2012 Benjamin Lewis <be...@gm...>
+##
+## This program 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 coeff = nuwaveletcoeff( X , Y , t , o , wgt = cubicwgt , wgtrad = 1 )
+  so = 0.05 .* o;
+  mask = abs( X - t ) * so < wgtrad;
+  mask = mask .* [ 1 : length(mask) ];
+  rx = X(mask);
+  ## This was the fastest way to extract a matching subset that I could think of, but it has a complexity O(2n).
+  ry = Y(mask);
+  ## Going by the R code, this can use the same mask.
+  s = sum( wgt( ( X - t ) .* so ) );
+  coeff = ifelse( s != 0 , sum( wgt( ( rx - t ) .* so) .* exp( i .* o .* ( rx - t ) ) .* ry ) ./ s , 0 );
+  
+endfunction
\ No newline at end of file

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





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

A collection of packages providing extra functionality for GNU Octave

octave-cvsupdate