3.10. Fundamental frequency (F0) — Introduction to Speech Processing

Excerpt

The fundamental frequency of a speech signal, often denoted by F0 or , refers to the approximate frequency of the (quasi-)periodic structure of voiced speech signals. The oscillation originates from the vocal folds, which oscillate in the airflow when appropriately tensed. The fundamental frequency is defined as the average number of oscillations per second and expressed in Hertz. Since the oscillation originates from an organic structure, it is not exactly periodic but contains significant fluctuations. In particular, amount of variation in period length and amplitude are known respectively as jitter and shimmer. Moreover, the F0 is typically not stationary, but changes constantly within a sentence. In fact, the F0 can be used for expressive purposes to signify, for example, emphasis and questions.

---

3.10.1. Introduction

The fundamental frequency of a speech signal, often denoted by F0 or , refers to the approximate frequency of the (quasi-)periodic structure of voiced speech signals. The oscillation originates from the vocal folds, which oscillate in the airflow when appropriately tensed. The fundamental frequency is defined as the average number of oscillations per second and expressed in Hertz. Since the oscillation originates from an organic structure, it is not exactly periodic but contains significant fluctuations. In particular, amount of variation in period length and amplitude are known respectively as jitter and shimmer. Moreover, the F0 is typically not stationary, but changes constantly within a sentence. In fact, the F0 can be used for expressive purposes to signify, for example, emphasis and questions.

Typically fundamental frequencies lie roughly in the range 80 to 450 Hz, where males have lower voices than females and children. The F0 of an individual speaker depends primarily on the length of the vocal folds, which is in turn correlated with overall body size. Cultural and stylistic aspects of speech naturally have also a large impact.

The fundamental frequency is closely related to pitch, which is defined as our perception of fundamental frequency. That is, the F0 describes the actual physical phenomenon, whereas pitch describes how our ears and brains interpret the signal, in terms of periodicity. For example, a voice signal could have an F0 of 100 Hz. If we then apply a high-pass filter to remove all signal components below 450 Hz, then that would remove the actual fundamental frequency. The lowest remaining periodic component would be 500 Hz, which correspond to the fifth harmonic of the original F0. However, a human listener would then typically still perceive a pitch of 100 Hz, even if it does not exist anymore. The brain somehow reconstructs the fundamental from the upper harmonics. This well-known phenomenon is however still not completely understood.

A speech signal with a fundamental frequency of approximately F0=93Hz.

The spectrum of a speech signal with a fundamental frequency of approximately F0=93Hz (original) and a high-pass filtered version of it such that the fundamental frequency has been removed (high-pass filtered).

A speech signal with a fundamental frequency of approximately F0=93Hz

Show code cell source

<span></span><span>import</span> <span>IPython.display</span> <span>as</span> <span>ipd</span>
<span>ipd</span><span>.</span><span>Audio</span><span>(</span><span>'attachments/175515683.wav'</span><span>)</span>

Your browser does not support the audio element.

A high-pass filtered version of it such that the fundamental frequency has been removed

Show code cell source

<span></span><span>import</span> <span>IPython.display</span> <span>as</span> <span>ipd</span>
<span>ipd</span><span>.</span><span>Audio</span><span>(</span><span>'attachments/175515684.wav'</span><span>)</span>

Your browser does not support the audio element.

If is the fundamental frequency, then the length of a single period in seconds is

The speech waveform thus repeats itself after every seconds.

A simple way of modelling the fundamental frequency is to repeat the signal after a delay of seconds. If a signal is sampled with a sampling rate of , then the signal repeats after a delay of samples where

A signal then approximately repeats itself such that

In the Z-domain this can be modelled by an IIR-filter as

where the scalar scales with the accuracy of the period. The Z-transform of the signal can then be written as where is the Z-transform of a single period.

Segment of a speech signal, with the period length , and fundamental frequency .

Spectrum of speech signal with the fundamental frequency and harmonics , as well as the formants F1, F2, F3… Notice how the harmonics form a regular comb-structure.

The magnitude spectrum of , has then a periodic comb-structure. That is, the magnitude spectrum has peaks at , for integer . For a discussion about the fundamental frequency in the cepstral domain, see Cepstrum and MFCC.

Spectrum of fundamental frequency model , showing the characteristic comb-structure with harmonic peaks appearing at integer multiples of .

3.10.2. Fundamental frequency estimation

The fundamental frequency (F0) is central in describing speech signals whereby we need methods for estimating the from speech signals. In speech analysis applications, it can be informative to study the absolute value of the fundamental frequency as such, but more commonly, extraction of the is usually a pre-processing step. For example, in recognition tasks, is often used as a feature for machine learning methods. A voice activity detector could, for instance, set a lower and higher threshold on the , such that sounds with an outside the valid range would be classified as non-speech.

The fundamental frequency is visible in multiple different domains:

• In an acoustic time signal, the is visible as a repetition after every samples.

• In the autocovariance or -correlation, the is visible as a peak at lag as well as its integer multiples .

• In the magnitude, power or log-magnitude spectrum, the is visible as a peak at the frequency , as well as its integer multiples, where  is the sampling frequency.

• In the cepstrum, the is visible as a peak at quefrency  as well as its integer multiples .

Consequently, we can use any of these domains to estimate the fundamental frequency . A typical approach applicable in all domains except the time domain is peak-picking. The fundamental frequency corresponds to a peak in each domain, such that we can determine the by finding the highest peak. For better robustness to spurious peaks and computational efficiency, we naturally limit our search to the range of valid ’s, such as 

The harmonic structure, however, poses a problem for peak-picking. Peaks appear at integer multiples of either or lag T, such that sometimes, by coincidence or due to estimation errors, the harmonic peaks can be higher than the primary peak. Such estimation errors are known as octave errors because the error in corresponds to the musical interval of an octave. A typical post-processing step is, therefore, to check for octave jumps. We can check whether or correspond to a sensible . Alternatively, we can check whether the previous analysis frame had an , an octave, or two octaves off. Depending on the application, we can then fix errors or label problematic zones for later use.

