annotate 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
rev   line source
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