function [cep,dctm] = spec2cep(spec, ncep, type) % [cep,dctm] = spec2cep(spec, ncep, type) % Calculate cepstra from spectral samples (in columns of spec) % Return ncep cepstral rows (defaults to 9) % This one does type II dct, or type I if type is specified as 1 % dctm returns the DCT matrix that spec was multiplied by to give cep. % 2005-04-19 dpwe@ee.columbia.edu for mfcc_dpwe if nargin < 2; ncep = 13; end if nargin < 3; type = 2; end % type of DCT [nrow, ncol] = size(spec); % Make the DCT matrix dctm = zeros(ncep, nrow); if type == 2 || type == 3 % this is the orthogonal one, the one you want for i = 1:ncep dctm(i,:) = cos((i-1)*[1:2:(2*nrow-1)]/(2*nrow)*pi) * sqrt(2/nrow); end if type == 2 % make it unitary! (but not for HTK type 3) dctm(1,:) = dctm(1,:)/sqrt(2); end elseif type == 4 % type 1 with implicit repeating of first, last bins % Deep in the heart of the rasta/feacalc code, there is the logic % that the first and last auditory bands extend beyond the edge of % the actual spectra, and they are thus copied from their neighbors. % Normally, we just ignore those bands and take the 19 in the middle, % but when feacalc calculates mfccs, it actually takes the cepstrum % over the spectrum *including* the repeated bins at each end. % Here, we simulate 'repeating' the bins and an nrow+2-length % spectrum by adding in extra DCT weight to the first and last % bins. for i = 1:ncep dctm(i,:) = cos((i-1)*[1:nrow]/(nrow+1)*pi) * 2; % Add in edge points at ends (includes fixup scale) dctm(i,1) = dctm(i,1) + 1; dctm(i,nrow) = dctm(i,nrow) + ((-1)^(i-1)); end dctm = dctm / (2*(nrow+1)); else % dpwe type 1 - same as old spec2cep that expanded & used fft for i = 1:ncep dctm(i,:) = cos((i-1)*[0:(nrow-1)]/(nrow-1)*pi) * 2 / (2*(nrow-1)); end % fixup 'non-repeated' points dctm(:,[1 nrow]) = dctm(:, [1 nrow])/2; end cep = dctm*log(spec);