Another problem with peak picking is that peak locations might not align with the samples. For example, in the autocorrelation domain, the actual length of the period could be 100.3 samples. However, the peak in the autocorrelation would then appear at lag 100. One approach would then be to use quadratic interpolation between samples near a peak and use the location of the maximum of the interpolated peak to estimate the peak location. Interpolation makes the estimate less sensitive to noise. For example, background noise could peak at lag 102, so the desired maximum of 100 is lower than the peak at 102. By using data from more samples, as in the interpolation approach, we can, therefore, reduce the likelihood that a single corrupted data point would cause an error.

An alternative approach to peak-picking is to use all distinctive peaks to estimate the jointly. That is, if you find  peaks at frequencies , which approximately correspond to harmonic peaks , then you can approximate Another alternative is to calculate the distance between consecutive peaks and estimate We can also combine these methods as we like.

A further potential improvement is to compare subsequent frames. The pitch in speech is fairly continuous over time, and subsequent frames should, therefore, have similar F0s. This can be used to reduce the danger of octave jumps.

In many other applications, we use discrete Fourier transform (DFT) or cosine transforms (DCT) to resolve frequency components. It would, therefore, be tempting to apply the same approach here. In such a domain, we would also already have a joint estimate which does not rely on single data points. However, note that the spectrum is already the DFT of the time signal, and the cepstrum is the DCT or DFT of the log-spectrum. Additional transforms, therefore, usually do not resolve the ambiguity between harmonics.

A classic algorithm for fundamental frequency estimation is YIN [De Cheveigné and Kawahara, 2002], which has been later improved with a probabilistic version pYIN [Mauch and Dixon, 2014]. YIN is based on peak picking in the autocorrelation domain while emphasizing continuity over time.

3.10.3. Examples

Show code cell source

<span></span><span># Initialization</span>
<span>import</span> <span>matplotlib.pyplot</span> <span>as</span> <span>plt</span>
<span>from</span> <span>scipy.io</span> <span>import</span> <span>wavfile</span>
<span>import</span> <span>scipy</span>
<span>import</span> <span>scipy.fft</span>
<span>import</span> <span>numpy</span> <span>as</span> <span>np</span>
<span>import</span> <span>IPython</span>

<span>filename</span> <span>=</span> <span>'sounds/f0sample.wav'</span>

<span>plt</span><span>.</span><span>rcParams</span><span>.</span><span>update</span><span>({</span>
    <span>"text.usetex"</span><span>:</span> <span>True</span><span>,</span>
    <span>"font.family"</span><span>:</span> <span>"Helvetica"</span>
<span>})</span>

<span># read from storage</span>
<span>fs</span><span>,</span> <span>data</span> <span>=</span> <span>wavfile</span><span>.</span><span>read</span><span>(</span><span>filename</span><span>)</span>
<span>data</span> <span>=</span> <span>data</span><span>[:,</span><span>0</span><span>]</span>
<span>time</span> <span>=</span> <span>np</span><span>.</span><span>arange</span><span>(</span><span>len</span><span>(</span><span>data</span><span>))</span><span>/</span><span>fs</span>

<span>IPython</span><span>.</span><span>display</span><span>.</span><span>display</span><span>(</span><span>IPython</span><span>.</span><span>display</span><span>.</span><span>Audio</span><span>(</span><span>data</span><span>,</span><span>rate</span><span>=</span><span>fs</span><span>))</span>
<span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>6</span><span>,</span><span>3</span><span>))</span>
<span>plt</span><span>.</span><span>plot</span><span>(</span><span>time</span><span>,</span><span>data</span><span>)</span>
<span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
<span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Amplitude'</span><span>)</span>
<span>plt</span><span>.</span><span>yticks</span><span>([])</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

Your browser does not support the audio element.

Excerpts of each tone (columns) with a square (upper row) and a Hann window (lower row).

Show code cell source

<span></span><span>window_length_ms</span> <span>=</span> <span>30</span>

<span>window_length</span> <span>=</span> <span>int</span><span>(</span><span>window_length_ms</span><span>*</span><span>fs</span><span>/</span><span>1000</span><span>)</span>
<span>window_time</span> <span>=</span> <span>np</span><span>.</span><span>arange</span><span>(</span><span>window_length</span><span>)</span>
<span>windowing_function</span> <span>=</span> <span>np</span><span>.</span><span>sin</span><span>(</span><span>np</span><span>.</span><span>pi</span><span>*</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>0.5</span><span>,</span><span>window_length</span><span>,</span><span>1</span><span>)</span><span>/</span><span>window_length</span><span>)</span><span>**</span><span>2</span>


<span>time_sec_list</span> <span>=</span> <span>[</span><span>0.7</span><span>,</span> <span>1.7</span><span>,</span> <span>3</span><span>,</span> <span>4.1</span><span>,</span> <span>5.4</span><span>]</span>
<span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>4</span><span>))</span>

<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>time_ix</span> <span>=</span> <span>int</span><span>(</span><span>time_sec_list</span><span>[</span><span>k</span><span>]</span><span>*</span><span>fs</span><span>)</span> <span>+</span> <span>window_time</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>251</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>time_ix</span><span>/</span><span>fs</span><span>,</span> <span>data</span><span>[</span><span>time_ix</span><span>])</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Amplitude'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>6</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>time_ix</span><span>/</span><span>fs</span><span>,</span> <span>data</span><span>[</span><span>time_ix</span><span>]</span><span>*</span><span>windowing_function</span><span>)</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Amplitude'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>

<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>

The (loosly) periodic structure is visible for each tone as a repeated waveform. Left to right, the pitch increases and the period length becomes shorter such that more periods fit the 30ms analysis window. We also note that for the lowest tone, the periods are nicely visible with the square window (upper row), but less so with the Hann window.

To estimate the fundamental frequency, recall that for a zero-mean signal with a period length , we should have or equivalently, zero expectation For a sample segment of length , the expectation has

Here is a vector containing the samples of to .

To find the fundamental frequency, we can then find that where the above norm is closest to zero.

