Skip to content

Commit c9fedc1

Browse files
feat(analysis): add mel filterbank and log-mel/MFCC front end (#345)
* feat(analysis): add mel filterbank and MFCC front end Implement MelFilterbank<T,FftSize,NumMelBands> and Mfcc<T,FftSize,NumMelBands,NumCoefficients> for keyword-spotting feature extraction on resource-constrained devices. MelFilterbank constructs HTK triangular filters at runtime from per-band edge frequencies (no full weight matrix stored), supports optional Slaney area normalisation, and exposes static HzToMel / MelToHz helpers. Mfcc composes RealFastFourierTransform → power spectrum → MelFilterbank → log (with ε floor) → direct orthonormal DCT-II (identical to scipy dct(x, type=2, norm='ortho')). Three tests: - mel/Hz HTK round-trip for 200, 1000, 4000 Hz - Slaney per-band discrete area ≈ 1.0 (tolerance 0.15) - Reference frame: unit impulse (flat spectrum), values generated with numpy+scipy using float32 arithmetic Reference values produced by: python3 -c "import numpy as np; from scipy.fft import dct; ..." with FftSize=32, sampleRate=8000, numMelBands=8, numCoefficients=4, fMin=125, fMax=3800, HTK scale, epsilon=1e-12, norm='ortho'. Closes #340 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01N4pngRiBqeb1xiseKXNLG6 * docs: align markdown tables per MegaLinter Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01N4pngRiBqeb1xiseKXNLG6 --------- Co-authored-by: Claude <noreply@anthropic.com>
1 parent e024946 commit c9fedc1

11 files changed

Lines changed: 631 additions & 12 deletions

File tree

‎README.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ Refer to the documentation to quickly integrate and utilize the library's signal
1616

1717
| Category | Description |
1818
|--------------------------------------------------------------------|----------------------------------------------------------------------|
19-
| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Discrete Wavelet Transform (Haar/Daubechies), Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm, Decibels, Hilbert Transform / Analytic Signal |
19+
| [Analysis](doc/analysis/README.md) | FFT, Real-Input FFT (RFFT), Power Spectral Density, DCT, Discrete Wavelet Transform (Haar/Daubechies), Window Functions, Signal Detectors, Convolution & Correlation, Goertzel Algorithm, Decibels, Hilbert Transform / Analytic Signal, Mel Filterbank / MFCC |
2020
| [Control Analysis](doc/control_analysis/README.md) | Frequency Response, Root Locus, Controllability/Observability Matrices & Gramians, Continuous-to-Discrete, Transfer Function ↔ State Space |
2121
| [Controllers](doc/controllers/README.md) | Bang-Bang/Hysteresis, Deadbeat Control, PID, LQR, LQI (Integral/Servo State Feedback), MPC, Saturation, Rate Limiter, Slew-Limited Saturation, Feedforward/2-DOF, Gain-Scheduled Controller, Lead-Lag Compensator, Luenberger Observer |
2222
| [Estimators](doc/estimators/README.md) | Linear Regression, Polynomial Fitting, Total Least Squares, Yule-Walker (offline), Recursive Least Squares, LMS / NLMS Adaptive Filter (online), Consistency Metrics / NEES / NIS |

‎doc/analysis/MelFilterbankMfcc.md‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
# Mel Filterbank and MFCC
2+
3+
## Overview & Motivation
4+
5+
Keyword-spotting and small speech models on microcontrollers need a compact, perceptually motivated
6+
description of a short audio frame. The log-mel spectrum and its decorrelated form, the
7+
Mel-Frequency Cepstral Coefficients (MFCC), compress a frame of $N$ samples into a few tens of
8+
values that follow the coarse spectral envelope on a frequency axis that matches human pitch
9+
perception. Since Davis and Mermelstein (1980) the pipeline has been the standard speech front end,
10+
and every stage has a fixed size, so it fits a no-heap, real-time budget.
11+
12+
The pipeline is: window → real FFT → power spectrum → triangular mel filterbank → floored natural
13+
log (log-mel) → orthonormal DCT-II (MFCC, first $C$ coefficients).
14+
15+
## Mathematical Theory
16+
17+
### Mel scales
18+
19+
Two mel scales are in common use.
20+
21+
**HTK** (O'Shaughnessy):
22+
23+
$$m = 2595 \log_{10}\!\left(1 + \frac{f}{700}\right), \qquad f = 700\left(10^{m/2595} - 1\right)$$
24+
25+
**Slaney** (Auditory Toolbox, librosa default): linear below 1 kHz and logarithmic above, with
26+
$f_{sp} = 200/3$ Hz per mel, a break point $f_b = 1000$ Hz ($m_b = f_b / f_{sp} = 15$) and step
27+
$\lambda = \ln(6.4)/27$:
28+
29+
$$m = \begin{cases} f / f_{sp} & f < f_b \\ m_b + \ln(f / f_b)/\lambda & f \ge f_b \end{cases}
30+
\qquad
31+
f = \begin{cases} f_{sp}\, m & m < m_b \\ f_b\, e^{\lambda (m - m_b)} & m \ge m_b \end{cases}$$
32+
33+
Both pairs are exact inverses, so a Hz → mel → Hz round trip reproduces the input up to rounding.
34+
35+
### Triangular filterbank
36+
37+
For $M$ bands over $[f_\text{min}, f_\text{max}]$, place $M+2$ points equally spaced in mel and map
38+
them back to Hz: $f_0 = f_\text{min} < f_1 < \dots < f_{M+1} = f_\text{max}$. Band $m$ rises from
39+
$f_m$ to its centre $f_{m+1}$ and falls to $f_{m+2}$. FFT bin $k$ sits at $\nu_k = k f_s / N$:
40+
41+
$$H_m(k) = \max\!\left(0,\; \min\!\left(\frac{\nu_k - f_m}{f_{m+1} - f_m},\; \frac{f_{m+2} - \nu_k}{f_{m+2} - f_{m+1}}\right)\right)$$
42+
43+
**Sparse structure.** If $\nu_k$ lies in segment $j$, $[f_j, f_{j+1}]$, only two triangles are
44+
non-zero there. Band $j$ rises with $r_k = (\nu_k - f_j)/(f_{j+1} - f_j)$ and band $j-1$ falls with
45+
$1 - r_k$. Storing one segment index and one weight per bin is therefore enough, so applying the
46+
filterbank costs one pass over the $N/2+1$ bins and needs no $M \times (N/2+1)$ weight matrix.
47+
48+
**Partition of unity.** Without normalisation, for every bin between the first centre $f_1$ and the
49+
last centre $f_M$, the two active weights add up to $r_k + (1 - r_k) = 1$. This holds per bin and
50+
exactly, not just in the continuous limit.
51+
52+
**Slaney (area) normalisation.** Each triangle is scaled by $2/(f_{m+2} - f_m)$, so its continuous
53+
area is 1 and wide high-frequency bands do not dominate. The scale and the normalisation are
54+
independent choices. With the Slaney scale, Slaney normalisation and the same $f_s, N, M,
55+
f_\text{min}, f_\text{max}$, the weights equal librosa's `filters.mel` defaults.
56+
57+
### Log-mel and cepstrum
58+
59+
$$S_m = \sum_k H_m(k)\,|X_k|^2, \qquad \tilde S_m = \ln \max(S_m, \varepsilon)$$
60+
61+
The floor $\varepsilon$ bounds silent or band-limited frames at $\ln \varepsilon$ instead of
62+
$-\infty$. The cepstrum is the orthonormal DCT-II, truncated to the first $C \le M$ terms:
63+
64+
$$c_k = s_k \sum_{m=0}^{M-1} \tilde S_m \cos\!\left(\frac{\pi k (2m+1)}{2M}\right), \qquad s_0 = \sqrt{1/M},\; s_{k>0} = \sqrt{2/M}$$
65+
66+
This is `scipy.fft.dct(·, type=2, norm='ortho')`. The $C \times M$ cosine basis and the $N$ window
67+
coefficients are computed once at construction.
68+
69+
## Complexity Analysis
70+
71+
| Stage | Time | Memory (words) | Notes |
72+
|-------------------|------------------|----------------|------------------------------------------------|
73+
| Construction | $O(N + M + C M)$ | — | Mel points, bin→segment map, window, DCT basis |
74+
| Window + real FFT | $O(N \log N)$ | $O(N)$ | Dominant cost |
75+
| Power spectrum | $O(N)$ | $N/2+1$ | |
76+
| Filterbank | $O(N)$ | $2(N/2+1) + M$ | Two weights per bin plus a gain per band |
77+
| Log | $O(M)$ | $M$ | |
78+
| DCT-II | $O(C M)$ | $C M$ | Direct product with the precomputed basis |
79+
80+
With $N=512$, $M=40$, $C=13$, the tables take about 1.6 k words, where a dense filterbank matrix
81+
would need 10 k.
82+
83+
## Step-by-Step Walkthrough
84+
85+
$f_s = 16$ kHz, $N = 256$ ($\Delta\nu = 62.5$ Hz), $M = 20$, $C = 13$, $[20, 8000]$ Hz,
86+
Slaney scale and normalisation, periodic Hann window, $\varepsilon = 10^{-10}$.
87+
88+
1. $m(20) = 0.3$ and $m(8000) = 45.25$, so the 22 mel points are spaced $2.14$ mel apart. The first
89+
band spans $20$–$305$ Hz with its centre at $162$ Hz.
90+
2. Bin 2 ($125$ Hz) lies in segment 0, giving $r = (125 - 20)/(162 - 20) = 0.74$ for band 0. Bin 3
91+
($187.5$ Hz) lies in segment 1: band 1 rises with $r = 0.17$ and band 0 falls with $0.83$.
92+
3. Band 0's gain is $2/(305 - 20) = 7.0\cdot10^{-3}$. Its librosa row sum is $0.015754$.
93+
4. The test frame (tones at 125, 440 and 1800 Hz plus a 100–7900 Hz chirp) gives
94+
$\tilde S \approx [-1.23, 0.30, 2.16, -0.32, -7.90, \dots, -5.76]$ and
95+
$c \approx [-13.82, 1.49, 2.08, 7.49, \dots]$, which match the float64 librosa + SciPy
96+
reference to within $5\cdot10^{-3}$.
97+
5. An all-zero frame gives $\tilde S_m = \ln 10^{-10} = -23.03$ for every band, so
98+
$c_0 = \sqrt{20}\cdot(-23.03)$ and $c_{k>0} = 0$.
99+
100+
## Pitfalls & Edge Cases
101+
102+
- **Narrow low bands.** When $\Delta\nu$ is wider than the lowest triangles, a band can catch zero
103+
or one bin, and its energy comes from a single bin or is 0, in which case it is floored. Raise $N$
104+
or $f_\text{min}$, or lower $M$.
105+
- **Floor choice.** $\varepsilon$ should sit below the quietest meaningful band energy, but above
106+
the float32 noise of the FFT, about $10^{-7}$ relative to the frame's peak power. Otherwise
107+
quiet bands report rounding noise instead of the floor.
108+
- **Fast-math.** The stages use plain sums and products. Reassociation changes results only at the
109+
rounding level, and the explicit floor (rather than relying on $\ln 0$) keeps the output finite.
110+
- **Scale and normalisation pairing.** Toolkits differ (HTK: HTK scale without normalisation;
111+
librosa: Slaney with Slaney). Match both choices when you compare against a reference.
112+
113+
## Variants & Generalizations
114+
115+
- **Log-mel / FBANK features** stop before the DCT. They are exposed alongside the cepstrum and are
116+
the usual input for CNN spectrogram models.
117+
- **dB scaling** ($10 \log_{10}$) differs from the natural log by a constant factor of
118+
$10/\ln 10$, which the first layer of a learned model absorbs.
119+
- **Liftering and deltas** (temporal derivatives over frames) are post-processing steps on the
120+
cepstrum sequence.
121+
122+
## Applications
123+
124+
Keyword spotting (DS-CNN / "Hello Edge" class models), speaker verification, audio event detection
125+
and voice activity detection. A typical MCU configuration is 25–40 ms frames, a 10–20 ms hop,
126+
$M = 40$ and $C = 10$–$13$.
127+
128+
## Connections to Other Algorithms
129+
130+
The stage consumes the half spectrum of the real-input FFT and the analysis windows. The cepstral
131+
stage is a truncated DCT-II. Power-spectral-density estimation shares the window → FFT → $|X|^2$
132+
front.
133+
134+
## References & Further Reading
135+
136+
- S. Davis, P. Mermelstein, "Comparison of parametric representations for monosyllabic word
137+
recognition in continuously spoken sentences," *IEEE Trans. ASSP* 28(4), 1980.
138+
- M. Slaney, "Auditory Toolbox, Version 2," Interval Research Tech. Report 1998-010, 1998.
139+
- S. Young et al., *The HTK Book*, ch. 5.
140+
- Y. Zhang et al., "Hello Edge: Keyword Spotting on Microcontrollers," arXiv:1711.07128, 2017.
141+
- B. McFee et al., "librosa: Audio and Music Signal Analysis in Python," *SciPy*, 2015.

‎doc/analysis/README.md‎

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -4,17 +4,18 @@ Signal analysis algorithms for frequency-domain decomposition and spectral estim
44

55
## Algorithms
66

7-
| Algorithm | Description |
8-
|---------------------------------------------------------|--------------------------------------------------------------------------------------------------|
9-
| [Fast Fourier Transform](FastFourierTransform.md) | Efficient computation of the Discrete Fourier Transform using the Cooley-Tukey radix-2 algorithm |
10-
| [Real-Input FFT](RealFastFourierTransform.md) | Length-N real FFT via even/odd split into two N/2-point complex DFTs, halving compute and memory |
11-
| [Power Spectral Density](PowerDensitySpectrum.md) | Estimation of signal power distribution across frequencies using Welch's method |
12-
| [Discrete Cosine Transform](DiscreteCosineTransform.md) | Real-valued frequency decomposition via cosine basis functions, computed through FFT |
13-
| [Signal Detectors](SignalDetectors.md) | Peak hold, zero-crossing counter, and RMS envelope detectors for real-time signal monitoring |
14-
| [Decibels](Decibels.md) | `ToDecibels` / `FromDecibels` conversion helpers with zero-floor guard, plus attenuation and ripple utilities |
15-
| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory |
16-
| [Discrete Wavelet Transform](DiscreteWaveletTransform.md) | Multilevel Haar / Daubechies filter bank for O(N) time-frequency decomposition with perfect reconstruction |
17-
| [Hilbert Transform](HilbertTransform.md) | Analytic signal and instantaneous amplitude/phase/frequency via FFT one-sided spectrum or FIR approximation |
7+
| Algorithm | Description |
8+
|-----------------------------------------------------------|----------------------------------------------------------------------------------------------------------------|
9+
| [Fast Fourier Transform](FastFourierTransform.md) | Efficient computation of the Discrete Fourier Transform using the Cooley-Tukey radix-2 algorithm |
10+
| [Real-Input FFT](RealFastFourierTransform.md) | Length-N real FFT via even/odd split into two N/2-point complex DFTs, halving compute and memory |
11+
| [Power Spectral Density](PowerDensitySpectrum.md) | Estimation of signal power distribution across frequencies using Welch's method |
12+
| [Discrete Cosine Transform](DiscreteCosineTransform.md) | Real-valued frequency decomposition via cosine basis functions, computed through FFT |
13+
| [Signal Detectors](SignalDetectors.md) | Peak hold, zero-crossing counter, and RMS envelope detectors for real-time signal monitoring |
14+
| [Decibels](Decibels.md) | `ToDecibels` / `FromDecibels` conversion helpers with zero-floor guard, plus attenuation and ripple utilities |
15+
| [Goertzel Algorithm](GoertzelAlgorithm.md) | Single-bin DFT via a second-order recurrence for O(N) tone detection with O(1) memory |
16+
| [Discrete Wavelet Transform](DiscreteWaveletTransform.md) | Multilevel Haar / Daubechies filter bank for O(N) time-frequency decomposition with perfect reconstruction |
17+
| [Hilbert Transform](HilbertTransform.md) | Analytic signal and instantaneous amplitude/phase/frequency via FFT one-sided spectrum or FIR approximation |
18+
| [Mel Filterbank / MFCC](MelFilterbankMfcc.md) | Sparse triangular mel filterbank (HTK/Slaney scale and normalisation), log-mel and orthonormal DCT-II cepstrum |
1819

1920
## Sub-domains
2021

‎numerical/analysis/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,8 @@ target_sources(numerical.analysis PRIVATE
2020
FastFourierTransformRadix2Impl.hpp
2121
GoertzelAlgorithm.hpp
2222
HilbertTransform.hpp
23+
MelFilterbank.hpp
24+
Mfcc.hpp
2325
PowerDensitySpectrum.hpp
2426
RealFastFourierTransform.hpp
2527
SignalDetectors.hpp
@@ -33,6 +35,8 @@ numerical_add_coverage_sources(numerical.analysis
3335
FastFourierTransformRadix2Impl.cpp
3436
GoertzelAlgorithm.cpp
3537
HilbertTransform.cpp
38+
MelFilterbank.cpp
39+
Mfcc.cpp
3640
PowerDensitySpectrum.cpp
3741
RealFastFourierTransform.cpp
3842
SignalDetectors.cpp
Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
#include "numerical/analysis/MelFilterbank.hpp"
2+
3+
namespace analysis
4+
{
5+
template float HzToMel<float>(float, MelScale);
6+
template float MelToHz<float>(float, MelScale);
7+
template class MelFilterbank<float, 256, 20>;
8+
}

0 commit comments

Comments
 (0)