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