where constants that do not depend on were omitted, and the minimum was replaced by the maximum by changing the sign. We thus obtained the maximization of the correlation between the signal and its delayed version. Note that this makes this approach equivalent with estimating the maximum of the autocorrelation function.

<span></span><span>f0_range_Hz</span> <span>=</span> <span>np</span><span>.</span><span>array</span><span>([</span><span>80</span><span>,</span> <span>450</span><span>])</span>
<span>f0_range</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>fs</span><span>/</span><span>np</span><span>.</span><span>flipud</span><span>(</span><span>f0_range_Hz</span><span>))</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>print</span><span>(</span><span>f</span><span>"Frequency range </span><span>{</span><span>f0_range_Hz</span><span>[</span><span>0</span><span>]</span><span>}</span><span> to </span><span>{</span><span>f0_range_Hz</span><span>[</span><span>1</span><span>]</span><span>}</span><span> Hz. "</span>
      <span>f</span><span>"Lag range </span><span>{</span><span>(</span><span>1000</span><span>/</span><span>f0_range_Hz</span><span>[</span><span>1</span><span>])</span><span>:</span><span>2.2f</span><span>}</span><span> to </span><span>{</span><span>(</span><span>1000</span><span>/</span><span>f0_range_Hz</span><span>[</span><span>0</span><span>])</span><span>:</span><span>2.2f</span><span>}</span><span> ms or</span><span>\n</span><span>{</span><span>f0_range</span><span>[</span><span>0</span><span>]</span><span>}</span><span> to </span><span>{</span><span>f0_range</span><span>[</span><span>1</span><span>]</span><span>}</span><span> samples with the sampling rate </span><span>{</span><span>fs</span><span>/</span><span>1000</span><span>}</span><span> kHz."</span><span>)</span>

<span>correlation</span> <span>=</span> <span>np</span><span>.</span><span>zeros</span><span>([</span><span>5</span><span>,</span><span>window_length</span><span>-</span><span>1</span><span>])</span>

<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>time_ix</span> <span>=</span> <span>int</span><span>(</span><span>time_sec_list</span><span>[</span><span>k</span><span>]</span><span>*</span><span>fs</span><span>)</span> <span>+</span> <span>window_time</span>
    <span>xwin</span> <span>=</span> <span>data</span><span>[</span><span>time_ix</span><span>]</span><span>*</span><span>windowing_function</span>    
    <span>for</span> <span>L</span> <span>in</span> <span>range</span><span>(</span><span>window_length</span><span>-</span><span>1</span><span>):</span>
        <span>x_k</span> <span>=</span> <span>xwin</span><span>[:</span><span>window_length</span><span>-</span><span>L</span><span>]</span>
        <span>x_kL</span> <span>=</span> <span>xwin</span><span>[</span><span>L</span><span>:</span><span>window_length</span><span>]</span>
        <span>correlation</span><span>[</span><span>k</span><span>,</span><span>L</span><span>]</span> <span>=</span> <span>np</span><span>.</span><span>sum</span><span>(</span><span>x_k</span><span>*</span><span>x_kL</span><span>)</span><span>/</span><span>(</span><span>window_length</span><span>-</span><span>L</span><span>)</span>

<span>max_lag</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span><span>correlation</span><span>[:,</span><span>f0_range</span><span>[</span><span>0</span><span>]:</span><span>f0_range</span><span>[</span><span>1</span><span>]],</span><span>axis</span><span>=</span><span>1</span><span>)</span><span>+</span><span>f0_range</span><span>[</span><span>0</span><span>]</span>
        

<span></span>Frequency range 80 to 450 Hz. Lag range 2.22 to 12.50 ms or
98 to 551 samples with the sampling rate 44.1 kHz.

Show code cell source

<span></span><span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>5</span><span>))</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>1</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>window_length</span><span>-</span><span>1</span><span>)</span><span>/</span><span>fs</span><span>,</span> <span>correlation</span><span>[</span><span>k</span><span>,:])</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Lag $L$ (ms)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Correlation'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>    
    <span>plt</span><span>.</span><span>xticks</span><span>([</span><span>0</span><span>,</span> <span>15</span><span>,</span><span>30</span><span>])</span> 

    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>6</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>window_length</span><span>-</span><span>1</span><span>)</span><span>/</span><span>fs</span><span>,</span> <span>correlation</span><span>[</span><span>k</span><span>,:])</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>]])</span><span>/</span><span>fs</span><span>,</span>
             <span>correlation</span><span>[</span><span>k</span><span>,</span><span>0</span><span>]</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>1.15</span><span>,</span> <span>1.05</span><span>,</span> <span>1.1</span><span>,</span> <span>1.1</span><span>,</span> <span>1.05</span><span>,</span> <span>1.15</span><span>]),</span><span>color</span><span>=</span><span>'red'</span><span>)</span>

    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>max_lag</span><span>[</span><span>k</span><span>]</span><span>/</span><span>fs</span><span>,</span> <span>correlation</span><span>[</span><span>k</span><span>,</span><span>max_lag</span><span>[</span><span>k</span><span>]],</span> <span>'rx'</span><span>)</span>
    <span>plt</span><span>.</span><span>title</span><span>(</span><span>f</span><span>"Peak at lag </span><span>{</span><span>max_lag</span><span>[</span><span>k</span><span>]</span><span>}</span><span>, </span><span>\n</span><span>or </span><span>{</span><span>fs</span><span>/</span><span>max_lag</span><span>[</span><span>k</span><span>]</span><span>:</span><span>2.1f</span><span>}</span><span> Hz"</span><span>)</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Lag $L$ (ms)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Correlation'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>
    <span>plt</span><span>.</span><span>xticks</span><span>([</span><span>0</span><span>,</span> <span>15</span><span>])</span> 
    <span>plt</span><span>.</span><span>xlim</span><span>([</span><span>-</span><span>1</span><span>,</span> <span>16</span><span>])</span>
    

<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

The upper row shows the whole correlation, while the lower row zooms near the valid range of lags (indicated by the red lines). Red crosses show the maximum peak.

