librosa的短时傅里叶实现librosa.stft()

目录

1. API参数解析

2. 分帧机制

2.1 C语言实现

2.1.1 测试代码

2.2 librosa的短时傅里叶变换实现

2.2.1 libsora.stft()的分帧机制

2.2.2 测试代码

2.2.3 运行结果分析

2.3 librosa的短时傅里叶逆变换实现

3. 附录：librosa官网

---

1. API参数解析

2. 分帧机制

Libsora的短时傅里叶变换，分帧有点“不常规”。以256个采样点的傅里叶变换为例。

实验数据：为了更直观、简洁、方便地进行试验和理解API的参数以及实现机制，预先从音频文件中读出256个采样点，用数组或者列表进行存储。同样，汉宁窗也可以实现计算出来，并存储，使用时直接查表，可避免运算量，这是算法处理中常见的“以牺牲空间换取时间效率”的做法。

2.1 C语言实现

直接对这256个点进行傅里叶变换，傅里叶变换的点数就是256，阶数就是8。采用复数的FFT计算算法，对于输入：实数部分就是音频数据，虚部全为0。

2.1.1 测试代码

#define FRAME_LEN 256
#define FFT_ORDER ((int)log2(FRAME_LEN))
float dat_i[FRAME_LEN], dat_a[FRAME_LEN];
float dat_r[FRAME_LEN]={
	0.000214,
	0.000122,
	-0.000061,
	0.000275,
	0.000092,
	0.000336,
	0.000153,
	-0.000092,
	0.000031,
	0.000305,
	0.000641,
	0.000458,
	0.000305,
	0.000153,
	0.000397,
	0.000732,
	0.000885,
	0.000580,
	0.000427,
	0.000397,
	0.000641,
	0.001190,
	0.001282,
	0.000977,
	0.000336,
	0.000458,
	0.001007,
	0.001617,
	0.001465,
	0.000824,
	0.000275,
	0.000366,
	0.001007,
	0.001282,
	0.000885,
	0.000122,
	-0.000336,
	-0.000122,
	0.000336,
	0.000275,
	-0.000153,
	-0.000824,
	-0.000916,
	-0.000580,
	-0.000519,
	-0.000763,
	-0.001099,
	-0.001221,
	-0.001007,
	-0.000580,
	-0.000793,
	-0.000885,
	-0.001038,
	-0.000580,
	-0.000275,
	-0.000214,
	-0.000549,
	-0.000610,
	-0.000061,
	0.000702,
	0.001190,
	0.000763,
	0.000305,
	0.000275,
	0.001068,
	0.001923,
	0.001953,
	0.001190,
	0.000488,
	0.000519,
	0.001251,
	0.001740,
	0.001282,
	0.000397,
	-0.000183,
	0.000000,
	0.000214,
	0.000305,
	-0.000092,
	-0.000458,
	-0.000732,
	-0.000610,
	-0.000671,
	-0.000916,
	-0.000946,
	-0.000824,
	-0.000702,
	-0.001068,
	-0.000977,
	-0.001190,
	-0.000671,
	-0.000244,
	-0.000122,
	-0.000580,
	-0.000854,
	-0.000519,
	0.000153,
	0.000641,
	0.000671,
	0.000214,
	-0.000122,
	0.000336,
	0.000732,
	0.000702,
	0.000397,
	0.000122,
	0.000122,
	0.000275,
	0.000000,
	-0.000458,
	-0.000580,
	-0.000275,
	0.000092,
	-0.000244,
	-0.001373,
	-0.001892,
	-0.001465,
	-0.000153,
	0.000214,
	-0.000641,
	-0.002533,
	-0.002808,
	-0.001526,
	0.000549,
	0.000946,
	-0.000824,
	-0.002594,
	-0.002625,
	-0.000427,
	0.001678,
	0.001587,
	-0.000397,
	-0.002075,
	-0.001526,
	0.000549,
	0.002014,
	0.001740,
	-0.000122,
	-0.001678,
	-0.000916,
	0.000610,
	0.001556,
	0.000977,
	-0.000580,
	-0.001251,
	-0.000854,
	-0.000183,
	0.000214,
	-0.000031,
	-0.000671,
	-0.001007,
	-0.001221,
	-0.000916,
	-0.000824,
	-0.000427,
	-0.000549,
	-0.000732,
	-0.000610,
	-0.000854,
	-0.000671,
	-0.000458,
	-0.000183,
	0.000031,
	0.000275,
	0.000427,
	-0.000092,
	-0.000397,
	0.000000,
	0.000793,
	0.001526,
	0.001465,
	0.000122,
	-0.001068,
	-0.000854,
	0.000946,
	0.002472,
	0.002014,
	-0.000488,
	-0.002350,
	-0.002075,
	0.000702,
	0.002686,
	0.001801,
	-0.001526,
	-0.003662,
	-0.002472,
	0.000793,
	0.002747,
	0.001740,
	-0.001251,
	-0.002991,
	-0.001648,
	0.001160,
	0.003143,
	0.002289,
	0.000244,
	-0.000916,
	-0.000305,
	0.001526,
	0.002960,
	0.002808,
	0.002289,
	0.001221,
	0.000092,
	0.000427,
	0.001404,
	0.003021,
	0.003479,
	0.001709,
	-0.001068,
	-0.002502,
	-0.001007,
	0.001831,
	0.002960,
	0.000427,
	-0.003510,
	-0.005157,
	-0.002930,
	0.000580,
	0.001648,
	-0.001038,
	-0.004486,
	-0.005463,
	-0.002991,
	-0.000122,
	0.000885,
	-0.000488,
	-0.002045,
	-0.002289,
	-0.001129,
	0.000214,
	0.001343,
	0.002441,
	0.003448,
	0.002991,
	0.001282,
	0.000153,
	0.002014,
	0.006134,
	0.008820,
	0.006439,
	0.000793,
	-0.002625,
	0.000488,
	0.007111,
	0.010132,
	0.004669,
	-0.004791,
	-0.009033,
	-0.004028,
	0.004089,
	0.006104,
	-0.002228,
	-0.013214,
	-0.016327,
	-0.009583,
};
 
