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