Get DeepFormants working again

* Minor syntax tweaks to make the code Python 3 compatible

* Fixes for various NumPy warnings/errors, either due to use of "float"
  where "int" is required, or domain errors on log functions

* Replaced the use of the obsolete Python-2-only scikits.talkbox
  library with a compatible LPC implementation from the Conch project

* Documentation update to indicate that an old version of "rnn" is
  required

* Invoke Lua scripts via "luajit" directly, instead of going through
  the "th" frontend (to reduce the dependency footprint)

README.md

@@ -24,9 +24,7 @@ Download the code. The code is based on signal processing package in Python call
 Dependencies:
 Run these lines in a terminal to install everything necessary for feature extraction.
 ```
-sudo apt-get install python-numpy python-scipy python-nose python-pip
+sudo apt-get install python3-numpy python3-scipy python3-nose
-
-sudo pip install scikits.talkbox 
 ```
 Next for the installation of Torch for loading the models run this.
 ```
@@ -37,30 +35,31 @@ cd ~/torch; bash install-deps;
 ./install.sh
 ```
 ```
-luarocks install rnn
+git clone https://github.com/Element-Research/rnn.git old-rnn
+cd old-rnn; luarocks make rocks/rnn-scm-1.rockspec
 ```
-The Estimation model can be downloaded here and because of size constraints the Tracking model can be abtained by download from this link
+The Estimation model can be downloaded here and because of size constraints the Tracking model can be obtained by download from this link:
-[tracking_model.mat] (https://drive.google.com/open?id=0Bxkc5_D0JjpiZWx4eTU1d0hsVXc)
+[tracking_model.mat](https://drive.google.com/open?id=0Bxkc5_D0JjpiZWx4eTU1d0hsVXc)
 
 ## How to use:
 
 For vowel formant estimation, call the main script in a terminal with the following inputs: wav file, formant output filename, and the vowel begin and end times:
 
 ```
-python formants.py data/Example.wav data/ExamplePredictions.csv --begin 1.2 --end 1.3
+python3 formants.py data/Example.wav data/ExamplePredictions.csv --begin 1.2 --end 1.3
 ```
 
 or the vowel begin and end times can be taken from a TextGrid file (here the name of the TextGrid is Example.TextGrid and the vowel is taken from a tier called "VOWEL"):
 
 ```
-python formants.py data/Example.wav data/examplePredictions.csv --textgrid_filename data/Example.TextGrid \
+python3 formants.py data/Example.wav data/examplePredictions.csv --textgrid_filename data/Example.TextGrid \
           --textgrid_tier VOWEL
 ```
 
 For formant tracking, just call the script with the wav file and output filename:
 
 ```
-python formants.py data/Example.wav data/ExamplePredictions.csv
+python3 formants.py data/Example.wav data/ExamplePredictions.csv
 ```
 
 

VTR_Results.py

@@ -0,0 +1,142 @@
+from keras.models import model_from_json
+import numpy as np
+import csv
+import math
+
+model = model_from_json(open('model.json').read())
+model.load_weights('weights.h5')
+data_dir = ""
+X_test = np.load(data_dir+'VTR_test_X.npy')
+Y = np.load(data_dir+'VTR_test_Y.npy')
+
+names = Y[:, :1]
+Y_test = Y[:,1:]
+predictions = []
+
+loss1 = 0.0
+loss2 = 0.0
+loss3 = 0.0
+loss4 = 0.0
+max_1 = 0.0
+max_2 = 0.0
+max_3 = 0.0
+max_4 = 0.0
+list_1 = []
+list_2 = []
+list_3 = []
+list_4 = []
+male = [0.0, 0.0, 0.0, 0.0, 0.0, [], [], [], []]
+female = [0.0, 0.0, 0.0, 0.0, 0.0, [], [], [], []]
+karma_list = [0, 0.0, 0.0, 0.0, 0.0]
+AVG_list = [0, 0.0, 0.0, 0.0, 0.0]
+
+y_hat = model.predict(X_test)
+for i in range(0,len(Y_test)):
+    l1 = np.abs(float(Y_test[i, 0]) - y_hat[i, 0])
+    l2 = np.abs(float(Y_test[i, 1]) - y_hat[i, 1])
+    l3 = np.abs(float(Y_test[i, 2]) - y_hat[i, 2])
+    l4 = np.abs(float(Y_test[i, 3]) - y_hat[i, 3])
+    pred = [names[i][0], float(Y_test[i, 0]), float(Y_test[i, 1]), float(Y_test[i, 2]), float(Y_test[i, 3])]
+
+    AVG_list[0] += 1
+    AVG_list[1] += float(Y_test[i, 0]) - y_hat[i, 0]
+    AVG_list[2] += float(Y_test[i, 1]) - y_hat[i, 1]
+    AVG_list[3] += float(Y_test[i, 2]) - y_hat[i, 2]
+    AVG_list[4] += float(Y_test[i, 3]) - y_hat[i, 3]
+
+    pred.extend([y_hat[i, 0], y_hat[i, 1], y_hat[i, 2], y_hat[i, 3]])
+
+    if names[i][0].split('_')[3][0] == 'f':
+        female[0] += 1
+        female[1] += l1
+        female[2] += l2
+        female[3] += l3
+        female[4] += l4
+        female[5].append(l1)
+        female[6].append(l2)
+        female[7].append(l3)
+        female[8].append(l4)
+    elif names[i][0].split('_')[3][0] == 'm':
+        male[0] += 1
+        male[1] += l1
+        male[2] += l2
+        male[3] += l3
+        male[4] += l4
+        male[5].append(l1)
+        male[6].append(l2)
+        male[7].append(l3)
+        male[8].append(l4)
+
+    predictions.append(pred)
+
+    list_1.append(l1)
+    list_2.append(l2)
+    list_3.append(l3)
+    list_4.append(l4)
+    max_1 = max(max_1,l1)
+    max_2 = max(max_2,l2)
+    max_3 = max(max_3,l3)
+    max_4 = max(max_4,l4)
+
+    loss1 += l1
+    loss2 += l2
+    loss3 += l3
+    loss4 += l4
+
+    karma_list[0] += 1
+    karma_list[1] += l1 * l1
+    karma_list[2] += l2 * l2
+    karma_list[3] += l3 * l3
+    karma_list[4] += l4 * l4
+loss1 /= len(Y_test)
+loss2 /= len(Y_test)
+loss3 /= len(Y_test)
+loss4 /= len(Y_test)
+total_loss = loss1+loss2+loss3+loss4
+total_loss /= 4.0
+print('standard deviation', round(np.std(list_1)*1000, 2), round(np.std(list_2)*1000, 2), round(np.std(list_3)*1000, 2), round(np.std(list_4)*1000, 2))
+print('median', round(np.median(list_1)*1000, 2), round(np.median(list_2)*1000, 2), round(np.median(list_3)*1000, 2), round(np.median(list_4)*1000, 2))
+print('max loss ', round(max_1*1000, 2), round(max_2*1000, 2), round(max_3*1000, 2), round(max_4*1000, 2))
+print('total loss ', round(total_loss*1000, 2))
+print('Real test score:', round(loss1*1000, 2), round(loss2*1000, 2), round(loss3*1000, 2), round(loss4*1000, 2))
+
+female[1] = round((female[1] / female[0])*1000, 2)
+female[2] = round((female[2] / female[0])*1000, 2)
+female[3] = round((female[3] / female[0])*1000, 2)
+female[4] = round((female[4] / female[0])*1000, 2)
+female[5] = round(np.std(female[5])*1000, 2)
+female[6] = round(np.std(female[6])*1000, 2)
+female[7] = round(np.std(female[7])*1000, 2)
+female[8] = round(np.std(female[8])*1000, 2)
+
+male[1] = round((male[1] / male[0])*1000, 2)
+male[2] = round((male[2] / male[0])*1000, 2)
+male[3] = round((male[3] / male[0])*1000, 2)
+male[4] = round((male[4] / male[0])*1000, 2)
+male[5] = round(np.std(male[5])*1000, 2)
+male[6] = round(np.std(male[6])*1000, 2)
+male[7] = round(np.std(male[7])*1000, 2)
+male[8] = round(np.std(male[8])*1000, 2)
+
+print("male: ", male)
+print("female: ", female)
+
+# karma
+
+karma_list[1] /= karma_list[0]
+karma_list[2] /= karma_list[0]
+karma_list[3] /= karma_list[0]
+karma_list[4] /= karma_list[0]
+print('root mean squared error ', round(math.sqrt(karma_list[1]) * 1000, 2), round(math.sqrt(karma_list[2]) * 1000, 2),
+                            round(math.sqrt(karma_list[3]) * 1000, 2), round(math.sqrt(karma_list[4]) * 1000, 2))
+
+AVG_list[1] /= AVG_list[0]
+AVG_list[2] /= AVG_list[0]
+AVG_list[3] /= AVG_list[0]
+AVG_list[4] /= AVG_list[0]
+print('AVG ', round(AVG_list[1] * 1000, 2), round(AVG_list[2] * 1000, 2), round(AVG_list[3] * 1000, 2), round(AVG_list[4] * 1000, 2))
+
+
+with open("results/VTR.csv", "wb") as f:
+    writer = csv.writer(f)
+    writer.writerows(predictions)

extract_features.py

@@ -9,9 +9,9 @@ from os.path import isfile, join
 import math
 from scipy.fftpack.realtransforms import dct
 from scipy.signal import lfilter, hamming
-from copy import deepcopy
 from scipy.fftpack import fft, ifft
-from scikits.talkbox.linpred import lpc
+#from scikits.talkbox.linpred import lpc  # obsolete
+from helpers.conch_lpc import lpc
 import shutil
 from helpers.utilities import *
 
@@ -26,7 +26,8 @@ def build_data(wav,begin=None,end=None):
     dstr = wav_in_file.readframes(N)
     data = np.fromstring(dstr, np.int16)
     if begin is not None and end is not None:
-        return data[begin*16000:end*16000]
+        #return data[begin*16000:end*16000] #numpy 1.11.0
+        return data[np.int(begin*16000):np.int(end*16000)] #numpy 1.14.0
     X = []
     l = len(data)
     for i in range(0, l-100, 160):
@@ -87,9 +88,9 @@ def periodogram(x, nfft=None, fs=1):
 
     pxx = np.abs(fft(x, nfft)) ** 2
     if nfft % 2 == 0:
-        pn = nfft / 2 + 1
+        pn = nfft // 2 + 1
     else:
-        pn = (nfft + 1 )/ 2
+        pn = (nfft + 1) // 2
 
     fgrid = np.linspace(0, fs * 0.5, pn)
     return pxx[:pn] / (n * fs), fgrid
@@ -136,9 +137,9 @@ def arspec(x, order, nfft=None, fs=1):
 
     # This is not enough to deal correctly with even/odd size
     if nfft % 2 == 0:
-        pn = nfft / 2 + 1
+        pn = nfft // 2 + 1
     else:
-        pn = (nfft + 1 )/ 2
+        pn = (nfft + 1) // 2
 
     px = 1 / np.fft.fft(a, nfft)[:pn]
     pxx = np.real(np.conj(px) * px)
@@ -199,7 +200,6 @@ def preemp(input, p):
 
 
 def arspecs(input_wav,order,Atal=False):
-    epsilon = 0.0000000001
     data = input_wav
     if Atal:
         ar = atal(data, order, 30)
@@ -210,8 +210,10 @@ def arspecs(input_wav,order,Atal=False):
         for k, l in zip(ars[0], ars[1]):
             ar.append(math.log(math.sqrt((k**2)+(l**2))))
         for val in range(0,len(ar)):
-            if ar[val] == 0.0:
+            if ar[val] < 0.0:
-                ar[val] = deepcopy(epsilon)
+                ar[val] = np.nan
+            elif ar[val] == 0.0:
+                ar[val] = epsilon
         mspec1 = np.log10(ar)
         # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)
         ar = dct(mspec1, type=2, norm='ortho', axis=-1)
@@ -220,10 +222,10 @@ def arspecs(input_wav,order,Atal=False):
 
 def specPS(input_wav,pitch):
         N = len(input_wav)
-        samps = N/pitch
+        samps = N // pitch
         if samps == 0:
             samps = 1
-        frames = N/samps
+        frames = N // samps
         data = input_wav[0:frames]
         specs = periodogram(data,nfft=4096)
         for i in range(1,int(samps)):
@@ -235,10 +237,11 @@ def specPS(input_wav,pitch):
             specs[0][s] /= float(samps)
         peri = []
         for k, l in zip(specs[0], specs[1]):
-            if k == 0 and l == 0:
+            m = math.sqrt((k ** 2) + (l ** 2))
-                peri.append(epsilon)
+            if m > 0: m = math.log(m)
-            else:
+            if m == 0: m = epsilon
-                peri.append(math.log(math.sqrt((k ** 2) + (l ** 2))))
+            elif m < 0: m = np.nan
+            peri.append(m)
         # Filter the spectrum through the triangle filterbank
         mspec = np.log10(peri)
         # Use the DCT to 'compress' the coefficients (spectrum -> cepstrum domain)

formants.py

@@ -9,19 +9,19 @@ import shutil
 
 def predict_from_times(wav_filename, preds_filename, begin, end):
     tmp_features_filename = tempfile._get_default_tempdir() + "/" + next(tempfile._get_candidate_names()) + ".txt"
-    print tmp_features_filename
+    print(tmp_features_filename)
 
     if begin > 0.0 or end > 0.0:
         features.create_features(wav_filename, tmp_features_filename, begin, end)
-        easy_call("th load_estimation_model.lua " + tmp_features_filename + ' ' + preds_filename)
+        easy_call("luajit load_estimation_model.lua " + tmp_features_filename + ' ' + preds_filename)
     else:
         features.create_features(wav_filename, tmp_features_filename)
-        easy_call("th load_tracking_model.lua " + tmp_features_filename + ' ' + preds_filename)
+        easy_call("luajit load_tracking_model.lua " + tmp_features_filename + ' ' + preds_filename)
 
 
 def predict_from_textgrid(wav_filename, preds_filename, textgrid_filename, textgrid_tier):
 
-    print wav_filename
+    print(wav_filename)
 
     if os.path.exists(preds_filename):
         os.remove(preds_filename)

formants.sh

@@ -4,12 +4,12 @@ if [ $# -eq 2 ]
 then
     tempfile=`mktemp -t txt`
     python extract_features.py $1 $tempfile
-    th load_estimation_model.lua $tempfile $2
+    luajit load_estimation_model.lua $tempfile $2
 elif [ $# -eq 4 ]
 then
     tempfile=`mktemp -t txt`
     python extract_features.py $1 $tempfile --begin $3 --end $4
-    th load_estimation_model.lua $tempfile $2
+    luajit load_estimation_model.lua $tempfile $2
 else
     echo "$0 wav_filename pred_csv_filename [begin_time end_time]"
 fi

helpers/conch_lpc.py

@@ -0,0 +1,286 @@
+# This file has been copied (with minor changes) from Michael
+# McAuliffe's Conch project, to provide a compatible replacement
+# implementation of the lpc() function from the obsolete Python-2-only
+# scikits.talkbox library.
+#
+# Conch repository: https://github.com/mmcauliffe/Conch-sounds
+# Source: https://github.com/mmcauliffe/Conch-sounds/blob/master/conch/analysis/formants/lpc.py
+
+# Copyright (c) 2015 Michael McAuliffe
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+
+#import librosa
+import numpy as np
+import scipy as sp
+from scipy.signal import lfilter
+
+from scipy.fftpack import fft, ifft
+from scipy.signal import gaussian
+
+#from ..helper import nextpow2
+#from ..functions import BaseAnalysisFunction
+
+# Source: https://github.com/mmcauliffe/Conch-sounds/blob/master/conch/analysis/helper.py
+def nextpow2(x):
+    """Return the first integer N such that 2**N >= abs(x)"""
+    return np.ceil(np.log2(np.abs(x)))
+
+def lpc_ref(signal, order):
+    """Compute the Linear Prediction Coefficients.
+
+    Return the order + 1 LPC coefficients for the signal. c = lpc(x, k) will
+    find the k+1 coefficients of a k order linear filter:
+
+      xp[n] = -c[1] * x[n-2] - ... - c[k-1] * x[n-k-1]
+
+    Such as the sum of the squared-error e[i] = xp[i] - x[i] is minimized.
+
+    Parameters
+    ----------
+    signal: array_like
+        input signal
+    order : int
+        LPC order (the output will have order + 1 items)
+
+    Notes
+    ----
+    This is just for reference, as it is using the direct inversion of the
+    toeplitz matrix, which is really slow"""
+    if signal.ndim > 1:
+        raise ValueError("Array of rank > 1 not supported yet")
+    if order > signal.size:
+        raise ValueError("Input signal must have a lenght >= lpc order")
+
+    if order > 0:
+        p = order + 1
+        r = np.zeros(p, 'float32')
+        # Number of non zero values in autocorrelation one needs for p LPC
+        # coefficients
+        nx = np.min([p, signal.size])
+        x = np.correlate(signal, signal, 'full')
+        r[:nx] = x[signal.size - 1:signal.size + order]
+        phi = np.dot(sp.linalg.inv(sp.linalg.toeplitz(r[:-1])), -r[1:])
+        return np.concatenate(([1.], phi))
+    else:
+        return np.ones(1, dtype='float32')
+
+
+# @jit
+def levinson_1d(r, order):
+    """Levinson-Durbin recursion, to efficiently solve symmetric linear systems
+    with toeplitz structure.
+
+    Parameters
+    ---------
+    r : array-like
+        input array to invert (since the matrix is symmetric Toeplitz, the
+        corresponding pxp matrix is defined by p items only). Generally the
+        autocorrelation of the signal for linear prediction coefficients
+        estimation. The first item must be a non zero real.
+
+    Notes
+    ----
+    This implementation is in python, hence unsuitable for any serious
+    computation. Use it as educational and reference purpose only.
+
+    Levinson is a well-known algorithm to solve the Hermitian toeplitz
+    equation:
+
+                       _          _
+        -R[1] = R[0]   R[1]   ... R[p-1]    a[1]
+         :      :      :          :      *  :
+         :      :      :          _      *  :
+        -R[p] = R[p-1] R[p-2] ... R[0]      a[p]
+                       _
+    with respect to a (  is the complex conjugate). Using the special symmetry
+    in the matrix, the inversion can be done in O(p^2) instead of O(p^3).
+    """
+    r = np.atleast_1d(r)
+    if r.ndim > 1:
+        raise ValueError("Only rank 1 are supported for now.")
+
+    n = r.size
+    if n < 1:
+        raise ValueError("Cannot operate on empty array !")
+    elif order > n - 1:
+        raise ValueError("Order should be <= size-1")
+
+    if not np.isreal(r[0]):
+        raise ValueError("First item of input must be real.")
+    elif not np.isfinite(1 / r[0]):
+        raise ValueError("First item should be != 0")
+
+    # Estimated coefficients
+    a = np.empty(order + 1, 'float32')
+    # temporary array
+    t = np.empty(order + 1, 'float32')
+    # Reflection coefficients
+    k = np.empty(order, 'float32')
+
+    a[0] = 1.
+    e = r[0]
+
+    for i in range(1, order + 1):
+        acc = r[i]
+        for j in range(1, i):
+            acc += a[j] * r[i - j]
+        k[i - 1] = -acc / e
+        a[i] = k[i - 1]
+
+        for j in range(order):
+            t[j] = a[j]
+
+        for j in range(1, i):
+            a[j] += k[i - 1] * np.conj(t[i - j])
+
+        e *= 1 - k[i - 1] * np.conj(k[i - 1])
+
+    return a, e, k
+
+
+# @jit
+def _acorr_last_axis(x, nfft, maxlag):
+    a = np.real(ifft(np.abs(fft(x, n=nfft) ** 2)))
+    return a[..., :maxlag + 1] / x.shape[-1]
+
+
+# @jit
+def acorr_lpc(x, axis=-1):
+    """Compute autocorrelation of x along the given axis.
+
+    This compute the biased autocorrelation estimator (divided by the size of
+    input signal)
+
+    Notes
+    -----
+        The reason why we do not use acorr directly is for speed issue."""
+    if not np.isrealobj(x):
+        raise ValueError("Complex input not supported yet")
+
+    maxlag = x.shape[axis]
+    nfft = int(2 ** nextpow2(2 * maxlag - 1))
+
+    if axis != -1:
+        x = np.swapaxes(x, -1, axis)
+    a = _acorr_last_axis(x, nfft, maxlag)
+    if axis != -1:
+        a = np.swapaxes(a, -1, axis)
+    return a
+
+
+# @jit
+def lpc(signal, order, axis=-1):
+    """Compute the Linear Prediction Coefficients.
+
+    Return the order + 1 LPC coefficients for the signal. c = lpc(x, k) will
+    find the k+1 coefficients of a k order linear filter:
+
+      xp[n] = -c[1] * x[n-2] - ... - c[k-1] * x[n-k-1]
+
+    Such as the sum of the squared-error e[i] = xp[i] - x[i] is minimized.
+
+    Parameters
+    ----------
+    signal: array_like
+        input signal
+    order : int
+        LPC order (the output will have order + 1 items)
+
+    Returns
+    -------
+    a : array-like
+        the solution of the inversion.
+    e : array-like
+        the prediction error.
+    k : array-like
+        reflection coefficients.
+
+    Notes
+    -----
+    This uses Levinson-Durbin recursion for the autocorrelation matrix
+    inversion, and fft for the autocorrelation computation.
+
+    For small order, particularly if order << signal size, direct computation
+    of the autocorrelation is faster: use levinson and correlate in this case."""
+    n = signal.shape[axis]
+    if order > n:
+        raise ValueError("Input signal must have length >= order")
+
+    r = acorr_lpc(signal, axis)
+    return levinson_1d(r, order)
+
+
+def process_frame(X, window, num_formants, new_sr):
+    X = X * window
+    A, e, k = lpc(X, num_formants * 2)
+
+    rts = np.roots(A)
+    rts = rts[np.where(np.imag(rts) >= 0)]
+    angz = np.arctan2(np.imag(rts), np.real(rts))
+    frqs = angz * (new_sr / (2 * np.pi))
+    frq_inds = np.argsort(frqs)
+    frqs = frqs[frq_inds]
+    bw = -1 / 2 * (new_sr / (2 * np.pi)) * np.log(np.abs(rts[frq_inds]))
+    return frqs, bw
+
+
+def lpc_formants(signal, sr, num_formants, max_freq, time_step,
+                 win_len, window_shape='gaussian'):
+    output = {}
+    new_sr = 2 * max_freq
+    alpha = np.exp(-2 * np.pi * 50 * (1 / new_sr))
+    proc = lfilter([1., -alpha], 1, signal)
+    if sr > new_sr:
+        proc = librosa.resample(proc, sr, new_sr)
+    nperseg = int(win_len * new_sr)
+    nperstep = int(time_step * new_sr)
+    if window_shape == 'gaussian':
+        window = gaussian(nperseg + 2, 0.45 * (nperseg - 1) / 2)[1:nperseg + 1]
+    else:
+        window = np.hanning(nperseg + 2)[1:nperseg + 1]
+    indices = np.arange(int(nperseg / 2), proc.shape[0] - int(nperseg / 2) + 1, nperstep)
+    num_frames = len(indices)
+    for i in range(num_frames):
+        if nperseg % 2 != 0:
+            X = proc[indices[i] - int(nperseg / 2):indices[i] + int(nperseg / 2) + 1]
+        else:
+            X = proc[indices[i] - int(nperseg / 2):indices[i] + int(nperseg / 2)]
+        frqs, bw = process_frame(X, window, num_formants, new_sr)
+        formants = []
+        for j, f in enumerate(frqs):
+            if f < 50:
+                continue
+            if f > max_freq - 50:
+                continue
+            formants.append((np.asscalar(f), np.asscalar(bw[j])))
+        missing = num_formants - len(formants)
+        if missing:
+            formants += [(None, None)] * missing
+        output[indices[i] / new_sr] = formants
+    return output
+
+
+#class FormantTrackFunction(BaseAnalysisFunction):
+#    def __init__(self, num_formants=5, max_frequency=5000,
+#                 time_step=0.01, window_length=0.025, window_shape='gaussian'):
+#        super(FormantTrackFunction, self).__init__()
+#        self.arguments = [num_formants, max_frequency, time_step, window_length, window_shape]
+#        self._function = lpc_formants
+#        self.requires_file = False

training_for_git/batch.lua

@@ -1,104 +0,0 @@
-require 'torch'   -- torch
-require 'optim'
-require 'nn'      -- provides a normalization operator
-local train_file_path = 'train.th7' 
-local test_file_path = 'test.th7'
-local train_data = torch.load(train_file_path)
-local test_data = torch.load(test_file_path)
-local Y = train_data[{{},{2,5}}]
-local X = train_data[{{},{6,-1}}]
-local test_labels = test_data[{{},{2,5}}]
-local test_X = test_data[{{},{6,-1}}]
-local batch_size = 30
-epochs = 3
-
-model = nn.Sequential()                 -- define the container
-ninputs = 350; noutputs = 4 ; nhiddens1 = 1024; nhiddens2 = 512; nhiddens3 = 256
-model:add(nn.Linear(ninputs,nhiddens1))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens1,nhiddens2))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens2,nhiddens3))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens3,noutputs))
-criterion = nn.AbsCriterion()--MSECriterion()
-x, dl_dx = model:getParameters()
-sgd_params = {
-   learningRate = 0.01,
-   learningRateDecay = 1e-08,
-   weightDecay = 0,
-   momentum = 0
-}
-
-function train(X,Y)
-	   
-   current_loss = 0
-   for batch = 1,(#train_data)[1], batch_size do
-
-      local inputs = {}
-      local targets = {} 
-      local x_start = batch
-      local x_end = math.min(batch + batch_size-1, (#train_data)[1])
-      for i = x_start,x_end do    
-      local target = Y[i]    
-      local input = X[i]  
-      table.insert(inputs, input)
-      table.insert(targets, target)
-      end
-	local feval = function(x_new)
-	   if x ~= x_new then
-	      x:copy(x_new)
-	   end
-	   dl_dx:zero()
-           local f=0
-	for i = 1, #inputs do
-	   local loss_x = criterion:forward(model:forward(inputs[i]), targets[i])
-	   model:backward(inputs[i], criterion:backward(model.output, targets[i]))
-	f = f+loss_x
-	end
-	   return f/#inputs, dl_dx:div(#inputs)
-	end
-      _,fs = optim.adagrad(feval,x,sgd_params)
-      current_loss = current_loss + fs[1]
-   end
-   current_loss = current_loss/( (#train_data)[1]/batch_size)
-   print('train loss = ' .. current_loss)
-   return current_loss
-end
-
-time = sys.clock()
-local cumm_loss = 0.
-for j = 1, epochs do
-    print(j)
-    cumm_loss = train( X, Y )
-print( 'Final loss = ' .. cumm_loss )
-if j%10 == 0 then
-print('id  approx   text')
-local loss1 = 0.0
-local loss2 = 0.0
-local loss3 = 0.0
-local loss4 = 0.0
-for i = 1,(#test_data)[1] do
-   local myPrediction = model:forward(test_X[i])
-   loss1 = loss1+math.abs(myPrediction[1] - test_labels[i][1])
-   loss2 = loss2+math.abs(myPrediction[2] - test_labels[i][2])
-   loss3 = loss3+math.abs(myPrediction[3] - test_labels[i][3])
-   loss4 = loss4+math.abs(myPrediction[4] - test_labels[i][4])
-end
-
-loss1 = loss1/(#test_data)[1]
-loss2 = loss2/(#test_data)[1]
-loss3 = loss3/(#test_data)[1]
-loss4 = loss4/(#test_data)[1]
-end
-end
-
-
--- time taken
-time = sys.clock() - time
-print( "Time per epoch = " .. (time / epochs) .. '[s]')
-
-
-
-print(loss1,loss2,loss3,loss4)
-torch.save('estimation_model.dat',model)

training_for_git/gpu_rec.lua

@@ -1,130 +0,0 @@
-require 'rnn'
-require 'optim'
-
-batchSize = 30
-rho = 10
-hiddenSize = 512
-hiddenSize1 = 256
-inputSize = 400
-outputSize = 3
-epochs = 100
-xStart = 6
-yStart = 2
-yEnd = 4
-
-
-local train_file_path = 'recurrent_train.th7' 
-local train_data = torch.load(train_file_path)
-local Y = train_data[{{},{yStart,yEnd}}]
-local X = train_data[{{},{xStart,-1}}]
-seriesSize = (#train_data)[1]
-print(seriesSize)
-local test_file_path = 'recurrent_test.th7'
-local test_data = torch.load(test_file_path)
-local test_labels = test_data[{{},{yStart,yEnd}}]
-local test_X = test_data[{{},{xStart,-1}}]
-
-model = nn.Sequential()
-model:add(nn.Sequencer(nn.FastLSTM(inputSize, hiddenSize, rho)))
-model:add(nn.Sequencer(nn.FastLSTM(hiddenSize, hiddenSize1, rho)))
-model:add(nn.Sequencer(nn.Linear(hiddenSize1, outputSize)))
-
-criterion = nn.SequencerCriterion(nn.AbsCriterion())
-
--- dummy dataset (task predict the next item)
---dataset = torch.randn(seriesSize, inputSize)
-
--- define the index of the batch elements
-offsets = {}
-for i= 1, batchSize do
-   table.insert(offsets, i)--math.ceil(math.random() * batchSize))
-end
-offsets = torch.LongTensor(offsets)
-
-function nextBatch()
-   local inputs, targets = {}, {}
-   for step = 1, rho do
-      --get a batch of inputs
-      table.insert(inputs, X:index(1, offsets))
-      -- shift of one batch indexes
-      offsets:add(1)
-      for j=1,batchSize do
-         if offsets[j] > seriesSize then
-            offsets[j] = 1
-         end
-      end
-      -- a batch of targets
-      table.insert(targets, Y[{{},{1,3}}]:index(1,offsets))
-   end
-   return inputs, targets
-end
-
--- get weights and loss wrt weights from the model
-x, dl_dx = model:getParameters()
-
-feval = function(x_new)
-	-- copy the weight if are changed
-	if x ~= x_new then
-		x:copy(x_new)
-	end
-
-	-- select a training batch
-	local inputs, targets = nextBatch()
-
-	-- reset gradients (gradients are always accumulated, to accommodate
-	-- batch methods)
-	dl_dx:zero()
-
-	-- evaluate the loss function and its derivative wrt x, given mini batch
-	local prediction = model:forward(inputs)
-	local loss_x = criterion:forward(prediction, targets)
-	model:backward(inputs, criterion:backward(prediction, targets))
-	
-	return loss_x, dl_dx
-end
-
-sgd_params = {
-   learningRate = 0.01,
-   learningRateDecay = 1e-08,
-   weightDecay = 0,
-   momentum = 0
-}
-
-time = sys.clock()
-for j = 1, epochs do
-   -- train a mini_batch of batchSize in parallel
-	 _, fs = optim.adagrad(feval,x, sgd_params)
-	print('error for iteration ' .. sgd_params.evalCounter  .. ' is ' .. fs[1])
-
-end
-
-
-print('id  approx   text')
-local loss1 = 0.0
-local loss2 = 0.0
-local loss3 = 0.0
-local loss4 = 0.0
-for i = 1,(#test_data)[1], 1 do
-   local inputs = {}
-   for step = 1, 1 do
-      --get a batch of inputs
-      table.insert(inputs, test_X[i])
-   end
-   local myPrediction = model:forward(inputs)
-   loss1 = loss1+math.abs(myPrediction[1][1] - test_labels[i][1])
-   loss2 = loss2+math.abs(myPrediction[1][2] - test_labels[i][2])
-   loss3 = loss3+math.abs(myPrediction[1][3] - test_labels[i][3])
-   --loss4 = loss4+math.abs(myPrediction[4] - test_labels[i][4])
-end
-
-loss1 = loss1/(#test_data)[1]
-loss2 = loss2/(#test_data)[1]
-loss3 = loss3/(#test_data)[1]
---loss4 = loss4/(#test_data)[1]
-
--- time taken
-time = sys.clock() - time
-print( "Time per epoch = " .. (time / epochs) .. '[s]')
-
-print(loss1,loss2,loss3,loss4)
-torch.save('recurrent.dat',model)

training_for_git/no_validation.lua

@@ -1,129 +0,0 @@
-require 'torch'   -- torch
-require 'optim'
-require 'nn'      -- provides a normalization operator
-local train_file_path = 'train.th7' 
-local test_file_path = 'test.th7'
-local train_data = torch.load(train_file_path)
-local test_data = torch.load(test_file_path)
-local train_labels = train_data[{{},{2,5}}]
-local train_X = train_data[{{},{6,-1}}]
-local test_labels = test_data[{{},{2,5}}]
-local test_X = test_data[{{},{6,-1}}]
-local batch_size = 30
-model = nn.Sequential()                 -- define the container
-ninputs = 350; noutputs = 4 ; nhiddens1 = 1024; nhiddens2 = 512; nhiddens3 = 256
---model:add(nn.Linear(ninputs, noutputs)) -- define the only module
-model:add(nn.Linear(ninputs,nhiddens1))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens1,nhiddens2))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens2,nhiddens3))
-model:add(nn.Sigmoid())
-model:add(nn.Linear(nhiddens3,noutputs))
-criterion = nn.AbsCriterion()--MSECriterion()
-x, dl_dx = model:getParameters()
-
-feval = function(x_new)
-   if x ~= x_new then
-      x:copy(x_new)
-   end
-   -- select a new training sample
-   _nidx_ = (_nidx_ or 0) + 1
-   if _nidx_ > (#train_data)[1] then _nidx_ = 1 end
-   --local sample = data[_nidx_]
-   local target = train_labels[_nidx_]      -- this funny looking syntax allows
-   local inputs = train_X[_nidx_]    -- slicing of arrays.
-   -- reset gradients (gradients are always accumulated, to accommodate 
-   -- batch methods)
-   dl_dx:zero()
-   -- evaluate the loss function and its derivative wrt x, for that sample
-   --print(inputs)
-   --print(target)
-   for i=1, 350 do
-  if type(inputs[i]) ~= 'number' then
-  print(i)
-  print(inputs[i])
-  print(type(inputs[i])) end
-  end
-   --io.write("continue with this operation (y/n)?")
-   --answer=io.read()
-   local loss_x = criterion:forward(model:forward(inputs), target)
-   model:backward(inputs, criterion:backward(model.output, target))
-   -- return loss(x) and dloss/dx
-   return loss_x, dl_dx
-end
--- Given the function above, we can now easily train the model using SGD.
--- For that, we need to define four key parameters:
---   + a learning rate: the size of the step taken at each stochastic 
---     estimate of the gradient
---   + a weight decay, to regularize the solution (L2 regularization)
---   + a momentum term, to average steps over time
---   + a learning rate decay, to let the algorithm converge more precisely
-sgd_params = {
-   learningRate = 0.01,
-   learningRateDecay = 1e-08,
-   weightDecay = 0,
-   momentum = 0
-}
--- We're now good to go... all we have left to do is run over the dataset
--- for a certain number of iterations, and perform a stochastic update 
--- at each iteration. The number of iterations is found empirically here,
--- but should typically be determinined using cross-validation.
--- we cycle 1e4 times over our training data
-for i = 1,1 do
-   print(i)
-   -- this variable is used to estimate the average loss
-   current_loss = 0
-   -- an epoch is a full loop over our training data
-   for i = 1,(#train_data)[1] do
-      -- optim contains several optimization algorithms. 
-      -- All of these algorithms assume the same parameters:
-      --   + a closure that computes the loss, and its gradient wrt to x, 
-      --     given a point x
-      --   + a point x
-      --   + some parameters, which are algorithm-specific      
-      _,fs = optim.adagrad(feval,x,sgd_params)
-      -- Functions in optim all return two things:
-      --   + the new x, found by the optimization method (here SGD)
-      --   + the value of the loss functions at all points that were used by
-      --     the algorithm. SGD only estimates the function once, so
-      --     that list just contains one value.
-      current_loss = current_loss + fs[1]
-   end
-   -- report average error on epoch
-   current_loss = current_loss / (#train_data)[1]
-   print('train loss = ' .. current_loss)
-   
-end
-----------------------------------------------------------------------
--- 5. Test the trained model.
-
--- Now that the model is trained, one can test it by evaluating it
--- on new samples.
-
--- The text solves the model exactly using matrix techniques and determines
--- that 
---   corn = 31.98 + 0.65 * fertilizer + 1.11 * insecticides
-
--- We compare our approximate results with the text's results.
-
-print('id  approx   text')
-local loss1 = 0.0
-local loss2 = 0.0
-local loss3 = 0.0
-local loss4 = 0.0
-for i = 1,(#test_data)[1] do
-   local myPrediction = model:forward(test_X[i])
-   loss1 = loss1+math.abs(myPrediction[1] - test_labels[i][1])
-   loss2 = loss2+math.abs(myPrediction[2] - test_labels[i][2])
-   loss3 = loss3+math.abs(myPrediction[3] - test_labels[i][3])
-   loss4 = loss4+math.abs(myPrediction[4] - test_labels[i][4])
-end
-
-loss1 = loss1/(#test_data)[1]
-loss2 = loss2/(#test_data)[1]
-loss3 = loss3/(#test_data)[1]
-loss4 = loss4/(#test_data)[1]
-
-print(loss1,loss2,loss3,loss4)
-torch.save('save.dat',model)

training_for_git/recurrent.lua

@@ -1,109 +0,0 @@
-require 'rnn'
-require 'optim'
-
-function range(from, to, step)
-  step = step or 1
-  return function(_, lastvalue)
-    local nextvalue = lastvalue + step
-    if step > 0 and nextvalue <= to or step < 0 and nextvalue >= to or
-       step == 0
-    then
-      return nextvalue
-    end
-  end, nil, from - step
-end
-
-local train_file_path = 'recurrent_train.th7' 
-local test_file_path = 'recurrent_test.th7'
-local train_data = torch.load(train_file_path)
-local test_data = torch.load(test_file_path)
-local Y = train_data[{{},{2,5}}]
-local X = train_data[{{},{6,-1}}]
-local test_labels = test_data[{{},{2,5}}]
-local test_X = test_data[{{},{6,-1}}]
-
-batchSize = 5
-rho = 10
-hiddenSize1 = 1024
-hiddenSize2 = 512
-hiddenSize3 = 256
-inputSize = 1
-outputSize = 1
-seriesSize = 100
-
-model = nn.Sequential()
-model:add(nn.Sequencer(nn.FastLSTM(inputSize, hiddenSize2, rho)))
-model:add(nn.Sequencer(nn.FastLSTM(hiddenSize2, hiddenSize3, rho)))
---model:add(nn.Sequencer(nn.Linear(hiddenSize2, hiddenSize3, rho)))
---model:add(nn.Sequencer(nn.Sigmoid()))
-model:add(nn.Sequencer(nn.Linear(hiddenSize3, outputSize)))
-
-criterion = nn.SequencerCriterion(nn.MSECriterion())
-
--- dummy dataset (task predict the next item)
---dataset = torch.randn(seriesSize, inputSize)
-
--- define the index of the batch elements
-offsets = {}
-for i= 1, batchSize do
-   table.insert(offsets,i)
-end
-offsets = torch.LongTensor(offsets)
-print(offsets)
-function nextBatch()
-   local inputs, targets = {}, {}
-   for step = 1, rho do
-      --get a batch of inputs
-      table.insert(inputs, X:index(1, offsets))
-      -- shift of one batch indexes
-      offsets:add(1)
-      for j=1,batchSize do
-         if offsets[j] > seriesSize then
-            offsets[j] = 1
-         end
-      end
-      -- a batch of targets
-      table.insert(targets, Y:index(1,offsets))
-   end
-   return inputs, targets
-end
-
--- get weights and loss wrt weights from the model
-x, dl_dx = model:getParameters()
-
-feval = function(x_new)
-  -- copy the weight if are changed
-  if x ~= x_new then
-    x:copy(x_new)
-  end
-
-  -- select a training batch
-  local inputs, targets = nextBatch()
-
-  -- reset gradients (gradients are always accumulated, to accommodate
-  -- batch methods)
-  dl_dx:zero()
-
-  -- evaluate the loss function and its derivative wrt x, given mini batch
-  local prediction = model:forward(inputs)
-  local loss_x = criterion:forward(prediction, targets)
-  model:backward(inputs, criterion:backward(prediction, targets))
-
-  return loss_x, dl_dx
-end
-
-sgd_params = {
-   learningRate = 0.01,
-   learningRateDecay = 1e-08,
-   weightDecay = 0,
-   momentum = 0
-}
-
-for i = 1, 2 do
-   -- train a mini_batch of batchSize in parallel
-   _, fs = optim.adagrad(feval,x, sgd_params)
-
-  if sgd_params.evalCounter % 100 == 0 then
-    print('error for iteration ' .. sgd_params.evalCounter  .. ' is ' .. fs[1] / rho)
-  end
-end

training_for_git/recurrent_training.lua

@@ -1,144 +0,0 @@
-require 'rnn'
-require 'optim'
-
-batchSize = 30
-rho = 20
-hiddenSize = 512
-hiddenSize1 = 256
-inputSize = 400
-outputSize = 4
-epochs = 10000
-xStart = 6
-yStart = 2
-yEnd = 5
-
-
-local train_file_path = 'recurrent_train.th7' 
-local train_data = torch.load(train_file_path)
-local Y = train_data[{{},{yStart,yEnd}}]
-local X = train_data[{{},{xStart,-1}}]
-local place = train_data[{{},{1}}]
-seriesSize = (#train_data)[1]
-print(seriesSize)
-local test_file_path = 'recurrent_test.th7'
-local test_data = torch.load(test_file_path)
-local test_labels = test_data[{{},{yStart,yEnd}}]
-local test_X = test_data[{{},{xStart,-1}}]
-
-model = nn.Sequential()
-model:add(nn.Sequencer(nn.FastLSTM(inputSize, hiddenSize, rho)))
-model:add(nn.Sequencer(nn.FastLSTM(hiddenSize, hiddenSize1, rho)))
-model:add(nn.Sequencer(nn.Linear(hiddenSize1, outputSize)))
-
-criterion = nn.SequencerCriterion(nn.AbsCriterion())
---local method = 'xavier'
---local model_new = require('weight-init')(model, method)
-
--- define the index of the batch elements
-offsets = {}
-function offset_(seed)
-offsets = {}
-math.randomseed(seed)
-for i= 1, batchSize do
-   table.insert(offsets, math.ceil(math.random() * batchSize))
-end
-offsets = torch.LongTensor(offsets)
-end
-function nextBatch()
-   local inputs, targets = {}, {}
-   local nums = {}
-   for step = 1, rho do
-      --get a batch of inputs
-      table.insert(inputs, X:index(1, offsets))
-      -- shift of one batch indexes
-      offsets:add(1)
-      for j=1,batchSize do
-         if offsets[j] > seriesSize then
-            offsets[j] = 1
-         end
-      end
-      -- a batch of targets
-      table.insert(targets, Y[{{},{1,4}}]:index(1,offsets))
-      table.insert(nums,place:index(1,offsets))
-   end
-   return inputs, targets
-end
-
--- get weights and loss wrt weights from the model
-x, dl_dx = model:getParameters()
-
-feval = function(x_new)
-	-- copy the weight if are changed
-	if x ~= x_new then
-		x:copy(x_new)
-	end
-
-	-- select a training batch
-	local inputs, targets = nextBatch()
-
-	-- reset gradients (gradients are always accumulated, to accommodate
-	-- batch methods)
-	dl_dx:zero()
-
-	-- evaluate the loss function and its derivative wrt x, given mini batch
-	local prediction = model:forward(inputs)
-	local loss_x = criterion:forward(prediction, targets)
-	model:backward(inputs, criterion:backward(prediction, targets))
-	
-	return loss_x, dl_dx
-end
-
-adagrad_params = {
-   learningRate = 0.01,
-   learningRateDecay = 1e-08,
-   weightDecay = 0,
-   momentum = 0
-}
-seed = 1
-offset_(seed)
-time = sys.clock()
-for j = 1, epochs do
-	if j%1000 == 0 then 
-		seed = seed + 1
-		offset_(seed)
-	end
-   -- train a mini_batch of batchSize in parallel
-	 _, fs = optim.adagrad(feval,x, adagrad_params)
-	print('error for iteration ' .. adagrad_params.evalCounter  .. ' is ' .. fs[1]/rho)
-
-end
-
-
-print('id  approx   text')
-local loss1 = 0.0
-local loss2 = 0.0
-local loss3 = 0.0
-local loss4 = 0.0
-predict_batch = 100
-for i = 1,(#test_data)[1], predict_batch do
-   local inputs = {}
-   for step = 0, predict_batch-1 do
-      --get a batch of inputs
-      table.insert(inputs, test_X[i+step])
-   end
-   local myPrediction = model:forward(inputs)
-   for step = 1, predict_batch do
-   loss1 = loss1+math.abs(myPrediction[step][1] - test_labels[i+step-1][1])
-   loss2 = loss2+math.abs(myPrediction[step][2] - test_labels[i+step-1][2])
-   loss3 = loss3+math.abs(myPrediction[step][3] - test_labels[i+step-1][3])
-   loss4 = loss4+math.abs(myPrediction[4] - test_labels[i][4])
-   end
-
-end
-
-loss1 = loss1/(#test_data)[1]
-loss2 = loss2/(#test_data)[1]
-loss3 = loss3/(#test_data)[1]
-loss4 = loss4/(#test_data)[1]
-
--- time taken
-time = sys.clock() - time
-print( "Time per epoch = " .. (time / epochs) .. '[s]')
-
-print(loss1,loss2,loss3,loss4)
-torch.save('recurrent3.dat',model)