1-convert-ej-f.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -*- coding: utf-8 -*-\n",
    "\n",
    "\"\"\"\n",
    "This program makes learning ev-gmm.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# __future__ module make compatible python2 and python3\n",
    "from __future__ import division, print_function\n",
    "\n",
    "# basic modules\n",
    "import os\n",
    "import os.path\n",
    "import time\n",
    "\n",
    "# for warning ignore\n",
    "import warnings\n",
    "#warning.filterwarnings('ignore')\n",
    "\n",
    "# for file system manupulation \n",
    "from shutil import rmtree \n",
    "import glob\n",
    "import argparse\n",
    "\n",
    "# for save object\n",
    "import pickle\n",
    "\n",
    "# for make glaph\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_style(\"whitegrid\")\n",
    "plt.rcParams['figure.figsize'] = (16, 5)\n",
    "import librosa.display\n",
    "\n",
    "# for scientific computing\n",
    "import numpy as np\n",
    "from numpy.linalg import norm \n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.mixture import GMM # GMM class cannot use after sklearn 0.20.0\n",
    "import sklearn.mixture\n",
    "#from sklearn.mixture.gaussian_mixture import _compute_precision_cholesky\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import scipy.sparse\n",
    "from scipy.signal import firwin, lfilter\n",
    "\n",
    "# for display audio controler\n",
    "from IPython.display import Audio\n",
    "\n",
    "# for manuplate audio data\n",
    "import soundfile as sf\n",
    "import pyworld as pw\n",
    "import pysptk\n",
    "from dtw import dtw\n",
    "from fastdtw import fastdtw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class WORLD(object):\n",
    "    \"\"\"\n",
    "    WORLD based speech analyzer and synthezer.\n",
    "    \n",
    "    Ref : https://github.com/k2kobayashi/sprocket/\n",
    "    \"\"\"\n",
    "    def __init__(self, fs=16000, fftl=1024, shiftms=5.0, minf0=40.0, maxf0=500.0):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "        fs : int \n",
    "            Sampling frequency \n",
    "        fftl : int\n",
    "            FFT length\n",
    "        shiftms : float\n",
    "            Shift length [ms]\n",
    "        minf0 : float\n",
    "            Floor in F0 estimation\n",
    "        maxf0 : float\n",
    "            Seli in F0 estimation\n",
    "        \"\"\"\n",
    "        \n",
    "        self.fs = fs\n",
    "        self.fftl = fftl\n",
    "        self.shiftms = shiftms\n",
    "        self.minf0 = minf0\n",
    "        self.maxf0 = maxf0\n",
    "        \n",
    "    def analyze(self, x):\n",
    "        \"\"\"\n",
    "        Analyze acoustic featueres.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        x : array, shape(`T`)\n",
    "            monoral speech signal in time domain\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        sp : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            Spectral envelope sequence\n",
    "        ap : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            aperiodicity sequence\n",
    "        \"\"\"\n",
    "        \n",
    "        f0, time_axis = pw.harvest(x, self.fs, f0_floor=self.minf0,\n",
    "                                   f0_ceil=self.maxf0, frame_period=self.shiftms)\n",
    "        sp = pw.cheaptrick(x, f0, time_axis, self.fs, fft_size=self.fftl)\n",
    "        ap = pw.d4c(x, f0, time_axis, self.fs, fft_size=self.fftl)\n",
    "        \n",
    "        assert sp.shape == ap.shape\n",
    "        \n",
    "        return f0, sp, ap\n",
    "    \n",
    "    def analyze_f0(self, x):\n",
    "        \"\"\"\n",
    "        Analyze f0.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        x : array, shape(`T`)\n",
    "            monoral speech signal in time domain\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        \"\"\"\n",
    "\n",
    "        f0, time_axis = pw.harvest(x, self.fs, f0_floor=self.minf0,\n",
    "                                   f0_ceil=self.maxf0, frame_period=self.shiftms)\n",
    "        \n",
    "        assert f0.shape == x.shape()\n",
    "        \n",
    "        return f0\n",
    "    \n",
    "    def synthesis(self, f0, sp, ap):\n",
    "        \"\"\"\n",
    "        Re-synthesizes a speech waveform from acoustic featueres.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        sp : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            Spectral envelope sequence\n",
    "        ap : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            aperiodicity sequence\n",
    "        \"\"\"\n",
    "\n",
    "        return pw.synthesize(f0, sp, ap, self.fs, frame_period=self.shiftms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FeatureExtractor(object):\n",
    "    \"\"\"\n",
    "    Analyze acoustic features from a waveform.\n",
    "    \n",
    "    This class may have several types of estimeter like WORLD or STRAIGHT.\n",
    "    Default type is WORLD.\n",
    "    \n",
    "    Ref : https://github.com/k2kobayashi/sprocket/\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, analyzer='world', fs=16000, fftl=1024, \n",
    "                 shiftms=5.0, minf0=50.0, maxf0=500.0):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "        analyzer : str\n",
    "            Analyzer\n",
    "        fs : int \n",
    "            Sampling frequency \n",
    "        fftl : int\n",
    "            FFT length\n",
    "        shiftms : float\n",
    "            Shift length [ms]\n",
    "        minf0 : float\n",
    "            Floor in F0 estimation\n",
    "        maxf0 : float\n",
    "            Seli in F0 estimation\n",
    "        \"\"\"\n",
    "        \n",
    "        self.analyzer = analyzer\n",
    "        self.fs = fs\n",
    "        self.fftl = fftl\n",
    "        self.shiftms = shiftms\n",
    "        self.minf0 = minf0\n",
    "        self.maxf0 = maxf0\n",
    "    \n",
    "        if self.analyzer == 'world':\n",
    "            self.analyzer = WORLD(fs=self.fs, fftl=self.fftl, \n",
    "                                  minf0=self.minf0, maxf0=self.maxf0, shiftms=self.shiftms)\n",
    "        else:\n",
    "            raise('Analyzer Error : not support type, see FeatureExtractor class.')\n",
    "        \n",
    "        self._f0 = None\n",
    "        self._sp = None\n",
    "        self._ap = None\n",
    "        \n",
    "    def analyze(self, x):\n",
    "        \"\"\"\n",
    "        Analyze acoustic featueres.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        x : array, shape(`T`)\n",
    "            monoral speech signal in time domain\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        sp : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            Spectral envelope sequence\n",
    "        ap : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            aperiodicity sequence\n",
    "        \"\"\"\n",
    "        \n",
    "        self.x = np.array(x, dtype=np.float)\n",
    "        self._f0, self._sp, self._ap = self.analyzer.analyze(self.x)\n",
    "        \n",
    "        # check f0 < 0\n",
    "        self._f0[self._f0 < 0] = 0\n",
    "        \n",
    "        if np.sum(self._f0) == 0.0:\n",
    "            print(\"Warning : F0 values are all zero.\")\n",
    "        \n",
    "        return self._f0, self._sp, self._ap\n",
    "    \n",
    "    def analyze_f0(self, x):\n",
    "        \"\"\"\n",
    "        Analyze f0.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        x : array, shape(`T`)\n",
    "            monoral speech signal in time domain\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        \"\"\"\n",
    "\n",
    "        self.x = np.array(x, dtype=np.float)\n",
    "        self._f0 = self.analyzer.analyze_f0(self.x)\n",
    "\n",
    "        # check f0 < 0\n",
    "        self._f0[self._f0 < 0] = 0\n",
    "        \n",
    "        if np.sum(self._f0) == 0.0:\n",
    "            print(\"Warning : F0 values are all zero.\")\n",
    "        \n",
    "        return self._f0\n",
    "    \n",
    "    def mcep(self, dim=24, alpha=0.42):\n",
    "        \"\"\"\n",
    "        Convert mel-cepstrum sequence from spectral envelope.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        dim : int\n",
    "            mel-cepstrum dimension\n",
    "        alpha : float\n",
    "            parameter of all-path filter\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        mcep : array, shape(`T`, `dim + 1`)\n",
    "            mel-cepstrum sequence\n",
    "        \"\"\"        \n",
    "        \n",
    "        self._analyzed_check()\n",
    "        \n",
    "        return pysptk.sp2mc(self._sp, dim, alpha)\n",
    "    \n",
    "    def codeap(self):\n",
    "        \"\"\"\n",
    "        \"\"\"\n",
    "        self._analyzed_check()\n",
    "        \n",
    "        return pw.code_aperiodicity(self._ap, self.fs)\n",
    "    \n",
    "    def npow(self):\n",
    "        \"\"\"\n",
    "        Normalized power sequence from spectral envelope.\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        npow : vector, shape(`T`, `1`)\n",
    "            Normalized power sequence of the given waveform\n",
    "        \"\"\"\n",
    "        \n",
    "        self._analyzed_check()\n",
    "        \n",
    "        npow = np.apply_along_axis(self._spvec2pow, 1, self._sp)\n",
    "        \n",
    "        meanpow = np.mean(npow)\n",
    "        npow = 10.0 * np.log10(npow / meanpow)\n",
    "        \n",
    "        return npow\n",
    "    \n",
    "    def _spvec2pow(self, specvec):\n",
    "        \"\"\"\n",
    "        \"\"\"\n",
    "        fftl2 = len(specvec) - 1\n",
    "        fftl = fftl2 * 2\n",
    "        \n",
    "        power = specvec[0] + specvec[fftl2]\n",
    "        for k in range(1, fftl2):\n",
    "            power += 2.0 * specvec[k]\n",
    "        power /= fftl\n",
    "        \n",
    "        return power\n",
    "        \n",
    "    def _analyzed_check(self):\n",
    "        if self._f0 is None and self._sp is None and self._ap is None:\n",
    "            raise('Call FeatureExtractor.analyze() before this method.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Synthesizer(object):\n",
    "    \"\"\"\n",
    "    Synthesize a waveform from acoustic features.\n",
    "    \n",
    "    Ref : https://github.com/k2kobayashi/sprocket/\n",
    "    \"\"\"\n",
    "    def __init__(self, fs=16000, fftl=1024, shiftms=5.0):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "        fs : int \n",
    "            Sampling frequency \n",
    "        fftl : int\n",
    "            FFT length\n",
    "        shiftms : float\n",
    "            Shift length [ms]\n",
    "        \"\"\"\n",
    "        \n",
    "        self.fs = fs\n",
    "        self.fftl = fftl\n",
    "        self.shiftms = shiftms\n",
    "    \n",
    "    def synthesis(self, f0, mcep, ap, rmcep=None, alpha=0.42):\n",
    "        \"\"\"\n",
    "        Re-synthesizes a speech waveform from acoustic featueres.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        mcep : array, shape(`T`, `dim`)\n",
    "            mel-cepstrum sequence\n",
    "        ap : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            aperiodicity sequence\n",
    "        rmcep : array, shape(`T`, `dim`)\n",
    "            array of reference mel-cepstrum sequence\n",
    "        alpha : float\n",
    "            parameter of all-path filter\n",
    "            \n",
    "        Returns\n",
    "        ----------\n",
    "        wav : array,\n",
    "            syntesized waveform\n",
    "        \"\"\"\n",
    "        \n",
    "        if rmcep is not None:\n",
    "            # power modification\n",
    "            mcep = mod_power(mcep, rmcep, alpha=alpha)\n",
    "        \n",
    "        sp = pysptk.mc2sp(mcep, alpha, self.fftl)\n",
    "        wav = pw.synthesize(f0, sp, ap, self.fs, frame_period=self.shiftms)\n",
    "        \n",
    "        return wav\n",
    "        \n",
    "    def synthesis_diff(self, x, diffmcep, rmcep=None, alpha=0.42):\n",
    "        \"\"\"\n",
    "        Re-synthesizes a speech waveform from acoustic featueres.\n",
    "        filtering with a differential mel-cepstrum.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        x : array, shape(`samples`)\n",
    "            array of waveform sequence\n",
    "        diffmcep : array, shape(`T`, `dim`)\n",
    "            array of differential mel-cepstrum sequence\n",
    "        rmcep : array, shape(`T`, `dim`)\n",
    "            array of reference mel-cepstrum sequence\n",
    "        alpha : float\n",
    "            parameter of all-path filter\n",
    "            \n",
    "        Returns\n",
    "        ----------\n",
    "        wav : array,\n",
    "            syntesized waveform\n",
    "        \"\"\"      \n",
    "        \n",
    "        x = x.astype(np.float64)\n",
    "        dim = diffmcep.shape[1] - 1 \n",
    "        shiftl = int(self.fs / 1000 * self.shiftms)\n",
    "        \n",
    "        if rmcep is not None:\n",
    "            # power modification\n",
    "            diffmcep = mod_power(rmcep + diffmcep, rmcep, alpha=alpha) - rmcep        \n",
    "        \n",
    "        # mc2b = transform mel-cepstrum to MLSA digital filter coefficients.\n",
    "        b = np.apply_along_axis(pysptk.mc2b, 1, diffmcep, alpha)\n",
    "        \n",
    "        mlsa_fil = pysptk.synthesis.Synthesizer(pysptk.synthesis.MLSADF(dim, alpha=alpha),\n",
    "                                                shiftl)\n",
    "        wav = mlsa_fil.synthesis(x, b)\n",
    "        \n",
    "        return wav\n",
    "    \n",
    "    def synthesis_sp(self, f0, sp, ap):\n",
    "        \"\"\"\n",
    "        Re-synthesizes a speech waveform from acoustic featueres.\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        f0 : array, shape(`T`)\n",
    "            F0 sequence\n",
    "        spc : array, shape(`T`, `dim`)\n",
    "            mel-cepstrum sequence\n",
    "        ap : array, shape(`T`, `fftl / 2 + 1`)\n",
    "            aperiodicity sequence\n",
    "            \n",
    "        Returns\n",
    "        ----------\n",
    "        wav : array,\n",
    "            syntesized waveform\n",
    "        \"\"\"      \n",
    "        \n",
    "        wav = pw.synthesize(f0, sp, ap, self.fs, frame_period=self.shiftms)\n",
    "        \n",
    "        return wav\n",
    "    \n",
    "def mod_power(cvmcep, rmcep, alpha=0.42, irlen=256):\n",
    "    \"\"\"\n",
    "    power modification based on inpuulse responce\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    cvmcep : array, shape(`T`, `dim`)\n",
    "        array of converted mel-cepstrum\n",
    "    rmcep : arraym shape(`T`, `dim`)\n",
    "        array of reference mel-cepstrum\n",
    "    alpha : float\n",
    "        parameter of all-path filter\n",
    "    irlen : int\n",
    "        Length for IIR filter\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    modified_cvmcep : array, shape(`T`, `dim`)\n",
    "        array of power modified converted mel-cepstrum\n",
    "    \"\"\"\n",
    "\n",
    "    if rmcep.shape != cvmcep.shape:\n",
    "        raise ValueError(\n",
    "            \"The shape of the converted and reference mel-cepstrum are different : {} / {}.format(cvmcep.shape, rmcep.shape)\"\n",
    "        )\n",
    "\n",
    "    # mc2e = Compute energy from mel-cepstrum. e-option\n",
    "    cv_e = pysptk.mc2e(cvmcep, alpha=alpha, irlen=irlen)\n",
    "    r_e = pysptk.mc2e(rmcep, alpha=alpha, irlen=irlen)\n",
    "\n",
    "    dpow = np.log(r_e / cv_e) / 2\n",
    "\n",
    "    modified_cvmcep = np.copy(cvmcep)\n",
    "    modified_cvmcep[:, 0] += dpow\n",
    "\n",
    "    return modified_cvmcep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def util methods\n",
    "def melcd(array1, array2):\n",
    "    \"\"\"\n",
    "    calculate mel-cepstrum distortion\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    array1, array2 : array, shape(`T`, `dim`) or shape(`dim`)\n",
    "        Array of original and target.\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    mcd : scala, number > 0\n",
    "        Scala of mel-cepstrum distoriton\n",
    "    \"\"\"\n",
    "    if array1.shape != array2.shape:\n",
    "        raise ValueError(\n",
    "            \"The shape of both array are different : {} / {}.format(array1.shape,array2.shape)\"\n",
    "        )    \n",
    "   \n",
    "    if array1.ndim == 2:\n",
    "        diff = array1 - array2\n",
    "        mcd = 10.0 / np.log(10) * np.mean(np.sqrt(2.0 * np.sum(diff ** 2, axis=1)))\n",
    "    elif array1.ndim == 1:\n",
    "        diff = array1 - array2\n",
    "        mcd = 10.0 / np.log(10) * np.sqrt(2.0 * np.sum(diff ** 2))\n",
    "    else:\n",
    "        raise ValueError(\"Dimension mismatch.\")\n",
    "        \n",
    "    return mcd\n",
    "\n",
    "def delta(data, win=[-1.0, 1.0, 0]):\n",
    "    \"\"\"\n",
    "    calculate delta component\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array, shape(`T`, `dim`)\n",
    "        Array of static matrix sequence.\n",
    "    win : array, shape(`3`)\n",
    "        The shape of window matrix.\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    delta : array, shape(`T`, `dim`)\n",
    "        Array of delta matrix sequence.\n",
    "    \"\"\"\n",
    "    \n",
    "    if data.ndim == 1:\n",
    "        # change vector into 1d-array\n",
    "        T = len(data)\n",
    "        dim = data.ndim\n",
    "        data = data.reshape(T, dim)\n",
    "    else:\n",
    "        T, dim = data.shape\n",
    "    \n",
    "    win = np.array(win, dtype=np.float64)\n",
    "    delta = np.zeros((T, dim))\n",
    "    \n",
    "    delta[0] = win[0] * data[0] + win[1] * data[1]\n",
    "    delta[-1] = win[0] * data[-2] + win[1] * data[-1]\n",
    "    \n",
    "    for i in range(len(win)):\n",
    "        delta[1:T - 1] += win[i] * delta[i:T - 2 + i]\n",
    "    \n",
    "    return delta\n",
    "\n",
    "def static_delta(data, win=[-1.0, 1.0, 0]):\n",
    "    \"\"\"\n",
    "    calculate static and delta component\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array, shape(`T`, `dim`)\n",
    "        Array of static matrix sequence.\n",
    "    win : array, shape(`3`)\n",
    "        The shape of window matrix.\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    sddata : array, shape(`T`, `dim * 2`)\n",
    "        Array of static and delta matrix sequence.\n",
    "    \"\"\"\n",
    "    \n",
    "    sddata = np.c_[data, delta(data, win)]\n",
    "    \n",
    "    assert sddata.shape[1] == data.shape[1] * 2\n",
    "    \n",
    "    return sddata\n",
    "\n",
    "def construct_static_and_delta_matrix(T, D, win=[-1.0, 1.0, 0]):\n",
    "    \"\"\"\n",
    "    calculate static and delta transformation matrix\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    T : scala, `T`\n",
    "        Scala of time length\n",
    "    D : scala, `D`\n",
    "        Scala of the number of dimension.\n",
    "    win : array, shape(`3`)\n",
    "        The shape of window matrix.\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    W : array, shape(`2 * D * T`, `D * T`)\n",
    "        Array of static and delta transformation matrix.\n",
    "    \"\"\"\n",
    "    \n",
    "    static = [0, 1, 0]\n",
    "    delta = win\n",
    "    assert len(static) == len(delta)\n",
    "    \n",
    "    # generate full W\n",
    "    DT = D * T\n",
    "    ones = np.ones(DT)\n",
    "    row = np.arange(2 * DT).reshape(2 * T, D) # generate serial numbers\n",
    "    static_row = row[::2] # [1,2,3,4,5] => [1,3,5]\n",
    "    delta_row = row[1::2] # [1,2,3,4,5] => [2,4]\n",
    "    col = np.arange(DT)\n",
    "    \n",
    "    data = np.array([ones * static[0], ones * static[1],\n",
    "                      ones * static[2], ones * delta[0],\n",
    "                      ones * delta[1], ones * delta[2]]).flatten()\n",
    "    row = np.array([[static_row] * 3, [delta_row] * 3]).flatten()\n",
    "    col = np.array([[col - D, col, col + D] * 2]).flatten()\n",
    "\n",
    "    # remove component at first and end frame\n",
    "    valid_idx = np.logical_not(np.logical_or(col < 0, col >= DT))\n",
    "    \n",
    "    W = scipy.sparse.csr_matrix(\n",
    "        (data[valid_idx], (row[valid_idx], col[valid_idx])), shape=(2 * DT, DT))\n",
    "    W.eliminate_zeros()\n",
    "    \n",
    "    return W\n",
    "    \n",
    "def extfrm(data, npow, power_threshold=-20):\n",
    "    \"\"\"\n",
    "    Extract frame over the power threshold\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array, shape(`T`, `dim`)\n",
    "        array of input data\n",
    "    npow : array, shape(`T`)\n",
    "        vector of normalized power sequence\n",
    "    threshold : scala\n",
    "        scala of power threshold [dB]\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    data : array, shape(`T_ext`, `dim`)\n",
    "        remaining data after extracting frame\n",
    "        `T_ext` <= `T`\n",
    "    \"\"\"\n",
    "    T = data.shape[0]\n",
    "    if T != len(npow):\n",
    "        raise(\"Length of two vectors is different.\")\n",
    "        \n",
    "    valid_index = np.where(npow > power_threshold)\n",
    "    extdata = data[valid_index]\n",
    "    assert extdata.shape[0] <= T\n",
    "    \n",
    "    return extdata\n",
    "\n",
    "def estimate_twf(orgdata, tardata, distance='melcd', fast=True, otflag=None):\n",
    "    \"\"\"\n",
    "    time warping function estimator\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    orgdata : array, shape(`T_org`, `dim`)\n",
    "        array of source feature\n",
    "    tardata : array, shape(`T_tar`, `dim`)\n",
    "        array of target feature\n",
    "    distance : str \n",
    "        distance function\n",
    "    fast : bool\n",
    "        use fastdtw instead of dtw\n",
    "    otflag : str\n",
    "        Alignment into the length of specification\n",
    "        'org' : alignment into original length\n",
    "        'tar' : alignment into target length   \n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    twf : array, shape(`2`, `T`)\n",
    "        time warping function between original and target\n",
    "    \"\"\"\n",
    "    \n",
    "    if distance == 'melcd':\n",
    "        def distance_func(x, y): return melcd(x, y)\n",
    "    else:\n",
    "        raise ValueError('this distance method is not support.')\n",
    "    \n",
    "    if fast:\n",
    "        _, path = fastdtw(orgdata, tardata, dist=distance_func)\n",
    "        twf = np.array(path).T\n",
    "    else:\n",
    "        _, _, _, twf = dtw(orgdata, tardata, distance_func) \n",
    "        \n",
    "    if otflag is not None:\n",
    "        twf = modify_twf(twf, otflag=otflag)\n",
    "    \n",
    "    return twf\n",
    "\n",
    "def align_data(org_data, tar_data, twf):\n",
    "    \"\"\"\n",
    "    get aligned joint feature vector\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    org_data : array, shape(`T_org`, `dim_org`)\n",
    "        Acoustic feature vector of original speaker\n",
    "    tar_data : array, shape(`T_tar`, `dim_tar`)\n",
    "        Acoustic feature vector of target speaker\n",
    "    twf : array, shape(`2`, `T`)\n",
    "        time warping function between original and target\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    jdata : array, shape(`T_new`, `dim_org + dim_tar`)\n",
    "        Joint feature vector between source and target\n",
    "    \"\"\"\n",
    "    \n",
    "    jdata = np.c_[org_data[twf[0]], tar_data[twf[1]]]\n",
    "    return jdata\n",
    "\n",
    "def modify_twf(twf, otflag=None):\n",
    "    \"\"\"\n",
    "    align specified length\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    twf : array, shape(`2`, `T`)\n",
    "        time warping function between original and target\n",
    "    otflag : str\n",
    "        Alignment into the length of specification\n",
    "        'org' : alignment into original length\n",
    "        'tar' : alignment into target length   \n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    mod_twf : array, shape(`2`, `T_new`)\n",
    "        time warping function of modified alignment\n",
    "    \"\"\"\n",
    "    \n",
    "    if otflag == 'org':\n",
    "        of, indice = np.unique(twf[0], return_index=True)\n",
    "        mod_twf = np.c_[of, twf[1][indice]].T\n",
    "    elif otflag == 'tar':\n",
    "        tf, indice = np.unique(twf[1], return_index=True)\n",
    "        mod_twf = np.c_[twf[0][indice], tf].T \n",
    "    \n",
    "    return mod_twf\n",
    "\n",
    "def low_cut_filter(x, fs, cutoff=70):\n",
    "    \"\"\"\n",
    "    low cut filter\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array, shape('samples')\n",
    "        waveform sequence\n",
    "    fs : array, int\n",
    "        Sampling frequency\n",
    "    cutoff : float\n",
    "        cutoff frequency of low cut filter\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    lct_x : array, shape('samples')\n",
    "        Low cut filtered waveform sequence\n",
    "    \"\"\"\n",
    "    \n",
    "    nyquist = fs // 2\n",
    "    norm_cutoff = cutoff / nyquist\n",
    "    \n",
    "    # low cut filter\n",
    "    fil = firwin(255, norm_cutoff, pass_zero=False)\n",
    "    lct_x = lfilter(fil, 1, x)\n",
    "    \n",
    "    return lct_x\n",
    "\n",
    "def extsddata(data, npow, power_threshold=-20):\n",
    "    \"\"\"\n",
    "    get power extract static and delta feature vector\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    data : array, shape(`T`, `dim`)\n",
    "        acoustic feature vector\n",
    "    npow : array, shape(`T`)\n",
    "        normalized power vector\n",
    "    power_threshold : float\n",
    "        power threshold\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    extsddata : array, shape(`T_new`, `dim * 2`)\n",
    "        silence remove static and delta feature vector\n",
    "    \"\"\"\n",
    "    \n",
    "    extsddata = extfrm(static_delta(data), npow, power_threshold=power_threshold)\n",
    "    return extsddata\n",
    "\n",
    "def transform_jnt(array_list):\n",
    "    num_files = len(array_list)\n",
    "    for i in range(num_files):\n",
    "        if i == 0:\n",
    "            jnt = array_list[i]\n",
    "        else:\n",
    "            jnt = np.r_[jnt, array_list[i]]\n",
    "    return jnt\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class F0statistics(object):\n",
    "    \"\"\"\n",
    "    Estimate F0 statistics and convert F0\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        pass\n",
    "    \n",
    "    def estimate(self, f0list):\n",
    "        \"\"\"\n",
    "        estimate F0 statistics from list of f0\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        f0list : list, shape(`f0num`)\n",
    "            List of several F0 sequence\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        f0stats : array, shape(`[mean, std]`)\n",
    "            values of mean and standard deviation for log f0\n",
    "        \"\"\"\n",
    "        \n",
    "        n_files = len(f0list)\n",
    "        assert n_files != 0\n",
    "\n",
    "        for i in range(n_files):\n",
    "            f0 = f0list[i]\n",
    "            nonzero_indices = np.nonzero(f0)\n",
    "            if i == 0:\n",
    "                f0s = np.log(f0[nonzero_indices])\n",
    "            else:\n",
    "                f0s = np.r_[f0s, np.log(f0[nonzero_indices])]\n",
    "        \n",
    "        f0stats = np.array([np.mean(f0s), np.std(f0s)])\n",
    "        \n",
    "        return f0stats\n",
    "\n",
    "    def convert(self, f0, orgf0stats, tarf0stats):\n",
    "        \"\"\"\n",
    "        convert F0 based on F0 statistics\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        f0 : array, shape(`T`, `1`)\n",
    "            array of F0 sequence\n",
    "        orgf0stats : array, shape(`[mean, std]`)\n",
    "            vectors of mean and standard deviation of log f0 for original speaker\n",
    "        tarf0stats : array, shape(`[mean, std]`)\n",
    "            vectors of mean and standard deviation of log f0 for target speaker\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        cvf0 : array, shape(`T`, `1`)\n",
    "            array of converted F0 sequence\n",
    "        \"\"\"\n",
    "        \n",
    "        # get length and dimension\n",
    "        T = len(f0)\n",
    "        \n",
    "        # perform f0 conversion\n",
    "        cvf0 = np.zeros(T)\n",
    "        \n",
    "        nonzero_indices  = f0 > 0\n",
    "        cvf0[nonzero_indices] = np.exp((tarf0stats[1] / orgf0stats[1]) * (np.log(f0[nonzero_indices]) - orgf0stats[0]) + tarf0stats[0])\n",
    "        \n",
    "        return cvf0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GV(object):\n",
    "    \"\"\"\n",
    "    Estimate statistics and perform postfilter based on the GV statistics.\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        pass\n",
    "    \n",
    "    def estimate(self, datalist):\n",
    "        \"\"\"\n",
    "        estimate GV statistics from list of data\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        datalist : list, shape(`num_data`)\n",
    "            List of several data ([T, dim]) sequence\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        gvstats : array, shape(`2`, `dim`)\n",
    "            array of mean and standard deviation for GV\n",
    "        \"\"\"\n",
    "        \n",
    "        n_files = len(datalist)\n",
    "        assert n_files != 0\n",
    "        \n",
    "        var = []\n",
    "        for i in range(n_files):\n",
    "            data = datalist[i]\n",
    "            var.append(np.var(data, axis=0))\n",
    "            \n",
    "        # calculate vm and vv\n",
    "        vm = np.mean(np.array(var), axis=0)\n",
    "        vv = np.var(np.array(var), axis=0)\n",
    "        gvstats = np.r_[vm, vv]\n",
    "        gvstats = gvstats.reshape(2, len(vm))\n",
    "        \n",
    "        return gvstats\n",
    "\n",
    "    def postfilter(self, data, gvstats, cvgvstats=None, alpha=1.0, startdim=1):\n",
    "        \"\"\"\n",
    "        perform postfilter based on GV statistics into data\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        data : array, shape(`T`, `dim`)\n",
    "            array of data sequence\n",
    "        gvstats : array, shape(`2`, `dim`)\n",
    "            array of mean and variance for target GV\n",
    "        cvgvstats : array, shape(`2`, `dim`)\n",
    "            array of mean and variance for converted GV\n",
    "        alpha : float\n",
    "            morphing coefficient between GV transformed data and data.\n",
    "            alpha * gvpf(data) + (1 - alpha) * data\n",
    "        startdim : int\n",
    "            start dimension to perform GV postfilter\n",
    "        \n",
    "        Returns\n",
    "        ----------\n",
    "        filtered_data : array, shape(`T`, `data`)\n",
    "            array of GV postfiltered data sequnece\n",
    "        \"\"\"\n",
    "        \n",
    "        # get length and dimension\n",
    "        T, dim = data.shape\n",
    "        assert gvstats is not None\n",
    "        assert dim == gvstats.shape[1]\n",
    "        \n",
    "        # calculate statics of input data\n",
    "        datamean = np.mean(data, axis=0)\n",
    "        \n",
    "        if cvgvstats is None:\n",
    "            # use variance of the given data\n",
    "            datavar = np.var(data, axis=0)\n",
    "        else:\n",
    "            # use variance of trained gv stats\n",
    "            datavar = cvgvstats[0]\n",
    "        \n",
    "        # perform GV postfilter\n",
    "        filterd = np.sqrt(gvstats[0, startdim:] / datavar[startdim:]) * (data[:, startdim:] - datamean[startdim:]) + datamean[startdim:]\n",
    "        \n",
    "        filterd_data = np.c_[data[:, :startdim], filterd]\n",
    "        \n",
    "        return alpha * filterd_data + (1 - alpha) * data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 0. config path\n",
    "__versions = \"pre-stored-ej\"\n",
    "__same_path = \"./utterance/\" + __versions + \"/\"\n",
    "pre_stored_source_list = __same_path + 'pre-source/**/V01/T01/**/*.wav'\n",
    "pre_stored_list = __same_path + \"pre/**/V01/T01/**/*.wav\"\n",
    "output_path = \"./utterance/old-test/orl/ej-f/1/adapt31/\"\n",
    "\n",
    "# 1. estimate features\n",
    "feat = FeatureExtractor()\n",
    "synthesizer = Synthesizer()\n",
    "\n",
    "org_f0list = None\n",
    "org_splist = None\n",
    "org_mceplist = None\n",
    "org_aplist = None\n",
    "org_npowlist = None\n",
    "org_codeaplist = None\n",
    "\n",
    "if os.path.exists(output_path + \"_org_f0.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_org_sp.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_org_ap.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_org_mcep.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_org_npow.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_org_codeap.pickle\"):\n",
    "        \n",
    "    with open(output_path + \"_org_f0.pickle\", 'rb') as f:    \n",
    "        org_f0list = pickle.load(f)\n",
    "    with open(output_path + \"_org_sp.pickle\", 'rb') as f:   \n",
    "        org_splist = pickle.load(f)\n",
    "    with open(output_path + \"_org_ap.pickle\", 'rb') as f:  \n",
    "        org_aplist = pickle.load(f)\n",
    "    with open(output_path + \"_org_mcep.pickle\", 'rb') as f:  \n",
    "        org_mceplist = pickle.load(f)\n",
    "    with open(output_path + \"_org_npow.pickle\", 'rb') as f:  \n",
    "        org_npowlist = pickle.load(f)\n",
    "    with open(output_path + \"_org_codeap.pickle\", 'rb') as f:  \n",
    "        org_codeaplist = pickle.load(f) \n",
    "else:\n",
    "    org_f0list = []\n",
    "    org_splist = []\n",
    "    org_mceplist = []\n",
    "    org_aplist = []\n",
    "    org_npowlist = []\n",
    "    org_codeaplist = []\n",
    "    ite = 0 \n",
    "    for files in sorted(glob.iglob(pre_stored_source_list, recursive=True)):\n",
    "        wavf = files\n",
    "        x, fs = sf.read(wavf)\n",
    "        x = np.array(x, dtype=np.float)\n",
    "        x = low_cut_filter(x, fs, cutoff=70)\n",
    "        assert fs == 16000\n",
    "\n",
    "        print(\"extract acoustic featuers: \" + wavf)\n",
    "\n",
    "        f0, sp, ap = feat.analyze(x)\n",
    "        mcep = feat.mcep()\n",
    "        npow = feat.npow()\n",
    "        codeap = feat.codeap()\n",
    "        wav = synthesizer.synthesis_sp(f0, sp, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"src_ansys_{}_.wav\".format(ite), wav, fs)\n",
    "\n",
    "        org_f0list.append(f0)\n",
    "        org_splist.append(sp)\n",
    "        org_mceplist.append(mcep)\n",
    "        org_aplist.append(ap)\n",
    "        org_npowlist.append(npow)\n",
    "        org_codeaplist.append(codeap)\n",
    "\n",
    "        wav = synthesizer.synthesis(f0, mcep, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"src_mcep_{}_.wav\".format(ite), wav, fs)\n",
    "        ite = ite + 1\n",
    "\n",
    "    with open(output_path + \"_org_f0.pickle\", 'wb') as f:    \n",
    "        pickle.dump(org_f0list, f)\n",
    "    with open(output_path + \"_org_sp.pickle\", 'wb') as f:   \n",
    "        pickle.dump(org_splist, f)\n",
    "    with open(output_path + \"_org_npow.pickle\", 'wb') as f:   \n",
    "        pickle.dump(org_npowlist, f)\n",
    "    with open(output_path + \"_org_ap.pickle\", 'wb') as f:  \n",
    "        pickle.dump(org_aplist, f)\n",
    "    with open(output_path + \"_org_mcep.pickle\", 'wb') as f:  \n",
    "        pickle.dump(org_mceplist, f)\n",
    "    with open(output_path + \"_org_codeap.pickle\", 'wb') as f:  \n",
    "        pickle.dump(org_codeaplist, f) \n",
    "\n",
    "mid_f0list = None\n",
    "mid_mceplist = None\n",
    "mid_aplist = None\n",
    "mid_npowlist = None\n",
    "mid_splist = None\n",
    "mid_codeaplist = None        \n",
    "\n",
    "if os.path.exists(output_path + \"_mid_f0.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_mid_sp_0_.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_mid_ap_0_.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_mid_mcep.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_mid_npow.pickle\") \\\n",
    "    and os.path.exists(output_path + \"_mid_codeap.pickle\"):\n",
    "        \n",
    "    with open(output_path + \"_mid_f0.pickle\", 'rb') as f:    \n",
    "        mid_f0list = pickle.load(f)\n",
    "    for i in range(0, len(org_splist)*21, len(org_splist)):  \n",
    "        with open(output_path + \"_mid_sp_{}_.pickle\".format(i), 'rb') as f:\n",
    "            temp_splist = pickle.load(f)\n",
    "            if mid_splist is None:\n",
    "                mid_splist = temp_splist\n",
    "            else:\n",
    "                mid_splist = mid_splist + temp_splist\n",
    "    for i in range(0, len(org_aplist)*21, len(org_aplist)):  \n",
    "        with open(output_path + \"_mid_ap_{}_.pickle\".format(i), 'rb') as f:\n",
    "            temp_aplist = pickle.load(f)\n",
    "            if mid_aplist is None:\n",
    "                mid_aplist = temp_aplist\n",
    "            else:\n",
    "                mid_aplist = mid_aplist + temp_aplist    \n",
    "    with open(output_path + \"_mid_mcep.pickle\", 'rb') as f:  \n",
    "        mid_mceplist = pickle.load(f)\n",
    "    with open(output_path + \"_mid_npow.pickle\", 'rb') as f:  \n",
    "        mid_npowlist = pickle.load(f)\n",
    "    with open(output_path + \"_mid_codeap.pickle\", 'rb') as f:  \n",
    "        mid_codeaplist = pickle.load(f) \n",
    "else:        \n",
    "    mid_f0list = []\n",
    "    mid_mceplist = []\n",
    "    mid_aplist = []\n",
    "    mid_npowlist = []\n",
    "    mid_splist = []\n",
    "    mid_codeaplist = []\n",
    "    ite = 0\n",
    "    for files in sorted(glob.iglob(pre_stored_list, recursive=True)):\n",
    "        wavf = files\n",
    "        x, fs = sf.read(wavf)\n",
    "        x = np.array(x, dtype=np.float)\n",
    "        x = low_cut_filter(x, fs, cutoff=70)\n",
    "        assert fs == 16000\n",
    "\n",
    "        print(\"extract acoustic featuers: \" + wavf)\n",
    "\n",
    "        f0, sp, ap = feat.analyze(x)\n",
    "        mcep = feat.mcep()\n",
    "        npow = feat.npow()\n",
    "        codeap = feat.codeap()\n",
    "        name, ext = os.path.splitext(wavf)\n",
    "        wav = synthesizer.synthesis_sp(f0, sp, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"mid_ansys_{}_.wav\".format(ite), wav, fs)\n",
    "\n",
    "        org_f0list.append(f0)\n",
    "        org_splist.append(sp)\n",
    "        org_mceplist.append(mcep)\n",
    "        org_aplist.append(ap)\n",
    "        org_npowlist.append(npow)\n",
    "        org_codeaplist.append(codeap)\n",
    "\n",
    "        wav = synthesizer.synthesis(f0, mcep, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"mid_mcep_{}_.wav\".format(ite), wav, fs)\n",
    "        ite = ite + 1\n",
    "        \n",
    "    with open(output_path + \"_mid_f0.pickle\", 'wb') as f:\n",
    "        print(f)\n",
    "        pickle.dump(mid_f0list, f)\n",
    "    with open(output_path + \"_mid_npow.pickle\", 'wb') as f:\n",
    "        print(f)\n",
    "        pickle.dump(mid_npowlist, f)\n",
    "    for i in range(0, len(mid_splist), len(org_splist)):\n",
    "        with open(output_path + \"_mid_sp_{}_.pickle\".format(i), 'wb') as f:   \n",
    "            print(f)\n",
    "            pickle.dump(mid_splist[i:i+len(org_splist)], f)\n",
    "    for i in range(0, len(mid_aplist), len(org_aplist)):\n",
    "        with open(output_path + \"_mid_ap_{}_.pickle\".format(i), 'wb') as f:\n",
    "            print(f)\n",
    "            pickle.dump(mid_aplist[i:i+len(org_aplist)], f)\n",
    "    with open(output_path + \"_mid_mcep.pickle\", 'wb') as f:\n",
    "        print(f)\n",
    "        pickle.dump(mid_mceplist, f)\n",
    "    with open(output_path + \"_mid_codeap.pickle\", 'wb') as f:\n",
    "        print(f)\n",
    "        pickle.dump(mid_codeaplist, f) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GMMTrainer(object):\n",
    "    \"\"\"\n",
    "    this class offers the training of GMM with several types of covariance matrix.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    n_mix : int \n",
    "        the number of mixture components of the GMM\n",
    "    n_iter : int\n",
    "        the number of iteration for EM algorithm\n",
    "    covtype : str\n",
    "        the type of covariance matrix of the GMM\n",
    "        'full': full-covariance matrix\n",
    "    \n",
    "    Attributes\n",
    "    ---------\n",
    "    param : \n",
    "        sklearn-based model parameters of the GMM\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, n_mix=64, n_iter=100, covtype='full', params='wmc'):\n",
    "        self.n_mix = n_mix\n",
    "        self.n_iter = n_iter\n",
    "        self.covtype = covtype\n",
    "        self.params = params\n",
    "        self.param = sklearn.mixture.GMM(n_components=self.n_mix,\n",
    "                                                     covariance_type=self.covtype,\n",
    "                                                     n_iter=self.n_iter, params=self.params)\n",
    "        \n",
    "    def train(self, jnt):\n",
    "        \"\"\"\n",
    "        fit GMM parameter from given joint feature vector\n",
    "        \n",
    "        Parametes\n",
    "        ---------\n",
    "        jnt : array, shape(`T`, `jnt.shape[0]`)\n",
    "            joint feature vector of original and target feature vector consisting of static and delta components\n",
    "        \"\"\"\n",
    "        \n",
    "        if self.covtype == 'full':\n",
    "            self.param.fit(jnt)\n",
    "            \n",
    "        return\n",
    "    \n",
    "class GMMConvertor(object):\n",
    "    \"\"\"\n",
    "    this class offers the several conversion techniques such as Maximum Likelihood Parameter Generation (MLPG)\n",
    "    and Minimum Mean Square Error (MMSE).\n",
    "    \n",
    "    Parametes\n",
    "    ---------\n",
    "    n_mix : int\n",
    "        the number of mixture components of the GMM\n",
    "    covtype : str\n",
    "        the type of covariance matrix of the GMM\n",
    "        'full': full-covariance matrix\n",
    "    gmmmode : str\n",
    "        the type of the GMM for opening\n",
    "        `None` : Normal Joint Density - GMM (JD-GMM)\n",
    "    \n",
    "    Attributes\n",
    "    ---------\n",
    "    param : \n",
    "        sklearn-based model parameters of the GMM\n",
    "    w : shape(`n_mix`)\n",
    "        vector of mixture component weight of the GMM\n",
    "    jmean : shape(`n_mix`, `jnt.shape[0]`)\n",
    "        Array of joint mean vector of the GMM\n",
    "    jcov : shape(`n_mix`, `jnt.shape[0]`, `jnt.shape[0]`)\n",
    "        array of joint covariance matrix of the GMM\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, n_mix=64, covtype='full', gmmmode=None):\n",
    "        self.n_mix = n_mix\n",
    "        self.covtype = covtype\n",
    "        self.gmmmode = gmmmode\n",
    "        \n",
    "    def open_from_param(self, param):\n",
    "        \"\"\"\n",
    "        open GMM from GMMTrainer\n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        param : GMMTrainer\n",
    "            GMMTrainer class\n",
    "        \"\"\"\n",
    "        \n",
    "        self.param = param\n",
    "        self._deploy_parameters()\n",
    "        \n",
    "        return\n",
    "    \n",
    "    def convert(self, data, cvtype='mlpg'):\n",
    "        \"\"\"\n",
    "        convert data based on conditional probability density function\n",
    "        \n",
    "        Parametes\n",
    "        ---------\n",
    "        data : array, shape(`T`, `dim`)\n",
    "            original data will be converted\n",
    "        cvtype : str\n",
    "            type of conversion technique\n",
    "            `mlpg` : maximum likelihood parameter generation\n",
    "            \n",
    "        Returns\n",
    "        ----------\n",
    "        odata : array, shape(`T`, `dim`)\n",
    "            converted data\n",
    "        \"\"\"\n",
    "        \n",
    "        # estimate parameter sequence\n",
    "        cseq, wseq, mseq, covseq = self._gmmmap(data)\n",
    "        \n",
    "        if cvtype == 'mlpg':\n",
    "            odata = self._mlpg(mseq, covseq)\n",
    "        else:\n",
    "            raise ValueError('please choose conversion mode in `mlpg`.')\n",
    "        \n",
    "        return odata\n",
    "    \n",
    "    def _gmmmap(self, sddata):\n",
    "        # paramete for sequencial data\n",
    "        T, sddim = sddata.shape\n",
    "        \n",
    "        # estimate posterior sequence\n",
    "        wseq = self.pX.predict_proba(sddata)\n",
    "        \n",
    "        # estimate mixture sequence\n",
    "        cseq = np.argmax(wseq, axis=1)\n",
    "        \n",
    "        mseq = np.zeros((T, sddim))\n",
    "        covseq = np.zeros((T, sddim, sddim))\n",
    "        for t in range(T):\n",
    "            # read maximum likelihood mixture component in frame t\n",
    "            m = cseq[t]\n",
    "            \n",
    "            # conditional mean vector sequence\n",
    "            mseq[t] = self.meanY[m] + self.A[m] @ (sddata[t] - self.meanX[m])\n",
    "            \n",
    "            # conditional covariance sequence\n",
    "            covseq[t] = self.cond_cov_inv[m]\n",
    "        \n",
    "        return cseq, wseq, mseq, covseq\n",
    "    \n",
    "    def _mlpg(self, mseq, covseq):\n",
    "        # parameter for sequencial data\n",
    "        T, sddim = mseq.shape\n",
    "        \n",
    "        # prepare W\n",
    "        W = construct_static_and_delta_matrix(T, sddim // 2)\n",
    "        \n",
    "        # prepare D\n",
    "        D = get_diagonal_precision_matrix(T, sddim, covseq)\n",
    "        \n",
    "        # calculate W'D\n",
    "        WD = W.T @ D\n",
    "        \n",
    "        # W'DW\n",
    "        WDW = WD @ W\n",
    "        \n",
    "        # W'Dm\n",
    "        WDM = WD @ mseq.flatten()\n",
    "        \n",
    "        # estimate y = (W'DW)^-1 * W'Dm\n",
    "        odata = scipy.sparse.linalg.spsolve(WDW, WDM, use_umfpack=False).reshape(T, sddim // 2)\n",
    "        \n",
    "        return odata\n",
    "    \n",
    "    def _deploy_parameters(self):\n",
    "        # read JD-GMM parameters from self.param\n",
    "        self.W = self.param.weights_\n",
    "        self.jmean = self.param.means_\n",
    "        self.jcov = self.param.covars_\n",
    "        \n",
    "        # devide GMM parameters into source and target parameters\n",
    "        sddim = self.jmean.shape[1] // 2\n",
    "        self.meanX = self.jmean[:, 0:sddim]\n",
    "        self.meanY = self.jmean[:, sddim:]\n",
    "        self.covXX = self.jcov[:, :sddim, :sddim]\n",
    "        self.covXY = self.jcov[:, :sddim, sddim:]\n",
    "        self.covYX = self.jcov[:, sddim:, :sddim]\n",
    "        self.covYY = self.jcov[:, sddim:, sddim:]\n",
    "        \n",
    "        # change model parameter of GMM into that of gmmmode\n",
    "        if self.gmmmode is None:\n",
    "            pass\n",
    "        else:\n",
    "            raise ValueError('please choose GMM mode in [None]')\n",
    "            \n",
    "        # estimate parameters for conversion\n",
    "        self._set_Ab()\n",
    "        self._set_pX()\n",
    "        \n",
    "        return\n",
    "    \n",
    "    def _set_Ab(self):\n",
    "        # calculate A and b from self.jmean, self.jcov\n",
    "        sddim = self.jmean.shape[1] // 2\n",
    "        \n",
    "        # calculate inverse covariance for covariance XX in each mixture\n",
    "        self.covXXinv = np.zeros((self.n_mix, sddim, sddim))\n",
    "        for m in range(self.n_mix):\n",
    "            self.covXXinv[m] = np.linalg.inv(self.covXX[m])\n",
    "            \n",
    "        # calculate A, b, and conditional covariance given X\n",
    "        self.A = np.zeros((self.n_mix, sddim, sddim))\n",
    "        self.b = np.zeros((self.n_mix, sddim))\n",
    "        self.cond_cov_inv = np.zeros((self.n_mix, sddim, sddim))\n",
    "        for m in range(self.n_mix):\n",
    "            # calculate A (A = yxcov_m * xxcov_m^-1)\n",
    "            self.A[m] = self.covYX[m] @ self.covXXinv[m]\n",
    "            \n",
    "            # calculate b (b = mean^Y - A * mean^X)\n",
    "            self.b[m] = self.meanY[m] - self.A[m] @ self.meanX[m]\n",
    "            \n",
    "            # calculate conditional covariance (cov^(Y|X)^-1 = (yycov - A * xycov)^-1)\n",
    "            self.cond_cov_inv[m] = np.linalg.inv(self.covYY[m] - self.A[m] @ self.covXY[m])\n",
    "            \n",
    "        return\n",
    "    \n",
    "    def _set_pX(self):\n",
    "        # probability density function of X \n",
    "        self.pX = sklearn.mixture.GMM(n_components=self.n_mix, covariance_type=self.covtype)\n",
    "        self.pX.weights_ = self.W\n",
    "        self.pX.means_ = self.meanX\n",
    "        self.pX.covars_ = self.covXX\n",
    "        \n",
    "        # following function is required to estimate porsterior\n",
    "        # p(x | \\lambda^(X))\n",
    "        #self.pX.precisions_cholesky_ = _compute_precision_cholesky(self.covXX, self.covtype)\n",
    "        \n",
    "        return\n",
    "    \n",
    "def get_diagonal_precision_matrix(T, D, covseq):\n",
    "    return scipy.sparse.block_diag(covseq, format='csr')        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_alignment(odata, onpow, tdata, tnpow, opow=-20, tpow=-20, sd=0, cvdata=None, given_twf=None, otflag=None, distance='melcd'):\n",
    "    \"\"\"\n",
    "    get alignment between original and target.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    odata : array, shape(`T`, `dim`)\n",
    "        acoustic feature vector of original\n",
    "    onpow : array, shape(`T`)\n",
    "        Normalized power vector of original\n",
    "    tdata : array, shape(`T`, `dim`)\n",
    "        acoustic feature vector of target\n",
    "    tnpow : array, shape(`T`)\n",
    "        Normalized power vector of target\n",
    "    opow : float\n",
    "        power threshold of original\n",
    "    tpow : float\n",
    "        power threshold of target\n",
    "    sd : int\n",
    "        start dimension to be used for alignment\n",
    "    cvdata : array, shape(`T`, `dim`)\n",
    "        converted original data\n",
    "    given_twf : array, shape(`T_new`, `dim * 2`)\n",
    "        Alignment given twf\n",
    "    otflag : str\n",
    "        Alignment into the length of specification\n",
    "        'org' : alignment into original length\n",
    "        'tar' : alignment into target length\n",
    "    distance : str\n",
    "        Distance function to be used\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    jdata : array, shape(`T_new`, `dim * 2`)\n",
    "        joint static and delta feature vector\n",
    "    twf : array, shape(`T_new`, `dim * 2`)\n",
    "        Time warping function\n",
    "    mcd : float\n",
    "        Mel-cepstrum distortion between arrays\n",
    "    \"\"\"\n",
    "    \n",
    "    oexdata = extsddata(odata[:, sd:], onpow, power_threshold=opow)\n",
    "    texdata = extsddata(tdata[:, sd:], tnpow, power_threshold=tpow)\n",
    "    \n",
    "    if cvdata is None:\n",
    "        align_odata = oexdata\n",
    "    else:\n",
    "        cvexdata = extsddata(cvdata, onpow, power_threshold=opow)\n",
    "        align_odata = cvexdata\n",
    "    \n",
    "    if given_twf is None: \n",
    "        twf = estimate_twf(align_odata, texdata, distance=distance, fast=False, otflag=otflag)\n",
    "    else:\n",
    "        twf = given_twf\n",
    "    \n",
    "    jdata = align_data(oexdata, texdata, twf)\n",
    "    mcd = melcd(align_odata[twf[0]], texdata[twf[1]])\n",
    "    \n",
    "    return jdata, twf, mcd\n",
    "\n",
    "def align_feature_vectors(odata, onpows, tdata, tnpows, opow=-100, tpow=-100, itnum=3, sd=0, given_twfs=None, otflag=None):\n",
    "    \"\"\"\n",
    "    get alignment to create joint feature vector\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    odata : list, (`num_files`)\n",
    "        List of original feature vectors\n",
    "    onpow : list, (`num_files`)\n",
    "        List of original npows\n",
    "    tdata : list, (`num_files`)\n",
    "        List of target feature vectors\n",
    "    tnpow : list, (`num_files`)\n",
    "        List of target npows\n",
    "    opow : float\n",
    "        power threshold of original\n",
    "    tpow : float\n",
    "        power threshold of target\n",
    "    itnum : int\n",
    "        the number of iteration\n",
    "    sd : int\n",
    "        start dimension of feature vector to be used for alignment\n",
    "    given_twf : array, shape(`T_new`, `dim * 2`)\n",
    "        use given alignment while 1st iteration\n",
    "    otflag : str\n",
    "        Alignment into the length of specification\n",
    "        'org' : alignment into original length\n",
    "        'tar' : alignment into target length\n",
    "    distance : str\n",
    "        Distance function to be used\n",
    "    \n",
    "    Returns\n",
    "    ----------\n",
    "    jdata : array, shape(`T_new`, `dim * 2`)\n",
    "        joint static and delta feature vector\n",
    "    twf : array, shape(`T_new`, `dim * 2`)\n",
    "        Time warping function\n",
    "    mcd : float\n",
    "        Mel-cepstrum distortion between arrays\n",
    "    \"\"\"\n",
    "    \n",
    "    it = 1\n",
    "    num_files = len(odata)\n",
    "    cvgmm, cvdata = None, None\n",
    "    for it in range(1, itnum+1):\n",
    "        print('{}-th joint feature extraction starts.'.format(it))\n",
    "        \n",
    "        # alignment\n",
    "        twfs, jfvs = [], []\n",
    "        for i in range(num_files):\n",
    "            if it == 1 and given_twfs is not None:\n",
    "                gtwf = given_twfs[i]\n",
    "            else:\n",
    "                gtwf = None\n",
    "            \n",
    "            if it > 1:\n",
    "                cvdata = cvgmm.convert(static_delta(odata[i][:, sd:]))\n",
    "            \n",
    "            jdata, twf, mcd = get_alignment(odata[i], onpows[i], tdata[i], tnpows[i], opow=opow, tpow=tpow,\n",
    "                                            sd=sd, cvdata=cvdata, given_twf=gtwf, otflag=otflag)\n",
    "            twfs.append(twf)\n",
    "            jfvs.append(jdata)\n",
    "            print('distortion [dB] for {}-th file: {}'.format(i+1, mcd))\n",
    "            \n",
    "        jnt_data = transform_jnt(jfvs)\n",
    "        \n",
    "        if it != itnum:\n",
    "            # train GMM, if not final iteration\n",
    "            datagmm = GMMTrainer()\n",
    "            datagmm.train(jnt_data)\n",
    "            cvgmm = GMMConvertor()\n",
    "            cvgmm.open_from_param(datagmm.param)\n",
    "        it += 1\n",
    "    return jfvs, twfs    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2. estimate twf and jnt\n",
    "if os.path.exists(output_path + \"_jnt_mcep_0_.pickle\"):\n",
    "    pass\n",
    "else:\n",
    "    for i in range(0, len(mid_mceplist), len(org_mceplist)):\n",
    "        org_mceps = org_mceplist\n",
    "        org_npows = org_npowlist\n",
    "        mid_mceps = mid_mceplist[i:i+len(org_mceps)]\n",
    "        mid_npows = mid_npowlist[i:i+len(org_npows)]\n",
    "        assert len(org_mceps) == len(mid_mceps)\n",
    "        assert len(org_npows) == len(mid_npows)\n",
    "        assert len(org_mceps) == len(org_npows)\n",
    "\n",
    "        # dtw between original and target 0-th and silence\n",
    "        print(\"## alignment mcep 0-th and silence ##\")\n",
    "        jmceps, twfs = align_feature_vectors(org_mceps, org_npows, mid_mceps, mid_npows, opow=-15, tpow=-15, sd=1)\n",
    "        jnt_mcep = transform_jnt(jmceps)\n",
    "\n",
    "        # save joint feature vectors\n",
    "        with open(output_path + \"_jnt_mcep_{}_.pickle\".format(i), 'wb') as f:   \n",
    "            print(f)\n",
    "            pickle.dump(jnt_mcep, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 3. make EV-GMM\n",
    "initgmm, initgmm_codeap = None, None\n",
    "if os.path.exists(output_path + \"initgmm.pickle\"):\n",
    "    with open(output_path + \"initgmm.pickle\".format(i), 'rb') as f:   \n",
    "        print(f)\n",
    "        initgmm = pickle.load(f)\n",
    "else:\n",
    "    jnt, jnt_codeap = None, []\n",
    "    for i in range(0, len(mid_mceplist), len(org_mceplist)):\n",
    "        with open(output_path + \"_jnt_mcep_{}_.pickle\".format(i), 'rb') as f:  \n",
    "            temp_jnt = pickle.load(f)\n",
    "            if jnt is None:\n",
    "                jnt = temp_jnt\n",
    "            else:\n",
    "                jnt = np.r_[jnt, temp_jnt]\n",
    "\n",
    "    # train initial gmm\n",
    "    initgmm = GMMTrainer()\n",
    "    initgmm.train(jnt)\n",
    "\n",
    "    with open(output_path + \"initgmm.pickle\", 'wb') as f:   \n",
    "        print(f)\n",
    "        pickle.dump(initgmm, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get initial gmm params\n",
    "init_W = initgmm.param.weights_\n",
    "init_jmean = initgmm.param.means_\n",
    "init_jcov = initgmm.param.covars_\n",
    "sddim = init_jmean.shape[1] // 2\n",
    "init_meanX = init_jmean[:, :sddim]\n",
    "init_meanY = init_jmean[:, sddim:]\n",
    "init_covXX = init_jcov[:, :sddim, :sddim]\n",
    "init_covXY = init_jcov[:, :sddim, sddim:]\n",
    "init_covYX = init_jcov[:, sddim:, :sddim]\n",
    "init_covYY = init_jcov[:, sddim:, sddim:]\n",
    "fitted_source = init_meanX\n",
    "fitted_target = init_meanY\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sv = None\n",
    "if os.path.exists(output_path + \"_sv.npy\"):\n",
    "    sv = np.array(sv)\n",
    "    sv = np.load(output_path + '_sv.npy')\n",
    "else:\n",
    "    depengmm, depenjnt = None, None\n",
    "    sv = []\n",
    "    for i in range(0, len(mid_mceplist), len(org_mceplist)):\n",
    "        with open(output_path + \"_jnt_mcep_{}_.pickle\".format(i), 'rb') as f:  \n",
    "            depenjnt = pickle.load(f)\n",
    "            depengmm = GMMTrainer(params='m')\n",
    "            depengmm.param.weights_ = init_W\n",
    "            depengmm.param.means_ = init_jmean\n",
    "            depengmm.param.covars_ = init_jcov\n",
    "            depengmm.train(depenjnt)\n",
    "            sv.append(depengmm.param.means_)\n",
    "    sv = np.array(sv)\n",
    "    np.save(output_path + \"_sv\", sv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_mix = 64\n",
    "S = int(len(mid_mceplist) / len(org_mceplist))\n",
    "assert S == 22\n",
    "\n",
    "source_pca = sklearn.decomposition.PCA()\n",
    "source_pca.fit(sv[:,:,:sddim].reshape((S, n_mix*sddim)))\n",
    "\n",
    "target_pca = sklearn.decomposition.PCA()\n",
    "target_pca.fit(sv[:,:,sddim:].reshape((S, n_mix*sddim)))\n",
    "\n",
    "eigenvectors = source_pca.components_.reshape((n_mix, sddim, S)), target_pca.components_.reshape((n_mix, sddim, S)) \n",
    "biasvectors = source_pca.mean_.reshape((n_mix, sddim)), target_pca.mean_.reshape((n_mix, sddim))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# estimate statistic features\n",
    "for_convert_source = __same_path + 'input/EJM10/V01/T01/TIMIT/000/*.wav'\n",
    "for_convert_target = __same_path + 'adaptation/EJF01/V01/T01/ATR503/A/*.wav'\n",
    "\n",
    "src_f0list = []\n",
    "src_splist = []\n",
    "src_mceplist = []\n",
    "src_aplist = []\n",
    "src_npowlist = []\n",
    "src_codeaplist = []\n",
    "if os.path.exists(__same_path + 'input/EJM10/V01/T01/TIMIT/000/A11.wav'):\n",
    "    ite = 0\n",
    "    for files in sorted(glob.iglob(for_convert_source, recursive=True)):\n",
    "        wavf = files\n",
    "        x, fs = sf.read(wavf)\n",
    "        x = np.array(x, dtype=np.float)\n",
    "        x = low_cut_filter(x, fs, cutoff=70)\n",
    "        assert fs == 16000\n",
    "\n",
    "        print(\"extract acoustic featuers: \" + wavf)\n",
    "\n",
    "        f0, sp, ap = feat.analyze(x)\n",
    "        mcep = feat.mcep()\n",
    "        npow = feat.npow()\n",
    "        codeap = feat.codeap()\n",
    "        wav = synthesizer.synthesis_sp(f0, sp, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"input_ansys_{}_.wav\".format(ite), wav, fs)\n",
    "\n",
    "        org_f0list.append(f0)\n",
    "        org_splist.append(sp)\n",
    "        org_mceplist.append(mcep)\n",
    "        org_aplist.append(ap)\n",
    "        org_npowlist.append(npow)\n",
    "        org_codeaplist.append(codeap)\n",
    "\n",
    "        wav = synthesizer.synthesis(f0, mcep, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"input_mcep_{}_.wav\".format(ite), wav, fs)\n",
    "        ite = ite + 1\n",
    "\n",
    "else:\n",
    "    raise ValueError(\"No such files.\")\n",
    "tar_f0list = []\n",
    "tar_mceplist = []\n",
    "tar_aplist = []\n",
    "tar_npowlist = []\n",
    "tar_splist = []\n",
    "tar_codeaplist = []\n",
    "\n",
    "if os.path.exists(__same_path + 'adaptation/EJF01/V01/T01/ATR503/A/A01.wav'):\n",
    "    ite = 0\n",
    "    for files in sorted(glob.iglob(for_convert_target, recursive=True)):\n",
    "        wavf = files\n",
    "        x, fs = sf.read(wavf)\n",
    "        x = np.array(x, dtype=np.float)\n",
    "        x = low_cut_filter(x, fs, cutoff=70)\n",
    "        assert fs == 16000\n",
    "\n",
    "        print(\"extract acoustic featuers: \" + wavf)\n",
    "\n",
    "        f0, sp, ap = feat.analyze(x)\n",
    "        mcep = feat.mcep()\n",
    "        npow = feat.npow()\n",
    "        codeap = feat.codeap()\n",
    "        name, ext = os.path.splitext(wavf)\n",
    "        wav = synthesizer.synthesis_sp(f0, sp, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"target_ansys_{}_.wav\".format(ite), wav, fs)\n",
    "\n",
    "        org_f0list.append(f0)\n",
    "        org_splist.append(sp)\n",
    "        org_mceplist.append(mcep)\n",
    "        org_aplist.append(ap)\n",
    "        org_npowlist.append(npow)\n",
    "        org_codeaplist.append(codeap)\n",
    "\n",
    "        wav = synthesizer.synthesis(f0, mcep, ap)\n",
    "        wav = np.clip(wav, -32768, 32767)\n",
    "        sf.write(output_path + \"target_mcep_{}_.wav\".format(ite), wav, fs)\n",
    "        ite = ite + 1\n",
    "else:\n",
    "    raise ValueError(\"No such files.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f0statis = F0statistics()\n",
    "tarf0stats = f0statis.estimate(tar_f0list)\n",
    "srcf0stats = f0statis.estimate(org_f0list)\n",
    "\n",
    "gv = GV()\n",
    "srcgvstats = gv.estimate(org_mceplist)\n",
    "targvstats = gv.estimate(tar_mceplist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 5. fitting target\n",
    "epoch = 100\n",
    "\n",
    "fitgmm = sklearn.mixture.GMM(n_components=n_mix,\n",
    "                             covariance_type='full',\n",
    "                             n_iter=100)\n",
    "fitgmm.weights_ = init_W\n",
    "fitgmm.means_ = init_meanY\n",
    "fitgmm.covars_ = init_covYY\n",
    "\n",
    "for i in range(len(tar_mceplist)):\n",
    "    print(\"adapt: \", i+1, \"/\", len(tar_mceplist))\n",
    "    target = tar_mceplist[i]\n",
    "    target_pow = target[:, 0]\n",
    "    target = target[:, 1:]\n",
    "    for x in range(epoch):\n",
    "        print(\"epoch = \", x)\n",
    "        predict = fitgmm.predict_proba(np.atleast_2d(static_delta(target)))\n",
    "        \n",
    "        y = np.sum([predict[:, k:k+1] * (static_delta(target) - biasvectors[1][k]) for k in range(n_mix)], axis=1)\n",
    "        gamma = np.sum(predict, axis=0)\n",
    "        \n",
    "        left = np.sum([gamma[k] * np.dot(eigenvectors[1][k].T,\n",
    "                                         np.linalg.solve(fitgmm.covars_, eigenvectors[1])[k])\n",
    "                                                        for k in range(n_mix)], axis=0)\n",
    "        right = np.sum([np.dot(eigenvectors[1][k].T,\n",
    "                               np.linalg.solve(fitgmm.covars_, y)[k])\n",
    "                        for k in range(n_mix)], axis=0)\n",
    "        weight = np.linalg.solve(left, right)\n",
    "        \n",
    "        fitted_target = np.dot(eigenvectors[1], weight) + biasvectors[1]\n",
    "        fitgmm.means_ = fitted_target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mcepconvert(source, weights, jmean, meanX, covarXX, covarXY, covarYX, covarYY,\n",
    "                fitted_source, fitted_target):\n",
    "    \n",
    "    M = 64\n",
    "    # set pX\n",
    "    px = sklearn.mixture.GMM(n_components=M, covariance_type='full', n_iter=100)\n",
    "    px.weights_ = weights\n",
    "    px.means_ = meanX\n",
    "    px.covars_ = covarXX\n",
    "    \n",
    "    # set Ab\n",
    "    sddim = jmean.shape[1] // 2\n",
    "    covXXinv = np.zeros((M, sddim, sddim))\n",
    "    for m in range(M):\n",
    "        covXXinv[m] = np.linalg.inv(covarXX[m])\n",
    "    A = np.zeros((M, sddim, sddim))\n",
    "    b = np.zeros((M, sddim))\n",
    "    cond_cov_inv = np.zeros((M, sddim, sddim))\n",
    "    for m in range(M):\n",
    "        A[m] = covarYX[m] @ covXXinv[m]\n",
    "        b[m] = fitted_target[m] - A[m] @ meanX[m]\n",
    "        cond_cov_inv[m] = np.linalg.inv(covarYY[m] - A[m] @ covarXY[m])\n",
    "        \n",
    "    # _gmmmap\n",
    "    T, sddim = source.shape\n",
    "    wseq = px.predict_proba(source)\n",
    "    cseq = np.argmax(wseq, axis=1)\n",
    "    mseq = np.zeros((T, sddim))\n",
    "    covseq = np.zeros((T, sddim, sddim))\n",
    "    for t in range(T):\n",
    "        m = cseq[t]\n",
    "        mseq[t] = fitted_target[m] + A[m] @ (source[t] - meanX[m])\n",
    "        covseq[t] = cond_cov_inv[m]\n",
    "        \n",
    "    # _mlpg\n",
    "    T, sddim = mseq.shape\n",
    "    W = construct_static_and_delta_matrix(T, sddim // 2)\n",
    "    D = get_diagonal_precision_matrix(T, sddim, covseq)\n",
    "    WD = W.T @ D\n",
    "    WDW = WD @ W\n",
    "    WDM = WD @ mseq.flatten()\n",
    "    \n",
    "    output = scipy.sparse.linalg.spsolve(WDW, WDM, use_umfpack=False).reshape(T, sddim // 2)\n",
    "    return output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# learn cvgvstats \n",
    "cv_mceps = []\n",
    "for i in range(len(src_mceplist)):\n",
    "    temp_mcep = src_mceplist[i]\n",
    "    temp_mcep_0th = temp_mcep[:, 0]\n",
    "    temp_mcep = temp_mcep[:, 1:]\n",
    "    sta_mcep = static_delta(temp_mcep)\n",
    "    cvmcep_wopow = np.array(mcepconvert(sta_mcep, init_W, init_jmean, init_meanX,\n",
    "                                        init_covXX, init_covXY, init_covYX, init_covYY,\n",
    "                                        fitted_source, fitted_target))\n",
    "    cvmcep = np.c_[temp_mcep_0th, cvmcep_wopow]\n",
    "    cv_mceps.append(cvmcep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvgvstats = gv.estimate(cv_mceps)\n",
    "\n",
    "for i in range(len(src_mceplist)):\n",
    "    cvmcep_wGV = gv.postfilter(cv_mceps[i], targvstats, cvgvstats=cvgvstats)\n",
    "    cvf0 = f0statis.convert(src_f0list[i], srcf0stats, tarf0stats)\n",
    "    wav = synthesizer.synthesis(cvf0, cvmcep_wGV, src_aplist[i], rmcep=src_mceplist[i])\n",
    "    wav = np.clip(wav, -32768, 32767)\n",
    "    sf.write(output_path + \"cv_{}_.wav\".format(i), wav, 16000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(src_mceplist)):\n",
    "    wav = synthesizer.synthesis(src_f0list[i], src_mceplist[i], src_aplist[i])\n",
    "    wav = np.clip(wav, -32768, 32767)\n",
    "    sf.write(output_path + \"mcep_{}_.wav\".format(i), wav, 16000)\n",
    "    wav = synthesizer.synthesis_sp(src_f0list[i], src_splist[i], src_aplist[i])\n",
    "    wav = np.clip(wav, -32768, 32767)\n",
    "    sf.write(output_path + \"ansys_{}_.wav\".format(i), wav, 16000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvf0 = f0statis.convert(src_f0list[0], srcf0stats, tarf0stats)\n",
    "plt.plot(cvf0)\n",
    "plt.plot(src_f0list[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cvmcep_wGV = gv.postfilter(cv_mceps[0], srcgvstats, cvgvstats=cvgvstats)\n",
    "cvf0 = f0statis.convert(src_f0list[0], srcf0stats, tarf0stats)\n",
    "wav = synthesizer.synthesis(cvf0, cvmcep_wGV, src_aplist[0], rmcep=src_mceplist[0])\n",
    "wav = np.clip(wav, -32768, 32767)\n",
    "sf.write(output_path + \"te.wav\", wav, 16000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}