本文实现了基于 librosa 的 LFCC 和 CQCC 特征提取,主要参考 librosa 中 MFCC 特征提取的过程,同时使用 torchaudio 来验证 LFCC 的正确性,使用 matlab 来验证 CQCC 的正确性。

LFCC

原理

LFCC 和 MFCC 的区别就是 fliter bank 的不同,MFCC 用的是 mel freq 的滤波器组,而 LFCC 用的是频率为线性分布的滤波器组,因此只要改变 MFCC 中滤波器组的获得方式保持其他代码不变即可。

实现

基于 librosa 库实现的 LFCC 如下:

import warnings
import numpy as np
import librosa
import scipydef linear(sr, n_fft, n_filters=128, fmin=0.0, fmax=None, dtype=np.float32):if fmax is None:fmax = float(sr) / 2# Initialize the weightsn_filters = int(n_filters)weights = np.zeros((n_filters, int(1 + n_fft // 2)), dtype=dtype)# Center freqs of each FFT binfftfreqs = librosa.fft_frequencies(sr=sr, n_fft=n_fft)# 'Center freqs' of liner bands - uniformly spaced between limits# * Need to validatelinear_f = np.linspace(fmin, fmax, n_filters + 2)fdiff = np.diff(linear_f)ramps = np.subtract.outer(linear_f, fftfreqs)for i in range(n_filters):# lower and upper slopes for all binslower = -ramps[i] / fdiff[i]upper = ramps[i + 2] / fdiff[i + 1]# .. then intersect them with each other and zeroweights[i] = np.maximum(0, np.minimum(lower, upper))# Only check weights if f_mel[0] is positiveif not np.all((linear_f[:-2] == 0) | (weights.max(axis=1) > 0)):# This means we have an empty channel somewherewarnings.warn("Empty filters detected in linear frequency basis. ""Some channels will produce empty responses. ""Try increasing your sampling rate (and fmax) or ""reducing n_filters.",stacklevel=2,)return weightsdef linear_spec(y=None,sr=22050,S=None,n_fft=2048,hop_length=512,win_length=None,window='hann',center=True,pad_mode='constant',power=2.0,**kwargs
):if S is not None:# Infer n_fft from spectrogram shape, but only if it mismatchesif n_fft // 2 + 1 != S.shape[-2]:n_fft = 2 * (S.shape[-2] - 1)else:S = (np.abs(librosa.stft(y,n_fft=n_fft,hop_length=hop_length,win_length=win_length,center=center,window=window,pad_mode=pad_mode,))** power)# Build a linear filterlinear_basis = linear(sr=sr, n_fft=n_fft, **kwargs)return np.einsum("...ft,mf->...mt", S, linear_basis, optimize=True)def expand_to(x, *, ndim, axes):try:axes = tuple(axes)except TypeError:axes = tuple([axes])if len(axes) != x.ndim:raise Exception("Shape mismatch between axes={} and input x.shape={}".format(axes, x.shape))if ndim < x.ndim:raise Exception("Cannot expand x.shape={} to fewer dimensions ndim={}".format(x.shape, ndim))shape = [1] * ndimfor i, axi in enumerate(axes):shape[axi] = x.shape[i]return x.reshape(shape)def lfcc(y=None,sr=22050,S=None,n_lfcc=20,dct_type=2,norm='ortho',lifter=0,**kwargs):if S is None:S = librosa.power_to_db(linear_spec(y=y, sr=sr, **kwargs))M = scipy.fftpack.dct(S, axis=-2, type=dct_type,norm=norm)[..., :n_lfcc, :]if lifter > 0:# shape lifter for broadcastingLI = np.sin(np.pi * np.arange(1, 1 + n_lfcc, dtype=M.dtype) / lifter)LI = expand_to(LI, ndim=S.ndim, axes=-2)M *= 1 + (lifter / 2) * LIreturn Melif lifter == 0:return Melse:raise Exception("LFCC lifter={} must be a non-negative number".format(lifter))

用 torchaudio 现有的函数进行验证:

