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  {} {} spectral_xdf (@var{x})
28## @deftypefnx {} {} spectral_xdf (@var{x}, @var{win})
29## @deftypefnx {} {} spectral_xdf (@var{x}, @var{win}, @var{b})
30## Return the spectral density estimator given a data vector @var{x}, window
31## name @var{win}, and bandwidth, @var{b}.
32##
33## The window name, e.g., @qcode{"triangle"} or @qcode{"rectangle"} is used to
34## search for a function called @code{@var{win}_sw}.
35##
36## If @var{win} is omitted, the triangle window is used.
37##
38## If @var{b} is omitted, @code{1 / sqrt (length (@var{x}))} is used.
39## @seealso{spectral_adf}
40## @end deftypefn
41
42function retval = spectral_xdf (x, win, b)
43
44  if (nargin < 1 || nargin > 3)
45    print_usage ();
46  endif
47
48  xr = length (x);
49
50  if (columns (x) > 1)
51    x = x';
52  endif
53
54  if (nargin < 3)
55    b = 1 / ceil (sqrt (xr));
56  endif
57
58  if (nargin == 1)
59    w = triangle_sw (xr, b);
60  elseif (! ischar (win))
61    error ("spectral_xdf: WIN must be a string");
62  else
63    win = str2func ([win "_sw"]);
64    w = feval (win, xr, b);
65  endif
66
67  x -= sum (x) / xr;
68
69  retval = (abs (fft (x)) / xr).^2;
70  retval = real (ifft (fft (retval) .* fft (w)));
71
72  retval = [(zeros (xr, 1)), retval];
73  retval(:, 1) = (0 : xr-1)' / xr;
74
75endfunction
76
77
78## Test input validation
79%!error spectral_xdf ()
80%!error spectral_xdf (1, 2, 3, 4)
81%!error spectral_xdf (1, 2)
82%!error spectral_xdf (1, "invalid")
83