octave_packages/signal-1.1.3/pwelch.m

   1 %% Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
   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 %% USAGE:
  17 %%   [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
  18 %%                           range,plot_type,detrend,sloppy)
  19 %%     Estimate power spectral density of data "x" by the Welch (1967)
  20 %%     periodogram/FFT method.  All arguments except "x" are optional.
  21 %%         The data is divided into segments.  If "window" is a vector, each
  22 %%     segment has the same length as "window" and is multiplied by "window"
  23 %%     before (optional) zero-padding and calculation of its periodogram. If
  24 %%     "window" is a scalar, each segment has a length of "window" and a
  25 %%     Hamming window is used.
  26 %%         The spectral density is the mean of the periodograms, scaled so that
  27 %%     area under the spectrum is the same as the mean square of the
  28 %%     data.  This equivalence is supposed to be exact, but in practice there
  29 %%     is a mismatch of up to 0.5% when comparing area under a periodogram
  30 %%     with the mean square of the data.
  31 %%
  32 %%  [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs,
  33 %%                          range,plot_type,detrend,sloppy,results)
  34 %%     Two-channel spectrum analyser.  Estimate power spectral density, cross-
  35 %%     spectral density, transfer function and/or coherence functions of time-
  36 %%     series input data "x" and output data "y" by the Welch (1967)
  37 %%     periodogram/FFT method.
  38 %%       pwelch treats the second argument as "y" if there is a control-string
  39 %%     argument "cross", "trans", "coher" or "ypower"; "power" does not force
  40 %%     the 2nd argument to be treated as "y".  All other arguments are
  41 %%     optional.  All spectra are returned in matrix "spectra".
  42 %%
  43 %%  [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
  44 %%                                 range,plot_type,detrend,sloppy)
  45 %%  [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
  46 %%                                 range,plot_type,detrend,sloppy,results)
  47 %%     Estimates confidence intervals for the spectral density.
  48 %%     See Hint (7) below for compatibility options.  Confidence level "conf"
  49 %%     is the 6th or 7th numeric argument.  If "results" control-string
  50 %%     arguments are used, one of them must be "power" when the "conf"
  51 %%     argument is present; pwelch can estimate confidence intervals only for
  52 %%     the power spectrum of the "x" data.  It does not know how to estimate
  53 %%     confidence intervals of the cross-power spectrum, transfer function or
  54 %%     coherence; if you can suggest a good method, please send a bug report.
  55 %%
  56 %% ARGUMENTS
  57 %% All but the first argument are optional and may be empty, except that
  58 %% the "results" argument may require the second argument to be "y".
  59 %%
  60 %% x           %% [non-empty vector] system-input time-series data
  61 %% y           %% [non-empty vector] system-output time-series data
  62 %%
  63 %% window      %% [real vector] of window-function values between 0 and 1; the
  64 %%             %%       data segment has the same length as the window.
  65 %%             %%       Default window shape is Hamming.
  66 %%             %% [integer scalar] length of each data segment.  The default
  67 %%             %%       value is window=sqrt(length(x)) rounded up to the
  68 %%             %%       nearest integer power of 2; see 'sloppy' argument.
  69 %%
  70 %% overlap     %% [real scalar] segment overlap expressed as a multiple of
  71 %%             %%       window or segment length.   0 <= overlap < 1,
  72 %%             %%       The default is overlap=0.5 .
  73 %%
  74 %% Nfft        %% [integer scalar] Length of FFT.  The default is the length
  75 %%             %%       of the "window" vector or has the same value as the
  76 %%             %%       scalar "window" argument.  If Nfft is larger than the
  77 %%             %%       segment length, "seg_len", the data segment is padded
  78 %%             %%       with "Nfft-seg_len" zeros.  The default is no padding.
  79 %%             %%       Nfft values smaller than the length of the data
  80 %%             %%       segment (or window) are ignored silently.
  81 %%
  82 %% Fs          %% [real scalar] sampling frequency (Hertz); default=1.0
  83 %%
  84 %% conf        %% [real scalar] confidence level between 0 and 1.  Confidence
  85 %%             %%       intervals of the spectral density are estimated from
  86 %%             %%       scatter in the periodograms and are returned as Pxx_ci.
  87 %%             %%       Pxx_ci(:,1) is the lower bound of the confidence
  88 %%             %%       interval and Pxx_ci(:,2) is the upper bound.  If there
  89 %%             %%       are three return values, or conf is an empty matrix,
  90 %%             %%       confidence intervals are calculated for conf=0.95 .
  91 %%             %%       If conf is zero or is not given, confidence intervals
  92 %%             %%       are not calculated. Confidence intervals can be
  93 %%             %%       obtained only for the power spectral density of x;
  94 %%             %%       nothing else.
  95 %%
  96 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
  97 %%   Control-string arguments must be after the other arguments but can be in
  98 %%   any order.
  99 %%
 100 %% range     %% 'half',  'onesided' : frequency range of the spectrum is
 101 %%           %%       zero up to but not including Fs/2.  Power from
 102 %%           %%       negative frequencies is added to the positive side of
 103 %%           %%       the spectrum, but not at zero or Nyquist (Fs/2)
 104 %%           %%       frequencies.  This keeps power equal in time and
 105 %%           %%       spectral domains.  See reference [2].
 106 %%           %% 'whole', 'twosided' : frequency range of the spectrum is
 107 %%           %%       -Fs/2 to Fs/2, with negative frequencies
 108 %%           %%       stored in "wrap around" order after the positive
 109 %%           %%       frequencies; e.g. frequencies for a 10-point 'twosided'
 110 %%           %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
 111 %%           %% 'shift', 'centerdc' : same as 'whole' but with the first half
 112 %%           %%       of the spectrum swapped with second half to put the
 113 %%           %%       zero-frequency value in the middle. (See "help
 114 %%           %%       fftshift".
 115 %%           %% If data (x and y) are real, the default range is 'half',
 116 %%           %% otherwise default range is 'whole'.
 117 %%
 118 %% plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
 119 %%           %% specifies the type of plot.  The default is 'plot', which
 120 %%           %% means linear-linear axes. 'squared' is the same as 'plot'.
 121 %%           %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
 122 %%           %% spectrum is not plotted if the caller requires a returned
 123 %%           %% value.
 124 %%
 125 %% detrend   %% 'no-strip', 'none' -- do NOT remove mean value from the data
 126 %%           %% 'short', 'mean' -- remove the mean value of each segment from
 127 %%           %%                    each segment of the data.
 128 %%           %% 'linear',       -- remove linear trend from each segment of
 129 %%           %%                    the data.
 130 %%           %% 'long-mean'     -- remove the mean value from the data before
 131 %%           %%              splitting it into segments.  This is the default.
 132 %%
 133 %%   sloppy  %% 'sloppy': FFT length is rounded up to the nearest integer
 134 %%           %%       power of 2 by zero padding.  FFT length is adjusted
 135 %%           %%       after addition of padding by explicit Nfft argument.
 136 %%           %%       The default is to use exactly the FFT and window/
 137
 138 %%           %%       segment lengths specified in argument list.
 139 %%
 140 %%   results %% specifies what results to return (in the order specified
 141 %%           %%   and as many as desired).
 142 %%           %% 'power' calculate power spectral density of "x"
 143 %%           %% 'cross' calculate cross spectral density of "x" and "y"
 144 %%           %% 'trans' calculate transfer function of a system with
 145 %%           %%         input "x" and output "y"
 146 %%           %% 'coher' calculate coherence function of "x" and "y"
 147 %%           %% 'ypower' calculate power spectral density of "y"
 148 %%           %%  The default is 'power', with argument "y" omitted.
 149 %%
 150 %% RETURNED VALUES:
 151 %%   If return values are not required by the caller, the results are
 152 %%     plotted and nothing is returned.
 153 %%
 154 %%   spectra %% [real-or-complex matrix] columns of the matrix contain results
 155 %%           %%        in the same order as specified by "results" arguments.
 156 %%           %%        Each column contains one of the result vectors.
 157 %%
 158 %%   Pxx_ci  %% [real matrix] estimate of confidence interval for power
 159 %%           %%        spectral density of x.  First column is the lower
 160 %%           %%        bound.  Second column is the upper bound.
 161 %%
 162 %%   freq    %% [real column vector] frequency values
 163 %%
 164 %% HINTS
 165 %% 1) EMPTY ARGS:
 166 %%    if you don't want to use an optional argument you can leave it empty
 167 %%    by writing its value as [].
 168 %% 2) FOR BEGINNERS:
 169 %%    The profusion of arguments may make pwelch difficult to use, and an
 170 %%    unskilled user can easily produce a meaningless result or can easily
 171 %%    mis-interpret the result.  With real data "x" and sampling frequency
 172 %%    "Fs", the easiest and best way for a beginner to use pwelch is
 173 %%    probably "pwelch(x,[],[],[],Fs)".  Use the "window" argument to
 174 %%    control the length of the spectrum vector.  For real data and integer
 175 %%    scalar M, "pwelch(x,2*M,[],[],Fs)" gives an M+1 point spectrum.
 176 %%    Run "demo pwelch" (octave only).
 177 %% 3) WINDOWING FUNCTIONS:
 178 %%    Without a window function, sharp spectral peaks can have strong
 179 %%    sidelobes because the FFT of a data in a segment is in effect convolved
 180 %%    with a rectangular window.  A window function which tapers off
 181 %%    (gradually) at the ends produces much weaker sidelobes in the FFT.
 182 %%    Hann (hanning), hamming, bartlett, blackman, flattopwin etc are
 183 %%    available as separate Matlab/sigproc or Octave functions.  The sidelobes
 184 %%    of the Hann window have a roll-off rate of 60dB/decade of frequency.
 185 %%    The first sidelobe of the Hamming window is suppressed and is about 12dB
 186 %%    lower than the first Hann sidelobe, but the roll-off rate is only
 187 %%    20dB/decade.  You can inspect the FFT of a Hann window by plotting
 188 %%    "abs(fft(postpad(hanning(256),4096,0)))".
 189 %%    The default window is Hamming.
 190 %% 4) ZERO PADDING:
 191 %%    Zero-padding reduces the frequency step in the
 192 %%    spectrum, and produces an artificially smoothed spectrum.  For example,
 193 %%    "Nfft=2*length(window)" gives twice as many frequency values, but
 194 %%    adjacent PSD (power spectral density) values are not independent;
 195 %%    adjacent PSD values are independent if "Nfft=length(window)", which is
 196 %%    the default value of Nfft.
 197 %% 5) REMOVING MEAN FROM SIGNAL:
 198 %%    If the mean is not removed from the signal there is a large spectral
 199 %%    peak at zero frequency and the sidelobes of this peak are likely to
 200 %%    swamp the rest of the spectrum.  For this reason, the default behaviour
 201 %%    is to remove the mean.  However, the matlab pwelch does not do this.
 202 %% 6) WARNING ON CONFIDENCE INTERVALS
 203 %%    Confidence intervals are obtained by measuring the sample variance of
 204 %%    the periodograms and assuming that the periodograms have a Gaussian
 205 %%    probability distribution.  This assumption is not accurate.  If, for
 206 %%    example, the data (x) is Gaussian, the periodogram has a Rayleigh
 207 %%    distribution.  However, the confidence intervals may still be  useful.
 208 %%
 209 %% 7) COMPATIBILITY WITH Matlab R11, R12, etc
 210 %%    When used without the second data (y) argument, arguments are compatible
 211 %%    with the pwelch of Matlab R12, R13, R14, 2006a and 2006b except that
 212 %%     1) overlap is expressed as a multiple of window length ---
 213 %%        effect of overlap scales with window length
 214 %%     2) default values of length(window), Nfft and Fs are more sensible, and
 215 %%     3) Goertzel algorithm is not available so Nfft cannot be an array of
 216 %%        frequencies as in Matlab 2006b.
 217 %%    Pwelch has four persistent Matlab-compatibility levels.  Calling pwelch
 218 %%    with an empty first argument sets the order of arguments and defaults
 219 %%    specified above in the USAGE and ARGUMENTS section of this documentation.
 220 %%          prev_compat=pwelch([]);
 221 %%          [Pxx,f]=pwelch(x,window,overlap,Nfft,Fs,conf,...);
 222 %%    Calling pwelch with a single string argument (as described below) gives
 223 %%    compatibility with Matlab R11 or R12, or the R14 spectrum.welch
 224 %%    defaults.  The returned value is the PREVIOUS compatibility string.
 225 %%
 226 %%    Matlab R11:  For compatibility with the Matlab R11 pwelch:
 227 %%          prev_compat=pwelch('R11-');
 228 %%          [Pxx,f]=pwelch(x,Nfft,Fs,window,overlap,conf,range,units);
 229 %%          %% units of overlap are "number of samples"
 230 %%          %% defaults: Nfft=min(length(x),256), Fs=2*pi, length(window)=Nfft,
 231 %%          %%           window=Hanning, do not detrend,
 232 %%          %% N.B.  "Sloppy" is not available.
 233 %%
 234 %%    Matlab R12:  For compatibility with Matlab R12 to 2006a pwelch:
 235 %%          prev_compat=pwelch('R12+');
 236 %%          [Pxx,f]=pwelch(x,window,overlap,nfft,Fs,...);
 237 %%          %% units of overlap are "number of samples"
 238 %%          %% defaults: length(window)==length(x)/8, window=Hamming,
 239 %%          %%           Nfft=max(256,NextPow2), Fs=2*pi, do not detrend
 240 %%          %% NextPow2 is the next power of 2 greater than or equal to the
 241 %%          %% window length. "Sloppy", "conf" are not available.  Default
 242 %%          %% window length gives very poor amplitude resolution.
 243 %%
 244 %%    To adopt defaults of the Matlab R14 "spectrum.welch" spectrum object
 245 %%    associated "psd" method.
 246 %%          prev_compat=pwelch('psd');
 247 %%          [Pxx,f] = pwelch(x,window,overlap,Nfft,Fs,conf,...);
 248 %%          %% overlap is expressed as a percentage of window length,
 249 %%          %% defaults: length(window)==64, Nfft=max(256,NextPow2), Fs=2*pi
 250 %%          %%           do not detrend
 251 %%          %% NextPow2 is the next power of 2 greater than or equal to the
 252 %%          %% window length. "Sloppy" is not available.
 253 %%          %% Default window length gives coarse frequency resolution.
 254 %%
 255 %%
 256 %% REFERENCES
 257 %%  [1] Peter D. Welch (June 1967):
 258 %%   "The use of fast Fourier transform for the estimation of power spectra:
 259 %%   a method based on time averaging over short, modified periodograms."
 260 %%   IEEE Transactions on Audio Electroacoustics, Vol AU-15(6), pp 70-73
 261 %%
 262 %%  [2] William H. Press and Saul A. Teukolsky and William T. Vetterling and
 263 %%               Brian P. Flannery",
 264 %%   "Numerical recipes in C, The art of scientific computing", 2nd edition,
 265 %%      Cambridge University Press, 2002 --- Section 13.7.
 266 %%  [3] Paul Kienzle (1999-2001): "pwelch", http://octave.sourceforge.net/
 267
 268 function [varargout] = pwelch(x,varargin)
 269   %%
 270   %% COMPATIBILITY LEVEL
 271   %% Argument positions and defaults depend on compatibility level selected
 272   %% by calling pwelch without arguments or with a single string argument.
 273   %%   native:      compatib=1; prev_compat=pwelch(); prev_compat=pwelch([]);
 274   %%   matlab R11:  compatib=2; prev_compat=pwelch('R11-');
 275   %%   matlab R12:  compatib=3; prev_compat=pwelch('R12+');
 276   %%   spectrum.welch defaults:  compatib=4; prev_compat=pwelch('psd');
 277   %% In each case, the returned value is the PREVIOUS compatibility string.
 278   %%
 279   compat_str = {[]; 'R11-'; 'R12+'; 'psd'};
 280   persistent compatib;
 281   if ( isempty(compatib) || compatib<=0 || compatib>4 )
 282     %% legal values are 1, 2, 3, 4
 283
 284     compatib = 1;
 285   end
 286   if ( nargin <= 0 )
 287     error( 'pwelch: Need at least 1 arg. Use "help pwelch".' );
 288   elseif ( nargin==1 && (ischar(x) || isempty(x)) )
 289     varargout{1} = compat_str{compatib};
 290     if ( isempty(x) ) % native
 291       compatib = 1;
 292     elseif ( strcmp(x,'R11-') )
 293       compatib = 2;
 294     elseif ( strcmp(x,'R12+') )
 295       compatib = 3;
 296     elseif ( strcmp(x,'psd') )
 297       compatib = 4;
 298     else
 299       error( 'pwelch: compatibility arg must be empty, R11-, R12+ or psd' );
 300     end
 301     %% return
 302   %%
 303   %% Check fixed argument
 304   elseif ( isempty(x) || ~isvector(x) )
 305     error( 'pwelch: arg 1 (x) must be vector.' );
 306   else
 307     %%  force x to be COLUMN vector
 308     if ( size(x,1)==1 )
 309       x=x(:);
 310     end
 311     %%
 312     %% Look through all args to check if  cross PSD, transfer function or
 313     %% coherence is required.  If yes, the second arg is data vector "y".
 314     arg2_is_y = 0;
 315     x_len = length(x);
 316     nvarargin = length(varargin);
 317     for iarg=1:nvarargin
 318       arg = varargin{iarg};
 319       if ( ~isempty(arg) && ischar(arg) && ...
 320            ( strcmp(arg,'cross') || strcmp(arg,'trans') || ...
 321              strcmp(arg,'coher') || strcmp(arg,'ypower') ))
 322         %% OK. Need "y". Grab it from 2nd arg.
 323         arg = varargin{1};
 324         if ( nargin<2 || isempty(arg) || ~isvector(arg) || length(arg)~=x_len )
 325           error( 'pwelch: arg 2 (y) must be vector, same length as x.' );
 326         end
 327         %% force  COLUMN vector
 328         y = varargin{1}(:);
 329         arg2_is_y = 1;
 330         break;
 331       end
 332     end
 333     %%
 334     %% COMPATIBILITY
 335     %% To select default argument values, "compatib" is used as an array index.
 336     %% Index values are   1=native,  2=R11,  3=R12,  4=spectrum.welch
 337     %%
 338     %%  argument positions:
 339     %%  arg_posn = varargin index of window, overlap, Nfft, Fs and conf
 340     %%             args respectively, a value of zero ==>> arg does not exist
 341     arg_posn = [1 2 3 4 5;  %% native
 342                 3 4 1 2 5;  %% Matlab R11- pwelch
 343                 1 2 3 4 0;  %% Matlab R12+ pwelch
 344                 1 2 3 4 5]; %% spectrum.welch defaults
 345     arg_posn  = arg_posn(compatib,:) + arg2_is_y;
 346     %%
 347     %%  SPECIFY SOME DEFAULT VALUES for (not all) optional arguments
 348     %%  Use compatib as array index.
 349     %%  Fs = sampling frequency
 350     Fs        = [ 1.0 2*pi 2*pi 2*pi ];
 351     Fs        = Fs(compatib);
 352     %%  plot_type: 1='plot'|'squared'; 5='db'|'dB'
 353     plot_type = [ 1 5 5 5 ];
 354     plot_type = plot_type(compatib);
 355     %%  rm_mean: 3='long-mean'; 0='no-strip'|'none'
 356     rm_mean   = [ 3 0 0 0 ];
 357     rm_mean   = rm_mean(compatib);
 358     %% use max_overlap=x_len-1 because seg_len is not available yet
 359     %% units of overlap are different for each version:
 360     %%    fraction, samples, or percent
 361     max_overlap = [ 0.95 x_len-1 x_len-1 95];
 362     max_overlap = max_overlap(compatib);
 363     %% default confidence interval
 364     %%  if there are more than 2 return values and if there is a "conf" arg
 365     conf      = 0.95 * (nargout>2) * (arg_posn(5)>0);
 366     %%
 367     is_win    = 0;    % =0 means valid window arg is not provided yet
 368     Nfft      = [];   % default depends on segment length
 369     overlap   = [];   % WARNING: units can be #samples, fraction or percentage
 370     range     = ~isreal(x) || ( arg2_is_y && ~isreal(y) );
 371     is_sloppy = 0;
 372     n_results = 0;
 373     do_power  = 0;
 374     do_cross  = 0;
 375     do_trans  = 0;
 376     do_coher  = 0;
 377     do_ypower = 0;
 378   %%
 379   %%  DECODE AND CHECK OPTIONAL ARGUMENTS
 380     end_numeric_args = 0;
 381     for iarg = 1+arg2_is_y:nvarargin
 382       arg = varargin{iarg};
 383       if ( ischar(arg) )
 384         %% first string arg ==> no more numeric args
 385         %% non-string args cannot follow a string arg
 386         end_numeric_args = 1;
 387         %%
 388         %% decode control-string arguments
 389         if ( strcmp(arg,'sloppy') )
 390           is_sloppy = ~is_win || is_win==1;
 391         elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
 392           plot_type = 1;
 393         elseif ( strcmp(arg,'semilogx') )
 394           plot_type = 2;
 395         elseif ( strcmp(arg,'semilogy') )
 396           plot_type = 3;
 397         elseif ( strcmp(arg,'loglog') )
 398           plot_type = 4;
 399         elseif ( strcmp(arg,'db') || strcmp(arg,'dB') )
 400           plot_type = 5;
 401         elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
 402           range = 0;
 403         elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
 404           range = 1;
 405         elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
 406           range = 2;
 407         elseif ( strcmp(arg,'long-mean') )
 408           rm_mean = 3;
 409         elseif ( strcmp(arg,'linear') )
 410           rm_mean = 2;
 411         elseif ( strcmp(arg,'short') || strcmp(arg,'mean') )
 412           rm_mean = 1;
 413         elseif ( strcmp(arg,'no-strip') || strcmp(arg,'none') )
 414           rm_mean = 0;
 415         elseif ( strcmp(arg, 'power' ) )
 416           if ( ~do_power )
 417             n_results = n_results+1;
 418             do_power = n_results;
 419           end
 420         elseif ( strcmp(arg, 'cross' ) )
 421           if ( ~do_cross )
 422             n_results = n_results+1;
 423             do_cross = n_results;
 424           end
 425         elseif ( strcmp(arg, 'trans' ) )
 426           if ( ~do_trans )
 427             n_results = n_results+1;
 428             do_trans = n_results;
 429           end
 430         elseif ( strcmp(arg, 'coher' ) )
 431           if ( ~do_coher )
 432             n_results = n_results+1;
 433             do_coher = n_results;
 434           end
 435         elseif ( strcmp(arg, 'ypower' ) )
 436           if ( ~do_ypower )
 437             n_results = n_results+1;
 438             do_ypower = n_results;
 439           end
 440         else
 441           error( 'pwelch: string arg %d illegal value: %s', iarg+1, arg );
 442         end
 443         %% end of processing string args
 444         %%
 445       elseif ( end_numeric_args )
 446         if ( ~isempty(arg) )
 447           %% found non-string arg after a string arg ... oops
 448           error( 'pwelch: control arg must be string' );
 449         end
 450       %%
 451       %% first 4 optional arguments are numeric -- in fixed order
 452       %%
 453       %% deal with "Fs" and "conf" first because empty arg is a special default
 454       %% -- "Fs" arg -- sampling frequency
 455       elseif ( iarg == arg_posn(4) )
 456         if ( isempty(arg) )
 457           Fs = 1;
 458         elseif ( ~isscalar(arg) || ~isreal(arg) || arg<0 )
 459           error( 'pwelch: arg %d (Fs) must be real scalar >0', iarg+1 );
 460         else
 461           Fs = arg;
 462         end
 463       %%
 464       %%  -- "conf" arg -- confidence level
 465       %%    guard against the "it cannot happen" iarg==0
 466       elseif ( arg_posn(5) && iarg == arg_posn(5) )
 467         if ( isempty(arg) )
 468           conf = 0.95;
 469         elseif ( ~isscalar(arg) || ~isreal(arg) || arg < 0.0 || arg >= 1.0 )
 470           error( 'pwelch: arg %d (conf) must be real scalar, >=0, <1',iarg+1 );
 471         else
 472           conf = arg;
 473         end
 474       %%
 475       %% skip all empty args from this point onward
 476       elseif ( isempty(arg) )
 477         1;
 478       %%
 479       %%  -- "window" arg -- window function
 480       elseif ( iarg == arg_posn(1) )
 481         if ( isscalar(arg) )
 482           is_win = 1;
 483         elseif ( isvector(arg) )
 484           is_win = length(arg);
 485           if ( size(arg,2)>1 )  %% vector must be COLUMN vector
 486             arg = arg(:);
 487           end
 488         else
 489           is_win = 0;
 490         end
 491         if ( ~is_win )
 492           error( 'pwelch: arg %d (window) must be scalar or vector', iarg+1 );
 493         elseif ( is_win==1 && ( ~isreal(arg) || fix(arg)~=arg || arg<=3 ) )
 494           error( 'pwelch: arg %d (window) must be integer >3', iarg+1 );
 495         elseif ( is_win>1 && ( ~isreal(arg) || any(arg<0) ) )
 496           error( 'pwelch: arg %d (window) vector must be real and >=0',iarg+1);
 497         end
 498         window = arg;
 499         is_sloppy = 0;
 500       %%
 501       %% -- "overlap" arg -- segment overlap
 502       elseif ( iarg == arg_posn(2) )
 503         if (~isscalar(arg) || ~isreal(arg) || arg<0 || arg>max_overlap )
 504           error( 'pwelch: arg %d (overlap) must be real from 0 to %f', ...
 505                  iarg+1, max_overlap );
 506         end
 507         overlap = arg;
 508       %%
 509       %% -- "Nfft" arg -- FFT length
 510       elseif ( iarg == arg_posn(3) )
 511         if ( ~isscalar(arg) || ~isreal(arg) || fix(arg)~=arg || arg<0 )
 512           error( 'pwelch: arg %d (Nfft) must be integer >=0', iarg+1 );
 513         end
 514         Nfft = arg;
 515       %%
 516       else
 517         error( 'pwelch: arg %d  must be string', iarg+1 );
 518       end
 519     end
 520     if ( conf>0 && (n_results && ~do_power ) )
 521      error('pwelch: can give confidence interval for x power spectrum only' );
 522       end
 523     %%
 524     %% end DECODE AND CHECK OPTIONAL ARGUMENTS.
 525     %%
 526     %% SETUP REMAINING PARAMETERS
 527     %% default action is to calculate power spectrum only
 528     if ( ~n_results )
 529       n_results = 1;
 530       do_power = 1;
 531     end
 532     need_Pxx = do_power || do_trans || do_coher;
 533     need_Pxy = do_cross || do_trans || do_coher;
 534     need_Pyy = do_coher || do_ypower;
 535     log_two = log(2);
 536     nearly_one = 0.99999999999;
 537     %%
 538     %% compatibility-options
 539     %% provides exact compatibility with Matlab R11 or R12
 540     %%
 541     %% Matlab R11 compatibility
 542     if ( compatib==2 )
 543       if ( isempty(Nfft) )
 544         Nfft = min( 256, x_len );
 545       end
 546       if ( is_win > 1 )
 547         seg_len = min( length(window), Nfft );
 548         window = window(1:seg_len);
 549       else
 550         if ( is_win )
 551           %% window arg is scalar
 552           seg_len = window;
 553         else
 554           seg_len = Nfft;
 555         end
 556         %% make Hann window (don't depend on sigproc)
 557         xx = seg_len - 1;
 558         window = 0.5 - 0.5 * cos( (2*pi/xx)*[0:xx].' );
 559       end
 560     %%
 561     %% Matlab R12 compatibility
 562     elseif ( compatib==3 )
 563       if ( is_win > 1 )
 564         %% window arg provides window function
 565         seg_len = length(window);
 566       else
 567         %% window arg does not provide window function; use Hamming
 568         if ( is_win )
 569           %% window arg is scalar
 570           seg_len = window;
 571         else
 572           %% window arg not available; use R12 default, 8 windows
 573           %% ignore overlap arg; use overlap=50% -- only choice that makes sense
 574           %% this is the magic formula for 8 segments with 50% overlap
 575           seg_len = fix( (x_len-3)*2/9 );
 576         end
 577         %% make Hamming window (don't depend on sigproc)
 578         xx = seg_len - 1;
 579         window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
 580         end
 581       if ( isempty(Nfft) )
 582         Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
 583       end
 584     %%
 585     %% Matlab R14 psd(spectrum.welch) defaults
 586     elseif ( compatib==4 )
 587       if ( is_win > 1 )
 588         %% window arg provides window function
 589         seg_len = length(window);
 590       else
 591         %% window arg does not provide window function; use Hamming
 592         if ( is_win )
 593           %% window arg is scalar
 594           seg_len = window;
 595         else
 596           %% window arg not available; use default seg_len = 64
 597           seg_len = 64;
 598         end
 599         %% make Hamming window (don't depend on sigproc)
 600         xx = seg_len - 1;
 601         window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
 602         end
 603       %% Now we know segment length,
 604       %% so we can set default overlap as number of samples
 605       if ( ~isempty(overlap) )
 606         overlap = fix(seg_len * overlap / 100 );
 607         end
 608       if ( isempty(Nfft) )
 609         Nfft = max( 256, 2^ceil(log(seg_len)*nearly_one/log_two) );
 610       end
 611     %%
 612     %% default compatibility level
 613     else %if ( compatib==1 )
 614       %% calculate/adjust segment length, window function
 615       if ( is_win > 1 )
 616         %% window arg provides window function
 617         seg_len = length(window);
 618       else
 619         %% window arg does not provide window function; use Hamming
 620         if ( is_win )       %% window arg is scalar
 621           seg_len = window;
 622         else
 623           %% window arg not available; use default length:
 624           %% = sqrt(length(x)) rounded up to nearest integer power of 2
 625           if ( isempty(overlap) )
 626             overlap=0.5;
 627           end
 628           seg_len = 2 ^ ceil( log(sqrt(x_len/(1-overlap)))*nearly_one/log_two );
 629         end
 630         %% make Hamming window (don't depend on sigproc)
 631         xx = seg_len - 1;
 632         window = 0.54 - 0.46 * cos( (2*pi/xx)*[0:xx].' );
 633         end
 634       %% Now we know segment length,
 635       %% so we can set default overlap as number of samples
 636       if ( ~isempty(overlap) )
 637         overlap = fix(seg_len * overlap);
 638         end
 639       %%
 640       %% calculate FFT length
 641       if ( isempty(Nfft) )
 642         Nfft = seg_len;
 643       end
 644       if ( is_sloppy )
 645         Nfft = 2 ^ ceil( log(Nfft) * nearly_one / log_two );
 646       end
 647     end
 648     %% end of compatibility options
 649     %%
 650     %% minimum FFT length is seg_len
 651     Nfft = max( Nfft, seg_len );
 652     %% Mean square of window is required for normalising PSD amplitude.
 653     win_meansq = (window.' * window) / seg_len;
 654     %%
 655     %% Set default or check overlap.
 656     if ( isempty(overlap) )
 657       overlap = fix(seg_len /2);
 658     elseif ( overlap >= seg_len )
 659       error( 'pwelch: arg (overlap=%d) too big. Must be <length(window)=%d',...
 660              overlap, seg_len );
 661     end
 662     %%
 663     %% Pad data with zeros if shorter than segment. This should not happen.
 664     if ( x_len < seg_len )
 665       x = [x; zeros(seg_len-x_len,1)];
 666       if ( arg2_is_y )
 667         y = [y; zeros(seg_len-x_len,1)];
 668         end
 669       x_len = seg_len;
 670       end
 671     %% end SETUP REMAINING PARAMETERS
 672     %%
 673     %%
 674     %% MAIN CALCULATIONS
 675     %% Remove mean from the data
 676     if ( rm_mean == 3 )
 677       n_ffts = max( 0, fix( (x_len-seg_len)/(seg_len-overlap) ) ) + 1;
 678       x_len  = min( x_len, (seg_len-overlap)*(n_ffts-1)+seg_len );
 679       if ( need_Pxx || need_Pxy )
 680         x = x - sum( x(1:x_len) ) / x_len;
 681       end
 682       if ( arg2_is_y || need_Pxy)
 683         y = y - sum( y(1:x_len) ) / x_len;
 684       end
 685     end
 686     %%
 687     %% Calculate and accumulate periodograms
 688     %%   xx and yy are padded data segments
 689     %%   Pxx, Pyy, Pyy are periodogram sums, Vxx is for confidence interval
 690     xx = zeros(Nfft,1);
 691     yy = xx;
 692     Pxx = xx;
 693     Pxy = xx;
 694     Pyy = xx;
 695     if ( conf>0 )
 696       Vxx = xx;
 697     else
 698       Vxx = [];
 699     end
 700     n_ffts = 0;
 701     for start_seg = [1:seg_len-overlap:x_len-seg_len+1]
 702       end_seg = start_seg+seg_len-1;
 703       %% Don't truncate/remove the zero padding in xx and yy
 704       if ( need_Pxx || need_Pxy )
 705         if ( rm_mean==1 ) % remove mean from segment
 706           xx(1:seg_len) = window .* ( ...
 707             x(start_seg:end_seg) - sum(x(start_seg:end_seg)) / seg_len);
 708         elseif ( rm_mean == 2 ) % remove linear trend from segment
 709           xx(1:seg_len) = window .* detrend( x(start_seg:end_seg) );
 710         else % rm_mean==0 or 3
 711           xx(1:seg_len) = window .* x(start_seg:end_seg);
 712         end
 713         fft_x = fft(xx);
 714       end
 715       if ( need_Pxy || need_Pyy )
 716         if ( rm_mean==1 ) % remove mean from segment
 717           yy(1:seg_len) = window .* ( ...
 718             y(start_seg:end_seg) - sum(y(start_seg:end_seg)) / seg_len);
 719         elseif ( rm_mean == 2 ) % remove linear trend from segment
 720           yy(1:seg_len) = window .* detrend( y(start_seg:end_seg) );
 721         else % rm_mean==0 or 3
 722           yy(1:seg_len) = window .* y(start_seg:end_seg);
 723         end
 724         fft_y = fft(yy);
 725       end
 726       if ( need_Pxx )
 727         %% force Pxx to be real; pgram = periodogram
 728         pgram = real(fft_x .* conj(fft_x));
 729         Pxx = Pxx + pgram;
 730         %% sum of squared periodograms is required for confidence interval
 731         if ( conf>0 )
 732           Vxx = Vxx + pgram .^2;
 733           end
 734       end
 735       if ( need_Pxy )
 736         %% Pxy (cross power spectrum) is complex. Do not force to be real.
 737         Pxy = Pxy + fft_y .* conj(fft_x);
 738       end
 739       if ( need_Pyy )
 740         %% force Pyy to be real
 741         Pyy = Pyy + real(fft_y .* conj(fft_y));
 742       end
 743       n_ffts = n_ffts +1;
 744     end
 745     %%
 746     %% Calculate confidence interval
 747     %%    -- incorrectly assumes that the periodogram has Gaussian probability
 748     %%       distribution (actually, it has a single-sided (e.g. exponential)
 749     %%       distribution.
 750     %% Sample variance of periodograms is (Vxx-Pxx.^2/n_ffts)/(n_ffts-1).
 751     %%    This method of calculating variance is more susceptible to round-off
 752     %%  error, but is quicker, and for double-precision arithmetic and the
 753     %%  inherently noisy periodogram (variance==mean^2), it should be OK.
 754     if ( conf>0 && need_Pxx )
 755       if ( n_ffts<2 )
 756         Vxx = zeros(Nfft,1);
 757       else
 758         %% Should use student distribution here (for unknown variance), but tinv
 759         %% is not a core Matlab function (is in statistics toolbox. Grrr)
 760         Vxx = (erfinv(conf)*sqrt(2*n_ffts/(n_ffts-1))) * sqrt(Vxx-Pxx.^2/n_ffts);
 761       end
 762     end
 763     %%
 764     %% Convert two-sided spectra to one-sided spectra (if range == 0).
 765     %% For one-sided spectra, contributions from negative frequencies are added
 766     %% to the positive side of the spectrum -- but not at zero or Nyquist
 767     %% (half sampling) frequencies.  This keeps power equal in time and spectral
 768     %% domains, as required by Parseval theorem.
 769     %%
 770     if ( range == 0 )
 771       if ( ~ rem(Nfft,2) )    %% one-sided, Nfft is even
 772         psd_len = Nfft/2+1;
 773         if ( need_Pxx )
 774           Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1); 0];
 775           if ( conf>0 )
 776             Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1); 0];
 777           end
 778         end
 779         if ( need_Pxy )
 780           Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1); 0]);
 781         end
 782         if ( need_Pyy )
 783           Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1); 0];
 784         end
 785       else                    %% one-sided, Nfft is odd
 786         psd_len = (Nfft+1)/2;
 787         if ( need_Pxx )
 788           Pxx = Pxx(1:psd_len) + [0; Pxx(Nfft:-1:psd_len+1)];
 789           if ( conf>0 )
 790             Vxx = Vxx(1:psd_len) + [0; Vxx(Nfft:-1:psd_len+1)];
 791           end
 792         end
 793         if ( need_Pxy )
 794           Pxy = Pxy(1:psd_len) + conj([0; Pxy(Nfft:-1:psd_len+1)]);
 795         end
 796         if ( need_Pyy )
 797           Pyy = Pyy(1:psd_len) + [0; Pyy(Nfft:-1:psd_len+1)];
 798         end
 799       end
 800     else                      %% two-sided (and shifted)
 801       psd_len = Nfft;
 802     end
 803     %% end MAIN CALCULATIONS
 804     %%
 805     %% SCALING AND OUTPUT
 806     %% Put all results in matrix, one row per spectrum
 807     %%   Pxx, Pxy, Pyy are sums of periodograms, so "n_ffts"
 808     %%   in the scale factor converts them into averages
 809     spectra    = zeros(psd_len,n_results);
 810     spect_type = zeros(n_results,1);
 811     scale = n_ffts * seg_len * Fs * win_meansq;
 812     if ( do_power )
 813       spectra(:,do_power) = Pxx / scale;
 814       spect_type(do_power) = 1;
 815       if ( conf>0 )
 816         Vxx = [Pxx-Vxx Pxx+Vxx]/scale;
 817       end
 818     end
 819     if ( do_cross )
 820       spectra(:,do_cross) = Pxy / scale;
 821       spect_type(do_cross) = 2;
 822     end
 823     if ( do_trans )
 824       spectra(:,do_trans) = Pxy ./ Pxx;
 825       spect_type(do_trans) = 3;
 826     end
 827     if ( do_coher )
 828       %% force coherence to be real
 829       spectra(:,do_coher) = real(Pxy .* conj(Pxy)) ./ Pxx ./ Pyy;
 830       spect_type(do_coher) = 4;
 831     end
 832     if ( do_ypower )
 833       spectra(:,do_ypower) = Pyy / scale;
 834       spect_type(do_ypower) = 5;
 835     end
 836     freq = [0:psd_len-1].' * ( Fs / Nfft );
 837     %%
 838     %% range='shift': Shift zero-frequency to the middle
 839     if ( range == 2 )
 840       len2 = fix((Nfft+1)/2);
 841       spectra = [ spectra(len2+1:Nfft,:); spectra(1:len2,:)];
 842       freq    = [ freq(len2+1:Nfft)-Fs; freq(1:len2)];
 843       if ( conf>0 )
 844         Vxx = [ Vxx(len2+1:Nfft,:); Vxx(1:len2,:)];
 845       end
 846     end
 847     %%
 848     %%  RETURN RESULTS or PLOT
 849     if ( nargout>=2 && conf>0 )
 850       varargout{2} = Vxx;
 851     end
 852     if ( nargout>=(2+(conf>0)) )
 853       %% frequency is 2nd or 3rd return value,
 854       %% depends on if 2nd is confidence interval
 855       varargout{2+(conf>0)} = freq;
 856     end
 857     if ( nargout>=1 )
 858       varargout{1} = spectra;
 859     else
 860       %%
 861       %% Plot the spectra if there are no return variables.
 862       plot_title=['power spectrum x ';
 863                   'cross spectrum   ';
 864                   'transfer function';
 865                   'coherence        ';
 866                   'power spectrum y ' ];
 867       for ii = 1: n_results
 868         if ( conf>0 && spect_type(ii)==1 )
 869           Vxxxx = Vxx;
 870         else
 871           Vxxxx = [];
 872         end
 873         if ( n_results > 1 )
 874           figure();
 875           end
 876         if ( plot_type == 1 )
 877           plot(freq,[abs(spectra(:,ii)) Vxxxx]);
 878         elseif ( plot_type == 2 )
 879           semilogx(freq,[abs(spectra(:,ii)) Vxxxx]);
 880         elseif ( plot_type == 3 )
 881           semilogy(freq,[abs(spectra(:,ii)) Vxxxx]);
 882         elseif ( plot_type == 4 )
 883           loglog(freq,[abs(spectra(:,ii)) Vxxxx]);
 884         elseif ( plot_type == 5 )  % db
 885           ylabel( 'amplitude (dB)' );
 886           plot(freq,[10*log10(abs(spectra(:,ii))) 10*log10(abs(Vxxxx))]);
 887         end
 888         title( char(plot_title(spect_type(ii),:)) );
 889         ylabel( 'amplitude' );
 890         %% Plot phase of cross spectrum and transfer function
 891         if ( spect_type(ii)==2 || spect_type(ii)==3 )
 892           figure();
 893           if ( plot_type==2 || plot_type==4 )
 894             semilogx(freq,180/pi*angle(spectra(:,ii)));
 895           else
 896             plot(freq,180/pi*angle(spectra(:,ii)));
 897           end
 898           title( char(plot_title(spect_type(ii),:)) );
 899           ylabel( 'phase' );
 900         end
 901       end %for
 902     end
 903   end
 904 end
 905
 906 %!demo
 907 %! fflush(stdout);
 908 %! rand('seed',2038014164);
 909 %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
 910 %! white = rand(1,16384);
 911 %! signal = detrend(filter(0.70181,a,white));
 912 %! % frequency shift by modulating with exp(j.omega.t)
 913 %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
 914 %! Fs = 25; % sampling frequency
 915 %! hold off
 916 %! pwelch([]);
 917 %! pwelch(signal);
 918 %! disp('Default settings: Fs=1Hz, overlap=0.5, no padding' )
 919 %! input('Onesided power spectral density (real data). Press ENTER', 's' );
 920 %! hold on
 921 %! pwelch(skewed);
 922 %! disp('Frequency-shifted complex data.  Twosided wrap-around spectrum.' );
 923 %! input('Area is same as one-sided spectrum. Press ENTER', 's' );
 924 %! pwelch(signal,'shift','semilogy');
 925 %! input('Twosided, centred zero-frequency, lin-log plot. Press ENTER', 's' );
 926 %! hold off
 927 %! figure();
 928 %! pwelch(skewed,[],[],[],Fs,'shift','semilogy');
 929 %! input('Actual Fs=25 Hz. Note change of scales. Press ENTER', 's' );
 930 %! pwelch(skewed,[],[],[],Fs,0.95,'shift','semilogy');
 931 %! input('Spectral density with 95% confidence interval. Press ENTER', 's' );
 932 %! pwelch('R12+');
 933 %! pwelch(signal,'squared');
 934 %! input('Spectral density with Matlab R12 defaults. Press ENTER', 's' );
 935 %! figure();
 936 %! pwelch([]);
 937 %! pwelch(signal,3640,[],4096,2*pi,[],'no-strip');
 938 %! input('Same spectrum with 95% confidence interval. Press ENTER', 's' );
 939 %! figure();
 940 %! pwelch(signal,[],[],[],2*pi,0.95,'no-strip');
 941 %! input('95% confidence interval with native defaults. Press ENTER', 's' );
 942 %! pwelch(signal,64,[],[],2*pi,'no-strip');
 943 %! input('Only 32 frequency values in this spectrum. Press ENTER', 's' );
 944 %! hold on
 945 %! pwelch(signal,64,[],256,2*pi,'no-strip');
 946 %! input('4:1 zero padding gives artificial smoothing. Press ENTER', 's' );
 947 %! figure();
 948 %! pwelch('psd');
 949 %! pwelch(signal,'squared');
 950 %! input('Just like Matlab spectrum.welch(...) defaults. Press ENTER', 's' );
 951 %! hold off
 952 %! pwelch({});
 953 %! pwelch(white,signal,'trans','coher','short')
 954 %! input('Transfer and coherence functions. Press ENTER', 's' );
 955 %! disp('Use "close all" to remove plotting windows.' );