With this easy sound example, the maxima are easy to discern. However, with both the lowest (leftmost) and highest (rightmost) pitches, other peaks are very close in magnitude to the true pitch. As this was an easy signal, it is thus not difficult to imagine more challenging sounds where peak picking would land at the incorrect peak.

We can repeat the experiment using the power spectrum (squared absolute spectrum).

<span></span><span>fft_length</span> <span>=</span> <span>window_length</span>
<span>spectrum_length</span> <span>=</span> <span>(</span><span>fft_length</span><span>//</span><span>2</span><span>)</span><span>+</span><span>1</span>
<span>f0_range_fft_indices</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>f0_range_Hz</span><span>*</span><span>fft_length</span><span>/</span><span>fs</span><span>)</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>power_spectrum</span> <span>=</span> <span>np</span><span>.</span><span>zeros</span><span>([</span><span>5</span><span>,</span><span>spectrum_length</span><span>])</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>time_ix</span> <span>=</span> <span>int</span><span>(</span><span>time_sec_list</span><span>[</span><span>k</span><span>]</span><span>*</span><span>fs</span><span>)</span> <span>+</span> <span>window_time</span>
    <span>xwin</span> <span>=</span> <span>data</span><span>[</span><span>time_ix</span><span>]</span><span>*</span><span>windowing_function</span>    
    <span>power_spectrum</span><span>[</span><span>k</span><span>,:]</span> <span>=</span> <span>np</span><span>.</span><span>abs</span><span>(</span><span>scipy</span><span>.</span><span>fft</span><span>.</span><span>rfft</span><span>(</span><span>xwin</span><span>,</span><span>n</span><span>=</span><span>fft_length</span><span>))</span><span>**</span><span>2</span>

<span>max_frequency_index</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span>
    <span>power_spectrum</span><span>[:,</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]:</span><span>f0_range_fft_indices</span><span>[</span><span>1</span><span>]],</span>
    <span>axis</span><span>=</span><span>1</span><span>)</span> <span>+</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]</span>
<span>max_frequency</span> <span>=</span> <span>max_frequency_index</span><span>*</span><span>fs</span><span>/</span><span>fft_length</span>        

Show code cell source

<span></span><span>fft_ix</span> <span>=</span> <span>np</span><span>.</span><span>linspace</span><span>(</span><span>0</span><span>,</span><span>fs</span><span>/</span><span>2</span><span>,</span><span>spectrum_length</span><span>)</span>
<span>plot_range_Hz</span> <span>=</span> <span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span> <span>1200</span><span>)</span>
<span>plot_range</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>plot_range_Hz</span><span>*</span><span>fft_length</span><span>/</span><span>fs</span><span>)</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>5</span><span>))</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>1</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>fft_ix</span><span>/</span><span>1000</span><span>,</span> <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,:]))</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Frequency (kHz)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Power (dB)'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>    

    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>6</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>fft_ix</span><span>[</span><span>plot_range</span><span>],</span> 
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,</span><span>plot_range</span><span>]))</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>array</span><span>([</span><span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>]]),</span>
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>np</span><span>.</span><span>max</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,:]))</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>1.15</span><span>,</span> <span>1.05</span><span>,</span> <span>1.1</span><span>,</span> <span>1.1</span><span>,</span> <span>1.05</span><span>,</span> <span>1.15</span><span>]),</span><span>color</span><span>=</span><span>'red'</span><span>)</span>

    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>max_frequency</span><span>[</span><span>k</span><span>],</span> 
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,</span><span>max_frequency_index</span><span>[</span><span>k</span><span>]]),</span> 
             <span>'rx'</span><span>)</span>
    <span>plt</span><span>.</span><span>title</span><span>(</span><span>f</span><span>"Peak at index </span><span>{</span><span>max_frequency_index</span><span>[</span><span>k</span><span>]</span><span>}</span><span>, </span><span>\n</span><span>or </span><span>{</span><span>max_frequency</span><span>[</span><span>k</span><span>]</span><span>:</span><span>2.1f</span><span>}</span><span> Hz"</span><span>)</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Frequency (Hz)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Power (dB)'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>
    <span>plt</span><span>.</span><span>xlim</span><span>([</span><span>-</span><span>1</span><span>,</span> <span>1000</span><span>])</span>
    

<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

Observe that the y-axis is now expressed in decibels. The upper row has the complete spectrum of each tone, while the lower row zooms in on the frequency range where the fundamental is visible.

We immediately observe that for the lowest (leftmost) tone, the pitch is incorrectly estimated as it picks up the second peak of the harmonic structure. In fact, for the two highest pitches, the second harmonic peak is higher than the first, but it is not incorrectly picked up only because it is outside the defined frequency analysis range. Furthermore, the accuracy of estimated frequencies is low because that the density of frequencies is not sufficient. The accuracy can however be improved by oversampling the spectrum, that is, by extending the analysis window by zeros before the Fourier transform.

In other words, let us extend the vector with zeros to a length equal to the sampling rate (44100).

Show code cell source

<span></span><span>fft_length</span> <span>=</span> <span>fs</span>
<span>spectrum_length</span> <span>=</span> <span>(</span><span>fft_length</span><span>//</span><span>2</span><span>)</span><span>+</span><span>1</span>
<span>f0_range_fft_indices</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>f0_range_Hz</span><span>*</span><span>fft_length</span><span>/</span><span>fs</span><span>)</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>power_spectrum</span> <span>=</span> <span>np</span><span>.</span><span>zeros</span><span>([</span><span>5</span><span>,</span><span>spectrum_length</span><span>])</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>time_ix</span> <span>=</span> <span>int</span><span>(</span><span>time_sec_list</span><span>[</span><span>k</span><span>]</span><span>*</span><span>fs</span><span>)</span> <span>+</span> <span>window_time</span>
    <span>xwin</span> <span>=</span> <span>data</span><span>[</span><span>time_ix</span><span>]</span><span>*</span><span>windowing_function</span>    
    <span>power_spectrum</span><span>[</span><span>k</span><span>,:]</span> <span>=</span> <span>np</span><span>.</span><span>abs</span><span>(</span><span>scipy</span><span>.</span><span>fft</span><span>.</span><span>rfft</span><span>(</span><span>xwin</span><span>,</span><span>n</span><span>=</span><span>fft_length</span><span>))</span><span>**</span><span>2</span>

