cython translation of the estimating alpha parameter for mel-cepstrum in go (https://bitbucket.org/happyalu/mcep_alpha_calc)

estimate_alpha.pyx

# -*- 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