<jupyter_start><jupyter_text>Plotting diphthongsFirst, we import some standard numeric and plotting libraries, and set some basic defaults to apply to all following plots.<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

%matplotlib inline
plt.rc('figure', figsize=(16, 9))  # make plots wider
plt.rc('font', size=16)            # make text bigger<jupyter_output><empty_output><jupyter_text>Since vowel plots are pretty specialized, let's write a reusable function that will put our axes on the top and right sides, and invert the axis directions.  We'll do this with a function called `garnish` that can be applied to any plot.<jupyter_code>def garnish(ax):
    ax.xaxis.tick_top()
    ax.yaxis.tick_right()
    ax.set_xlabel('F2')
    ax.set_ylabel('F1')
    ax.xaxis.set_label_position('top')
    ax.yaxis.set_label_position('right')
    ax.set_xlim(2800, 600)
    ax.set_ylim(1100, 100)
    return ax<jupyter_output><empty_output><jupyter_text>Lit reviewHere's what some studies in the literature have done:<jupyter_code>pd.read_csv('lit-review.csv')<jupyter_output><empty_output><jupyter_text>Load the data<jupyter_code>formant_data = pd.read_csv('diphthongs-bVd.csv')
formant_data.head()<jupyter_output><empty_output><jupyter_text>Data cleaningAfter loading the data, we notice that the `pct` column is strings. It would be easier if they were numbers (for sorting purposes), so let's convert them. Let's also reshape the data from "long format" (1 formant measurement per line) into "wide format" (1 record per line).  The way we do this will have the side effect of dropping the `time` values (which we don't need anyway).  The resulting table will have "heirarchical indexes" on both the columns and rows, which will make it easy to select whole spans of rows or columns at a time.<jupyter_code>formant_data['pct'] = formant_data['pct'].map(lambda x: int(x.strip('%')))
wide_data = pd.pivot_table(formant_data, index=['ipa', 'filename'],
                           columns=['formant', 'pct'], values='formant_value')
wide_data.head(10)<jupyter_output><empty_output><jupyter_text>it will be useful for later to have the IPA labels in both the table index *and* in their own column, so let's do that now. While we're at it, let's remove the heirarchical index names (`formant`, `pct`, etc), so our software doesn't get confused between the `ipa` column and the `ipa` sequence of row labels. Finally, for obscure reasons, indexing works best if the top level of heirarchical column names is in alphabetic order, so we'll do that now too.<jupyter_code>idx = wide_data.index
wide_data.reset_index(0, inplace=True, col_level=0)
wide_data.index = idx
wide_data.columns.rename(['', ''], inplace=True)
wide_data.index.rename(['', ''], inplace=True)
wide_data = wide_data[['f1', 'f2', 'ipa']]  # sort top-level
wide_data.head()<jupyter_output><empty_output><jupyter_text>While we're munging data, let's add a column for color, and make the color the same for each instance of a given diphthong. Then we'll show an example of heirarchical indexing at work:<jupyter_code>colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
color_mapping = dict(zip(wide_data['ipa'].unique(), colors))
wide_data['color'] = wide_data['ipa'].replace(color_mapping)

# example of heirarchical indexing:
wide_data.loc[['aʊ', 'oʊ'], [('ipa', ''), ('color', ''), ('f1', 50), ('f2', 50)]]<jupyter_output><empty_output><jupyter_text>PlottingLet's write a reusable function that will plot the vowels' IPA symbols instead of points.<jupyter_code>def plot_vowel(f1, f2, vowel, ax=None, **kwargs):
    if ax is None:
        fig, ax = plt.subplots()
    ax.text(f2, f1, vowel, ha='center', va='center', **kwargs)
    ax = garnish(ax)
    return ax<jupyter_output><empty_output><jupyter_text>For starters, let's plot just the midpoint of each vowel, like we would do with monophthongs. We could do it all at once with a command like this:```pythonplot_vowels(wide_data[('f1', 50)], wide_data[('f2', 50)], wide_data['ipa'])```but it would be nicer if the color was the same for each token of the same diphthong.  So we'll use a for-loop to do that:<jupyter_code>fig, ax = plt.subplots()
for row in wide_data[[('f1', 50), ('f2', 50),
                      ('ipa', ''), ('color', '')]].itertuples(index=False):
    plot_vowel(row[0], row[1], row[2], color=row[3], ax=ax);<jupyter_output><empty_output><jupyter_text>Plotting 2 points (but which two?)Now let's pick 2 timepoints in each diphthong. We'll want to connect the start- and end-points with a line, so we'll need a new plotting function.<jupyter_code>def plot_diphthong(f1, f2, vowel, ax=None, **kwargs):
    '''f1 and f2 should be matrices, where each column is one diphthong,
       with later timepoints at higher row numbers'''
    f1 = np.array(f1)
    f2 = np.array(f2)
    if ax is None:
        fig, ax = plt.subplots()
    ax.plot(f2, f1, **kwargs, linestyle='-', marker='o', markersize=3)
    ax = plot_vowel(f1[0], f2[0], vowel, ax=ax)
    ax = garnish(ax)
    return ax<jupyter_output><empty_output><jupyter_text>We could pick the `20%` and `80%` points, or `30%`-`70%`, or `10%`-`90%`... which is best?<jupyter_code>pairs = [[5, 95], [10, 90], [15, 85], [20, 80], [25, 75], [30, 70]]

fig, axs = plt.subplots(2, 3, figsize=(16, 15))
for pair, ax in zip(pairs, axs.ravel()):
    f1 = wide_data[[('f1', pair[0]), ('f1', pair[1])]]
    f2 = wide_data[[('f2', pair[0]), ('f2', pair[1])]]
    for this_f1, this_f2, ipa, color in zip(f1.values, f2.values,
                                            wide_data['ipa'],
                                            wide_data['color']):
        ax = plot_diphthong(this_f1, this_f2, ipa, color=color, ax=ax)
    ax.set_title('{} to {}'.format(*pair), y=1.1);<jupyter_output><empty_output><jupyter_text>Looking at the red traces for oʊ, notice that the `30-70%` lines go up and to the right, whereas the `5-95%` lines go up and to the left.  Let's plot several spans of just that vowel to get a closer look:<jupyter_code>ou = wide_data.loc['oʊ', ['f1', 'f2', 'ipa', 'color']]

fig, ax = plt.subplots()
for pair in pairs:
    for f1, f2, color in zip(ou[[('f1', pair[0]), ('f1', pair[1])]].values,
                             ou[[('f2', pair[0]), ('f2', pair[1])]].values,
                             colors):
        ax = plot_diphthong(f1, f2, pair[0], ax=ax, color=color, linewidth=0.7)
        ax = plot_vowel(f1[-1], f2[-1], pair[1], ax=ax);<jupyter_output><empty_output><jupyter_text>This plot is a little crowded, but should make clear that you would have a very different impression of the direction of formant movement if you picked `30-70%` vs. `10-90%` or `5-95%`.  We can plot the whole time course to get the full picture:<jupyter_code>fig, ax = plt.subplots()
for f1, f2, ipa, color in zip(ou['f1'].values, ou['f2'].values,
                              ou['ipa'], colors):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color);<jupyter_output><empty_output><jupyter_text>Presumably all that fronting toward the end of the vowel is due to the following `/d/` consonant, but we can check the spectrogram just to be sure. The `scipy` and `matplotlib` modules both have spectrogram functions built in, but here we'll define a wrapper function that lets us use praat-style parameters like dynamic range in dB and window size in seconds (rather than in samples).<jupyter_code>from scipy.io import wavfile
