Commit 937783f4 authored by Rafael Laboissiere's avatar Rafael Laboissiere

Imported Upstream version 0.1.2

parents
This diff is collapsed.
Name: lssa
Version: 0.1.2
Date: 2012-08-17
Author: Ben Lewis <benjf5@gmail.com>
Maintainer: Ben Lewis <benjf5@gmail.com>
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: octave (>= 3.6.0)
Autoload: no
License: GPLv3+
lssa >> Least Squares Spectral Analysis
Windowing
cubicwgt
Periodogram
lombcoeff lombnormcoeff
Accelerated time-series functions
fastlscomplex
Complex time-series functions
lscomplex
Real time-series functions
lsreal
Correlation
lscorrcoeff
Wavelet Transform
lswaveletcoeff
# lscomplexwavelet lsrealwavelet
## The wavelet functions are unavailable until I can get them working.
Summary of changes in lssa 0.1.2:
** All functions now have input checks in place to return useful errors as
opposed to division by zero, etc. Documentation has also been improved.
Summary of status of the intial lssa release, 0.1.1:
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.
** There are two wavelet functions under development, but they are not included
in this release as they are currently not functional. For all your wavelet
needs, the specific transformation used is available in the lswaveletcoeff
function, and will generate a single cosine/sine magnitude pair (as a
complex number) for a complex-valued series (this function may be joined by
a companion for real-valued series) and can be looped to simulate a full
wavelet transform.
** For all the working functions, tests have been written and formatted to
Octave coding standards. These tests should pass on any given architecture
(there was some question about that previously) and often provide examples
of how the function operates. For a few functions, there are demo scripts.
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{a} =} cubicwgt (@var{series})
##
## Returns the input series, windowed by a polynomial similar to a Hanning
## window. To window an arbitrary section of the series, subtract or add an
## offset to it to adjust the centre of the window; for an offset of k, the call
## would be cubicwgt (@var{s} - k). Similarly, the radius of the window is 1;
## if an arbitrary radius r is desired, dividing the series by the radius after
## centering is the best way to adjust to fit the window: cubicwgt ((@var{s} -
## k) / r).
##
## The windowing function itself is:
## w = 1 + ( x ^ 2 * ( 2 x - 3 ) ), x in [-1,1], else w = 0.
##
## @end deftypefn
function a = cubicwgt (s)
if (nargin != 1)
print_usage ();
endif
## s is the value/vector/matrix to be windowed
a = abs (s);
a = ifelse ((a < 1), 1 + ((a .^ 2) .* (2 .* a - 3)), 0);
endfunction
%!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.
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lombcoeff (@var{time}, @var{mag}, @var{freq})
##
## Return the Lomb Periodogram value at one frequency for a time series.
##
## @seealso{lombnormcoeff}
## @end deftypefn
function coeff = lombcoeff (T, X, o)
if (nargin != 3)
print_usage ();
endif
if (! all (size (T) == size (X)))
error ("lombcoeff: Time series vectors of uneven size.\n");
endif
if (! isscalar (o))
error ("lombcoeff: Supplied frequency is not a scalar.\n");
endif
if (o == 0)
error ("lombcoeff: Supplied frequency is not a frequency.\n");
endif
oT = o .* T;
theta = atan2 (sum (sin (2 * oT)),
sum (cos (2 * oT))) ./ (2 * o);
coeff = (sum (X .* cos (oT - theta)) ^2 /
sum (cos (oT - theta) .^2) +
sum (X .* sin (oT - theta)) ^2 /
sum (sin (oT - theta) .^2));
endfunction
%!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, maxfreq), 1076.77574184435, 5e-10);
%!assert (lombcoeff (t, x, 3/4*maxfreq), 1226.53572492183, 5e-10);
%!assert (lombcoeff (t, x, 1/4*maxfreq), 1341.63962181896, 5e-10);
## Copyright (c) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## 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} {@var{c} =} lombnormcoeff (@var{time}, @var{mag}, @var{freq})
##
## Return the normalized Lomb Periodogram value at one frequency for a time
## series.
##
## @seealso{lombcoeff}
##
## @end deftypefn
function coeff = lombnormcoeff (T, X, omega)
if (nargin != 3)
print_usage ();
endif
if (! all (size (T) == size (X)))
error ("lombnormcoeff: Time series vectors of uneven size.\n");
endif
if (! isscalar (omega))
error ("lombnormcoeff: Supplied frequency is not a scalar.\n");
endif
if (omega == 0)
error ("lombnormcoeff: Supplied frequency is not a frequency.\n");
endif
xmean = mean (X);
theta = atan2 (sum (sin (2 .* omega .*T)),
sum (cos (2 .* omega .* T))) / (2*omega);
coeff = ((sum ((X-xmean) .* cos (omega .* T - theta)) .^ 2 /
sum (cos (omega .* T - theta) .^ 2) +
sum ((X-xmean) .* sin (omega .* T - theta)) .^ 2 /
sum (sin (omega .* T - theta) .^ 2 )) /
(2 * var(X)));
endfunction
%!shared t, x, o, maxfreq
%! 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 (lombnormcoeff (t,x,o(1)), 44.7068607258824, 5e-10);
%!assert (lombnormcoeff (t,x,o(2)), 35.7769955188467, 5e-10);
%!assert (lombnormcoeff (t,x,o(3)), 20.7577786183241, 5e-10);
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{t} =} lscomplex (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
##
## Return a series of least-squares transforms of a complex-valued time series.
## Each transform is minimized independently at each frequency. @var{numcoeff}
## frequencies are tested for each of @var{numoctaves} octaves, starting from
## @var{maxfreq}.
##
## Each result (a + bi) at a given frequency, o, defines the real and imaginary
## coefficients for a sum of cosine and sine functions: a cos(ot) + b i
## sin(ot). The specific frequency can be determined by its index in @var{t},
## @var{ind}, as @var{maxfreq} * 2 ^ (- (@var{ind} - 1) / @var{numcoeff}).
##
## @seealso{lsreal}
## @end deftypefn
function transform = lscomplex (t, x, omegamax, ncoeff, noctave)
if (nargin != 5)
print_usage ();
endif
if (! isvector (t))
error ("lscomplex: Time values are not a vector.\n");
endif
if (! isvector (x))
error ("lscomplex: Magnitude values are not a vector.\n");
endif
if (! all (size (t) == size (x)))
error ("lscomplex: Size of time vector, magnitude vector unequal.\n");
endif
if (! isscalar (omegamax))
error ("lscomplex: More than one value for maximum frequency specified.\n");
endif
if (! isscalar (ncoeff))
error ("lscomplex: More than one number of frequencies per octave specified.\n");
endif
if (! isscalar (noctave))
error ("lscomplex: More than one number of octaves to traverse specified.\n");
endif
if (omegamax == 0)
error ("lscomplex: Specified maximum frequency is not a frequency.\n");
endif
if (noctave == 0)
error ("lscomplex: No octaves of results requested.\n");
endif
if (ncoeff == 0)
error ("lscomplex: No frequencies per octave requested.\n");
endif
if (ncoeff != floor (ncoeff))
error ("lscomplex: Specified number of frequencies per octave is not integral.\n");
endif
if (noctave != floor (noctave))
error ("lscomplex: Specified number of octaves of results is not integral.\n");
endif
n = numel (t);
iter = 0 : (ncoeff * noctave - 1);
omul = (2 .^ (- iter / ncoeff));
ot = t(:) * (omul * omegamax);
transform = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
endfunction
%!test
%! maxfreq = 4 / ( 2 * pi );
%! t = [0:0.008:8];
%! x = ( 2 .* sin (maxfreq .* t) +
%! 3 .* sin ( (3 / 4) * maxfreq .* t)-
%! 0.5 .* sin ((1/4) * maxfreq .* t) -
%! 0.2 .* cos (maxfreq .* t) +
%! cos ((1/4) * maxfreq .* t));
%! assert (lscomplex (t, x, maxfreq, 2, 2),
%! [(-0.400924546169395 - 2.371555305867469i), ...
%! (1.218065147708429 - 2.256125004156890i), ...
%! (1.935428592212907 - 1.539488163739336i), ...
%! (2.136692292751917 - 0.980532175174563i)], 5e-10);
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq})
## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} = @var{cubicwgt})
## @deftypefnx {Function File} {@var{c} =} lscorrcoeff (@var{time1}, @var{mag1}, @var{time2}, @var{mag2}, @var{time}, @var{freq}, @var{window} = @var{cubicwgt}, @var{winradius} = 1)
##
## Return the coefficient of the wavelet correlation of two complex time
## series. The correlation is only effective at a given time and frequency.
## The windowing function applied by default is cubicwgt, this can be changed by
## passing a different function handle to @var{window}, while the radius applied
## is set by @var{winradius}. Note that this will be most effective when both
## series have had their mean value (if it is not zero) subtracted (and stored
## separately); this reduces the constant-offset error further, and allows the
## functions to be compared on their periodic features rather than their
## constant features.
##
## @seealso{lswaveletcoeff, lscomplexwavelet, lsrealwavelet}
##
## @end deftypefn
function coeff = lscorrcoeff (x1, y1, x2, y2, t, o, wgt = @cubicwgt, wgtrad = 1)
## Input checking is absolutely necessary.
if (!((nargin >= 6) && (nargin <= 8)))
print_usage ();
endif
## Test to be sure x1, y1, x2, y2 are all vectors, and that t and o are
## scalars.
if (! isvector (x1))
error ("lscorrcoeff: First time series time values are not a vector.\n");
endif
if (! isvector (y1))
error ("lscorrcoeff: First time series magnitude values are not a vector.\n");
endif
if (! isvector (x2))
error ("lscorrcoeff: Second time series time values are not a vector.\n");
endif
if (! isvector (y2))
error ("lscorrcoeff: Second time series magnitude values are not a vector.\n");
endif
if (! isscalar (t))
error ("lscorrcoeff: Window centre is not a scalar.\n");
endif
if (! isscalar (o))
error ("lscorrcoeff: Specified frequency is not a scalar.\n");
endif
if (! isscalar (wgtrad))
error ("lscorrcoeff: Window radius is not a scalar.\n");
endif
if (! all (size (x1) == size (y1)))
error ("lscorrcoeff: First time series vectors not of matching size.\n");
endif
if (! all (size (x2) == size (y2)))
error ("lscorrcoeff: Second time series vectors not of matching size.\n");
endif
## How to determine if a weight function has been assigned or not? (Possible
## to get name of function?)
so = 0.05 * o;
## The first solution that comes to mind is admittedly slightly
## ugly and has a data footprint of O(2n) but it is vectorised.
mask = (abs (x1 - t) * so) < wgtrad;
rx1 = x1(mask);
## FIXME : Needs to have a noisy error if length(y1) != length(x1) -- add this!
ry1 = y1(mask);
mask = (abs (x2 - t) * so ) < wgtrad;
rx2 = x2(mask);
ry2 = y2(mask);
windowed_element_count = length (rx1);
if (windowed_element_count == 0)
error("lscorrcoeff: No time-series elements contained in window.\n");
endif
s = sum (wgt ((rx1 - t) .* so)) * sum (wgt ((rx2 - t ) .* so ));
if (s != 0)
coeff = sum (wgt ((rx1 - t) .* so) .* exp (i * o .* rx1) .* ry1) * ...
sum (wgt ((rx2 - t) .* so) .* exp (i * o .* rx2) .* conj (ry2)) / s;
else
coeff = 0;
endif
endfunction
%!shared t, p, x, y, z, o, maxfreq
%! maxfreq = 4 / (2 * pi);
%! t = linspace (0, 8);
%! x = (2 .* sin (maxfreq .* t) +
%! 3 .* sin ((3/4) * maxfreq .* t) -
%! 0.5 .* sin ((1/4) * maxfreq .* t) -
%! 0.2 .* cos (maxfreq .* t) +
%! cos ((1/4) * maxfreq .* t));
%! y = - x;
%! p = linspace (0, 8, 500);
%! z = (2 .* sin (maxfreq .* p) +
%! 3 .* sin ((3/4) * maxfreq .* p) -
%! 0.5 .* sin ((1/4) * maxfreq .* p) -
%! 0.2 .* cos (maxfreq .* p) +
%! cos ((1/4) * maxfreq .* p));
%! o = [maxfreq , (3/4 * maxfreq) , (1/4 * maxfreq)];
%!assert (lscorrcoeff (t, x, t, x, 0.5, maxfreq),
%! -5.54390340863576 - 1.82439880893383i, 5e-10);
%!assert (lscorrcoeff (t, x, t, y, 0.5, maxfreq),
%! 5.54390340863576 + 1.82439880893383i, 5e-10);
%!assert (lscorrcoeff (t, x, p, z, 0.5, maxfreq),
%! -5.55636741054624 - 1.82803733863170i, 5e-10);
## Demo with sin, cos as Nir suggested.
%!demo
%! ## This generates the correlation coefficient at time 0.5 and circular freq. 0.9
%! x = 1:10;
%! y = sin (x);
%! z = cos (x);
%! a = lscorrcoeff (x, y, x, z, 0.5, 0.9)
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{t} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
##
## Return a series of least-squares transforms of a real-valued time series.
## Each transform is minimized independently for each frequency. The method
## used is a Lomb-Scargle transform of the real-valued (@var{time}, @var{mag})
## series, starting from frequency @var{maxfreq} and descending @var{numoctaves}
## octaves with @var{numcoeff} coefficients per octave.
##
## The result of the transform for each frequency is the coefficient of a sum of
## sine and cosine functions modified by that frequency, in the form of a
## complex number—where the cosine coefficient is encoded in the real term, and
## the sine coefficient is encoded in the imaginary term. Each frequency is fit
## independently from the others, and to minimize very low frequency error,
## consider storing the mean of a dataset with a constant or near-constant
## offset separately, and subtracting it from the dataset.
##
## @seealso{lscomplex}
## @end deftypefn
function transform = lsreal (t, x, omegamax, ncoeff, noctave)
## Sanity checks to make sure that the user can get meaningful errors.
if (nargin != 5)
print_usage ();
endif
if (! isvector (t))
error ("lsreal: Time values are not a vector.\n");
endif
if (! isvector (x))
error ("lsreal: Magnitude values are not a vector.\n");
endif
if (! all (size (t) == size (x)))
error ("lsreal: Size of time vector, magnitude vector unequal.\n");
endif
if (! isscalar (omegamax))
error ("lsreal: More than one value for maximum frequency specified.\n");
endif
if (! isscalar (ncoeff))
error ("lsreal: More than one number of frequencies per octave specified.\n");
endif
if (! isscalar (noctave))
error ("lsreal: More than one number of octaves to traverse specified.\n");
endif
if (omegamax == 0)
error ("lsreal: Specified maximum frequency is not a frequency.\n");
endif
if (noctave == 0)
error ("lsreal: No octaves of results requested.\n");
endif
if (ncoeff == 0)
error ("lsreal: No frequencies per octave requested.\n");
endif
if (ncoeff != floor (ncoeff))
error ("lsreal: Specified number of frequencies per octave is not integral.\n");
endif
if (noctave != floor (noctave))
error ("lsreal: Specified number of octaves of results is not integral.\n");
endif
n = numel (t);
iter = 0 : (ncoeff * noctave - 1);
omul = (2 .^ (- iter / ncoeff));
## For a given frequency, the iota term is taken at twice the frequency of the
## zeta term.
ot = t(:) * (omul * omegamax);
oit = t(:) * (omul * omegamax * 2);
zeta = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
iota = sum ((cos (oit) - (sin (oit) .* i)), 1) / n;
transform = 2 .* (conj (zeta) - conj (iota) .* zeta) ./ (1 - abs (iota) .^ 2);
endfunction
%!test
%! maxfreq = 4 / ( 2 * pi );
%! t = linspace(0,8);
%! x = ( 2 .* sin ( maxfreq .* t ) +
%! 3 .* sin ( (3/4) * maxfreq .* t ) -
%! 0.5 .* sin ( (1/4) * maxfreq .* t ) -
%! 0.2 .* cos ( maxfreq .* t ) +
%! cos ( (1/4) * maxfreq .* t ) );
%! # In the assert here, I've got an error bound large enough to catch
%! # individual system errors which would present no real issue.
%! assert (lsreal (t,x,maxfreq,2,2),
%! [(-1.68275915310663 + 4.70126183846743i), ...
%! (1.93821553170889 + 4.95660209883437i), ...
%! (4.38145452686697 + 2.14403733658600i), ...
%! (5.27425332281147 - 0.73933440226597i)],
%! 5e-10)
## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq})
## @deftypefnx {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq}, @var{window}=cubicwgt)
## @deftypefnx {Function File} {@var{c} =} lswaveletcoeff (@var{t}, @var{x}, @var{time}, @var{freq}, @var{window}=cubicwgt, @var{winradius}=1)
##
## Return the wavelet transform of a complex time series in a given window. The
## transform takes a complex time series (@var{t}, @var{x}) at time @var{time}
## and frequency @var{freq}, then applies a windowing function to it; the
## default is cubicwgt, however by providing a function handle for the optional
## variable @var{window}, the user may select their own function; to determine
## the radius of the interval around the @var{time} selected, set
## @var{winradius} to some value other than 1.
##
## This transform operates identically to the transform at the heart of
## lscomplexwavelet, however for one window only.
##
## @seealso{lscorrcoeff, lscomplexwavelet, lsrealwavelet}
##
## @end deftypefn
function coeff = lswaveletcoeff (x, y, t, o, wgt = @cubicwgt, wgtrad = 1)
if (! (nargin >= 4) && (nargin <= 6))
print_usage ();
endif
if (! isvector (x))
error ("lswaveletcoeff: Time values are not a vector.\n");
endif
if (! isvector (y))
error ("lswaveletcoeff: Magnitude values are not a vector.\n");
endif
if (! all (size (x) == size (y)))
error ("lswaveletcoeff: Time series vectors of uneven size.\n");
endif
if (! isscalar (t))
error ("lswaveletcoeff: Window centre specified is not scalar.\n");
endif
if (! isscalar (o))
error ("lswaveletcoeff: Frequency specified is not scalar.\n");
endif
if (! isscalar (wgtrad))
error ("lswaveletcoeff: Window radius specified is not scalar.\n");
endif
so = 0.05 .* o;
if ((ndims (x) == 2) && ! (rows (x) == 1))
x = reshape (x, 1, length (x));
y = reshape (y, 1, length (y));
endif
mask = (abs (x - t) * so < wgtrad);
rx = x(mask);
ry = y(mask);
## Going by the R code, this can use the same mask.
s = sum (wgt ((x - t) .* so));
coeff = ifelse (s != 0,
sum (wgt ((rx - t) .* so) .*
exp (i .* o .* (rx - t)) .* ry) ./ s,
0);
endfunction
%!shared t, p, x, y, maxfreq
%! maxfreq = 4 / (2 * pi);
%! t = linspace (0, 8);
%! x = (2 .* sin (maxfreq .* t) +
%! 3 .* sin ((3/4) * maxfreq .* t) -
%! 0.5 .* sin ((1/4) * maxfreq .* t) -
%! 0.2 .* cos (maxfreq .* t) +
%! cos ((1/4) * maxfreq .* t));
%! y = - x;
%! p = linspace (0, 8, 500);
%!assert (lswaveletcoeff (t, x, 0.5, maxfreq),
%! 0.383340407638780 + 2.385251997545446i, 5e-10);
%!assert (lswaveletcoeff (t, y, 3.3, 3/4 * maxfreq),
%! -2.35465091096084 + 1.01892561714824i, 5e-10);
%!demo
%! ## Generates the wavelet transform coefficient for time 0.5 and circ. freq. 0.9, for row & column vectors.
%! x = 1:10;
%! y = sin (x);
%! xt = x';
%! yt = y';
%! a = lswaveletcoeff (x, y, 0.5, 0.9)
%! b = lswaveletcoeff (xt, yt, 0.5, 0.9)
MKOCTFILE ?= mkoctfile
all: fastlscomplex.oct #\
# fastlsreal.oct
fastlscomplex.oct: fastlscomplex.cc
$(MKOCTFILE) fastlscomplex.cc
# fastlsreal compilation is disabled for the time being
#fastlsreal.oct: fastlsreal.cc
# $(MKOCTFILE) fastlsreal.cc
# helper function just in case
clean:
rm *.o *.oct *~ octave-core
This diff is collapsed.
src/fastlsreal.cc 0 → 100644