def test_definition_ortho(self):
     """Test orthornomal mode."""
     for i in range(len(X)):
         x = np.array(X[i], dtype=self.rdt)
         y = dct(x, norm="ortho", type=2)
         xi = dct(y, norm="ortho", type=3)
         self.assertTrue(xi.dtype == self.rdt, "Output dtype is %s, expected %s" % (xi.dtype, self.rdt))
         assert_array_almost_equal(xi, x, decimal=self.dec)

def dct_2d_ref(x, **kwargs):
    """Calculate reference values for testing dct2."""
    x = np.array(x, copy=True)
    for row in range(x.shape[0]):
        x[row, :] = dct(x[row, :], **kwargs)
    for col in range(x.shape[1]):
        x[:, col] = dct(x[:, col], **kwargs)
    return x

 def test_definition_ortho(self):
     # Test orthornomal mode.
     for i in range(len(X)):
         x = np.array(X[i], dtype=self.rdt)
         y = dct(x, norm='ortho', type=2)
         xi = dct(y, norm="ortho", type=3)
         assert_equal(xi.dtype, self.rdt)
         assert_array_almost_equal(xi, x, decimal=self.dec)

    def test_axis(self):
        nt = 2
        for i in [7, 8, 9, 16, 32, 64]:
            x = np.random.randn(nt, i)
            y = dct(x, type=self.type)
            for j in range(nt):
                assert_array_almost_equal(y[j], dct(x[j], type=self.type), decimal=self.dec)

            x = x.T
            y = dct(x, axis=0, type=self.type)
            for j in range(nt):
                assert_array_almost_equal(y[:, j], dct(x[:, j], type=self.type), decimal=self.dec)

def st_mfcc(cur_pos_signal, fbank, nceps):
    """
    短时mfcc
    """
    mspec = numpy.log10(numpy.dot(cur_pos_signal, fbank.T) + EPS)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:nceps]
    return ceps

 def smoothData(self,x,y,weight,nMiss=0):
   '''
   smooth data
   '''
   import scipy.optimize.lbfgsb as lbfgsb
   from scipy.fftpack.realtransforms import dct,idct
   n0 = len(x)
   #x = np.array([x,x,x]).flatten()
   #y = np.array([y,y,y]).flatten()
   #weight = np.array([weight,weight,weight]).flatten()
   n = len(x)
   weight = 1./weight
   # scale 0 to 1
   weight = weight/np.max(weight)
   i = np.arange(1,n+1)
   eigenvalues = -2. + 2.*np.cos((i-1)*np.pi/n)
   DCTy = dct(y,norm='ortho',type=2)
   dcty2 = DCTy**2
   eigenvalues2 = eigenvalues**2
   x0 = np.atleast_1d(1.)
   y_hat = np.zeros_like(y)
   xpost,f,d = lbfgsb.fmin_l_bfgs_b(gcv,x0,fprime=None,factr=10.,\
          approx_grad=True,args=(y,weight,eigenvalues2,n,nMiss,y_hat))
   solvedGamma = np.exp(xpost)[0]
   return y_hat,solvedGamma 

