diff functions/dceps.m @ 0:b4e26b53072f tip

Initial commit.

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