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