def mfcc(input, nceps=13):
    """Compute Mel Frequency Cepstral Coefficients.

    Parameters
    ----------
    input: ndarray
        input spectrogram from which the coefficients are computed

    Returns
    -------
    ceps: ndarray
        Mel-cepstrum coefficients
    mspec: ndarray
        Log-spectrum in the mel-domain.

    Notes
    -----
    MFCC are computed as follows:
        * Pre-processing in time-domain (pre-emphasizing)
        * Compute the spectrum amplitude by windowing with a Hamming window
        * Filter the signal in the spectral domain with a triangular
        filter-bank, whose filters are approximatively linearly spaced on the
        mel scale, and have equal bandwith in the mel scale
        * Compute the DCT of the log-spectrum

    This is based on the talkbox module:
    http://pydoc.net/Python/scikits.talkbox/0.2.4.dev/scikits.talkbox.features.mfcc/

    References
    ----------
    .. [1] S.B. Davis and P. Mermelstein, "Comparison of parametric
           representations for monosyllabic word recognition in continuously
           spoken sentences", IEEE Trans. Acoustics. Speech, Signal Proc.
           ASSP-28 (4): 357-366, August 1980."""


    nfft = input.metadata.sampling_configuration.dft_length
    fs = input.metadata.sampling_configuration.fs
    over = input.metadata.sampling_configuration.window_length \
        - input.metadata.sampling_configuration.window_step

    #lowfreq = 400 / 3.
    lowfreq = 133.33
    #highfreq = 6855.4976
    linsc = 200/3.
    logsc = 1.0711703

    nlinfil = 13
    nlogfil = 27
    nlinfil + nlogfil

    fbank = trfbank(fs, nfft, lowfreq, linsc, logsc, nlinfil, nlogfil)[0]
    fbank = fbank.T[0:input.data.shape[0], :]

    mspec = np.log10(np.maximum(np.dot(fbank.T, input.data), 0.0000001)).T

    # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:, :nceps]

    return ceps

def get_mfcc(path):
    """Finds the MFCCs and FFTs of a WAVE file.

  Args:
    path: The path to a WAVE file.

  Returns:
    A tuple of two iterables, the FFTs and MFCCs of the frames of the
    WAVE file.
  """
    global COMP_FRAME_SIZE
    # Read the file, and determine its length in frames
    (sample, data) = utils.read_wave_from_file(path)
    total_frames = (data.size / sample) / COMP_FRAME_SIZE

    step = COMP_FRAME_SIZE * sample
    window = hamming(step)

    # Allocate space for the FFT decompositions of each frame of sound data
    fft_out = []
    mfcc_out = []

    # Loop invariant:
    #   0 <= frame_index <= total_frames
    #   results in an array (fft_out) of FFTs that correspond to the
    #    frames of the WAVE file
    filterbank_cache = {}
    frame_index = 0

    while frame_index + (1 - FRAME_OVERLAP_FACTOR) < total_frames:
        # Obtain the frame_indexth frame from the data
        frame = data[frame_index * step : (frame_index + 1) * step]

        # Generate the FFT of the frame windowed by the hamming window
        frame_fft = numpy.fft.rfft(frame * window, n=256)
        frame_fft[frame_fft == 0] = 0.000003
        nfft = len(frame_fft)

        # Compute the mel triangular filterbank or get a cached version
        fb_key = (sample, nfft)
        if fb_key in filterbank_cache:
            filterbank = filterbank_cache[fb_key]
        else:
            filterbank = triangular_filters(sample, nfft).T
            filterbank[filterbank == 0] = 0.00003
            filterbank_cache[fb_key] = filterbank

        # The power spectrum of the frame
        power_spectrum = numpy.abs(frame_fft)
        # Filtered by the mel filterbank
        mel_power_spectrum = numpy.log10(numpy.dot(power_spectrum, filterbank))
        # With the discrete cosine transform to find the cepstrum
        cepstrum = dct(mel_power_spectrum, type=2, norm="ortho", axis=-1)

        fft_out.append(frame_fft)
        mfcc_out.append(cepstrum[: int(len(cepstrum) * SIGNIFICANT_MFCC)])
        frame_index = frame_index + FRAME_OVERLAP_FACTOR

    return numpy.array(mfcc_out)

 def test_definition_matlab(self):
     """Test correspondance with matlab (orthornomal mode)."""
     for i in range(len(X)):
         x = np.array(X[i], dtype=self.rdt)
         yr = Y[i]
         y = dct(x, norm="ortho", type=2)
         self.assertTrue(y.dtype == self.rdt, "Output dtype is %s, expected %s" % (y.dtype, self.rdt))
         assert_array_almost_equal(y, yr, decimal=self.dec)

 def test_definition_matlab(self):
     # Test correspondance with matlab (orthornomal mode).
     for i in range(len(X)):
         x = np.array(X[i], dtype=self.rdt)
         yr = Y[i]
         y = dct(x, norm="ortho", type=2)
         assert_equal(y.dtype, self.rdt)
         assert_array_almost_equal(y, yr, decimal=self.dec)

 def test_definition_ortho(self):
     # Test orthornomal mode.
     for i in range(len(X)):
         x = np.array(X[i], dtype=self.rdt)
         dt = np.result_type(np.float32, self.rdt)
         y = dct(x, norm='ortho', type=4)
         y2 = naive_dct4(x, norm='ortho')
         assert_equal(y.dtype, dt)
         assert_array_almost_equal(y / np.max(y), y2 / np.max(y), decimal=self.dec)

 def test_definition(self):
     for i in FFTWDATA_SIZES:
         x, yr = fftw_ref(self.type, i, self.rdt)
         y = dct(x, type=self.type)
         self.assertTrue(y.dtype == self.rdt, "Output dtype is %s, expected %s" % (y.dtype, self.rdt))
         # XXX: we divide by np.max(y) because the tests fail otherwise. We
         # should really use something like assert_array_approx_equal. The
         # difference is due to fftw using a better algorithm w.r.t error
         # propagation compared to the ones from fftpack.
         assert_array_almost_equal(y / np.max(y), yr / np.max(y), decimal=self.dec, err_msg="Size %d failed" % i)

