---
title: Vowel Formant Track Normalization Using Discrete Cosine Transform Coefficients
author:
  - name: 
      given: Josef 
      family: Fruehwald
    affiliation: University of Kentucky
    email: josef.fruehwald@uky.edu
    roles: writing
    corresponding: true
abstract: >
  This paper provides an overview of the Discrete Cosine Transform (DCT)
  as a method for smoothing vowel formant tracks, as well as a procedure
  to take any speaker normalization method that has been defined for
  formant point measurements and define an equivalent method to be applied
  directly to DCT coefficients. This procedure is followed for three established 
  normalization methods, and the difference between DCT normalization and 
  formant point normalization is found to be marginal.
keywords:
  - vowel formants
  - speaker normalization
  - vowel intrinsic spectral change
  - discrete cosine transform
acknowledgements: > 
  I would like to thank Santiago Barreda, Kevin McGowan,
  Dan Villarreal, Jack Rechsteiner, 
  the journal's anonymous reviewers, 
  and the audience at NWAV 52 for 
  feedback on this work.
date-modified: last-modified
bibliography: references.bib
citeproc: true
csl: unified-style-sheet-for-linguistics.csl
execute: 
  freeze: auto
keep-typst: true
knitr:
  opts_chunk: 
    warning: false
    message: false
    fig.width: 6
    fig.asp: 0.6
    out.width: 70%
    fig.align: center
    dev: svglite
---

```{r}
library(tidyverse)
library(reticulate)
library(ggdist)
library(showtext)
library(systemfonts)
library(geomtextpath)
library(tidynorm)
library(gt)
library(patchwork)

ptol_red = "#EE6677"
ptol_blue = "#4477AA"
ptol_lblue = "#66CCEE"
ptol_green =  "#117733"
ptol_yellow =  "#CCBB44"
ptol_purple =  "#AA3377"

# custom fonts
font_add_google("Public Sans", "Public Sans")
font_add_google("Voces", "Voces")
font_add_google("Fira Code", "Fira Code")
font_add_google("Roboto Mono", "Roboto Mono")
showtext_auto()


theme_set(
  theme_ggdist(
     base_family = "Public Sans",
     base_size = 16
  )
)

theme_no_y <- function(slab = T){

  out_theme <- theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank()
  )

  out_theme  
}

scale_y_tight <- function(...) {
  scale_y_continuous(expand = expansion(0), ...)
}

options(
  mc.cores = parallel::detectCores(),
  ggplot2.discrete.colour = c(
    lapply(
      1:6, 
      \(x) c("#4477AA", "#EE6677", "#228833", "#CCBB44", "#66CCEE", "#AA3377")[1:x]
    ),
    lapply(
      7:9,
      \(x) c("#CC5577", "#332288", "#DDCC77", "#117733","#88CCEE", "#882255", "#44AA99", "#999933", "#AA4499")[1:x]
    )
  ),
  ggplot2.discrete.fill =  c(
    lapply(
      1:6, 
      \(x) c("#4477AA", "#EE6677", "#228833", "#CCBB44", "#66CCEE", "#AA3377")[1:x]
    ),
    lapply(
      7:9,
      \(x) c("#CC5577", "#332288", "#DDCC77", "#117733","#88CCEE", "#882255", "#44AA99", "#999933", "#AA4499")[1:x]
    )
  ),
  ggplot2.ordinal.colour = \(x, ...) scico::scale_fill_scico_d(
    palette = "devon", 
    direction = -1, 
    begin = 0.1,
    end = 0.85, 
    ...
    ),
  ggplot2.ordinal.fill = \(x, ...) scico::scale_fill_scico_d(
    palette = "devon", 
    direction = -1,
    begin = 0.1,
    end = 0.85, 
    ...
  ),
  ggplot2.continuous.colour = \(x,...) scico::scale_color_scico(palette = "batlow", ...),
  ggplot2.continuous.fill = \(x,...) scico::scale_fill_scico(palette = "batlow", ...)
)


out2fig = function(out.width, out.width.default = 0.7, fig.width.default = 6) {
    fig.width.default * out.width / out.width.default 
}

```

# Introduction

Some landmark research in North American dialectology and sociolinguistics of vowels has relied on single point measurements of F1 and F2 to characterize the entire vowel [@Labov2001; @ANAE].
However, there is a wealth of work that makes it clear that full formant dynamics, or Vowel Intrinsic Spectral Change [@nearey1986], have important dialectal and sociolinguistic properties [e.g. @fox2009; @risdal2014; @docherty2015; @tanner2022].
Additionally, vowel formant extraction tools such as the FAVE suite [@rosenfelder2024], fasttrack [@barredaFastTrackFast2021; @fruehwald2024], and new-fave [@fruehwaldNewfave2025] make full formant tracks available to researchers, and curve fitting techniques are becoming more widely adopted [@soskuthyEvaluatingGeneralisedAdditive2021].

One open question about the use of full vowel formant tracks is how normalization for differences in vocal tract length between speakers ought to be carried out.
There is an *extensive* literature on the topic in the use of point measurements [see @Adank2004], but how these methods should (or could) be applied to full formant tracks has not been addressed as fully.

In this paper, I will be exploring the use the Discrete Cosine Transform (DCT) to smooth and normalize formant tracks.
The DCT was chosen, in part, because it is used in the fasttrack method to optimize the Linear Predictive Coding parameters in formant extraction, and the fasttrack python implementation allows users to save the DCT coefficients.
I will demonstrate that it is possible to carry out speaker normalization *directly* on the DCT coefficients, which can then themselves be used for quantitative analysis.

The goal of this paper is not to compare normalization techniques and evaluate their benefits and drawbacks.
Rather, it is to demonstrate that if we were to construe any given formant normalization technique as a function $f(x)$ that is applied to measured formant tracks, there is an approximately equivalent function $g(y)$ that can applied to the DCT coefficients of formant tracks.
The specific procedure for defining $g(y)$ given $f(x)$ is as follows:

::: {.block fill="luma(230)" inset="8pt" radius="4pt"}
1.  Any *non-linear* transformations (log, exponentiation, Mel, Bark, etc) must be applied to formant track data *before* the DCT is applied.
2.  For *location* parameters (values added or subtracted in the normalization):
    -   They should be estimated using the 0^th^ DCT parameter.
    -   They should be applied to just the 0^th^ DCT parameter.
3.  For *scale* parameters (values divided or multiplied by in the normalization):
    -   They should be estimated using the 0^th^ DCT parameter and multiplied by $\sqrt{2}$.[^1]
    -   They should be applied to all DCT parameters.
:::

[^1]: This $\sqrt{2}$ constant presumes the "orthogonalization" of the DCT.
    See Appendix C for more information.

::: {.content-visible when-format="html"}
:::

::: {.content-visible when-format="typst"}
:::

## Notation conventions

There are a number of mathematical formulae in this paper and mathematical notation used throughout.
I have chosen to follow the following notation conventions.

$F$, $x$

:   A formant track.
    $F$ is used when discussing formants specifically, and $x$ when discussing a signal in general.

::: {.content-visible when-format="html"}

$\bar{F}_i$

:   For a single vowel token, the average across its entire formant track.
:::

::: {.content-visible when-format="typst"}
```{=typst}
/ $macron(F)_i$: #block[
For a single vowel token, the average across its entire formant track.
]
```
:::

$y$

:   A vector of Discrete Cosine Transform coefficients (described below).

$i$

:   A variable for formant numbers.
    For example $F_1$ is the first formant, and $F_i$ is the $i^{th}$ formant.

$j$

:   A variable for an observation number.
    For example, $x_3$ would be the third value from a signal $x$.
    This will sometimes be used interchangeably to refer to the $j^{th}$ value in a signal, and sometimes to the $j^{th}$ vowel token in a data set.

$k$

:   A variable for a Discrete Cosine Transform coefficient.
    For example $y_0$ will be the 0^th^ coefficient, and $y_k$ will be the $k^{th}$ coefficient.

# The Discrete Cosine Transform

The Discrete Cosine Transform (DCT) is a general purpose curve fitting method (it is the basis of jpeg image compression [@wallace1992] ).
It is sometimes used directly on audio spectra [e.g. @zahorian1991; @zahorian1993; @guzik2007; @jannedy2017] and has also been used to model and classify vowel formant trajectories [@watson1999; @hillenbrand2001; @morrison2009; @williams2014a; @williams2015; @cox2019].
Versions of the DCT are available in the fftw R package [@fftw] (which itself calls the fftw C library) and in the scipy python package [@virtanen2020].
It has also been implemented in the emuR R package [@emuR] and the fasttrackpy python package [@fruehwald2024].