import torchaudio
import librosa
from torchaudio.transforms import LFCC
path = ''# torchaudio
waveform, sample_rate = torchaudio.load(full_path)
print(waveform.shape)
lfcc_transform = LFCC(sample_rate=sample_rate,n_lfcc=13,n_filter=128,speckwargs={"n_fft": 512, "hop_length": 160, },
)
torch_lfcc = lfcc_transform(waveform)
print(torch_lfcc[0, :, 0])
# librosa
wave, sr = librosa.load(full_path, sr=16000)
print(wave.shape)
librosa_lfcc = lfcc(y=wave, sr=sr, n_lfcc=13, n_fft=512, hop_length=240, n_filters=128, pad_mode='reflect')
print(librosa_lfcc[:, 0])

运行结果为:

两者得到的结果基本一致。

CQCC

原理

STFT 和 CQT 可以看成是时间信号频域分析的两个同级的方法。

在 STFT 中,从长的序列中提取固定长度的片段与窗函数相乘后进行 FFT,然后重复应用滑动窗口得到最终的谱图。

STFT 实际上是一个滤波器组,定义 Q 因子为滤波器的中心频率和滤波器带宽之比: Q = f k δ f Q=\frac{f_k}{\delta f} Q=δffk​​在STFT 中,每个滤波器的带宽是恒定的,当频率从低频到高频时,Q 因子增加。

但是人类的感知系统的 Q 因子在 500Hz-20KHz 之间近似常数,因此,至少从感知角度来看,STFT对于语音信号的时频分析可能并不理想。

于是提出 CQT 变换,第一种算法由Youngberg和Boll于1978年提出的(Youngberg-Boll),而另一种算法则是由 Kashima 和MontReynaud-Kashima在1986年提出(Mont Reynaud)。在这些方法中,倍频程(octaves)是几何分布的,如下图:

CQT 计算

离散时域信号 x ( n ) x(n) x(n) 的 CQT 计算如下:
X C Q ( k , n ) = ∑ j = n − ⌊ N k / 2 ⌋ n + ⌊ N k / 2 ⌋ x ( j ) a k ∗ ( j − n + N k / 2 ) X^{C Q}(k, n)=\sum_{j=n-\left\lfloor N_k / 2\right\rfloor}^{n+\left\lfloor N_k / 2\right\rfloor} x(j) a_k^*\left(j-n+N_k / 2\right) XCQ(k,n)=j=n−⌊Nk​/2⌋∑n+⌊Nk​/2⌋​x(j)ak∗​(j−n+Nk​/2)
其中, k = 1 , 2 , ⋯ , K k = 1,2,\cdots, K k=1,2,⋯,K 为频率索引, a k ∗ ( n ) a_k^*(n) ak∗​(n) 为 a k ( n ) a_k(n) ak​(n) 的复共轭, N k N_k Nk​ 为可变的窗口大小, ⌊ ⋅ ⌋ \left\lfloor \cdot \right\rfloor ⌊⋅⌋ 代表向下取整,基函数 a k ( n ) a_k(n) ak​(n) 定义如下:
a k ( n ) = 1 C ( n N k ) exp ⁡ [ i ( 2 π n f k f s + Φ k ) ] a_k(n)=\frac{1}{C}\left(\frac{n}{N_k}\right) \exp \left[i\left(2 \pi n \frac{f_k}{f_s}+\Phi_k\right)\right] ak​(n)=C1​(Nk​n​)exp[i(2πnfs​fk​​+Φk​)]
其中, f k f_k fk​ 是第 k k k 个频带的中心频率, f s f_s fs​ 为采样率, w ( t ) w(t) w(t) 为窗函数, Φ k \Phi_k Φk​ 为相位偏置,缩放因子 C C C 由下式给出:
C = ∑ l = − ⌊ N k / 2 ⌋ ⌊ N k / 2 ⌋ w ( l + N k / 2 N k ) C=\sum_{l=-\left\lfloor N_k / 2\right\rfloor}^{\left\lfloor N_k / 2\right\rfloor} w\left(\frac{l+N_k / 2}{N_k}\right) C=l=−⌊Nk​/2⌋∑⌊Nk​/2⌋​w(Nk​l+Nk​/2​)
为了满足恒Q变换,中心频率必须满足:
f k = f 1 2 k − 1 B f_k=f_1 2^{\frac{k-1}{B}} fk​=f1​2Bk−1​
其中, f 1 f_1 f1​ 是最低频率带的中心频率, B B B 确定每个倍频程的频带数。

