changeset 10299:4f08f44a2705 octave-forge

Added cubicwgt and the wavelet coefficient functions.
author benjf5
date Tue, 22 May 2012 02:02:32 +0000
parents a187606d8d02
children 8ffb34d12d44
files extra/lssa/cubicwgt.m extra/lssa/nucorrcoeff.m extra/lssa/nuwaveletcoeff.m
diffstat 3 files changed, 87 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/lssa/cubicwgt.m	Tue May 22 02:02:32 2012 +0000
@@ -0,0 +1,22 @@
+## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
+##
+## 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/lssa/nucorrcoeff.m	Tue May 22 02:02:32 2012 +0000
@@ -0,0 +1,38 @@
+## Copyright (C) 2012 Benjamin Lewis  <benjf5@gmail.com>
+## 
+## 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/lssa/nuwaveletcoeff.m	Tue May 22 02:02:32 2012 +0000
@@ -0,0 +1,27 @@
+## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
+##
+## 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