octave_packages/signal-1.1.3/cceps.m

   1 ## Copyright (C) 1994 Dept of Probability Theory and Statistics TU Wien <Andreas.Weingessel@ci.tuwien.ac.at>
   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:  cceps (x [, correct])
  17 ##
  18 ## Returns the complex cepstrum of the vector x.
  19 ## If the optional argument correct has the value 1, a correction
  20 ## method is applied.  The default is not to do this.
  21
  22 ## Author: Andreas Weingessel <Andreas.Weingessel@ci.tuwien.ac.at>
  23 ## Apr 1, 1994
  24 ## Last modifified by AW on Nov 8, 1994
  25
  26 function cep = cceps (x, c)
  27
  28   if (nargin == 1)
  29     c = 0;
  30   elseif (nargin != 2)
  31     print_usage;
  32   endif
  33
  34   [nr, nc] = size (x);
  35   if (nc != 1)
  36     if (nr == 1)
  37       x = x';
  38       nr = nc;
  39     else
  40       error ("cceps: x must be a vector");
  41     endif
  42   endif
  43
  44   bad_signal_message = ["cceps:  bad signal x, ", ...
  45       "some Fourier coefficients are zero."];
  46
  47   F = fft (x);
  48   if (min (abs (F)) == 0)
  49     error (bad_signal_message);
  50   endif
  51
  52   # determine if correction necessary
  53   half = fix (nr / 2);
  54   cor = 0;
  55   if (2 * half == nr)
  56     cor = (c && (real (F (half + 1)) < 0));
  57     if (cor)
  58       F = fft (x(1:nr-1))
  59       if (min (abs (F)) == 0)
  60         error (bad_signal_message);
  61       endif
  62     endif
  63   endif
  64
  65   cep = fftshift (ifft (log (F)));
  66
  67   # make result real
  68   if (c)
  69     cep = real (cep);
  70     if (cor)
  71       # make cepstrum of same length as input vector
  72       cep (nr) = 0;
  73     endif
  74   endif
  75
  76 endfunction