# LibROSA库提取MFCC特征的过程解析

1 篇文章 1 订阅

LibROSA（本文使用的版本是0.6.3）中的mfcc函数可以用来提取音频的梅尔频率倒谱系数（Mel-Frequency Cepstral Coefficients，MFCCs）特征，MFCC被广泛应用于语音识别。LibROSA的mfcc函数源码如下：

# -- Mel spectrogram and MFCCs -- #
def mfcc(y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm='ortho', **kwargs):
if S is None:
S = power_to_db(melspectrogram(y=y, sr=sr, **kwargs))

return scipy.fftpack.dct(S, axis=0, type=dct_type, norm=norm)[:n_mfcc]

# 源码解析

## 获取梅尔频谱

### 分帧

LibROSA提取MFCC的过程没有预加重步骤，而是直接进行了分帧。分帧函数为librosa.util.frame()，源码如下：

def frame(y, frame_length=2048, hop_length=512):
'''Slice a time series into overlapping frames.

This implementation uses low-level stride manipulation to avoid
redundant copies of the time series data.

Parameters
----------
y : np.ndarray [shape=(n,)]
Time series to frame. Must be one-dimensional and contiguous
in memory.

frame_length : int > 0 [scalar]
Length of the frame in samples

hop_length : int > 0 [scalar]
Number of samples to hop between frames

Returns
-------
y_frames : np.ndarray [shape=(frame_length, N_FRAMES)]
An array of frames sampled from y:
y_frames[i, j] == y[j * hop_length + i]

Raises
------
ParameterError
If y is not contiguous in memory, not an np.ndarray, or
not one-dimensional.  See np.ascontiguous() for details.

If hop_length < 1, frames cannot advance.

If len(y) < frame_length.

'''

if not isinstance(y, np.ndarray):
raise ParameterError('Input must be of type numpy.ndarray, '
'given type(y)={}'.format(type(y)))

if y.ndim != 1:
raise ParameterError('Input must be one-dimensional, '
'given y.ndim={}'.format(y.ndim))

if len(y) < frame_length:
raise ParameterError('Buffer is too short (n={:d})'
' for frame_length={:d}'.format(len(y), frame_length))

if hop_length < 1:
raise ParameterError('Invalid hop_length: {:d}'.format(hop_length))

if not y.flags['C_CONTIGUOUS']:
raise ParameterError('Input buffer must be contiguous.')

# Compute the number of frames that will fit. The end may get truncated.
n_frames = 1 + int((len(y) - frame_length) / hop_length)

# Vertical stride is one sample
# Horizontal stride is hop_length samples
y_frames = as_strided(y, shape=(frame_length, n_frames),
strides=(y.itemsize, hop_length * y.itemsize))
return y_frames

### 加窗

def get_window(window, Nx, fftbins=True):
'''Compute a window function.

This is a wrapper for scipy.signal.get_window that additionally
supports callable or pre-computed windows.

Parameters
----------
window : string, tuple, number, callable, or list-like
The window specification:

- If string, it's the name of the window function (e.g., 'hann')
- If tuple, it's the name of the window function and any parameters
(e.g., ('kaiser', 4.0))
- If numeric, it is treated as the beta parameter of the 'kaiser'
window, as in scipy.signal.get_window.
- If callable, it's a function that accepts one integer argument
(the window length)
- If list-like, it's a pre-computed window of the correct length Nx

Nx : int > 0
The length of the window

fftbins : bool, optional
If True (default), create a periodic window for use with FFT
If False, create a symmetric window for filter design applications.

Returns
-------
get_window : np.ndarray
A window of length Nx and type window

--------
scipy.signal.get_window

Notes
-----
This function caches at level 10.

Raises
------
ParameterError
If window is supplied as a vector of length != n_fft,
or is otherwise mis-specified.
'''
if six.callable(window):
return window(Nx)

elif (isinstance(window, (six.string_types, tuple)) or
np.isscalar(window)):
# TODO: if we add custom window functions in librosa, call them here

return scipy.signal.get_window(window, Nx, fftbins=fftbins)

elif isinstance(window, (np.ndarray, list)):
if len(window) == Nx:
return np.asarray(window)

raise ParameterError('Window size mismatch: '
'{:d} != {:d}'.format(len(window), Nx))
else:
raise ParameterError('Invalid window specification: {}'.format(window))

### 快速傅里叶变换