例如,假设 B = 1 B=1 B=1,则 f k = f 1 2 k − 1 f_k = f_1 2^{k-1} fk​=f1​2k−1 ,取 f s = 8000 H z , f 1 = 500 H z f_s=8000Hz, f_1 = 500Hz fs​=8000Hz,f1​=500Hz,则 f 2 = 1000 H z , f 3 = 2000 H z , f 4 = 4000 H z , f 5 = 8000 H z f_2=1000Hz, f_3=2000Hz, f_4=4000Hz, f_5=8000Hz f2​=1000Hz,f3​=2000Hz,f4​=4000Hz,f5​=8000Hz,不能再大了。而在DFT中,这几个频率呈线性变化。

则 Q 因子计算如下:
Q = f k f k + 1 − f k = ( 2 1 / B − 1 ) − 1 Q=\frac{f_k}{f_{k+1}-f_k}=\left(2^{1 / B}-1\right)^{-1} Q=fk+1​−fk​fk​​=(21/B−1)−1
同时窗函数长度 N k N_k Nk​ 满足:
N k = f s f k Q N_k=\frac{f_s}{f_k} Q Nk​=fk​fs​​Q

CQCC 提取

首先 MFCC 计算如下:
MFCC ⁡ ( q ) = ∑ m = 1 M log ⁡ [ M F ( m ) ] cos ⁡ [ q ( m − 1 2 ) π M ] \operatorname{MFCC}(q)=\sum_{m=1}^M \log [M F(m)] \cos \left[\frac{q\left(m-\frac{1}{2}\right) \pi}{M}\right] MFCC(q)=m=1∑M​log[MF(m)]cos[Mq(m−21​)π​]
其中,Mel 谱计算如下:
M F ( m ) = ∑ k = 1 K ∣ X D F T ( k ) ∣ 2 H m ( k ) M F(m)=\sum_{k=1}^K\left|X^{D F T}(k)\right|^2 H_m(k) MF(m)=k=1∑K​ ​XDFT(k) ​2Hm​(k)
其中, k k k 是DFT之后的频率索引系数, H m ( k ) H_m(k) Hm​(k) 是第 m m m 个Mel 尺度下的三角加权函数带通滤波器。这里, M M M 代表滤波器的个数( M F ( m ) M F(m) MF(m) 为 M M M 点的序列), q q q 代表离散余弦变换的点数。

倒谱分析不能直接被用于 CQT,因为 X C Q ( k ) X^{C Q}(k) XCQ(k) 和 DCT 的余弦函数的尺度不同(一个是几何一个是线性)。可以通过将几何空间转换为线性空间来解决这个问题。