int main(void)
{
	for (int i=0; i
		dat_i[i]=0.0;
		dat_r[i]=dat_r[i]*hanning[i];
	}
 
	FFT(dat_r, dat_i, dat_a, FRAME_LEN, FFT_ORDER);
//	for (int i=0; i
//		printf("[%d] (%f+%fj)\n", i, dat_r[i], dat_i[i]);
//		printf("%d %f\n", i, dat_r[i]);
//	}
#if 1
	FILE* fp_dbg=NULL;
	fp_dbg = fopen("res_fft.txt", "wb");
	for (int i=0; i
		fprintf(fp_dbg, "[%d]:%f\n", i, dat_r[i]);
	}
	fclose(fp_dbg);
	printf("data output done\n");
#endif
	return 0;
}

2.1.2 运行结果分析

注：

（1）实数FFT的运算结果得到复数序列

（2）第0个和第FRAME_LEN/2个是直流分量，其余的结果具有“共轭对称性”，即DFT的共轭对称性：X(m)=(N-m)

（3）注意运算结果的前129个数据，要与下文的libsora.stft()的结果做对比

根据实数的FFT结果具有“共轭对称性”的性质，可以优化FFT计算算法，将大大减小运算量和存储空间，一般来说像Matlab或者Python都有实数FFT计算的API，比如做256个点的FFT，只输出129个点的数据即可(没必要返回冗余数据)，。为方便理解，这里采用的是复数的FFT算法。测试代码中，个人实现的API简单说明如下：

2.2 librosa的短时傅里叶变换实现

用相同的数据源，由于librosa.stft()默认加的是汉宁窗，这里就没必要先进行计算。

2.2.1 libsora.stft()的分帧机制

(1) librosa.stft输出的帧数比正常计算的帧数多一帧，它分帧策略的准则是中心对齐的，即分帧点位于该帧的中心，首尾帧之外都有半个帧长度的padding。而padding部分使用0数据填充。

(2) 输出帧数的计算公式：signal_length/hop_length + 1

例：比如，一个语音端共有256个点。分帧的帧长为256，帧移取一半，也就是128，那么折叠的部分就是256-128=128。

正常的分帧方法可分为两个帧，第二帧(最后一帧)的后半部分用0填充，或者一般来说，最后一帧可以直接舍弃【但是在实时语音处理中，分段处理，不能这么做！】。如下图：

 而librosa的做法是会在数据的头部和尾部增加半个帧长度的“填充”，再进行分帧，如下图：

输入的数据：data=dat1+dat2

第一帧数据：128个0+dat1

第二帧数据：data=dat1+dat2

第三帧数据：dat2+128个0

每次的傅里叶变换都是256个点的输入。

2.2.2 测试代码

文件：libsora_test.py

import librosa
from librosa.core.spectrum import amplitude_to_db
import numpy as np
import soundfile as sf
import matplotlib.pyplot as plt
 