def stft(y, n_fft=2048, hop_length=None, win_length=None, window='hann',
"""Short-time Fourier transform (STFT)
Returns a complex-valued matrix D such that
np.abs(D[f, t]) is the magnitude of frequency bin f
at frame t
np.angle(D[f, t]) is the phase of frequency bin f
at frame t
Parameters
----------
y : np.ndarray [shape=(n,)], real-valued
the input signal (audio time series)
n_fft : int > 0 [scalar]
FFT window size
hop_length : int > 0 [scalar]
number audio of frames between STFT columns.
If unspecified, defaults win_length / 4.
win_length : int <= n_fft [scalar]
Each frame of audio is windowed by window().
The window will be of length win_length and then padded
with zeros to match n_fft.
If unspecified, defaults to win_length = n_fft.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, or number);
see scipy.signal.get_window
- a window function, such as scipy.signal.hanning
- a vector or array of length n_fft
.. see also:: filters.get_window
center : boolean
- If True, the signal y is padded so that frame
D[:, t] is centered at y[t * hop_length].
- If False, then D[:, t] begins at y[t * hop_length]
dtype : numeric type
Complex numeric type for D.  Default is 64-bit complex.
If center=True, the padding mode to use at the edges of the signal.
By default, STFT uses reflection padding.
Returns
-------
D : np.ndarray [shape=(1 + n_fft/2, t), dtype=dtype]
STFT matrix
--------
istft : Inverse STFT
ifgram : Instantaneous frequency spectrogram
Notes
-----
This function caches at level 20.
"""

# By default, use the entire frame
if win_length is None:
win_length = n_fft

# Set the default hop, if it's not already specified
if hop_length is None:
hop_length = int(win_length // 4)

fft_window = get_window(window, win_length, fftbins=True)

# Pad the window out to n_fft size

# Reshape so that the window can be broadcast
fft_window = fft_window.reshape((-1, 1))

# Check audio is valid
util.valid_audio(y)

# Pad the time series so that frames are centered
if center:

# Window the time series.
y_frames = util.frame(y, frame_length=n_fft, hop_length=hop_length)

# Pre-allocate the STFT matrix
stft_matrix = np.empty((int(1 + n_fft // 2), y_frames.shape[1]),
dtype=dtype,
order='F')

fft = get_fftlib()

# how many columns can we fit within MAX_MEM_BLOCK?
n_columns = int(util.MAX_MEM_BLOCK / (stft_matrix.shape[0] *
stft_matrix.itemsize))

for bl_s in range(0, stft_matrix.shape[1], n_columns):
bl_t = min(bl_s + n_columns, stft_matrix.shape[1])

stft_matrix[:, bl_s:bl_t] = fft.rfft(fft_window *
y_frames[:, bl_s:bl_t],
axis=0)
return stft_matrix

def _spectrogram(y=None, S=None, n_fft=2048, hop_length=512, power=1,
'''Helper function to retrieve a magnitude spectrogram.
This is primarily used in feature extraction functions that can operate on
either audio time-series or spectrogram input.
Parameters
----------
y : None or np.ndarray [ndim=1]
If provided, an audio time series
S : None or np.ndarray
Spectrogram input, optional
n_fft : int > 0
STFT window size
hop_length : int > 0
STFT hop length
power : float > 0
Exponent for the magnitude spectrogram,
e.g., 1 for energy, 2 for power, etc.
win_length : int <= n_fft [scalar]
Each frame of audio is windowed by window().
The window will be of length win_length and then padded
with zeros to match n_fft.
If unspecified, defaults to win_length = n_fft.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, or number);
see scipy.signal.get_window
- a window function, such as scipy.signal.hanning
- a vector or array of length n_fft
.. see also:: filters.get_window
center : boolean
- If True, the signal y is padded so that frame
t is centered at y[t * hop_length].
- If False, then frame t begins at y[t * hop_length]
If center=True, the padding mode to use at the edges of the signal.
By default, STFT uses reflection padding.
Returns
-------
S_out : np.ndarray [dtype=np.float32]
- If S is provided as input, then S_out == S
- Else, S_out = |stft(y, ...)|**power
n_fft : int > 0
- If S is provided, then n_fft is inferred from S
- Else, copied from input
'''

if S is not None:
# Infer n_fft from spectrogram shape
n_fft = 2 * (S.shape[0] - 1)
else:
# Otherwise, compute a magnitude spectrogram from input
S = np.abs(stft(y, n_fft=n_fft, hop_length=hop_length,
win_length=win_length, center=center,

return S, n_fft

_spectrogram函数的power虽然默认值为1，但被上层函数melspectrogram调用时传入的值为2（后面讲到）。

## 梅尔滤波器

def mel(sr, n_fft, n_mels=128, fmin=0.0, fmax=None, htk=False,
norm=1, dtype=np.float32):
"""Create a Filterbank matrix to combine FFT bins into Mel-frequency bins
Parameters
----------
sr        : number > 0 [scalar]
sampling rate of the incoming signal
n_fft     : int > 0 [scalar]
number of FFT components
n_mels    : int > 0 [scalar]
number of Mel bands to generate
fmin      : float >= 0 [scalar]
lowest frequency (in Hz)
fmax      : float >= 0 [scalar]
highest frequency (in Hz).
If None, use fmax = sr / 2.0
htk       : bool [scalar]
use HTK formula instead of Slaney
norm : {None, 1, np.inf} [scalar]
if 1, divide the triangular mel weights by the width of the mel band
(area normalization).  Otherwise, leave all the triangles aiming for
a peak value of 1.0
dtype : np.dtype
The data type of the output basis.
By default, uses 32-bit (single-precision) floating point.
Returns
-------
M         : np.ndarray [shape=(n_mels, 1 + n_fft/2)]
Mel transform matrix
Notes
-----
This function caches at level 10,
"""

if fmax is None:
fmax = float(sr) / 2

if norm is not None and norm != 1 and norm != np.inf:
raise ParameterError('Unsupported norm: {}'.format(repr(norm)))

# Initialize the weights
n_mels = int(n_mels)
weights = np.zeros((n_mels, int(1 + n_fft // 2)), dtype=dtype)

# Center freqs of each FFT bin
fftfreqs = fft_frequencies(sr=sr, n_fft=n_fft)

# 'Center freqs' of mel bands - uniformly spaced between limits
mel_f = mel_frequencies(n_mels + 2, fmin=fmin, fmax=fmax, htk=htk)

fdiff = np.diff(mel_f)
ramps = np.subtract.outer(mel_f, fftfreqs)

for i in range(n_mels):
# lower and upper slopes for all bins
lower = -ramps[i] / fdiff[i]
upper = ramps[i+2] / fdiff[i+1]

# .. then intersect them with each other and zero
weights[i] = np.maximum(0, np.minimum(lower, upper))

if norm == 1:
# Slaney-style mel is scaled to be approx constant energy per channel
enorm = 2.0 / (mel_f[2:n_mels+2] - mel_f[:n_mels])
weights *= enorm[:, np.newaxis]

# Only check weights if f_mel[0] is positive
if not np.all((mel_f[:-2] == 0) | (weights.max(axis=1) > 0)):
# This means we have an empty channel somewhere
warnings.warn('Empty filters detected in mel frequency basis. '
'Some channels will produce empty responses. '
'Try increasing your sampling rate (and fmax) or '
'reducing n_mels.')

return weights

def melspectrogram(y=None, sr=22050, S=None, n_fft=2048, hop_length=512,
power=2.0, **kwargs):
"""Compute a mel-scaled spectrogram.
If a spectrogram input S is provided, then it is mapped directly onto
the mel basis mel_f by mel_f.dot(S).
If a time-series input y, sr is provided, then its magnitude spectrogram
S is first computed, and then mapped onto the mel scale by
mel_f.dot(S**power).  By default, power=2 operates on a power spectrum.
Parameters
----------
y : np.ndarray [shape=(n,)] or None
audio time-series
sr : number > 0 [scalar]
sampling rate of y
S : np.ndarray [shape=(d, t)]
spectrogram
n_fft : int > 0 [scalar]
length of the FFT window
hop_length : int > 0 [scalar]
number of samples between successive frames.
See librosa.core.stft
win_length : int <= n_fft [scalar]
Each frame of audio is windowed by window().
The window will be of length win_length and then padded
with zeros to match n_fft.
If unspecified, defaults to win_length = n_fft.
window : string, tuple, number, function, or np.ndarray [shape=(n_fft,)]
- a window specification (string, tuple, or number);
see scipy.signal.get_window
- a window function, such as scipy.signal.hanning
- a vector or array of length n_fft
.. see also:: filters.get_window
center : boolean
- If True, the signal y is padded so that frame
t is centered at y[t * hop_length].
- If False, then frame t begins at y[t * hop_length]
If center=True, the padding mode to use at the edges of the signal.
By default, STFT uses reflection padding.
power : float > 0 [scalar]
Exponent for the magnitude melspectrogram.
e.g., 1 for energy, 2 for power, etc.
Mel filter bank parameters.
See librosa.filters.mel for details.
Returns
-------
S : np.ndarray [shape=(n_mels, t)]
Mel spectrogram
--------
librosa.filters.mel
Mel filter bank construction
librosa.core.stft
Short-time Fourier Transform

"""

S, n_fft = _spectrogram(y=y, S=S, n_fft=n_fft, hop_length=hop_length, power=power,
win_length=win_length, window=window, center=center,

# Build a Mel filter
mel_basis = filters.mel(sr, n_fft, **kwargs)

return np.dot(mel_basis, S)

## 取对数

LibROSA将这一步骤实现在librosa.core.power_to_db()中：

def power_to_db(S, ref=1.0, amin=1e-10, top_db=80.0):
"""Convert a power spectrogram (amplitude squared) to decibel (dB) units
This computes the scaling 10 * log10(S / ref) in a numerically
stable way.
Parameters
----------
S : np.ndarray
input power
ref : scalar or callable
If scalar, the amplitude abs(S) is scaled relative to ref:
10 * log10(S / ref).
Zeros in the output correspond to positions where S == ref.
If callable, the reference value is computed as ref(S).
amin : float > 0 [scalar]
minimum threshold for abs(S) and ref
top_db : float >= 0 [scalar]
threshold the output at top_db below the peak:
max(10 * log10(S)) - top_db
Returns
-------
S_db : np.ndarray
S_db ~= 10 * log10(S) - 10 * log10(ref)
--------
perceptual_weighting
db_to_power
amplitude_to_db
db_to_amplitude
Notes
-----
This function caches at level 30.

"""

S = np.asarray(S)

if amin <= 0:
raise ParameterError('amin must be strictly positive')

if np.issubdtype(S.dtype, np.complexfloating):
warnings.warn('power_to_db was called on complex input so phase '
'information will be discarded. To suppress this warning, '
magnitude = np.abs(S)
else:
magnitude = S

if six.callable(ref):
# User supplied a function to calculate reference power
ref_value = ref(magnitude)
else:
ref_value = np.abs(ref)

log_spec = 10.0 * np.log10(np.maximum(amin, magnitude))
log_spec -= 10.0 * np.log10(np.maximum(amin, ref_value))

if top_db is not None:
if top_db < 0:
raise ParameterError('top_db must be non-negative')
log_spec = np.maximum(log_spec, log_spec.max() - top_db)

return log_spec

## 离散余弦变换

# -- Mel spectrogram and MFCCs -- #
def mfcc(y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm='ortho', **kwargs):
"""Mel-frequency cepstral coefficients (MFCCs)
Parameters
----------
y : np.ndarray [shape=(n,)] or None
audio time series
sr : number > 0 [scalar]
sampling rate of y
S : np.ndarray [shape=(d, t)] or None
log-power Mel spectrogram
n_mfcc: int > 0 [scalar]
number of MFCCs to return
dct_type : None, or {1, 2, 3}
Discrete cosine transform (DCT) type.
By default, DCT type-2 is used.
norm : None or 'ortho'
If dct_type is 2 or 3, setting norm='ortho' uses an ortho-normal
DCT basis.
Normalization is not supported for dct_type=1.
Arguments to melspectrogram, if operating
on time series input
Returns
-------
M : np.ndarray [shape=(n_mfcc, t)]
MFCC sequence
--------
melspectrogram
scipy.fftpack.dct

"""

if S is None:
S = power_to_db(melspectrogram(y=y, sr=sr, **kwargs))

return scipy.fftpack.dct(S, axis=0, type=dct_type, norm=norm)[:n_mfcc]

# 总结

LibROSA提取MFCC的函数调用链如下图所示：

本文的代码摘自LibROSA的官方GitHub：https://github.com/librosa/librosa

• 40
点赞
• 122
收藏
• 打赏
• 17
评论
11-30 435
04-04 5317
11-09 401
10-04
05-08 2959
06-13 3万+
07-01 3312
09-23 9009
12-12 7611

### “相关推荐”对你有帮助么？

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

___Blue_H

¥2 ¥4 ¥6 ¥10 ¥20

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、C币套餐、付费专栏及课程。