def stMFCC(X, fbank, nceps):
    """
    Computes the MFCCs of a frame, given the fft mag
    ARGUMENTS:
        X:        fft magnitude abs(FFT)
        fbank:    filter bank (see mfccInitFilterBanks)
    RETURN
        ceps:     MFCCs (13 element vector)
    Note:    MFCC calculation is, in general, taken from the scikits.talkbox library (MIT Licence),
    #    with a small number of modifications to make it more compact and suitable for the pyAudioAnalysis Lib
    """

    mspec = numpy.log10(numpy.dot(X, fbank.T)+eps)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:nceps]
    return ceps

def gcv(gamma_,y,weight,eigenvalues2,n,nMiss,y_hat_final):
  # a GCV function for the smoother
  from scipy.fftpack.realtransforms import dct,idct
  gamma = np.exp((gamma_))
  G = 1./(1+gamma*eigenvalues2)
  y0 = y.copy()
  e = 1e20
  while (e > 1e-10):
    y_hat = idct(G*dct(weight*weight*(y-y0)+y0,norm='ortho',type=2),norm='ortho',type=2)
    dy = y_hat - y0
    e = np.mean(dy*dy)
    y0 = y_hat
  y_hat_final[:] = y_hat
  d = weight*(y_hat-y)
  numerator = np.dot(d,d)/(n-nMiss)
  traceH = (1./(1 + gamma*eigenvalues2)).sum()
  denominator = (1 - traceH/n)**2
  return numerator/denominator

