やりたいこと

pydubで読み込んだ音声信号に対し、Sliding Windowしながら周波数スペクトラムを表示したい。

周波数について

人間は音を聞き取る際、音声信号をなす波形情報の周波数を使用しています。
f:id:nsr_9:20210927141024p:plain

世の中で聞こえてくる美しい音色や人の声、騒音などのすべての音は、これらの周波数を持つ信号源が合成して出来ています。
f:id:nsr_9:20210927141528p:plain

ある信号に含まれる周波数成分を解析する手法としてフーリエ変換があります。
フーリエ変換は、その理論の成り立ちと無限に応用範囲があることから、情報工学の世界では最終奥義みたいな存在になっています。

やっている事はとてもシンプルで、連続した周期信号は正弦波の合成に分解（近似）出来るというものです。
f:id:nsr_9:20210927142328p:plain

先程のすべての音はあらゆる音の合成信号であるという話と対応している話なので、直感的に理解することが出来ますね。
フーリエ変換は、理論自体面白さや高速化技術（バタフライ演算）、画像や音声への応用事例など、非常に興味深いトピックスが沢山あるのですが、今回はnumpyで簡単に触れる所にとどめておきます。
フーリエ変換の基本は以下の書籍がとっつきやすいです。高専時代に大変お世話になりました。

マンガでわかるフーリエ解析 | 渋谷道雄, 晴瀬ひろき, 晴瀬ひろき, トレンド・プロ, トレンドプロ | 数学 | Kindleストア | Amazon

numpyである信号をフーリエ変換する際は次のように実装します。

import numpy as np
data = np.random.radnom([1024])   # 長さ1024の適当な信号源
F = np.fft.fft(data)

np.fft.fftで返ってくるFは信号長が入力と同じ1024の複素配列です。
複素配列のreal成分に信号の振幅成分、imag成分に振幅成分が含まれています。

Fから振幅スペクトラム、周波数スペクトラムを得るには次のように実装します。

amp = np.abs(F)
phase = np.angle(F)

この際、フーリエ解析できる周波数の最小単位は入力信号長に依存します。(今回のケースでは1024)
Fの各要素に対応する実際の周波数は、np.fft.fftfreqで求める事が出来ます。
f:id:nsr_9:20210927151810p:plain

window_len = 1024 # FFTに与えるベクトルの長さ
dt = 1/44100    # サンプリング周波数の逆数（波長）
freq = np.fft.fftfreq(window_len, dt)  # 各要素に対応する周波数

実際の実装

前回実装したらsliding windowのプログラムと組み合わせてみました。

nsr-9.hatenablog.jp

import pydub
import numpy as np
from pydub import AudioSegment
import matplotlib.pyplot as plt


if __name__ == "__main__":

    song = AudioSegment.from_wav("audio.wav")
    song_array = np.array(song.get_array_of_samples())
    song_ch1 = song_array[::2]
    song_ch2 = song_array[1::2]

    wsize = 1024
    frame_rate = song.frame_rate

    vmax, vmin = 0, 0
    for index in range(song_ch1.shape[0]//wsize):
        fig, ax = plt.subplots(1, 1)
        data = song_ch1[index*wsize: (index+1)*wsize]
        F = np.fft.fft(data)
        freq = np.fft.fftfreq(wsize, 1/frame_rate)


        amp = np.abs(F/(wsize/2))
        ax.plot(freq[1:int(wsize/2)], amp[1:int(wsize/2)])
        ax.set_xlabel("Freqency [Hz]")
        ax.set_ylabel("Amplitude")
        ax.grid()

        plt.savefig("./out/fig_{0:06d}.jpg".format(index))