python+numpy+scipyの組み合わせにより、DFTで簡単に周波数スペクトラムを求める事ができます。
ここでは一歩進んで、スペクトログラムを作成してみます。
スペクトログラムとは、STFT（短時間フーリエ変換）によって得られた周波数スペクトラムの振幅成分を画像化したものです。要は、数サンプルステップずつDFTの窓をずらしながら周波数スペクトルを求め、時系列に並べただけの画像です。
理屈は簡単ですが、Ｗaveファイルの入力、入力信号を指定ステップ、指定窓幅で切り出す処理、画像化等、本質ではない部分でコーディング量が多くなってしまいました。しかし、一旦ツール化しておくと、音声を周波数成分で分析する際に役に立ちます。
import sys
import os.path
import numpy
import scipy.fftpack
import scipy.signal
import scipy.io.wavfile
import argparse
def lpcms_int16_to_float(lpcms):
lpcm_array = []
for index, lpcm in enumerate(lpcms):
lpcm_array.append(numpy.divide(lpcm.astype(float), 32767.0))
lpcms = numpy.array(lpcm_array)
return lpcms
def read_wave_file_to_float(wave_file_name):
fs, data = scipy.io.wavfile.read(args.input)
nch = 1
if type(data[0]) == numpy.ndarray:
lpcms = data.transpose()
nch = len(lpcms)
else:
lpcms = numpy.array([data])
pcm_type_from_file = type(lpcms[0][0])
if pcm_type_from_file == numpy.int16:
lpcms = lpcms_int16_to_float(lpcms)
return fs, nch, lpcms
def split_lpcm(lpcm, span, width):
n = int(numpy.floor((len(lpcm) - width) / span))
extended_length = n * span + width
if len(lpcm) < extended_length:
lpcm = numpy.r_[lpcm, numpy.zeros(extended_length - len(lpcm))]
splited_lpcms = []
for index in range(0, n):
offset = index * span
splited_lpcms.append(lpcm[offset:offset + width])
return numpy.array(splited_lpcms)
def db_array_to_grayscale_bytes(db_array, lower_db, upper_db):
db_array = numpy.clip(db_array, lower_db, upper_db)
db_depth = upper_db - lower_db
grayscales = numpy.multiply(numpy.divide(numpy.subtract(db_array, lower_db), db_depth), 255)
return bytes(grayscales.astype(numpy.int8))
if __name__ == "__main__":
argment_parser = argparse.ArgumentParser()
argment_parser.add_argument("--input", "-i", type=str)
argment_parser.add_argument("--window_size", type=int, default=512)
argment_parser.add_argument("--span", type=int, default=128)
args = argment_parser.parse_args()
fs, nch, lpcms = read_wave_file_to_float(args.input)
window_half_size = int(args.window_size / 2)
fragments = split_lpcm(lpcms[0], args.span, args.window_size)
grayscale_bytes_array = []
window = scipy.signal.boxcar(args.window_size)
for fragment in fragments:
windowed = numpy.multiply(fragment, window)
spectrum = scipy.fftpack.fft(windowed)
amp = numpy.absolute(spectrum) / numpy.sqrt(args.window_size)
amp_db = 20 * numpy.log(amp + sys.float_info.min)
grayscale_bytes_array.append(db_array_to_grayscale_bytes(amp_db[:window_half_size], -120, 0))
basename = os.path.basename(args.input)
spectrogram_file_name = basename + "_spectrogram.pgm"
with open(spectrogram_file_name, "wb") as pgm_file:
pgm_file.write("P5\n".encode())
pgm_file.write(("{:d} {:d}\n".format(window_half_size, len(grayscale_bytes_array))).encode())
pgm_file.write("255\n".encode())
for grayscale_bytes in grayscale_bytes_array:
pgm_file.write(grayscale_bytes)
exit(0)
入力ファイル指定、切り出し間隔指定、窓幅指定ができる簡単なツールです。
適当にＷebから拾ってきたＷaveファイルを入力に与えて生成したスペクトログラムを以下に示します。
実際に出力する画像はPGM画像です。Mac OS上で作成しているので、ToyViewerで左９０度回転し、PNG保存しています。
画像例は、上下方向が周波数、左右方向が時間です。上が高周波成分、右が未来方向です。
DFTを施す前にかけている窓関数は矩形窓（boxcar）です。そのため、サイドローブ成分が多く出ており、タイミングによっては全周波数に成分が出ているように見えています。しかし、実際に全周波数に成分が強く存在するわけではありません。
この辺り、どういった分析をしたいかによって、DFTを施す前の窓関数を慎重に決める必要があります。
矩形窓は全窓関数中、メインローブが一番狭く周波数分解能は最高ですが、メインローブとサイドローブのダイナミックレンジが狭いため、上図のように、存在しない周波数成分が出てしまったり、強く出ている特定の周波数成分に他の周波数成分が強く影響を受けてしまったりします。
スペクトログラム生成処理の窓関数を他の窓関数に差し替えられるようにする事で、窓関数の違いを確認しやすくなるので、一度はスペクトログラムを生成するツールを作成してみることをお勧めします。