<span>max_frequency_index</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span>
    <span>power_spectrum</span><span>[:,</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]:</span><span>f0_range_fft_indices</span><span>[</span><span>1</span><span>]],</span>
    <span>axis</span><span>=</span><span>1</span><span>)</span> <span>+</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]</span>
<span>max_frequency</span> <span>=</span> <span>max_frequency_index</span><span>*</span><span>fs</span><span>/</span><span>fft_length</span>        

Show code cell source

<span></span><span>fft_ix</span> <span>=</span> <span>np</span><span>.</span><span>linspace</span><span>(</span><span>0</span><span>,</span><span>fs</span><span>/</span><span>2</span><span>,</span><span>spectrum_length</span><span>)</span>
<span>plot_range_Hz</span> <span>=</span> <span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span> <span>1200</span><span>)</span>
<span>plot_range</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>plot_range_Hz</span><span>*</span><span>fft_length</span><span>/</span><span>fs</span><span>)</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>3</span><span>))</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>2</span><span>,</span><span>5</span><span>,</span><span>1</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>fft_ix</span><span>[</span><span>plot_range</span><span>],</span> 
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,</span><span>plot_range</span><span>]))</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>array</span><span>([</span><span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>],</span> 
                            <span>f0_range_Hz</span><span>[</span><span>1</span><span>]]),</span>
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>np</span><span>.</span><span>max</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,:]))</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>1.15</span><span>,</span> <span>1.05</span><span>,</span> <span>1.1</span><span>,</span> <span>1.1</span><span>,</span> <span>1.05</span><span>,</span> <span>1.15</span><span>]),</span><span>color</span><span>=</span><span>'red'</span><span>)</span>

    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>max_frequency</span><span>[</span><span>k</span><span>],</span> 
             <span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>[</span><span>k</span><span>,</span><span>max_frequency_index</span><span>[</span><span>k</span><span>]]),</span> 
             <span>'rx'</span><span>)</span>
    <span>plt</span><span>.</span><span>title</span><span>(</span><span>f</span><span>"Peak at index </span><span>{</span><span>max_frequency_index</span><span>[</span><span>k</span><span>]</span><span>}</span><span>, </span><span>\n</span><span>or </span><span>{</span><span>max_frequency</span><span>[</span><span>k</span><span>]</span><span>:</span><span>2.1f</span><span>}</span><span> Hz"</span><span>)</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Frequency (Hz)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Power (dB)'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>
    <span>plt</span><span>.</span><span>xlim</span><span>([</span><span>-</span><span>1</span><span>,</span> <span>1000</span><span>])</span>
    

<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

It is now much easier to determine the exact location of each peak and thus the frequency estimates are also more accurate. However, this did not solve the octave jump in the lowest (leftmost) tone, where the second harmonic is still larger than the first, true peak.

The third main feature often used for F0 estimation is the cepstrum.

Show code cell source

<span></span><span>fft_length</span> <span>=</span> <span>window_length</span>
<span>spectrum_length</span> <span>=</span> <span>(</span><span>fft_length</span><span>//</span><span>2</span><span>)</span><span>+</span><span>1</span>
<span>f0_range_fft_indices</span> <span>=</span> <span>np</span><span>.</span><span>round</span><span>(</span><span>f0_range_Hz</span><span>*</span><span>fft_length</span><span>/</span><span>fs</span><span>)</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>)</span>

<span>power_spectrum</span> <span>=</span> <span>np</span><span>.</span><span>zeros</span><span>([</span><span>5</span><span>,</span><span>spectrum_length</span><span>])</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>time_ix</span> <span>=</span> <span>int</span><span>(</span><span>time_sec_list</span><span>[</span><span>k</span><span>]</span><span>*</span><span>fs</span><span>)</span> <span>+</span> <span>window_time</span>
    <span>xwin</span> <span>=</span> <span>data</span><span>[</span><span>time_ix</span><span>]</span><span>*</span><span>windowing_function</span>    
    <span>power_spectrum</span><span>[</span><span>k</span><span>,:]</span> <span>=</span> <span>np</span><span>.</span><span>abs</span><span>(</span><span>scipy</span><span>.</span><span>fft</span><span>.</span><span>rfft</span><span>(</span><span>xwin</span><span>,</span><span>n</span><span>=</span><span>fft_length</span><span>))</span><span>**</span><span>2</span>

<span></span><span>cepstrum</span> <span>=</span> <span>scipy</span><span>.</span><span>fft</span><span>.</span><span>irfft</span><span>(</span><span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>),</span><span>axis</span><span>=</span><span>1</span><span>)</span>

<span>max_cepstral_peak</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span><span>cepstrum</span><span>[:,</span><span>f0_range</span><span>[</span><span>0</span><span>]:</span><span>f0_range</span><span>[</span><span>1</span><span>]],</span><span>axis</span><span>=</span><span>1</span><span>)</span><span>+</span><span>f0_range</span><span>[</span><span>0</span><span>]</span>

Show code cell source

