octave_packages/signal-1.1.3/rceps.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: [y, xm] = rceps(x)
  17 ##   Produce the cepstrum of the signal x, and if desired, the minimum
  18 ##   phase reconstruction of the signal x.  If x is a matrix, do so
  19 ##   for each column of the matrix.
  20 ##
  21 ## Example
  22 ##   f0=70; Fs=10000;           # 100 Hz fundamental, 10kHz sampling rate
  23 ##   a=poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3])); # two formants
  24 ##   s=0.005*randn(1024,1);      # Noise excitation signal
  25 ##   s(1:Fs/f0:length(s)) = 1;   # Impulse glottal wave
  26 ##   x=filter(1,a,s);            # Speech signal in x
  27 ##   [y, xm] = rceps(x.*hanning(1024)); # cepstrum and min phase reconstruction
  28 ##
  29 ## Reference
  30 ##    Programs for digital signal processing. IEEE Press.
  31 ##    New York: John Wiley & Sons. 1979.
  32
  33 function [y, ym] = rceps(x)
  34   if (nargin != 1)
  35     print_usage;
  36   end
  37   y = real(ifft(log(abs(fft(x)))));
  38   if nargout == 2
  39     n=length(x);
  40     if rows(x)==1
  41       if rem(n,2)==1
  42         ym = [y(1), 2*y(2:n/2), zeros(1,n/2-1)];
  43       else
  44         ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
  45       endif
  46     else
  47       if rem(n,2)==1
  48         ym = [y(1,:); 2*y(2:n/2,:); zeros(n/2-1,columns(y))];
  49       else
  50         ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,columns(y))];
  51       endif
  52     endif
  53     ym = real(ifft(exp(fft(ym))));
  54   endif
  55 endfunction
  56
  57 %!error rceps
  58 %!error rceps(1,2)  # too many arguments
  59
  60 %!test
  61 %! ## accepts matrices
  62 %! x=randn(32,3);
  63 %! [y, xm] = rceps(x);
  64 %! ## check the mag-phase response of the reproduction
  65 %! hx = fft(x);
  66 %! hxm = fft(xm);
  67 %! assert(abs(hx), abs(hxm), 200*eps); # good magnitude response match
  68 %! #XXX FIXME XXX test for minimum phase?  Stop using random datasets!
  69 %! #assert(arg(hx) != arg(hxm));        # phase mismatch
  70
  71 %!test
  72 %! ## accepts column and row vectors
  73 %! x=randn(256,1);
  74 %! [y, xm] = rceps(x);
  75 %! [yt, xmt] = rceps(x.');
  76 %! tol = 1e-14;
  77 %! assert(yt.', y, tol);
  78 %! assert(xmt.', xm, tol);
  79
  80 %!demo
  81 %! f0=70; Fs=10000;           # 100 Hz fundamental, 10kHz sampling rate
  82 %! a=real(poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3]))); # two formants
  83 %! s=0.05*randn(1024,1);      # Noise excitation signal
  84 %! s(floor(1:Fs/f0:length(s))) = 1;  # Impulse glottal wave
  85 %! x=filter(1,a,s);           # Speech signal in x
  86 %! [y, xm] = rceps(x);        # cepstrum and minimum phase x
  87 %! [hx, w] = freqz(x,1,[],Fs); hxm = freqz(xm);
  88 %! figure(1);
  89 %! subplot(311);
  90 %!    auplot(x,Fs,'b',';signal;');
  91 %!    hold on; auplot(xm,Fs,'g',';reconstruction;');
  92 %!    hold off;
  93 %! subplot(312);
  94 %!    axis("ticy");
  95 %!    plot(w,log(abs(hx)), ";magnitude;", ...
  96 %!         w,log(abs(hxm)),";reconstruction;");
  97 %! subplot(313);
  98 %!    axis("on");
  99 %!    plot(w,unwrap(arg(hx))/(2*pi), ";phase;",...
 100 %!         w,unwrap(arg(hxm))/(2*pi),";reconstruction;");
 101 %! figure(2); auplot(y,Fs,';cepstrum;');
 102 %! %-------------------------------------------------------------
 103 %! % confirm the magnitude spectrum is identical in the signal
 104 %! % and the reconstruction and that there are peaks in the
 105 %! % cepstrum at 14 ms intervals corresponding to an F0 of 70 Hz.