def mfcc(input, nwin=256, nfft=512, fs=16000, nceps=13):
    import numpy as numpy
    from scipy.io import loadmat
    from scipy.signal import lfilter, hamming
    from scipy.fftpack import fft
    from scipy.fftpack.realtransforms import dct
    over = nwin - 160
    prefac = 0.97

    #lowfreq = 400 / 3.
    lowfreq = 133.33

    linsc = 200/3.
    logsc = 1.0711703

    nlinfil = 13
    nlogfil = 27
    nfil = nlinfil + nlogfil

    w = hamming(nwin, sym=0)

    fbank = trfbank(fs, nfft, lowfreq, linsc, logsc, nlinfil, nlogfil)[0]
    if fbank<=0:
        fbank=0.0001
    #------------------
    # Compute the MFCC
    #------------------
    extract = preemp(input, prefac)
    framed = segment_axis(extract, nwin, over) * w

    # Compute the spectrum magnitude
    spec = numpy.abs(fft(framed, nfft, axis=-1))

    if spec<=0:
        spec=0.00001

    # Filter the spectrum through the triangle filterbank
    arr= numpy.dot(spec,fbank.T)                                         ##CHANGED CODE
    print "LOG ARRAy =",arr
    mspec=numpy.log10(arr)
    # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:, :nceps]

    return ceps

  def __init_mfcc(self,
      num_mel_bands = DEFAULT_MFCC_BANDS,
      num_mfcc = DEFAULT_NUM_MFCC_COEFFICIENTS,
      delta_N = DEFAULT_MFCC_DELTA_N):
    mel_bin_matrix, freqs = self.get_mel_binning_matrix(num_mel_bands)
    Pxx2 = np.dot(self.specgram.T, mel_bin_matrix)

    # Unlike the mlab implementation, we threshold and log our FFT magnitudes
    # before returning
    Pxx2[Pxx2 < 1e-10] = 1e-10
    Pxx2 = 10. * np.log10(Pxx2)
    Pxx2[Pxx2 <= 0.0] = 0.0

    # http://pydoc.net/Python/scikits.talkbox/0.2.4.dev/scikits.talkbox.features.mfcc/
    ceps = dct(Pxx2, type=2, norm='ortho', axis=-1)[:, :num_mfcc]
    ceps = np.flipud(ceps)

    deltas = np.zeros(ceps.shape)
    delta_deltas = np.zeros(ceps.shape)
    for cep_frame_i in xrange(len(ceps)):
      if cep_frame_i < delta_N:
        del_N = cep_frame_i
      elif cep_frame_i > len(ceps) - delta_N - 1:
        del_N = len(ceps) - cep_frame_i - 1
      else:
        del_N = delta_N
      if del_N == 0:
        continue
      deltas[cep_frame_i] = sum([n*(ceps[cep_frame_i + n] - ceps[cep_frame_i - n]) for n in xrange(1,del_N+1)]) / (2.0*sum([n**2 for n in xrange(1, del_N + 1)]))
    for cep_frame_i in xrange(len(deltas)):
      if cep_frame_i < delta_N:
        del_N = cep_frame_i
      elif cep_frame_i > len(ceps) - delta_N - 1:
        del_N = len(ceps) - cep_frame_i - 1
      else:
        del_N = delta_N
      if del_N == 0:
        continue
      delta_deltas[cep_frame_i] = sum([n*(deltas[cep_frame_i + n] - deltas[cep_frame_i - n]) for n in xrange(1,del_N+1)]) / (2.0*sum([n**2 for n in xrange(1, del_N + 1)]))
    ceps = np.fliplr(ceps.T[1:])
    deltas = np.fliplr(deltas.T[1:])
    delta_deltas = np.fliplr(delta_deltas.T[1:])
    return ceps, deltas, delta_deltas

def mfcc(input, nwin=256, nfft=512, fs=16000, nceps=13):
    """Compute Mel Frequency Cepstral Coefficients.

    Parameters
    ----------
    input: ndarray
        input from which the coefficients are computed

    Returns
    -------
    ceps: ndarray
        Mel-cepstrum coefficients
    mspec: ndarray
        Log-spectrum in the mel-domain.

    Notes
    -----
    MFCC are computed as follows:
        * Pre-processing in time-domain (pre-emphasizing)
        * Compute the spectrum amplitude by windowing with a Hamming window
        * Filter the signal in the spectral domain with a triangular
        filter-bank, whose filters are approximatively linearly spaced on the
        mel scale, and have equal bandwith in the mel scale
        * Compute the DCT of the log-spectrum

    References
    ----------
    .. [1] S.B. Davis and P. Mermelstein, "Comparison of parametric
           representations for monosyllabic word recognition in continuously
           spoken sentences", IEEE Trans. Acoustics. Speech, Signal Proc.
           ASSP-28 (4): 357-366, August 1980."""

    # Number of overlapping samples in each frame
    t_overlap = 10*10**(-3) # Time in seconds of overlapping between frames
    over = int(t_overlap*fs)
