Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10851] trunk/octave-forge/extra/lssa
|
From: <be...@us...> - 2012-08-09 23:49:28
|
Revision: 10851
http://octave.svn.sourceforge.net/octave/?rev=10851&view=rev
Author: benjf5
Date: 2012-08-09 23:49:21 +0000 (Thu, 09 Aug 2012)
Log Message:
-----------
Fixed INDEX, NEWS for lssa 0.1.1 release. Also updated formatting.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/DESCRIPTION
trunk/octave-forge/extra/lssa/INDEX
trunk/octave-forge/extra/lssa/NEWS
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/lscomplex.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
trunk/octave-forge/extra/lssa/inst/lsreal.m
trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
Modified: trunk/octave-forge/extra/lssa/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/lssa/DESCRIPTION 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/DESCRIPTION 2012-08-09 23:49:21 UTC (rev 10851)
@@ -5,9 +5,9 @@
Maintainer: Ben Lewis <be...@gm...>
Title: Least squares spectral analysis
Description: A package implementing tools to compute spectral decompositions of
-irregularly-spaced time series. Currently includes functions based off the
-Lomb-Scargle periodogram and Adolf Mathias' implementation for R and C (see
-URLs).
+ irregularly-spaced time series. Currently includes functions based off the
+ Lomb-Scargle periodogram and Adolf Mathias' implementation for R and C (see
+ URLs).
Url: http://www.jstatsoft.org/v11/i02
Problems: fast implementations, wavelet functions are currently not functional.
Depends: octave (>= 3.4.3)
Modified: trunk/octave-forge/extra/lssa/INDEX
===================================================================
--- trunk/octave-forge/extra/lssa/INDEX 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/INDEX 2012-08-09 23:49:21 UTC (rev 10851)
@@ -1,10 +1,19 @@
-cubicwgt
-fastlscomplex
-lombcoeff
-lombnormcoeff
-lscomplex
-lscomplexwavelet
-lscorrcoeff
-lsreal
-lsrealwavelet
-lswaveletcoeff
\ No newline at end of file
+lssa >> Least Squares Spectral Analysis
+Windowing
+ cubicwgt
+Periodogram
+ lombcoeff lombnormcoeff
+Accelerated time-series functions
+ fastlscomplex
+Complex time-series functions
+ lscomplex fastlscomplex lswaveletcoeff
+# lscomplexwavelet
+Real time-series functions
+ lsreal
+# lsrealwavelet
+Correlation
+ lscorrcoeff
+Wavelet Transform
+ lswaveletcoeff
+# lscomplexwavelet lsrealwavelet
+## The wavelet functions are unavailable until I can get them working.
Modified: trunk/octave-forge/extra/lssa/NEWS
===================================================================
--- trunk/octave-forge/extra/lssa/NEWS 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/NEWS 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,9 +14,13 @@
author's reasoning. Its generated results are relevant to any given data
set, but it does not appear to be normalized to any great extent.
- ** lsrealwavelet and lscomplexwavelet are not currently correct or even really
- functional; this is an ongoing struggle for me, and if anyone can suggest
- what I might do to improve them, please! By all means, send me an email.
+ ** There are two wavelet functions under development, but they are not included
+ in this release as they are currently not functional. For all your wavelet
+ needs, the specific transformation used is available in the lswaveletcoeff
+ function, and will generate a single cosine/sine magnitude pair (as a
+ complex number) for a complex-valued series (this function may be joined by
+ a companion for real-valued series) and can be looped to simulate a full
+ wavelet transform.
** For all the working functions, tests have been written and formatted to
Octave coding standards. These tests should pass on any given architecture
Modified: trunk/octave-forge/extra/lssa/inst/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,8 +14,8 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*-texinfo-*-
-## @deftypefn {Function File} {a =} cubicwgt {series}
-## Return the series as windowed by a cubic polynomial,
+## @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
@@ -36,4 +36,4 @@
%!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.
\ No newline at end of file
+%! ## on any array input.
Modified: trunk/octave-forge/extra/lssa/inst/lombcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,7 +14,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {c =} lombcoeff (time, mag, freq)
+## @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.
Modified: trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,7 +14,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {c =} lombnormcoeff (time, mag, freq)
+## @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
@@ -25,11 +25,17 @@
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) +
- sum (X .* sin (omega .* T - tau)) .^ 2 / sum (sin (omega .* T - tau) .^ 2 )) /
- (2 * var(X)));
+ xmean = mean (X);
+ theta = atan2 (sum (sin (2 .* omega .*T)),
+ sum (cos (2 .* omega .* T))) / (2*omega);
+
+ coeff = ((sum ((X-xmean) .* cos (omega .* T - theta)) .^ 2 /
+ sum (cos (omega .* T - theta) .^ 2) +
+ sum ((X-xmean) .* sin (omega .* T - theta)) .^ 2 /
+ sum (sin (omega .* T - theta) .^ 2 )) /
+ (2 * var(X)));
+
endfunction
%!shared t, x, o, maxfreq
@@ -41,6 +47,6 @@
%! 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);
+%!assert (lombnormcoeff (t,x,o(1)), 44.7068607258824, 5e-10);
+%!assert (lombnormcoeff (t,x,o(2)), 35.7769955188467, 5e-10);
+%!assert (lombnormcoeff (t,x,o(3)), 20.7577786183241, 5e-10);
Modified: trunk/octave-forge/extra/lssa/inst/lscomplex.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,7 +14,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {t =} lscomplex ( time, mag, maxfreq, numcoeff, numoctaves)
+## @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}
@@ -25,14 +25,25 @@
function transform = lscomplex( t , x , omegamax , ncoeff , noctave )
- n = length(t); ## VECTOR ONLY, and since t and x have the same number of entries, there's no problem.
+
+ n = length(t); ## VECTOR ONLY, and since t and x have the same number of
+ ## entries, there's no problem.
+
transform = zeros(1,ncoeff*noctave);
+
o = omegamax;
- omul = 2 ^ ( - 1 / ncoeff );
+
+ omul = 2 ^ (- 1 / ncoeff);
+
for iter = 1:ncoeff*noctave
+
ot = o .* t;
- transform(iter) = sum( ( cos(ot) - ( sin(ot) .* i ) ) .* x ) / n; ## See the paper for the expression
+
+ transform(iter) = sum ((cos (ot) - (sin (ot) .* i)) .* x) / n; ## See the
+ ## paper for the expression
+
o *= omul; ## To advance the transform to the next coefficient in the octave
+
endfor
endfunction
@@ -48,4 +59,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 );
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,12 +14,15 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {c =} lscorrcoeff (abc1, ord1, abc2, ord2, time, freq)
-## @deftypefnx {Function File} {c =} lscorrcoeff (abc1, ord1, abc2, ord2, time, freq, window)
-## @deftypefnx {Function File} {c =} lscorrcoeff (abc1, ord1, abc2, ord2, time, freq, window, winradius)
+## @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)
##
## Return the coefficient of the wavelet correlation of time
-## series (@var{abc1}, @var{ord1}) and (@var{abc2}, @var{ord2}).
+## 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}.
Modified: trunk/octave-forge/extra/lssa/inst/lsreal.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lsreal.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,7 +14,7 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {transform =} lsreal ( time, mag, maxfreq, numcoeff, numoctaves)
+## @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
@@ -26,35 +26,39 @@
## @seealso{lscomplex}
## @end deftypefn
-function transform = lsreal( 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.
+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)); ## 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?
+ transform = zeros(1,(noctave * ncoeff));
+
+ od = 2 ^ (- 1 / ncoeff);
+
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 ...
+
+ n1 = 1 / n;
+
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 ) );
+
+ zeta = sum ((cos (ot) .* x) -
+ (sin (ot) .* x .* i));
+
ot = ot .* 2;
- iota = sum( cos(ot) - ( sin(ot) .* i ) );
+
+ 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 ));
+
+ transform(iter) = (2 / (1 -(real (iota) ^ 2) - (imag(iota) ^ 2)) * (conj(zeta) - (conj(iota) * zeta)));
o = o .* od;
endfor
@@ -73,4 +77,4 @@
%! [-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
+%! # individual system errors which would present no real issue.
Modified: trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-09 16:45:30 UTC (rev 10850)
+++ trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-09 23:49:21 UTC (rev 10851)
@@ -14,12 +14,12 @@
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {c =} lswaveletcoeff (abc, ord, time, freq)
-## @deftypefnx {Function File} {c =} lswaveletcoeff (abc, ord, time, freq, window)
-## @deftypefnx {Function File} {c =} lswaveletcoeff (abc, ord, time, freq, window, winradius)
+## @deftypefn {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq})
+## @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{abc}, @var{ord}) at time @var{time}
+## 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
@@ -74,4 +74,4 @@
%!assert (lswaveletcoeff (t, x, 0.5, maxfreq), 0.383340407638780 +
%! 2.385251997545446i , 5e-10);
%!assert (lswaveletcoeff (t, y, 3.3, 3/4 * maxfreq), -2.35465091096084 +
-%! 1.01892561714824i, 5e-10);
\ No newline at end of file
+%! 1.01892561714824i, 5e-10);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
