Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10852] trunk/octave-forge/extra/lssa/inst
|
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.
|
