刚开始学习MFCC，从网上找了两种方法，求MFCC，试用了下，发现结果完全不同，请高手帮忙解释，或能给出正确结果：
代码如下 ：
import numpy as np
from scipy import signal
from scipy.fftpack import dct
import pylab as plt
import librosa
def enframe(wave_data, nw, inc, winfunc):
'''将音频信号分帧。
参数含义：
wave_data:原始音频型号
nw:每一帧的长度(这里指采样点的长度，即采样频率乘以时间间隔)
inc:相邻帧的间隔
'''
wlen=len(wave_data) #信号总长度
if wlen<=nw:
nf=1
else:
nf=int(np.ceil((1.0*wlen-nw+inc)/inc))
pad_length=int((nf-1)*inc+nw) #所有帧加起来总的铺平后的长度
print(wlen,nf, pad_length)
zeros=np.zeros((pad_length-wlen,)) #不够的长度使用0填补
pad_signal=np.concatenate((wave_data,zeros)) #填补后的信号
indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T #相当于对所有帧的时间点进行抽取，得到nf*nw长度的矩阵
indices=np.array(indices,dtype=np.int32 #indices 为pad_signal的位置
frames=pad_signal[indices] #得到帧信号
win=np.tile(winfunc,(nf,1)) #window窗函数，这里默认取1
return frames*win #返回帧信号矩阵
Df=5 #采样点时间间隔
fs=8000 #采样频率
N=fs/Df #采样点数
t = np.arange(0,(N-1)/fs,1/fs) #取样时间
wave_data=np.sin(2*np.pi*200*t) # 待处理的信号
#预加重
b,a = signal.butter(1,1-0.97,'high')
emphasized_signal = signal.filtfilt(b,a,wave_data)
#归一化
lifts=[]
for n in range(1,13):
lift =1 + 6 * np.sin(np.pi * n / 12)
lifts.append(lift)
#分帧、加窗
winfunc = signal.hamming(256) #汉明窗
#分帧函数：每帧长度256 ，不重叠长度80：
X=enframe(wave_data, 256, 80, winfunc)
frameNum =X.shape[0] #获取分帧后的帧数
下面是三组代码及结果
（1）第一种情况：
for i in range(frameNum):
y=X[i,:]
yf = np.abs(np.fft.rfft(y))
melM = librosa.feature.mfcc(y,fs,S=yf,n_mfcc=12)
print(melM)
输出结果：
[12.01187176 16.44704344 15.61216672 14.38118246 12.77777051 10.84203157 8.62860918 6.19697322 3.61622903 0.95595607 -1.70837193 -4.30547829]
（2）第二种情况：
for i in range(frameNum):
y=X[i,:]
yf = np.abs(np.fft.rfft(y))
yf=yf**2
melM = librosa.feature.mfcc(y,fs,S=yf,n_mfcc=12)
print(melM)
输出结果：
[ 570.37065963 795.19627788 761.25293536 705.78283689 630.39580985 537.27662246 429.12048164 309.05311082 180.53860667 47.27705996 -86.90415637 -218.16302183]
这两种情况，MFCC对应的图的形状是相同
（3）第三种情况：
for i in range(frameNum):
y=X[i,:]
melM = librosa.feature.mfcc(y,fs,S=None,n_mfcc=12)
print(melM)
结果：
[[-278.17089678] [ 144.78769227] [ 56.44110736] [ 40.11968822]
[ 25.1757355 ] [ 14.67492614] [ 5.77245624] [ -0.96605092]
[ -7.23574856] [ -12.32471135] [ -17.40294169] [ -20.9566881 ]]
（4）第四种情况：
for i in range(frameNum):
y=X[i,:]
yf = np.abs(np.fft.rfft(y)) #频谱取模
yf = yf**2 #谱线能量
nfilt = 24 #Mel滤波器数量
low_freq_mel = 0
NFFT=256
high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700)) # 把 Hz 变成 Mel
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # 将梅尔刻度等间隔
hz_points = (700 * (10**(mel_points / 2595) - 1)) # 把 Mel 变成 Hz
bin = np.floor((NFFT + 1) * hz_points / fs)
fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))
for m in range(1, nfilt + 1):
f_m_minus = int(bin[m - 1]) # left
f_m = int(bin[m]) # center
f_m_plus = int(bin[m + 1]) # right
for k in range(f_m_minus, f_m):
fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
for k in range(f_m, f_m_plus):
fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])
filter_banks = np.dot(yf[0:129], fbank.T)
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks) # 数值稳定性
filter_banks = 10 * np.log10(filter_banks) # dB 此为分贝标准公式
filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8)
#DCT系数
num_ceps = 12
c2 = dct(filter_banks, type=2, axis=-1, norm='ortho')[ 1 : (num_ceps + 1)] # Keep 2-13
c2 *= lifts
print(c2)
结果是：
[ 152.16939829 103.18780826 53.15848861 -3.0957145 -67.3315891
-123.00839875 -157.48879067 -150.09115643 -104.25335171 -47.27345817 -6.68998885 1.39787529]
