Menu
▾
▴
octave-cvsupdate
|
From: <be...@us...> - 2012-08-03 15:19:47
|
Revision: 10804
http://octave.svn.sourceforge.net/octave/?rev=10804&view=rev
Author: benjf5
Date: 2012-08-03 15:19:41 +0000 (Fri, 03 Aug 2012)
Log Message:
-----------
Tests, having been rewritten to match the manual, still broken.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/inst/lombcoeff.m
trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
trunk/octave-forge/extra/lssa/inst/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
trunk/octave-forge/extra/lssa/inst/lsreal.m
trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m
Modified: trunk/octave-forge/extra/lssa/inst/lombcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -23,14 +23,14 @@
## @end deftypefn
%!test
-%! shared t, x, o, maxfreq
+%!shared t, x, o, maxfreq;
%! maxfreq = 4 / ( 2 * pi );
-%! t = linspace(0,8); x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t)
-%! - 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t)
-%! + cos((1/4)*maxfreq.*t)); o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
-%!assert( lombcoeff(t,x,o(1)),10788.9848389923,5e-10 );
-%!assert( lombcoeff(t,x,o(2)),12352.6413413457,5e-10 );
-%!assert( lombcoeff(t,x,o(3)),13673.4098969780,5e-10 );
+%! t = linspace(0,8);
+%! x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t) - 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t) + cos((1/4)*maxfreq.*t));
+%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
+%!assert( lombcoeff(t,x,maxfreq),10788.9848389923,5e-10 );
+%!assert( lombcoeff(t,x,3/4*maxfreq),12352.6413413457,5e-10 );
+%!assert( lombcoeff(t,x,1/4*maxfreq),13673.4098969780,5e-10 );
function coeff = lombcoeff(T, X, o)
Modified: trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -22,6 +22,17 @@
##
## @end deftypefn
+%!test
+%!shared t, x, o, maxfreq;
+%! maxfreq = 4 / ( 2 * pi );
+%! t = linspace(0,8);
+%! x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t) - 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t) + cos((1/4)*maxfreq.*t));
+%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
+%!assert( lombnormcoeff(t,x,o(1)),63.3294946603949,5e-10 );
+%!assert( lombnormcoeff(t,x,o(2)),73.2601360674868,5e-10 );
+%!assert( lombnormcoeff(t,x,o(3)),53.0799752083903,5e-10 );
+
+
function coeff = lombnormcoeff(T,X,omega)
tau = atan2( sum( sin( 2.*omega.*T)), sum(cos(2.*omega.*T))) / 2;
coeff = ( ( sum ( X .* cos( omega .* T - tau ) ) .^ 2 ./ sum ( cos ( omega .* T - tau ) .^ 2 )
Modified: trunk/octave-forge/extra/lssa/inst/lscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -24,9 +24,12 @@
## @end deftypefn
%!test
-%! shared t, x, o, maxfreq
-%! maxfreq = 4 / ( 2 * pi ); t = [0:0.008:8]; x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t)- 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t) + cos((1/4)*maxfreq.*t)); o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
-%! assert( lscomplex(t,x,maxfreq,2,2), [-0.400754376933531 - 2.366871097665244i, 1.226663545950135 - 2.243899314661490i, 1.936433327880238 - 1.515538553198501i, 2.125045509991203 - 0.954100898917708i ], 6e-14 );
+%!shared t, x, o, maxfreq
+%! maxfreq = 4 / ( 2 * pi );
+%! t = [0:0.008:8];
+%! x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t)- 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t) + cos((1/4)*maxfreq.*t));
+%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
+%!assert( lscomplex(t,x,maxfreq,2,2), [-0.400754376933531 - 2.366871097665244i, 1.226663545950135 - 2.243899314661490i, 1.936433327880238 - 1.515538553198501i, 2.125045509991203 - 0.954100898917708i ], 6e-14 );
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -50,7 +50,6 @@
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 = ;
mask = find( ( abs( x1 - t ) * so ) < wgtrad );
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!
@@ -62,8 +61,7 @@
## in particular to maintain an exact duplicate of the R function.
s = sum( wgt( ( rx1 - t ) .* so ) ) * sum( wgt( ( rx2 - t ) .* so ) );
if s != 0
- coeff = sum( wgt((rx1-t).*so).*exp(i*o.*rx1).*ry1) *
- sum(wgt((rx2-t).*so).*exp(i*o.*rx2).*conj(ry2)) / s;
+ coeff = sum( wgt((rx1-t).*so).*exp(i*o.*rx1).*ry1) * sum(wgt((rx2-t).*so).*exp(i*o.*rx2).*conj(ry2)) / s;
else
coeff = 0;
endif
Modified: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -29,7 +29,7 @@
%! maxfreq = 4 / ( 2 * pi );
%! t = linspace(0,8);
%! x = ( 2.*sin(maxfreq.*t) + 3.*sin((3/4)*maxfreq.*t) - 0.5 .* sin((1/4)*maxfreq.*t) - 0.2 .* cos(maxfreq .* t) + cos((1/4)*maxfreq.*t));
-%! assert(lsreal(t,x,maxfreq,2,2),[-1.68275915310663 + 4.70126183846743i, 1.93821553170889 + 4.95660209883437i, 4.38145452686697 + 2.14403733658600i, 5.27425332281147 - 0.73933440226597i],6e-14)
+%!assert(lsreal(t,x,maxfreq,2,2),[-1.68275915310663 + 4.70126183846743i, 1.93821553170889 + 4.95660209883437i, 4.38145452686697 + 2.14403733658600i, 5.27425332281147 - 0.73933440226597i],6e-14)
%! #In the assert here, I've got an error bound large enough to catch individual system errors which would present no real issue.
function transform = lsreal( t, x, omegamax, ncoeff, noctave)
Modified: trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m 2012-08-03 12:22:15 UTC (rev 10803)
+++ trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m 2012-08-03 15:19:41 UTC (rev 10804)
@@ -14,7 +14,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
##-*- texinfo -*-
-## @deffun {Function File} {t =} lsrealwavelet( @var{time}, @var{mag},
+## @deftypefn {Function File} {t =} lsrealwavelet( @var{time}, @var{mag},
## @var{maxfreq}, @var{coefficients}, @var{octaves}, @var{time_min},
## @var{time_max}, @var{min_window_count} )
##
@@ -26,7 +26,7 @@
##
## @seealso lscomplexwavelet lswaveletcoeff lscorrcoeff
##
-## @end deffun
+## @end deftypefn
function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number )
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-08-09 12:59:01
|
Revision: 10848
http://octave.svn.sourceforge.net/octave/?rev=10848&view=rev
Author: benjf5
Date: 2012-08-09 12:58:55 +0000 (Thu, 09 Aug 2012)
Log Message:
-----------
A few fixes and some tests for lscorrcoeff.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
Modified: trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-09 02:11:43 UTC (rev 10847)
+++ trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-09 12:58:55 UTC (rev 10848)
@@ -15,39 +15,47 @@
-function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, minimum_window_count, sigma = 0.05)
+function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, tmin, tmax, tstep, sigma = 0.05)
-## 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.
+ ## 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.
+
+ transform = cell(noctave*ncoeff,1);
+
+ for octave_iter = 1:noctave
+ ## In fastnu.c, winrad is set as π/(sigma*omegaoct); I suppose this is
+ ## ... feasible, although it will need to be noted that if sigma is set too
+ ## large, the windows will exclude data. I can work with that.
+ ##
+ ## An additional consideration is that
+
+ for coeff_iter = 1:ncoeff
+
+ ## in this, win_t is the centre of the window in question
+ ## Although that will vary depending on the window. This is just an
+ ## implementation for the first window.
+
+ current_iteration = (octave_iter-1)*ncoeff+coeff_iter;
+ window_radius = pi / ( sigma * omegamax * ( 2 ^ ( current_iteration - 1 ) ) );
+ window_count = 2 * ceil ( ( tmax - tmin ) / window_step ) - 1;
+ omega = current_frequency = maxfreq * 2 ^ ( - octave_iter*coeff_iter / ncoeff );
+
+
+
+ transform{current_iteration}=zeros(1,window_count);
-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.
- omega = 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
- ## Computes the transform as stated in the paper for each given frequency.
- zeta = sum ( cubicwgt ( sigma .* omega .* ( T - win_t ) ) .* exp ( -i .* omega .* ( T - win_t ) ) .* X ) / sum ( cubicwgt ( sigma .* omega .* ( T - win_t ) ) .* exp ( -i .* omega .* ( T - win_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?
+ ## win_t is the centre of the current window.
+ win_t = tmin + window_radius;
+ for iter_window = 1:window_count
+ ## Computes the transform as stated in the paper for each given frequency.
+ zeta = sum ( cubicwgt ( sigma .* omega .* ( T - win_t ) ) .* exp ( -i .* omega .* ( T - win_t ) ) .* X ) / sum ( cubicwgt ( sigma .* omega .* ( T - win_t ) ) .* exp ( -i .* omega .* ( T - win_t ) ) );
+ transform{current_iteration}(iter_window) = zeta;
+ window_min += window_radius ;
endfor
endfor
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-09 02:11:43 UTC (rev 10847)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-09 12:58:55 UTC (rev 10848)
@@ -59,11 +59,38 @@
## 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 included it
## in particular to maintain an exact duplicate of the R function.
- s = sum( wgt( ( rx1 - t ) .* so ) ) * sum( wgt( ( rx2 - t ) .* so ) );
+ s = sum ( wgt ( ( rx1 - t ) .* so ) ) * sum ( wgt ( ( rx2 - t ) .* so ) );
if s != 0
- coeff = sum( wgt((rx1-t).*so).*exp(i*o.*rx1).*ry1) * sum(wgt((rx2-t).*so).*exp(i*o.*rx2).*conj(ry2)) / s;
+ coeff = sum ( wgt ((rx1-t).*so).*exp(i*o.*rx1).*ry1) * sum(wgt((rx2-t).*so).*exp(i*o.*rx2).*conj(ry2)) / s;
else
coeff = 0;
endif
endfunction
+
+
+%!shared t, p, x, y, z, o, maxfreq
+%! maxfreq = 4 / (2 * pi);
+%! t = linspace (0, 8);
+%! x = (2 .* sin (maxfreq .* t) +
+%! 3 .* sin ((3/4) * maxfreq .* t) -
+%! 0.5 .* sin ((1/4) * maxfreq .* t) -
+%! 0.2 .* cos (maxfreq .* t) +
+%! cos ((1/4) * maxfreq .* t));
+%! y = - x;
+%! p = linspace (0, 8, 500);
+%! z = (2 .* sin (maxfreq .* p) +
+%! 3 .* sin ((3/4) * maxfreq .* p) -
+%! 0.5 .* sin ((1/4) * maxfreq .* p) -
+%! 0.2 .* cos (maxfreq .* p) +
+%! cos ((1/4) * maxfreq .* p));
+%! o = [maxfreq , (3/4 * maxfreq) , (1/4 * maxfreq)];
+%!assert (lscorrcoeff (t, x, t, x, 0.5, maxfreq), -5.54390340863576 -
+%! 1.82439880893383i, 5e-10);
+%!assert (lscorrcoeff (t, x, t, y, 0.5, maxfreq), 5.54390340863576 +
+%! 1.82439880893383i, 5e-10);
+%!assert (lscorrcoeff (t, x, p, z, 0.5, maxfreq), -5.55636741054624 -
+%! 1.82803733863170i, 5e-10);
+
+
+
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-08-09 16:45:39
|
Revision: 10850
http://octave.svn.sourceforge.net/octave/?rev=10850&view=rev
Author: benjf5
Date: 2012-08-09 16:45:30 +0000 (Thu, 09 Aug 2012)
Log Message:
-----------
Really the last commit for lssa release 0.1.1.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/inst/cubicwgt.m
trunk/octave-forge/extra/lssa/inst/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m
trunk/octave-forge/extra/lssa/inst/lsreal.m
Modified: trunk/octave-forge/extra/lssa/inst/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-09 16:08:18 UTC (rev 10849)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-09 16:45:30 UTC (rev 10850)
@@ -34,7 +34,6 @@
%! k = [ 0 , 3 , 1.5, -1, -0.5, -0.25, 0.75 ];
%!assert( cubicwgt(h), 0 );
%!assert( cubicwgt(m), 1 + m ^ 2 * ( 2 * m - 3 ));
-%!assert( cubicwgt(k), [ 1.00000, 0.00000, 0.00000, 0.00000, 0.50000,
-%! 0.84375, 0.15625], 1e-6);
+%!assert( cubicwgt(k), [ 1.00000, 0.00000, 0.00000, 0.00000, 0.50000, 0.84375, 0.15625], 1e-6);
%! ## Tests cubicwgt on two scalars and two vectors; cubicwgt will work
%! ## on any array input.
\ No newline at end of file
Modified: trunk/octave-forge/extra/lssa/inst/lscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-09 16:08:18 UTC (rev 10849)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-09 16:45:30 UTC (rev 10850)
@@ -48,7 +48,4 @@
%! cos ((1/4)*maxfreq.*t));
%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
%!assert( lscomplex(t,x,maxfreq,2,2),
-%! [-0.400924546169395 - 2.371555305867469i,
-%! 1.218065147708429 - 2.256125004156890i,
-%! 1.935428592212907 - 1.539488163739336i,
-%! 2.136692292751917 - 0.980532175174563i ], 5e-10 );
\ No newline at end of file
+%! [-0.400924546169395 - 2.371555305867469i, 1.218065147708429 - 2.256125004156890i, 1.935428592212907 - 1.539488163739336i, 2.136692292751917 - 0.980532175174563i ], 5e-10 );
\ No newline at end of file
Modified: trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-09 16:08:18 UTC (rev 10849)
+++ trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-09 16:45:30 UTC (rev 10850)
@@ -15,7 +15,7 @@
-function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, tmin, tmax, tstep, sigma = 0.05)
+function transform = lscomplexwavelet( T, X, omegamax, ncoeff, noctave, tmin, tmax, tstep, sigma = 0.05)
## 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
@@ -42,8 +42,8 @@
current_iteration = (octave_iter-1)*ncoeff+coeff_iter;
window_radius = pi / ( sigma * omegamax * ( 2 ^ ( current_iteration - 1 ) ) );
- window_count = 2 * ceil ( ( tmax - tmin ) / window_step ) - 1;
- omega = current_frequency = maxfreq * 2 ^ ( - octave_iter*coeff_iter / ncoeff );
+ window_count = 2 * ceil ( ( tmax - tmin ) / window_radius ) - 1;
+ omega = current_frequency = omegamax * 2 ^ ( - octave_iter*coeff_iter / ncoeff );
Modified: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-09 16:08:18 UTC (rev 10849)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-09 16:45:30 UTC (rev 10850)
@@ -70,10 +70,7 @@
%! 0.2 .* cos ( maxfreq .* t ) +
%! cos ( (1/4) * maxfreq .* t ) );
%!assert (lsreal (t,x,maxfreq,2,2),
-%! [-1.68275915310663 + 4.70126183846743i,
-%! 1.93821553170889 + 4.95660209883437i,
-%! 4.38145452686697 + 2.14403733658600i,
-%! 5.27425332281147 - 0.73933440226597i ],
+%! [-1.68275915310663 + 4.70126183846743i, 1.93821553170889 + 4.95660209883437i, 4.38145452686697 + 2.14403733658600i, 5.27425332281147 - 0.73933440226597i ],
%! 5e-10)
%! # In the assert here, I've got an error bound large enough to catch
%! # individual system errors which would present no real issue.
\ No newline at end of file
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <cd...@us...> - 2012-08-10 09:30:20
|
Revision: 10852
http://octave.svn.sourceforge.net/octave/?rev=10852&view=rev
Author: cdf
Date: 2012-08-10 09:30:12 +0000 (Fri, 10 Aug 2012)
Log Message:
-----------
code review part 1
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/inst/cubicwgt.m
trunk/octave-forge/extra/lssa/inst/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
trunk/octave-forge/extra/lssa/inst/lsreal.m
Modified: trunk/octave-forge/extra/lssa/inst/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-09 23:49:21 UTC (rev 10851)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-10 09:30:12 UTC (rev 10852)
@@ -13,27 +13,29 @@
## 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} {@var{a} =} cubicwgt (@var{series})
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{a} =} cubicwgt (@var{series})
## Return @var{series} as windowed by a cubic polynomial,
## 1 + ( x ^ 2 * ( 2 x - 3 ) ), assuming x is in [-1,1].
-## @end deftypefn
-
## This function implements the windowing function on page 10 of the paper.
## 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(s) ## where s is the value to be windowed
- a = abs(s);
- a = ifelse( ( a < 1 ), 1 + ( ( a .^ 2 ) .* ( 2 .* a - 3 ) ), a = 0);
+## @end deftypefn
+
+function a = cubicwgt (s)
+ ## s is the value to be windowed
+ a = abs (s);
+ a = ifelse ((a < 1), 1 + ((a .^ 2) .* (2 .* a - 3)), 0);
endfunction
-%!test
+
%!shared h, m, k
%! h = 2;
%! m = 0.01;
-%! k = [ 0 , 3 , 1.5, -1, -0.5, -0.25, 0.75 ];
-%!assert( cubicwgt(h), 0 );
-%!assert( cubicwgt(m), 1 + m ^ 2 * ( 2 * m - 3 ));
-%!assert( cubicwgt(k), [ 1.00000, 0.00000, 0.00000, 0.00000, 0.50000, 0.84375, 0.15625], 1e-6);
+%! k = [0, 3, 1.5, -1, -0.5, -0.25, 0.75];
+%!assert (cubicwgt (h), 0 );
+%!assert (cubicwgt (m), 1 + m ^ 2 * (2 * m - 3));
+%!assert (cubicwgt (k), [1.00000, 0.00000, 0.00000, 0.00000, ...
+%! 0.50000, 0.84375, 0.15625], 1e-6);
%! ## Tests cubicwgt on two scalars and two vectors; cubicwgt will work
%! ## on any array input.
Modified: trunk/octave-forge/extra/lssa/inst/lscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-09 23:49:21 UTC (rev 10851)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-10 09:30:12 UTC (rev 10852)
@@ -24,39 +24,45 @@
## @end deftypefn
-function transform = lscomplex( t , x , omegamax , ncoeff , noctave )
+function transform = lscomplex (t, x, omegamax, ncoeff, noctave)
- n = length(t); ## VECTOR ONLY, and since t and x have the same number of
+ ## VECTOR ONLY, and since t and x have the same number of
## entries, there's no problem.
+ n = length (t);
+
- transform = zeros(1,ncoeff*noctave);
+ transform = zeros (1, ncoeff * noctave);
o = omegamax;
omul = 2 ^ (- 1 / ncoeff);
- for iter = 1:ncoeff*noctave
+ for iter = 1:ncoeff * noctave
ot = o .* t;
- transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n; ## See the
- ## paper for the expression
+ ## See the paper for the expression below
+ transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n;
+
- o *= omul; ## To advance the transform to the next coefficient in the octave
+ ## Advance the transform to the next coefficient in the octave
+ o *= omul;
endfor
endfunction
%!test
-%!shared t, x, o, maxfreq
%! maxfreq = 4 / ( 2 * pi );
%! t = [0:0.008:8];
%! x = ( 2 .* sin (maxfreq .* t) +
%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
%! 0.5 .* sin ((1/4) * maxfreq .* t) -
%! 0.2 .* cos (maxfreq .* t) +
-%! cos ((1/4)*maxfreq.*t));
+%! cos ((1/4) * maxfreq .* t));
%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
-%!assert( lscomplex(t,x,maxfreq,2,2),
-%! [-0.400924546169395 - 2.371555305867469i, 1.218065147708429 - 2.256125004156890i, 1.935428592212907 - 1.539488163739336i, 2.136692292751917 - 0.980532175174563i ], 5e-10 );
+%! assert (lscomplex (t, x, maxfreq, 2, 2),
+%! [(-0.400924546169395 - 2.371555305867469i), ...
+%! (1.218065147708429 - 2.256125004156890i), ...
+%! (1.935428592212907 - 1.539488163739336i), ...
+%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-09 23:49:21 UTC (rev 10851)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-10 09:30:12 UTC (rev 10852)
@@ -31,47 +31,51 @@
##
## @end deftypefn
-## Demo with sin, cos as Nir suggested.
-%!demo
-%! x = 1:10;
-%! y = sin(x);
-%! z = cos(x);
-%! a = lscorrcoeff(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 = lscorrcoeff (x1, y1, x2, y2, t, o, wgt = @cubicwgt, wgtrad = 1)
-## nucorrcoeff, computes a coefficient of the wavelet correlation of two time series
+ so = 0.05 * o;
-function coeff = lscorrcoeff(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));
+ ## 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 = find( ( abs( x1 - t ) * so ) < wgtrad );
- 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 = find( ( abs( x2 - t ) * so ) < wgtrad );
+
+ ## 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;
+ ## I've kept the variable names from the R function here
+ rx1 = x1(mask);
+ ## FIXME : Needs to have a noisy error if length(y1) != length(x1) -- add this!
+ ry1 = y1(mask);
+
+ ## I've used the same mask for all of these as it's an
+ ## otherwise unimportant variable ... can this leak memory?
+ mask = (abs (x2 - t) * so ) < wgtrad;
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 included it
+
+ ## printing this length is probably used as a warning if 0 is returned;
+ ## I included it
+ length (rx1)
+
## in particular to maintain an exact duplicate of the R function.
- s = sum ( wgt ( ( rx1 - t ) .* so ) ) * sum ( wgt ( ( rx2 - t ) .* so ) );
- if s != 0
- coeff = sum ( wgt ((rx1-t).*so).*exp(i*o.*rx1).*ry1) * sum(wgt((rx2-t).*so).*exp(i*o.*rx2).*conj(ry2)) / s;
+ s = sum (wgt ((rx1 - t) .* so)) * sum (wgt ((rx2 - t ) .* so ));
+ if (s != 0)
+ coeff = sum (wgt ((rx1 - t) .* so) .* exp (i * o .* rx1) .* ry1) * ...
+ sum (wgt ((rx2 - t) .* so) .* exp (i * o .* rx2) .* conj (ry2)) / s;
else
coeff = 0;
endif
endfunction
-%!test
+
%!shared t, p, x, y, z, o, maxfreq
%! maxfreq = 4 / (2 * pi);
%! t = linspace (0, 8);
@@ -88,12 +92,18 @@
%! 0.2 .* cos (maxfreq .* p) +
%! cos ((1/4) * maxfreq .* p));
%! o = [maxfreq , (3/4 * maxfreq) , (1/4 * maxfreq)];
-%!assert (lscorrcoeff (t, x, t, x, 0.5, maxfreq), -5.54390340863576 -
-%! 1.82439880893383i, 5e-10);
-%!assert (lscorrcoeff (t, x, t, y, 0.5, maxfreq), 5.54390340863576 +
-%! 1.82439880893383i, 5e-10);
-%!assert (lscorrcoeff (t, x, p, z, 0.5, maxfreq), -5.55636741054624 -
-%! 1.82803733863170i, 5e-10);
+%!assert (lscorrcoeff (t, x, t, x, 0.5, maxfreq),
+%! -5.54390340863576 - 1.82439880893383i, 5e-10);
+%!assert (lscorrcoeff (t, x, t, y, 0.5, maxfreq),
+%! 5.54390340863576 + 1.82439880893383i, 5e-10);
+%!assert (lscorrcoeff (t, x, p, z, 0.5, maxfreq),
+%! -5.55636741054624 - 1.82803733863170i, 5e-10);
+## Demo with sin, cos as Nir suggested.
+%!demo
+%! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9
+%! x = 1:10;
+%! y = sin (x);
+%! z = cos (x);
+%! a = lscorrcoeff (x, y, x, z, 0.5, 0.9)
-
Modified: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-09 23:49:21 UTC (rev 10851)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-10 09:30:12 UTC (rev 10852)
@@ -28,44 +28,40 @@
function transform = lsreal (t, x, omegamax, ncoeff, noctave)
- k = n = 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));
+ ## FIXME : 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.
+ k = n = length (t);
- od = 2 ^ (- 1 / ncoeff);
+ transform = zeros (1, (noctave * ncoeff));
+ od = 2 ^ (- 1 / ncoeff);
o = omegamax;
-
n1 = 1 / n;
- ncoeffp = ncoeff;
+ ncoeffp = ncoeff * noctave;
- 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
+ ## 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));
+ zeta = n1 * sum ((cos (ot) - i * sin (ot)) .* x);
- ot = ot .* 2;
+ ot *= 2;
- iota = sum (cos (ot) -
- (sin (ot) .* i));
+ iota = n1 * sum (cos (ot) - i * sin (ot));
- zeta = zeta .* n1;
- iota = iota .* n1;
+ transform(iter) = (2 * (conj (zeta) - (conj (iota) * zeta)) /
+ (1 - (real (iota) ^ 2) - (imag (iota) ^ 2)));
- transform(iter) = (2 / (1 -(real (iota) ^ 2) - (imag(iota) ^ 2)) * (conj(zeta) - (conj(iota) * zeta)));
- o = o .* od;
+ o *= od;
endfor
endfunction
%!test
-%!shared t, x, o, maxfreq
%! maxfreq = 4 / ( 2 * pi );
%! t = linspace(0,8);
%! x = ( 2 .* sin ( maxfreq .* t ) +
@@ -73,8 +69,11 @@
%! 0.5 .* sin ( (1/4) * maxfreq .* t ) -
%! 0.2 .* cos ( maxfreq .* t ) +
%! cos ( (1/4) * maxfreq .* t ) );
-%!assert (lsreal (t,x,maxfreq,2,2),
-%! [-1.68275915310663 + 4.70126183846743i, 1.93821553170889 + 4.95660209883437i, 4.38145452686697 + 2.14403733658600i, 5.27425332281147 - 0.73933440226597i ],
-%! 5e-10)
%! # In the assert here, I've got an error bound large enough to catch
%! # individual system errors which would present no real issue.
+%! assert (lsreal (t,x,maxfreq,2,2),
+%! [(-1.68275915310663 + 4.70126183846743i), ...
+%! (1.93821553170889 + 4.95660209883437i), ...
+%! (4.38145452686697 + 2.14403733658600i), ...
+%! (5.27425332281147 - 0.73933440226597i)],
+%! 5e-10)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-08-14 13:42:59
|
Revision: 10860
http://octave.svn.sourceforge.net/octave/?rev=10860&view=rev
Author: benjf5
Date: 2012-08-14 13:42:51 +0000 (Tue, 14 Aug 2012)
Log Message:
-----------
Moved to loopless LS transforms, rewrote and improved help strings for all functions in lssa.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/inst/cubicwgt.m
trunk/octave-forge/extra/lssa/inst/lombcoeff.m
trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
Added Paths:
-----------
trunk/octave-forge/extra/lssa/inst/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lscomplex.m.old
trunk/octave-forge/extra/lssa/inst/lsreal.m
trunk/octave-forge/extra/lssa/inst/lsreal.m.old
Removed Paths:
-------------
trunk/octave-forge/extra/lssa/inst/fastlscomplex.m
trunk/octave-forge/extra/lssa/inst/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lsreal.m
Modified: trunk/octave-forge/extra/lssa/inst/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -15,11 +15,19 @@
## -*- texinfo -*-
## @deftypefn {Function File} {@var{a} =} cubicwgt (@var{series})
-## Return @var{series} as windowed by a cubic polynomial,
-## 1 + ( x ^ 2 * ( 2 x - 3 ) ), assuming x is in [-1,1].
+##
+## Returns the input series, windowed by a polynomial similar to a Hanning
+## window. To window an arbitrary section of the series, subtract or add an
+## offset to it to adjust the centre of the window; for an offset of k, the call
+## would be cubicwgt (@var{s} - k). Similarly, the radius of the window is 1;
+## if an arbitrary radius r is desired, dividing the series by the radius after
+## centering is the best way to adjust to fit the window: cubicwgt ((@var{s} -
+## k) / r).
+##
+## The windowing function itself is:
+## w = 1 + ( x ^ 2 * ( 2 x - 3 ) ), x in [-1,1], else w = 0.
## This function implements the windowing function on page 10 of the paper.
-## if t is in [-1,1] then the windowed term is a = 1 + ( |t|^2 * ( 2|t| - 3 )
-## else the windowed term is 0.
+##
## @end deftypefn
function a = cubicwgt (s)
Deleted: trunk/octave-forge/extra/lssa/inst/fastlscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/fastlscomplex.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/fastlscomplex.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -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 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} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
-##
-## Return the complex least-squares transform of the (@var{time},@var{mag})
-## series, considering frequencies up to @var{maxfreq}, over @var{numoctaves}
-## octaves and @var{numcoeff} coefficients.
-##
-## @seealso{lsreal}
-## @end deftypefn
-
-
-function transform = fastlscomplex (t, x, omegamax, ncoeff, noctave)
-
- ## t will be unrolled to a column vector below
- ## no metter what its original shape is
- n = numel (t);
-
- iter = 0 : (ncoeff * noctave - 1);
- omul = (2 .^ (- iter / ncoeff));
-
- ot = t(:) * (omul * omegamax);
-
- ## See the paper for the expression below
- transform = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
-
-endfunction
-
-%!test
-%! maxfreq = 4 / ( 2 * pi );
-%! t = [0:0.008:8];
-%! x = ( 2 .* sin (maxfreq .* t) +
-%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
-%! 0.5 .* sin ((1/4) * maxfreq .* t) -
-%! 0.2 .* cos (maxfreq .* t) +
-%! cos ((1/4) * maxfreq .* t));
-%! assert (fastlscomplex (t, x, maxfreq, 2, 2),
-%! [(-0.400924546169395 - 2.371555305867469i), ...
-%! (1.218065147708429 - 2.256125004156890i), ...
-%! (1.935428592212907 - 1.539488163739336i), ...
-%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
Modified: trunk/octave-forge/extra/lssa/inst/lombcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -16,8 +16,7 @@
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lombcoeff (@var{time}, @var{mag}, @var{freq})
##
-## Return the coefficient of the Lomb periodogram (unnormalized) for the
-## (@var{time},@var{mag}) series for the @var{freq} provided.
+## Return the Lomb Periodogram value at one frequency for a time series.
##
## @seealso{lombnormcoeff}
## @end deftypefn
Modified: trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -16,10 +16,11 @@
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lombnormcoeff (@var{time}, @var{mag}, @var{freq})
##
-## Return the coefficient of the Lomb Normalised Periodogram at the
-## specified @var{frequency} of the periodogram applied to the
-## (@var{time}, @var{mag}) series.
+## Return the normalized Lomb Periodogram value at one frequency for a time
+## series.
##
+## @seealso{lombcoeff}
+##
## @end deftypefn
Deleted: trunk/octave-forge/extra/lssa/inst/lscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -1,68 +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} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
-##
-## Return the complex least-squares transform of the (@var{time},@var{mag})
-## series, considering frequencies up to @var{maxfreq}, over @var{numoctaves}
-## octaves and @var{numcoeff} coefficients.
-##
-## @seealso{lsreal}
-## @end deftypefn
-
-
-function transform = lscomplex (t, x, omegamax, ncoeff, noctave)
-
- ## VECTOR ONLY, and since t and x have the same number of
- ## entries, there's no problem.
- n = length (t);
-
-
- transform = zeros (1, ncoeff * noctave);
-
- o = omegamax;
-
- omul = 2 ^ (- 1 / ncoeff);
-
- for iter = 1:ncoeff * noctave
-
- ot = o .* t;
-
- ## See the paper for the expression below
- transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n;
-
-
- ## Advance the transform to the next coefficient in the octave
- o *= omul;
-
- endfor
-
-endfunction
-
-%!test
-%! maxfreq = 4 / ( 2 * pi );
-%! t = [0:0.008:8];
-%! x = ( 2 .* sin (maxfreq .* t) +
-%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
-%! 0.5 .* sin ((1/4) * maxfreq .* t) -
-%! 0.2 .* cos (maxfreq .* t) +
-%! cos ((1/4) * maxfreq .* t));
-%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
-%! assert (lscomplex (t, x, maxfreq, 2, 2),
-%! [(-0.400924546169395 - 2.371555305867469i), ...
-%! (1.218065147708429 - 2.256125004156890i), ...
-%! (1.935428592212907 - 1.539488163739336i), ...
-%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
Copied: trunk/octave-forge/extra/lssa/inst/lscomplex.m (from rev 10856, trunk/octave-forge/extra/lssa/inst/fastlscomplex.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -0,0 +1,61 @@
+## 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} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
+##
+## Return a series of least-squares transforms of a complex-valued time series.
+## Each transform is minimized independently at each frequency. @var{numcoeff}
+## frequencies are tested for each of @var{numoctaves} octaves, starting from
+## @var{maxfreq}.
+##
+## Each result (a + bi) at a given frequency, o, defines the real and imaginary
+## coefficients for a sum of cosine and sine functions: a cos(ot) + b i
+## sin(ot). The specific frequency can be determined by its index in @var{t},
+## @var{ind}, as @var{maxfreq} * 2 ^ (- (@var{ind} - 1) / @var{numcoeff}).
+##
+## @seealso{lsreal}
+## @end deftypefn
+
+
+function transform = lscomplex (t, x, omegamax, ncoeff, noctave)
+
+ ## t will be unrolled to a column vector below
+ ## no metter what its original shape is
+ n = numel (t);
+
+ iter = 0 : (ncoeff * noctave - 1);
+ omul = (2 .^ (- iter / ncoeff));
+
+ ot = t(:) * (omul * omegamax);
+
+ ## See the paper for the expression below
+ transform = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
+
+endfunction
+
+%!test
+%! maxfreq = 4 / ( 2 * pi );
+%! t = [0:0.008:8];
+%! x = ( 2 .* sin (maxfreq .* t) +
+%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
+%! 0.5 .* sin ((1/4) * maxfreq .* t) -
+%! 0.2 .* cos (maxfreq .* t) +
+%! cos ((1/4) * maxfreq .* t));
+%! assert (fastlscomplex (t, x, maxfreq, 2, 2),
+%! [(-0.400924546169395 - 2.371555305867469i), ...
+%! (1.218065147708429 - 2.256125004156890i), ...
+%! (1.935428592212907 - 1.539488163739336i), ...
+%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
Copied: trunk/octave-forge/extra/lssa/inst/lscomplex.m.old (from rev 10856, trunk/octave-forge/extra/lssa/inst/lscomplex.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m.old (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m.old 2012-08-14 13:42:51 UTC (rev 10860)
@@ -0,0 +1,68 @@
+## 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} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
+##
+## Return the complex least-squares transform of the (@var{time},@var{mag})
+## series, considering frequencies up to @var{maxfreq}, over @var{numoctaves}
+## octaves and @var{numcoeff} coefficients.
+##
+## @seealso{lsreal}
+## @end deftypefn
+
+
+function transform = lscomplex (t, x, omegamax, ncoeff, noctave)
+
+ ## VECTOR ONLY, and since t and x have the same number of
+ ## entries, there's no problem.
+ n = length (t);
+
+
+ transform = zeros (1, ncoeff * noctave);
+
+ o = omegamax;
+
+ omul = 2 ^ (- 1 / ncoeff);
+
+ for iter = 1:ncoeff * noctave
+
+ ot = o .* t;
+
+ ## See the paper for the expression below
+ transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n;
+
+
+ ## Advance the transform to the next coefficient in the octave
+ o *= omul;
+
+ endfor
+
+endfunction
+
+%!test
+%! maxfreq = 4 / ( 2 * pi );
+%! t = [0:0.008:8];
+%! x = ( 2 .* sin (maxfreq .* t) +
+%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
+%! 0.5 .* sin ((1/4) * maxfreq .* t) -
+%! 0.2 .* cos (maxfreq .* t) +
+%! cos ((1/4) * maxfreq .* t));
+%! o = [ maxfreq , 3 / 4 * maxfreq , 1 / 4 * maxfreq ];
+%! assert (lscomplex (t, x, maxfreq, 2, 2),
+%! [(-0.400924546169395 - 2.371555305867469i), ...
+%! (1.218065147708429 - 2.256125004156890i), ...
+%! (1.935428592212907 - 1.539488163739336i), ...
+%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -15,24 +15,22 @@
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq})
-## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1},
-## @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} = @var{cubicwgt})
-## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1},
-## @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} =
-## @var{cubicwgt}, @var{winradius} = 1)
+## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} = @var{cubicwgt})
+## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} = @var{cubicwgt}, @var{winradius} = 1)
##
-## Return the coefficient of the wavelet correlation of time
-## series (@var{time1}, @var{mag1}) and (@var{time2}, @var{mag2}).
-## @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}.
+## Return the coefficient of the wavelet correlation of two complex-valued time
+## series at a given time and frequency. The windowing function applied by
+## default is cubicwgt, this can be changed by passing a different function
+## handle to @var{window}, while the radius applied is set by @var{winradius}.
+## Note that this will be most effective when both series have had their mean
+## value (if it is not zero) subtracted (and stored separately); this reduces
+## the constant-offset error further, and allows the functions to be compared on
+## their periodic features rather than their constant features.
##
## @seealso{lswaveletcoeff, lscomplexwavelet, lsrealwavelet}
##
## @end deftypefn
-## nucorrcoeff, computes a coefficient of the wavelet correlation of two time series
-
function coeff = lscorrcoeff (x1, y1, x2, y2, t, o, wgt = @cubicwgt, wgtrad = 1)
so = 0.05 * o;
Deleted: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -1,79 +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 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} {@var{transform} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{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}. Each complex-valued result is the
-## pair (c_o, s_o) defining the coefficients which best fit the
-## function y = c_o * cos(ot) + s_o * sin(ot) to the (@var{time}, @var{mag}) data.
-##
-## @seealso{lscomplex}
-## @end deftypefn
-
-function transform = lsreal (t, x, omegamax, ncoeff, noctave)
-
- ## FIXME : 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.
- k = n = length (t);
-
- transform = zeros (1, (noctave * ncoeff));
-
- od = 2 ^ (- 1 / ncoeff);
- o = omegamax;
- n1 = 1 / n;
-
- ncoeffp = ncoeff * 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 = n1 * sum ((cos (ot) - i * sin (ot)) .* x);
-
- ot *= 2;
-
- iota = n1 * sum (cos (ot) - i * sin (ot));
-
-
- transform(iter) = (2 * (conj (zeta) - (conj (iota) * zeta)) /
- (1 - (real (iota) ^ 2) - (imag (iota) ^ 2)));
-
- o *= od;
- endfor
-
-endfunction
-
-%!test
-%! maxfreq = 4 / ( 2 * pi );
-%! t = linspace(0,8);
-%! x = ( 2 .* sin ( maxfreq .* t ) +
-%! 3 .* sin ( (3/4) * maxfreq .* t ) -
-%! 0.5 .* sin ( (1/4) * maxfreq .* t ) -
-%! 0.2 .* cos ( maxfreq .* t ) +
-%! cos ( (1/4) * maxfreq .* t ) );
-%! # In the assert here, I've got an error bound large enough to catch
-%! # individual system errors which would present no real issue.
-%! assert (lsreal (t,x,maxfreq,2,2),
-%! [(-1.68275915310663 + 4.70126183846743i), ...
-%! (1.93821553170889 + 4.95660209883437i), ...
-%! (4.38145452686697 + 2.14403733658600i), ...
-%! (5.27425332281147 - 0.73933440226597i)],
-%! 5e-10)
Added: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -0,0 +1,73 @@
+## 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} {@var{t} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
+##
+## Return a series of least-squares transforms of a real-valued time series.
+## Each transform is minimized independently for each frequency. The method
+## used is a Lomb-Scargle transform of the real-valued (@var{time}, @var{mag})
+## series, starting from frequency @var{maxfreq} and descending @var{numoctaves}
+## octaves with @var{numcoeff} coefficients per octave.
+##
+## The result of the transform for each frequency is the coefficient of a sum of
+## sine and cosine functions modified by that frequency, in the form of a
+## complex number—where the cosine coefficient is encoded in the real term, and
+## the sine coefficient is encoded in the imaginary term. Each frequency is fit
+## independently from the others, and to minimize very low frequency error,
+## consider storing the mean of a dataset with a constant or near-constant
+## offset separately, and subtracting it from the dataset.
+##
+## @seealso{lscomplex}
+## @end deftypefn
+
+
+
+function transform = lsreal (t, x, omegamax, ncoeff, noctave)
+
+ n = numel (t);
+
+ iter = 0 : (ncoeff * noctave - 1);
+ omul = (2 .^ (- iter / ncoeff));
+
+ ## For a given frequency, the iota term is taken at twice the frequency of the
+ ## zeta term.
+ ot = t(:) * (omul * omegamax);
+ oit = t(:) * (omul * omegamax * 2);
+
+ zeta = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
+ iota = sum ((cos (oit) - (sin (oit) .* i)), 1) / n;
+
+ transform = 2 .* (conj (zeta) - conj (iota) .* zeta) ./ (1 - abs (iota) .^ 2);
+
+endfunction
+
+%!test
+%! maxfreq = 4 / ( 2 * pi );
+%! t = linspace(0,8);
+%! x = ( 2 .* sin ( maxfreq .* t ) +
+%! 3 .* sin ( (3/4) * maxfreq .* t ) -
+%! 0.5 .* sin ( (1/4) * maxfreq .* t ) -
+%! 0.2 .* cos ( maxfreq .* t ) +
+%! cos ( (1/4) * maxfreq .* t ) );
+%! # In the assert here, I've got an error bound large enough to catch
+%! # individual system errors which would present no real issue.
+%! assert (lsreal (t,x,maxfreq,2,2),
+%! [(-1.68275915310663 + 4.70126183846743i), ...
+%! (1.93821553170889 + 4.95660209883437i), ...
+%! (4.38145452686697 + 2.14403733658600i), ...
+%! (5.27425332281147 - 0.73933440226597i)],
+%! 5e-10)
+
Copied: trunk/octave-forge/extra/lssa/inst/lsreal.m.old (from rev 10856, trunk/octave-forge/extra/lssa/inst/lsreal.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m.old (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m.old 2012-08-14 13:42:51 UTC (rev 10860)
@@ -0,0 +1,79 @@
+## 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 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} {@var{transform} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{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}. Each complex-valued result is the
+## pair (c_o, s_o) defining the coefficients which best fit the
+## function y = c_o * cos(ot) + s_o * sin(ot) to the (@var{time}, @var{mag}) data.
+##
+## @seealso{lscomplex}
+## @end deftypefn
+
+function transform = lsreal (t, x, omegamax, ncoeff, noctave)
+
+ ## FIXME : 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.
+ k = n = length (t);
+
+ transform = zeros (1, (noctave * ncoeff));
+
+ od = 2 ^ (- 1 / ncoeff);
+ o = omegamax;
+ n1 = 1 / n;
+
+ ncoeffp = ncoeff * 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 = n1 * sum ((cos (ot) - i * sin (ot)) .* x);
+
+ ot *= 2;
+
+ iota = n1 * sum (cos (ot) - i * sin (ot));
+
+
+ transform(iter) = (2 * (conj (zeta) - (conj (iota) * zeta)) /
+ (1 - (real (iota) ^ 2) - (imag (iota) ^ 2)));
+
+ o *= od;
+ endfor
+
+endfunction
+
+%!test
+%! maxfreq = 4 / ( 2 * pi );
+%! t = linspace(0,8);
+%! x = ( 2 .* sin ( maxfreq .* t ) +
+%! 3 .* sin ( (3/4) * maxfreq .* t ) -
+%! 0.5 .* sin ( (1/4) * maxfreq .* t ) -
+%! 0.2 .* cos ( maxfreq .* t ) +
+%! cos ( (1/4) * maxfreq .* t ) );
+%! # In the assert here, I've got an error bound large enough to catch
+%! # individual system errors which would present no real issue.
+%! assert (lsreal (t,x,maxfreq,2,2),
+%! [(-1.68275915310663 + 4.70126183846743i), ...
+%! (1.93821553170889 + 4.95660209883437i), ...
+%! (4.38145452686697 + 2.14403733658600i), ...
+%! (5.27425332281147 - 0.73933440226597i)],
+%! 5e-10)
Modified: trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-13 02:06:19 UTC (rev 10859)
+++ trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-14 13:42:51 UTC (rev 10860)
@@ -18,13 +18,17 @@
## @deftypefnx {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq}, @var{window}=cubicwgt)
## @deftypefnx {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq}, @var{window}=cubicwgt, @var{winradius}=1)
##
-## Return the coefficient of the wavelet transform of the
-## complex time series (@var{t}, @var{x}) 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.)
+## Return the wavelet transform of a complex time series in a given window. The
+## transform takes a complex time series (@var{t}, @var{x}) at time @var{time}
+## and frequency @var{freq}, then applies a windowing function to it; the
+## default is cubicwgt, however by providing a function handle for the optional
+## variable @var{window}, the user may select their own function; to determine
+## the radius of the interval around the @var{time} selected, set
+## @var{winradius} to some value other than 1.
##
+## This transform operates identically to the transform at the heart of
+## lscomplexwavelet, however for one window only.
+##
## @seealso{lscorrcoeff, lscomplexwavelet, lsrealwavelet}
##
## @end deftypefn
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
×
Want the latest updates on software, tech news, and AI?
Get latest updates about software, tech news, and AI from SourceForge directly in your inbox once a month.
Thanks for helping keep SourceForge clean.
X