numpyでスペクトル包絡の算出（ケプストラム分析）

・引用/データ

メル周波数ケプストラム（MFCC） - Miyazawa’s Pukiwiki

・スペクトル

import numpy as np
from scipy.io import wavfile
import matplotlib.pyplot as plt

# waveファイル入力
sample_rate, wave_data_orig = wavfile.read('a.wav')

# 切り抜き
cuttime = 0.04  # 切り出し長さ0.04[s]
N = int(sample_rate * cuttime)  # サンプル数1764
center = int(wave_data_orig.shape[0] / 2)
wave_data = wave_data_orig[center - int(N / 2): center + int(N / 2)]

# 離散フーリエ変換
fourier = np.fft.fft(wave_data)
amp = np.abs(fourier)

# 表示
freq = np.linspace(0, sample_rate, N)
plt.plot(freq[:int(N/10)], amp[:int(N/10)])
plt.show()

・対数スペクトル

解析対象の音声が、音源と声道フィルタのたたみ込みで生成されているとする。

時間領域のたたみ込みは、周波数領域では積になる。

対数をとると積は和になり、音源特性と声道特性の重ね合わせを可視化できる。

# スペクトルのlog
log_amp = np.log10(amp)
plt.plot(freq[:int(N/10)], log_amp[:int(N/10)])
plt.show()

・ケプストラム/スペクトル包絡

対数スペクトルをフーリエ変換したものがケプストラム。

　※スペクトルのフーリエ変換＝逆フーリエ変換

ケプストラムの高周波成分をカットして周波数成分に戻し、対数を外すと、声道特性（スペクトル包絡）が抽出できる。

# ケプストラム/スペクトル包絡
cps = np.fft.ifft(log_amp).real
cps[120:cps.shape[0] - 119] = 0
envelop = abs(10 ** np.fft.fft(cps))
plt.plot(freq[:int(N/10)], envelop[:int(N/10)])
plt.show()