/usr/share/octave/packages/signal-1.3.0/rceps.m is in octave-signal 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 | ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{y}, @var{xm}] =} rceps (@var{x})
## Produce the cepstrum of the signal x, and if desired, the minimum
## phase reconstruction of the signal x. If x is a matrix, do so
## for each column of the matrix.
##
## Example:
## @example
## @group
## f0 = 70; Fs = 10000; # 100 Hz fundamental, 10kHz sampling rate
## a = poly (0.985 * exp (1i*pi*[0.1, -0.1, 0.3, -0.3])); # two formants
## s = 0.005 * randn (1024, 1); # Noise excitation signal
## s(1:Fs/f0:length(s)) = 1; # Impulse glottal wave
## x = filter (1, a, s); # Speech signal in x
## [y, xm] = rceps (x .* hanning (1024)); # cepstrum and min phase reconstruction
## @end group
## @end example
##
## Reference
## Programs for digital signal processing. IEEE Press.
## New York: John Wiley & Sons. 1979.
## @end deftypefn
function [y, ym] = rceps(x)
if (nargin != 1)
print_usage;
endif
f = abs(fft(x));
if (any (f == 0))
error ("rceps: the spectrum of x contains zeros, unable to compute real cepstrum");
endif
y = real(ifft(log(f)));
if nargout == 2
n=length(x);
if rows(x)==1
if rem(n,2)==1
ym = [y(1), 2*y(2:n/2+1), zeros(1,n/2)];
else
ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
endif
else
if rem(n,2)==1
ym = [y(1,:); 2*y(2:n/2+1,:); zeros(n/2,columns(y))];
else
ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,columns(y))];
endif
endif
ym = real(ifft(exp(fft(ym))));
endif
endfunction
%!error rceps
%!error rceps(1,2) # too many arguments
%!test
%! ## accepts matrices
%! x=randn(32,3);
%! [y, xm] = rceps(x);
%! ## check the mag-phase response of the reproduction
%! hx = fft(x);
%! hxm = fft(xm);
%! assert(abs(hx), abs(hxm), 200*eps); # good magnitude response match
%! ## FIXME: test for minimum phase? Stop using random datasets!
%! #assert(arg(hx) != arg(hxm)); # phase mismatch
%!test
%! ## accepts column and row vectors
%! x=randn(256,1);
%! [y, xm] = rceps(x);
%! [yt, xmt] = rceps(x.');
%! tol = 1e-14;
%! assert(yt.', y, tol);
%! assert(xmt.', xm, tol);
%% Test that an odd-length input produces an odd-length output
%!test
%! x = randn(33, 4);
%! [y, xm] = rceps(x);
%! assert(size(y) == size(x));
%! assert(size(xm) == size(x));
%!demo
%! f0=70; Fs=10000; # 100 Hz fundamental, 10kHz sampling rate
%! a=real(poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3]))); # two formants
%! s=0.05*randn(1024,1); # Noise excitation signal
%! s(floor(1:Fs/f0:length(s))) = 1; # Impulse glottal wave
%! x=filter(1,a,s); # Speech signal in x
%! [y, xm] = rceps(x); # cepstrum and minimum phase x
%! [hx, w] = freqz(x,1,[],Fs); hxm = freqz(xm);
%! figure(1);
%! subplot(311);
%! len = 1000 * fix (min (length (x), length (xm)) / 1000);
%! plot([0:len-1]*1000/Fs,x(1:len),'b;signal;');
%! hold on; plot([0:len-1]*1000/Fs,xm(1:len),'g;reconstruction;');
%! ylabel("amplitude"); xlabel("time (ms)");
%! hold off;
%! subplot(312);
%! axis("ticy");
%! plot(w,log(abs(hx)), ";magnitude;", ...
%! w,log(abs(hxm)),";reconstruction;");
%! xlabel("frequency (Hz)");
%! subplot(313);
%! axis("on");
%! plot(w,unwrap(arg(hx))/(2*pi), ";phase;",...
%! w,unwrap(arg(hxm))/(2*pi),";reconstruction;");
%! xlabel("frequency (Hz)");
%! len = 1000 * fix (length (y) / 1000);
%! figure(2); plot([0:len-1]*1000/Fs,y(1:len),';cepstrum;');
%! ylabel("amplitude"); xlabel("quefrency (ms)");
%! %-------------------------------------------------------------
%! % confirm the magnitude spectrum is identical in the signal
%! % and the reconstruction and that there are peaks in the
%! % cepstrum at 14 ms intervals corresponding to an F0 of 70 Hz.
|