#     over = nwin - 160

    # Pre-emphasis factor (to take into account the -6dB/octave rolloff of the
    # radiation at the lips level)
    prefac = 0.97

    #lowfreq = 400 / 3.
    lowfreq = 133.33
    #highfreq = 6855.4976
    linsc = 200/3.
    logsc = 1.0711703

    nlinfil = 13
    nlogfil = 27
    nfil = nlinfil + nlogfil

    w = hamming(nwin, sym=0)
    
    [fbank, freqs] = trfbank(fs, nfft, lowfreq, linsc, logsc, nlinfil, nlogfil) # "fbank" is a nfil-by-nfft Numpy 2D array.
    
    '''
    # Visualizando o banco de filtros:
    plt.figure()
    nfiltros,lenfiltros = fbank.shape
    for i in range(nfiltros):
        plt.plot(range(lenfiltros),fbank[i,:])
    plt.axis([0, lenfiltros, 0, np.max(fbank)])
    plt.show()
    '''
     
    #------------------
    # Compute the MFCC
    #------------------
    extract = preemp(input, prefac)
    framed = segment_axis(extract, nwin, over) * w


    # Compute the spectrum magnitude
    spec = np.abs(fft(framed, nfft, axis=-1))
    # Filter the spectrum through the triangle filterbank
    mspec = np.log10(np.dot(spec, fbank.T))
    # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:, :nceps]
    nframes = ceps.shape[0]
    print 'nframes: ', nframes
	    
#     -----------------------------------------
#     Cepstrum mean subtraction
#     mean_along_frames = np.mean(mspec,axis=0) # Mean along the vertical dimension of the mel-spectrum
#     
#     mean_along_frames_stack = mean_along_frames
#     for i in range(nframes-1):
#         mean_along_frames_stack = np.vstack((mean_along_frames_stack, mean_along_frames))
#     ceps = ceps - mean_along_frames_stack[:,0:nceps]
    
    return ceps, mspec, spec

def mfcc(input, nwin=256, nfft=512, fs=16000, nceps=13):
    """Compute Mel Frequency Cepstral Coefficients.

    Parameters
    ----------
    input: ndarray
        input from which the coefficients are computed

    Returns
    -------
    ceps: ndarray
        Mel-cepstrum coefficients
    mspec: ndarray
        Log-spectrum in the mel-domain.

    Notes
    -----
    MFCC are computed as follows:
        * Pre-processing in time-domain (pre-emphasizing)
        * Compute the spectrum amplitude by windowing with a Hamming window
        * Filter the signal in the spectral domain with a triangular
        filter-bank, whose filters are approximatively linearly spaced on the
        mel scale, and have equal bandwith in the mel scale
        * Compute the DCT of the log-spectrum

    References
    ----------
    .. [1] S.B. Davis and P. Mermelstein, "Comparison of parametric
           representations for monosyllabic word recognition in continuously
           spoken sentences", IEEE Trans. Acoustics. Speech, Signal Proc.
           ASSP-28 (4): 357-366, August 1980."""

    # MFCC parameters: taken from auditory toolbox
    over = nwin - 160
    # Pre-emphasis factor (to take into account the -6dB/octave rolloff of the
    # radiation at the lips level)
    prefac = 0.97

    #lowfreq = 400 / 3.
    lowfreq = 133.33
    #highfreq = 6855.4976
    linsc = 200/3.
    logsc = 1.0711703

    nlinfil = 13
    nlogfil = 27
    nfil = nlinfil + nlogfil

    w = hamming(nwin, sym=0)

    fbank = trfbank(fs, nfft, lowfreq, linsc, logsc, nlinfil, nlogfil)[0]

    #------------------
    # Compute the MFCC
    #------------------
    extract = preemp(input, prefac)
    framed = segment_axis(extract, nwin, over) * w

    # Compute the spectrum magnitude
    spec = np.abs(fft(framed, nfft, axis=-1))
    # Filter the spectrum through the triangle filterbank
    mspec = np.log10(np.dot(spec, fbank.T))
    # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:, :nceps]

    return ceps, mspec, spec

def MFCC(X, fbank, nceps):
    mspec = np.log10(np.dot(X, fbank.T)+eps)
    ceps = dct(mspec, type=2, norm='ortho', axis=-1)[:nceps]
    return ceps   

 def test_dct_complex(self):
     y = dct(np.arange(5) * 1j)
     x = 1j * dct(np.arange(5))
     assert_array_almost_equal(x, y)