由于在 X C Q ( k ) X^{C Q}(k) XCQ(k) 中, k k k 个频带几何分布,信号重构的过程可以看成是前 k k k 个频带(低频部分)进行下采样,后 K − k K-k K−k 个频带(高频部分)进行上采样得到的,将第 k k k 个频带的中心频率 f k f_k fk​ 和 第一个频带的中心频率 f 1 = f m i n f_1=f_{min} f1​=fmin​ 的频率差记为:
Δ f k ↔ 1 = f k − f 1 = f 1 ( 2 k − 1 B − 1 ) \Delta f^{k \leftrightarrow 1}=f_k-f_1=f_1\left(2^{\frac{k-1}{B}}-1\right) Δfk↔1=fk​−f1​=f1​(2Bk−1​−1)
其中, k = 1 , 2 , ⋯ , K k = 1,2,\cdots, K k=1,2,⋯,K 为频率索引,距离 Δ f k ↔ 1 \Delta f^{k \leftrightarrow 1} Δfk↔1 为 k k k 的递增函数。我们以周期 T l T_l Tl​ 进行重采样,这等效于确定一个关于 k l k_l kl​ 的函数使得:
T l = Δ f k l ↔ 1 T_l=\Delta f^{k_l \leftrightarrow 1} Tl​=Δfkl​↔1
以下关注第一个倍频程,通过将第一个倍频程以周期 T l T_l Tl​ 进行 d d d 等分,解得 k l k_l kl​ 为:
f 1 d = f 1 ( 2 k l − 1 B − 1 ) → k l = B log ⁡ 2 ( 1 + 1 d ) \frac{f_1}{d}=f_1\left(2^{\frac{k_l-1}{B}}-1\right) \rightarrow k_l=B \log _2\left(1+\frac{1}{d}\right) df1​​=f1​(2Bkl​−1​−1)→kl​=Blog2​(1+d1​)
新的频率计算为:
F l = 1 T l = [ f 1 ( 2 k l − 1 B − 1 ) ] − 1 F_l=\frac{1}{T_l}=\left[f_1\left(2^{\frac{k_l-1}{B}}-1\right)\right]^{-1} Fl​=Tl​1​=[f1​(2Bkl​−1​−1)]−1
因此,在第一个倍频程中有 d d d 个均匀采样,第二个倍频程中有 2 d 2d 2d 个,第 j − 1 j-1 j−1 个倍频程中有 2 j d 2^jd 2jd 个。

信号重构方法采用了多相抗混叠滤波器和样条插值方法,以均匀的采样率 F l F_l Fl​ 对信号进行重新采样。

最后,CQCC 系数计算如下:
C Q C C ( p ) = ∑ l = 1 L log ⁡ ∣ X C Q ( l ) ∣ 2 cos ⁡ [ p ( l − 1 2 ) π L ] C Q C C(p)=\sum_{l=1}^L \log \left|X^{C Q}(l)\right|^2 \cos \left[\frac{p\left(l-\frac{1}{2}\right) \pi}{L}\right] CQCC(p)=l=1∑L​log ​XCQ(l) ​2cos[Lp(l−21​)π​]
其中, p = 0 , 1 , ⋯ , L − 1 p=0,1,\cdots,L-1 p=0,1,⋯,L−1 是重采样之后得频带索引。

我理解是,比如说设置 d = 5 d=5 d=5 ,频率 500 − 1000 H z 500-1000Hz 500−1000Hz 之间采样就可以分成五个 500 , 600 , 700 , 800 , 900 500,600,700,800,900 500,600,700,800,900,而频率 1000 − 2000 H z 1000-2000Hz 1000−2000Hz 之间采样分成 2 d = 10 2d=10 2d=10 个频带,即 1000 , 1100 , 1200 , ⋯ , 1900 1000,1100,1200,\cdots,1900 1000,1100,1200,⋯,1900,同理频率范围在 2000 − 4000 H z 2000-4000Hz 2000−4000Hz 之间的频率可以分成 2 2 ∗ 5 = 20 2^2*5=20 22∗5=20 个频带,以此类推。这样虽然两两中心频率之间是几何间隔,但是采样的数量也是呈几何增长,所以最终重采样得到的频率还是线性的,从而可以进行 DCT 变换。

