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