octave_packages/signal-1.1.3/__power.m

   1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
   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:  [P, w] = __power (b, a, [, nfft [, Fs]] [, range] [, units])
  17 ##
  18 ## Plot the power spectrum of the given ARMA model.
  19 ##
  20 ## b, a: filter coefficients (b=numerator, a=denominator)
  21 ## nfft is number of points at which to sample the power spectrum
  22 ## Fs is the sampling frequency of x
  23 ## range is 'half' (default) or 'whole'
  24 ## units is  'squared' or 'db' (default)
  25 ## range and units may be specified any time after the filter, in either
  26 ## order
  27 ##
  28 ## Returns P, the magnitude vector, and w, the frequencies at which it
  29 ## is sampled.  If there are no return values requested, then plot the power
  30 ## spectrum and don't return anything.
  31
  32 ## TODO: consider folding this into freqz --- just one more parameter to
  33 ## TODO:    distinguish between 'linear', 'log', 'logsquared' and 'squared'
  34
  35 function [varargout] = __power (b, a, varargin)
  36   if (nargin < 2 || nargin > 6) print_usage; endif
  37
  38   nfft = [];
  39   Fs = [];
  40   range = [];
  41   range_fact = 1.0;
  42   units = [];
  43
  44   pos = 0;
  45   for i=1:length(varargin)
  46     arg = varargin{i};
  47     if strcmp(arg, 'squared') || strcmp(arg, 'db')
  48       units = arg;
  49     elseif strcmp(arg, 'whole')
  50       range = arg;
  51       range_fact = 1.0;
  52     elseif strcmp(arg, 'half')
  53       range = arg;
  54       range_fact = 2.0;
  55     elseif ischar(arg)
  56       usage(usagestr);
  57     elseif pos == 0
  58       nfft = arg;
  59       pos++;
  60     elseif pos == 1
  61       Fs = arg;
  62       pos++;
  63     else
  64       usage(usagestr);
  65     endif
  66   endfor
  67
  68   if isempty(nfft); nfft = 256; endif
  69   if isempty(Fs); Fs = 2; endif
  70   if isempty(range)
  71     range = 'half';
  72     range_fact = 2.0;
  73     endif
  74
  75   [P, w] = freqz(b, a, nfft, range, Fs);
  76
  77   P = (range_fact/Fs)*(P.*conj(P));
  78   if nargout == 0,
  79     if strcmp(units, 'squared')
  80       plot(w, P, ";;");
  81     else
  82       plot(w, 10.0*log10(abs(P)), ";;");
  83     endif
  84   endif
  85   if nargout >= 1, varargout{1} = P; endif
  86   if nargout >= 2, varargout{2} = w; endif
  87
  88 endfunction
  89