1 ## Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
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
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
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/>.
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.
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
30 ## Programs for digital signal processing. IEEE Press.
31 ## New York: John Wiley & Sons. 1979.
33 function [y, ym] = rceps(x)
37 y = real(ifft(log(abs(fft(x)))));
42 ym = [y(1), 2*y(2:n/2), zeros(1,n/2-1)];
44 ym = [y(1), 2*y(2:n/2), y(n/2+1), zeros(1,n/2-1)];
48 ym = [y(1,:); 2*y(2:n/2,:); zeros(n/2-1,columns(y))];
50 ym = [y(1,:); 2*y(2:n/2,:); y(n/2+1,:); zeros(n/2-1,columns(y))];
53 ym = real(ifft(exp(fft(ym))));
58 %!error rceps(1,2) # too many arguments
61 %! ## accepts matrices
63 %! [y, xm] = rceps(x);
64 %! ## check the mag-phase response of the reproduction
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
72 %! ## accepts column and row vectors
74 %! [y, xm] = rceps(x);
75 %! [yt, xmt] = rceps(x.');
77 %! assert(yt.', y, tol);
78 %! assert(xmt.', xm, tol);
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);
90 %! auplot(x,Fs,'b',';signal;');
91 %! hold on; auplot(xm,Fs,'g',';reconstruction;');
95 %! plot(w,log(abs(hx)), ";magnitude;", ...
96 %! w,log(abs(hxm)),";reconstruction;");
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.