from matplotlib.colors import LogNorm

def draw_spectrogram(x, fs, window_dur=0.005, step_dur=None, dyn_range=70, 
                     cmap=None, ax=None):
    from scipy.signal import spectrogram, gaussian
    # set default for step_dur, if unspecified
    if step_dur is None:
        step_dur = window_dur / np.sqrt(np.pi) / 8.
    # convert window & step durations to numbers of samples
    window_nsamp = int(window_dur * fs * 2)
    step_nsamp = int(step_dur * fs)
    # compute the window
    window_sigma = (window_nsamp + 1) // 6
    window = gaussian(window_nsamp, window_sigma)
    # other arguments to compute spectrogram
    noverlap = window_nsamp - step_nsamp
    kwargs = dict(noverlap=noverlap, window=window, nperseg=window_nsamp, fs=fs)
    freqs, times, power = spectrogram(x, detrend=False, mode='psd', 
                                      scaling='density', **kwargs)
    # create the figure
    if ax is None:
        fig, ax = plt.subplots()
    # default colormap
    if cmap is None:
        cmap = plt.get_cmap('Greys')
    # other arguments to the figure
    extent = (times.min(), times.max(), freqs.min(), freqs.max())
    p_ref = 20e-3  # 20 micropascals
    dB_max = 10 * np.log10(power.max() / p_ref)
    vmin = p_ref * 10 ** ((dB_max - dyn_range) / 10)
    ax.imshow(power, origin='lower', aspect='auto', cmap=cmap,
              norm=LogNorm(), extent=extent, vmin=vmin, vmax=None)
    return ax<jupyter_output><empty_output><jupyter_text>We'll use the `wavio` module to import the audio file (`scipy` also has a function `scipy.io.wavfile.read()`, but it can't handle audio with 24 bits per sample).<jupyter_code>import wavio
