Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10521] trunk/octave-forge/extra/lssa
|
From: <be...@us...> - 2012-05-25 14:49:59
|
Revision: 10521
http://octave.svn.sourceforge.net/octave/?rev=10521&view=rev
Author: benjf5
Date: 2012-05-25 14:49:45 +0000 (Fri, 25 May 2012)
Log Message:
-----------
Fixed lssa filenames, blatant error in lsreal.
Added Paths:
-----------
trunk/octave-forge/extra/lssa/lscorrcoeff.m
trunk/octave-forge/extra/lssa/lsreal.m
trunk/octave-forge/extra/lssa/lswaveletcoeff.m
Removed Paths:
-------------
trunk/octave-forge/extra/lssa/nucorrcoeff.m
trunk/octave-forge/extra/lssa/nureal.m
trunk/octave-forge/extra/lssa/nuwaveletcoeff.m
Copied: trunk/octave-forge/extra/lssa/lscorrcoeff.m (from rev 10520, trunk/octave-forge/extra/lssa/nucorrcoeff.m)
===================================================================
--- trunk/octave-forge/extra/lssa/lscorrcoeff.m (rev 0)
+++ trunk/octave-forge/extra/lssa/lscorrcoeff.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -0,0 +1,67 @@
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq)
+## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window)
+## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window, winradius)
+##
+## Return the coefficient of the wavelet correlation of time
+## series (@var{abc1}, @var{ord1}) and (@var{abc2}, @var{ord2}).
+## @var{window} is used to apply a windowing function, its
+## default is cubicwgt if left blank, and its radius is 1,
+## as defined as the default for @var{winradius}.
+##
+## @end deftypefn
+
+## Demo with sin, cos as Nir suggested.
+%!demo
+%! x = 1:10;
+%! y = sin(x);
+%! z = cos(x);
+%! a = nucorrcoeff(x,y,x,z,0.5,0.9)
+%! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9
+
+
+## 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;
+ ## This code can only, as of currently, work on vectors; I haven't figured out a way to make it work on a matrix.
+ if( ( ndims(x1) == 2 ) && ! ( rows(x1) == 1 ) )
+ x1 = reshape(x1,1,length(x1));
+ y1 = reshape(y1,1,length(y1));
+ x2 = reshape(x2,1,length(x2));
+ y2 = reshape(y2,1,length(y2));
+ endif
+ ## 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?
+ lnength(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
Copied: trunk/octave-forge/extra/lssa/lsreal.m (from rev 10520, trunk/octave-forge/extra/lssa/nureal.m)
===================================================================
--- trunk/octave-forge/extra/lssa/lsreal.m (rev 0)
+++ trunk/octave-forge/extra/lssa/lsreal.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -0,0 +1,60 @@
+## Copyright (C) 2012 Benjamin Lewis
+##
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {transform =} nureal ( time, mag, maxfreq, numcoeff, numoctaves)
+##
+## Return the real least-squares transform of the time series
+## defined, based on the maximal frequency @var{maxfreq}, the
+## number of coefficients @var{numcoeff}, and the number of
+## octaves @var{numoctaves}. It works only for vectors currently.
+##
+## @end deftypefn
+
+function transform = nureal( t, x, omegamax, ncoeff, noctave)
+ ## the R function runs the following command:
+ ## nureal( double X, double Y, int min(X,Y), int ncoeff, int noctave, double omegamax, complex rp)
+ ## this means that in the C, *tptr is X and *xptr is Y. Yes, I know. I think I'll rename them.
+ ## n is the min of X and Y ... (as is k) and ncoeff is ... ncoeff, while noctave is noctave and
+ ## o is omegamax.
+ ## where rp = complex(noctave*ncoeff) so ... I can just store that as noctave*ncoeff and have no
+ ## problems, I guess.
+ ## Possibly throw an error if ncoeff <= 0.
+ k = n = min(length(x),length(t)); ## THIS IS VECTOR-ONLY. I'd need to add another bit of code to
+ ## make it array-safe, and that's not knowing right now what else will be necessary.
+ transform = zeros(1,(noctave * ncoeff)); ## In the C code, this is rendered as a Complex, but Octave doesn't care.
+ od = 2 ^ ( - 1 / ncoeff ); ## this will cause a crash if ncoeff=0; prefer error & quit?
+ o = omegamax;
+ ## ot is just defined as a Real here, I'm leaving it until it needs to be called.
+ n1 = 1 / n; ## written n_1 in the C, but I'd prefer to not get into underscores here.
+ ## zeta and iota are defined as Complex here, leaving them until they need to be defined.
+ ## I'm not convinced I won't want ncoeff again, so ...
+ ncoeffp = ncoeff;
+ ncoeffp *= noctave;
+ for iter = 1:ncoeffp
+ ## This method is an application of Eq. 8 on page 6 of the text, as well as Eq. 7
+ ot = o .* t;
+ zeta = sum( ( cos(ot) .* x ) - ( sin(ot) .* x .* i ) );
+ ot = ot .* 2;
+ iota = sum( cos(ot) - ( sin(ot) .* i ) );
+ zeta = zeta .* n1;
+ iota = iota .* n1;
+ transform(iter) = 2 / ( 1 - ( real(iota) ^ 2 ) - ( imag(iota) ^ 2 ) ) * ( conj(zeta) - (conj(iota) * zeta ));
+ o = o .* od;
+ endfor
+
+ ## transform = rp;
+
+endfunction
\ No newline at end of file
Copied: trunk/octave-forge/extra/lssa/lswaveletcoeff.m (from rev 10520, trunk/octave-forge/extra/lssa/nuwaveletcoeff.m)
===================================================================
--- trunk/octave-forge/extra/lssa/lswaveletcoeff.m (rev 0)
+++ trunk/octave-forge/extra/lssa/lswaveletcoeff.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -0,0 +1,55 @@
+## 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/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq)
+## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window)
+## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window, winradius)
+##
+## Return the coefficient of the wavelet transform of the
+## time series (@var{abc}, @var{ord}) at time @var{time}
+## and frequency @var{freq}; optional variable @var{window}
+## provides a windowing function and defaults to cubicwgt,
+## while @var{winradius} is the windowing radius, and defaults
+## to 1 (the radius of cubicwgt.)
+##
+## @end deftypefn
+
+%!demo
+%! x = 1:10;
+%! y = sin(x);
+%! xt = x';
+%! yt = y';
+%! a = nuwaveletcoeff(x,y,0.5,0.9)
+%! b = nuwaveletcoeff(xt,yt,0.5,0.9)
+%! ## Generates the wavelet transform coefficient for time 0.5 and circ. freq. 0.9, for row & column vectors.
+
+
+function coeff = nuwaveletcoeff( x , y , t , o , wgt = @cubicwgt , wgtrad = 1 )
+ so = 0.05 .* o;
+ if ( ( ndims(x) == 2 ) && ! ( rows(x) == 1 ) )
+ x = reshape(x,1,length(x));
+ y = reshape(y,1,length(y));
+ endif
+ 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
Deleted: trunk/octave-forge/extra/lssa/nucorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/nucorrcoeff.m 2012-05-24 20:41:47 UTC (rev 10520)
+++ trunk/octave-forge/extra/lssa/nucorrcoeff.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -1,67 +0,0 @@
-## 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/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq)
-## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window)
-## @deftypefnx {Function File} {c =} nucorrcoeff (abc1, ord1, abc2, ord2, time, freq, window, winradius)
-##
-## Return the coefficient of the wavelet correlation of time
-## series (@var{abc1}, @var{ord1}) and (@var{abc2}, @var{ord2}).
-## @var{window} is used to apply a windowing function, its
-## default is cubicwgt if left blank, and its radius is 1,
-## as defined as the default for @var{winradius}.
-##
-## @end deftypefn
-
-## Demo with sin, cos as Nir suggested.
-%!demo
-%! x = 1:10;
-%! y = sin(x);
-%! z = cos(x);
-%! a = nucorrcoeff(x,y,x,z,0.5,0.9)
-%! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9
-
-
-## 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;
- ## This code can only, as of currently, work on vectors; I haven't figured out a way to make it work on a matrix.
- if( ( ndims(x1) == 2 ) && ! ( rows(x1) == 1 ) )
- x1 = reshape(x1,1,length(x1));
- y1 = reshape(y1,1,length(y1));
- x2 = reshape(x2,1,length(x2));
- y2 = reshape(y2,1,length(y2));
- endif
- ## 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?
- lnength(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
Deleted: trunk/octave-forge/extra/lssa/nureal.m
===================================================================
--- trunk/octave-forge/extra/lssa/nureal.m 2012-05-24 20:41:47 UTC (rev 10520)
+++ trunk/octave-forge/extra/lssa/nureal.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -1,72 +0,0 @@
-## Copyright (C) 2012 Benjamin Lewis
-##
-## 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/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {transform =} nureal ( mag, time, maxfreq, numcoeff, numoctaves)
-##
-## Return the real least-squares transform of the time series
-## defined, based on the maximal frequency @var{maxfreq}, the
-## number of coefficients @var{numcoeff}, and the number of
-## octaves @var{numoctaves}. It works only for vectors currently.
-##
-## @end deftypefn
-
-function transform = nureal( x, t, omegamax, ncoeff, noctave)
- ## the R function runs the following command:
- ## nureal( double X, double Y, int min(X,Y), int ncoeff, int noctave, double omegamax, complex rp)
- ## this means that in the C, *tptr is X and *xptr is Y. Yes, I know. I think I'll rename them.
- ## n is the min of X and Y ... (as is k) and ncoeff is ... ncoeff, while noctave is noctave and
- ## o is omegamax.
- ## where rp = complex(noctave*ncoeff) so ... I can just store that as noctave*ncoeff and have no
- ## problems, I guess.
- ## Possibly throw an error if ncoeff <= 0.
- k = n = min ( min ( x , t ) ); ## THIS IS VECTOR-ONLY. I'd need to add another bit of code to
- ## make it array-safe, and that's not knowing right now what else will be necessary.
- rp = zeros(1,(noctave * ncoeff)); ## In the C code, this is rendered as a Complex, but Octave doesn't care.
- od = 2 ^ ( - 1 / ncoeff ); ## this will cause a crash if ncoeff=0; prefer error & quit?
- o = omegamax;
- ## ot is just defined as a Real here, I'm leaving it until it needs to be called.
- n1 = 1 / n; ## written n_1 in the C, but I'd prefer to not get into underscores here.
- ## zeta and iota are defined as Complex here, leaving them until they need to be defined.
- ## I'm not convinced I won't want ncoeff again, so ...
- ncoeffp = ncoeff;
- for ( ncoeffp *= noctave , iter = 1 ; sign(ncoeffp--) ; o *= od )
- #{
- zeta = iota = 0;
- for ( SRCFIRST ; SRCAVAIL ; SRCNEXT ) ##This is going to be vectorised shortly.
- ## This code can't run yet ... I'm going to work out what SRCFIRST, SRCAVAIL, SRCNEXT are and
- ## replace them with what they should be in this context. I've kept them as reminders.
- ot = o * SRCT; ## Same with SRCT ... I think it means what was originally set as Y. Macros?
- zeta += cos(ot) * SRCX;
- zeta -= sin(ot) * SRCX * i; ## More sure now. I don't think I can vectorise this ... who am I kidding?
- ot *= 2;
- iota += cos(ot);
- iota -= sin(ot) * i;
- endfor
- }#
- ## Commented out the converted-from-C code because it's replaced by the four lines below.
- ## This method is an application of Eq. 8 on page 6 of the text, as well as Eq. 7
- ot = o .* t;
- zeta = sum( ( cos(ot) .* x ) - ( sin(ot) .* x .* i ) );
- ot = ot .* 2;
- iota = sum( cos(ot) - ( sin(ot) .* i ) );
- zeta *= n1;
- iota *= n1;
- rp(iter++) = 2 / ( 1 - ( real(iota) ^ 2 ) - ( imag(iota) ^ 2 ) ) * ( conj(zeta) - (conj(iota) * zeta ));
- endfor
-
- transform = rp;
-
-endfunction
\ No newline at end of file
Deleted: trunk/octave-forge/extra/lssa/nuwaveletcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/nuwaveletcoeff.m 2012-05-24 20:41:47 UTC (rev 10520)
+++ trunk/octave-forge/extra/lssa/nuwaveletcoeff.m 2012-05-25 14:49:45 UTC (rev 10521)
@@ -1,55 +0,0 @@
-## 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/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq)
-## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window)
-## @deftypefnx {Function File} {c =} nuwaveletcoeff (abc, ord, time, freq, window, winradius)
-##
-## Return the coefficient of the wavelet transform of the
-## time series (@var{abc}, @var{ord}) at time @var{time}
-## and frequency @var{freq}; optional variable @var{window}
-## provides a windowing function and defaults to cubicwgt,
-## while @var{winradius} is the windowing radius, and defaults
-## to 1 (the radius of cubicwgt.)
-##
-## @end deftypefn
-
-%!demo
-%! x = 1:10;
-%! y = sin(x);
-%! xt = x';
-%! yt = y';
-%! a = nuwaveletcoeff(x,y,0.5,0.9)
-%! b = nuwaveletcoeff(xt,yt,0.5,0.9)
-%! ## Generates the wavelet transform coefficient for time 0.5 and circ. freq. 0.9, for row & column vectors.
-
-
-function coeff = nuwaveletcoeff( x , y , t , o , wgt = @cubicwgt , wgtrad = 1 )
- so = 0.05 .* o;
- if ( ( ndims(x) == 2 ) && ! ( rows(x) == 1 ) )
- x = reshape(x,1,length(x));
- y = reshape(y,1,length(y));
- endif
- 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.
|
