octave_packages/signal-1.1.3/ar_psd.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] = ar_psd(a,v,freq,Fs,range,method,plot_type)
  18 %%
  19 %%  Calculate the power spectrum of the autoregressive model
  20 %%
  21 %%                         M
  22 %%  x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
  23 %%                        k=1
  24 %%  where x(n) is the output of the model and e(n) is white noise.
  25 %%  This function is intended for use with
  26 %%    [a,v,k] = arburg(x,poles,criterion)
  27 %%  which use the Burg (1968) method to calculate a "maximum entropy"
  28 %%  autoregressive model of "x".  This function runs on octave and matlab.
  29 %%
  30 %%  If the "freq" argument is a vector (of frequencies) the spectrum is
  31 %%  calculated using the polynomial method and the "method" argument is
  32 %%  ignored.  For scalar "freq", an integer power of 2, or "method='FFT'",
  33 %%  causes the spectrum to be calculated by FFT.  Otherwise, the spectrum
  34 %%  is calculated as a polynomial.  It may be computationally more
  35 %%  efficient to use the FFT method if length of the model is not much
  36 %%  smaller than the number of frequency values. The spectrum is scaled so
  37 %%  that spectral energy (area under spectrum) is the same as the
  38 %%  time-domain energy (mean square of the signal).
  39 %%
  40 %% ARGUMENTS:
  41 %%     All but the first two arguments are optional and may be empty.
  42 %%
  43 %%   a      %% [vector] list of M=(order+1) autoregressive model
  44 %%          %%      coefficients.  The first element of "ar_coeffs" is the
  45 %%          %%      zero-lag coefficient, which always has a value of 1.
  46 %%
  47 %%   v      %% [real scalar] square of the moving-average coefficient of
  48 %%          %%               the AR model.
  49 %%
  50 %%   freq   %% [real vector] frequencies at which power spectral density
  51 %%          %%               is calculated
  52 %%          %% [integer scalar] number of uniformly distributed frequency
  53 %%          %%          values at which spectral density is calculated.
  54 %%          %%          [default=256]
  55 %%
  56 %%   Fs     %% [real scalar] sampling frequency (Hertz) [default=1]
  57 %%
  58 %% CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
  59 %%   Control-string arguments can be in any order after the other arguments.
  60 %%
  61 %%   range  %% 'half',  'onesided' : frequency range of the spectrum is
  62 %%          %%       from zero up to but not including sample_f/2.  Power
  63 %%          %%       from negative frequencies is added to the positive
  64 %%          %%       side of the spectrum.
  65 %%          %% 'whole', 'twosided' : frequency range of the spectrum is
  66 %%          %%       -sample_f/2 to sample_f/2, with negative frequencies
  67 %%          %%       stored in "wrap around" order after the positive
  68 %%          %%       frequencies; e.g. frequencies for a 10-point 'twosided'
  69 %%          %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
  70 %%          %% 'shift', 'centerdc' : same as 'whole' but with the first half
  71 %%          %%       of the spectrum swapped with second half to put the
  72 %%          %%       zero-frequency value in the middle. (See "help
  73 %%          %%       fftshift". If "freq" is vector, 'shift' is ignored.
  74 %%          %% If model coefficients "ar_coeffs" are real, the default
  75 %%          %% range is 'half', otherwise default range is 'whole'.
  76 %%
  77 %%   method %% 'fft':  use FFT to calculate power spectrum.
  78 %%          %% 'poly': calculate power spectrum as a polynomial of 1/z
  79 %%          %% N.B. this argument is ignored if the "freq" argument is a
  80 %%          %%      vector.  The default is 'poly' unless the "freq"
  81 %%          %%      argument is an integer power of 2.
  82 %%
  83 %% plot_type%% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
  84 %%          %% specifies the type of plot.  The default is 'plot', which
  85 %%          %% means linear-linear axes. 'squared' is the same as 'plot'.
  86 %%          %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
  87 %%          %% spectrum is not plotted if the caller requires a returned
  88 %%          %% value.
  89 %%
  90 %% RETURNED VALUES:
  91 %%     If returned values are not required by the caller, the spectrum
  92 %%     is plotted and nothing is returned.
  93 %%   psd    %% [real vector] estimate of power-spectral density
  94 %%   f_out  %% [real vector] frequency values
  95 %%
  96 %% N.B. arburg runs in octave and matlab, and does not depend on octave-forge
  97 %%      or signal-processing-toolbox functions.
  98 %%
  99 %% REFERENCE
 100 %% [1] Equation 2.28 from Steven M. Kay and Stanley Lawrence Marple Jr.:
 101 %%   "Spectrum analysis -- a modern perspective",
 102 %%   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981
 103 %%
 104
 105 function [varargout]=ar_psd(a,v,varargin)
 106 %%
 107 %% Check fixed arguments
 108 if ( nargin < 2 )
 109   error( 'ar_psd: needs at least 2 args. Use help ar_psd.' );
 110 elseif ( ~isvector(a) || length(a)<2 )
 111   error( 'ar_psd: arg 1 (a) must be vector, length>=2.' );
 112 elseif ( ~isscalar(v) )
 113   error( 'ar_psd: arg 2 (v) must be real scalar >0.' );
 114 else
 115   real_model = isreal(a);
 116 %%
 117 %%  default values for optional areguments
 118   freq = 256;
 119   user_freqs = 0;    %% boolean: true for user-specified frequencies
 120   Fs   = 1.0;
 121   %%  FFT padding factor (is also frequency range divisor): 1=whole, 2=half.
 122   pad_fact = 1 + real_model;
 123   do_shift   = 0;
 124   force_FFT  = 0;
 125   force_poly = 0;
 126   plot_type  = 1;
 127 %%
 128 %%  decode and check optional arguments
 129 %%  end_numeric_args is boolean; becomes true at 1st string arg
 130   end_numeric_args = 0;
 131   for iarg = 1:length(varargin)
 132     arg = varargin{iarg};
 133     end_numeric_args = end_numeric_args || ischar(arg);
 134     %% skip empty arguments
 135     if ( isempty(arg) )
 136       1;
 137     %% numeric optional arguments must be first, cannot follow string args
 138     %% N.B. older versions of matlab may not have the function "error" so
 139     %% the user writes "function error(msg); disp(msg); end" and we need
 140     %% a "return" here.
 141     elseif ( ~ischar(arg) )
 142       if ( end_numeric_args )
 143         error( 'ar_psd: control arg must be string.' );
 144       %%
 145       %% first optional numeric arg is "freq"
 146       elseif ( iarg == 1 )
 147         user_freqs = isvector(arg) && length(arg)>1;
 148         if ( ~isscalar(arg) && ~user_freqs )
 149           error( 'ar_psd: arg 3 (freq) must be vector or scalar.' );
 150         elseif ( ~user_freqs && ( ~isreal(arg) || ...
 151                  fix(arg)~=arg || arg <= 2 || arg >= 1048576 ) )
 152           error('ar_psd: arg 3 (freq) must be integer >=2, <=1048576' );
 153         elseif ( user_freqs && ~isreal(arg) )
 154           error( 'ar_psd: arg 3 (freq) vector must be real.' );
 155         end
 156         freq = arg(:); % -> column vector
 157       %%
 158       %% second optional numeric arg is  "Fs" - sampling frequency
 159       elseif ( iarg == 2 )
 160         if ( ~isscalar(arg) || ~isreal(arg) || arg<=0 )
 161           error( 'ar_psd: arg 4 (Fs) must be real positive scalar.' );
 162         end
 163         Fs = arg;
 164       %%
 165       else
 166         error( 'ar_psd: control arg must be string.' );
 167       end
 168   %%
 169   %% decode control-string arguments
 170     elseif ( strcmp(arg,'plot') || strcmp(arg,'squared') )
 171       plot_type = 1;
 172     elseif ( strcmp(arg,'semilogx') )
 173       plot_type = 2;
 174     elseif ( strcmp(arg,'semilogy') )
 175       plot_type = 3;
 176     elseif ( strcmp(arg,'loglog') )
 177       plot_type = 4;
 178     elseif ( strcmp(arg,'dB') )
 179       plot_type = 5;
 180     elseif ( strcmp(arg,'fft') )
 181       force_FFT  = 1;
 182       force_poly = 0;
 183     elseif ( strcmp(arg,'poly') )
 184       force_FFT  = 0;
 185       force_poly = 1;
 186     elseif ( strcmp(arg,'half') || strcmp(arg,'onesided') )
 187       pad_fact = 2;    % FFT zero-padding factor (pad FFT to double length)
 188       do_shift = 0;
 189     elseif ( strcmp(arg,'whole') || strcmp(arg,'twosided') )
 190       pad_fact = 1;    % FFT zero-padding factor (do not pad)
 191       do_shift = 0;
 192     elseif ( strcmp(arg,'shift') || strcmp(arg,'centerdc') )
 193       pad_fact = 1;
 194       do_shift = 1;
 195     else
 196       error( 'ar_psd: string arg: illegal value: %s', arg );
 197     end
 198   end
 199 %%  end of decoding and checking args
 200 %%
 201   if ( user_freqs )
 202     %% user provides (column) vector of frequencies
 203     if ( any(abs(freq)>Fs/2) )
 204       error( 'ar_psd: arg 3 (freq) cannot exceed half sampling frequency.' );
 205     elseif ( pad_fact==2 && any(freq<0) )
 206       error( 'ar_psd: arg 3 (freq) must be positive in onesided spectrum' );
 207     end
 208     freq_len = length(freq);
 209     fft_len  = freq_len;
 210     use_FFT  = 0;
 211     do_shift = 0;
 212   else
 213     %% internally generated frequencies
 214     freq_len = freq;
 215     freq = (Fs/pad_fact/freq_len) * [0:freq_len-1]';
 216     %% decide which method to use (poly or FFT)
 217     is_power_of_2 = rem(log(freq_len),log(2))<10.*eps;
 218     use_FFT = ( ~ force_poly && is_power_of_2 ) || force_FFT;
 219     fft_len = freq_len * pad_fact;
 220     end
 221   end
 222   %%
 223   %% calculate denominator of Equation 2.28, Kay and Marple, ref [1]Jr.:
 224   len_coeffs = length(a);
 225   if ( use_FFT )
 226     %% FFT method
 227     fft_out = fft( [ a(:); zeros(fft_len-len_coeffs,1) ] );
 228   else
 229     %% polynomial method
 230     %% complex data on "half" frequency range needs -ve frequency values
 231     if ( pad_fact==2 && ~real_model )
 232       freq = [freq; -freq(freq_len:-1:1)];
 233       fft_len = 2*freq_len;
 234     end
 235     fft_out = polyval( a(len_coeffs:-1:1), exp( (-i*2*pi/Fs) * freq ) );
 236   end
 237   %%
 238   %% The power spectrum (PSD) is the scaled squared reciprocal of amplitude
 239   %% of the FFT/polynomial. This is NOT the reciprocal of the periodogram.
 240   %% The PSD is a continuous function of frequency.  For uniformly
 241   %% distributed frequency values, the FFT algorithm might be the most
 242   %% efficient way of calculating it.
 243   %%
 244   psd = ( v / Fs ) ./ ( fft_out .* conj(fft_out) );
 245   %%
 246   %% range='half' or 'onesided',
 247   %%   add PSD at -ve frequencies to PSD at +ve frequencies
 248   %% N.B. unlike periodogram, PSD at zero frequency _is_ doubled.
 249   if ( pad_fact==2 )
 250     freq = freq(1:freq_len);
 251     if ( real_model )
 252       %% real data, double the psd
 253       psd = 2 * psd(1:freq_len);
 254     elseif ( use_FFT )
 255       %% complex data, FFT method, internally-generated frequencies
 256       psd = psd(1:freq_len)+[psd(1); psd(fft_len:-1:freq_len+2)];
 257     else
 258       %% complex data, polynomial method
 259       %%  user-defined and internally-generated frequencies
 260       psd = psd(1:freq_len)+psd(fft_len:-1:freq_len+1);
 261     end
 262   %%
 263   %% range='shift'
 264   %%   disabled for user-supplied frequencies
 265   %%   Shift zero-frequency to the middle (pad_fact==1)
 266   elseif ( do_shift )
 267     len2 = fix((fft_len+1)/2);
 268     psd  = [psd(len2+1:fft_len); psd(1:len2)];
 269     freq = [freq(len2+1:fft_len)-Fs; freq(1:len2)];
 270   end
 271   %%
 272   %% Plot the spectrum if there are no return variables.
 273   if ( nargout >= 2 )
 274      varargout{1} = psd;
 275      varargout{2} = freq;
 276   elseif ( nargout == 1 )
 277      varargout{1} = psd;
 278   else
 279     if ( plot_type == 1 )
 280       plot(freq,psd);
 281     elseif ( plot_type == 2 )
 282       semilogx(freq,psd);
 283     elseif ( plot_type == 3 )
 284       semilogy(freq,psd);
 285     elseif ( plot_type == 4 )
 286       loglog(freq,psd);
 287     elseif ( plot_type == 5 )
 288       plot(freq,10*log10(psd));
 289     end
 290   end
 291 end