wav = wavio.read('audio/bVd/bowed1_11.wav')
wav_data = np.squeeze(wav.data)
sampling_frequency = wav.rate<jupyter_output><empty_output><jupyter_text>Now to plot the spectrogram:<jupyter_code>fig, ax = plt.subplots(figsize=(16, 6))
ax = draw_spectrogram(wav_data, sampling_frequency, ax=ax)
# draw an arrow
ax.annotate('F2 transition', xy=(0.45, 1650), xytext=(0.55, 2200),
            color='r', fontsize=24, fontweight='bold',
            arrowprops=dict(facecolor='r', edgecolor='r'));<jupyter_output><empty_output><jupyter_text>Yep, looks like that raising of F2 is all right at the end, where we'd expect to see the effect of the following consonant.Let's see what happens when we plot all the points for all the vowels.  For good measure, let's average the tokens together, and add a mean track for each vowel type:<jupyter_code>means = wide_data[['f1', 'f2', 'ipa']].groupby('ipa').aggregate('mean')
blank_labels = [''] * wide_data.shape[0]

fig, ax = plt.subplots()
for f1, f2, ipa, color in zip(wide_data['f1'].values, wide_data['f2'].values,
                              blank_labels, wide_data['color']):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color, alpha=0.6, linewidth=0.7);

for f1, f2, ipa, color in zip(means['f1'].values, means['f2'].values,
                              means.index, colors):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color, linewidth=2.5);<jupyter_output><empty_output><jupyter_text>Nice, but a little messy...  let's try smoothing with a moving average of length 3. When we do this, we lose the endpoints (5% and 95%) because there's nothing before/after them to average them with, so we'll remove those points.<jupyter_code>rolled = wide_data[['f1', 'f2']].rolling(axis=1, window=3, min_periods=3,
                                         center=True, win_type='boxcar').mean()
rolled['ipa'] = wide_data['ipa']
nan_columns = [(f, p) for p in (5, 95) for f in ('f1', 'f2')]
rolled.drop(nan_columns, axis=1, inplace=True)

fig, ax = plt.subplots()
for f1, f2, ipa, color in zip(rolled['f1'].values, rolled['f2'].values,
                              rolled['ipa'], wide_data['color']):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color, alpha=0.6, linewidth=0.7);<jupyter_output><empty_output><jupyter_text>Issues to watch out forOur data was fairly clean, but that won't always be the case.  Let's load some data from /hVd/ tokens from the same talker.<jupyter_code># read in the data
formant_data_hvd = pd.read_csv('diphthongs-hVd.csv')
# convert percentage strings to numbers
formant_data_hvd['pct'] = formant_data_hvd['pct'].map(lambda x: int(x.strip('%')))
# convert long to wide format
wide_data_hvd = pd.pivot_table(formant_data_hvd, index=['ipa', 'filename'],
                               columns=['formant', 'pct'], values='formant_value')
