holger@0
|
1 function c = dceps(X, w, numCoeffs, lambda)
|
holger@0
|
2 % computes the discrete cepstral coefficients of a sequence of partials X.
|
holger@0
|
3 %
|
holger@0
|
4 % c = dceps(X, w, numCoeffs)
|
holger@0
|
5 %
|
holger@0
|
6 % w contains the normalized frequencies of the partials. Normalized
|
holger@0
|
7 % frequencies should be in the range [0 .. pi]
|
holger@0
|
8 % X and w must be column vectors.
|
holger@0
|
9 % numCoeffs specifies the number of cepstral coefficients.
|
holger@0
|
10 %
|
holger@0
|
11 % Further details can be found in:
|
holger@0
|
12 % [1] Diemo Schwarz. Spectral envelopes in sound analysis and synthesis.
|
holger@0
|
13 % Master's thesis, Universitaet Stuttgart, 1998, pp. 36-38.
|
holger@0
|
14 % [2] T. Galas and X. Rodet. An improved cepstral method for deconvolution
|
holger@0
|
15 % of source filter systems with discrete spectra: Application to musical
|
holger@0
|
16 % sound signals. In International Computer Music Conference, 1990.
|
holger@0
|
17
|
holger@0
|
18
|
holger@0
|
19 p = numCoeffs-1;
|
holger@0
|
20 n = length(X); % number of partials
|
holger@0
|
21 S = ones(n,1); % amplitude of partials of source spectrum
|
holger@0
|
22 %lambda = 0.0005; % regularization term
|
holger@0
|
23
|
holger@0
|
24 if sum(X == 0) %if one or more element is zero
|
holger@0
|
25 error('All elements of X must be nonzero');
|
holger@0
|
26 end
|
holger@0
|
27
|
holger@0
|
28 % compute all r_i (i from 0 to 2p)
|
holger@0
|
29 r = .5 * sum(cos(w * (0:2*p)));
|
holger@0
|
30
|
holger@0
|
31 % compute matrix A
|
holger@0
|
32 A = zeros(p+1);
|
holger@0
|
33 B = zeros(p+1);
|
holger@0
|
34 for i = 0:p
|
holger@0
|
35 for j = 0:p
|
holger@0
|
36 A(i+1,j+1) = r(i+j+1) + r(abs(i-j)+1); % +1 because 1st idx in r is 0!
|
holger@0
|
37 end
|
holger@0
|
38 B(i+1,i+1) = 8 * (pi*(i-1))^2;
|
holger@0
|
39 end
|
holger@0
|
40
|
holger@0
|
41 A = A + lambda * B;
|
holger@0
|
42
|
holger@0
|
43 % compute b
|
holger@0
|
44 b = sum( repmat(log(X./S),1,p+1) .* cos(w *(0:p)) )';
|
holger@0
|
45
|
holger@0
|
46 % compute c
|
holger@0
|
47 c = A \ b; % equals inv(A) * b |