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