<span></span><span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>3</span><span>))</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>5</span><span>):</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>1</span><span>,</span><span>5</span><span>,</span><span>1</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>window_length</span><span>-</span><span>1</span><span>)</span><span>/</span><span>fs</span><span>,</span> <span>cepstrum</span><span>[</span><span>k</span><span>,:])</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>0</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>],</span> <span>f0_range</span><span>[</span><span>1</span><span>]])</span><span>/</span><span>fs</span><span>,</span>
             <span>cepstrum</span><span>[</span><span>k</span><span>,</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]]</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>1.55</span><span>,</span> <span>1.35</span><span>,</span> <span>1.4</span><span>,</span> <span>1.4</span><span>,</span> <span>1.35</span><span>,</span> <span>1.55</span><span>]),</span><span>color</span><span>=</span><span>'red'</span><span>)</span>

    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>1000</span><span>*</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]</span><span>/</span><span>fs</span><span>,</span> <span>cepstrum</span><span>[</span><span>k</span><span>,</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]],</span> <span>'rx'</span><span>)</span>
    <span>plt</span><span>.</span><span>title</span><span>(</span><span>f</span><span>"Peak at lag </span><span>{</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]</span><span>}</span><span>, </span><span>\n</span><span>or </span><span>{</span><span>fs</span><span>/</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]</span><span>:</span><span>2.1f</span><span>}</span><span> Hz"</span><span>)</span>
    <span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Quefrency $L$ (ms)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span> <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Cepstral magnitude'</span><span>)</span>
    <span>plt</span><span>.</span><span>yticks</span><span>([])</span>    
    <span>plt</span><span>.</span><span>xlim</span><span>([</span><span>0</span><span>,</span> <span>15</span><span>])</span> 
    <span>plt</span><span>.</span><span>ylim</span><span>(</span><span>cepstrum</span><span>[</span><span>k</span><span>,</span><span>max_cepstral_peak</span><span>[</span><span>k</span><span>]]</span><span>*</span><span>np</span><span>.</span><span>array</span><span>([</span><span>-</span><span>1</span><span>,</span><span>1.5</span><span>]))</span>
    

<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

The peaks are much more pronounced here in the cepstrum in the sense that they are sharp and easily distinguishable from the noise. Still, the same danger of octave jumps persists in the lowest and highest tones, as there are other high peaks in the search range (red lines). However, since the peaks are now narrow, we cannot easily interpolate between points to find a more accurate location of the peak.

3.10.4. Pitch contours and -tracking

Pitch is used in speech for expressive purposes. In particular, relatively high pitch (and intensity) segments typically indicate an emphasis. It is therefore important to observe changes in pitch over time. Such pitch-over-time plots are known as pitch contours. The challenge of finding “good” pitch contours, that avoids intermittent octave-jumps, is known as pitch tracking.

Here is a sample sentence that can be used to demonstrate pitch contours and -tracking.

Show code cell source

<span></span><span># Initialization</span>
<span>import</span> <span>matplotlib.pyplot</span> <span>as</span> <span>plt</span>
<span>from</span> <span>scipy.io</span> <span>import</span> <span>wavfile</span>
<span>import</span> <span>scipy</span>
<span>import</span> <span>scipy.fft</span>
<span>import</span> <span>numpy</span> <span>as</span> <span>np</span>
<span>import</span> <span>IPython</span>
<span>import</span> <span>librosa</span>

<span>filename</span> <span>=</span> <span>'sounds/f0speechsample.wav'</span>

<span>plt</span><span>.</span><span>rcParams</span><span>.</span><span>update</span><span>({</span>
    <span>"text.usetex"</span><span>:</span> <span>True</span><span>,</span>
    <span>"font.family"</span><span>:</span> <span>"Helvetica"</span>
<span>})</span>

<span># read from storage</span>
<span>fs</span><span>,</span> <span>data</span> <span>=</span> <span>wavfile</span><span>.</span><span>read</span><span>(</span><span>filename</span><span>)</span>
<span>datalength</span> <span>=</span> <span>len</span><span>(</span><span>data</span><span>)</span>
<span>time</span> <span>=</span> <span>np</span><span>.</span><span>arange</span><span>(</span><span>datalength</span><span>)</span><span>/</span><span>fs</span>

<span>IPython</span><span>.</span><span>display</span><span>.</span><span>display</span><span>(</span><span>IPython</span><span>.</span><span>display</span><span>.</span><span>Audio</span><span>(</span><span>data</span><span>,</span><span>rate</span><span>=</span><span>fs</span><span>))</span>
<span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>6</span><span>,</span><span>3</span><span>))</span>
<span>plt</span><span>.</span><span>plot</span><span>(</span><span>time</span><span>,</span><span>data</span><span>)</span>
<span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
<span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Amplitude'</span><span>)</span>
<span>plt</span><span>.</span><span>yticks</span><span>([])</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

Your browser does not support the audio element.

Show code cell source

<span></span><span>window_step</span> <span>=</span> <span>window_length</span> <span>//</span> <span>2</span>
<span>window_count</span> <span>=</span> <span>(</span><span>datalength</span> <span>-</span> <span>window_length</span><span>)</span><span>//</span><span>window_step</span> <span>-</span> <span>1</span>

<span>f0estimates</span> <span>=</span> <span>np</span><span>.</span><span>zeros</span><span>([</span><span>window_count</span><span>,</span> <span>4</span><span>])</span>

<span>f0pyin</span><span>,</span> <span>voiced_flag</span><span>,</span> <span>voiced_prob</span> <span>=</span> <span>librosa</span><span>.</span><span>pyin</span><span>(</span><span>data</span><span>.</span><span>astype</span><span>(</span><span>float</span><span>),</span> 
                                    <span>sr</span> <span>=</span> <span>fs</span><span>,</span> <span># sampling frequency</span>
                                    <span>fmin</span><span>=</span><span>f0_range_Hz</span><span>[</span><span>0</span><span>],</span> 
                                    <span>fmax</span><span>=</span><span>f0_range_Hz</span><span>[</span><span>1</span><span>],</span> 
                                    <span>frame_length</span><span>=</span><span>window_length</span><span>,</span> 
                                    <span>hop_length</span><span>=</span><span>window_step</span><span>)</span>    

<span>f0estimates</span><span>[:,</span><span>3</span><span>]</span> <span>=</span> <span>f0pyin</span><span>[</span><span>1</span><span>:</span><span>window_count</span><span>+</span><span>1</span><span>]</span>

