function y = hmmfeatures(s,N,deltaN,M,Q)
%  hmmfeatures --> Feature extraction for HMM recognizer.
%
%  <Synopsis>
%    y = hmmfeatures(s,N,deltaN,M,Q)
%
%  <Description>
%    A frame based analysis of the speech signal, s, is performed to
%    give observation vectors (columns of y), which can be used to train
%    HMMs for speech recognition.
%
%    The speech signal is blocked into frames of N samples, and
%    consecutive frames are spaced deltaN samples apart. Each frame is
%    multiplied by an N-sample Hamming window, and Mth-order LP analysis
%    is performed. The LPC coefficients are then converted to Q cepstral
%    coefficients, which are weighted by a raised sine window. The result
%    is the first half of an observation vector, the second half is the
%    differenced cepstral coefficients used to add dynamic information.
%    Thus, the returned argument y is an 2Q-by-T matrix, where T is the
%    number of frames.
%
%  <See Also>
%    hmmcodebook --> Codebook generation for HMM recognizer.

%  <References>
%  [1] J.R Deller, J.G. Proakis and F.H.L. Hansen, "Discrete-Time
%      Processing of Speech Signals", IEEE Press, chapter 12, (2000).
%
%  <Revision>
%    Peter S.K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: September 30, 2000
%-----------------------------------------------------------------------

Ns = length(s);                         % Signal length.
T  = 1 + fix((Ns-N)/deltaN);            % No. of frames.

a       = zeros(Q,1);
gamma   = zeros(Q,1);
gamma_w = zeros(Q,T);

win_gamma = 1 + (Q/2)*sin(pi/Q*(1:Q)'); % Cepstral window function.

for (t = 1:T)                           % Loop frames.
  % Block into frames.
  idx = (deltaN*(t-1)+1):(deltaN*(t-1)+N);

  % Window frame.
  sw = s(idx).*hamming(N);

  % Short-term autocorrelation.
  [rs,eta] = xcorr(sw,M,'biased');

  % LP analysis based on Levinson-Durbin recursion.
  [a(1:M),xi,kappa] = durbin(rs(M+1:2*M+1),M);

  % Cepstral coefficients.
  gamma(1) = a(1);
  for (i = 2:Q)
    gamma(i) = a(i) + (1:i-1)*(gamma(1:i-1).*a(i-1:-1:1))/i;
  end

  % Weighted cepstral sequence for frame t.
  gamma_w(:,t) = gamma.*win_gamma;
end

% Time differenced weighted cepstral sequence.
delta_gamma_w = gradient(gamma_w);

% Observation vectors.
y = [gamma_w; delta_gamma_w];

%-----------------------------------------------------------------------
% End of function hmmfeatures
%-----------------------------------------------------------------------
