comparison functions/dceps.m @ 0:b4e26b53072f tip

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