Menu
▾
▴
octave-cvsupdate — CVS update notifications
This list is closed, nobody may subscribe to it.
| 2001 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(10) |
Aug
(5) |
Sep
(3) |
Oct
(41) |
Nov
(41) |
Dec
(33) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2002 |
Jan
(75) |
Feb
(10) |
Mar
(170) |
Apr
(174) |
May
(66) |
Jun
(11) |
Jul
(10) |
Aug
(44) |
Sep
(73) |
Oct
(28) |
Nov
(139) |
Dec
(52) |
| 2003 |
Jan
(35) |
Feb
(93) |
Mar
(62) |
Apr
(10) |
May
(55) |
Jun
(70) |
Jul
(37) |
Aug
(16) |
Sep
(56) |
Oct
(31) |
Nov
(57) |
Dec
(83) |
| 2004 |
Jan
(85) |
Feb
(67) |
Mar
(27) |
Apr
(37) |
May
(75) |
Jun
(85) |
Jul
(160) |
Aug
(68) |
Sep
(104) |
Oct
(25) |
Nov
(39) |
Dec
(23) |
| 2005 |
Jan
(10) |
Feb
(45) |
Mar
(43) |
Apr
(19) |
May
(108) |
Jun
(31) |
Jul
(41) |
Aug
(23) |
Sep
(65) |
Oct
(58) |
Nov
(44) |
Dec
(54) |
| 2006 |
Jan
(96) |
Feb
(27) |
Mar
(69) |
Apr
(59) |
May
(67) |
Jun
(35) |
Jul
(13) |
Aug
(461) |
Sep
(160) |
Oct
(399) |
Nov
(32) |
Dec
(72) |
| 2007 |
Jan
(316) |
Feb
(305) |
Mar
(318) |
Apr
(54) |
May
(194) |
Jun
(173) |
Jul
(282) |
Aug
(91) |
Sep
(227) |
Oct
(365) |
Nov
(168) |
Dec
(18) |
| 2008 |
Jan
(71) |
Feb
(111) |
Mar
(155) |
Apr
(173) |
May
(70) |
Jun
(67) |
Jul
(55) |
Aug
(83) |
Sep
(32) |
Oct
(68) |
Nov
(80) |
Dec
(29) |
| 2009 |
Jan
(46) |
Feb
(18) |
Mar
(95) |
Apr
(76) |
May
(140) |
Jun
(98) |
Jul
(84) |
Aug
(123) |
Sep
(94) |
Oct
(131) |
Nov
(142) |
Dec
(125) |
| 2010 |
Jan
(128) |
Feb
(158) |
Mar
(172) |
Apr
(134) |
May
(94) |
Jun
(84) |
Jul
(32) |
Aug
(127) |
Sep
(167) |
Oct
(109) |
Nov
(69) |
Dec
(78) |
| 2011 |
Jan
(39) |
Feb
(58) |
Mar
(52) |
Apr
(47) |
May
(56) |
Jun
(76) |
Jul
(55) |
Aug
(54) |
Sep
(165) |
Oct
(255) |
Nov
(328) |
Dec
(263) |
| 2012 |
Jan
(82) |
Feb
(147) |
Mar
(400) |
Apr
(216) |
May
(209) |
Jun
(160) |
Jul
(86) |
Aug
(141) |
Sep
(156) |
Oct
(6) |
Nov
|
Dec
|
| 2015 |
Jan
|
Feb
|
Mar
|
Apr
|
May
|
Jun
|
Jul
(1) |
Aug
|
Sep
(1) |
Oct
|
Nov
(1) |
Dec
(2) |
| 2016 |
Jan
|
Feb
(2) |
Mar
(2) |
Apr
(1) |
May
(1) |
Jun
(2) |
Jul
(1) |
Aug
(1) |
Sep
|
Oct
|
Nov
(1) |
Dec
|
| 2019 |
Jan
|
Feb
|
Mar
(1) |
Apr
|
May
|
Jun
|
Jul
|
Aug
|
Sep
|
Oct
|
Nov
|
Dec
|
| 2021 |
Jan
|
Feb
|
Mar
|
Apr
(3) |
May
(4) |
Jun
(8) |
Jul
(2) |
Aug
(5) |
Sep
(9) |
Oct
|
Nov
|
Dec
|
|
From: <par...@us...> - 2012-08-01 19:44:39
|
Revision: 10798
http://octave.svn.sourceforge.net/octave/?rev=10798&view=rev
Author: paramaniac
Date: 2012-08-01 19:44:32 +0000 (Wed, 01 Aug 2012)
Log Message:
-----------
control-devel: show singular value plot if no return arguments are specified
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
trunk/octave-forge/extra/control-devel/inst/moen4.m
trunk/octave-forge/extra/control-devel/inst/moesp.m
trunk/octave-forge/extra/control-devel/inst/n4sid.m
Modified: trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-08-01 18:53:17 UTC (rev 10797)
+++ trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-08-01 19:44:32 UTC (rev 10798)
@@ -24,7 +24,7 @@
## Created: May 2012
## Version: 0.1
-function [sys, x0, info] = __slicot_identification__ (method, dat, varargin)
+function [sys, x0, info] = __slicot_identification__ (method, nout, dat, varargin)
## determine identification method
switch (method)
@@ -42,7 +42,7 @@
error ("%s: first argument must be a time-domain 'iddata' dataset", method);
endif
- if (nargin > 2) # ident (dat, ...)
+ if (nargin > 3) # ident (dat, ...)
if (is_real_scalar (varargin{1})) # ident (dat, n, ...)
varargin = horzcat (varargin(2:end), {"order"}, varargin(1));
endif
@@ -143,53 +143,72 @@
if (! isempty (conf))
ctrl = ! conf;
endif
-
- ## perform system identification
- [a, b, c, d, q, ry, s, k, x0] = slident (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
- ## compute noise variance matrix factor L
- ## L L' = Ry, e = L v
- ## v becomes white noise with identity covariance matrix
- l = chol (ry, "lower");
+ if (nout == 0)
+ ## compute singular values
+ [sv, nrec] = slib01ad (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
- ## assemble model
- [inname, outname] = get (dat, "inname", "outname");
- if (strncmpi (noise, "e", 1)) # add error inputs e, not normalized
- sys = ss (a, [b, k], c, [d, eye(p)], tsam);
- in_u = __labels__ (inname, "u");
- in_e = __labels__ (outname, "y");
- in_e = cellfun (@(x) ["e@", x], in_e, "uniformoutput", false);
- inname = [in_u; in_e];
- elseif (strncmpi (noise, "v", 1)) # add error inputs v, normalized
- sys = ss (a, [b, k*l], c, [d, l], tsam);
- in_u = __labels__ (inname, "u");
- in_v = __labels__ (outname, "y");
- in_v = cellfun (@(x) ["v@", x], in_v, "uniformoutput", false);
- inname = [in_u; in_v];
- elseif (strncmpi (noise, "k", 1)) # Kalman predictor
- sys = ss ([a-k*c], [b-k*d, k], c, [d, zeros(p)], tsam);
- in_u = __labels__ (inname, "u");
- in_y = __labels__ (outname, "y");
- inname = [in_u; in_y];
- else # no error inputs, default
- sys = ss (a, b, c, d, tsam);
- endif
+ ## there is no 'logbar' function
+ svl = log10 (sv);
+ base = floor (min (svl));
- sys = set (sys, "inname", inname, "outname", outname);
+ clf
+ bar (svl, "basevalue", base)
+ xlim ([0, length(sv)+1])
+ yl = ylim;
+ ylim ([base, yl(2)])
+ title ("Singular Values")
+ ylabel ("Logarithm of Singular Values")
+ xlabel (sprintf ("Estimated System Order with current Tolerance: %d", nrec))
+ grid on
+ else
+ ## perform system identification
+ [a, b, c, d, q, ry, s, k, x0] = slident (dat.y, dat.u, nobr, n, meth, alg, conct, ctrl, rcond, tol);
- ## return x0 as vector for single-experiment data
- ## instead of a cell containing one vector
- if (numel (x0) == 1)
- x0 = x0{1};
+ ## compute noise variance matrix factor L
+ ## L L' = Ry, e = L v
+ ## v becomes white noise with identity covariance matrix
+ l = chol (ry, "lower");
+
+ ## assemble model
+ [inname, outname] = get (dat, "inname", "outname");
+ if (strncmpi (noise, "e", 1)) # add error inputs e, not normalized
+ sys = ss (a, [b, k], c, [d, eye(p)], tsam);
+ in_u = __labels__ (inname, "u");
+ in_e = __labels__ (outname, "y");
+ in_e = cellfun (@(x) ["e@", x], in_e, "uniformoutput", false);
+ inname = [in_u; in_e];
+ elseif (strncmpi (noise, "v", 1)) # add error inputs v, normalized
+ sys = ss (a, [b, k*l], c, [d, l], tsam);
+ in_u = __labels__ (inname, "u");
+ in_v = __labels__ (outname, "y");
+ in_v = cellfun (@(x) ["v@", x], in_v, "uniformoutput", false);
+ inname = [in_u; in_v];
+ elseif (strncmpi (noise, "k", 1)) # Kalman predictor
+ sys = ss ([a-k*c], [b-k*d, k], c, [d, zeros(p)], tsam);
+ in_u = __labels__ (inname, "u");
+ in_y = __labels__ (outname, "y");
+ inname = [in_u; in_y];
+ else # no error inputs, default
+ sys = ss (a, b, c, d, tsam);
+ endif
+
+ sys = set (sys, "inname", inname, "outname", outname);
+
+ ## return x0 as vector for single-experiment data
+ ## instead of a cell containing one vector
+ if (numel (x0) == 1)
+ x0 = x0{1};
+ endif
+
+ ## assemble info struct
+ ## Kalman gain matrix K
+ ## state covariance matrix Q
+ ## output covariance matrix Ry
+ ## state-output cross-covariance matrix S
+ ## noise variance matrix factor L
+ info = struct ("K", k, "Q", q, "Ry", ry, "S", s, "L", l);
endif
-
- ## assemble info struct
- ## Kalman gain matrix K
- ## state covariance matrix Q
- ## output covariance matrix Ry
- ## state-output cross-covariance matrix S
- ## noise variance matrix factor L
- info = struct ("K", k, "Q", q, "Ry", ry, "S", s, "L", l);
endfunction
Modified: trunk/octave-forge/extra/control-devel/inst/moen4.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/moen4.m 2012-08-01 18:53:17 UTC (rev 10797)
+++ trunk/octave-forge/extra/control-devel/inst/moen4.m 2012-08-01 19:44:32 UTC (rev 10798)
@@ -232,7 +232,11 @@
print_usage ();
endif
- [sys, x0, info] = __slicot_identification__ ("moen4", varargin{:});
+ if (nargout == 0)
+ __slicot_identification__ ("moen4", nargout, varargin{:});
+ else
+ [sys, x0, info] = __slicot_identification__ ("moen4", nargout, varargin{:});
+ endif
endfunction
Modified: trunk/octave-forge/extra/control-devel/inst/moesp.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/moesp.m 2012-08-01 18:53:17 UTC (rev 10797)
+++ trunk/octave-forge/extra/control-devel/inst/moesp.m 2012-08-01 19:44:32 UTC (rev 10798)
@@ -34,6 +34,10 @@
print_usage ();
endif
- [sys, x0, info] = __slicot_identification__ ("moesp", varargin{:});
+ if (nargout == 0)
+ __slicot_identification__ ("moesp", nargout, varargin{:});
+ else
+ [sys, x0, info] = __slicot_identification__ ("moesp", nargout, varargin{:});
+ endif
-endfunction
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/extra/control-devel/inst/n4sid.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/n4sid.m 2012-08-01 18:53:17 UTC (rev 10797)
+++ trunk/octave-forge/extra/control-devel/inst/n4sid.m 2012-08-01 19:44:32 UTC (rev 10798)
@@ -34,6 +34,10 @@
print_usage ();
endif
- [sys, x0, info] = __slicot_identification__ ("n4sid", varargin{:});
+ if (nargout == 0)
+ __slicot_identification__ ("n4sid", nargout, varargin{:});
+ else
+ [sys, x0, info] = __slicot_identification__ ("n4sid", nargout, varargin{:});
+ endif
-endfunction
\ No newline at end of file
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-08-01 18:53:27
|
Revision: 10797
http://octave.svn.sourceforge.net/octave/?rev=10797&view=rev
Author: benjf5
Date: 2012-08-01 18:53:17 +0000 (Wed, 01 Aug 2012)
Log Message:
-----------
Fixing a few errors for packaging.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/INDEX
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/lscomplexwavelet.m
trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m
trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
trunk/octave-forge/extra/lssa/src/fastlscomplex.cc
trunk/octave-forge/extra/lssa/src/fastlsreal.cc
Modified: trunk/octave-forge/extra/lssa/INDEX
===================================================================
--- trunk/octave-forge/extra/lssa/INDEX 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/INDEX 2012-08-01 18:53:17 UTC (rev 10797)
@@ -1,7 +1,7 @@
cubicwgt
fastlscomplex
lombcoeff
-lomgnormcoeff
+lombnormcoeff
lscomplex
lscomplexwavelet
lscorrcoeff
Modified: trunk/octave-forge/extra/lssa/inst/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -2,7 +2,7 @@
##
## 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
+## Foundation; either version 2 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
Modified: trunk/octave-forge/extra/lssa/inst/lombcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -36,5 +36,5 @@
function coeff = lombcoeff(T, X, o)
theta = atan2(sum(sin(2 .* o .* T )), sum(cos(2.*o.*T)))/ (2 * o );
coeff = ( sum(X .* cos(o .* T - theta))**2)/(sum(cos(o.*T-theta).**2)) + ( sum(X .* sin(o .* T - theta))**2)/(sum(sin(o.*T-theta).**2));
-end function
+endfunction
Modified: trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -27,4 +27,4 @@
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) ) );
-end
\ No newline at end of file
+endfunction
Modified: trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -1,7 +1,5 @@
## Copyright (C) 2012 Ben Lewis <be...@gm...>
##
-## This code is part of GNU Octave.
-##
## This software 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 2 of the License, or (at your option) any later
Modified: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -2,7 +2,7 @@
##
## 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
+## Foundation; either version 2 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
@@ -24,6 +24,8 @@
## default is cubicwgt if left blank, and its radius is 1,
## as defined as the default for @var{winradius}.
##
+## @seealso lswaveletcoeff lscomplexwavelet lsrealwavelet
+##
## @end deftypefn
## Demo with sin, cos as Nir suggested.
Modified: trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lsrealwavelet.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -1,9 +1,30 @@
-# Copyright (C) 2012 Benjamin Lewis <be...@gm...>
+## 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 2 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 -*-
-## @deffun {Function File} {t =}
+## @deffun {Function File} {t =} lsrealwavelet( @var{time}, @var{mag},
+## @var{maxfreq}, @var{coefficients}, @var{octaves}, @var{time_min},
+## @var{time_max}, @var{min_window_count} )
##
+## Computes a windowed transform of the supplied (@var{time}, @var{mag}) series
+## of real-valued magnitudes, applying progressively wider windows as the
+## frequencies tested decline from the maximum frequency.
+##
+## Currently non-functional.
##
+## @seealso lscomplexwavelet lswaveletcoeff lscorrcoeff
##
## @end deffun
@@ -65,6 +86,4 @@
o *= omegamult;
endfor
-
-
endfunction
Modified: trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m 2012-08-01 18:53:17 UTC (rev 10797)
@@ -25,6 +25,8 @@
## while @var{winradius} is the windowing radius, and defaults
## to 1 (the radius of cubicwgt.)
##
+## @seealso lscorrcoeff lscomplexwavelet lsrealwavelet
+##
## @end deftypefn
%!demo
Modified: trunk/octave-forge/extra/lssa/src/fastlscomplex.cc
===================================================================
--- trunk/octave-forge/extra/lssa/src/fastlscomplex.cc 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/src/fastlscomplex.cc 2012-08-01 18:53:17 UTC (rev 10797)
@@ -1,5 +1,18 @@
/* Copyright (C) 2012 Benjamin Lewis <be...@gm...>
- * GNU GPLv2
+ *
+ * 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 2 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, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
Modified: trunk/octave-forge/extra/lssa/src/fastlsreal.cc
===================================================================
--- trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-08-01 15:45:50 UTC (rev 10796)
+++ trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-08-01 18:53:17 UTC (rev 10797)
@@ -1,5 +1,18 @@
/* Copyright (C) 2012 Benjamin Lewis <be...@gm...>
- * Licensed under the GNU GPLv2
+ *
+ * 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 2 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, write to the Free Software Foundation, Inc.,
+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
@@ -285,13 +298,13 @@
// In this case, there is no next record, but this record has data.
if ( record_current->stored_data ) {
int p = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1 , temp_ptr_alpha = temp_alpha ; p < 12 ; p++ , exp_pse_ptr++ ) {
+ for( exp_pse_ptr = exp_power_series_elements , temp_ptr_alpha = temp_alpha ; ; ) {
tpra = temp_ptr_alpha;
temp_ptr_alpha->x = record_current->power_series[p]->x;
(temp_ptr_alpha++)->t = record_current->power_series[p]->t;
temp_ptr_beta->x = -record_current->power_series[p]->x;
(temp_ptr_beta++)->t = -record_current->power_series[p]->t;
- for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *temp_ptr_alpha * *exp_ptr; ; ) {
+ for( exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->t * *exp_ptr ; ; ) {
/* This next block is from Mathias' code, and it does a few
* ... unsavoury things. First off, it uses conditionals with
* break in order to avoid potentially accessing null regions
@@ -299,87 +312,82 @@
* numbers. However, remembering that most of these will not
* actually be accessed for the first iterations, it's easier.
*/
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = *tpra * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_current->power_series[p] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = -*tpra * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_current->power_series[p] += *tpra * *exp_ptr;
+ if ( --exp_ptr < exp_power_series_elements ) break;
+ ++tpra;
+ record_current->power_series[p]->x -= tpra->x * *exp_ptr;
+ record_current->power_series[p]->t -= tpra->t * *exp_ptr;
+ if ( --exp_ptr < exp_power_series_elements ) break;
+ ++tpra;
+ record_current->power_series[p]->x += tpra->x * *exp_ptr;
+ record_current->power_series[p]->t += tpra->x * *exp_ptr;
}
+ if ( ++p >= 12 ) break;
+ temp_ptr_alpha->x = -record_current->power_series[p]->x;
+ (temp_ptr_alpha++)->t = -record_current->power_series[p]->t;
+ temp_ptr_beta->x = record_current->power_series[p]->x;
+ (temp_ptr_beta++)->t = record_current->power_series[p]->t;
+ for( tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->t = tprb->t * *exp_ptr; exp_ptr > exp_power_series_elements ; ) {
+ ++tprb;
+ --exp_ptr;
+ record_current->power_series[p]->t += tprb->t * *exp_ptr;
+ }
+ if ( ++p >= 12 ) break;
}
- if ( ! record_ref ) break; // Last record was reached
}
- else {
- record_next = record_ref;
- if ( record_current->stored_data ) {
- int p = 0, q = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; p < 12 ; p++, q++, exp_pse_ptr++ ) {
- tpra = temp_ptr_alpha;
- *temp_ptr_alpha++ = record_current->power_series[p] + record_next->power_series[q];
- *temp_ptr_beta++ = record_current->power_series[p] - record_next->power_series[1];
- tprb = temp_ptr_beta;
- for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *tpra * *exp_ptr; ; ) {
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tprb -= 2;
- h = *tprb * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tpra -= 2;
- record_current->power_series[p] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tprb -= 2;
- h = - *tprb * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tpra -= 2;
- record_current->power_series[p] += *tpra * *exp_ptr;
- }
+ if ( ! record_ref ) break; // Last record was reached
+ }
+ else {
+ record_next = record_ref;
+ if ( record_current->stored_data ) {
+ int p = 0;
+ for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) {
+ temp_ptr_alpha->x = record_current->power_series[p]->x + record_next->power_series[p]->x;
+ (temp_ptr_alpha++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t;
+ temp_ptr_beta->x = record_ref->power_series[p]->x - record_current->power_series[p]->x;
+ (temp_ptr_beta++)->t = record_ref->power_series[p]->t - record_current->power_series[p]->t;
+ for( tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tpra->x * *exp_ptr, record_current->power_series[p]->t = tpra->x * *exp_ptr; ; ) {
+ if ( --exp_ptr < exp_pse_ptr ) break;
+ ++tpra;
+ record_current->power_series[p]->x -= tpra->x * *exp_ptr;
+ record_current->power_series[p]->t -= tpra->t * *exp_ptr;
+ if ( --exp_ptr < exp_pse_ptr ) break;
+ ++tpra;
+ record_current->power_series[p]->x += tpra->x * *exp_ptr;
+ record_current->power_series[p]->t += tpra->t * *exp_ptr;
}
- } else {
- int q = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; q < 12; q++, exp_pse_ptr++ ) {
- tpra = temp_ptr_alpha;
- *temp_ptr_alpha++ = std::complex<double>(record_next->power_series[q]);
- for ( exp_ptr = exp_power_series_elements, record_next->power_series[q] = *tpra * *exp_ptr; ; ) {
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = *tpra * *exp_ptr;
- record_next->power_series[q].real() -= h.imag();
- record_next->power_series[q].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_next->power_series[q] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = -*tpra * *exp_ptr;
- record_next->power_series[q].real() -= h.imag();
- record_next->power_series[q].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_next->power_series[q] += *tpra * *exp_ptr;
- }
+ if ( ++p >= 12 ) break;
+ temp_ptr_alpha->x = record_next->power_series[p]->x - record_current->power_series[p]->x;
+ (temp_ptr_alpha++)->t = record_next->power_series[p]->t - record_current->power_series[p]->t;
+ temp_ptr_beta->x = record_current->power_series[p]->x + record_next->power_series[p]->x;
+ (temp_ptr_beta++)->t = record_current->power_series[p]->t + record_next->power_series[p]->t;
+ for(tprb = temp_beta, exp_ptr = exp_pse_ptr++, record_current->power_series[p]->x = tprb->x * *exp_ptr, record_current->power_series[p]->t = tprb->x * *exp_ptr; exp_ptr > exp_power_series_elements; ) {
+ ++tprb;
+ --exp_ptr;
+ record_current->power_series[p]->x += tprb->x * *exp_ptr;
+ record_current->power_series[p]->t += tprb->t * *exp_ptr;
}
+ if ( ++p >= 12 ) break;
}
- record_current->stored_data = true;
- record_ref = record_next;
- record_current->next = record_ref->next;
- delete record_ref;
- }
+ } else {
+ int q = 0;
+ for( exp_pse_ptr = exp_power_series_elements, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; ; ) {
+ temp_ptr_alpha->x = record_next->power_series[q]->x;
+ temp_ptr_alpha->t = record_next->power_series[q]->t;
+ for(tpra = temp_alpha, exp_ptr = exp_pse_ptr++, record_next->power_series[q]->x = tpra->x * *exp_ptr, record_next->power_series[q]->t = tpra->t * *exp_ptr; exp_ptr > exp_power_series_elements; ) {
+ ++tpra;
+ --exp_ptr;
+ record_next->power_series[q]->x += tpra->x * *exp_ptr;
+ record_next->power_series[q]->t += tpra->t * *exp_ptr;
+ }
+ if ( ++q >= 12 ) break;
+ }
+ record_current->stored_data = true;
+ record_ref = record_next;
+ record_current->next = record_ref->next;
+ record_next = 0;
+ delete record_ref;
}
}
-
-
+ }
return results;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-08-01 15:45:56
|
Revision: 10796
http://octave.svn.sourceforge.net/octave/?rev=10796&view=rev
Author: paramaniac
Date: 2012-08-01 15:45:50 +0000 (Wed, 01 Aug 2012)
Log Message:
-----------
control-devel: improve argument checking
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
Modified: trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m
===================================================================
--- trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-07-31 15:34:36 UTC (rev 10795)
+++ trunk/octave-forge/extra/control-devel/inst/__slicot_identification__.m 2012-08-01 15:45:50 UTC (rev 10796)
@@ -72,7 +72,7 @@
conct = 1; # no connection between experiments
ctrl = 1; # don't confirm order n
rcond = 0.0;
- tol = -1.0; % 0;
+ tol = -1.0;
s = [];
n = [];
conf = [];
@@ -83,16 +83,27 @@
key = lower (varargin{k});
val = varargin{k+1};
switch (key)
- ## TODO: proper argument checking
case {"n", "order"}
+ if (! issample (val, 0) || val != round (val))
+ error ("%s: 'n' must be a positive integer", method);
+ endif
n = val;
case "s"
+ if (! issample (val, 0) || val != round (val))
+ error ("%s: 's' must be a positive integer", method);
+ endif
s = val;
case {"alg", "algorithm"}
error ("alg");
case "tol"
+ if (! is_real_scalar (val))
+ error ("%s: tolerance 'tol' must be a real scalar", method);
+ endif
tol = val;
case "rcond"
+ if (! is_real_scalar (val))
+ error ("%s: 'rcond' must be a real scalar", method);
+ endif
rcond = val;
case "confirm"
conf = logical (val);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-07-31 15:34:47
|
Revision: 10795
http://octave.svn.sourceforge.net/octave/?rev=10795&view=rev
Author: benjf5
Date: 2012-07-31 15:34:36 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
Ready for an alpha release. The data files will be improved shortly.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/DESCRIPTION
trunk/octave-forge/extra/lssa/src/Makefile
Added Paths:
-----------
trunk/octave-forge/extra/lssa/INDEX
trunk/octave-forge/extra/lssa/NEWS
trunk/octave-forge/extra/lssa/samples/
trunk/octave-forge/extra/lssa/samples/SampleDataSet.m
trunk/octave-forge/extra/lssa/samples/SampleScriptWithVostokData
Removed Paths:
-------------
trunk/octave-forge/extra/lssa/SampleScriptWithVostokData
trunk/octave-forge/extra/lssa/inst/SampleDataSet.m
Modified: trunk/octave-forge/extra/lssa/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/lssa/DESCRIPTION 2012-07-31 14:43:11 UTC (rev 10794)
+++ trunk/octave-forge/extra/lssa/DESCRIPTION 2012-07-31 15:34:36 UTC (rev 10795)
@@ -1,5 +1,5 @@
Name: lssa
-Version:
+Version: 0.1.0
Date: 2012-07-28
Author: Ben Lewis
Maintainer: Ben Lewis
@@ -10,5 +10,6 @@
URLs).
Url: http://www.jstatsoft.org/v11/i02
Problems: fast implementations, wavelet functions are currently not functional.
-Depends:
-License: GPL version 2 or later
+Depends: octave (>= 3.4.3)
+Autoload: no
+License: GPLv2+
Added: trunk/octave-forge/extra/lssa/INDEX
===================================================================
--- trunk/octave-forge/extra/lssa/INDEX (rev 0)
+++ trunk/octave-forge/extra/lssa/INDEX 2012-07-31 15:34:36 UTC (rev 10795)
@@ -0,0 +1,10 @@
+cubicwgt
+fastlscomplex
+lombcoeff
+lomgnormcoeff
+lscomplex
+lscomplexwavelet
+lscorrcoeff
+lsreal
+lsrealwavelet
+lswaveletcoeff
\ No newline at end of file
Added: trunk/octave-forge/extra/lssa/NEWS
===================================================================
--- trunk/octave-forge/extra/lssa/NEWS (rev 0)
+++ trunk/octave-forge/extra/lssa/NEWS 2012-07-31 15:34:36 UTC (rev 10795)
@@ -0,0 +1,20 @@
+Welcome to the first release of lssa, 0.1.0
+
+Current status:
+
+ ** lscomplex and lsreal both produce accurate results; they can be slow for
+ very large datasets.
+
+ ** fastlscomplex is accurate for the first octave of results; there is still an
+ error I need to pin down in the merging for additional octaves. fastlsreal
+ is disabled at the moment as I move to an implementation based on the new
+ fastlscomplex.
+
+ ** lscorrcoeff works, although I'm still attempting to understand the initial
+ 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.
+
Deleted: trunk/octave-forge/extra/lssa/SampleScriptWithVostokData
===================================================================
--- trunk/octave-forge/extra/lssa/SampleScriptWithVostokData 2012-07-31 14:43:11 UTC (rev 10794)
+++ trunk/octave-forge/extra/lssa/SampleScriptWithVostokData 2012-07-31 15:34:36 UTC (rev 10795)
@@ -1,226 +0,0 @@
-## Copyright (C) 2012 Benjamin Lewis
-## Licensed under the GNU GPL v2
-
-## This is just a sample script to introduce the purpose and usage of
-## the Lomb-Scargle Least Squares method with experimental data, here
-## using the Vostok ice core data collected and measured by J.R. Petit
-## et. al. and published in Nature; also available from the NOAA's
-## Paleoclimatology pages here:
-## <http://www.ncdc.noaa.gov/paleo/icecore/antarctica/vostok/vostok_data.html>.
-
-co2 = csvread("./data/co2.csv")(2:end,2:end);
-ch4 = csvread("./data/ch4.csv")(2:end,2:end);
-o18 = csvread("./data/o18.csv")(2:end,2:end);
-deut = csvread("./data/deut.csv")(2:end,2:end);
-dust = csvread("./data/dust.csv")(2:end,2:end);
-## The limited ranges are to deal with artifacts from the extraction of
-## the R data, notably that it leaves an extra column on the front of 0s
-## and the first row is text that Octave refuses to process.
-
-## Columns in co2 are Depth, Ice Age, Gas Age, and CO2 Concentration.
-## Columns in ch4 are Depth, Ice Age, Gas Age, and CH4 Concentration.
-## Columns in o18 are Depth, Ice Age, Gas Age, and Atmospheric O18.
-## Columns in dust are Depth, Ice Age, Dust Concentration.
-## Columns in deut are Depth, Ice Age, D concentration, and DeltaTS.
-
-co2_fig = figure("visible","off","name","CO2");
-ch4_fig = figure("visible","off","name","CH4");
-o18_fig = figure("visible","off","name","O18");
-deut_fig = figure("visible","off","name","Deuterium");
-dust_fig = figure("visible","off","name","Dust");
-## Generates figures and attaches handles to them for easy access.
-
-## Now we need some data to display; I'll run a few functions.
-ls_complex_co2_ice_age = lscomplex(co2(:,2),co2(:,4),1,100,20);
-ls_complex_co2_gas_age = lscomplex(co2(:,3),co2(:,4),1,100,20);
-ls_real_co2_ice_age = lsreal(co2(:,2),co2(:,4),1,100,20);
-ls_real_co2_gas_age = lsreal(co2(:,3),co2(:,4),1,100,20);
-ls_complex_ch4_ice_age = lscomplex(ch4(:,2),ch4(:,4),1,100,20);
-ls_complex_ch4_gas_age = lscomplex(ch4(:,3),ch4(:,4),1,100,20);
-ls_real_ch4_ice_age = lsreal(ch4(:,2),ch4(:,4),1,100,20);
-ls_real_ch4_gas_age = lsreal(ch4(:,3),ch4(:,4),1,100,20);
-ls_complex_o18_ice_age = lscomplex(o18(:,2),o18(:,4),1,100,20);
-ls_complex_o18_gas_age = lscomplex(o18(:,3),o18(:,4),1,100,20);
-ls_real_o18_ice_age = lsreal(o18(:,2),o18(:,4),1,100,20);
-ls_real_o18_gas_age = lsreal(o18(:,3),o18(:,4),1,100,20);
-ls_complex_deut = lscomplex(deut(:,2),deut(:,3),1,100,20);
-ls_real_deut = lsreal(deut(:,2),deut(:,3),1,100,20);
-ls_complex_dust = lscomplex(dust(:,2),dust(:,3),1,100,20);
-ls_real_dust = lsreal(dust(:,2),dust(:,3),1,100,20);
-
-x_data_axis_vector = [ -430000, 0 ];
-## Useful because all of the data extends over 430 000 years up to the
-## present.
-
-## Setting up the CO2 plots:
-figure(co2_fig);
-subplot(4,2,1);
-axis(x_data_axis_vector);
-plot(-(co2(:,2)),co2(:,4));
-title("Gas levels over ice age");
-subplot(4,2,2);
-axis(x_data_axis_vector);
-plot(-(co2(:,3),co2(:,4));
-title("Gas levels over gas age");
-subplot(4,2,3);
-plot(real(ls_complex_co2_ice_age));
-hold on;
-plot(imag(ls_complex_co2_ice_age),'r');
-title("Complex L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,4);
-plot(real(ls_complex_co2_gas_age));
-hold on;
-plot(imag(ls_complex_co2_gas_age),'r');
-title("Complex L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-subplot(4,2,5);
-plot(real(ls_real_co2_ice_age));
-hold on;
-plot(imag(ls_real_co2_ice_age),'r');
-title("Real L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,6);
-plot(real(ls_real_co2_gas_age));
-hold on;
-plot(imag(ls_real_co2_gas_age));
-title("Real L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-## At this point, we have transforms of both datasets, real and complex,
-## and just need to figure out what cool thing to do with the remaining slot.
-
-## Setting up the CH4 plots
-figure(ch4_fig);
-subplot(4,2,1);
-axis(x_data_axis_vector);
-plot(-(ch4(:,2)),ch4(:,4));
-title("Gas levels over ice age");
-subplot(4,2,2);
-axis(x_data_axis_vector);
-plot(-(ch4(:,3),ch4(:,4));
-title("Gas levels over gas age");
-subplot(4,2,3);
-plot(real(ls_complex_ch4_ice_age));
-hold on;
-plot(imag(ls_complex_ch4_ice_age),'r');
-title("Complex L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,4);
-plot(real(ls_complex_ch4_gas_age));
-hold on;
-plot(imag(ls_complex_ch4_gas_age),'r');
-title("Complex L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-subplot(4,2,5);
-plot(real(ls_real_ch4_ice_age));
-hold on;
-plot(imag(ls_real_ch4_ice_age),'r');
-title("Real L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,6);
-plot(real(ls_real_ch4_gas_age));
-hold on;
-plot(imag(ls_real_ch4_gas_age));
-title("Real L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-
-## Setting up the O18 plots:
-figure(o18_fig);
-subplot(4,2,1);
-axis(x_data_axis_vector);
-plot(-(o18(:,2)),o18(:,4));
-title("Gas levels over ice age");
-subplot(4,2,2);
-axis(x_data_axis_vector);
-plot(-(o18(:,3),o18(:,4));
-title("Gas levels over gas age");
-subplot(4,2,3);
-plot(real(ls_complex_o18_ice_age));
-hold on;
-plot(imag(ls_complex_o18_ice_age),'r');
-title("Complex L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,4);
-plot(real(ls_complex_o18_gas_age));
-hold on;
-plot(imag(ls_complex_o18_gas_age),'r');
-title("Complex L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-subplot(4,2,5);
-plot(real(ls_real_o18_ice_age));
-hold on;
-plot(imag(ls_real_o18_ice_age),'r');
-title("Real L-S transform of Gas/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,2,6);
-plot(real(ls_real_o18_gas_age));
-hold on;
-plot(imag(ls_real_o18_gas_age));
-title("Real L-S transform of Gas/gas age data");
-legend("Real part","Imaginary part");
-
-## Setting up Dust plots:
-figure(dust_fig);
-subplot(4,1,1);
-axis(x_data_axis_vector);
-plot(-(dust(:,2)),dust(:,3));
-title("Dust levels over ice age");
-subplot(4,1,2);
-plot(real(ls_complex_dust_ice_age));
-hold on;
-plot(imag(ls_complex_dust_ice_age),'r');
-title("Complex L-S transform of Dust/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,1,3);
-plot(real(ls_real_dust_ice_age));
-hold on;
-plot(imag(ls_real_dust_ice_age),'r');
-title("Real L-S transform of Dust/ice age data");
-legend("Real part","Imaginary part");
-
-## Setting up Deuterium plots:
-figure(deut_fig);
-subplot(4,1,1);
-axis(x_data_axis_vector);
-plot(-(deut(:,2)),deut(:,3));
-title("Deuterium levels over ice age");
-subplot(4,1,2);
-plot(real(ls_complex_deut_ice_age));
-hold on;
-plot(imag(ls_complex_deut_ice_age),'r');
-title("Complex L-S transform of Deuterium/ice age data");
-legend("Real part","Imaginary part");
-subplot(4,1,3);
-plot(real(ls_real_deut_ice_age));
-hold on;
-plot(imag(ls_real_deut_ice_age),'r');
-title("Real L-S transform of Deuterium/ice age data");
-legend("Real part","Imaginary part");
-
-co2_ch4_comparison_figure = figure("visible","off","name","CO2/CH4
-comparison");
-subplot(4,1,1);
-axes(x_data_axis_vector);
-plot(-(co2(:,2)),co2(:,4));
-hold on;
-plot(-(ch4(:,2)),ch4(:,4),'g');
-title("CO2 and CH4 data");
-legend("CO2","CH4");
-
-subplot(4,1,2);
-plot(abs(ls_complex_co2_ice_age));
-hold on;
-plot(abs(ls_complex_ch4_gas_age),'g');
-title("Abs. values of CO2 and CH4 L-S complex transforms");
-legend("CO2,CH4");
-
-
-
-
-
-
-## to implement:
-## - displays of all the data and flaws in trying to model with just
- ## using L-S data
-## - correlations of every data set with every other data set
-## - Comparing ls* results to periodogram results
\ No newline at end of file
Deleted: trunk/octave-forge/extra/lssa/inst/SampleDataSet.m
===================================================================
--- trunk/octave-forge/extra/lssa/inst/SampleDataSet.m 2012-07-31 14:43:11 UTC (rev 10794)
+++ trunk/octave-forge/extra/lssa/inst/SampleDataSet.m 2012-07-31 15:34:36 UTC (rev 10795)
@@ -1,25 +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 2 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/>.
-
-## No function structure to this, I just want to use it to store the
-## sums of sines and cosines I'll use for testing.
-
-xvec = linspace(0,8,1000);
-maxfreq = 4 / ( 2 * pi );
-
-yvec = ( 2.*sin(maxfreq.*xvec) + 3.*sin((3/4)*maxfreq.*xvec)
- - 0.5 .* sin((1/4)*maxfreq.*xvec) - 0.2 .* cos(maxfreq .* xvec)
- + cos((1/4)*maxfreq.*xvec));
-
Copied: trunk/octave-forge/extra/lssa/samples/SampleDataSet.m (from rev 10785, trunk/octave-forge/extra/lssa/inst/SampleDataSet.m)
===================================================================
--- trunk/octave-forge/extra/lssa/samples/SampleDataSet.m (rev 0)
+++ trunk/octave-forge/extra/lssa/samples/SampleDataSet.m 2012-07-31 15:34:36 UTC (rev 10795)
@@ -0,0 +1,25 @@
+## 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 2 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/>.
+
+## No function structure to this, I just want to use it to store the
+## sums of sines and cosines I'll use for testing.
+
+xvec = linspace(0,8,1000);
+maxfreq = 4 / ( 2 * pi );
+
+yvec = ( 2.*sin(maxfreq.*xvec) + 3.*sin((3/4)*maxfreq.*xvec)
+ - 0.5 .* sin((1/4)*maxfreq.*xvec) - 0.2 .* cos(maxfreq .* xvec)
+ + cos((1/4)*maxfreq.*xvec));
+
Copied: trunk/octave-forge/extra/lssa/samples/SampleScriptWithVostokData (from rev 10785, trunk/octave-forge/extra/lssa/SampleScriptWithVostokData)
===================================================================
--- trunk/octave-forge/extra/lssa/samples/SampleScriptWithVostokData (rev 0)
+++ trunk/octave-forge/extra/lssa/samples/SampleScriptWithVostokData 2012-07-31 15:34:36 UTC (rev 10795)
@@ -0,0 +1,226 @@
+## Copyright (C) 2012 Benjamin Lewis
+## Licensed under the GNU GPL v2
+
+## This is just a sample script to introduce the purpose and usage of
+## the Lomb-Scargle Least Squares method with experimental data, here
+## using the Vostok ice core data collected and measured by J.R. Petit
+## et. al. and published in Nature; also available from the NOAA's
+## Paleoclimatology pages here:
+## <http://www.ncdc.noaa.gov/paleo/icecore/antarctica/vostok/vostok_data.html>.
+
+co2 = csvread("./data/co2.csv")(2:end,2:end);
+ch4 = csvread("./data/ch4.csv")(2:end,2:end);
+o18 = csvread("./data/o18.csv")(2:end,2:end);
+deut = csvread("./data/deut.csv")(2:end,2:end);
+dust = csvread("./data/dust.csv")(2:end,2:end);
+## The limited ranges are to deal with artifacts from the extraction of
+## the R data, notably that it leaves an extra column on the front of 0s
+## and the first row is text that Octave refuses to process.
+
+## Columns in co2 are Depth, Ice Age, Gas Age, and CO2 Concentration.
+## Columns in ch4 are Depth, Ice Age, Gas Age, and CH4 Concentration.
+## Columns in o18 are Depth, Ice Age, Gas Age, and Atmospheric O18.
+## Columns in dust are Depth, Ice Age, Dust Concentration.
+## Columns in deut are Depth, Ice Age, D concentration, and DeltaTS.
+
+co2_fig = figure("visible","off","name","CO2");
+ch4_fig = figure("visible","off","name","CH4");
+o18_fig = figure("visible","off","name","O18");
+deut_fig = figure("visible","off","name","Deuterium");
+dust_fig = figure("visible","off","name","Dust");
+## Generates figures and attaches handles to them for easy access.
+
+## Now we need some data to display; I'll run a few functions.
+ls_complex_co2_ice_age = lscomplex(co2(:,2),co2(:,4),1,100,20);
+ls_complex_co2_gas_age = lscomplex(co2(:,3),co2(:,4),1,100,20);
+ls_real_co2_ice_age = lsreal(co2(:,2),co2(:,4),1,100,20);
+ls_real_co2_gas_age = lsreal(co2(:,3),co2(:,4),1,100,20);
+ls_complex_ch4_ice_age = lscomplex(ch4(:,2),ch4(:,4),1,100,20);
+ls_complex_ch4_gas_age = lscomplex(ch4(:,3),ch4(:,4),1,100,20);
+ls_real_ch4_ice_age = lsreal(ch4(:,2),ch4(:,4),1,100,20);
+ls_real_ch4_gas_age = lsreal(ch4(:,3),ch4(:,4),1,100,20);
+ls_complex_o18_ice_age = lscomplex(o18(:,2),o18(:,4),1,100,20);
+ls_complex_o18_gas_age = lscomplex(o18(:,3),o18(:,4),1,100,20);
+ls_real_o18_ice_age = lsreal(o18(:,2),o18(:,4),1,100,20);
+ls_real_o18_gas_age = lsreal(o18(:,3),o18(:,4),1,100,20);
+ls_complex_deut = lscomplex(deut(:,2),deut(:,3),1,100,20);
+ls_real_deut = lsreal(deut(:,2),deut(:,3),1,100,20);
+ls_complex_dust = lscomplex(dust(:,2),dust(:,3),1,100,20);
+ls_real_dust = lsreal(dust(:,2),dust(:,3),1,100,20);
+
+x_data_axis_vector = [ -430000, 0 ];
+## Useful because all of the data extends over 430 000 years up to the
+## present.
+
+## Setting up the CO2 plots:
+figure(co2_fig);
+subplot(4,2,1);
+axis(x_data_axis_vector);
+plot(-(co2(:,2)),co2(:,4));
+title("Gas levels over ice age");
+subplot(4,2,2);
+axis(x_data_axis_vector);
+plot(-(co2(:,3),co2(:,4));
+title("Gas levels over gas age");
+subplot(4,2,3);
+plot(real(ls_complex_co2_ice_age));
+hold on;
+plot(imag(ls_complex_co2_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_co2_gas_age));
+hold on;
+plot(imag(ls_complex_co2_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_co2_ice_age));
+hold on;
+plot(imag(ls_real_co2_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_co2_gas_age));
+hold on;
+plot(imag(ls_real_co2_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+## At this point, we have transforms of both datasets, real and complex,
+## and just need to figure out what cool thing to do with the remaining slot.
+
+## Setting up the CH4 plots
+figure(ch4_fig);
+subplot(4,2,1);
+axis(x_data_axis_vector);
+plot(-(ch4(:,2)),ch4(:,4));
+title("Gas levels over ice age");
+subplot(4,2,2);
+axis(x_data_axis_vector);
+plot(-(ch4(:,3),ch4(:,4));
+title("Gas levels over gas age");
+subplot(4,2,3);
+plot(real(ls_complex_ch4_ice_age));
+hold on;
+plot(imag(ls_complex_ch4_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_ch4_gas_age));
+hold on;
+plot(imag(ls_complex_ch4_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_ch4_ice_age));
+hold on;
+plot(imag(ls_real_ch4_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_ch4_gas_age));
+hold on;
+plot(imag(ls_real_ch4_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+
+## Setting up the O18 plots:
+figure(o18_fig);
+subplot(4,2,1);
+axis(x_data_axis_vector);
+plot(-(o18(:,2)),o18(:,4));
+title("Gas levels over ice age");
+subplot(4,2,2);
+axis(x_data_axis_vector);
+plot(-(o18(:,3),o18(:,4));
+title("Gas levels over gas age");
+subplot(4,2,3);
+plot(real(ls_complex_o18_ice_age));
+hold on;
+plot(imag(ls_complex_o18_ice_age),'r');
+title("Complex L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,4);
+plot(real(ls_complex_o18_gas_age));
+hold on;
+plot(imag(ls_complex_o18_gas_age),'r');
+title("Complex L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+subplot(4,2,5);
+plot(real(ls_real_o18_ice_age));
+hold on;
+plot(imag(ls_real_o18_ice_age),'r');
+title("Real L-S transform of Gas/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,2,6);
+plot(real(ls_real_o18_gas_age));
+hold on;
+plot(imag(ls_real_o18_gas_age));
+title("Real L-S transform of Gas/gas age data");
+legend("Real part","Imaginary part");
+
+## Setting up Dust plots:
+figure(dust_fig);
+subplot(4,1,1);
+axis(x_data_axis_vector);
+plot(-(dust(:,2)),dust(:,3));
+title("Dust levels over ice age");
+subplot(4,1,2);
+plot(real(ls_complex_dust_ice_age));
+hold on;
+plot(imag(ls_complex_dust_ice_age),'r');
+title("Complex L-S transform of Dust/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,1,3);
+plot(real(ls_real_dust_ice_age));
+hold on;
+plot(imag(ls_real_dust_ice_age),'r');
+title("Real L-S transform of Dust/ice age data");
+legend("Real part","Imaginary part");
+
+## Setting up Deuterium plots:
+figure(deut_fig);
+subplot(4,1,1);
+axis(x_data_axis_vector);
+plot(-(deut(:,2)),deut(:,3));
+title("Deuterium levels over ice age");
+subplot(4,1,2);
+plot(real(ls_complex_deut_ice_age));
+hold on;
+plot(imag(ls_complex_deut_ice_age),'r');
+title("Complex L-S transform of Deuterium/ice age data");
+legend("Real part","Imaginary part");
+subplot(4,1,3);
+plot(real(ls_real_deut_ice_age));
+hold on;
+plot(imag(ls_real_deut_ice_age),'r');
+title("Real L-S transform of Deuterium/ice age data");
+legend("Real part","Imaginary part");
+
+co2_ch4_comparison_figure = figure("visible","off","name","CO2/CH4
+comparison");
+subplot(4,1,1);
+axes(x_data_axis_vector);
+plot(-(co2(:,2)),co2(:,4));
+hold on;
+plot(-(ch4(:,2)),ch4(:,4),'g');
+title("CO2 and CH4 data");
+legend("CO2","CH4");
+
+subplot(4,1,2);
+plot(abs(ls_complex_co2_ice_age));
+hold on;
+plot(abs(ls_complex_ch4_gas_age),'g');
+title("Abs. values of CO2 and CH4 L-S complex transforms");
+legend("CO2,CH4");
+
+
+
+
+
+
+## to implement:
+## - displays of all the data and flaws in trying to model with just
+ ## using L-S data
+## - correlations of every data set with every other data set
+## - Comparing ls* results to periodogram results
\ No newline at end of file
Modified: trunk/octave-forge/extra/lssa/src/Makefile
===================================================================
--- trunk/octave-forge/extra/lssa/src/Makefile 2012-07-31 14:43:11 UTC (rev 10794)
+++ trunk/octave-forge/extra/lssa/src/Makefile 2012-07-31 15:34:36 UTC (rev 10795)
@@ -1,13 +1,14 @@
MKOCTFILE ?= mkoctfile
-all: fastlscomplex.oct \
- fastlsreal.oct
+all: fastlscomplex.oct #\
+# fastlsreal.oct
fastlscomplex.oct: fastlscomplex.cc
$(MKOCTFILE) fastlscomplex.cc
-fastlsreal.oct: fastlsreal.cc
- $(MKOCTFILE) fastlsreal.cc
+# fastlsreal compilation is disabled for the time being
+#fastlsreal.oct: fastlsreal.cc
+# $(MKOCTFILE) fastlsreal.cc
# helper function just in case
clean:
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 14:43:21
|
Revision: 10794
http://octave.svn.sourceforge.net/octave/?rev=10794&view=rev
Author: abarth93
Date: 2012-07-31 14:43:11 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/nccoord.m
Modified: trunk/octave-forge/extra/ncArray/inst/nccoord.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-31 14:33:36 UTC (rev 10793)
+++ trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-31 14:43:11 UTC (rev 10794)
@@ -38,7 +38,7 @@
% check for coordinate dimensions
for i=1:length(dims)
% check if variable with the same name than the dimension exist
- index = find(strcmp(dims{i},{finfo.Variables(:).Name}));
+ index = find(strcmp(dims{i},{finfo.Variables(:).Name}),1);
if ~isempty(index)
coord = addcoord(coord,dims{i},finfo);
end
@@ -50,7 +50,7 @@
function coord = addcoord(coord,name,finfo)
% check if coordinate is aleady in the list
-if isempty(find(strcmp(name,{coord(:).name})))
+if isempty(find(strcmp(name,{coord(:).name}),1))
% check if name is variable
index = find(strcmp(name,{finfo.Variables(:).Name}));
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 14:33:47
|
Revision: 10793
http://octave.svn.sourceforge.net/octave/?rev=10793&view=rev
Author: abarth93
Date: 2012-07-31 14:33:36 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/main/octcdf/inst/ncinfo.m
trunk/octave-forge/main/octcdf/inst/ncread.m
Modified: trunk/octave-forge/main/octcdf/inst/ncinfo.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 14:32:53 UTC (rev 10792)
+++ trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 14:33:36 UTC (rev 10793)
@@ -41,7 +41,12 @@
tmp.Length = dims{i}(:);
% requires octcdf 1.1.6
%tmp.Unlimited = isrecdim(dims{i});
- vinfo.Dimensions = [vinfo.Dimensions; tmp];
+
+ if isempty(vinfo.Dimensions)
+ vinfo.Dimensions = [tmp];
+ else
+ vinfo.Dimensions(i) = tmp;
+ end
end
@@ -55,7 +60,12 @@
tmp.Name = nm;
tmp.Value = na{j}(:);
- vinfo.Attributes = [vinfo.Attributes; tmp];
+
+ if isempty(vinfo.Attributes)
+ vinfo.Attributes = [tmp];
+ else
+ vinfo.Attributes(j) = tmp;
+ end
end
vinfo.FillValue = fillval(nv);
Modified: trunk/octave-forge/main/octcdf/inst/ncread.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncread.m 2012-07-31 14:32:53 UTC (rev 10792)
+++ trunk/octave-forge/main/octcdf/inst/ncread.m 2012-07-31 14:33:36 UTC (rev 10793)
@@ -77,6 +77,11 @@
end
x = permute(x,[ndims(x):-1:1]);
+
+if length(count) < 2
+ count(2) = 1;
+end
+
x = reshape(x,count);
close(nc)
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 14:33:04
|
Revision: 10792
http://octave.svn.sourceforge.net/octave/?rev=10792&view=rev
Author: abarth93
Date: 2012-07-31 14:32:53 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/test_ncarray.m
Modified: trunk/octave-forge/extra/ncArray/inst/test_ncarray.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/test_ncarray.m 2012-07-31 13:45:19 UTC (rev 10791)
+++ trunk/octave-forge/extra/ncArray/inst/test_ncarray.m 2012-07-31 14:32:53 UTC (rev 10792)
@@ -314,7 +314,7 @@
assert(isequalwithequalnans(x,lon_ref(1:3:end,:)))
- assert(strcmp(SST.units,'degC'))
+ %assert(strcmp(SST.units,'degC'))
assert(strcmp(SST.('units'),'degC'))
end
@@ -335,7 +335,7 @@
SST_ref = ncread(files{2},'SST');
assert(isequalwithequalnans(SST_test,SST_ref))
-assert(strcmp(CA2.units,'degC'));
+assert(strcmp(CA2.('units'),'degC'));
disp('All tests passed.')
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 13:45:25
|
Revision: 10791
http://octave.svn.sourceforge.net/octave/?rev=10791&view=rev
Author: abarth93
Date: 2012-07-31 13:45:19 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/main/octcdf/inst/@ncdim/isrecord.m
trunk/octave-forge/main/octcdf/inst/ncinfo.m
Modified: trunk/octave-forge/main/octcdf/inst/@ncdim/isrecord.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/@ncdim/isrecord.m 2012-07-31 11:44:13 UTC (rev 10790)
+++ trunk/octave-forge/main/octcdf/inst/@ncdim/isrecord.m 2012-07-31 13:45:19 UTC (rev 10791)
@@ -1,3 +1,4 @@
+% depreciated: use isrecdim instead
function isr = isrecord(self)
isr = ncisrecord(self);
end
Modified: trunk/octave-forge/main/octcdf/inst/ncinfo.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 11:44:13 UTC (rev 10790)
+++ trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 13:45:19 UTC (rev 10791)
@@ -30,7 +30,7 @@
vinfo.Size = fliplr(size(nv));
vinfo.Filename = filename;
-vinfo.Dimensions = {};
+vinfo.Dimensions = [];
vinfo.Name = varname;
dims = fliplr(dim(nv));
@@ -39,8 +39,9 @@
tmp = struct();
tmp.Name = name(dims{i});
tmp.Length = dims{i}(:);
- tmp.Unlimited = isrecord(dims{i});
- vinfo.Dimensions(i) = tmp;
+ % requires octcdf 1.1.6
+ %tmp.Unlimited = isrecdim(dims{i});
+ vinfo.Dimensions = [vinfo.Dimensions; tmp];
end
@@ -54,7 +55,7 @@
tmp.Name = nm;
tmp.Value = na{j}(:);
- vinfo.Attributes(j) = tmp;
+ vinfo.Attributes = [vinfo.Attributes; tmp];
end
vinfo.FillValue = fillval(nv);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 11:44:19
|
Revision: 10790
http://octave.svn.sourceforge.net/octave/?rev=10790&view=rev
Author: abarth93
Date: 2012-07-31 11:44:13 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
matlab compat
Modified Paths:
--------------
trunk/octave-forge/main/octcdf/inst/ncinfo.m
Modified: trunk/octave-forge/main/octcdf/inst/ncinfo.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 11:34:27 UTC (rev 10789)
+++ trunk/octave-forge/main/octcdf/inst/ncinfo.m 2012-07-31 11:44:13 UTC (rev 10790)
@@ -5,6 +5,7 @@
% return information about complete NetCDF file (filename) or about
% the specific variable varname.
+%
function info = ncinfo(filename,varname)
nc = netcdf(filename,'r');
@@ -35,13 +36,17 @@
dims = fliplr(dim(nv));
for i=1:length(dims)
- vinfo.Dimensions{i} = name(dims{i});
+ tmp = struct();
+ tmp.Name = name(dims{i});
+ tmp.Length = dims{i}(:);
+ tmp.Unlimited = isrecord(dims{i});
+ vinfo.Dimensions(i) = tmp;
end
na = att(nv);
-%vinfo.Attributes = [];
+vinfo.Attributes = [];
for j=1:length(na)
tmp = struct();
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 11:34:34
|
Revision: 10789
http://octave.svn.sourceforge.net/octave/?rev=10789&view=rev
Author: abarth93
Date: 2012-07-31 11:34:27 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
bug fix: elements of struct arrays are found by strcmp instead of strmatch
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/@ncBaseArray/subsref.m
trunk/octave-forge/extra/ncArray/inst/ncCatArray.m
trunk/octave-forge/extra/ncArray/inst/nccoord.m
Modified: trunk/octave-forge/extra/ncArray/inst/@ncBaseArray/subsref.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/@ncBaseArray/subsref.m 2012-07-31 08:46:57 UTC (rev 10788)
+++ trunk/octave-forge/extra/ncArray/inst/@ncBaseArray/subsref.m 2012-07-31 11:34:27 UTC (rev 10789)
@@ -52,7 +52,7 @@
elseif strcmp(idx.type,'.')
% load attribute
name = idx.subs;
- index = strmatch(name,{self.vinfo.Attributes(:).Name});
+ index = find(strcmp(name,{self.vinfo.Attributes(:).Name}));
if isempty(index)
error('variable %s has no attribute called %s',self.varname,name);
Modified: trunk/octave-forge/extra/ncArray/inst/ncCatArray.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/ncCatArray.m 2012-07-31 08:46:57 UTC (rev 10788)
+++ trunk/octave-forge/extra/ncArray/inst/ncCatArray.m 2012-07-31 11:34:27 UTC (rev 10789)
@@ -77,7 +77,7 @@
for i=1:length(coord)
- % coordinates do also depend on the dimension only which we concatenate
+ % coordinates do also depend on the dimension on which we concatenate
coord(i).val = arr(dim,filenames,coord(i).name);
if dim > length(coord(i).dims)
coord(i).dims{dim} = catdimname;
Modified: trunk/octave-forge/extra/ncArray/inst/nccoord.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-31 08:46:57 UTC (rev 10788)
+++ trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-31 11:34:27 UTC (rev 10789)
@@ -24,19 +24,21 @@
coord = struct('name',{},'dims',{});
% check the coordinate attribute
-index = strmatch('coordinates',{vinfo.Attributes(:).Name});
-if ~isempty(index)
+if ~isempty(vinfo.Attributes)
+ index = find(strcmp('coordinates',{vinfo.Attributes(:).Name}));
+ if ~isempty(index)
tmp = strsplit(vinfo.Attributes(index).Value,' ');
for i=1:length(tmp)
coord = addcoord(coord,tmp{i},finfo);
end
+ end
end
% check for coordinate dimensions
for i=1:length(dims)
% check if variable with the same name than the dimension exist
- index = strmatch(dims{i},{finfo.Variables(:).Name});
+ index = find(strcmp(dims{i},{finfo.Variables(:).Name}));
if ~isempty(index)
coord = addcoord(coord,dims{i},finfo);
end
@@ -48,14 +50,17 @@
function coord = addcoord(coord,name,finfo)
% check if coordinate is aleady in the list
-if isempty(strmatch(name,{coord(:).name}))
+if isempty(find(strcmp(name,{coord(:).name})))
% check if name is variable
- index = strmatch(name,{finfo.Variables(:).Name});
+ index = find(strcmp(name,{finfo.Variables(:).Name}));
if ~isempty(index)
c.name = name;
- c.dims = {finfo.Variables(index).Dimensions(:).Name};
-
+ d = finfo.Variables(index).Dimensions;
+ c.dims = {d(:).Name};
+ % does not work in octave
+ %c.dims = {finfo.Variables(index).Dimensions(:).Name};
+ %keyboard
coord(end+1) = c;
end
end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 08:47:08
|
Revision: 10788
http://octave.svn.sourceforge.net/octave/?rev=10788&view=rev
Author: abarth93
Date: 2012-07-31 08:46:57 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
bug fix for compressed files
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
Modified: trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/cached_decompress.m 2012-07-31 08:46:24 UTC (rev 10787)
+++ trunk/octave-forge/extra/ncArray/inst/cached_decompress.m 2012-07-31 08:46:57 UTC (rev 10788)
@@ -19,7 +19,6 @@
cache_dir = CACHED_DECOMPRESS_DIR;
if isempty(cache_dir)
-% cache_dir = fullfile(getenv('HOME'),'tmp','Cache');
cache_dir = fullfile(getenv('HOME'),'tmp','Cache');
end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 08:46:30
|
Revision: 10787
http://octave.svn.sourceforge.net/octave/?rev=10787&view=rev
Author: abarth93
Date: 2012-07-31 08:46:24 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
bug fix for compressed files
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
Modified: trunk/octave-forge/extra/ncArray/inst/cached_decompress.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/cached_decompress.m 2012-07-31 08:29:25 UTC (rev 10786)
+++ trunk/octave-forge/extra/ncArray/inst/cached_decompress.m 2012-07-31 08:46:24 UTC (rev 10787)
@@ -10,7 +10,7 @@
% Alexander Barth, 2012-06-13
%
-function [fname]=cached_decompress(url)
+function [fname] = cached_decompress(url)
global CACHED_DECOMPRESS_DIR
@@ -23,7 +23,11 @@
cache_dir = fullfile(getenv('HOME'),'tmp','Cache');
end
-if beginswith(url,'http:') || ~endswith(url,'.gz') || ~endswith(url,'.bz2')
+% do nothing if
+% file is a a remote url (begins with http:)
+% or file does not end with .gz or .bz2
+
+if beginswith(url,'http:') || ~(endswith(url,'.gz') || endswith(url,'.bz2'))
% opendap url or not compressed file
fname = url;
return
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <aba...@us...> - 2012-07-31 08:29:35
|
Revision: 10786
http://octave.svn.sourceforge.net/octave/?rev=10786&view=rev
Author: abarth93
Date: 2012-07-31 08:29:25 +0000 (Tue, 31 Jul 2012)
Log Message:
-----------
matlab compat
Modified Paths:
--------------
trunk/octave-forge/extra/ncArray/inst/nccoord.m
Modified: trunk/octave-forge/extra/ncArray/inst/nccoord.m
===================================================================
--- trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-30 18:00:45 UTC (rev 10785)
+++ trunk/octave-forge/extra/ncArray/inst/nccoord.m 2012-07-31 08:29:25 UTC (rev 10786)
@@ -18,7 +18,7 @@
% determine coordinates
% using CF convention
-dims = vinfo.Dimensions;
+dims = {vinfo.Dimensions(:).Name};
% create empty coord array with the fields name and dims
coord = struct('name',{},'dims',{});
@@ -54,7 +54,7 @@
index = strmatch(name,{finfo.Variables(:).Name});
if ~isempty(index)
c.name = name;
- c.dims = finfo.Variables(index).Dimensions;
+ c.dims = {finfo.Variables(index).Dimensions(:).Name};
coord(end+1) = c;
end
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-07-30 18:00:55
|
Revision: 10785
http://octave.svn.sourceforge.net/octave/?rev=10785&view=rev
Author: benjf5
Date: 2012-07-30 18:00:45 +0000 (Mon, 30 Jul 2012)
Log Message:
-----------
Some modifications made, still getting everything ready.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/src/fastlsreal.cc
Added Paths:
-----------
trunk/octave-forge/extra/lssa/src/Makefile
Added: trunk/octave-forge/extra/lssa/src/Makefile
===================================================================
--- trunk/octave-forge/extra/lssa/src/Makefile (rev 0)
+++ trunk/octave-forge/extra/lssa/src/Makefile 2012-07-30 18:00:45 UTC (rev 10785)
@@ -0,0 +1,14 @@
+MKOCTFILE ?= mkoctfile
+
+all: fastlscomplex.oct \
+ fastlsreal.oct
+
+fastlscomplex.oct: fastlscomplex.cc
+ $(MKOCTFILE) fastlscomplex.cc
+
+fastlsreal.oct: fastlsreal.cc
+ $(MKOCTFILE) fastlsreal.cc
+
+# helper function just in case
+clean:
+ rm *.o *.oct *~ octave-core
Modified: trunk/octave-forge/extra/lssa/src/fastlsreal.cc
===================================================================
--- trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-07-29 21:51:22 UTC (rev 10784)
+++ trunk/octave-forge/extra/lssa/src/fastlsreal.cc 2012-07-30 18:00:45 UTC (rev 10785)
@@ -55,11 +55,14 @@
}
-ComplexRowVector flsreal( RowVector tvec , ComplexRowVector xvec ,
+ComplexRowVector flsreal( RowVector tvec , RowVector xvec ,
double maxfreq, int octaves, int coefficients ) {
+ struct XTElem {
+ double x, t;
+ };
struct Precomputation_Record {
Precomputation_Record *next;
- std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only.
+ XTElem power_series[12]; // I'm using 12 as a matter of compatibility, only.
bool stored_data;
};
@@ -67,16 +70,15 @@
double tau, delta_tau, tau_0, tau_h, n_inv, mu,
omega_oct, omega_multiplier, octavemax, omega_working,
- loop_tau_0, loop_delta_tau;
+ loop_tau_0, loop_delta_tau, x;
double length = ( tvec((tvec.numel()-1)) - tvec( octave_idx_type (0)));
int octave_iter, coeff_iter;
std::complex<double> zeta, z_accumulator, zeta_exp_term, zeta_exp_multiplier, alpha,
- iota, i_accumulator, iota_exp_term, iota_exp_multiplier;
+ iota, i_accumulator, iota_exp_term, iota_exp_multiplier, exp_squared, exp_squared_multiplier;
octave_idx_type n = tvec.numel();
- std::complex<double> temp_array[12];
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- temp_array[array_iter] = std::complex<double> ( 0 , 0 );
- }
+ XTElem *tpra, *temp_ptr_alpha, temp_alpha[12], *tprb, *temp_ptr_beta, temp_beta[12], temp_array[12];
+
+
int factorial_array[12];
factorial_array[0] = 1;
for ( int i = 1 ; i < 12 ; i++ ) {
@@ -95,104 +97,127 @@
const std::complex<double> im = std::complex<double> ( 0 , 1 ); //I seriously prefer C99's complex numbers.
octave_idx_type k ( 0 ); // Iterator for accessing xvec, tvec.
-
- Precomputation_Record * complex_precomp_records_head, *complex_record_current,
- *complex_record_tail, *complex_record_ref, *complex_record_next, *iota_precomp_records_head,
- *iota_record_current, *iota_record_tail, *iota_record_ref, *iota_record_next;
- complex_record_current = complex_precomp_records_head = new Precomputation_Record;
- iota_record_current = iota_precomp_records_head = new Precomputation_Record;
- for ( size_t p_r_iter = 1 ; p_r_iter < precomp_subset_count ; p_r_iter++ ) {
- complex_record_current->next = new Precomputation_Record;
- iota_record_current->next = new Precomputation_Record;
- complex_record_current = complex_record_current->next;
- iota_record_current = iota_record_current->next;
- }
- // Error's past this point
- complex_record_tail = complex_record_current;
- iota_record_tail = iota_record_current;
- complex_record_current = complex_precomp_records_head;
- iota_record_current = iota_precomp_records_head;
- complex_record_tail->next = 0;
- iota_record_tail->next = 0;
- /* A test needs to be included for if there was a failure, but since
- * precomp_subset_count is of type size_t, it should be okay. */
- for( ; complex_record_current != 0 ; complex_record_current = complex_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- complex_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- } // To avoid any trouble down the line, although it is an annoyance.
- while ( (k < n) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double p;
- for ( int j = 0 ; j < 12 ; j++ ) {
- alpha.real() = xvec(k).real();
- alpha.imag() = xvec(k).imag();
- // Change to switches for easier manipulation, fewer tests. This is two tests where one will do.
- if ( !( j % 2 ) ) {
- if ( ! ( j % 4 ) ) {
- alpha.real() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- } else {
- if ( ! ( j % 3 ) ) {
- alpha.real() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- }
- complex_record_current->power_series[j].real() += alpha.real();
- complex_record_current->power_series[j].imag() += alpha.imag();
- }
- // Computes each next step of the power series for the given power series element.
- // j was reused since it was a handy inner-loop variable, even though I used it twice here.
- complex_record_current->stored_data = true;
- k++;
+
+ Precomputation_Record * precomp_records_head, *record_current, *record_tail, *record_ref, *record_next;
+ record_current = precomp_records_head = new Precomputation_Record;
+ for ( te = tvec(k) + (2 * delta_tau) ; ; ) {
+ x = xvec(k);
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ p->x = x;
+ (p++)->t = 1;
+ // p = 1
+ tt = -t;
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 2
+ tt *= t*(1.0/2.0);
+ p->x = x*tt;
+ (p++)->t = tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 4
+ tt *= t*(1.0/4.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 6
+ tt *= t*(1.0/6.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 8
+ tt *= t*(1.0/8.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 10
+ tt *= t*(1.0/10.0);
+ p->x = x * tt;
+ (p++)->t = tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ p->x = x * tt;
+ (p++)->t = tt;
}
- tau_h += ( 2 * delta_tau );
- }
- // At this point all precomputation records have been
- // exhausted for complex records. Short-circuit is abused
- // to avoid overflow errors.
- // Reset k, tau_h to reset the process. I may rewrite
- // these loops to be one, since running twice as long to
- // do the same thing is painful. May also move to a switch
- // in the previous section too.
- k = 0;
- tau_h = tau_0;
- for( ; iota_record_current != 0 ; iota_record_current = iota_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- iota_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- }
- while( ( k < n ) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double comps[12];
- iota_record_current->power_series[0].real() = 1;
- comps[0] = 1;
- for ( int j = 1 ; j < 12 ; j++ ) {
- comps[j] = comps[j-1] * mu * (tvec(k)-tau_h);
- switch ( j % 4 ) {
- case 0 :
- iota_record_current->power_series[j].real() += comps[j] / factorial_array[j] ;
- break;
- case 1:
- iota_record_current->power_series[j].imag() += comps[j] / factorial_array[j] ;
- break;
- case 2:
- iota_record_current->power_series[j].real() -= comps[j] / factorial_array[j] ;
- break;
- case 3:
- iota_record_current->power_series[j].imag() -= comps[j] / factorial_array[j] ;
- break;
- }
+ record_current->stored_data = true;
+ for(k++; ( k < n ) && tvec(k) < te ; k++ ) {
+ x = xvec(k);
+ {
+ double t = mu*(tvec(k)-tau_h), tt;
+ p = record_current->power_series;
+ // p = 0
+ p->x += x;
+ (p++)->t += 1;
+ // p = 1
+ tt = -t;
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 2
+ tt *= t*(1.0/2.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 3
+ tt *= t*(-1.0/3.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 4
+ tt *= t*(1.0/4.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 5
+ tt *= t*(-1.0/5.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 6
+ tt *= t*(1.0/6.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 7
+ tt *= t*(-1.0/7.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 8
+ tt *= t*(1.0/8.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 9
+ tt *= t*(-1.0/9.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 10
+ tt *= t*(1.0/10.0);
+ p->x += x * tt;
+ (p++)->t += tt;
+ // p = 11
+ tt *= t*(-1.0/11.0);
+ p->x += x * tt;
+ (p++)->t += tt;
}
- iota_record_current->stored_data = true;
- k++;
+ record_current->stored_data = true;
}
- tau_h += ( 2 * delta_tau );
+ if( k >= n ) break;
+ tau_h = te + delta_tau;
+ te = tau_h + delta_tau;
+ record_current->next = new Precomputation_Record;
+ record_current = record_current->next;
}
-
+ record_tail = record_current;
+ record_current = precomp_records_head;
+ record_tail->next = 0;
/* Summation of coefficients for each frequency. As we have ncoeffs * noctaves elements,
* it makes sense to work from the top down, as we have omega_max by default (maxfreq)
@@ -206,142 +231,155 @@
octave_idx_type iter ( 0 );
- double zeta_real_part = 0, zeta_imag_part = 0, zeta_real_part_accumulator = 0, zeta_imag_part_accumulator = 0,
- iota_real_part = 0, iota_imag_part = 0, iota_real_part_accumulator = 0, iota_imag_part_accumulator = 0;
-
// Loops need to first travel over octaves, then coefficients;
for ( octave_iter = octaves ; ; omega_oct *= 0.5 , octavemax *= 0.5 , loop_tau_0 += loop_delta_tau , loop_delta_tau *= 2 ) {
- omega_working = omega_oct;
- zeta_exp_term = std::complex<double> ( cos ( - omega_working * loop_tau_0 ) ,
- sin ( - omega_working * loop_tau_0 ) );
- zeta_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- iota_exp_term = std::complex<double> ( cos ( - 2 * omega_working * loop_tau_0 ) ,
- sin ( - 2 * omega_working * loop_tau_0 ) );
- iota_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, omega_working *= omega_multiplier){
- zeta_real_part_accumulator = 0;
- zeta_imag_part_accumulator = 0;
- zeta_real_part = 0;
- zeta_imag_part = 0;
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next, zeta_exp_term *= zeta_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- z_accumulator = ( pow(omega_working,array_iter) * complex_record_current->power_series[array_iter] );
- zeta_real_part_accumulator += z_accumulator.real();
- zeta_imag_part_accumulator += z_accumulator.imag();
+ o = omega_oct;
+ omega_working = octavemax;
+ for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, o *= omega_multiplier, omega_working *= omega_multiplier){
+ exp_term = std::complex<double> ( cos( - omega_working * loop_tau_0 ) ,
+ sin ( - omega_working * loop_tau_0 ) );
+ exp_squared = exp_term * exp_term;
+ exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
+ sin ( - 2 * omega_working * loop_delta_tau ) );
+ exp_squared_multiplier = exp_multiplier * exp_multiplier;
+ for ( zeta = iota = 0, record_current = precomp_records_head ; record_current ;
+ record_current = record_current->next, exp_term *= exp_multiplier,
+ exp_squared *= exp_squared_multiplier ) {
+ if ( record_current->stored_data ) {
+ int p;
+ for ( zz = ii = 0 , p = 0, on_1 = n_1 ; p < 12 ; ) {
+ zz.real() += record_current->power_series[p]->x * on_1;
+ ii.real() += record_current->power_series[p++]-> t * o2n_1;
+ on_1 *= o;
+ o2n_1 *= o2;
+ zz.imag() += record_current->power_series[p]->x * on_1;
+ ii.imag() += record_current->power_series[p++]-> t * o2n_1;
+ on_1 *= o;
+ o2n_1 *= o2;
+ }
+ zeta += exp_term * zz;
+ iota += exp_squared * ii;
}
- zeta_real_part = zeta_real_part + ( zeta_exp_term.real() * zeta_real_part_accumulator - zeta_exp_term.imag() * zeta_imag_part_accumulator );
- zeta_imag_part = zeta_imag_part + ( zeta_exp_term.imag() * zeta_real_part_accumulator + zeta_exp_term.real() * zeta_imag_part_accumulator );
}
- for ( iota_record_current = iota_precomp_records_head; iota_record_current ;
- iota_record_current = iota_record_current->next, iota_exp_term *= iota_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- i_accumulator = ( pow(omega_working,array_iter) * iota_record_current->power_series[array_iter] );
- iota_real_part_accumulator += i_accumulator.real();
- iota_imag_part_accumulator += i_accumulator.imag();
- }
- iota_real_part = iota_real_part + ( iota_exp_term.real() * iota_real_part_accumulator - iota_exp_term.imag() * iota_imag_part_accumulator );
- iota_imag_part = iota_imag_part + ( iota_exp_term.imag() * iota_real_part_accumulator + iota_exp_term.real() * iota_imag_part_accumulator );
- }
- // c + i s = 2 * ( zeta-omega-conjugate - iota-2omega-conjuage * zeta-omega ) / ( 1 - abs(iota-2omega) ^ 2 )
- // (is what the results will be.)
- zeta_real_part *= n_inv;
- zeta_imag_part *= n_inv;
- iota_real_part *= n_inv;
- iota_imag_part *= n_inv;
- double result_real_part = 2 * ( zeta_real_part - iota_real_part * zeta_real_part - iota_imag_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- double result_imag_part = 2 * ( - zeta_imag_part + iota_imag_part * zeta_real_part - iota_real_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- results(iter) = std::complex<double> ( result_real_part , result_imag_part );
+ results(iter) = 2 / ( 1 - ( iota.real() * iota.real() ) - (iota.imag() *
+ iota.imag() )
+ * ( conj(zeta) - conj(iota) * zeta );
iter++;
}
if ( !(--octave_iter) ) break;
/* If we've already reached the lowest value, stop.
* Otherwise, merge with the next computation range.
*/
- double exp_power_series_elements[12];
+ double *exp_pse_ptr, *exp_ptr, exp_power_series_elements[12];
exp_power_series_elements[0] = 1;
+ exp_pse_ptr = exp_ptr = exp_power_series_elements;
for ( int r_iter = 1 ; r_iter < 12 ; r_iter++ ) {
exp_power_series_elements[r_iter] = exp_power_series_elements[r_iter-1]
* ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) );
}
try{
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next ) {
- if ( ! ( complex_record_ref = complex_record_current->next ) || ! complex_record_ref->stored_data ) {
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q )];
+ for ( record_current = precomp_records_head ; record_current ;
+ record_current = record_current->next ) {
+ if ( ! ( record_ref = record_current->next ) || ! record_ref->stored_data ) {
+ // In this case, there is no next record, but this record has data.
+ if ( record_current->stored_data ) {
+ int p = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1 , temp_ptr_alpha = temp_alpha ; p < 12 ; p++ , exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ temp_ptr_alpha->x = record_current->power_series[p]->x;
+ (temp_ptr_alpha++)->t = record_current->power_series[p]->t;
+ temp_ptr_beta->x = -record_current->power_series[p]->x;
+ (temp_ptr_beta++)->t = -record_current->power_series[p]->t;
+ for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *temp_ptr_alpha * *exp_ptr; ; ) {
+ /* This next block is from Mathias' code, and it does a few
+ * ... unsavoury things. First off, it uses conditionals with
+ * break in order to avoid potentially accessing null regions
+ * of memory, and then it does ... painful things with a few
+ * numbers. However, remembering that most of these will not
+ * actually be accessed for the first iterations, it's easier.
+ */
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = *tpra * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_current->power_series[p] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = -*tpra * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_current->power_series[p] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
- if ( ! complex_record_ref ) break; // Last record was reached
+ if ( ! record_ref ) break; // Last record was reached
}
else {
- complex_record_next = complex_record_ref;
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q )] - complex_record_next->power_series[p - (2*q)] );
+ record_next = record_ref;
+ if ( record_current->stored_data ) {
+ int p = 0, q = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; p < 12 ; p++, q++, exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ *temp_ptr_alpha++ = record_current->power_series[p] + record_next->power_series[q];
+ *temp_ptr_beta++ = record_current->power_series[p] - record_next->power_series[1];
+ tprb = temp_ptr_beta;
+ for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *tpra * *exp_ptr; ; ) {
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tprb -= 2;
+ h = *tprb * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tpra -= 2;
+ record_current->power_series[p] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tprb -= 2;
+ h = - *tprb * *exp_ptr;
+ record_current->power_series[p].real() -= h.imag();
+ record_current->power_series[p].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ tpra -= 2;
+ record_current->power_series[p] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q ) - 1] - complex_record_next->power_series[p - ( 2 * q ) - 1 ] );
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
} else {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_next->power_series[p - ( 2 * q )];
+ int q = 0;
+ for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; q < 12; q++, exp_pse_ptr++ ) {
+ tpra = temp_ptr_alpha;
+ *temp_ptr_alpha++ = std::complex<double>(record_next->power_series[q]);
+ for ( exp_ptr = exp_power_series_elements, record_next->power_series[q] = *tpra * *exp_ptr; ; ) {
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = *tpra * *exp_ptr;
+ record_next->power_series[q].real() -= h.imag();
+ record_next->power_series[q].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_next->power_series[q] -= *tpra * *exp_ptr;
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ h = -*tpra * *exp_ptr;
+ record_next->power_series[q].real() -= h.imag();
+ record_next->power_series[q].imag() += h.real();
+ if ( ++exp_ptr >= exp_pse_ptr ) break;
+ --tpra;
+ record_next->power_series[q] += *tpra * *exp_ptr;
}
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,(q+1)) * complex_record_next->power_series[p - ( 2 * q ) - 1];
- }
}
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
}
- complex_record_current->stored_data = true;
- complex_record_ref = complex_record_next;
- complex_record_current->next = complex_record_ref->next;
- delete complex_record_ref;
+ record_current->stored_data = true;
+ record_ref = record_next;
+ record_current->next = record_ref->next;
+ delete record_ref;
}
}
}
- } catch (std::exception& e) { //This section was part of my debugging, and may be removed.
- std::cout << "Exception thrown: " << e.what() << std::endl;
- ComplexRowVector exception_result (1);
- exception_result(0) = std::complex<double> ( 0,0);
- return exception_result;
- }
- }
+
+
return results;
}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-07-29 21:51:28
|
Revision: 10784
http://octave.svn.sourceforge.net/octave/?rev=10784&view=rev
Author: paramaniac
Date: 2012-07-29 21:51:22 +0000 (Sun, 29 Jul 2012)
Log Message:
-----------
control: fix copy-paste error in error message
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/lsim.m
Modified: trunk/octave-forge/main/control/inst/lsim.m
===================================================================
--- trunk/octave-forge/main/control/inst/lsim.m 2012-07-29 13:17:39 UTC (rev 10783)
+++ trunk/octave-forge/main/control/inst/lsim.m 2012-07-29 21:51:22 UTC (rev 10784)
@@ -144,7 +144,7 @@
if (isempty (x0))
x0 = zeros (n, 1);
elseif (n != length (x0) || ! is_real_vector (x0))
- error ("initial: x0 must be a vector with %d elements", n);
+ error ("lsim: x0 must be a vector with %d elements", n);
endif
x = reshape (x0, [], 1); # make sure that x is a column vector
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-07-29 13:17:45
|
Revision: 10783
http://octave.svn.sourceforge.net/octave/?rev=10783&view=rev
Author: paramaniac
Date: 2012-07-29 13:17:39 +0000 (Sun, 29 Jul 2012)
Log Message:
-----------
control-devel: minor changes
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/INDEX
trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m
trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
Modified: trunk/octave-forge/extra/control-devel/INDEX
===================================================================
--- trunk/octave-forge/extra/control-devel/INDEX 2012-07-29 12:09:22 UTC (rev 10782)
+++ trunk/octave-forge/extra/control-devel/INDEX 2012-07-29 13:17:39 UTC (rev 10783)
@@ -26,4 +26,5 @@
@iddata/subsref
@iddata/vertcat
Miscellaneous
+ options
test_devel
\ No newline at end of file
Modified: trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-07-29 12:09:22 UTC (rev 10782)
+++ trunk/octave-forge/extra/control-devel/devel/HeatingSystem.m 2012-07-29 13:17:39 UTC (rev 10783)
@@ -72,7 +72,7 @@
dat = iddata (Y, U, 2.0, 'inname', 'input drive voltage', \
'inunit', 'Volt', \
'outname', 'temperature', \
- 'outunit', '°C')
+ 'outunit', 'Degree Celsius')
% s=15, n=7
[sys1, x0] = moen4 (dat, 's', 15, 'n', 7)
Modified: trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m 2012-07-29 12:09:22 UTC (rev 10782)
+++ trunk/octave-forge/extra/control-devel/devel/LakeErieARX.m 2012-07-29 13:17:39 UTC (rev 10783)
@@ -150,6 +150,6 @@
xlabel ('Time [months]')
xlim ([0, 56])
-legend ('measurement DaISy', 'Kalman Predictor', 'Kalman Predictor weak', 'location', 'northeast')
+legend ('Measurement DaISy', 'MOEN4 Kalman Predictor', 'MOEN4 Kalman Predictor (weak)', 'location', 'northeast')
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-07-29 12:09:29
|
Revision: 10782
http://octave.svn.sourceforge.net/octave/?rev=10782&view=rev
Author: paramaniac
Date: 2012-07-29 12:09:22 +0000 (Sun, 29 Jul 2012)
Log Message:
-----------
control: improve argument checks
Modified Paths:
--------------
trunk/octave-forge/main/control/inst/@lti/feedback.m
Modified: trunk/octave-forge/main/control/inst/@lti/feedback.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/feedback.m 2012-07-28 22:00:12 UTC (rev 10781)
+++ trunk/octave-forge/main/control/inst/@lti/feedback.m 2012-07-29 12:09:22 UTC (rev 10782)
@@ -110,6 +110,14 @@
endswitch
+ if (! is_real_vector (feedin) || ! isequal (feedin, abs (fix (feedin))))
+ error ("feedback: require 'feedin' to be a vector of integers");
+ endif
+
+ if (! is_real_vector (feedout) || ! isequal (feedout, abs (fix (feedout))))
+ error ("feedback: require 'feedout' to be a vector of integers");
+ endif
+
[p2, m2] = size (sys2);
l_feedin = length (feedin);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-07-28 22:00:22
|
Revision: 10781
http://octave.svn.sourceforge.net/octave/?rev=10781&view=rev
Author: benjf5
Date: 2012-07-28 22:00:12 +0000 (Sat, 28 Jul 2012)
Log Message:
-----------
Moving files to prepare for a first release. Not yet ready; lacks makefile and other files.
Added Paths:
-----------
trunk/octave-forge/extra/lssa/DESCRIPTION
trunk/octave-forge/extra/lssa/inst/
trunk/octave-forge/extra/lssa/inst/SampleDataSet.m
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/lscomplexwavelet.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
trunk/octave-forge/extra/lssa/inst/lswaveletcoeff.m
trunk/octave-forge/extra/lssa/src/
trunk/octave-forge/extra/lssa/src/fastlscomplex.cc
trunk/octave-forge/extra/lssa/src/fastlsreal.cc
Removed Paths:
-------------
trunk/octave-forge/extra/lssa/SampleDataSet.m
trunk/octave-forge/extra/lssa/cubicwgt.m
trunk/octave-forge/extra/lssa/fastlscomplex.cc
trunk/octave-forge/extra/lssa/fastlscomplex.cpp
trunk/octave-forge/extra/lssa/fastlsreal.cc
trunk/octave-forge/extra/lssa/lombcoeff.m
trunk/octave-forge/extra/lssa/lombnormcoeff.m
trunk/octave-forge/extra/lssa/lscomplex.m
trunk/octave-forge/extra/lssa/lscomplexwavelet.m
trunk/octave-forge/extra/lssa/lscorrcoeff.m
trunk/octave-forge/extra/lssa/lsreal.m
trunk/octave-forge/extra/lssa/lsrealwavelet.m
trunk/octave-forge/extra/lssa/lswaveletcoeff.m
Added: trunk/octave-forge/extra/lssa/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/lssa/DESCRIPTION (rev 0)
+++ trunk/octave-forge/extra/lssa/DESCRIPTION 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,14 @@
+Name: lssa
+Version:
+Date: 2012-07-28
+Author: Ben Lewis
+Maintainer: Ben Lewis
+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).
+Url: http://www.jstatsoft.org/v11/i02
+Problems: fast implementations, wavelet functions are currently not functional.
+Depends:
+License: GPL version 2 or later
Deleted: trunk/octave-forge/extra/lssa/SampleDataSet.m
===================================================================
--- trunk/octave-forge/extra/lssa/SampleDataSet.m 2012-07-28 10:41:08 UTC (rev 10780)
+++ trunk/octave-forge/extra/lssa/SampleDataSet.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -1,25 +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 2 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/>.
-
-## No function structure to this, I just want to use it to store the
-## sums of sines and cosines I'll use for testing.
-
-xvec = linspace(0,8,1000);
-maxfreq = 4 / ( 2 * pi );
-
-yvec = ( 2.*sin(maxfreq.*xvec) + 3.*sin((3/4)*maxfreq.*xvec)
- - 0.5 .* sin((1/4)*maxfreq.*xvec) - 0.2 .* cos(maxfreq .* xvec)
- + cos((1/4)*maxfreq.*xvec));
-
Deleted: trunk/octave-forge/extra/lssa/cubicwgt.m
===================================================================
--- trunk/octave-forge/extra/lssa/cubicwgt.m 2012-07-28 10:41:08 UTC (rev 10780)
+++ trunk/octave-forge/extra/lssa/cubicwgt.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -1,37 +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} {a =} cubicwgt {series}
-## Return the series as windowed by a cubic polynomial,
-## 1 + ( x ^ 2 * ( 2 x - 3 ) ), assuming x is in [-1,1].
-## @end deftypefn
-
-%!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);
-%! ## Tests cubicwgt on two scalars and two vectors; cubicwgt will work on any array input.
-
-
-## This function implements the windowing function on page 10 of the doc.
-## 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);
-endfunction
Deleted: trunk/octave-forge/extra/lssa/fastlscomplex.cc
===================================================================
--- trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-07-28 10:41:08 UTC (rev 10780)
+++ trunk/octave-forge/extra/lssa/fastlscomplex.cc 2012-07-28 22:00:12 UTC (rev 10781)
@@ -1,361 +0,0 @@
-/* Copyright (C) 2012 Benjamin Lewis <be...@gm...>
- * GNU GPLv2
- */
-
-
-#include <octave/oct.h>
-#include <octave/unwind-prot.h>
-#include <complex>
-#include <string>
-#include <math.h>
-#include <iostream>
-#include <exception>
-
-
-
-ComplexRowVector flscomplex( RowVector tvec , ComplexRowVector xvec ,
- double maxfreq , int octaves , int coefficients);
-
-
-DEFUN_DLD(fastlscomplex,args,nargout,
- "-*- texinfo -*-\n"
-"@deftypefn {Function File} { C = } fastlscomplex"
- "(@var{time},@var{magnitude},@var{maximum_frequency},@var{octaves},@var{coefficients})\n"
-"\n"
-"Return the complex least squares transform of the (@var{time},@var{magnitude}) series\n\
-supplied, using the fast algorithm.\n"
-"\n"
-"@seealso{lscomplex}\n"
-"@seealso{fastlsreal}\n"
-"\n"
-"@end deftypefn") {
- if ( args.length() != 5 ) {
- print_usage();
- return octave_value_list ();
- }
- RowVector tvals = args(0).row_vector_value();
- ComplexRowVector xvals = args(1).complex_row_vector_value();
- double omegamax = args(2).double_value();
- int noctaves = args(3).int_value();
- int ncoeff = args(4).int_value();
- if ( tvals.numel() != xvals.numel() ){
- if ( tvals.numel() > xvals.numel() ) {
- error("More time values than magnitude values.");
- } else {
- error("More magnitude values than time values.");
- }
- }
- if ( ncoeff == 0 ) error("No coefficients to compute.");
- if ( noctaves == 0 ) error("No octaves to compute over.");
- if ( omegamax == 0 ) error("No difference between minimal and maximal frequency.");
- octave_value_list retval;
- if ( !error_state) {
- ComplexRowVector results = flscomplex(tvals,xvals,omegamax,noctaves,ncoeff);
- retval(0) = octave_value(results);
- } else {
- return octave_value_list ();
- }
- return retval;
-
-}
-
-ComplexRowVector flscomplex( RowVector tvec , ComplexRowVector xvec ,
- double maxfreq, int octaves, int coefficients ) {
- struct Precomputation_Record {
- Precomputation_Record *next;
- std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only.
- bool stored_data;
- };
-
- ComplexRowVector results = ComplexRowVector (coefficients * octaves );
-
- double tau, delta_tau, tau_0, tau_h, n_inv, mu, te,
- omega_oct, omega_multiplier, octavemax, omega_working,
- loop_tau_0, loop_delta_tau, on_1, n_1, o;
- double length = ( tvec((tvec.numel()-1)) - tvec( octave_idx_type (0)));
- int octave_iter, coeff_iter;
- std::complex<double> zeta, zz, z_accumulator, exp_term, exp_multiplier, alpha, h, *tpra, *temp_ptr_alpha, temp_alpha[12], *tprb, *temp_ptr_beta, temp_beta[12], temp_array[12], *p, x;
- octave_idx_type n = tvec.numel();
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- temp_array[array_iter] = std::complex<double> ( 0 , 0 );
- }
- int factorial_array[12];
- factorial_array[0] = 1;
- for ( int i = 1 ; i < 12 ; i++ ) {
- factorial_array[i] = factorial_array[i-1] * i;
- }
- n_1 = n_inv = 1.0 / n;
- mu = (0.5 * M_PI)/length; // Per the article; this is in place to improve numerical accuracy if desired.
- /* Viz. the paper, in which Dtau = c / omega_max, and c is stated as pi/2 for floating point processors,
- * In the case of this computation, I'll go by the recommendation.
- */
- delta_tau = M_PI / ( 2 * maxfreq );
- tau_0 = tvec(0) + delta_tau;
- tau_h = tau_0;
- te = tau_h + delta_tau;
-
- octave_idx_type k ( 0 ); // Iterator for accessing xvec, tvec.
-
- Precomputation_Record * precomp_records_head, *record_current, *record_tail, *record_ref, *record_next;
- record_current = precomp_records_head = new Precomputation_Record;
- for ( te = tvec(k) + (2 * delta_tau) ; ; ) {
- x = std::complex<double>(xvec(k));
- {
- double t = mu*(tvec(k)-tau_h), tt;
- p = record_current->power_series;
- // p = 0
- *p++ = std::complex<double>(x);
- // p = 1
- tt = -t;
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- // p = 2
- tt *= t*(1.0/2.0);
- *p++ = x*tt;
- // p = 3
- tt *= t*(-1.0/3.0);
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- // p = 4
- tt *= t*(1.0/4.0);
- *p++ = x*tt;
- // p = 5
- tt *= t*(-1.0/5.0);
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- // p = 6
- tt *= t*(1.0/6.0);
- *p++ = x*tt;
- // p = 7
- tt *= t*(-1.0/7.0);
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- // p = 8
- tt *= t*(1.0/8.0);
- *p++ = x*tt;
- // p = 9
- tt *= t*(-1.0/9.0);
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- // p = 10
- tt *= t*(1.0/10.0);
- *p++ = x*tt;
- // p = 11
- tt *= t*(-1.0/11.0);
- h = x * tt;
- *p++ = std::complex<double>(-h.imag(),h.real());
- }
- record_current->stored_data = true;
- for(k++; ( k < n ) && tvec(k) < te ; k++ ) {
- x = std::complex<double>(xvec(k));
- {
- double t = mu*(tvec(k)-tau_h), tt;
- p = record_current->power_series;
- // p = 0
- *p++ += std::complex<double>(x);
- // p = 1
- tt = -t;
- h = x * tt;
- *p++ += std::complex<double>(- h.imag(), h.real());
- // p = 2
- tt *= t*(1.0/2.0);
- *p++ += x*tt;
- // p = 3
- tt *= t*(-1.0/3.0);
- h = x * tt;
- *p++ += std::complex<double>(-h.imag(),h.real());
- // p = 4
- tt *= t*(1.0/4.0);
- *p++ += x*tt;
- // p = 5
- tt *= t*(-1.0/5.0);
- h = x * tt;
- *p++ += std::complex<double>(-h.imag(),h.real());
- // p = 6
- tt *= t*(1.0/6.0);
- *p++ += x*tt;
- // p = 7
- tt *= t*(-1.0/7.0);
- h = x * tt;
- *p++ += std::complex<double>(-h.imag(),h.real());
- // p = 8
- tt *= t*(1.0/8.0);
- *p++ += x*tt;
- // p = 9
- tt *= t*(-1.0/9.0);
- h = x * tt;
- *p++ += std::complex<double>(-h.imag(),h.real());
- // p = 10
- tt *= t*(1.0/10.0);
- *p++ += x*tt;
- // p = 11
- tt *= t*(-1.0/11.0);
- h = x * tt;
- *p++ += std::complex<double>(-h.imag(),h.real());
- }
- record_current->stored_data = true;
- }
- if( k >= n ) break;
- tau_h = te + delta_tau;
- te = tau_h + delta_tau;
- record_current->next = new Precomputation_Record;
- record_current = record_current->next;
- }
- record_tail = record_current;
- record_current = precomp_records_head;
- record_tail->next = 0;
-
- /* Summation of coefficients for each frequency. As we have ncoeffs * noctaves elements,
- * it makes sense to work from the top down, as we have omega_max by default (maxfreq)
- */
-
- omega_oct = maxfreq / mu;
- omega_multiplier = exp(-log(2)/coefficients);
- octavemax = maxfreq;
- loop_tau_0 = tau_0;
- loop_delta_tau = delta_tau;
-
- octave_idx_type iter ( 0 );
-
- // Loops need to first travel over octaves, then coefficients;
-
- for ( octave_iter = octaves ; ; omega_oct *= 0.5 , octavemax *= 0.5 , loop_tau_0 += loop_delta_tau , loop_delta_tau *= 2 ) {
- o = omega_oct;
- omega_working = octavemax;
- for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, o *= omega_multiplier, omega_working *= omega_multiplier){
- exp_term = std::complex<double> ( cos( - omega_working * loop_tau_0 ) ,
- sin ( - omega_working * loop_tau_0 ) );
- exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- for ( zeta = 0, record_current = precomp_records_head ; record_current ;
- record_current = record_current->next, exp_term *= exp_multiplier ) {
- if ( record_current->stored_data ) {
- int p;
- for ( zz = 0 , p = 0, on_1 = n_1 ; p < 12 ; p++ ) {
- zz += record_current->power_series[p] * on_1 ;
- on_1 *= o;
- }
- zeta += exp_term * zz;
- }
- }
- results(iter) = std::complex<double> (zeta);
- iter++;
- }
- if ( !(--octave_iter) ) break;
- /* If we've already reached the lowest value, stop.
- * Otherwise, merge with the next computation range.
- */
- double *exp_pse_ptr, *exp_ptr, exp_power_series_elements[12];
- exp_power_series_elements[0] = 1;
- exp_pse_ptr = exp_ptr = exp_power_series_elements;
- for ( int r_iter = 1 ; r_iter < 12 ; r_iter++ ) {
- exp_power_series_elements[r_iter] = exp_power_series_elements[r_iter-1]
- * ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) );
- }
- try{
- for ( record_current = precomp_records_head ; record_current ;
- record_current = record_current->next ) {
- if ( ! ( record_ref = record_current->next ) || ! record_ref->stored_data ) {
- // In this case, there is no next record, but this record has data.
- if ( record_current->stored_data ) {
- int p = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1 , temp_ptr_alpha = temp_alpha ; p < 12 ; p++ , exp_pse_ptr++ ) {
- tpra = temp_ptr_alpha;
- *(temp_ptr_alpha++) = std::complex<double>(record_current->power_series[p]);
- for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *temp_ptr_alpha * *exp_ptr; ; ) {
- /* This next block is from Mathias' code, and it does a few
- * ... unsavoury things. First off, it uses conditionals with
- * break in order to avoid potentially accessing null regions
- * of memory, and then it does ... painful things with a few
- * numbers. However, remembering that most of these will not
- * actually be accessed for the first iterations, it's easier.
- */
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = *tpra * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_current->power_series[p] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = -*tpra * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_current->power_series[p] += *tpra * *exp_ptr;
- }
- }
- if ( ! record_ref ) break; // Last record was reached
- }
- else {
- record_next = record_ref;
- if ( record_current->stored_data ) {
- int p = 0, q = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; p < 12 ; p++, q++, exp_pse_ptr++ ) {
- tpra = temp_ptr_alpha;
- *temp_ptr_alpha++ = record_current->power_series[p] + record_next->power_series[q];
- *temp_ptr_beta++ = record_current->power_series[p] - record_next->power_series[1];
- tprb = temp_ptr_beta;
- for( exp_ptr = exp_power_series_elements, record_current->power_series[p] = *tpra * *exp_ptr; ; ) {
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tprb -= 2;
- h = *tprb * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tpra -= 2;
- record_current->power_series[p] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tprb -= 2;
- h = - *tprb * *exp_ptr;
- record_current->power_series[p].real() -= h.imag();
- record_current->power_series[p].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- tpra -= 2;
- record_current->power_series[p] += *tpra * *exp_ptr;
- }
- }
- } else {
- int q = 0;
- for( exp_pse_ptr = exp_power_series_elements + 1, temp_ptr_alpha = temp_alpha, temp_ptr_beta = temp_beta; q < 12; q++, exp_pse_ptr++ ) {
- tpra = temp_ptr_alpha;
- *temp_ptr_alpha++ = std::complex<double>(record_next->power_series[q]);
- for ( exp_ptr = exp_power_series_elements, record_next->power_series[q] = *tpra * *exp_ptr; ; ) {
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = *tpra * *exp_ptr;
- record_next->power_series[q].real() -= h.imag();
- record_next->power_series[q].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_next->power_series[q] -= *tpra * *exp_ptr;
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- h = -*tpra * *exp_ptr;
- record_next->power_series[q].real() -= h.imag();
- record_next->power_series[q].imag() += h.real();
- if ( ++exp_ptr >= exp_pse_ptr ) break;
- --tpra;
- record_next->power_series[q] += *tpra * *exp_ptr;
- }
- }
- }
- record_current->stored_data = true;
- record_ref = record_next;
- record_current->next = record_ref->next;
- delete record_ref;
- }
- }
- }
- } catch (std::exception& e) { //This section was part of my debugging, and may be removed.
- std::cout << "Exception thrown: " << e.what() << std::endl;
- ComplexRowVector exception_result (1);
- exception_result(0) = std::complex<double> ( 0,0);
- return exception_result;
- }
- }
- return results;
-}
Deleted: trunk/octave-forge/extra/lssa/fastlscomplex.cpp
===================================================================
--- trunk/octave-forge/extra/lssa/fastlscomplex.cpp 2012-07-28 10:41:08 UTC (rev 10780)
+++ trunk/octave-forge/extra/lssa/fastlscomplex.cpp 2012-07-28 22:00:12 UTC (rev 10781)
@@ -1,240 +0,0 @@
-/*
- * fastlscomplex.cpp, compiles to fastlscomplex.oct
- * Conversion to C++, with wrapper for Octave of the code from
- * A. Mathias' nuspectral package.
- */
-
-#include <octave/oct.h>
-#include <math.h>
-#include <complex>
-#include <string>
-#include <iostream>
-
-ComplexRowVector fastlscomplex ( RowVector tvals, ComplexRowVector xvals, octave_idx_type n,
- double length, int ncoeff, int noctaves, double omegamax );
-
-inline double sqr(double x)
-{ return x*x; }
-
-
-#define SETXT(p_, op_, x_, t_) (p_)->x op_ x_; (p_++)->t op_ t_;
-#define SETT(p_, op_, x_, t_) *p_++ op_ t_;
-#define SETX(p_, op_, x_, t_) *p_++ op_ x_;
-/* h is a complex aux. variable; it is used for assignment times I everywhere */
-#define SETIX(p_, op_, x_, t_) h = x_; (*(p_)).real() op_ -(h.imag()); (*(p_)).imag() op_ h.real(); p_++;
-
- /* Macro that sums up the power series terms into the power series
- * element record pointed to by p_.
- * By using = and += for op_, initial setting and accumulation can be selected.
- * t_ is the expression specifying the abscissa value. set_ can be either
- * SETXT to set the x and t fields of an XTElem record, or SETT/SETX to set
- * the elements of a Real array representing alternately real and imaginary
- * values.
- */
- // -10 points, comments don't match method.
-#define EXP_IOT_SERIES(p_, el_, t_, op_, setr_, seti_) \
-{ double t = t_, tt; p_ = el_; setr_(p_, op_, x, 1) \
- tt = -t; seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/2.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(-1.0/3.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/4.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(-1.0/5.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/6.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(-1.0/7.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/8.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(-1.0/9.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/10.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(-1.0/11.0); seti_(p_, op_, x*tt, tt) \
-}
-
-/* same as the above, but without alternating signs */
-#define EXPIOT_SERIES(p_, el_, t_, op_, setr_, seti_) \
-{ double t = t_, tt; p_ = el_; setr_(p_, op_, x, 1) \
- seti_(p_, op_, x*t, t ) \
- tt = t*t*(1.0/2.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/3.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/4.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/5.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/6.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/7.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/8.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/9.0); seti_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/10.0); setr_(p_, op_, x*tt, tt) \
- tt *= t*(1.0/11.0); seti_(p_, op_, x*tt, tt) \
-}
-
-# define SRCARG double *tptr, double *xptr, int n, double *lengthptr
-# define SRCVAR int k; double length = *lengthptr;
-
-//I'll remove these very shortly
-# define SRCFIRST k = 0
-# define SRCAVAIL (k<n)
-# define SRCNEXT k++
-
-DEFUN_DLD(fastlscomplex, args, nargout, "Takes the fast complex least-squares transform of irregularly sampled data.\n\
-When called, takes a time series with associated x-values, a number of octaves,\n\
-a number of coefficients per octave, the maximum frequency, and a length term.\n\
-It returns a complex row vector of the transformed series." ){
- RowVector tvals = args(0).row_vector_value();
- ComplexRowVector xvals = args(1).complex_row_vector_value();
- double length = ( tvals((tvals.numel()-1)) - tvals( octave_idx_type (0)));
- int ncoeff = args(4).int_value();
- int noctaves = args(5).int_value();
- double omegamax = args(6).double_value();
-
- if ( error_state ) {
- //print an error statement, see if error_state can be printed or interpreted.
- return octave_value_list ();
- }
- if ( tvals.length() != xvals.length() ) {
- std::cout << "Unmatched time series elements." << std::endl;
- // return error state if possible?
- return octave_value_list ();
- }
- octave_idx_type n = tvals.numel(); /* I think this will actually return octave_idx_type.
- * In that case I'll change the signature of fastlscomplex. */
- if ( ( ncoeff == 0 ) || ( noctaves == 0 ) ) {
- std::cout << "Cannot operate over either zero octaves or zero coefficients." << std::endl;
- // return error state if possible
- return octave_value_list ();
- }
- // possibly include error for when omegamax isn't right?
- ComplexRowVector results = fastlscomplex(tvals, xvals, n, length, ncoeff, noctaves, omegamax);
- return octave_value_list ( (octave_value) results );
-}
-
-ComplexRowVector fastlscomplex ( RowVector tvals, ComplexRowVector xvals, octave_idx_type n,
- double length, int ncoeff, int noctaves, double omegamax ) {
- /* Singly-linked list which contains each precomputation record
- * count stores the number of samples for which power series elements
- * were added. This will be useful for accelerating computation by ignoring
- * unused entries.
- */
- struct SumVec {
- SumVec *next;
- std::complex<double> elements[12];
- int count; };
- SumVec *shead, *stail, *sp, *sq;
- double dtelems[12], /* power series elements of exp(-i dtau) */
- *dte, *r, /* Pointers into dtelems */
- tau, tau0, te, /* Precomputation range centers and range end */
- tau0d, /* tau_h of first summand range at level d */
- dtau = (0.5*M_PI)/omegamax,/* initial precomputation interval radius */
- dtaud, /* precomputation interval radius at d'th merging step */
- n_1 = 1.0/n, /* reciprocal of sample count */ // n is implicitly cast to double.
- ooct, o, omul, /* omega/mu for octave's top omega and per band, mult. factor */
- omegaoct, omega, /* Max. frequency of octave and current frequency */
- on_1, /* n_1*(omega/mu)^p, n_1*(2*omega/mu)^p */
- mu = (0.5*M_PI)/length, /* Frequency shift: a quarter period of exp(i mu t) on length */
- tmp;
- std::complex<double> zeta, zz, // Accumulators for spectral coefficients to place in complex<double>
- e, emul,
- x,
- h, eoelems[12], oeelems[12],
- *eop, *oep,
- *ep, *op,
- *p, *q, *pe;
- ComplexRowVector results (ncoeff*noctaves);
-
- int i , j; // Counters; coefficient and octave, respectively.
-
- octave_idx_type k = 0;
-
- // Subdivision and precomputation, reinvisioned in an OOWorld.
- tau = tvals(k) + dtau;
- te = tau + dtau;
- tau0 = tau;
- shead = stail = sp = (SumVec*) operator new (sizeof(SumVec));
- sp->next = 0;
- { te = tvals(k) + ( 2 * dtau );
- while ( k < n ) {
- EXP_IOT_SERIES(p, sp->elements, mu*(tvals(k)-tau), =, SETX, SETIX);
- sp->count = 1;
- //Sum terms and show that there has been at least one precomputation.
- // I will probably replace the macro with a better explanation.
- for(SRCNEXT; SRCAVAIL && tvals(k) < te; SRCNEXT) {
- x = xvals(k);
- EXP_IOT_SERIES(p,sp->elements,mu*(tvals(k)-tau), +=, SETX, SETIX);
- sp->count++;
- }
- if ( k >= n ) break;
- tau = te + dtau;
- te = tau + dtau;
- sp = (SumVec*) operator new (sizeof(*sp));
- stail->next = sp;
- stail = sp;
- sp->next = 0;
- sp->count = 0;
- } }
- // Defining starting values for the loops over octaves:
- ooct = omegamax / mu;
- omul = exp(-log(2)/ncoeff);
- omegaoct = omegamax;
- tau0d = tau0;
- dtaud = dtau;
- octave_idx_type iter ( 0 );
- // Looping over octaves
- for ( j = noctaves ; ; ooct *= 0.5 , omegaoct *= 0.5 , tau0d += dtaud , dtaud *= 2 ) {
- // Looping over&results per frequency
- for ( i = ncoeff, o = ooct, omega = omegaoct; i-- ; o *= omul, omega *= omul ) {
- e.real() = cos( - ( omega * tau0d ) ); e.imag() = sin( - ( omega * tau0d ) );
- // sets real, imag parts of e
- emul.real() = cos( - 2 * ( omega * dtaud ) ); emul.imag() = sin( - 2 * ( omega * dtaud ) );
- // sets real, imag parts of emul
- for( zeta = 0 , sp = shead; sp; sp = sp->next, e *= emul ) {
- if ( sp->count ) {
- zz = std::complex<double>(0.0,0.0);
- octave_idx_type i ( 0 );
- for ( p = sp->elements , on_1 = n_1 ; i < (octave_idx_type) 12 ; i++ ) {
- zz += *p++ * on_1;
- on_1 *= 0;
- }
- zeta += e * zz;
- }
- results(iter) = std::complex<double>(zeta.real(), zeta.imag());
- iter++;
- }
- if ( --j <= 0 ) break; //Avoids excess merging
-
- EXPIOT_SERIES(r, dtelems, mu*dtaud, =, SETT, SETT);
- for(sp = shead; sp; sp = sp->next){
- if(!(sq = sp->next) || !sq->count ) {
- for(p = sp->elements, eop = eoelems, dte = dtelems+1, pe = p+12; p < pe; p++, dte++)
- { ep = eop; *eop++ = *p;
- for(r = dtelems, *p = *ep * *r; ; )
- { if(++r>=dte) break; --ep; h = *ep * *r; (*p).real() -= h.imag(); (*p).imag() += h.real();
- if(++r>=dte) break; --ep; *p -= *ep * *r;
- if(++r>=dte) break; --ep; h = -*ep * *r; (*p).real() -= h.imag(); (*p).imag() += h.real();
- if(++r>=dte) break; --ep; *p += *ep * *r;
- }
- }
- if(!sq) break; /* reached the last precomputation range */
- }
- else
- if(sp->count)
- for(p = sp->elements, q = sq->elements, eop = eoelems, oep = oeelems, dte = dtelems+1, pe = p+12;
- p < pe; p++, q++, dte++)
- { ep = eop; *eop++ = *p+*q; *oep++ = *p-*q; op = oep;
- for(r = dtelems, *p = *ep * *r; ; )
- { if(++r>=dte) break; op -= 2; h = *op * *r; (*p).real() -= h.imag(); (*p).imag() += h.real();
- if(++r>=dte) break; ep -= 2; *p -= *ep * *r;
- if(++r>=dte) break; op -= 2; h = -*op * *r; (*p).real() -= h.imag(); (*p).imag() += h.real();
- if(++r>=dte) break; ep -= 2; *p += *ep * *r;
- }
- }
- else
- for(q = sq->elements, eop = eoelems, oep = oeelems, dte = dtelems+1, pe = q+12; q<pe; q++, dte++)
- { ep = eop; *eop++ = *q;
- for(r = dtelems, *q = *ep * *r; ; )
- { if(++r>=dte) break; --ep; h = *ep * *r; (*q).real() -= h.imag(); (*q).imag() += h.real();
- if(++r>=dte) break; --ep; *p -= *ep * *r;
- if(++r>=dte) break; --ep; h = -*ep * *r; (*q).real() -= h.imag(); (*q).imag() += h.real();
- if(++r>=dte) break; --ep; *q += *ep * *r;
- }
- }
-
- sp->count += sq->count; sp->next = sq->next; /* free(sq) if malloc'ed */
- }
- }
- }
-}
Deleted: trunk/octave-forge/extra/lssa/fastlsreal.cc
===================================================================
--- trunk/octave-forge/extra/lssa/fastlsreal.cc 2012-07-28 10:41:08 UTC (rev 10780)
+++ trunk/octave-forge/extra/lssa/fastlsreal.cc 2012-07-28 22:00:12 UTC (rev 10781)
@@ -1,347 +0,0 @@
-/* Copyright (C) 2012 Benjamin Lewis <be...@gm...>
- * Licensed under the GNU GPLv2
- */
-
-
-#include <octave/oct.h>
-#include <octave/unwind-prot.h>
-#include <complex>
-#include <string>
-#include <math.h>
-#include <iostream>
-#include <exception>
-
-ComplexRowVector flsreal( RowVector tvec , ComplexRowVector xvec ,
- double maxfreq , int octaves , int coefficients);
-
-
-DEFUN_DLD(fastlsreal,args,nargout,
- "-*- texinfo -*-\n\
-@deftypefn {Function File} { C = } fastlsreal(@var{time},@var{magnitude},@var{maximum_frequency},@var{octaves},@var{coefficients})\n\
-\n\
-Return the real least-sqaures spectral fit to the (@var{time},@var{magnitude})\n\
-data supplied, using the fast algorithm.\n\
-\n\
-@seealso{fastlscomplex}\n\
-@seealso{lsreal}\n\
-@end deftypefn") {
- if ( args.length() != 5 ) {
- print_usage();
- return octave_value_list ();
- }
- RowVector tvals = args(0).row_vector_value();
- ComplexRowVector xvals = args(1).complex_row_vector_value();
- double omegamax = args(2).double_value();
- int noctaves = args(3).int_value();
- int ncoeff = args(4).int_value();
- if ( tvals.numel() != xvals.numel() ){
- if ( tvals.numel() > xvals.numel() ) {
- error("More time values than magnitude values.");
- } else {
- error("More magnitude values than time values.");
- }
- }
- if ( ncoeff == 0 ) error("No coefficients to compute.");
- if ( noctaves == 0 ) error("No octaves to compute over.");
- if ( omegamax == 0 ) error("No difference between minimal and maximal frequency.");
- octave_value_list retval;
- if ( !error_state) {
- ComplexRowVector results = flsreal(tvals,xvals,omegamax,noctaves,ncoeff);
- retval(0) = octave_value(results);
- } else {
- return octave_value_list ();
- }
- return retval;
-
-}
-
-ComplexRowVector flsreal( RowVector tvec , ComplexRowVector xvec ,
- double maxfreq, int octaves, int coefficients ) {
- struct Precomputation_Record {
- Precomputation_Record *next;
- std::complex<double> power_series[12]; // I'm using 12 as a matter of compatibility, only.
- bool stored_data;
- };
-
- ComplexRowVector results = ComplexRowVector (coefficients * octaves );
-
- double tau, delta_tau, tau_0, tau_h, n_inv, mu,
- omega_oct, omega_multiplier, octavemax, omega_working,
- loop_tau_0, loop_delta_tau;
- double length = ( tvec((tvec.numel()-1)) - tvec( octave_idx_type (0)));
- int octave_iter, coeff_iter;
- std::complex<double> zeta, z_accumulator, zeta_exp_term, zeta_exp_multiplier, alpha,
- iota, i_accumulator, iota_exp_term, iota_exp_multiplier;
- octave_idx_type n = tvec.numel();
- std::complex<double> temp_array[12];
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- temp_array[array_iter] = std::complex<double> ( 0 , 0 );
- }
- int factorial_array[12];
- factorial_array[0] = 1;
- for ( int i = 1 ; i < 12 ; i++ ) {
- factorial_array[i] = factorial_array[i-1] * i;
- }
- n_inv = 1.0 / n;
- mu = (0.5 * M_PI)/length; // Per the article; this is in place to improve numerical accuracy if desired.
- /* Viz. the paper, in which Dtau = c / omega_max, and c is stated as pi/2 for floating point processors,
- * In the case of this computation, I'll go by the recommendation.
- */
- delta_tau = M_PI / ( 2 * maxfreq );
- tau_0 = tvec(0) + delta_tau;
- tau_h = tau_0;
- size_t precomp_subset_count = (size_t) ceil( ( tvec(tvec.numel()-1) - tvec(0) ) / ( 2 * delta_tau ) );
- // I've used size_t because it will work for my purposes without threatening undefined behaviour.
- const std::complex<double> im = std::complex<double> ( 0 , 1 ); //I seriously prefer C99's complex numbers.
-
- octave_idx_type k ( 0 ); // Iterator for accessing xvec, tvec.
-
- Precomputation_Record * complex_precomp_records_head, *complex_record_current,
- *complex_record_tail, *complex_record_ref, *complex_record_next, *iota_precomp_records_head,
- *iota_record_current, *iota_record_tail, *iota_record_ref, *iota_record_next;
- complex_record_current = complex_precomp_records_head = new Precomputation_Record;
- iota_record_current = iota_precomp_records_head = new Precomputation_Record;
- for ( size_t p_r_iter = 1 ; p_r_iter < precomp_subset_count ; p_r_iter++ ) {
- complex_record_current->next = new Precomputation_Record;
- iota_record_current->next = new Precomputation_Record;
- complex_record_current = complex_record_current->next;
- iota_record_current = iota_record_current->next;
- }
- // Error's past this point
- complex_record_tail = complex_record_current;
- iota_record_tail = iota_record_current;
- complex_record_current = complex_precomp_records_head;
- iota_record_current = iota_precomp_records_head;
- complex_record_tail->next = 0;
- iota_record_tail->next = 0;
- /* A test needs to be included for if there was a failure, but since
- * precomp_subset_count is of type size_t, it should be okay. */
- for( ; complex_record_current != 0 ; complex_record_current = complex_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- complex_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- } // To avoid any trouble down the line, although it is an annoyance.
- while ( (k < n) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double p;
- for ( int j = 0 ; j < 12 ; j++ ) {
- alpha.real() = xvec(k).real();
- alpha.imag() = xvec(k).imag();
- // Change to switches for easier manipulation, fewer tests. This is two tests where one will do.
- if ( !( j % 2 ) ) {
- if ( ! ( j % 4 ) ) {
- alpha.real() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- } else {
- if ( ! ( j % 3 ) ) {
- alpha.real() = -1 * xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = -1 * xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- } else {
- alpha.real() = xvec(k).imag() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- alpha.imag() = xvec(k).real() * pow(mu,j) * pow(tvec(k)-tau_h,j) / factorial_array[j];
- }
- }
- complex_record_current->power_series[j].real() += alpha.real();
- complex_record_current->power_series[j].imag() += alpha.imag();
- }
- // Computes each next step of the power series for the given power series element.
- // j was reused since it was a handy inner-loop variable, even though I used it twice here.
- complex_record_current->stored_data = true;
- k++;
- }
- tau_h += ( 2 * delta_tau );
- }
- // At this point all precomputation records have been
- // exhausted for complex records. Short-circuit is abused
- // to avoid overflow errors.
- // Reset k, tau_h to reset the process. I may rewrite
- // these loops to be one, since running twice as long to
- // do the same thing is painful. May also move to a switch
- // in the previous section too.
- k = 0;
- tau_h = tau_0;
- for( ; iota_record_current != 0 ; iota_record_current = iota_record_current->next ) {
- for ( int j = 0 ; j < 12 ; j++ ) {
- iota_record_current->power_series[j] = std::complex<double> ( 0 , 0 );
- }
- while( ( k < n ) && (abs(tvec(k)-tau_h) <= delta_tau) ) {
- double comps[12];
- iota_record_current->power_series[0].real() = 1;
- comps[0] = 1;
- for ( int j = 1 ; j < 12 ; j++ ) {
- comps[j] = comps[j-1] * mu * (tvec(k)-tau_h);
- switch ( j % 4 ) {
- case 0 :
- iota_record_current->power_series[j].real() += comps[j] / factorial_array[j] ;
- break;
- case 1:
- iota_record_current->power_series[j].imag() += comps[j] / factorial_array[j] ;
- break;
- case 2:
- iota_record_current->power_series[j].real() -= comps[j] / factorial_array[j] ;
- break;
- case 3:
- iota_record_current->power_series[j].imag() -= comps[j] / factorial_array[j] ;
- break;
- }
- }
- iota_record_current->stored_data = true;
- k++;
- }
- tau_h += ( 2 * delta_tau );
- }
-
-
- /* Summation of coefficients for each frequency. As we have ncoeffs * noctaves elements,
- * it makes sense to work from the top down, as we have omega_max by default (maxfreq)
- */
-
- omega_oct = maxfreq / mu;
- omega_multiplier = exp(-log(2)/coefficients);
- octavemax = maxfreq;
- loop_tau_0 = tau_0;
- loop_delta_tau = delta_tau;
-
- octave_idx_type iter ( 0 );
-
- double zeta_real_part = 0, zeta_imag_part = 0, zeta_real_part_accumulator = 0, zeta_imag_part_accumulator = 0,
- iota_real_part = 0, iota_imag_part = 0, iota_real_part_accumulator = 0, iota_imag_part_accumulator = 0;
-
- // Loops need to first travel over octaves, then coefficients;
-
- for ( octave_iter = octaves ; ; omega_oct *= 0.5 , octavemax *= 0.5 , loop_tau_0 += loop_delta_tau , loop_delta_tau *= 2 ) {
- omega_working = omega_oct;
- zeta_exp_term = std::complex<double> ( cos ( - omega_working * loop_tau_0 ) ,
- sin ( - omega_working * loop_tau_0 ) );
- zeta_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- iota_exp_term = std::complex<double> ( cos ( - 2 * omega_working * loop_tau_0 ) ,
- sin ( - 2 * omega_working * loop_tau_0 ) );
- iota_exp_multiplier = std::complex<double> ( cos ( - 2 * omega_working * loop_delta_tau ) ,
- sin ( - 2 * omega_working * loop_delta_tau ) );
- for ( coeff_iter = 0 ; coeff_iter < coefficients ; coeff_iter++, omega_working *= omega_multiplier){
- zeta_real_part_accumulator = 0;
- zeta_imag_part_accumulator = 0;
- zeta_real_part = 0;
- zeta_imag_part = 0;
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next, zeta_exp_term *= zeta_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- z_accumulator = ( pow(omega_working,array_iter) * complex_record_current->power_series[array_iter] );
- zeta_real_part_accumulator += z_accumulator.real();
- zeta_imag_part_accumulator += z_accumulator.imag();
- }
- zeta_real_part = zeta_real_part + ( zeta_exp_term.real() * zeta_real_part_accumulator - zeta_exp_term.imag() * zeta_imag_part_accumulator );
- zeta_imag_part = zeta_imag_part + ( zeta_exp_term.imag() * zeta_real_part_accumulator + zeta_exp_term.real() * zeta_imag_part_accumulator );
- }
- for ( iota_record_current = iota_precomp_records_head; iota_record_current ;
- iota_record_current = iota_record_current->next, iota_exp_term *= iota_exp_multiplier ) {
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- i_accumulator = ( pow(omega_working,array_iter) * iota_record_current->power_series[array_iter] );
- iota_real_part_accumulator += i_accumulator.real();
- iota_imag_part_accumulator += i_accumulator.imag();
- }
- iota_real_part = iota_real_part + ( iota_exp_term.real() * iota_real_part_accumulator - iota_exp_term.imag() * iota_imag_part_accumulator );
- iota_imag_part = iota_imag_part + ( iota_exp_term.imag() * iota_real_part_accumulator + iota_exp_term.real() * iota_imag_part_accumulator );
- }
- // c + i s = 2 * ( zeta-omega-conjugate - iota-2omega-conjuage * zeta-omega ) / ( 1 - abs(iota-2omega) ^ 2 )
- // (is what the results will be.)
- zeta_real_part *= n_inv;
- zeta_imag_part *= n_inv;
- iota_real_part *= n_inv;
- iota_imag_part *= n_inv;
- double result_real_part = 2 * ( zeta_real_part - iota_real_part * zeta_real_part - iota_imag_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- double result_imag_part = 2 * ( - zeta_imag_part + iota_imag_part * zeta_real_part - iota_real_part * zeta_imag_part )
- / ( 1 - iota_real_part * iota_real_part - iota_imag_part * iota_imag_part );
- results(iter) = std::complex<double> ( result_real_part , result_imag_part );
- iter++;
- }
- if ( !(--octave_iter) ) break;
- /* If we've already reached the lowest value, stop.
- * Otherwise, merge with the next computation range.
- */
- double exp_power_series_elements[12];
- exp_power_series_elements[0] = 1;
- for ( int r_iter = 1 ; r_iter < 12 ; r_iter++ ) {
- exp_power_series_elements[r_iter] = exp_power_series_elements[r_iter-1]
- * ( mu * loop_delta_tau) * ( 1.0 / ( (double) r_iter ) );
- }
- try{
- for ( complex_record_current = complex_precomp_records_head ; complex_record_current ;
- complex_record_current = complex_record_current->next ) {
- if ( ! ( complex_record_ref = complex_record_current->next ) || ! complex_record_ref->stored_data ) {
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q )];
- }
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * complex_record_current->power_series[p - ( 2 * q ) - 1];
- }
- }
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
- if ( ! complex_record_ref ) break; // Last record was reached
- }
- else {
- complex_record_next = complex_record_ref;
- if ( complex_record_current->stored_data ) {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q )] - complex_record_next->power_series[p - (2*q)] );
- }
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,q) * ( complex_record_current->power_series[p - ( 2 * q ) - 1] - complex_record_next->power_series[p - ( 2 * q ) - 1 ] );
- }
- }
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
- } else {
- std::complex<double> temp[12];
- for( int array_init = 0 ; array_init < 12 ; array_init++ ) { temp[array_init] = std::complex<double>(0,0); }
- for( int p = 0 ; p < 12 ; p ++ ) {
- double step_floor_r = floor( ( (double) p ) / 2.0 );
- double step_floor_i = floor( ( (double) ( p - 1 ) ) / 2.0 );
- for( int q = 0 ; q < step_floor_r ; q++ ) {
- temp[p] += exp_power_series_elements[2*q] * pow((double)-1,q) * complex_record_next->power_series[p - ( 2 * q )];
- }
- for( int q = 0 ; q <= step_floor_i ; q++ ) {
- temp[p] += im * exp_power_series_elements[2 * q + 1] * pow((double)-1,(q+1)) * complex_record_next->power_series[p - ( 2 * q ) - 1];
- }
- }
- for ( int array_iter = 0 ; array_iter < 12 ; array_iter++ ) {
- complex_record_current->power_series[array_iter].real() = temp[array_iter].real();
- complex_record_current->power_series[array_iter].imag() = temp[array_iter].imag();
- }
- }
- complex_record_current->stored_data = true;
- complex_record_ref = complex_record_next;
- complex_record_current->next = complex_record_ref->next;
- delete complex_record_ref;
- }
- }
- }
- } catch (std::exception& e) { //This section was part of my debugging, and may be removed.
- std::cout << "Exception thrown: " << e.what() << std::endl;
- ComplexRowVector exception_result (1);
- exception_result(0) = std::complex<double> ( 0,0);
- return exception_result;
- }
- }
- return results;
-}
Copied: trunk/octave-forge/extra/lssa/inst/SampleDataSet.m (from rev 10774, trunk/octave-forge/extra/lssa/SampleDataSet.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/SampleDataSet.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/SampleDataSet.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,25 @@
+## 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 2 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/>.
+
+## No function structure to this, I just want to use it to store the
+## sums of sines and cosines I'll use for testing.
+
+xvec = linspace(0,8,1000);
+maxfreq = 4 / ( 2 * pi );
+
+yvec = ( 2.*sin(maxfreq.*xvec) + 3.*sin((3/4)*maxfreq.*xvec)
+ - 0.5 .* sin((1/4)*maxfreq.*xvec) - 0.2 .* cos(maxfreq .* xvec)
+ + cos((1/4)*maxfreq.*xvec));
+
Copied: trunk/octave-forge/extra/lssa/inst/cubicwgt.m (from rev 10774, trunk/octave-forge/extra/lssa/cubicwgt.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/cubicwgt.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/cubicwgt.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,37 @@
+## 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} {a =} cubicwgt {series}
+## Return the series as windowed by a cubic polynomial,
+## 1 + ( x ^ 2 * ( 2 x - 3 ) ), assuming x is in [-1,1].
+## @end deftypefn
+
+%!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);
+%! ## Tests cubicwgt on two scalars and two vectors; cubicwgt will work on any array input.
+
+
+## This function implements the windowing function on page 10 of the doc.
+## 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);
+endfunction
Copied: trunk/octave-forge/extra/lssa/inst/lombcoeff.m (from rev 10774, trunk/octave-forge/extra/lssa/lombcoeff.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombcoeff.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lombcoeff.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,40 @@
+## 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 2 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} {c =} lombcoeff (time, mag, freq)
+##
+## Return the coefficient of the Lomb periodogram (unnormalized) for the
+## (@var{time},@var{mag}) series for the @var{freq} provided.
+##
+## @seealso{lombnormcoeff}
+## @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( 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 );
+
+
+function coeff = lombcoeff(T, X, o)
+ theta = atan2(sum(sin(2 .* o .* T )), sum(cos(2.*o.*T)))/ (2 * o );
+ coeff = ( sum(X .* cos(o .* T - theta))**2)/(sum(cos(o.*T-theta).**2)) + ( sum(X .* sin(o .* T - theta))**2)/(sum(sin(o.*T-theta).**2));
+end function
+
Copied: trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m (from rev 10774, trunk/octave-forge/extra/lssa/lombnormcoeff.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lombnormcoeff.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,30 @@
+## 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 2 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} {c =} lombnormcoeff (time, mag, 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.
+##
+## @end deftypefn
+
+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) ) );
+end
\ No newline at end of file
Copied: trunk/octave-forge/extra/lssa/inst/lscomplex.m (from rev 10774, trunk/octave-forge/extra/lssa/lscomplex.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplex.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lscomplex.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,44 @@
+## 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 2 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} {t =} lscomplex ( time, mag, maxfreq, numcoeff, 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
+
+%!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 );
+
+
+
+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.
+ transform = zeros(1,ncoeff*noctave);
+ o = omegamax;
+ 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
+ o *= omul; ## To advance the transform to the next coefficient in the octave
+ endfor
+
+endfunction
Copied: trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m (from rev 10777, trunk/octave-forge/extra/lssa/lscomplexwavelet.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lscomplexwavelet.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,66 @@
+## Copyright (C) 2012 Ben Lewis <be...@gm...>
+##
+## This code is part of GNU Octave.
+##
+## This software 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 2 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/>.
+
+
+
+function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, minimum_window_count, sigma)
+
+## 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.
+
+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.
+ o = 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
+ ## the beautiful part of this code is that if parts end up falling outside the
+ ## vector, it won't matter (although it's wasted computations.)
+ ## I may add a trap to avoid too many wasted cycles.
+ windowed_t = ((abs (T-win_t) < current_radius) .* T);
+ ## this will of course involve an additional large memory allocation, at least in the short term,
+ ## but it's faster to do the operation once on the time series, then multiply it by the data series.
+ zeta = sum( cubicwgt ((windowed_t - win_t) ./ current_radius) .* exp( - i *
+ o .* windowed_t ) .* X ) / sum( cubicwgt ((windowed_t - win_t ) ./
+ current_radius) .* exp ( -i * o .* windowed_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?
+ endfor
+ endfor
+
+endfor
+
+endfunction
+
Copied: trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m (from rev 10774, trunk/octave-forge/extra/lssa/lscorrcoeff.m)
===================================================================
--- trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m (rev 0)
+++ trunk/octave-forge/extra/lssa/inst/lscorrcoeff.m 2012-07-28 22:00:12 UTC (rev 10781)
@@ -0,0 +1,69 @@
+## Copyright (C) 2012 Benjamin Lewis <ben...
[truncated message content] |
|
From: <aba...@us...> - 2012-07-28 10:41:14
|
Revision: 10780
http://octave.svn.sourceforge.net/octave/?rev=10780&view=rev
Author: abarth93
Date: 2012-07-28 10:41:08 +0000 (Sat, 28 Jul 2012)
Log Message:
-----------
Modified Paths:
--------------
trunk/octave-forge/main/octcdf/inst/example_opendap.m
Modified: trunk/octave-forge/main/octcdf/inst/example_opendap.m
===================================================================
--- trunk/octave-forge/main/octcdf/inst/example_opendap.m 2012-07-28 07:34:47 UTC (rev 10779)
+++ trunk/octave-forge/main/octcdf/inst/example_opendap.m 2012-07-28 10:41:08 UTC (rev 10780)
@@ -48,10 +48,4 @@
close(nc);
-colormap(hsv);
-axis xy
-iamgesc(ssh);
-
-% or with yapso
-% pcolor(ssh);
-
+pcolor(x,y,ssh), shading flat,colorbar
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <par...@us...> - 2012-07-28 07:34:54
|
Revision: 10779
http://octave.svn.sourceforge.net/octave/?rev=10779&view=rev
Author: paramaniac
Date: 2012-07-28 07:34:47 +0000 (Sat, 28 Jul 2012)
Log Message:
-----------
control: Fixed an argument check in feedback.m which caused false positive error messages. It was a copy-paste mistake affecting non-square systems. (Thanks to Tony Olivo)
Modified Paths:
--------------
trunk/octave-forge/main/control/NEWS
trunk/octave-forge/main/control/inst/@lti/feedback.m
Modified: trunk/octave-forge/main/control/NEWS
===================================================================
--- trunk/octave-forge/main/control/NEWS 2012-07-28 00:40:20 UTC (rev 10778)
+++ trunk/octave-forge/main/control/NEWS 2012-07-28 07:34:47 UTC (rev 10779)
@@ -1,6 +1,16 @@
Summary of important user-visible changes for releases of the control package
===============================================================================
+control-2.3.53 Release Date: 2012-xx-yy Release Manager: Lukas Reichlin
+===============================================================================
+
+** feedback
+ Fixed an argument check which caused false positive error messages.
+ It was a copy-paste mistake affecting non-square systems.
+ (Thanks to Tony Olivo)
+
+
+===============================================================================
control-2.3.52 Release Date: 2012-06-25 Release Manager: Lukas Reichlin
===============================================================================
Modified: trunk/octave-forge/main/control/inst/@lti/feedback.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/feedback.m 2012-07-28 00:40:20 UTC (rev 10778)
+++ trunk/octave-forge/main/control/inst/@lti/feedback.m 2012-07-28 07:34:47 UTC (rev 10779)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
+## Copyright (C) 2009, 2010, 2011, 2012 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -67,7 +67,7 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.5
+## Version: 0.6
function sys = feedback (sys1, sys2, feedin, feedout, fbsign = -1)
@@ -127,7 +127,7 @@
error ("feedback: range of feedin indices exceeds dimensions of sys1");
endif
- if (any (feedin > p1 | feedin < 1))
+ if (any (feedout > p1 | feedout < 1))
error ("feedback: range of feedout indices exceeds dimensions of sys1");
endif
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <mtm...@us...> - 2012-07-28 00:40:26
|
Revision: 10778
http://octave.svn.sourceforge.net/octave/?rev=10778&view=rev
Author: mtmiller
Date: 2012-07-28 00:40:20 +0000 (Sat, 28 Jul 2012)
Log Message:
-----------
Fix rceps to work correctly on odd-length inputs
Reported by Sergei Steshenko on the mailing list.
Modified Paths:
--------------
trunk/octave-forge/main/signal/NEWS
trunk/octave-forge/main/signal/inst/rceps.m
Modified: trunk/octave-forge/main/signal/NEWS
===================================================================
--- trunk/octave-forge/main/signal/NEWS 2012-07-26 23:07:34 UTC (rev 10777)
+++ trunk/octave-forge/main/signal/NEWS 2012-07-28 00:40:20 UTC (rev 10778)
@@ -4,6 +4,8 @@
signal-1.1.4 Release Date: 2012-XX-XX Release Manager:
===============================================================================
+ ** The function `rceps' was fixed to work correctly with odd-length inputs.
+
** The following functions are new:
movingrms schtrig clustersegment
Modified: trunk/octave-forge/main/signal/inst/rceps.m
===================================================================
--- trunk/octave-forge/main/signal/inst/rceps.m 2012-07-26 23:07:34 UTC (rev 10777)
+++ trunk/octave-forge/main/signal/inst/rceps.m 2012-07-28 00:40:20 UTC (rev 10778)
@@ -34,18 +34,22 @@
if (nargin != 1)
print_usage;
end
- y = real(ifft(log(abs(fft(x)))));
+ f = abs(fft(x));
+ if (any (f == 0))
+ error ("rceps: the spectrum of x contains zeros, unable to compute real cepstrum");
+ endif
+ y = real(ifft(log(f)));
if nargout == 2
n=length(x);
if rows(x)==1
if rem(n,2)==1
- ym = [y(1), 2*y(2:n/2), zeros(1,n/2-1)];
+ ym = [y(1), 2*y(2:n/2+1), zeros(1,n/2)];
else
ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
endif
else
if rem(n,2)==1
- ym = [y(1,:); 2*y(2:n/2,:); zeros(n/2-1,columns(y))];
+ ym = [y(1,:); 2*y(2:n/2+1,:); zeros(n/2,columns(y))];
else
ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,columns(y))];
endif
@@ -77,6 +81,13 @@
%! assert(yt.', y, tol);
%! assert(xmt.', xm, tol);
+%% Test that an odd-length input produces an odd-length output
+%!test
+%! x = randn(33, 4);
+%! [y, xm] = rceps(x);
+%! assert(size(y) == size(x));
+%! assert(size(xm) == size(x));
+
%!demo
%! f0=70; Fs=10000; # 100 Hz fundamental, 10kHz sampling rate
%! a=real(poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3]))); # two formants
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-07-26 23:07:40
|
Revision: 10777
http://octave.svn.sourceforge.net/octave/?rev=10777&view=rev
Author: benjf5
Date: 2012-07-26 23:07:34 +0000 (Thu, 26 Jul 2012)
Log Message:
-----------
Added lscomplexwavelet, still working on it.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/lsrealwavelet.m
Added Paths:
-----------
trunk/octave-forge/extra/lssa/lscomplexwavelet.m
Added: trunk/octave-forge/extra/lssa/lscomplexwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/lscomplexwavelet.m (rev 0)
+++ trunk/octave-forge/extra/lssa/lscomplexwavelet.m 2012-07-26 23:07:34 UTC (rev 10777)
@@ -0,0 +1,66 @@
+## Copyright (C) 2012 Ben Lewis <be...@gm...>
+##
+## This code is part of GNU Octave.
+##
+## This software 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 2 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/>.
+
+
+
+function transform = lscomplexwavelet( t , x, omegamax, ncoeff, noctave, minimum_window_count, sigma)
+
+## 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.
+
+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.
+ o = 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
+ ## the beautiful part of this code is that if parts end up falling outside the
+ ## vector, it won't matter (although it's wasted computations.)
+ ## I may add a trap to avoid too many wasted cycles.
+ windowed_t = ((abs (T-win_t) < current_radius) .* T);
+ ## this will of course involve an additional large memory allocation, at least in the short term,
+ ## but it's faster to do the operation once on the time series, then multiply it by the data series.
+ zeta = sum( cubicwgt ((windowed_t - win_t) ./ current_radius) .* exp( - i *
+ o .* windowed_t ) .* X ) / sum( cubicwgt ((windowed_t - win_t ) ./
+ current_radius) .* exp ( -i * o .* windowed_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?
+ endfor
+ endfor
+
+endfor
+
+endfunction
+
Modified: trunk/octave-forge/extra/lssa/lsrealwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-26 13:15:46 UTC (rev 10776)
+++ trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-26 23:07:34 UTC (rev 10777)
@@ -10,7 +10,6 @@
function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number )
-len_tx = length (T)
omegamult = 2 ^ ( - 1/ ncoeff )
omegamult_inv = 1 / omegamult
@@ -63,7 +62,7 @@
window_min += 2 * current_radius;
## I remain hesitant about this value, since it is entirely possible necessary precision will be lost. Should I try to reduce that?
endfor
-
+ o *= omegamult;
endfor
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <nir...@us...> - 2012-07-26 13:15:57
|
Revision: 10776
http://octave.svn.sourceforge.net/octave/?rev=10776&view=rev
Author: nir-krakauer
Date: 2012-07-26 13:15:46 +0000 (Thu, 26 Jul 2012)
Log Message:
-----------
suppressed warning about automatic broadcasting in csaps
Modified Paths:
--------------
trunk/octave-forge/main/splines/inst/csaps.m
Modified: trunk/octave-forge/main/splines/inst/csaps.m
===================================================================
--- trunk/octave-forge/main/splines/inst/csaps.m 2012-07-25 17:11:56 UTC (rev 10775)
+++ trunk/octave-forge/main/splines/inst/csaps.m 2012-07-26 13:15:46 UTC (rev 10776)
@@ -1,17 +1,17 @@
## Copyright (C) 2012 Nir Krakauer
##
-## 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 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 2 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.
+## 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/>.
+## 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{yi} @var{p}] =} csaps(@var{x}, @var{y}, @var{p}, @var{xi}, @var{w}=[])
@@ -94,8 +94,14 @@
aa = bb = cc = dd = zeros (n+1, m);
aa(2:end, :) = a;
cc(3:n, :) = 6*p*u; #second derivative at endpoints is 0 [natural spline]
+ warn_state = warning ("query", "Octave:broadcast").state;
+ warning ("off", "Octave:broadcast"); #turn off warning message for automatic broadcasting
+ unwind_protect
dd(2:n, :) = diff(cc(2:(n+1), :)) ./ h;
bb(2:n, :) = diff(a) ./ h - (cc(2:n, :)/2).*h - (dd(2:n, :)/6).*(h.^2);
+ unwind_protect_cleanup
+ warning (warn_state, "Octave:broadcast");
+ end_unwind_protect
bb(1, :) = bb(2, :); #linear extension of splines
bb(n + 1, :) = bb(n, :);
aa(1, :) = a(1, :) - eps(x(1))*bb(1, :);
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <prn...@us...> - 2012-07-25 17:12:33
|
Revision: 10775
http://octave.svn.sourceforge.net/octave/?rev=10775&view=rev
Author: prnienhuis
Date: 2012-07-25 17:11:56 +0000 (Wed, 25 Jul 2012)
Log Message:
-----------
Added NEWS file + updated/reconstructed old news to it
Added Paths:
-----------
trunk/octave-forge/extra/java/NEWS
Added: trunk/octave-forge/extra/java/NEWS
===================================================================
--- trunk/octave-forge/extra/java/NEWS (rev 0)
+++ trunk/octave-forge/extra/java/NEWS 2012-07-25 17:11:56 UTC (rev 10775)
@@ -0,0 +1,55 @@
+Summary of important user-visible changes for releases of the Java package
+
+===============================================================================
+java-1.2.9 Release Date: 2012-07-25 Release Manager: Philip Nienhuis
+===============================================================================
+
+** Bug fixes & improvements:
+--- configure.base: properly enable java if JAVA_ARCH was set beforehand
+ (Thomas Weber)
+--- __java__.cc: restore the locale after initializing the JVM (S\xE9bastien
+ Villemot)
+
+--- All dialog functions: accept cell array for dialog text caption (ML-
+ compatible) (not for dialog title!) bug #36311
+
+** Added pre_install.m, which should warn before, and if needed abort, actual
+ installation step if JAVA_HOME wasn't properly set. If needed, pre_install.m
+ also adds the directory containing the Java executables to the PATH,
+ provided this resides in the expected place (i.e., JAVA_HOME/bin)
+
+===============================================================================
+java-1.2.8 Release Date: 2011-06-08 Release Manager: Philip Nienhuis
+===============================================================================
+
+** Bug fixes & improvements:
+--- Octave.java: Support for java.lang.short & java.lang.long
+
+--- __java__.cc: Fixed wrong input arg order for javaMethod()
+ Removed extra ")" in getcwd call
+
+** Updated docs and added doc generation scripts(Martin Hepperle)
+
+===============================================================================
+java-1.2.7 Release Date: 2010-03-04 Release Managers: Martin Hepperle &
+ Philip Nienhuis
+===============================================================================
+
+** Bug fixes & improvements
+--- Many many improvements by Martin Hepperle
+
+** Added extensive docs in pdf, html and texinfo format (Martin Hepperle)
+
+** Added TeX support (Martin Hepperle)
+
+** Added dialog methods: msgbox, errordlg, inputdlg, warndlg, listdlg, questdlg,
+ helpdlg + test script (Martin Hepperle)
+
+** Implemented static classpath (Martin Hepperle)
+
+** Added Matlab-compatible javaMethod and javaObject as wrappers for
+ java_invoke() and java_new() (Martin Hepperle)
+
+** Added javarmpath(), javamethods() (Martin Hepperle)
+
+** Added javamem()
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
|
From: <be...@us...> - 2012-07-25 15:06:01
|
Revision: 10774
http://octave.svn.sourceforge.net/octave/?rev=10774&view=rev
Author: benjf5
Date: 2012-07-25 15:05:47 +0000 (Wed, 25 Jul 2012)
Log Message:
-----------
lsrealwavelet written; not fully tested, but performs as expected thus far.
Modified Paths:
--------------
trunk/octave-forge/extra/lssa/lsrealwavelet.m
Modified: trunk/octave-forge/extra/lssa/lsrealwavelet.m
===================================================================
--- trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-25 04:48:43 UTC (rev 10773)
+++ trunk/octave-forge/extra/lssa/lsrealwavelet.m 2012-07-25 15:05:47 UTC (rev 10774)
@@ -8,18 +8,23 @@
## @end deffun
-function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number, sigma=0.1 )
+function transform = lsrealwavelet(T, X, maxfreq, ncoeff, noctave, t_min, t_max, minimum_window_number )
-delta_t = (t_max - t_min)/(t_subdiv-1)
-# fails if you're trying to use less than 2 windows. Why use the wavelet tool at that point?
-
len_tx = length (T)
omegamult = 2 ^ ( - 1/ ncoeff )
omegamult_inv = 1 / omegamult
-winradius = pi / ( sigma * maxfreq )
-current_radius = winradius;
+minimum_window_width = ( t_max - t_min ) / minimum_window_number;
+minimum_window_radius = minimum_window_width / 2;
+sigma = maxfreq * 2 ^ ( noctave ) / minimum_window_radius;
+## sigma needs to be such that | t _ k - t | = minimum_window_radius implies
+## that sigma * maxfrequency * 2 ^ ( - noctave ) * minimum_window_radius = 1
+
+## for a specific other frequency, sigma * frequency * window_radius = 1 means
+## window_radius = 1 / ( frequency * sigma )
+
+
o = maxfreq;
# zeta _ ( t , omega ) = sum(w(sigma omega (t_k - t )e^(-i omega (t_k - t))xi_k)
@@ -33,57 +38,34 @@
transform = cell(noctave*ncoeff,1);
-## 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.
-win_t = window_min + window_radius;
-
-windowed_t = ((abs (T-win_t) < window_radius) .* T);
-zeta = sum( cubicwgt ((windowed_t - win_t) ./ window_radius) .* exp( - i * o .* windowed_t ) .* X ) / sum ( cubicwgt( windowed_t ./ window_radius ) );
-iota = sum( cubicwgt ((windowed_t - win_t) ./ window_radius) .* exp( - i * 2 * o .* windowed_t) .* X ) / sum ( cubicwgt( windowed_t .* 2 * o ) )
-## this will of course involve an additional large memory allocation, at least in the short term,
-## but it's faster to do the operation once on the time series, then multiply it by the data series.
-
-
+for iter = 1:(noctave*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.
+ current_frequency = maxfreq * 2 ^ ( - iter / ncoeff );
+ current_radius = 1 / ( current_frequency * sigma );
+ current_window_number = ceil( ( t_max - t_min ) / current_radius);
-for iter_freq = 1:( noctave * ncoeff )
- t = t_min;
- tx_point = 1;
- for iter_elem = 1:t_subdiv
- if tx_point > len_tx
- break
- endif
- while ( tx_point+1 <= len_tx) && ((t-T(tx_point+1))<winradius)
- tx_point++;
- endwhile
- zeta = 0;
- iota = 0;
- iota0 = 0;
- count = 0;
- ot = o * ( T(tx_point) - t);
- sot = sigma * ot;
- while ( tx_point <= len_tx) && ( (sot = sigma * (ot = o*(T(tx_point)-t))) < pi )
- w = 0.5 * ( 1 + cos (sot));
- iota0 += w;
- zeta += w * ( cos (ot) + sin (ot) * i ) * X(tx_point);
- ot *= 2;
- iota += w * ( cos(ot) + sin (ot) * i );
- count++;
- tx_point++;
- endwhile
- if count > 0
- result = 2 * ( conj (zeta) * iota0 + zeta * conj (iota) ) / ( ( count + iota0 ) ** 2 - real(iota) ** 2 - imag(iota) ** 2 )
- printf("%d\n",result);
- transform(iter_freq,iter_elem) = result;
- else
- transform(iter_freq,iter_elem) = 0;
- endif
- t += delta_t
+ transform{iter} = zeros(1,current_window_number);
+ win_t = window_min + current_radius;
+ for iter_window = 1:current_window_number
+ ## the beautiful part of this code is that if parts end up falling outside the
+ ## vector, it won't matter (although it's wasted computations.)
+ ## I may add a trap to avoid too many wasted cycles.
+ windowed_t = ((abs (T-win_t) < current_radius) .* T);
+ ## this will of course involve an additional large memory allocation, at least in the short term,
+ ## but it's faster to do the operation once on the time series, then multiply it by the data series.
+ iota0 = sum ( cubicwgt (windowed_t ./ current_radius ) );
+ zeta = sum( cubicwgt ((windowed_t - win_t) ./ current_radius) .* exp( - i * o .* windowed_t ) .* X ) / iota0;
+ iota = sum( cubicwgt ((windowed_t - win_t) ./ current_radius) .* exp( - i * 2 * o .* windowed_t) .* X ) / sum ( cubicwgt( windowed_t .* 2 * o ) );
+ transform{iter}(iter_window) = 2 * ( conj(zeta) * iota0 + zeta * conj(iota) ) / ( ( length ( find (windowed_t)) + iota0 ) ^ 2 - real(iota) ^ 2 - imag(iota) ^ 2 );
+ window_min += 2 * current_radius;
+ ## I remain hesitant about this value, since it is entirely possible necessary precision will be lost. Should I try to reduce that?
endfor
- o *= omegamult
- winradius *= omegamult_inv
- iter_freq
+
endfor
+
endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
48 messages has been excluded from this view by a project administrator.