基于 librosa 的 LFCC 和 CQCC 特征提取相关推荐

  1. 基于Python的卷积神经网络和特征提取

     基于Python的卷积神经网络和特征提取 发表于2015-08-27 21:39| 4577次阅读| 来源blog.christianperone.com/| 13 条评论| 作者Christi ...

  2. python提取人物特征_基于图像人物面部表情识别的特征提取优化方法与流程

    本发明涉及一种基于图像人物面部表情识别的特征提取优化方法,主要利用基于统计特征提取的二维主成分分析法和改进的粒子群算法优化图像矩阵的解,属于图像处理.模式识别和计算机视觉交叉技术应用领域. 背景技术: ...

  3. vlad用python实现_HF-Net(一)基于NetVLAD的global descriptor的特征提取

    0.数据准备及预训练权重 以Aachen Day-Night dataset为例,该数据集目录结构如下:aachen存放在编译HF-Net时设置的DATA_PATH下 基于NetVLAD的global ...

  4. 前景提取 matlab,基于MATLAB的动态前景目标特征提取与运动跟踪

    吴晶鑫 仲梁维 摘 要:动态前景目标识别和提取是计算机视觉领域的重要内容.对动态图像进行前景目标提取与运动跟踪,通过改进高斯混合背景模型,提出一种基于自适应特征加权的前景目标提取算法,目的是对动态画面 ...

  5. HF-Net(一)基于NetVLAD的global descriptor的特征提取

    参考:HF-Net git地址 0.数据准备及预训练权重 以Aachen Day-Night dataset为例,该数据集目录结构如下:aachen存放在编译HF-Net时设置的DATA_PATH下 ...

  6. 基于MATLAB的动态前景目标特征提取与运动跟踪

    基于MATLAB的动态前景目标特征提取与运动跟踪 摘 要:动态前景目标识别和提取是计算机视觉领域的重要内容.对动态图像进行前景目标提取与运动跟踪,通过改进高斯混合背景模型,提出一种基于自适应特征加权的 ...

  7. 基于MATLAB GUI的LMD的特征提取和分析计算

    基于MATLAB  GUI的LMD的特征提取和分析计算 LMD 方法实质上是把任意信号分解成不同尺度的包络信号和纯调频信号, 并定义包络信号和纯调频信号的积为瞬时频率具有物理意义的 PF 分量,不断迭 ...

  8. 基于并联SVM支持向量机训练HOG特征提取的人员目标提取

    1.问题描述: 一:基于视频流的兴趣HOG特征提取 由于,这里所涉及到的场景比较多,角度也比较多,而且你也没有提供合适的样本,所以采用单一的模板库无法提取到能够识别各个情况下的特征数据,针对这个情况, ...

  9. LE-MSFE-DDNet:基于微光增强和多尺度特征提取的缺陷检测网络--论文笔记

    论文的英文名称为:LE–MSFE–DDNet: a defect detection network based on low‑light enhancement and multi‑scale fe ...

最新文章

  1. python 类的封装、继承、重写方法
  2. 软件安装——internal error2503/2502
  3. 【转】Redis安装整理(window平台和Linux平台)
  4. DirectShowPlayerService::doRender: Unresolved error code 0x80040266 (IDispatch error #102)
  5. JavaScript 闭包的详细分享(三种创建方式)(附小实例)
  6. tcp/ip 协议栈Linux内核源码分析八 路由子系统分析三 路由表
  7. [转]那些年我还不懂:IList,ICollection,IEnumerable,IEnumerator,IQueryable
  8. Java 将Excel转为OFD
  9. unity 使用粒子系统 实现一个火焰燃烧效果
  10. 在线富文本html编辑,html编辑器 - 经典富文本网页在线编辑器 - HtmlEditor
  11. 设计基于HTML5的APP登录功能及安全调用接口的方式(原理篇)
  12. 2022餐饮加盟3大核心,让赚钱变得简单
  13. RFSoC应用笔记 - RF数据转换器 -10- RFSoC关键配置之其他功能(一)
  14. 2019北邮网安院机试真题(回忆版)@lantin
  15. 问题: Mac外联硬盘不能更改“-”中的一个或多个项目,因为它们正在使用中
  16. binutils学习笔记
  17. 锐迪科为何选择NASDAQ上市?
  18. 关于算法设计与分析学习报告
  19. 网页不能复制文本的解决办法
  20. Matlab R2016b安装教程

热门文章

  1. 通信协议以及protobuf使用、语法指南一
  2. java jndi使用_JNDI简单入门
  3. 简单方便的去水印(使用python)
  4. egret引入外部字体
  5. 下列计算机设备中属于多媒体输入设备的是,属于多媒体输入设备的是()。
  6. [CF1517G]Starry Night Camping
  7. CSI笔记【7】:Crowd Vetting: Rejecting Adversaries via Collaboration with Application to......论文阅读
  8. 使用autoIt 上传文件
  9. org.apache.catalina.startup.Catalina start之过程分析
  10. 乐视云视频(视频托管)