holger@0: function c = dceps(X, w, numCoeffs, lambda)
holger@0: % computes the discrete cepstral coefficients of a sequence of partials X.
holger@0: %
holger@0: % c = dceps(X, w, numCoeffs)
holger@0: %
holger@0: % w contains the normalized frequencies of the partials. Normalized
holger@0: % frequencies should be in the range [0 .. pi]
holger@0: % X and w must be column vectors.
holger@0: % numCoeffs specifies the number of cepstral coefficients.
holger@0: %
holger@0: % Further details can be found in:
holger@0: % [1] Diemo Schwarz. Spectral envelopes in sound analysis and synthesis.
holger@0: % Master's thesis, Universitaet Stuttgart, 1998, pp. 36-38.
holger@0: % [2] T. Galas and X. Rodet. An improved cepstral method for deconvolution
holger@0: % of source filter systems with discrete spectra: Application to musical
holger@0: % sound signals. In International Computer Music Conference, 1990.
holger@0: 
holger@0: 
holger@0: p = numCoeffs-1;
holger@0: n = length(X); % number of partials
holger@0: S = ones(n,1); % amplitude of partials of source spectrum
holger@0: %lambda = 0.0005; % regularization term
holger@0: 
holger@0: if sum(X == 0) %if one or more element is zero
holger@0:     error('All elements of X must be nonzero');
holger@0: end
holger@0: 
holger@0: % compute all r_i (i from 0 to 2p)
holger@0: r = .5 * sum(cos(w * (0:2*p)));
holger@0: 
holger@0: % compute matrix A
holger@0: A = zeros(p+1);
holger@0: B = zeros(p+1);
holger@0: for i = 0:p
holger@0:     for j = 0:p
holger@0:         A(i+1,j+1) = r(i+j+1) + r(abs(i-j)+1); % +1 because 1st idx in r is 0!
holger@0:     end
holger@0:     B(i+1,i+1) = 8 * (pi*(i-1))^2;
holger@0: end
holger@0: 
holger@0: A = A + lambda * B;
holger@0: 
holger@0: % compute b
holger@0: b = sum( repmat(log(X./S),1,p+1) .* cos(w *(0:p)) )';
holger@0: 
holger@0: % compute c
holger@0: c = A \ b; % equals inv(A) * b