Measure Pitch, HNR, Jitter, Shimmer, and Formants.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Measure Pitch, HNR, Jitter, Shimmer, Formants, and Estimate VTL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import the external modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import parselmouth \n",
    "import statistics\n",
    "\n",
    "\n",
    "from parselmouth.praat import call\n",
    "from scipy.stats.mstats import zscore\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This function measures duration, pitch, HNR, jitter, and shimmer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is the function to measure source acoustics using default male parameters.\n",
    "\n",
    "def measurePitch(voiceID, f0min, f0max, unit):\n",
    "    sound = parselmouth.Sound(voiceID) # read the sound\n",
    "    duration = call(sound, \"Get total duration\") # duration\n",
    "    pitch = call(sound, \"To Pitch\", 0.0, f0min, f0max) #create a praat pitch object\n",
    "    meanF0 = call(pitch, \"Get mean\", 0, 0, unit) # get mean pitch\n",
    "    stdevF0 = call(pitch, \"Get standard deviation\", 0 ,0, unit) # get standard deviation\n",
    "    harmonicity = call(sound, \"To Harmonicity (cc)\", 0.01, f0min, 0.1, 1.0)\n",
    "    hnr = call(harmonicity, \"Get mean\", 0, 0)\n",
    "    pointProcess = call(sound, \"To PointProcess (periodic, cc)\", f0min, f0max)\n",
    "    localJitter = call(pointProcess, \"Get jitter (local)\", 0, 0, 0.0001, 0.02, 1.3)\n",
    "    localabsoluteJitter = call(pointProcess, \"Get jitter (local, absolute)\", 0, 0, 0.0001, 0.02, 1.3)\n",
    "    rapJitter = call(pointProcess, \"Get jitter (rap)\", 0, 0, 0.0001, 0.02, 1.3)\n",
    "    ppq5Jitter = call(pointProcess, \"Get jitter (ppq5)\", 0, 0, 0.0001, 0.02, 1.3)\n",
    "    ddpJitter = call(pointProcess, \"Get jitter (ddp)\", 0, 0, 0.0001, 0.02, 1.3)\n",
    "    localShimmer =  call([sound, pointProcess], \"Get shimmer (local)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    localdbShimmer = call([sound, pointProcess], \"Get shimmer (local_dB)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    apq3Shimmer = call([sound, pointProcess], \"Get shimmer (apq3)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    aqpq5Shimmer = call([sound, pointProcess], \"Get shimmer (apq5)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    apq11Shimmer =  call([sound, pointProcess], \"Get shimmer (apq11)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    ddaShimmer = call([sound, pointProcess], \"Get shimmer (dda)\", 0, 0, 0.0001, 0.02, 1.3, 1.6)\n",
    "    \n",
    "    return duration, meanF0, stdevF0, hnr, localJitter, localabsoluteJitter, rapJitter, ppq5Jitter, ddpJitter, localShimmer, localdbShimmer, apq3Shimmer, aqpq5Shimmer, apq11Shimmer, ddaShimmer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This function measures formants at each glottal pulse\n",
    "\n",
    "Puts, D. A., Apicella, C. L., & Cárdenas, R. A. (2012). Masculine voices signal men's threat potential in forager and industrial societies. Proceedings of the Royal Society of London B: Biological Sciences, 279(1728), 601-609.\n",
    "\n",
    "Adapted from: DOI 10.17605/OSF.IO/K2BHS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function measures formants using Formant Position formula\n",
    "def measureFormants(sound, wave_file, f0min,f0max):\n",
    "    sound = parselmouth.Sound(sound) # read the sound\n",
    "    pitch = call(sound, \"To Pitch (cc)\", 0, f0min, 15, 'no', 0.03, 0.45, 0.01, 0.35, 0.14, f0max)\n",
    "    pointProcess = call(sound, \"To PointProcess (periodic, cc)\", f0min, f0max)\n",
    "    \n",
    "    formants = call(sound, \"To Formant (burg)\", 0.0025, 5, 5000, 0.025, 50)\n",
    "    numPoints = call(pointProcess, \"Get number of points\")\n",
    "\n",
    "    f1_list = []\n",
    "    f2_list = []\n",
    "    f3_list = []\n",
    "    f4_list = []\n",
    "    \n",
    "    # Measure formants only at glottal pulses\n",
    "    for point in range(0, numPoints):\n",
    "        point += 1\n",
    "        t = call(pointProcess, \"Get time from index\", point)\n",
    "        f1 = call(formants, \"Get value at time\", 1, t, 'Hertz', 'Linear')\n",
    "        f2 = call(formants, \"Get value at time\", 2, t, 'Hertz', 'Linear')\n",
    "        f3 = call(formants, \"Get value at time\", 3, t, 'Hertz', 'Linear')\n",
    "        f4 = call(formants, \"Get value at time\", 4, t, 'Hertz', 'Linear')\n",
    "        f1_list.append(f1)\n",
    "        f2_list.append(f2)\n",
    "        f3_list.append(f3)\n",
    "        f4_list.append(f4)\n",
    "    \n",
    "    f1_list = [f1 for f1 in f1_list if str(f1) != 'nan']\n",
    "    f2_list = [f2 for f2 in f2_list if str(f2) != 'nan']\n",
    "    f3_list = [f3 for f3 in f3_list if str(f3) != 'nan']\n",
    "    f4_list = [f4 for f4 in f4_list if str(f4) != 'nan']\n",
    "    \n",
    "    # calculate mean formants across pulses\n",
    "    f1_mean = statistics.mean(f1_list)\n",
    "    f2_mean = statistics.mean(f2_list)\n",
    "    f3_mean = statistics.mean(f3_list)\n",
    "    f4_mean = statistics.mean(f4_list)\n",
    "    \n",
    "    # calculate median formants across pulses, this is what is used in all subsequent calcualtions\n",
    "    # you can use mean if you want, just edit the code in the boxes below to replace median with mean\n",
    "    f1_median = statistics.median(f1_list)\n",
    "    f2_median = statistics.median(f2_list)\n",
    "    f3_median = statistics.median(f3_list)\n",
    "    f4_median = statistics.median(f4_list)\n",
    "    \n",
    "    return f1_mean, f2_mean, f3_mean, f4_mean, f1_median, f2_median, f3_median, f4_median\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This function runs a 2-factor Principle Components Analysis (PCA) on Jitter and Shimmer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "def runPCA(df):\n",
    "    # z-score the Jitter and Shimmer measurements\n",
    "    measures = ['localJitter', 'localabsoluteJitter', 'rapJitter', 'ppq5Jitter', 'ddpJitter',\n",
    "                'localShimmer', 'localdbShimmer', 'apq3Shimmer', 'apq5Shimmer', 'apq11Shimmer', 'ddaShimmer']\n",
    "    x = df.loc[:, measures].values\n",
    "    x = StandardScaler().fit_transform(x)\n",
    "    # PCA\n",
    "    pca = PCA(n_components=2)\n",
    "    principalComponents = pca.fit_transform(x)\n",
    "    principalDf = pd.DataFrame(data = principalComponents, columns = ['JitterPCA', 'ShimmerPCA'])\n",
    "    principalDf\n",
    "    return principalDf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This block of code runs the above functions on all of the '.wav' files in the /audio folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create lists to put the results\n",
    "file_list = []\n",
    "duration_list = []\n",
    "mean_F0_list = []\n",
    "sd_F0_list = []\n",
    "hnr_list = []\n",
    "localJitter_list = []\n",
    "localabsoluteJitter_list = []\n",
    "rapJitter_list = []\n",
    "ppq5Jitter_list = []\n",
    "ddpJitter_list = []\n",
    "localShimmer_list = []\n",
    "localdbShimmer_list = []\n",
    "apq3Shimmer_list = []\n",
    "aqpq5Shimmer_list = []\n",
    "apq11Shimmer_list = []\n",
    "ddaShimmer_list = []\n",
    "f1_mean_list = []\n",
    "f2_mean_list = []\n",
    "f3_mean_list = []\n",
    "f4_mean_list = []\n",
    "f1_median_list = []\n",
    "f2_median_list = []\n",
    "f3_median_list = []\n",
    "f4_median_list = []\n",
    "\n",
    "# Go through all the wave files in the folder and measure all the acoustics\n",
    "for wave_file in glob.glob(\"audio/*.wav\"):\n",
    "    sound = parselmouth.Sound(wave_file)\n",
    "    (duration, meanF0, stdevF0, hnr, localJitter, localabsoluteJitter, rapJitter, ppq5Jitter, ddpJitter, \n",
    "     localShimmer, localdbShimmer, apq3Shimmer, aqpq5Shimmer, apq11Shimmer, ddaShimmer) = measurePitch(\n",
    "        sound, 75, 300, \"Hertz\")\n",
    "    (f1_mean, f2_mean, f3_mean, f4_mean, f1_median, f2_median, f3_median, f4_median) = measureFormants(\n",
    "        sound, wave_file, 75, 300)\n",
    "    file_list.append(wave_file) # make an ID list\n",
    "    duration_list.append(duration) # make duration list\n",
    "    mean_F0_list.append(meanF0) # make a mean F0 list\n",
    "    sd_F0_list.append(stdevF0) # make a sd F0 list\n",
    "    hnr_list.append(hnr) #add HNR data\n",
    "    \n",
    "    # add raw jitter and shimmer measures\n",
    "    localJitter_list.append(localJitter)\n",
    "    localabsoluteJitter_list.append(localabsoluteJitter)\n",
    "    rapJitter_list.append(rapJitter)\n",
    "    ppq5Jitter_list.append(ppq5Jitter)\n",
    "    ddpJitter_list.append(ddpJitter)\n",
    "    localShimmer_list.append(localShimmer)\n",
    "    localdbShimmer_list.append(localdbShimmer)\n",
    "    apq3Shimmer_list.append(apq3Shimmer)\n",
    "    aqpq5Shimmer_list.append(aqpq5Shimmer)\n",
    "    apq11Shimmer_list.append(apq11Shimmer)\n",
    "    ddaShimmer_list.append(ddaShimmer)\n",
    "    \n",
    "    # add the formant data\n",
    "    f1_mean_list.append(f1_mean)\n",
    "    f2_mean_list.append(f2_mean)\n",
    "    f3_mean_list.append(f3_mean)\n",
    "    f4_mean_list.append(f4_mean)\n",
    "    f1_median_list.append(f1_median)\n",
    "    f2_median_list.append(f2_median)\n",
    "    f3_median_list.append(f3_median)\n",
    "    f4_median_list.append(f4_median)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This block of code adds all of that data we just generated to a Pandas data frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the data to Pandas\n",
    "df = pd.DataFrame(np.column_stack([file_list, duration_list, mean_F0_list, sd_F0_list, hnr_list, \n",
    "                                   localJitter_list, localabsoluteJitter_list, rapJitter_list, \n",
    "                                   ppq5Jitter_list, ddpJitter_list, localShimmer_list, \n",
    "                                   localdbShimmer_list, apq3Shimmer_list, aqpq5Shimmer_list, \n",
    "                                   apq11Shimmer_list, ddaShimmer_list, f1_mean_list, \n",
    "                                   f2_mean_list, f3_mean_list, f4_mean_list, \n",
    "                                   f1_median_list, f2_median_list, f3_median_list, \n",
    "                                   f4_median_list]),\n",
    "                                   columns=['voiceID', 'duration', 'meanF0Hz', 'stdevF0Hz', 'HNR', \n",
    "                                            'localJitter', 'localabsoluteJitter', 'rapJitter', \n",
    "                                            'ppq5Jitter', 'ddpJitter', 'localShimmer', \n",
    "                                            'localdbShimmer', 'apq3Shimmer', 'apq5Shimmer', \n",
    "                                            'apq11Shimmer', 'ddaShimmer', 'f1_mean', 'f2_mean', \n",
    "                                            'f3_mean', 'f4_mean', 'f1_median', \n",
    "                                            'f2_median', 'f3_median', 'f4_median'])\n",
    "\n",
    "pcaData = runPCA(df) # Run jitter and shimmer PCA\n",
    "df = pd.concat([df, pcaData], axis=1) # Add PCA data\n",
    "# reload the data so it's all numbers\n",
    "df.to_csv(\"processed_results.csv\", index=False)\n",
    "df = pd.read_csv('processed_results.csv', header=0)\n",
    "df.sort_values('voiceID').head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Next we calculate the vocal-tract length estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Formant position\n",
    " Puts, D. A., Apicella, C. L., & Cárdenas, R. A. (2012). Masculine voices signal men's threat potential in forager and industrial societies. Proceedings of the Royal Society of London B: Biological Sciences, 279(1728), 601-609."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/david/anaconda3/lib/python3.6/site-packages/scipy/stats/stats.py:2242: RuntimeWarning: Mean of empty slice.\n",
      "  mns = a.mean(axis=axis)\n",
      "/home/david/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n",
      "/home/david/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:135: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  keepdims=keepdims)\n",
      "/home/david/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:105: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)\n",
      "/home/david/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py:127: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  ret = ret.dtype.type(ret / rcount)\n"
     ]
    }
   ],
   "source": [
    "df['pF'] = (zscore(df.f1_median) + zscore(df.f2_median) + zscore(df.f3_median) + zscore(df.f4_median)) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Formant Dispersion\n",
    "Fitch, W. T. (1997). Vocal tract length and formant frequency dispersion correlate with body size in rhesus macaques. The Journal of the Acoustical Society of America, 102(2), 1213-1222."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['fdisp'] = (df['f4_median'] - df['f1_median']) / 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fn (Average Formant)\n",
    "Pisanski, K., & Rendall, D. (2011). The prioritization of voice fundamental frequency or formants in listeners’ assessments of speaker size, masculinity, and attractiveness. The Journal of the Acoustical Society of America, 129(4), 2201-2212."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['avgFormant'] = (df['f1_median'] + df['f2_median'] + df['f3_median'] + df['f4_median']) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MFF \n",
    "Smith, D. R., & Patterson, R. D. (2005). The interaction of glottal-pulse rate and vocal-tract length in judgements of speaker size, sex, and age. The Journal of the Acoustical Society of America, 118(5), 3177-3186."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['mff'] = (df['f1_median'] * df['f2_median'] * df['f3_median'] * df['f4_median']) ** 0.25"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitch VTL\n",
    "Fitch, W. T. (1997). Vocal tract length and formant frequency dispersion correlate with body size in rhesus macaques. The Journal of the Acoustical Society of America, 102(2), 1213-1222."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reload the data again\n",
    "df.to_csv(\"processed_results.csv\", index=False)\n",
    "df = pd.read_csv('processed_results.csv', header=0)\n",
    "\n",
    "df['fitch_vtl'] = ((1 * (35000 / (4 * df['f1_median']))) +\n",
    "                   (3 * (35000 / (4 * df['f2_median']))) + \n",
    "                   (5 * (35000 / (4 * df['f3_median']))) + \n",
    "                   (7 * (35000 / (4 * df['f4_median'])))) / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $\\Delta$F \n",
    "Reby,D.,& McComb,K.(2003). Anatomical constraints generate honesty: acoustic cues to age and weight in the roars of red deer stags. Animal Behaviour, 65, 519e-530."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "xysum = (0.5 * df['f1_median']) + (1.5 * df['f2_median']) + (2.5 * df['f3_median']) + (3.5 * df['f4_median'])\n",
    "xsquaredsum = (0.5 ** 2) + (1.5 ** 2) + (2.5 ** 2) + (3.5 ** 2)\n",
    "df['delta_f'] = xysum / xsquaredsum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VTL($\\Delta$F)\n",
    "Reby,D.,&McComb,K.(2003).Anatomical constraints generate honesty: acoustic cues to age and weight in the roars of red deer stags. Animal Behaviour, 65, 519e-530."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['vtl_delta_f'] = 35000 / (2 * df['delta_f'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the final data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write out the final dataframe\n",
    "df.to_csv(\"processed_results2.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run this to tell you when it's done"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "finished\n"
     ]
    }
   ],
   "source": [
    "print(\"finished\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}