dat_r = np.array([
0.000214,
0.000122,
-0.000061,
0.000275,
0.000092,
0.000336,
0.000153,
-0.000092,
0.000031,
0.000305,
0.000641,
0.000458,
0.000305,
0.000153,
0.000397,
0.000732,
0.000885,
0.000580,
0.000427,
0.000397,
0.000641,
0.001190,
0.001282,
0.000977,
0.000336,
0.000458,
0.001007,
0.001617,
0.001465,
0.000824,
0.000275,
0.000366,
0.001007,
0.001282,
0.000885,
0.000122,
-0.000336,
-0.000122,
0.000336,
0.000275,
-0.000153,
-0.000824,
-0.000916,
-0.000580,
-0.000519,
-0.000763,
-0.001099,
-0.001221,
-0.001007,
-0.000580,
-0.000793,
-0.000885,
-0.001038,
-0.000580,
-0.000275,
-0.000214,
-0.000549,
-0.000610,
-0.000061,
0.000702,
0.001190,
0.000763,
0.000305,
0.000275,
0.001068,
0.001923,
0.001953,
0.001190,
0.000488,
0.000519,
0.001251,
0.001740,
0.001282,
0.000397,
-0.000183,
0.000000,
0.000214,
0.000305,
-0.000092,
-0.000458,
-0.000732,
-0.000610,
-0.000671,
-0.000916,
-0.000946,
-0.000824,
-0.000702,
-0.001068,
-0.000977,
-0.001190,
-0.000671,
-0.000244,
-0.000122,
-0.000580,
-0.000854,
-0.000519,
0.000153,
0.000641,
0.000671,
0.000214,
-0.000122,
0.000336,
0.000732,
0.000702,
0.000397,
0.000122,
0.000122,
0.000275,
0.000000,
-0.000458,
-0.000580,
-0.000275,
0.000092,
-0.000244,
-0.001373,
-0.001892,
-0.001465,
-0.000153,
0.000214,
-0.000641,
-0.002533,
-0.002808,
-0.001526,
0.000549,
0.000946,
-0.000824,
-0.002594,
-0.002625,
-0.000427,
0.001678,
0.001587,
-0.000397,
-0.002075,
-0.001526,
0.000549,
0.002014,
0.001740,
-0.000122,
-0.001678,
-0.000916,
0.000610,
0.001556,
0.000977,
-0.000580,
-0.001251,
-0.000854,
-0.000183,
0.000214,
-0.000031,
-0.000671,
-0.001007,
-0.001221,
-0.000916,
-0.000824,
-0.000427,
-0.000549,
-0.000732,
-0.000610,
-0.000854,
-0.000671,
-0.000458,
-0.000183,
0.000031,
0.000275,
0.000427,
-0.000092,
-0.000397,
0.000000,
0.000793,
0.001526,
0.001465,
0.000122,
-0.001068,
-0.000854,
0.000946,
0.002472,
0.002014,
-0.000488,
-0.002350,
-0.002075,
0.000702,
0.002686,
0.001801,
-0.001526,
-0.003662,
-0.002472,
0.000793,
0.002747,
0.001740,
-0.001251,
-0.002991,
-0.001648,
0.001160,
0.003143,
0.002289,
0.000244,
-0.000916,
-0.000305,
0.001526,
0.002960,
0.002808,
0.002289,
0.001221,
0.000092,
0.000427,
0.001404,
0.003021,
0.003479,
0.001709,
-0.001068,
-0.002502,
-0.001007,
0.001831,
0.002960,
0.000427,
-0.003510,
-0.005157,
-0.002930,
0.000580,
0.001648,
-0.001038,
-0.004486,
-0.005463,
-0.002991,
-0.000122,
0.000885,
-0.000488,
-0.002045,
-0.002289,
-0.001129,
0.000214,
0.001343,
0.002441,
0.003448,
0.002991,
0.001282,
0.000153,
0.002014,
0.006134,
0.008820,
0.006439,
0.000793,
-0.002625,
0.000488,
0.007111,
0.010132,
0.004669,
-0.004791,
-0.009033,
-0.004028,
0.004089,
0.006104,
-0.002228,
-0.013214,
-0.016327,
-0.009583,
])
 
res = librosa.stft(dat_r, n_fft=256, hop_length=128, win_length=256)
print("shape-res:", res.shape)  
print("row-res:", res.shape[0]) #行数
print("col-res:", res.shape[1]) #列数
#np.savetxt('dat_r_fft.txt', res, fmt="%f") #保存结果矩阵(复数矩阵，a+bi的形式)
#np.savetxt('dat_r_fft.txt', np.real(res), fmt="%f") #保存每一帧的实部
#np.savetxt('dat_r_fft.txt', np.imag(res), fmt="%f") #保存每一帧的虚部
 
lis0 = []
lis1 = []
lis2 = []
for i in range(res.shape[0]):
	lis0.append(np.real(res[i][0])) #取第一帧的实数部分(结果矩阵的第一列数据)
	lis1.append(np.real(res[i][1])) #取第二帧的实数部分(结果矩阵的第二列数据)
	lis2.append(np.real(res[i][2])) #取第三帧的实数部分(结果矩阵的第三列数据)
	
np.savetxt('res_fft_0.txt', lis0, fmt="%f")
np.savetxt('res_fft_1.txt', lis1, fmt="%f")
np.savetxt('res_fft_2.txt', lis2, fmt="%f")

2.2.3 运行结果分析

为方便对比数据，将每一帧的运算结果的实数部分保存到txt文本中。

对比第二帧的结果，如下：可发现前129个数据近似相等[自己写的复数FFT算法有计算精度上的差异]。

 同样可以对比验证，第一和第三帧的FFT结果，从而证实librosa的分帧机制。

2.3 librosa的短时傅里叶逆变换实现

对应的API是：librosa.istft

参数一致即可。逆变换会根据正变换的规则，处理重叠部分，还原数据帧。测试代码如下：

对比结果，和输入序列数据是一致的。

3. 附录：librosa官网

Core IO and DSP — librosa 0.9.2 documentation