<span>for</span> <span>window_ix</span> <span>in</span> <span>range</span><span>(</span><span>window_count</span><span>):</span>
    <span>window</span> <span>=</span> <span>data</span><span>[(</span><span>window_ix</span><span>*</span><span>window_step</span><span>):(</span><span>window_ix</span><span>*</span><span>window_step</span><span>+</span><span>window_length</span><span>)]</span> <span>*</span> <span>windowing_function</span>

    <span>correlation</span> <span>=</span> <span>np</span><span>.</span><span>correlate</span><span>(</span><span>window</span><span>,</span><span>window</span><span>,</span><span>mode</span><span>=</span><span>'full'</span><span>)[</span><span>window_length</span><span>-</span><span>1</span><span>:]</span>
    <span>power_spectrum</span> <span>=</span> <span>np</span><span>.</span><span>abs</span><span>(</span><span>scipy</span><span>.</span><span>fft</span><span>.</span><span>rfft</span><span>(</span><span>window</span><span>,</span><span>n</span><span>=</span><span>fft_length</span><span>))</span><span>**</span><span>2</span>
    <span>cepstrum</span> <span>=</span> <span>scipy</span><span>.</span><span>fft</span><span>.</span><span>irfft</span><span>(</span><span>10</span><span>*</span><span>np</span><span>.</span><span>log10</span><span>(</span><span>power_spectrum</span><span>))</span>

    <span>max_corr_index</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span><span>correlation</span><span>[</span><span>f0_range</span><span>[</span><span>0</span><span>]:</span><span>f0_range</span><span>[</span><span>1</span><span>]])</span><span>+</span><span>f0_range</span><span>[</span><span>0</span><span>]</span>
    <span>max_frequency_index</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span>
        <span>power_spectrum</span><span>[</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]:</span><span>f0_range_fft_indices</span><span>[</span><span>1</span><span>]])</span> <span>+</span><span>f0_range_fft_indices</span><span>[</span><span>0</span><span>]</span>
    <span>max_cepstral_index</span> <span>=</span> <span>np</span><span>.</span><span>argmax</span><span>(</span><span>cepstrum</span><span>[</span><span>f0_range</span><span>[</span><span>0</span><span>]:</span><span>f0_range</span><span>[</span><span>1</span><span>]])</span><span>+</span><span>f0_range</span><span>[</span><span>0</span><span>]</span>

    <span>f0estimates</span><span>[</span><span>window_ix</span><span>,</span> <span>0</span><span>]</span> <span>=</span> <span>fs</span><span>/</span><span>max_corr_index</span>
    <span>f0estimates</span><span>[</span><span>window_ix</span><span>,</span> <span>1</span><span>]</span> <span>=</span> <span>max_frequency_index</span><span>*</span><span>fs</span><span>/</span><span>fft_length</span>        
    <span>f0estimates</span><span>[</span><span>window_ix</span><span>,</span> <span>2</span><span>]</span> <span>=</span> <span>fs</span><span>/</span><span>max_cepstral_index</span>

Show code cell source

<span></span><span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>4</span><span>))</span>
<span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span><span>window_count</span><span>)</span><span>*</span><span>fs</span><span>/</span><span>(</span><span>1000</span><span>*</span><span>window_step</span><span>),</span><span>f0estimates</span><span>)</span>
<span>plt</span><span>.</span><span>legend</span><span>([</span><span>'Correlation'</span><span>,</span><span>'Spectrum'</span><span>,</span><span>'Cepstrum'</span><span>,</span><span>'pYIN'</span><span>],</span><span>bbox_to_anchor</span><span>=</span><span>(</span><span>1.05</span><span>,</span> <span>1.0</span><span>))</span>
<span>plt</span><span>.</span><span>title</span><span>(</span><span>'Pitch contour estimates'</span><span>)</span>
<span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
<span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Pitch (Hz)'</span><span>)</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

The figure illustrates four different estimates of the pitch contour; the above presented three methods as well as pYIN, which is an improvement of the classic YIN method [De Cheveigné and Kawahara, 2002, Mauch and Dixon, 2014].

We can observe that near 35s, all four methods are nicely aligned, giving a continuous pitch contour. However, elsewhere, the three homebrew methods are very noisy. The first obvious reason is that the presented methods estimate a pitch in all frames, also when no voiced signal is present. We can use voice activity detection to limit analysis to only speech frames. This still leaves unvoiced frames for which we need some additional detection. We can, for example, check whether the detected peak is “prominent” or large enough in its surroundings such that it really is a peak and not noise. Also, if the detected maximum is at the border of the analysis range, that often indicates an anomaly, e.g. the maximum is likely outside the analysis range. Such heuristics should then be carefully tuned.

The pYIN library function used here (from Librosa) provides such a measure of voicing activity. By filtering the pitch contours to include only voiced frames, we obtain the following.

Show code cell source

<span></span><span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>4</span><span>))</span>
<span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span><span>window_count</span><span>)</span><span>*</span><span>fs</span><span>/</span><span>(</span><span>1000</span><span>*</span><span>window_step</span><span>),</span><span>f0estimates</span><span>*</span><span>voiced_flag</span><span>[</span><span>1</span><span>:</span><span>window_count</span><span>+</span><span>1</span><span>]</span><span>.</span><span>reshape</span><span>([</span><span>window_count</span><span>,</span><span>1</span><span>])</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>))</span>
<span>plt</span><span>.</span><span>legend</span><span>([</span><span>'Correlation'</span><span>,</span><span>'Spectrum'</span><span>,</span><span>'Cepstrum'</span><span>,</span><span>'pYIN'</span><span>],</span><span>bbox_to_anchor</span><span>=</span><span>(</span><span>1.01</span><span>,</span> <span>1.0</span><span>))</span>
<span>plt</span><span>.</span><span>title</span><span>(</span><span>'Pitch contour estimates for voiced frames'</span><span>)</span>
<span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
<span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Pitch (Hz)'</span><span>)</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

This seems vaguely better as a large portion of the noise has been omitted. A better visualization could be to compare the three first estimates with pYIN individually and with discrete points for each frame rather than a line.

Show code cell source

<span></span><span>plt</span><span>.</span><span>figure</span><span>(</span><span>figsize</span><span>=</span><span>(</span><span>8</span><span>,</span><span>6</span><span>))</span>
<span>for</span> <span>k</span> <span>in</span> <span>range</span><span>(</span><span>3</span><span>):</span>
    <span>plt</span><span>.</span><span>subplot</span><span>(</span><span>311</span><span>+</span><span>k</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span><span>window_count</span><span>)</span><span>*</span><span>fs</span><span>/</span><span>(</span><span>1000</span><span>*</span><span>window_step</span><span>),</span><span>f0estimates</span><span>[:,</span><span>3</span><span>]</span><span>*</span><span>voiced_flag</span><span>[</span><span>1</span><span>:</span><span>window_count</span><span>+</span><span>1</span><span>]</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>),</span>
            <span>'r-'</span><span>,</span><span>linewidth</span><span>=</span><span>2</span><span>)</span>
    <span>plt</span><span>.</span><span>plot</span><span>(</span><span>np</span><span>.</span><span>arange</span><span>(</span><span>0</span><span>,</span><span>window_count</span><span>)</span><span>*</span><span>fs</span><span>/</span><span>(</span><span>1000</span><span>*</span><span>window_step</span><span>),</span><span>f0estimates</span><span>[:,</span><span>k</span><span>]</span><span>*</span><span>voiced_flag</span><span>[</span><span>1</span><span>:</span><span>window_count</span><span>+</span><span>1</span><span>]</span><span>.</span><span>astype</span><span>(</span><span>int</span><span>),</span>
            <span>'kx'</span><span>,</span><span>markersize</span><span>=</span><span>1</span><span>)</span>
    <span>plt</span><span>.</span><span>ylabel</span><span>(</span><span>'Pitch (Hz)'</span><span>)</span>
    <span>if</span> <span>k</span><span>==</span><span>0</span><span>:</span>
        <span>plt</span><span>.</span><span>title</span><span>(</span><span>'Autocorrelation estimate of pitch contour'</span><span>)</span>
        <span>plt</span><span>.</span><span>legend</span><span>([</span><span>'pYIN'</span><span>,</span><span>'Correlation'</span><span>],</span><span>bbox_to_anchor</span><span>=</span><span>(</span><span>1.01</span><span>,</span> <span>1.0</span><span>))</span>
    <span>elif</span> <span>k</span><span>==</span><span>1</span><span>:</span>
        <span>plt</span><span>.</span><span>title</span><span>(</span><span>'Spectral estimate of pitch contour'</span><span>)</span>
        <span>plt</span><span>.</span><span>legend</span><span>([</span><span>'pYIN'</span><span>,</span><span>'Spectrum'</span><span>],</span><span>bbox_to_anchor</span><span>=</span><span>(</span><span>1.01</span><span>,</span> <span>1.0</span><span>))</span>
    <span>elif</span> <span>k</span><span>==</span><span>2</span><span>:</span>
        <span>plt</span><span>.</span><span>title</span><span>(</span><span>'Cepstral estimate of pitch contour'</span><span>)</span>
        <span>plt</span><span>.</span><span>legend</span><span>([</span><span>'pYIN'</span><span>,</span><span>'Cepstrum'</span><span>],</span><span>bbox_to_anchor</span><span>=</span><span>(</span><span>1.01</span><span>,</span> <span>1.0</span><span>))</span>
<span>plt</span><span>.</span><span>xlabel</span><span>(</span><span>'Time (s)'</span><span>)</span>
<span>plt</span><span>.</span><span>tight_layout</span><span>()</span>
<span>plt</span><span>.</span><span>show</span><span>()</span>

The (auto)correlation method clearly yields estimates which are closest to the high-quality estimates of pYIN. Estimates extracted from the power spectrum are quantized to discrete levels, which makes them less accurate (with small errors). It also features many large errors, where often the estimate has clearly an octave error, such as at the very first segment where the red line is near 100 Hz, and some points are at the double frequency (an octave higher) near 200 Hz.

The cepstral method was the worst in this comparison as it has a large number of large errors, most of which are not even octave errors but (seemingly) random values. This is a problem since octave errors can be fixed in post-processing by checking whether the result could be multiples of the true value and correcting when multiples are found. Thus, the likelihood that the cepstral estimates could be remedied by post-processing is low.

3.10.5. Discussion

Analysis of fundamental frequency often makes sense as it is a prominent feature of speech and has an expressive function. Classical DSP methods yield good results, though octave errors will always mar estimation. Here, we presented examples only with clean voices without background noises or reverberation. More realistic scenarios obviously bring more problems.

More accurate and robust results can potentially be achieved by deep neural networks such as Crepe or DeepF0 [Kim et al., 2018, Singh et al., 2021]. A challenge with a data-driven approach is, however, to define the ground truth. What is the true pitch? Who defines it and how is it determined? The ground truth is required by most data-driven approaches, yet both manual and automatic pitch estimation are susceptible to errors. The danger is that the model then approximates the method which was used to create the ground-truth, rather than the underlying phenomenon.

3.10.6. References

[DCheveigneK02] (1,2)

Alain De Cheveigné and Hideki Kawahara. Yin, a fundamental frequency estimator for speech and music. The Journal of the Acoustical Society of America, 111(4):1917–1930, 2002. URL: https://doi.org/10.1121/1.1458024.

[KSLB18]

Jong Wook Kim, Justin Salamon, Peter Li, and Juan Pablo Bello. Crepe: a convolutional representation for pitch estimation. In 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 161–165. IEEE, 2018. URL: https://doi.org/10.1109/ICASSP.2018.8461329.

[MD14] (1,2)

Matthias Mauch and Simon Dixon. Pyin: a fundamental frequency estimator using probabilistic threshold distributions. In 2014 ieee international conference on acoustics, speech and signal processing (icassp), 659–663. IEEE, 2014. URL: https://doi.org/10.1109/ICASSP.2014.6853678.

[SWQ21]

Satwinder Singh, Ruili Wang, and Yuanhang Qiu. Deepf0: end-to-end fundamental frequency estimation for music and speech signals. In ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 61–65. IEEE, 2021. URL: https://doi.org/10.1109/ICASSP39728.2021.9414050.