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  {} {[@var{Pxx}, @var{w}] =} periodogram (@var{x})
28## @deftypefnx {} {[@var{Pxx}, @var{w}] =} periodogram (@var{x}, @var{win})
29## @deftypefnx {} {[@var{Pxx}, @var{w}] =} periodogram (@var{x}, @var{win}, @var{nfft})
30## @deftypefnx {} {[@var{Pxx}, @var{f}] =} periodogram (@var{x}, @var{win}, @var{nfft}, @var{Fs})
31## @deftypefnx {} {[@var{Pxx}, @var{f}] =} periodogram (@dots{}, "@var{range}")
32## @deftypefnx {} {} periodogram (@dots{})
33## Return the periodogram (Power Spectral Density) of @var{x}.
34##
35## The possible inputs are:
36##
37## @table @var
38## @item x
39##
40## data vector.  If @var{x} is real-valued a one-sided spectrum is estimated.
41## If @var{x} is complex-valued, or @qcode{"@var{range}"} specifies
42## @qcode{"@nospell{twosided}"}, the full spectrum is estimated.
43##
44## @item win
45## window weight data.  If window is empty or unspecified a default rectangular
46## window is used.  Otherwise, the window is applied to the signal
47## (@code{@var{x} .* @var{win}}) before computing the periodogram.  The window
48## data must be a vector of the same length as @var{x}.
49##
50## @item nfft
51## number of frequency bins.  The default is 256 or the next higher power of
52## 2 greater than the length of @var{x}
53## (@code{max (256, 2.^nextpow2 (length (x)))}).  If @var{nfft} is greater
54## than the length of the input then @var{x} will be zero-padded to the length
55## of @var{nfft}.
56##
57## @item Fs
58## sampling rate.  The default is 1.
59##
60## @item range
61## range of spectrum.  @qcode{"@nospell{onesided}"} computes spectrum from
62## [0:nfft/2+1].  @qcode{"@nospell{twosided}"} computes spectrum from
63## [0:nfft-1].
64## @end table
65##
66## The optional second output @var{w} are the normalized angular frequencies.
67## For a one-sided calculation @var{w} is in the range [0, pi] if @var{nfft}
68## is even and [0, pi) if @var{nfft} is odd.  Similarly, for a two-sided
69## calculation @var{w} is in the range [0, 2*pi] or [0, 2*pi) depending on
70## @var{nfft}.
71##
72## If a sampling frequency is specified, @var{Fs}, then the output frequencies
73## @var{f} will be in the range [0, @var{Fs}/2] or [0, @var{Fs}/2) for
74## one-sided calculations.  For two-sided calculations the range will be
75## [0, @var{Fs}).
76##
77## When called with no outputs the periodogram is immediately plotted in the
78## current figure window.
79## @seealso{fft}
80## @end deftypefn
81
82function [pxx, f] = periodogram (x, varargin)
83
84  ## check input arguments
85  if (nargin < 1 || nargin > 5)
86    print_usage ();
87  endif
88
89  nfft = fs = range = window = [];
90  j = 2;
91  for k = 1:length (varargin)
92    if (ischar (varargin{k}))
93      range = varargin{k};
94    else
95      switch (j)
96        case 2
97          window = varargin{k};
98        case 3
99          nfft   = varargin{k};
100        case 4
101          fs     = varargin{k};
102      endswitch
103      j += 1;
104    endif
105  endfor
106
107  if (! isvector (x))
108    error ("periodogram: X must be a real or complex vector");
109  endif
110  x = x(:);  # Use column vectors from now on
111
112  n = rows (x);
113
114  if (! isempty (window))
115    if (! isvector (window) || length (window) != n)
116      error ("periodogram: WIN must be a vector of the same length as X");
117    endif
118    window = window(:);
119    x .*= window;
120  endif
121
122  if (isempty (nfft))
123    nfft = max (256, 2.^nextpow2 (n));
124  elseif (! isscalar (nfft))
125    error ("periodogram: NFFT must be a scalar");
126  endif
127
128  use_w_freq = isempty (fs);
129  if (! use_w_freq && ! isscalar (fs))
130    error ("periodogram: FS must be a scalar");
131  endif
132
133  if (strcmpi (range, "onesided"))
134    range = 1;
135  elseif (strcmpi (range, "twosided"))
136    range = 2;
137  elseif (strcmpi (range, "centered"))
138    error ('periodogram: "centered" range type is not implemented');
139  else
140    range = 2-isreal (x);
141  endif
142
143  ## compute periodogram
144
145  if (n > nfft)
146    Pxx = 0;
147    rr = rem (length (x), nfft);
148    if (rr)
149      x = [x(:); zeros(nfft-rr, 1)];
150    endif
151    x = sum (reshape (x, nfft, []), 2);
152  endif
153
154  if (! isempty (window))
155    n = sumsq (window);
156  endif
157  Pxx = (abs (fft (x, nfft))) .^ 2 / n;
158
159  if (use_w_freq)
160    Pxx /= 2*pi;
161  else
162    Pxx /= fs;
163  endif
164
165  ## generate output arguments
166
167  if (range == 1)  # onesided
168    if (! rem (nfft,2))  # nfft is even
169      psd_len = nfft/2+1;
170      Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1); 0];
171    else                 # nfft is odd
172      psd_len = (nfft+1)/2;
173      Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1)];
174    endif
175  endif
176
177  if (nargout != 1)
178    if (range == 1)
179      f = (0:nfft/2)' / nfft;
180    elseif (range == 2)
181      f = (0:nfft-1)' / nfft;
182    endif
183    if (use_w_freq)
184      f *= 2*pi;  # generate w=2*pi*f
185    else
186      f *= fs;
187    endif
188  endif
189
190  if (nargout == 0)
191    if (use_w_freq)
192      plot (f/(2*pi), 10*log10 (Pxx));
193      xlabel ("normalized frequency [x pi rad]");
194      ylabel ("Power density [dB/rad/sample]");
195    else
196      plot (f, 10*log10 (Pxx));
197      xlabel ("frequency [Hz]");
198      ylabel ("Power density [dB/Hz]");
199    endif
200    grid on;
201    title ("Periodogram Power Spectral Density Estimate");
202  else
203    pxx = Pxx;
204  endif
205
206endfunction
207
208
209## FIXME: Need some functional tests
210
211
212## Test input validation
213%!error periodogram ()
214%!error periodogram (1,2,3,4,5,6)
215%!error <X must be a real or complex vector> periodogram (ones (2,2))
216%!error <WIN must be a vector.*same length> periodogram (1:5, ones (2,2))
217%!error <WIN must be a vector.*same length> periodogram (1:5, 1:6)
218%!error <NFFT must be a scalar> periodogram (1:5, 1:5, 1:5)
219%!error <FS must be a scalar> periodogram (1:5, [], [], 1:5)
220%!error <"centered" range type is not implemented> periodogram (1:5, "centered")
221