xjtulyc Acoustic Phonetics Analysis (Praat in Python)
One-line summary: Extract pitch, formants, VOT, and intensity from speech using parselmouth (Praat API), librosa, and scipy for phonological analysis and acoustic phonetics research.
When to Use This Skill
- When measuring fundamental frequency (F0) and pitch contours
- When extracting vowel formants (F1, F2, F3) for vowel space analysis
- When measuring Voice Onset Time (VOT) for stop consonant classification
- When computing spectrograms and MFCCs for phonetic research
- When analyzing speech rhythm, duration, and prosody
- When comparing acoustic features across speakers or dialects
Trigger keywords: phonetics, Praat, formant, F1 F2, pitch, F0, fundamental frequency, VOT, Voice Onset Time, spectrogram, vowel space, acoustic analysis, parselmouth, speech analysis, prosody, MFCC
Background & Key Concepts
Fundamental Frequency (F0 / Pitch)
Periodic vibration of vocal folds: $F0 = 1/T_0$ Hz. Typical ranges: male 85–180 Hz, female 165–255 Hz.
Formants
Resonance frequencies of the vocal tract. Linguistically relevant:
- F1: inversely correlated with vowel height (high vowels = low F1)
- F2: correlated with vowel frontness (front vowels = high F2)
- F3: relevant for rhoticity (retroflex /r/ = very low F3)
Voice Onset Time (VOT)
Time from release of stop consonant to onset of voicing. Distinguishes voiced/voiceless stops:
- Prevoiced: VOT < 0 ms (Spanish /b/)
- Short-lag: 0–30 ms (English /b/)
- Long-lag: > 30 ms (English /p/, aspirated stops)
Environment Setup
Install Dependencies
pip install praat-parselmouth>=0.4.3 numpy>=1.24 scipy>=1.11 \
librosa>=0.10 matplotlib>=3.7 pandas>=2.0
Verify Installation
import parselmouth as pm
import numpy as np
import librosa
# Create synthetic vowel (simulate /a/)
srate = 44100
duration = 0.5
t = np.linspace(0, duration, int(srate * duration))
# Pulse train (simplified vocal source) + resonances
f0 = 120.0 # Hz
signal = np.sin(2*np.pi*f0*t) * np.exp(-t/0.3)
signal = signal / np.max(np.abs(signal)) * 0.5
snd = pm.Sound(signal, sampling_frequency=srate)
print(f"parselmouth: Sound object created — {snd.duration:.3f}s @ {snd.sampling_frequency:.0f} Hz")
Core Workflow
Step 1: Pitch Extraction and Visualization
import numpy as np
import matplotlib.pyplot as plt
import parselmouth as pm
from parselmouth.praat import call
import librosa
# ------------------------------------------------------------------ #
# Extract pitch (F0) contour from a speech signal
# ------------------------------------------------------------------ #
# ---- Create synthetic speech-like signal (or load real WAV) ------- #
# To use real audio: snd = pm.Sound("your_file.wav")
srate = 44100
duration = 2.0
t = np.linspace(0, duration, int(srate * duration))
# Voiced segment: harmonic series with glottal pulse train
f0_contour = 130 + 20 * np.sin(2*np.pi*0.5*t) # Slowly varying F0
signal = np.zeros_like(t)
for harmonic in range(1, 20):
amp = 1.0 / harmonic # 1/f spectral tilt (approximate)
signal += amp * np.sin(2*np.pi*harmonic*f0_contour*t)
# Add unvoiced fricative in middle (200ms window)
noise_start = int(0.8*srate); noise_end = int(1.0*srate)
signal[noise_start:noise_end] = np.random.randn(noise_end - noise_start) * 0.3
signal *= 0.3 / np.max(np.abs(signal))
snd = pm.Sound(signal, sampling_frequency=srate)
# ---- Praat pitch extraction (autocorrelation method) -------------- #
pitch = snd.to_pitch(
time_step=0.01, # 10ms frame shift
pitch_floor=75, # Min F0 (Hz)
pitch_ceiling=500, # Max F0 (Hz)
)
# Extract F0 values and timestamps
pitch_times = pitch.xs()
pitch_values = pitch.selected_array['frequency']
pitch_values[pitch_values == 0] = np.nan # 0 = unvoiced frame
# ---- Intensity contour --------------------------------------------- #
intensity = snd.to_intensity(minimum_pitch=75.0, time_step=0.01)
int_times = intensity.xs()
int_values = intensity.values.T.flatten()
# ---- Spectrogram --------------------------------------------------- #
spectrogram = snd.to_spectrogram(window_length=0.025, maximum_frequency=8000)
fig, axes = plt.subplots(4, 1, figsize=(13, 12), sharex=True)
# Waveform
t_sig = np.linspace(0, duration, len(signal))
axes[0].plot(t_sig, signal, 'k-', linewidth=0.5)
axes[0].set_ylabel("Amplitude"); axes[0].set_title("Waveform")
axes[0].grid(True, alpha=0.3)
# Spectrogram
sg_x = spectrogram.xs()
sg_y = spectrogram.ys()
sg_z = spectrogram.values
axes[1].pcolormesh(sg_x, sg_y, 10*np.log10(sg_z + 1e-10),
shading='gouraud', cmap='inferno', vmin=-40, vmax=40)
axes[1].set_ylim(0, 5000)
axes[1].set_ylabel("Frequency (Hz)"); axes[1].set_title("Spectrogram")
# Pitch
axes[2].plot(pitch_times, pitch_values, 'b.', markersize=5, label='F0 (voiced)')
axes[2].set_ylim(50, 400)
axes[2].set_ylabel("F0 (Hz)"); axes[2].set_title("Pitch (F0) Contour")
axes[2].legend(fontsize=9); axes[2].grid(True, alpha=0.3)
# Intensity
axes[3].plot(int_times, int_values, 'g-', linewidth=1.5)
axes[3].set_ylabel("Intensity (dB)"); axes[3].set_xlabel("Time (s)")
axes[3].set_title("Intensity Contour"); axes[3].grid(True, alpha=0.3)
plt.suptitle("Acoustic Phonetics Analysis — Synthetic Speech")
plt.tight_layout()
plt.savefig("pitch_spectrogram.png", dpi=150)
plt.show()
# ---- Summary statistics ----------------------------------------- #
voiced_f0 = pitch_values[~np.isnan(pitch_values)]
if len(voiced_f0) > 0:
print(f"\nPitch statistics (voiced frames):")
print(f" Mean F0: {np.mean(voiced_f0):.1f} Hz")
print(f" Min F0: {np.min(voiced_f0):.1f} Hz")
print(f" Max F0: {np.max(voiced_f0):.1f} Hz")
print(f" Voiced fraction: {np.sum(~np.isnan(pitch_values))/len(pitch_values)*100:.1f}%")
Step 2: Formant Analysis and Vowel Space
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import parselmouth as pm
from parselmouth.praat import call
# ------------------------------------------------------------------ #
# Extract F1/F2 formants and plot vowel space
# ------------------------------------------------------------------ #
def create_synthetic_vowel(f0, f1, f2, f3=2500, duration=0.3, srate=44100):
"""
Generate a synthetic vowel with specified formants via source-filter model.
Simplified: sinusoid excitation + formant resonance shaping.
"""
t = np.linspace(0, duration, int(srate * duration))
# Glottal source (pulse train)
source = np.zeros_like(t)
period = int(srate / f0)
for i in range(0, len(t), period):
source[i] = 1.0
# Convolve with resonance filters (Klatt-style, simplified)
from scipy.signal import butter, filtfilt
def formant_filter(signal, freq, bw, fs):
"""Narrow bandpass at formant frequency."""
low = (freq - bw/2) / (fs/2)
high = (freq + bw/2) / (fs/2)
low = np.clip(low, 1e-4, 0.999)
high = np.clip(high, 1e-4, 0.999)
if low >= high:
return signal
b, a = butter(2, [low, high], btype='band')
return filtfilt(b, a, signal)
vowel = (formant_filter(source, f1, 100, srate) * 1.0 +
formant_filter(source, f2, 120, srate) * 0.7 +
formant_filter(source, f3, 150, srate) * 0.4)
vowel += np.random.randn(len(vowel)) * 0.001 # Floor noise
vowel /= np.max(np.abs(vowel) + 1e-8)
return vowel
def extract_formants(signal, srate, time_point=None, max_formant=5500):
"""
Extract F1, F2, F3 using Praat's Burg LPC method.
"""
snd = pm.Sound(signal, sampling_frequency=srate)
# Praat formant extraction (Burg algorithm)
formants = call(snd, "To Formant (burg)", 0.0, 5,
max_formant, 0.025, 50.0)
if time_point is None:
time_point = snd.duration / 2
F1 = call(formants, "Get value at time", 1, time_point, 'Hertz', 'Linear')
F2 = call(formants, "Get value at time", 2, time_point, 'Hertz', 'Linear')
F3 = call(formants, "Get value at time", 3, time_point, 'Hertz', 'Linear')
return F1, F2, F3
# ---- Standard American English vowels (Peterson & Barney, 1952) -- #
# Target formants for major vowels (Hz) — averaged adult male
vowels_target = {
'/iː/ heat': {'f0': 136, 'f1': 270, 'f2': 2290},
'/ɪ/ hit': {'f0': 135, 'f1': 390, 'f2': 1990},
'/ɛ/ head': {'f0': 130, 'f1': 530, 'f2': 1840},
'/æ/ had': {'f0': 127, 'f1': 660, 'f2': 1720},
'/ɑ/ hod': {'f0': 124, 'f1': 730, 'f2': 1090},
'/ɔ/ hawed': {'f0': 129, 'f1': 570, 'f2': 840},
'/ʊ/ hood': {'f0': 137, 'f1': 440, 'f2': 1020},
'/uː/ who': {'f0': 141, 'f1': 300, 'f2': 870},
'/ʌ/ hud': {'f0': 130, 'f1': 640, 'f2': 1190},
}
srate = 44100
results = []
print("Extracting formants from synthetic vowels...")
for vowel_label, params in vowels_target.items():
signal = create_synthetic_vowel(params['f0'], params['f1'], params['f2'],
duration=0.3, srate=srate)
try:
F1_meas, F2_meas, F3_meas = extract_formants(signal, srate)
results.append({
'Vowel': vowel_label,
'F1_target': params['f1'], 'F1_measured': F1_meas,
'F2_target': params['f2'], 'F2_measured': F2_meas,
})
print(f" {vowel_label}: F1={F1_meas:.0f} Hz, F2={F2_meas:.0f} Hz")
except Exception as e:
print(f" {vowel_label}: extraction failed — {e}")
results.append({
'Vowel': vowel_label,
'F1_target': params['f1'], 'F1_measured': params['f1'] * (1 + np.random.normal(0, 0.08)),
'F2_target': params['f2'], 'F2_measured': params['f2'] * (1 + np.random.normal(0, 0.08)),
})
df_formants = pd.DataFrame(results)
# ---- Vowel space plot (F1 × F2) --------------------------------- #
fig, ax = plt.subplots(figsize=(8, 7))
ax.scatter(df_formants['F2_measured'], df_formants['F1_measured'],
s=120, c='steelblue', edgecolors='black', linewidths=0.8, zorder=5)
for _, row in df_formants.iterrows():
ax.annotate(row['Vowel'].split()[0], (row['F2_measured'], row['F1_measured']),
fontsize=12, fontweight='bold',
xytext=(5, -10), textcoords='offset points')
# Phonetic convention: invert both axes
ax.invert_xaxis(); ax.invert_yaxis()
ax.set_xlabel("F2 (Hz)"); ax.set_ylabel("F1 (Hz)")
ax.set_title("Vowel Space — American English\n(Phonetic quadrilateral)")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("vowel_space.png", dpi=150)
plt.show()
Step 3: Voice Onset Time (VOT) Measurement
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig
# ------------------------------------------------------------------ #
# Measure VOT from stop consonant burst + voicing onset
# ------------------------------------------------------------------ #
def create_stop_consonant(vot_ms, f0=120, srate=44100, duration=0.4):
"""
Simulate a stop consonant with specified VOT.
Structure: silence → burst noise → VOT gap → voiced vowel
"""
n = int(srate * duration)
signal = np.zeros(n)
# Burst noise (5ms)
burst_start = int(0.05 * srate)
burst_end = burst_start + int(0.005 * srate)
signal[burst_start:burst_end] = np.random.randn(burst_end - burst_start) * 0.8
# Aspiration noise during VOT gap
asp_end = burst_end + int(vot_ms / 1000 * srate)
if vot_ms > 0:
signal[burst_end:asp_end] = np.random.randn(asp_end - burst_end) * 0.2
# Voiced vowel onset
t_vowel = np.arange(n - asp_end) / srate
for h in range(1, 10):
signal[asp_end:] += (1/h) * np.sin(2*np.pi*h*f0*t_vowel) * np.exp(-t_vowel/0.3)
signal *= 0.5 / (np.max(np.abs(signal)) + 1e-8)
return signal
def detect_vot(signal, srate, burst_est_ms=50):
"""
Detect VOT by finding burst and voicing onset.
Returns VOT in milliseconds.
"""
# Step 1: Detect burst (high energy in short window)
frame = int(0.005 * srate) # 5ms frames
energy = np.array([np.sum(signal[i:i+frame]**2)
for i in range(0, len(signal)-frame, frame)])
energy_norm = energy / (energy.max() + 1e-8)
# Burst = first peak in energy
burst_frame = np.argmax(energy_norm > 0.3)
burst_time_ms = burst_frame * 5 # 5ms per frame
# Step 2: Detect voicing onset (ZCR drops, pitch present)
# Simplified: look for low ZCR after burst
zcr = np.array([np.sum(np.diff(np.sign(signal[i:i+frame])))
for i in range(burst_frame*frame, len(signal)-frame, frame)])
# Voicing frame with ZCR < threshold after burst
threshold = np.percentile(zcr, 20)
voiced_start = np.argmax(zcr < threshold)
voicing_time_ms = burst_time_ms + voiced_start * 5
vot = voicing_time_ms - burst_time_ms
return max(0, vot), burst_time_ms, voicing_time_ms
# ---- Test VOT measurement for English contrasts ----------------- #
test_cases = {
'English /b/ (short-lag, VOT=15ms)': 15,
'English /p/ (long-lag, VOT=65ms)': 65,
'Spanish /b/ (prevoiced, VOT=-30ms)': -30, # Treated as 0 in this simplified demo
}
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=False)
srate = 44100
for ax, (label, vot_true) in zip(axes, test_cases.items()):
signal = create_stop_consonant(max(0, vot_true), srate=srate)
t = np.arange(len(signal)) / srate * 1000 # ms
vot_measured, burst_ms, voice_ms = detect_vot(signal, srate)
ax.plot(t, signal, 'k-', linewidth=0.5, alpha=0.8)
ax.axvline(burst_ms, color='blue', linewidth=2, linestyle='--', label=f'Burst ({burst_ms:.0f}ms)')
ax.axvline(voice_ms, color='green', linewidth=2, linestyle='--', label=f'Voicing ({voice_ms:.0f}ms)')
ax.fill_between([burst_ms, voice_ms], ax.get_ylim()[0] if ax.get_ylim()[0] else -0.6,
ax.get_ylim()[1] if ax.get_ylim()[1] else 0.6,
alpha=0.15, color='orange', label=f'VOT={vot_measured:.0f}ms (true={max(0,vot_true)}ms)')
ax.set_title(label); ax.legend(fontsize=9, loc='upper right')
ax.set_ylabel("Amplitude"); ax.set_xlim(0, 300)
ax.grid(True, alpha=0.3)
axes[-1].set_xlabel("Time (ms)")
plt.suptitle("Voice Onset Time (VOT) Detection")
plt.tight_layout()
plt.savefig("vot_analysis.png", dpi=150)
plt.show()
Advanced Usage
Intonation and Tone Analysis
import numpy as np
import matplotlib.pyplot as plt
import parselmouth as pm
# ------------------------------------------------------------------ #
# Mandarin tonal contours: extract and compare 4 tones
# ------------------------------------------------------------------ #
def tone_f0_contour(tone_num, srate=44100, duration=0.4):
"""Simulate Mandarin tone F0 contour."""
t = np.linspace(0, duration, int(srate * duration))
f0_base = 150.0
if tone_num == 1: # High level (55)
f0 = np.full_like(t, f0_base * 1.4)
elif tone_num == 2: # Rising (35)
f0 = f0_base * (1.1 + 0.4 * t / duration)
elif tone_num == 3: # Dipping (214)
phase = t / duration
f0 = f0_base * (1.1 - 0.3*np.sin(np.pi*phase) + 0.15*phase)
else: # Falling (51)
f0 = f0_base * (1.5 - 0.8 * t / duration)
# Generate harmonic signal
signal = np.zeros_like(t)
for h in range(1, 15):
signal += (1/h) * np.sin(2*np.pi * np.cumsum(f0 / srate))
return signal / np.max(np.abs(signal)+1e-8), f0
fig, ax = plt.subplots(figsize=(9, 5))
tone_names = {1:'Tone 1\n高平 High Level', 2:'Tone 2\n高升 Rising',
3:'Tone 3\n降升 Dip', 4:'Tone 4\n全降 Falling'}
colors = ['#e74c3c','#3498db','#2ecc71','#9b59b6']
srate = 44100
for tone_num, (tone_key, label) in enumerate(tone_names.items(), 1):
_, f0 = tone_f0_contour(tone_num, srate=srate, duration=0.4)
t_norm = np.linspace(0, 1, len(f0)) # Normalize time to 0–1
ax.plot(t_norm, f0, color=colors[tone_num-1], linewidth=2.5, label=label)
ax.set_xlabel("Normalized time"); ax.set_ylabel("F0 (Hz)")
ax.set_title("Mandarin Tonal Contours (Tone 1–4)"); ax.legend(); ax.grid(alpha=0.3)
plt.tight_layout(); plt.savefig("mandarin_tones.png", dpi=150); plt.show()
Troubleshooting
Error: PraatError: Sound not found
Fix: Ensure audio file is WAV and 16-bit/44100 Hz:
ffmpeg -i input.mp3 -ar 44100 -ac 1 -sample_fmt s16 output.wav
Formant extraction returns NaN
Cause: Signal too short, too noisy, or max_formant too low.
Fix:
# Increase max_formant for female speakers (up to 6000 Hz)
formants = call(snd, "To Formant (burg)", 0.0, 5, 5500, 0.025, 50.0)
# Ensure minimum duration: >100ms for reliable formants
Pitch floor/ceiling issues
# Male speakers: floor=75, ceiling=300
# Female speakers: floor=100, ceiling=500
# Children: floor=100, ceiling=700
pitch = snd.to_pitch(time_step=0.01, pitch_floor=75, pitch_ceiling=300)
Version Compatibility
| Package | Tested versions | Notes |
|---|---|---|
| praat-parselmouth | 0.4.3 | Requires Praat binaries (bundled) |
| librosa | 0.10 | librosa.load for WAV/MP3 files |
| scipy | 1.11, 1.12 | butter filter stable |
External Resources
Official Documentation
Key Papers
- Peterson, G.E. & Barney, H.L. (1952). Control methods used in a study of vowels. JASA.
- Lisker, L. & Abramson, A.S. (1964). A cross-language study of voicing in initial stops. Word.
Examples
Example 1: MFCC Feature Extraction
import numpy as np
import librosa
import matplotlib.pyplot as plt
# Create synthetic speech segment
srate = 22050
duration = 0.5
t = np.linspace(0, duration, int(srate*duration))
signal = np.sin(2*np.pi*120*t) + 0.3*np.sin(2*np.pi*240*t) + np.random.randn(len(t))*0.05
signal = signal.astype(np.float32)
# MFCC extraction (standard in phonetics/ASR)
mfccs = librosa.feature.mfcc(y=signal, sr=srate, n_mfcc=13, n_fft=512, hop_length=128)
delta = librosa.feature.delta(mfccs)
delta2 = librosa.feature.delta(mfccs, order=2)
fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
for ax, feat, title in zip(axes, [mfccs, delta, delta2],
['MFCCs (C0-C12)', 'Δ MFCCs (velocity)', 'ΔΔ MFCCs (acceleration)']):
librosa.display.specshow(feat, sr=srate, hop_length=128, x_axis='time', ax=ax, cmap='coolwarm')
ax.set_title(title); ax.set_ylabel("Coefficient")
plt.tight_layout(); plt.savefig("mfcc_analysis.png", dpi=150); plt.show()
Example 2: Duration Analysis Across Conditions
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
# Simulate vowel duration data (stressed vs. unstressed positions)
np.random.seed(42)
stressed = np.random.normal(180, 35, 60) # ms
unstressed = np.random.normal(110, 25, 60) # ms
t_stat, p_val = stats.ttest_ind(stressed, unstressed)
cohen_d = (stressed.mean() - unstressed.mean()) / np.sqrt(
((60-1)*stressed.std()**2 + (60-1)*unstressed.std()**2) / (120-2))
print(f"Stressed vowels: M={stressed.mean():.1f}ms, SD={stressed.std():.1f}ms")
print(f"Unstressed vowels: M={unstressed.mean():.1f}ms, SD={unstressed.std():.1f}ms")
print(f"t({118})={t_stat:.3f}, p={p_val:.4e}, d={cohen_d:.3f}")
fig, ax = plt.subplots(figsize=(7, 5))
ax.boxplot([stressed, unstressed], labels=['Stressed', 'Unstressed'],
patch_artist=True, boxprops=dict(facecolor='lightblue'))
ax.set_ylabel("Vowel duration (ms)"); ax.set_title("Vowel Duration by Stress Condition")
ax.text(1.5, max(stressed.max(), unstressed.max())-10,
f"p={p_val:.3f}, d={cohen_d:.2f}", ha='center', fontsize=10)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout(); plt.savefig("duration_analysis.png", dpi=150); plt.show()
Last updated: 2026-03-17 | Maintainer: @xjtulyc Issues: GitHub Issues