Skip to content

Instantly share code, notes, and snippets.

@shunsukeaihara
Last active August 29, 2015 14:27
Show Gist options
  • Save shunsukeaihara/f5672e0a8eae5fb55fa3 to your computer and use it in GitHub Desktop.
Save shunsukeaihara/f5672e0a8eae5fb55fa3 to your computer and use it in GitHub Desktop.
cython translation of the estimating alpha parameter for mel-cepstrum in go (https://bitbucket.org/happyalu/mcep_alpha_calc)
# -*- coding: utf-8 -*-
# distutils: language = c++
import numpy as np
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free
from libc.math cimport log, M_PI, sin, cos, atan
from libc.float cimport DBL_MAX
@cython.boundscheck(False)
@cython.wraparound(False)
def estimate_alpha(int sampfreq, double start=0, double end=1.0, double step=0.001, int size=1000):
"""
estimate alpha parameter for mel-cepstrum from sampling rate.
original imprementation https://bitbucket.org/happyalu/mcep_alpha_calc
"""
cdef double *melscale_vector = create_melscale_vector(sampfreq, size)
cdef double *warping_vector = <double *>malloc(size * sizeof(double))
cdef double i,
cdef double min_dist, best_alpha, dist
best_alpha = 0.0
min_dist = DBL_MAX
for a from start <= a <= end by step:
calc_warping_vector(warping_vector, a, size)
dist = rms_distance_like(melscale_vector, warping_vector, size)
if dist < min_dist:
min_dist = dist
best_alpha = a
free(melscale_vector)
free(warping_vector)
return best_alpha
cdef double* create_melscale_vector(int sampfreq, int size) nogil:
cdef double *vec = <double *>malloc(size * sizeof(double))
cdef double step = (sampfreq / 2.0) / size
cdef int i
for i in range(size):
vec[i] = (1000.0 / log(2)) * log(1.0+((step * i)/1000.0))
cdef double last = vec[size - 1]
for i in range(size):
vec[i] /= last
return vec
cdef void calc_warping_vector(double *vec, double alpha, int size) nogil:
cdef double step = M_PI / size
cdef int i
cdef double omega, warpfreq, num, den
for i in range(size):
omega = step * i
num = (1 - alpha * alpha) * sin(omega)
den = (1 + alpha * alpha) * cos(omega) - 2 * alpha
warpfreq = atan(num / den)
if warpfreq < 0:
warpfreq += M_PI
vec[i] = warpfreq
cdef double last = vec[size - 1]
for i in range(size):
vec[i] /= last
cdef double rms_distance_like(double *a, double *b, int size) nogil:
cdef int i
cdef double s = 0.0
for i in range(size):
s += (a[i] - b[i]) ** 2
return s
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment