1########################################################################
2##
3## Copyright (C) 1995-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn {} {} detrend (@var{x}, @var{p})
28## If @var{x} is a vector, @code{detrend (@var{x}, @var{p})} removes the
29## best fit of a polynomial of order @var{p} from the data @var{x}.
30##
31## If @var{x} is a matrix, @code{detrend (@var{x}, @var{p})} does the same
32## for each column in @var{x}.
33##
34## The second argument @var{p} is optional.  If it is not specified, a value of
35## 1 is assumed.  This corresponds to removing a linear trend.
36##
37## The order of the polynomial can also be given as a string, in which case
38## @var{p} must be either @qcode{"constant"} (corresponds to @code{@var{p}=0})
39## or @qcode{"linear"} (corresponds to @code{@var{p}=1}).
40## @seealso{polyfit}
41## @end deftypefn
42
43function y = detrend (x, p = 1)
44
45  if (nargin < 1 || nargin > 2)
46    print_usage ();
47  endif
48
49  if (! isnumeric (x) || ndims (x) > 2)
50    error ("detrend: X must be a numeric vector or matrix");
51  endif
52
53  if (ischar (p) && strcmpi (p, "constant"))
54    p = 0;
55  elseif (ischar (p) && strcmpi (p, "linear"))
56    p = 1;
57  elseif (! isscalar (p) || p < 0 || p != fix (p))
58    error ("detrend: P must be \"constant\", \"linear\", or a positive integer");
59  endif
60
61  [m, n] = size (x);
62  if (m == 1)
63    x = x.';
64  endif
65
66  r = rows (x);
67  b = ((1 : r).' * ones (1, p + 1)) .^ (ones (r, 1) * (0 : p));
68  y = x - b * (b \ x);
69
70  if (m == 1)
71    y = y.';
72  endif
73
74endfunction
75
76
77%!test
78%! N = 32;
79%! x = (0:1:N-1)/N + 2;
80%! y = detrend (x);
81%! assert (abs (y(:)) < 20*eps);
82
83%!test
84%! N = 32;
85%! t = (0:1:N-1)/N;
86%! x = t .* t + 2;
87%! y = detrend (x,2);
88%! assert (abs (y(:)) < 30*eps);
89
90%!test
91%! N = 32;
92%! t = (0:1:N-1)/N;
93%! x = [t;4*t-3].';
94%! y = detrend (x);
95%! assert (abs (y(:)) < 20*eps);
96
97%!test
98%! N = 32;
99%! x = ((0:1:N-1)/N + 2) * 1i;
100%! y = detrend (x);
101%! assert (abs (y(:)) < 20*eps);
102
103## Test input validation
104%!error detrend ()
105%!error detrend (1, 2, 3)
106%!error detrend ("a")
107%!error detrend (true)
108%!error detrend (1, "invalid")
109%!error detrend (1, -1)
110%!error detrend (1, 1.25)
111
112