octave_packages/signal-1.1.3/pburg.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 %%    [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)
  18 %%
  19 %% Calculate Burg maximum-entropy power spectral density.
  20 %% The functions "arburg" and "ar_psd" do all the work.
  21 %% See "help arburg" and "help ar_psd" for further details.
  22 %%
  23 %% ARGUMENTS:
  24 %%     All but the first two arguments are optional and may be empty.
  25 %%   x       %% [vector] sampled data
  26 %%
  27 %%   poles   %% [integer scalar] required number of poles of the AR model
  28 %%
  29 %%   freq    %% [real vector] frequencies at which power spectral density
  30 %%           %%               is calculated
  31 %%           %% [integer scalar] number of uniformly distributed frequency
  32 %%           %%          values at which spectral density is calculated.
  33 %%           %%          [default=256]
  34 %%
  35 %%   Fs      %% [real scalar] sampling frequency (Hertz) [default=1]
  36 %%
  37 %%
  38 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
  39 %%   Control-string arguments can be in any order after the other arguments.
  40 %%
  41 %%
  42 %%   range   %% 'half',  'onesided' : frequency range of the spectrum is
  43 %%           %%       from zero up to but not including sample_f/2.  Power
  44 %%           %%       from negative frequencies is added to the positive
  45 %%           %%       side of the spectrum.
  46 %%           %% 'whole', 'twosided' : frequency range of the spectrum is
  47 %%           %%       -sample_f/2 to sample_f/2, with negative frequencies
  48 %%           %%       stored in "wrap around" order after the positive
  49 %%           %%       frequencies; e.g. frequencies for a 10-point 'twosided'
  50 %%           %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
  51 %%           %% 'shift', 'centerdc' : same as 'whole' but with the first half
  52 %%           %%       of the spectrum swapped with second half to put the
  53 %%           %%       zero-frequency value in the middle. (See "help
  54 %%           %%       fftshift". If "freq" is vector, 'shift' is ignored.
  55 %%           %% If model coefficients "ar_coeffs" are real, the default
  56 %%           %% range is 'half', otherwise default range is 'whole'.
  57 %%
  58 %%   method  %% 'fft':  use FFT to calculate power spectral density.
  59 %%           %% 'poly': calculate spectral density as a polynomial of 1/z
  60 %%           %% N.B. this argument is ignored if the "freq" argument is a
  61 %%           %%      vector.  The default is 'poly' unless the "freq"
  62 %%           %%      argument is an integer power of 2.
  63 %%
  64 %% plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
  65 %%           %% specifies the type of plot.  The default is 'plot', which
  66 %%           %% means linear-linear axes. 'squared' is the same as 'plot'.
  67 %%           %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
  68 %%           %% spectrum is not plotted if the caller requires a returned
  69 %%           %% value.
  70 %%
  71 %% criterion %% [optional string arg]  model-selection criterion.  Limits
  72 %%           %%       the number of poles so that spurious poles are not
  73 %%           %%       added when the whitened data has no more information
  74 %%           %%       in it (see Kay & Marple, 1981). Recognised values are
  75 %%           %%  'AKICc' -- approximate corrected Kullback information
  76 %%           %%             criterion (recommended),
  77 %%           %%   'KIC'  -- Kullback information criterion
  78 %%           %%   'AICc' -- corrected Akaike information criterion
  79 %%           %%   'AIC'  -- Akaike information criterion
  80 %%           %%   'FPE'  -- final prediction error" criterion
  81 %%           %% The default is to NOT use a model-selection criterion
  82 %%
  83 %% RETURNED VALUES:
  84 %%     If return values are not required by the caller, the spectrum
  85 %%     is plotted and nothing is returned.
  86 %%   psd       %% [real vector] power-spectral density estimate
  87 %%   f_out     %% [real vector] frequency values
  88 %%
  89 %% HINTS
  90 %%   This function is a wrapper for arburg and ar_psd.
  91 %%   See "help arburg", "help ar_psd".
  92
  93 function [psd,f_out]=pburg(x,poles,varargin)
  94   %%
  95   if ( nargin<2 )
  96     error( 'pburg: need at least 2 args. Use "help pburg"' );
  97   end
  98   nvarargin=length(varargin);
  99   criterion=[];
 100   %%
 101   %% Search for a "criterion" arg. If found, remove it
 102   %% from "varargin" list and feed it to arburg instead.
 103   for iarg = 1: nvarargin
 104     arrgh = varargin{iarg};
 105     if ( ischar(arrgh) && ( strcmp(arrgh,'AKICc') ||...
 106          strcmp(arrgh,'KIC') || strcmp(arrgh,'AICc') ||...
 107          strcmp(arrgh,'AIC') || strcmp(arrgh,'FPE') ) )
 108       criterion=arrgh;
 109       if ( nvarargin>1 )
 110         varargin{iarg}= [];
 111       else
 112         varargin={};
 113         end
 114       end
 115     end
 116   %%
 117   [ar_coeffs,residual]=arburg(x,poles,criterion);
 118   if ( nargout==0 )
 119     ar_psd(ar_coeffs,residual,varargin{:});
 120   elseif ( nargout==1 )
 121     psd = ar_psd(ar_coeffs,residual,varargin{:});
 122   elseif ( nargout>=2 )
 123     [psd,f_out] = ar_psd(ar_coeffs,residual,varargin{:});
 124   end
 125 end
 126
 127 %!demo
 128 %! fflush(stdout);
 129 %! rand('seed',2038014164);
 130 %! a = [ 1.0 -1.6216505 1.1102795 -0.4621741 0.2075552 -0.018756746 ];
 131 %! signal = detrend(filter(0.70181,a,rand(1,16384)));
 132 %! % frequency shift by modulating with exp(j.omega.t)
 133 %! skewed = signal.*exp(2*pi*i*2/25*[1:16384]);
 134 %! Fs = 25;
 135 %! hold on
 136 %! pburg(signal,3,[],Fs);
 137 %! input('Onesided 3-pole spectrum. Press ENTER', 's' );
 138 %! pburg(signal,4,[],Fs,'whole');
 139 %! input('Twosided 4-pole spectrum of same data. Press ENTER', 's' );
 140 %! pburg(signal,5,128,Fs,'shift', 'semilogy');
 141 %! input('Twosided, centred zero-frequency, 5-pole. Press ENTER', 's' );
 142 %! pburg(skewed,7,128,Fs,'AKICc','shift','semilogy');
 143 %! input('Complex data, AKICc chooses no. of poles. Press ENTER', 's' );
 144 %! user_freq=[-0.2:0.02:0.2]*Fs;
 145 %! pburg(skewed,7,user_freq,Fs,'AKICc','semilogy');
 146 %! input('User-specified frequency values. Press ENTER', 's' );
 147 %! hold off
 148 %! clf