idx = wide_data_hvd.index
wide_data_hvd.reset_index(0, inplace=True, col_level=0)
wide_data_hvd.index = idx
wide_data_hvd.columns.rename(['', ''], inplace=True)
wide_data_hvd.index.rename(['', ''], inplace=True)
wide_data_hvd = wide_data_hvd[['f1', 'f2', 'ipa']]  # sort top-level
wide_data_hvd['color'] = wide_data_hvd['ipa'].replace(color_mapping)
# suppress labels for each token
blank_labels = [''] * wide_data_hvd.shape[0]

# plot tokens
fig, ax = plt.subplots()
for f1, f2, color, ipa in zip(wide_data_hvd['f1'].values,
                              wide_data_hvd['f2'].values,
                              wide_data_hvd['color'],
                              blank_labels):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color, alpha=0.5, linewidth=0.6)

# plot means
means_hvd = wide_data_hvd[['f1', 'f2', 'ipa']].groupby('ipa').aggregate('mean')
for f1, f2, ipa, color in zip(means_hvd['f1'].values,
                              means_hvd['f2'].values,
                              means_hvd.index, colors):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color, linewidth=2);<jupyter_output><empty_output><jupyter_text>ExperimentalAnother popular way of summarizing lots of tokens is to plot a ribbon showing the standard deviation around each data point.  Since these plots are not true functions, mathematically speaking (i.e., they turn back on themselves in the horizontal direction), this is actually rather trickier than we might think. Here instead we show ellipses illustrating the variation at each timepoint.<jupyter_code>def ellipse(x, y, ax, quant=0.6827, n_points=50, plot_eigvecs=False, **kwargs):
    from scipy.stats import f
    from matplotlib.patches import Polygon
    mu = (x.mean(), y.mean())
    # to reflect greater uncertainty with smaller numbers of observations,
    # we use Hotelling's T-squared instead of chi-squared
    ddf = len(x) - len(mu)
    scale_factor = len(mu) * (len(x) - 1) / ddf
    crit_rad = np.sqrt(scale_factor * f.ppf(quant, dfn=len(mu), dfd=ddf))
    # compute eigens
    cov = np.cov(x, y)
    eigval, eigvec = np.linalg.eig(cov)
    ixs = np.argsort(eigval)[::-1]
    eigval = eigval[ixs]
    eigvec = eigvec[:, ixs]
    if plot_eigvecs:
        # plot eigenvectors (major/minor axes of ellipses)
        for val, vec in zip(eigval, eigvec.T):
            ax.plot(*zip(mu, mu + vec * np.sqrt(val)), **kwargs)
            ax.plot(*zip(mu, mu - vec * np.sqrt(val)), **kwargs)
    else:
        # plot ellipses
        es = eigvec @ np.diag(np.sqrt(eigval))
        theta = np.linspace(0, 2 * np.pi, n_points)
        unit_circ = np.c_[np.cos(theta), np.sin(theta)].T
        points = (np.array(mu)[:, np.newaxis] - es @ (crit_rad * unit_circ)).T
        pgon = Polygon(xy=points, **kwargs)
        ax.add_artist(pgon)
    return ax

fig, ax = plt.subplots()

# plot ellipses
for _, df in wide_data[['f1', 'f2', 'color', 'ipa']].groupby('ipa'):
    color = df[('color', '')].unique()[0]
    for pct in np.arange(5, 100, 5):
        f1 = df[('f1', pct)].values
        f2 = df[('f2', pct)].values
        ax = ellipse(f2, f1, ax, color=color, alpha=0.1, plot_eigvecs=False)

# plot mean lines
for f1, f2, ipa, color in zip(means['f1'].values, means['f2'].values,
                              means.index, colors):
    ax = plot_diphthong(f1, f2, ipa, ax=ax, color=color);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13722923892415417, 'fraction_numerical': 0.02153000458085204, 'mean_word_length': 3.8391386953768207, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 81, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '145132', 'license': 'MIT', 'n_tokens_mistral': 5567, 'n_tokens_neox': 5327, 'n_words': 1544, 'path': 'diphthong-plotting.ipynb', 'repo_name': 'mwinn83/diphthong-plotting-tutorial', 'size': '15281'}