The DCT re–describes a continuous function in terms of weights on a bank of cosine functions oscillating at increasing frequencies.
There are several definitions of the DCT; for the purpose of this paper, we'll be using the definition of the DCT as used in the python scipy library [@virtanen2020] both because this is a standard and replicable method, and because this is how the fasttrackpy library implements and saves these parameters (see @sec-the-scipy-implementation-of-the-dct for more information).

The mathematical definition for the $k$th DCT is given in @eq-dctii.

$$
y_k = \frac{1}{oN} \sum_{j=0}^{N-1} x_j\cos\left(\frac{\pi k(2j + 1)}{2N} \right)
$$ {#eq-dctii}

$$
o = \begin{cases}
\sqrt{2}& \text{for }k=0\\
1 & \text{for }k>0
\end{cases}
$$

For the purpose of this paper, it's not necessary to understand every component of @eq-dctii in detail.
@fig-dct-diagram annotates the formula with the important pieces that will come into play.
A signal $x$ (in our case, a formant track) is multiplied by a cosine function of the same length, then summed.
Then, this sum is normalized (divided by the length of the signal), resulting in the DCT coefficient.
The cosine function changes depending on the number of the DCT coefficient.

![Annotated DCT Definition](assets/dct-diagram.svg){#fig-dct-diagram fig-align="center"}

@fig-dct-basis plots the 0^th^ through 2^nd^ cosine functions that the coefficients calculated in @eq-dctii weight.

```{r}
#| label: fig-dct-basis
#| fig-cap: "DCT basis functions"
scipy = import("scipy")
spdct = scipy$fft$dct

dct_basis = spdct(
  diag(100),
  norm = "backward",
  orthogonalize = TRUE
)

dct_basis[,1:3] |> 
  as_tibble() |> 
  mutate(
    x = row_number()
  ) |> 
  pivot_longer(
    -x
  ) |> 
  mutate(
    coef = (name |> 
      str_remove("V") |> 
      as.numeric())-1,
    coef_f = factor(coef)
  ) ->basis_df

basis_df |> 
  ggplot(
    aes(x, value, color = coef_f)
  )+
  geom_textline(
    linewidth = 1,
    aes(label = coef_f),
    lineend = "round"
  )+
  guides(
    color = "none"
  )+
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )
```

The inverse DCT takes the coefficients generated by @eq-dctii and returns the original signal.
The Inverse DCT equation is given in @eq-idct.

$$
x_j = \sqrt{2}y_0 + 2\sum_{k=1}^{N-1}y_k\cos\left(\frac{\pi k(2j+1)}{2N}\right)
$$ {#eq-idct}

To get the *rate* of change of a signal from its DCT coefficients, we need the first derivative of @eq-idct, given in @eq-dct-rate

$$
\frac{\delta x_j}{\delta j} = -2\sum_{k=1}^{N-1}y_k\frac{\pi k}{N}\sin\left(\frac{\pi k(2j+1)}{2N}\right)
$$ {#eq-dct-rate}

To get the *acceleration* in a signal from its DCT coefficients, we need the second derivative of @eq-idct, given in

$$
\frac{\delta^2 x_j}{\delta j^2} = -2\sum_{k=1}^{N-1}y_k\left(\frac{\pi k}{N}\right)^2\cos\left(\frac{\pi k(2j+1)}{2N}\right)
$$ {#eq-dct-accel}

# Benefits of the DCT

There are a number of general benefits to using DCT coefficients for analyzing vowel formant tracks which I will outline here.
It's worth noting at this point that several of the benefits only really hold for within-speaker analyses.
In order to apply some of these methods to multi-speaker analyses, the DCT coefficients will need to be speaker normalized (the subject of the remainder of this paper).

## Smoothing by truncation

```{r}
jean_tracks <- read_csv("data/csvs/PH06-2-1-AB-Jean_tracks.csv")

jean_tracks |> 
  filter(
    label == "ay0"
  ) -> ay0

to_plot = "0-0-10040-1"

ay0 |> 
  filter(
    id == to_plot
  ) ->ay_to_plot

f1 <- ay_to_plot$F1
```

```{r}
ay_to_plot |> 
  reframe_with_dct(
    .token_id_col = id,
    F1,
    .order = nrow(ay_to_plot)
  )->
  ay_dct

ay_dct |> 
  reframe_with_idct(
    F1, 
    .param_col = .param,
    .token_id_col = id,
    .n = .n
  )->
  ay_idct
```

The DCT returns the same number of coefficients as there are data points in the original signal.
With all of these DCT coefficients, the original signal can be fully reconstructed by applying the inverse DCT.
For example, the F1 formant track in @fig-dct-recreate consists of `r length(f1)` measurement points (the plotted points), so the application of the DCT returns 63 DCT coefficients.
When applying the *inverse* DCT to all 63 coefficients, the resulting curve (plotted in the red line) recreates the original signal.

```{r}
#| label: fig-dct-recreate
#| fig-cap: |
#|   Recreation of the original signal by applying
#|   the inverse DCT.
#| out-width: 70%
#| fig-align: center
ay_to_plot |> 
  mutate(
    x = row_number()
  ) |> 
  ggplot(
    aes(
      x, F1
    )
  )+
  geom_point(
    size = 3,
    color = "grey40"
  ) +
  geom_line(
    data = ay_idct,
    aes(
      x = .time
    ),
    linewidth = 1.5,
    color = ptol_red,
    lineend = "round"
  )+
  labs(
    x = "Sampling point"
  )
```

The first 3 DCT coefficients have relatively interpretable meanings which can be seen in their shapes in @fig-dct-basis.

-   0^th^: The overall level of the curve

-   1^st^: The amount of fall to the curve

-   2^nd^: The amount of fall-rise to the curve

In fact, several papers have been successful classifying vowel phonemes using just these first three coefficients [@watson1999; @williams2015].
Higher DCT coefficients capture higher frequency oscillations of a signal.
For example, the cosine function that the 9^th^ DCT coefficient weights is plotted in @fig-dct-ninth.

```{r}
#| label: fig-dct-ninth
#| fig-cap: |
#|   The 9^th^ dct basis function.
#| out-width: 70%
#| fig-align: center
tibble(
  x = 1:100, 
  basis = dct_basis[,10]
) |> 
  ggplot(
    aes(
      x = x,
      y = basis
    )
  )+
  geom_line(lineend = "round")+
  labs(y = "value")+
  theme(
    axis.title.x = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank()
  )
```

While these higher frequency oscillations are important for perfectly recreating the original signal, they are probably not linguistically meaningful with respect to a single vowel's formant track.
By truncating the number of DCT coefficients used to describe a formant track to the first 3 (as done in prior work) or 5 (as done in fasttrackpy), we get back a smoothed version of the formant track when applying the inverse DCT.

```{r}
#| label: fig-dct-smoothing
#| out-width: 90%
#| fig-asp: 0.5
#| fig-width: !expr out2fig(0.9)
#| fig-cap: |
#|   Greater smoothing achieved by using 
#|   fewer high frequency coefficients
c(3, 5) |> 
  map(
    ~ay_dct |> 
      filter(.param <= .x) |> 
      mutate(
        coefs = .x
      )
  ) |> 
  list_rbind() |> 
  reframe_with_idct(
    F1,
    .param_col = .param,
    .token_id_col = coefs,
    .n = 63
  ) |> 
  ggplot(
    aes(
      .time,
      F1
    )
  )+
  geom_point(
    data = ay_to_plot |> 
      mutate(
        coefs = NULL,
        .time = row_number()
      ),
    color = "grey",
    size = 1
  )+
  geom_line(
    color = ptol_red,
    lineend = "round",
    linewidth = 1
  )+
  facet_wrap(~coefs, label = "label_both") +
  labs(
    x = "sampling point"
  )
```

The benefit here is that the degree of desired smoothing can be arbitrarily specified *after* the DCT has been applied by simply using just the first $n$ coefficients.
This is in contrast with some other commonly used curve fitting methods for vowel formants, like cubic regression splines [@soskuthyGeneralisedAdditiveMixed2017a; @soskuthyEvaluatingGeneralisedAdditive2021], where the number of basis functions (and therefore coefficients) must be specified *before* the curve is fit, and cannot be altered after the fact.

## Data Compression

While this may be a diminishing concern given increasing computational power, by using just the first 5 DCT coefficients for analysis, instead of full formant track data, the volume of data that needs to be stored or loaded into memory for quantitative analysis is dramatically decreased.
As shown in @tbl-data-size, the size of the data in the Philadelphia Neighborhood Corpus (PNC) [@philly] is about an order of magnitude smaller when using DCT coefficients than when using the raw formant track data.

```{r}
#| eval: false
library(fs)

dir_ls("~/OneDrive - University of Kentucky/Corpora/PNC/derived_data/logparam_results/", glob = "*_tracks.csv") -> 
  track_files

dir_ls("~/OneDrive - University of Kentucky/Corpora/PNC/derived_data/logparam_results/", glob = "*_param.csv") -> 
  param_files


track_files |> 
  file.size() |> 
  sum()
```

```{r}
#| label: tbl-data-size
#| tbl-cap: |
#|  Comparison of data size for formant track data
#|  from the Philadelphia Neighborhood Corpus. Both
#|  number of lines in comma separated files and 
#|  total size on disk are shown.
#| 
tibble(
  data = c("tracks", "dct coefs"),
  lines = c(46973521, 4943165),
  size = c(16519237968, 999306425)
) |> 
  gt(
    rowname_col = "data"
  )  |> 
  fmt_number(
    columns = lines,
    decimals = 0
  ) |> 
  fmt_bytes(
    columns = size,
    decimals = 0
  ) |>   
  grand_summary_rows(
    fns = list(
      ratio ~ max(.)/min(.)
    ),
    fmt = ~ fmt_number(., decimals = 1)
  ) |> 
  opt_table_font(
    font = "Public Sans"
  )
```

## Averaging Formant Tracks

Averaging over formant tracks of different durations is not straightforward, as they will have a different number of sampled values.
@fig-diff-formant illustrates this problem with three example formant tracks for the vowel /ay/ from one speaker in the Philadelphia Neighborhood Corpus.
It would perhaps be inappropriate to average the final sample from a short token with the aligned midpoint sample from a longer token.
Aligning the beginnings and ends tokens in, for example, proportional time, creates a new problem that the samples in between are no longer aligned, making point-wise averaging impossible.

```{r}
#| label: fig-diff-formant
#| fig-cap: |
#|   Three F1 formant tracks of different duration.
#|   a. Formant values against actual time.
#|   b. Formant values against proportional time
#| fig-asp: 0.4
#| out-width: 100%
#| fig-width: !expr out2fig(1)
ids <- c("0-0-9693-1", "0-0-1523-1", "0-0-3868-5")

ay0 |> 
  filter(
    id %in% ids
  ) |> 
  mutate(
    .by = id,
    prop_time = (time-min(time))/diff(range(time))
  ) ->
  ay0_diff_len

ay0_diff_len|> 
  ggplot(
    aes(time, F1_s)
  )+
  geom_point(aes(color= id))+
  guides(
    color = "none"
  )+
  labs(
    y = "F1"
  )->
  p_time


ay0_diff_len|> 
  ggplot(
    aes(prop_time, F1_s)
  )+
  geom_point(aes(color= id))+
  guides(
    color = "none"
  )+
  labs(
    y = "F1",
    x = "proportional time"
  )+
  theme_no_y()->
  p_proptime

p_time + p_proptime + plot_annotation(tag_levels = 'a')
```

Instead of point-wise averaging, or fitting a smoothing model, we can instead average over the DCT coefficients for these three tokens (@tbl-dct-mean), then use these averages as DCT coefficients themselves to get the inverse DCT.
The result of this averaging process is plotted in @fig-dct-mean.

```{r}
#| label: tbl-dct-mean
#| tbl-cap: Averaging over DCT coefficients
ay0_diff_len |> 
  reframe_with_dct(
    .by = id,
    F1
  ) |> 
  select(
    id, .param, F1
  ) |> 
  pivot_wider(
    names_from = .param,
    values_from = F1
  ) |> 
  mutate(
    token = row_number(),
    .before = 1
  ) |> 
  select(-id) |> 
  gt(
    rowname_col = "token"
  ) |> 
  fmt_number(
    columns = `0`:`4`,
    decimals = 1
  ) |> 
  grand_summary_rows(
    fns = list(
      average ~ mean(.)
    ),
    fmt = ~fmt_number(., decimals = 1)
  ) |> 
  # tab_spanner(
  #   columns = `0`:`4`,label = "DCT coefficient"
  # ) |>
  tab_row_group(
    label = "token",
    rows =matches("^[0-9]")
  ) |>  
  opt_table_font(
    font = "Public Sans"
  ) 
```

```{r}
#| label: fig-dct-mean
#| fig-cap: |
#|   The result of the inverse DCT when applied to the
#|   average DCT coefficients of the three formant tracks.
ay0_diff_len |> 
  reframe_with_dct(
    .by = id,
    F1
  ) |> 
  select(
    id, .param, F1
  )  |> 
  summarise(
    .by = .param,
    F1 = mean(F1)
  ) |> 
  reframe_with_idct(
    F1,
    .param_col = .param,
    .n = 100
  ) |> 
  mutate(
    prop_time = (.time-1)/99
  )->
  ay_diff_mean

ay0_diff_len|> 
  ggplot(
    aes(prop_time, F1_s)
  )+
  geom_point(
    aes(color = id),
    alpha = 0.6
  )+
  guides(
    color = "none"
  )+
  labs(
    y = "F1",
    x = "proportional time"
  )+
  geom_textline(
    aes(
      prop_time,
      F1
    ),
    label = "average",
    data = ay_diff_mean,
    color = "grey20",
    linewidth = 1.5,
    lineend = "round"
  )
```

## Modelling with DCT coefficients

Similarly to how we can calculate averages more easily with DCT coefficients, we can also fit regression models more easily.
Let’s say, for example, we wanted to model the effect of following voicing on the entire F1 trajectory of /ay/ (plotted in @fig-ay-voicing).
To do this with the formant tracks themselves, we would need to fit a generalized additive model with an interaction between the smoothing term and voicing context, taking into account autocorrelation and, perhaps, random smooths by token [@soskuthyEvaluatingGeneralisedAdditive2021].

```{r}
jean_tracks |> 
  filter(
    label %in% c("ay", "ay0")
  ) |> 
  reframe_with_dct(
    .token_id_col = id,
    F1
  ) |> 
  select(
    .param, 
    F1, 
    label,
    id
  ) |> 
  mutate(
    .param = str_glue("y{.param}"),
    context = case_match(
      label,
      "ay0" ~ "pre-voiceless",
      "ay" ~ "elsewhere"
    )
  )  |> 
  pivot_wider(
    names_from = .param,
    values_from = F1
  )->
  ay_data
```

```{r}
#| label: fig-ay-voicing
#| fig-asp: 0.5
#| out-width: 90%
#| fig-width: !expr out2fig(0.9)
#| fig-cap: |
#|   The effect of voicing on F1 formant tracks
#|   for one speaker.
jean_tracks |> 
  filter(
    label %in% c("ay", "ay0")
  ) |> 
  mutate(
    .by = id,
    context = case_match(
      label,
      "ay0" ~ "pre-voiceless",
      "ay" ~ "elsewhere"
    ),
    prop_time = (time-min(time))/diff(range(time))
  ) |> 
  ggplot(
    aes(
      prop_time,
      F1,
      color = context
    )
  )+
  stat_smooth(lineend = "round")+
  labs(
    x = "proportional time"
  )
```

With DCT coefficients, on the other hand, we could fit a linear model, using each individual DCT coefficient as the outcome, and the voicing effect as the predictor.
This is a form of functional regression [@Ramsay2006; @gubian2015].
@lst-dctreg shows the R code involved for fitting such a model, with the model results summarized in @tbl-dctreg.

```{r}
#| echo: true
#| lst-label: lst-dctreg
#| lst-cap: A model
ay_model <- lm(
  cbind(y0, y1, y2, y3, y4) ~ context,
  data = ay_data
)
```

```{r}
#| label: tbl-dctreg
#| tbl-cap: |
#|    Regression coefficients over DCTs
ay_model |> 
  broom::tidy() |> 
  select(
    response, term, estimate
  ) |> 
  mutate(
    term = str_remove(term, "context")
  ) |> 
  pivot_wider(
    names_from = response,
    values_from = estimate
  ) |> 
  gt(
    rowname_col = "term"
  )  |> 
  fmt_number(
    decimals = 1
  ) |> 
  # tab_spanner(
  #   columns = starts_with("y"),
  #   label = "DCT Coefficient"
  # ) |> 
  opt_table_font(
    font = "Public Sans"
  )
```

In @tbl-dctreg, the intercept terms for each DCT coefficient correspond to the average of that coefficient across the reference level (the elsewhere context), and the pre-voiceless terms correspond to the *difference* of the pre-voiceless context from the reference level for each coefficient.

Interestingly, if we apply the inverse DCT to these difference terms, the result is the difference *curve* between the the two voicing contexts, plotted in @fig-diffcurve.

```{r}
#| label: fig-diffcurve
#| fig-cap: |
#|   The fitted difference curve between pre-voiceless
#|   /ay/ and /ay/ elsewhere.

set.seed(2024)
cov_mat <- vcov(ay_model)
keeper <- str_detect(
  rownames(cov_mat),
  "context"
)

#context_cov <- cov_mat[keeper, keeper]
#intercept_cov <- cov_mat[!keeper, !keeper]

sims <- MASS::mvrnorm(n = 5000, mu=as.vector(coef(ay_model)), Sigma = cov_mat)


sims |> 
  as_tibble() |> 
  mutate(
    rep = row_number()
  ) |> 
  pivot_longer(
    starts_with("y"),
    values_to = "F1"
  ) |> 
  separate_wider_delim(
    name, 
    delim = ":",
    names = c("param", "term")
  ) |> 
  filter(
    str_detect(term, "context")
  ) |> 
  reframe_with_idct(
    F1,
    .token_id_col = rep,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    prop_time = (.time - 1)/49
  )->
  fit_sim

fit_sim |> 
  ggplot(
    aes(
      prop_time,
      F1
    )
  )+
  stat_lineribbon(.width = 0.95, lineend = "round")+
  labs(
    x = "proportional time",
    y = "pre-voiceless difference (F1)"
  )+
  guides(
    fill = "none"
  )
```

More complex models exploring the difference in formant dynamics between these two allophones could be specified, but by using DCT coefficients as the outcome variables, we don't need to include time across the formant track in our predictors.

## Rate of change analysis

Given the definition of the first derivative of the signal in @eq-dct-rate, it's also possible to analyze the rate of change in formant dynamics.
@fig-rate plots the predicted rate of change for these two allophones of /ay/ based on the model fit above.
These rate of change curves were calculated directly from the DCT coefficients by plugging them into @eq-dct-rate.
This kind of analysis could be useful in research on hypo/hyper-articulation, how coarticulation influences sound change, etc.

```{r}
diff_idct <- function(y, n = 100){
  N = length(y)
  j = seq(0, n - 1, length = n)
  
  sin_fun <- function(j, k, n){
    sin(
      (pi*k*((2*j)+1))/(2*n)
    )
  }
  
  sum_fun <- function(j, y, n){
    y1 <- y[-1]
    k = seq_along(y1)

    midterm <- ((pi * k)/n)
     
    sum(y1 * midterm * sin_fun(j, k, n)) * -2
  }
  
  map_vec(
    j,
    ~sum_fun(.x, y, n)
  ) * (n/100)
}
```

```{r}
sims |> 
  as_tibble() |> 
  mutate(
    rep = row_number()
  ) |> 
  pivot_longer(
    starts_with("y"),
    values_to = "F1"
  ) |> 
  separate_wider_delim(
    name, 
    delim = ":",
    names = c("param", "term")
  )  |> 
  reframe(
    .by = c(rep, param),
    F1_p = c(F1[1], sum(F1)),
    label = c("ay", "ay0")
  ) |> 
  reframe(
    .by = c(rep, label),
    F1 = idct(F1_p, n = 50),
    F1_rate = diff_idct(F1_p, n = 50),
    time = (0:49)/49
  ) ->
  rate_df
```

```{r}
#| label: fig-rate
#| fig-cap: |
#|   The rate of change of /ay/ allophones. The rate
#|   is in expressed in terms of Hz per 0.01 change in
#|   proportional time.
#| fig-asp: 0.5
#| out-width: 90%
#| fig-width: !expr out2fig(0.9)
rate_df  |> 
  mutate(
    context = case_match(
      label,
      "ay" ~ "elsewhere",
      "ay0" ~ "pre-voiceless"
    )
  ) |> 
  ggplot(
    aes(time, F1_rate)
  )+
  geom_hline(
    yintercept = 0
  )+
  stat_summary(
    fun = mean,
    geom = "line",
    aes(color = context),
    linewidth = 1,
    lineend = "round"
  )+  
  stat_ribbon(
    aes(
      fill = context,
      color = context
    ),
    alpha = 0.3,
    .width = 0.95 
  )+
  labs(
    x = "proportional time",
    y = "F1 change (Hz/1%)"
  )
```

## Multi-speaker analysis

The ability to regress over DCT coefficients could be invaluable for research into, for example, diachronic changes in formant dynamics.
However this modeling approach, along with the other benefits of DCT coefficients discussed above, is really only well defined for within-speaker analyses.
For example, @fig-mean-unnorm plots formant trajectories for the vowels /iy/, /ey/, /ay/, /uw/, /ow/ and /aw/ from two speakers in the PNC.
Specifically, DCT coefficients from each vowel class for each speaker were averaged, and then the inverse DCT applied to result in the plotted trajectories.

```{r}
fs::dir_ls("data/csvs/", glob = "*_param.csv") |> 
  map(
    read_csv
  ) |> 
  list_rbind()->
  all_param
```

```{r}
#| label: fig-mean-unnorm
#| fig-asp: 0.5
#| out-width: 90%
#| fig-width: !expr out2fig(0.9)
#| fig-cap: |
#|   Average vowel trajectories from unnormalized
#|   DCT coefficients.
all_param |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw")
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1:F2,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1:F2,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ]
  ) |> 
  ggplot(
    aes(F2, F1, color = speaker)
  ) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.2, "cm")
    ),
    size = 5,
    linewidth = 1,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse()+
  scale_x_reverse()
```

The consistent differences in vowel space size and location between these two speakers are also reflected in their DCT coefficients, such that any analysis combining their data will face the same problems as using unnormalized formant point measurements.
For example, Speaker B's /iy/ and /ey/ appear to be backer than Speaker A's, but it's unclear how much of this is due to their overall smaller vowel space size versus having a different vowel space target.
In order to conduct multi-speaker analyses using DCT coefficients, the DCT coefficients themselves will need to be normalized.

# Transforming DCT Coefficients

Nearly all speaker normalization methods involve:

-   Shifting the *location* of formant values by a constant $l$.

-   *Scaling* formant values by a constant $s$.

$$
F' = s(F+l)
$$ {#eq-norm-abstract}

If there are coefficients $y_k$ that when the inverse DCT is applied to them return $F$, then there are coefficients $y'_k$ that when the inverse DCT is applied to them return $F'$.
The goal of this section is to identify whether and how it is possible to directly transform $y_k$ to $y'_k$.
To achieve this, we'll begin by identifying how shifting and scaling the input signal by our desired $l$ and $s$ affect the resulting DCT coefficients.

## Useful Equivalencies

In order to simplify the notation in the following section, we'll define some equivalencies.
[^2] The cosine term in @eq-dctii will be replaced with $\psi_j$, a function over the signal index $j$ and the DCT coefficient index $k$.

[^2]: I would like to thank a reviewer for the suggested notation for this section.

$$
\psi_j(k) = \cos\left(\frac{\pi k(2j + 1)}{2N}\right)
$$ {#eq-psi}

The normalizing constant will be replaced with $c$.

$$
c = \frac{1}{oN}
$$ {#eq-norm}

With these definitions, the DCT can be expressed as

$$
y_k = c\sum_{j=0}^{N-1}x_j\psi_j(k)
$$ {#eq-dct-phi}

There are some additional important equivalencies we'll leverage below.
The first is for when $k=0$.

$$
1 = \psi_j(0)
$$ {#eq-psi1}

Therefore, summing $\psi_j(0)$ over $j$ will equal $N$.

$$
N = \sum_{j=0}^{N-1}\psi_j(0)
$$ {#eq-psi-k0}

Second, when $k>0$, summing $\psi_j(k)$ over $j$ will equal 0.
$$
0 = \sum_{j=0}^{N-1}\psi_j(k);~ k > 0
$$ {#eq-psi-k1}

## Scaling the input

Let's say we have a signal $x$ that we applied the DCT to to get $y_k$.

$$
y_k = c\sum_{j=0}^{N-1}x_j\psi_j(k)
$$ {#eq-dct-orig}

If our normalization method called for scaling $x$ by $s$, our input to the DCT would now be $sx_j$.
The new DCT coefficients $y'_k$ would be

$$
y'_k = c\sum_{j=0}^{N-1}(sx_j)\psi_j(k)
$$ {#eq-dct-scale0}

As a constant, $s$ can be brought outside the summation.

$$
y'_k = sc\sum_{j=0}^{N-1}x_j\psi_j(k)
$$

Now, everything to the right of $s$ is the same as @eq-dct-orig, so we can substitute in $y_k$,

$$
y'_k = sy_k
$$ {#eq-dct-scaling}

### Scaling conclusion

The result here is that any scaling term applied to an input signal results in the same scaling term being applied to all of its DCT coefficients.

## Shifting the input

Again, let's say we have a signal $x$ that we applied the DCT to to get $y_k$.
If our normalization method called for adding the constant $l$ to $x$, our input to the DCT would now be $(x_j + l)$.
The new DCT coefficients $y'_k$ would now be

$$
y'_k = c\sum_{j=0}^{N-1}(x_j+l)\psi_j(k)
$$ {#eq-shift0}

The summation term can be rewritten as the addition of two summation terms.

$$
y'_k = c(\sum_{j=0}^{N-1}x_j\psi_j(k) +\sum_{j=0}^{N-1}l\psi_j(k))
$$ {#eq-shift1}

Then, we can distribute the normalizing term $c$ and move the constant $l$ outside of its summation.

$$
y'_k = c\sum_{j=0}^{N-1}x_j\psi_j(k) + lc\sum_{j=0}^{N-1}\psi_j(k))
$$ {#eq-shift2}

Everything to the left of the addition is the same as @eq-dct-orig, so we can substitute in $y_k$.

$$
y'_k = y_k + lc\sum_{j=0}^{N-1}\psi_j(k)
$$ {#eq-shift3}

We can now use the equivalencies in @eq-psi-k0 and @eq-psi-k1 to evaluate @eq-shift3 when $k=0$ and when $k>0$.

For $k=0$, we'll substitute both the sum of $\psi_j(0)$ (@eq-psi-k0) and the normalizing constant $c$ (@eq-norm) to arrive at

$$
y'_0 = y_0 + l\frac{1}{\sqrt{2}N}N = y_0 + \frac{l}{\sqrt{2}}
$$ {#eq-shift-k0}

When $k>0$, the sum of $\psi_j(k)$ is 0 (@eq-psi-k1), which eliminates the second term entirely.

$$
y'_{k>0} = y_{k>0} + 0lc = y_{k>0}
$$ {#eq-shift-k1}

### Shifting conclusion

If an input signal is shifted by $l$, its 0^th^ DCT coefficient will be shifted by $\frac{l}{\sqrt{2}}$, and its remaining DCT coefficients will remain the same.

## Non-linear transforms of the input {#sec-non-linear-transforms-of-the-input}

Some speaker normalization procedures involve transformations of the signal that are non-linear.
The most common one is taking the logarithm of formant values.
If we logged an input signal $x$, its DCT coefficients would be given as

$$
y'_k = c\sum_{j=0}^{N-1}\log(x_j)\psi_j(k)
$$ {#eq-dct-log}

There are no additional equivalences we can leverage here to re-express $y'_k$ as a transformation of $y_k$.
In order to calculate $y'_k$, we must directly transform the original signal, then apply the DCT.

## DCT Coefficient Transformation Summary

In summary, if we want to shift our formant values by some value $l$ and scale them by some value $s$:

-   We need to add $\frac{l}{\sqrt{2}}$ to the 0^th^ DCT coefficient (@eq-shift-k0).

-   We don't add anything to the remaining DCT coefficients (@eq-shift-k1).

-   We need to scale all DCT coefficients by $s$ (@eq-dct-scaling).

Additionally, for non-linear transformations of formant frequencies, there is no way to transform the DCT coefficients of the original values directly into DCT coefficients of the transformed values.
Instead, the entire DCT needs to be applied to the transformed formants directly.

# Estimating normalization parameters

Now that we understand how to shift DCT coefficients by $l$ and scale them by $s$, we need to establish how to get these $l$ and $s$ normalization parameters from DCT coefficients.
The two most common normalization parameters are the mean and standard deviation over vowels, so we will focus on those.

## Calculating means

Let's begin with $F$ as a formant track for a single vowel token, and $y_k$ as the coefficients we get from applying the DCT to $F$.
As we already established in @eq-psi1, $\psi_j(0) =1$, which means the 0^th^ DCT coefficient can be given as

$$
y_0 = \frac{1}{\sqrt{2}}\frac{1}{N}\sum_{j=0}^{N-1}F_j
$$ {#eq-zero-mean0}

::: {.content-visible when-format="html"}
The sum of all values in a vector divided by the number of values in the vector is the mean of that vector, which we can represent with $\bar{F}$.
So, the 0^th^ coefficient is equal to the mean of $F$, divided by $\sqrt{2}$.
:::

::: {.content-visible when-format="typst"}
```{=typst}
The sum of all values in a vector divided by the number of values in the vector is the mean of that vector, which we can represent with $macron(F)$. So, the 0#super[th] coefficient is equal to the mean of $F$, divided by $sqrt(2)$.
```
:::

::: {.content-visible when-format="html"}
$$
y_0 = \frac{1}{\sqrt{2}}\bar{F}
$$
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ y_0 = 1 / sqrt(2) macron(F) $ ])<eq-zeroth-mean>
```
:::

And therefore

::: {.content-visible when-format="html"}
$$
\bar{F} = \sqrt{2}y_0
$$ {#eq-mean-zeroth}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ macron(F) = sqrt(2) y_0 $ ])<eq-mean-zeroth>
```
:::

The $\sqrt{2}$ constant will carry through any functions that are weighted sums.
So if the function $f(x)$ is a weighted sum,

::: {.content-visible when-format="html"}
$$
f(\bar{F}) = \sqrt{2}f(y_0)
$$ {#eq-y0-linear}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ f(macron(F)) = sqrt(2)f(y_0) $ ])<eq-y0-linear>
```
:::

The mean is a weighted sum, so

::: {.content-visible when-format="html"}
$$
\text{mean}(\bar{F}) = \sqrt{2}~\text{mean}(y_0)
$$ {#eq-mean-y0}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ upright("mean") (macron(F)) = sqrt(2) med upright("mean") (y_0) $ ])<eq-mean-y0>
```
:::

### Shifting by the mean

::: {.content-visible when-format="html"}
Let's use $\text{mean}(\bar{F})$ as our desired location parameter $l$ by which to shift all vowel tokens.
We've already established in @eq-shift-k0 and @eq-shift-k1 that this will involve only changes to the 0^th^ DCT coefficient for the $j^{th}$ token.
If we substitute @eq-mean-y0 into @eq-shift-k0, we get
:::

::: {.content-visible when-format="typst"}
```{=typst}
Let’s use $upright("mean")(macron(F))$ as our desired location parameter $l$ by which to shift all vowel tokens. We’ve already established in @eq-shift-k0 and @eq-shift-k1 that this will involve only changes to the 0#super[th] DCT coefficient for the $j^(t h)$ token. If we substitute @eq-mean-y0 into @eq-shift-k0, we get
```
:::

$$
y'_{0j} = y_{0j} - \frac{\sqrt{2}~\text{mean}(y_0)}{\sqrt{2}} = y_{0j}-\text{mean}(y_0)
$$ {#eq-y0-meancenter}

Because $\sqrt{2}$ cancels out in @eq-y0-meancenter, to center all DCT coefficients on the mean across tokens, we just need to center $y_0$ on the mean across $y_0$.

## Calculating standard deviations

::: {.content-visible when-format="html"}
The standard deviation over $\bar{F}$ is *not* a weighted sum, but as it happens for the standard deviation specifically, there is still a straightforward relationship between formant values and $y_0$.
:::

::: {.content-visible when-format="typst"}
```{=typst}
The standard deviation over $macron(F)$ is #emph[not] a weighted sum, but as it happens for the standard deviation specifically, there is still a straightforward relationship between formant values and $y_0$.
```
:::

::: {.content-visible when-format="html"}
Firstly, the standard deviation across $\bar{F}$ can be given as
:::

::: {.content-visible when-format="typst"}
```{=typst}
Firstly, the standard deviation across $macron(F)$ can be given as
```
:::

::: {.content-visible when-format="html"}
$$
\text{sd}(\bar{F}) = \sqrt{
  \frac{
  \sum(\bar{F}_j - \text{mean}(\bar{F}))^2
  }{
    N
  }
}
$$ {#eq-sd-basic}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ upright("sd")(macron(F)) = sqrt(frac(sum (macron(F)_j - upright("mean")(macron(F)))^2, N)) $ ])<eq-sd-basic>
```
:::

::: {.content-visible when-format="html"}
Substituting in @eq-mean-zeroth and @eq-mean-y0, we have
:::

::: {.content-visible when-format="typst"}
```{=typst}
Substituting in @eq-mean-zeroth and @eq-mean-y0, we have
```
:::

::: {.content-visible when-format="html"}
$$
\text{sd}(\bar{F}) = \sqrt{
  \frac{
  \sum(\sqrt{2}y_{0j} - \sqrt{2}~\text{mean}(y_0))^2
  }{
    N
  }
}
$$ {#eq-y0_sd0}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ upright("sd")(macron(F)) = sqrt(frac(sum (sqrt(2) y_(0 j) - sqrt(2) med upright("mean")(y_0))^2, N)) $ ])<eq-y0_sd0>
```
:::

This can be simplified in a few steps.

::: {.content-visible when-format="html"}
$$
\text{sd}(\bar{F}) = \sqrt{
  \frac{
  \sum\sqrt{2}^2(y_{0j} - \text{mean}(y_0))^2
  }{
    N
  }
}
$$ {#eq-y0-simplify}

$$
\text{sd}(\bar{F}) = \sqrt{
  \frac{
  \sum2(y_{0j} - \text{mean}(y_0))^2
  }{
    N
  }
}
$$

$$
\text{sd}(\bar{F}) = \sqrt{
  2\frac{
  \sum(y_{0j} - \text{mean}(y_0))^2
  }{
    N
  }
}
$$

$$
\text{sd}(\bar{F}) = \sqrt{2}\sqrt{
  \frac{
  \sum(y_{0j} - \text{mean}(y_0))
  }{
    N
  }
}
$$
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ upright("sd")(macron(F)) = sqrt(frac(sum sqrt(2)^2 (y_(0 j) - upright("mean")(y_0))^2, N)) $ ])<eq-y0-simplify>

$ upright("sd")(macron(F)) = sqrt(frac(sum 2 (y_(0 j) - upright("mean")(y_0))^2, N)) $

$ upright("sd")(macron(F)) = sqrt(2 frac(sum (y_(0 j) - upright("mean")(y_0))^2, N)) $

$ upright("sd")(macron(F)) = sqrt(2) sqrt(frac(sum (y_(0 j) - upright("mean")(y_0)), N)) $
```
:::

Resulting in

::: {.content-visible when-format="html"}
$$
\text{sd}(\bar{F}) = \sqrt{2}~\text{sd}(y_0)
$$ {#eq-y0-sd-final}
:::

::: {.content-visible when-format="typst"}
```{=typst}
#math.equation(block: true, numbering: "(1)", [ $ upright("sd")(macron(F)) = sqrt(2) med upright("sd")(y_0) $ ])<eq-y0-sd-final>
```
:::

### Scaling by the standard deviation

::: {.content-visible when-format="html"}
Let's use $\frac{1}{\text{sd}(\bar{F})}$ as our scaling parameter $s$ by which to scale all vowel tokens.
By substituting @eq-y0-sd-final into @eq-dct-scaling, we get
:::

::: {.content-visible when-format="typst"}
```{=typst}
Let’s use $frac(1, upright("sd")(macron(F)))$ as our scaling parameter $s$ by which to scale all vowel tokens. By substituting @eq-y0-sd-final into @eq-dct-scaling, we get
```
:::

$$
y'_{kj} = \frac{1}{\sqrt{2}~\text{sd}(y_0)}y_{kj}
$$ {#eq-sd-shift}

So, when scaling by the standard deviation, or by a linear function over $y_0$, we need to multiply the result by $\sqrt{2}$, and then use it as a scaling parameter over all DCT coefficients.

# Example Application

In this section, we'll take three normalization methods that were defined over formant point measurements, and redefine DCT normalization methods as outlined above.
These three methods were chosen because they cover all three possibilities of shifting and scaling.

| Method                   |            $l$ (shift) |                $s$ (scale) |
|------------------------|---------------------:|-------------------------:|
| $\Delta$F [@johnson2020] |                    $0$ |      $\frac{1}{\delta(F)}$ |
| Nearey [@nearey1978]     | $-\text{mean}(\log F)$ |                        $1$ |
| Lobanov [@lobanov1971]   |    $-\text{mean}(F_i)$ | $\frac{1}{\text{sd}(F_i)}$ |

: Example Normalization Methods {#tbl-methods}

We'll be using @eq-norm-abstract (repeated here) as the generalized normalization equation that all three methods have in common.

$$
F' = s(F+l)
$$ {#eq-norm2}

## $\Delta$F

@johnson2020 proposes a vocal tract length normalization method that scales all formants, but does not shift them.
In the original paper, a quantity $\Delta$F is defined, but for the purpose of this paper, we'll define it as a function $\delta(F)$.

$$
\delta(F) =\sum_{i=1}^{M}\sum_{j=1}^{N} \frac{F_{ij}}{i-0.5}
$$ {#eq-deltaf-def}

Where $i$ is the formant number.
The scaling and shifting parameters $s$ and $l$ for this method are $$
s = \frac{1}{\delta(F)}; ~~l=0
$$ {#eq-deltaf-sl}

Plugging these parameters in to @eq-norm2, we get\
\
$$
F'_{ij} = \frac{1}{\delta(F)}F_{ij}
$$ {#eq-deltaf-norm}

::: {.content-visible when-format="html"}
The $\delta(F)$ function is a weighted sum, so by @eq-y0-linear we get
:::

::: {.content-visible when-format="typst"}
```{=typst}
The $delta (F)$ function is a weighted sum, so by @eq-y0-linear we get
```
:::

$$
\delta(F) = \sqrt{2}\delta(y_0)
$$ {#eq-deltaf-y0}

Since this method involves only scaling, we can plug @eq-deltaf-y0 into @eq-dct-scaling to get the DCT coefficients of the normalized formants as

$$
y'_{kij} = \frac{1}{\sqrt{2}\delta(y_{0})}y_{kij}
$$ {#eq-deltaf-dct}

I used @eq-deltaf-dct to normalize the DCT coefficients for F1 through F3 for two speakers from the PNC.
@fig-deltaf compares the resulting inverse DCT on the unnormalized vs $\Delta$F normalized formant tracks, focusing on just the long vowels.

```{r}
all_param |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw")
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1:F2,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1:F2,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ],
    norm = "unnormalized"
  ) ->
  dct_unnorm
```

```{r}
all_param |> 
  norm_dct_generic(
    F1:F3,
    .by = file_name,
    .token_id_col = id,
    .param_col = param, 
    .S = mean(.formant/(.formant_num-0.5))
  ) |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw")
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1_n:F2_n,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1_n:F2_n,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ],
    norm = "delta F"
  )  ->
  dct_deltaf
```

```{r}
ggplot(
  dct_unnorm,
  aes(F2, F1, color = speaker)
) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "cm")
    ),
    linewidth = 1,
    size = 3,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse()+
  scale_x_reverse()+
  labs(
    title = "unnormalized"
  )->
  p_unnorm


ggplot(
  dct_deltaf,
  aes(F2_n, F1_n, color = speaker)
) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "cm")
    ),
    linewidth = 1,
    size = 3,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse("F1 (normalized)")+
  scale_x_reverse("F2 (normalized)")+
  labs(
    title = "delta F normalized"
  )->
  p_deltaf
```

```{r}
#| label: fig-deltaf
#| fig-cap: |
#|   Unnormalized vs Delta F normalized formant tracks.
#| out-width: 100%
#| fig-asp: 0.4
#| fig-width: !expr out2fig(1)

p_unnorm + plot_spacer() + p_deltaf + plot_layout(guides = "collect", widths = c(5, 1, 5))
```

## Nearey

Nearey normalization [@nearey1978] involves a single shift in location, and has been recently argued to more closely reflect listeners perceptual categorization [@barreda2021].

Importantly, the first step of Nearey normalization involves taking the logarithm of formant values.
As discussed in @sec-non-linear-transforms-of-the-input, there is no way to directly transform DCT coefficients of formants in Hz to DCT coefficients in log(Hz).
Instead, the formants must first be logged, and then have the DCT applied.
For the sake of notational brevity, for this section assume that any $F$ refers to a log transformed formant, and any $y_k$ refers to the DCT coefficients of that log transformed formant.

The scaling ($s$) and shifting ($l$) parameters for Nearey normalization are:

$$
s = 1; ~~l = -\text{mean}(F)
$$ {#eq-nearey-param}

Plugging this into @eq-norm2, we get

$$
F'_{ij} = F_{ij} - \text{mean}(F)
$$ {#eq-nearey-norm}

Since this involves shifting by the mean, we can use @eq-y0-meancenter to define the DCT normalization like so:

$$
y'_{0ij} = y_{0ij} - \text{mean}(y_{0})
$$ {#eq-nearey-dct}

$$
y'_{k>0ij} = y_{k>0ij}
$$

The 0^th^ DCT parameter of every token will be shifted by the mean across $y_0$, and all other DCT parameters will stay the same.

Using @eq-nearey-dct, I Nearey normalized the DCT coefficients for F1 through F3 for the same two speakers from the PNC as before.
@fig-nearey compares the resulting inverse DCT on the unnormalized vs normalized formant tracks.

```{r}
fs::dir_ls("data/csvs/", glob = "*_logparam.csv") |> 
  map(
    read_csv
  ) |> 
  list_rbind()->
  all_logparam
```

```{r}
all_logparam |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw"),
    fol_seg != "L"
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1:F2,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1:F2,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ],
    norm = "unnormalized"
  ) ->
  dct_log_unnorm

all_logparam |> 
  norm_dct_generic(
    F1:F3,
    .by = file_name,
    .token_id_col = id,
    .param_col = param,
    .L = mean(.formant)
  ) |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw"),
    fol_seg != "L"
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1_n:F2_n,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1_n:F2_n,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ],
    norm = "unnormalized"
  ) ->
  dct_nearey
```

```{r}
ggplot(
  dct_log_unnorm,
  aes(F2, F1, color = speaker)
) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "cm")
    ),
    linewidth = 1,
    size = 3,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse("log(F1)")+
  scale_x_reverse("log(F2)")+
  labs(
    title = "unnormalized"
  )->
  p_log_unnorm


ggplot(
  dct_nearey,
  aes(F2_n, F1_n, color = speaker)
) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "cm")
    ),
    linewidth = 1,
    size = 3,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse("log(F1) (normalized)")+
  scale_x_reverse("log(F2) (normalized)")+
  labs(
    title = "Nearey normalized"
  )->
  p_nearey
```

```{r}
#| label: fig-nearey
#| fig-cap: |
#|   Unnormalized vs Nearey normalized formant tracks.
#| out-width: 100%
#| fig-asp: 0.4
#| fig-width: !expr out2fig(1)
p_log_unnorm + plot_spacer() + p_nearey + plot_layout(guides = "collect", widths = c(5, 1, 5))
```

## Lobanov

When applied to other domains of data analysis, Lobanov normalization [@lobanov1971] is called z-scoring.
The scaling ($s$) and shifting ($l$) parameters for Lobanov normalization are:

$$
s = \frac{1}{\text{sd}(F_i)}; ~~ l=-\text{mean}(F_i)
$$ {#eq-lobanov-param}

Plugging these parameters into @eq-norm2, we get the following definition of Lobanov normalization.

$$
F'_{ij} = \frac{1}{\text{sd}(F_i)}(F_{ij}-\text{mean}(F_i))
$$ {#eq-lobanov-norm}

::: {.content-visible when-format="html"}
The shifting parameter will only affect the 0^th^ coefficient, and the scaling parameter will affect all coefficients.
Plugging the equivalencies from @eq-y0-meancenter and @eq-y0-sd-final into @eq-lobanov-norm, we get the following definition for Lobanov normalizing the DCT coefficients.
:::

::: {.content-visible when-format="typst"}
```{=typst}
The shifting parameter will only affect the 0#super[th] coefficient, and the scaling parameter will affect all coefficients. Plugging the equivalencies from @eq-y0-meancenter and @eq-y0-sd-final into @eq-lobanov-norm, we get the following definition for Lobanov normalizing the DCT coefficients.
```
:::

$$
y_{0ij} = \frac{1}{\sqrt{2}~\text{sd}(y_{0i})}(y_{0ij} - \text{mean}(y_{0i}))
$$ {#eq-lobanov-dct}

$$
y_{k>0;ij} = \frac{1}{\sqrt{2}~\text{sd}(y_{0i})}y_{k>0;ij}
$$

Again, using @eq-lobanov-dct, I Lobanov normalized the DCT coefficients for F1 and F2 for the same two speakers from the Philadelphia neighborhood corpus.
@fig-lobanov compares the resulting average formant tracks when applying the inverse DCT to the normalized coefficients.

```{r}
all_param |> 
  norm_dct_generic(
    F1:F3,
    .by = file_name,
    .token_id_col = id,
    .param_col = param,
    .L = mean(.formant),
    .S = sd(.formant),
    .by_formant = TRUE
  ) |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw")
  ) |> 
  summarise(
    .by = c(param, file_name, label),
    across(
      F1_n:F2_n,
      mean
    )
  ) |> 
  reframe_with_idct(
    F1_n:F2_n,
    .by = file_name,
    .token_id_col = label,
    .param_col = param,
    .n = 50
  ) |> 
  mutate(
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ],
    norm = "unnormalized"
  ) ->
  dct_lobanov
```

```{r}
ggplot(
  dct_lobanov,
  aes(F2_n, F1_n, color = speaker)
) +
  geom_textpath(
    aes(
      label = label,
      group = interaction(label, speaker)
    ),
    arrow = arrow(
      type = "closed",
      length = unit(0.1, "cm")
    ),
    linewidth = 1,
    size = 3,
    lineend = "round",
    vjust = "auto", hjust = 0.2,
    key_glyph = draw_key_rect
  )+
  scale_y_reverse("F1 (normalized)")+
  scale_x_reverse("F2 (normalized)")+
  labs(
    title = "Lobanov normalized"
  )->
  p_lobanov
```

```{r}
#| label: fig-lobanov
#| fig-cap: |
#|   Unnormalized vs Lobanov normalized formant tracks.
#| out-width: 100%
#| fig-asp: 0.4
#| fig-width: !expr out2fig(1)
p_unnorm + plot_spacer() + p_lobanov + plot_layout(guides = "collect", widths = c(5, 1, 5))
```

```{r}
fs::dir_ls("data/csvs/", glob = "*_tracks.csv") |> 
  map(
    read_csv
  ) |> 
  list_rbind()->
  all_tracks
```

```{r}
all_tracks |> 
  slice(
    .by = c(file_name, id),
    round(n()/2)
  )->
  track_means
```

## Comparison to point-based normalization

It's visually clear that none of the normalization methods resulted in complete overlap of the formant trajectories for these two speakers.
However, the goal of speaker normalization is to eliminate (socio)linguistically irrelevant differences, while preserving meaningful differences.
While different approaches have been taken to evaluating how well any given normalization method achieves this goal [@Adank2004; @barreda2021] the goal of this paper is just to define formant *track* normalization based on formant *point* normalization methods.
As such, we'll briefly compare the formant track results above to the same point-based methods.

For this comparison, the same two speakers' formant track data was used, and formant values at 50% of the duration taken as the point measurement.
These point measures were then normalized using the $\Delta$F, Nearey, and Lobanov methods, as defined above.
Then, for the same set of vowels, the average F1 and F2 was estimated.
@fig-point-norm plots the result of all three methods as well as the unnormalized data.

```{r}
#| label: fig-point-norm
#| fig-cap: |
#|   Comparison of three different formant point normalization
#|   methods to unnormalized formants.
#| fig-asp: 0.83
#| out-width: 90%
#| fig-width: !expr out2fig(0.9)
track_means |> 
  select(
    -(F1_s:F3_s)
  ) |> 
  norm_lobanov(
    c(F1, F2),
    .by = file_name
  ) |> 
  norm_deltaF(
    c(F1, F2, F3),
    .by = file_name
  ) |> 
  norm_nearey(
    c(F1, F2, F3),
    .by = file_name
  ) |> 
  mutate(
    F1_unnorm = exp(F1),
    F2_unnorm = exp(F2),    
    F3_unnorm = exp(F3)
  ) |> 
  select(
    -(F1:F3)
  ) |> 
  filter(
    label %in% c("iy", "ey", "ay", "ow", "aw", "uw")
  ) |> 
  pivot_longer(
    matches("^F\\d_")
  ) |> 
  separate_wider_delim(
    name,
    delim = "_",
    names = c("formant", "method")
  ) |> 
  pivot_wider(
    names_from = formant,
    values_from = value
  ) |> 
  summarise(
    .by = c(file_name, label, method),
    across(
      F1:F2,
      mean
    )
  ) |> 
  mutate(
    method = case_match(
      method,
      "unnorm" ~ "unnormalized",
      "df" ~ "delta F",
      "lm" ~ "Nearey",
      "z" ~ "Lobanov"
    ) |> 
      factor(
        levels = c(
          "unnormalized",
          "delta F",
          "Nearey",
          "Lobanov"
        )
      ),
    speaker = c("A", "B")[
      as.numeric(as.factor(file_name))
    ]
  ) |> 
  ggplot(
    aes(
      F2, F1
    )
  )+
  geom_text(
    aes(
      label = label,
      color = speaker
    ),
    size = 6,
    key_glyph = draw_key_rect,
    fontface = "bold"
  )+
  facet_wrap(
    ~method, scales = "free"
  )+
  scale_x_reverse()+
  scale_y_reverse()+
  theme(
    aspect.ratio = 0.9
  )
```

Broadly similar patterns emerge in the normalized point values as did in the normalized formant tracks (e.g. Speaker A's fronter /uw/ and/ ow/ in all normalizations, and speaker B's slightly backer /iy/ and /ey/ in the delta F and Nearey normalizations).

## Discussion

In this section, we've seen that normalization methods originally defined for formant point measurements can be directly translated into DCT coefficient normalization methods, and that the results on formant tracks are broadly similar to the results on formant points.

To briefly reiterate and refine the process of defining a DCT normalization procedure:

1.  Any non-linear transformations must be applied to the formant tracks *before* the DCT is applied.
2.  Any normalization parameters based on weighted sums (like the mean) or the standard deviation can be estimated using the 0^th^ DCT coefficient.
3.  Location ($l$) parameters can be directly added or subtracted from the 0^th^ DCT coefficient only.
4.  Scale parameters ($s$) need to be multiplied by $\sqrt{2}$, and applied to all DCT coefficients.

# Conclusion

The Discrete Cosine Transform is a formant track smoothing method that has already found usage in phonetic research on vowel intrinsic spectral change.
The use of DCT coefficients comes with a number of concrete benefits to research on vowel formant dynamics.
The large body of research on speaker normalization can be drawn upon to directly normalize DCT coefficients, as approximately equivalent normalization functions can be defined for DCT coefficients based on the definition of the formant normalization functions.

```{=typst}
#set heading(numbering: "A")
#counter(heading).update(0)
```
# Appendices {.appendix .unnumbered}

# The scipy implementation of the DCT {#sec-the-scipy-implementation-of-the-dct .appendix}

The scipy documentation for the DCT describes three ways the DCT can be "normalized", and two ways the DCT can be "orthogonalized" or "non-orthogonalized".
All of these options on the DCT alter the terms to the left of the sum in the DCT formula.
Let's define and simplify these components.

I'll use $S$ to indicate the sum function, which is defined as

$$
S_k(x) = \sum_{j=0}^{N-1} x_j \cos\left(\frac{\pi k(2j + 1)}{2N} \right)
$$

This term is unaltered by any of the different options scipy offers.
Any given DCT implementation can be given as

$$
y_k = ocS_k(x)
$$ Where $o$ is the orthogonalization term and $c$ is the normalization constant.

The orthogonalization term is the easiest to define.

$$
o = \begin{cases}
  1 & \text{if orth = False}\\
  \frac{1}{\sqrt{2}} & \text{if orth = True}
\end{cases}
$$

The scipy documentation provides the mathematical definition for "backward" normalization constant only, but the "forward" normalization can be inferred from its output.

$$
c = \begin{cases}
  2 & \text{if norm = backward} \\
  \frac{1}{N} & \text{if norm = forward}
\end{cases}
$$

As a demonstration by example, we can define a python function for just the sum function, then apply it to the formant track in @fig-fmt_appendix.

```{r}
#| label: fig-fmt_appendix
#| fig-width: 5
#| fig-height: 3
#| fig-align: center
#| out-width: 70%
#| fig-cap: Demonstration formant track

ay_to_plot |> 
  ggplot(
    aes(time, F1)
  )+
    geom_point()
```

```{python}
#| lst-label: lst-listing_S
#| lst-cap: Definition of the DCT sum function
#| echo: true
import numpy as np
from scipy.fft import dct

def S(x, k):
  j = np.arange(len(x))
  N = len(x)
  
  result = sum(
      x * np.cos(
        (np.pi * k * ((2*j)+1))/
        (2*N)
      )
    )
      
  return result
```

```{python}
np.set_printoptions(precision=2)
```

```{python}
f1 = np.array(r.f1)
```

We can get the result of the sum function for the 0^th^ and 1^st^ DCT coefficients to then examine the outcome of the different normalizations.

```{python}
#| lst-label: lst-listing_ys
#| lst-cap: Sum terms of the DCT
#| echo: true
s_0 = S(f1, 0)
s_1 = S(f1, 1)
```

At this point, we can also get the 0^th^ and 1^st^ DCT coefficients from the scipy implementation.

```{python}
#| lst-label: lst-listing_scipy-dct
#| lst-cap: Application of the scipy DCT
#| echo: true

dct_backward = dct(f1, norm = "backward")
dct_forward  = dct(f1, norm = "forward")
```

The normalizing constant for `norm = "backward"` is documented to be 2, so multiplying `s_0` and `s_1` by 2 should be equal to the 0^th^ and 1^st^ coefficients in `dct_backward`.

```{python}
#| lst-label: lst-listing_norm-backward
#| lst-cap: Backward DCT normalization
#| echo: true
np.array([
  2 * np.array([s_0, s_1]),
  dct_backward[0:2]
])
```

If the normalizing constant for `norm="forward"` is $\frac{1}{N}$, dividing `s_0` and `s_1` by the length of the input vector should be equal to the 0^th^ and 1^st^ coefficients in `dct_forward`.

```{python}
#| lst-label: lst-listing_norm-forward
#| lst-cap: forward DCT normalization
#| echo: true
N = len(f1)

np.array([
  np.array([s_0, s_1]) / N,
  dct_forward[0:2]
])
```

Admittedly, it would be more ideal to be able to reference the actual forward normalization constant from the scipy documentation, but it is not provided.

# The DCT Basis {.appendix}

While the formula in @eq-dctii can be used to calculate the DCT coefficients, the formula to calculate the DCT basis functions in @fig-dct-basis is different.
If $B$ is a matrix of the basis functions, the $k^{th}$ basis function will be in its columns.
To get $B$, we apply the DCT with backward normalization to an identity matrix $I$ (that is, a matrix with 1s along the diagonal, and 0s elsewhere).
The orthogonalization term $o$ is included in @eq-dct-basis-fun.

$$
B_{jk} = 2o \sum_{j=0}^{N-1} I_{jk} \cos\left(\frac{\pi k(2j + 1)}{2N} \right)
$$ {#eq-dct-basis-fun}

This can be quickly implemented using the scipy DCT implementation like so:

```{python}
#| lst-label: lst-listing_dct-basis
#| lst-cap: Getting the DCT basis functions
#| echo: true
N = len(f1)
dct_basis = dct(np.eye(N), norm = "backward", orthogonalize = True)
```

# The choice of orthogonalization {#sec-the-choice-of-orthogonalization .appendix}

The choice of "orthogonalizing" the DCT coefficients, that is, dividing the 0^th^ coefficient by $\sqrt{2}$, does introduce some awkwardness into the normalization procedures described here.
Using the orthogonalized DCT was a design decision within fasttrackpy due to its reliance on regression-based DCT coefficients.

As a practical issue, formant tracking sometimes returns missing, or NA values for some, but not all, time points along a formant track.
With missing values, the DCT cannot be directly applied.
However, the DCT coefficients can be approximated by linear regression, using the DCT basis as the "predictors".

```{python}
#| lst-label: lst-listing_dct-reg
#| lst-cap: Comparison of direct vs regression-based DCT
#| echo: true
dct_reg = np.linalg.lstsq(dct_basis, f1, rcond=None)[0]
dct_direct = dct(f1, norm = "forward", orthogonalize = True)

np.array([
  dct_reg[0:3],
  dct_direct[0:3]
])
```

Orthogonalizing the first coefficient was the only option that resulted in the same coefficients for both regression and direct DCT within the scipy implementation.
Without orthogonalizing the first coefficient, the 0^th^ coefficient is not equal between the regression based DCT and direct DCT.

```{python}
#| lst-label: lst-listing_dct-reg
#| lst-cap: |
#|   Comparison of direct vs regression-based 
#|   non-orthogonalized DCT
#| echo: true

dct_nonorth_basis = dct(np.eye(N), norm = "backward", orthogonalize = False)

dct_nonorth_reg = np.linalg.lstsq(dct_nonorth_basis, f1, rcond=None)[0]
dct_nonorth_direct = dct(f1, norm="forward", orthogonalize=False)

np.array([
  dct_nonorth_reg[0:3],
  dct_nonorth_direct[0:3]
])
```

Since the design decision to orthogonalize the DCT coefficients was made within fasttrackpy, which was the tool used to arrive at these DCT coefficients in this paper, this was also the version of the DCT used here.

::: {.content-visible when-format="html"}
# Acknowledgements {.appendix .unnumbered}

```{r}
#| results: asis
cat(str_squish(rmarkdown::metadata$acknowledgements))
```
:::

# References {.appendix .unnumbered}
