1## Copyright (C) 2012 Benjamin Lewis <benjf5@gmail.com>
2##
3## This program is free software; you can redistribute it and/or modify it under
4## the terms of the GNU General Public License as published by the Free Software
5## Foundation; either version 3 of the License, or (at your option) any later
6## version.
7##
8## This program is distributed in the hope that it will be useful, but WITHOUT
9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
11## details.
12##
13## You should have received a copy of the GNU General Public License along with
14## this program; if not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @deftypefn {Function File} {@var{t} =} lsreal (@var{time}, @var{mag}, @var{maxfreq}, @var{numcoeff}, @var{numoctaves})
18##
19## Return a series of least-squares transforms of a real-valued time series.
20## Each transform is minimized independently for each frequency.  The method
21## used is a Lomb-Scargle transform of the real-valued (@var{time}, @var{mag})
22## series, starting from frequency @var{maxfreq} and descending @var{numoctaves}
23## octaves with @var{numcoeff} coefficients per octave.
24##
25## The result of the transform for each frequency is the coefficient of a sum of
26## sine and cosine functions modified by that frequency, in the form of a
27## complex number—where the cosine coefficient is encoded in the real term, and
28## the sine coefficient is encoded in the imaginary term. Each frequency is fit
29## independently from the others, and to minimize very low frequency error,
30## consider storing the mean of a dataset with a constant or near-constant
31## offset separately, and subtracting it from the dataset.
32##
33## @seealso{lscomplex}
34## @end deftypefn
35
36
37
38function transform = lsreal (t, x, omegamax, ncoeff, noctave)
39
40  ## Sanity checks to make sure that the user can get meaningful errors.
41  if (nargin != 5)
42    print_usage ();
43  elseif (! isvector (t))
44    error ("lsreal: Time values are not a vector");
45  elseif (! isvector (x))
46    error ("lsreal: Magnitude values are not a vector");
47  elseif (! all (size (t) == size (x)))
48    error ("lsreal: Size of time vector, magnitude vector unequal");
49  elseif (! isscalar (omegamax))
50    error ("lsreal: More than one value for maximum frequency specified");
51  elseif (! isscalar (ncoeff))
52    error ("lsreal: More than one number of frequencies per octave specified");
53  elseif (! isscalar (noctave))
54    error ("lsreal: More than one number of octaves to traverse specified");
55  elseif (omegamax == 0)
56    error ("lsreal: Specified maximum frequency is not a frequency");
57  elseif (noctave == 0)
58    error ("lsreal: No octaves of results requested");
59  elseif (ncoeff == 0)
60    error ("lsreal: No frequencies per octave requested");
61  elseif (ncoeff != floor (ncoeff))
62    error ("lsreal: Specified number of frequencies per octave is not integral");
63  elseif (noctave != floor (noctave))
64    error ("lsreal: Specified number of octaves of results is not integral");
65  endif
66
67  n = numel (t);
68
69  iter = 0 : (ncoeff * noctave - 1);
70  omul = (2 .^ (- iter / ncoeff));
71
72  ## For a given frequency, the iota term is taken at twice the frequency of the
73  ## zeta term.
74  ot = t(:) * (omul * omegamax);
75  oit = t(:) * (omul * omegamax * 2);
76
77  zeta = sum ((cos (ot) - (sin (ot) .* i)) .* x(:), 1) / n;
78  iota = sum ((cos (oit) - (sin (oit) .* i)), 1) / n;
79
80  transform = 2 .* (conj (zeta) - conj (iota) .* zeta) ./ (1 - abs (iota) .^ 2);
81
82endfunction
83
84%!test
85%! maxfreq = 4 / ( 2 * pi );
86%! t = linspace(0,8);
87%! x = ( 2 .* sin ( maxfreq .* t ) +
88%!       3 .* sin ( (3/4) * maxfreq .* t ) -
89%!       0.5 .* sin ( (1/4) * maxfreq .* t ) -
90%!       0.2 .* cos ( maxfreq .* t ) +
91%!       cos ( (1/4) * maxfreq .* t ) );
92%! # In the assert here, I've got an error bound large enough to catch
93%! # individual system errors which would present no real issue.
94%! assert (lsreal (t,x,maxfreq,2,2),
95%!       [(-1.68275915310663 + 4.70126183846743i), ...
96%!        (1.93821553170889 + 4.95660209883437i), ...
97%!        (4.38145452686697 + 2.14403733658600i), ...
98%!        (5.27425332281147 - 0.73933440226597i)],
99%!         5e-10)
100
101