Pythonで音のSTFT計算を自作!スペクトログラム表示する方法

  • このエントリーをはてなブックマークに追加
Advertisements

FFTは時系列データの周波数分析に重宝しますが、時間の変化はわかりません。スペクトログラムは周波数と時間変化の両方を分析可能な手法です。しかし計算方法は少々手順を踏みます。ここではPythonによるコーディングを1つずつ解説して行きます。

こんにちは。wat(@watlablog)です。
スペクトログラムは特に音声データ分析で絶大な効果を発揮する手法です。
ここではPythonによるスペクトログラムコードを自作したので紹介します

ここでは周波数、時間、振幅の変化がひと目でわかるスペクトログラムを学びます。

しかし、Pythonの基本的なプログラミング手法を細かく説明しているわけではありません。もしPythonの基本文法等の学習をしたい場合はPython特化型学習サービスの「PyQ(パイキュー)」が大変オススメです!

スペクトログラムは時系列データから多彩な情報を得ることができる!

スペクトログラムとは?

スペクトログラムとは、音声をはじめとする様々な時系列データの分析に使われる可視化手法のことです。

このスペクトログラムは時間×周波数×振幅レベルという3つの要素で構成された3次元グラフです。

このページの後半ではPythonでスペクトログラムを計算するコードを載せますが、以下にそのプログラムで得られるグラフの例を先に紹介します。

スペクトログラム表示する前の入力データはこちら。
この波形はグランドピアノの「A(ラ)」の音を最低音から最高音までオクターブ毎に順番に打鍵していったものです。

分析する時間波形

wavファイルの音声をお聞きしたい方は、以下のメディアプレーヤーから再生してみて下さい。
※生のグランドピアノの音が再生されます。

そして以下が上の時系列データをスペクトログラム表示した結果です。

スペクトログラムの例

このように、スペクトログラム表示をすることで時系列データだけでは時間と振幅の関係しかわからなかったものに、周波数軸という奥行きを持たせることができます。

キレイだ!

この時系列データの処理に今回はSTFT(短時間フーリエ変換)という手法を用いました。

STFT(短時間フーリエ変換)とは?

STFTShort-Time Fourier Transform : 短時間フーリエ変換)とは、時系列データを細かく窓(フレーム)で区切ってFFTし、その窓を順次横にずらしながらFFTを繰り返していく解析方法です。

以下は時系列データからフレームでデータを切り出し、順次FFTしていくSTFTからスペクトログラムまでの過程を示す図です。

STFTとスペクトログラムの概要

このように時系列データを抜き出す際にフレームをオーバーラップさせて計算するのが一般的です。

フレームで抜き出してきた短時間のフーリエ変換の結果は、時系列順に並べることで3次元の配列データになりますが、FFTの振幅レベルを色に変換することで横軸時間、縦軸周波数の2次元表示が可能になります。

スペクトログラムから得られる情報

スペクトログラムは世間で広く時系列データの分析に使われていますが、一般的にどんな情報を読み取っているのかを紹介します。

通常のFFTが周波数と振幅レベル(+位相)という情報を持っているのに対し、スペクトログラムは時間の情報がプラスされています。

そのためスペクトログラムで主に分析する観点は「時間変化」になります。

以下は振幅レベルと周波数の時間変化の例です。

スペクトログラムで分析する内容

左の図は各ピーク(振幅レベルが卓越している部分)の振幅レベルが小さくなっていく(色が青くなっていく)時間がそれぞれ異なることを示しています。

右の図は時間が経つにつれピークが周波数方向に変化していくことを示しています。

このようにスペクトログラムではFFTで得られる情報に加えて、時間変化の情報を得ることができます。

Pythonでスペクトログラムを計算する方法

Advertisements

それではいよいよPythonによるSTFT計算とスペクトログラム表示までのプログラムを紹介します。

…とは言っても、実はWATLABブログでこれまで習得してきた内容の集大成となっていますので、それぞれ関連の記事を紹介していきます!

例により、ここに記載のプログラムは僕個人が学習して作った自作品です。間違いがある可能性もありますので使用には注意して下さい。

importするパッケージ

今回のSTFTとスペクトログラム表示のために、Pythonの外部ライブラリパッケージを宣言します。

インストールをしていない場合はpipで簡単にインストール可能なので、「Pythonのパッケージ管理ツール pipの使い方とコマンド集」の記事を確認して下さい。

配列処理のためのNumPyと科学技術計算のためのSciPy、wavファイルを使うためのPySoundFileを用意します。

時間波形データの用意(wavファイル読み込み)

まずは分析する時系列データを準備します。

具体的な方法は「PythonのPySoundFileでwavファイル波形表示」の記事で紹介した通りですが、コードを以下に示します。

オーバーラップ切り出し処理

時間波形を切り出す時にはオーバーラップ処理をできるようにしておきます。

オーバーラップ切り出し処理については「PythonでFFTをする前にオーバーラップ処理をしよう!」の記事に詳細を書いていますので参照下さい。

但し、上記記事とは戻り値に「最終時刻:final_time」を設定している部分に違いがあります。

最終時刻とは、以下の図に示すように、フレームをオーバーラップして切り出してきた時の最終フレームの終端における時刻です。

この時刻はスペクトログラムの時間軸を作成する時に使います

フレームサイズを決めてデータを切り出していくと、データ長とフレームサイズ、オーバーラップ率によってはフレーム切り出しできない余った部分がでます。

スペクトログラム表示の正確な時間軸を作るためには、この余った部分を除外する必要があるのでこのような計算をしています。

最終時刻を算出

STFT計算

上記までで時間波形切り出しの関数は準備できたので、続いて周波数分析としてFFTを行うための関数を用意します。

窓関数と振幅補正係数

FFTをかけるために窓関数としてハニング窓の関数を用意します。詳細は「PythonでFFT!SciPyで窓関数をかける」をご確認下さい。

そして窓関数を使って振幅レベルを適正に評価するためには、振幅補正係数acf(Amplitude Correction Factor)が必要です。acfの詳しい内容は「窓関数使用時の補正!FFTの時に忘れがちな計算とは?」に書きました。

デシベル(dB)変換

スペクトログラム表示をする場合は色成分をログスケールやdBスケールで表示させないと、音のような広い周波数帯域のデータを包括的に分析するには不向きになります(ピークの振幅レベルが10倍、100倍といった度合いで変化するのが普通なので)。

dBについての詳細は「Pythonで音圧のデシベル(dB)変換と逆変換式!」に記載しました。

聴感補正(A特性補正)

聴感補正は音のデータを分析する時に通常使用します。人間の耳の特性に近づけた分析をすることで、実際の活用に役立てます。

聴感補正に関しては、「Pythonで聴感補正(A特性)の曲線を作る!」で式と実装を解説しました。

FFT処理

FFTの関数は「PythonでFFT!SciPyのFFTまとめ」に記載した内容とほぼ同じですが、今回は聴感補正を考慮するので、最初からdB形式で計算するため、以下の順番で演算します。

①振幅補正係数acfをかける。
②dB変換する(聴感補正がdB値なので、dB同士のデータにする)。
③聴感補正曲線を足す。

ndarray形式であればforループを使わなくても配列演算ができるので、効率的ですね。

メイン実行文とグラフ表示

以上までのdef文で作成したコードを「function.py」というファイルで用意し、以下に示すメインコード(main_spectrogram.py)を実行することで冒頭に紹介したスペクトログラムを得ることができます(import文で関数をまとめたファイルを宣言する必要はあります)。

関数の使い方は「Pythonの関数 def文の使い方!引数や別ファイル式も解説」で説明しています。

wavファイルが無くても、自分で波形を作成(samplerate, x, dataを決めるだけ)して使っても動作しますので、是非お試し下さい。

追記:wavファイルが無いけどすぐ試したい方は自分で時間波形を作ろう!

メインのコードを変更

スペクトログラムはwavファイルでなくても描くことができます。上記のコードはモノラルのwavファイルを使うことが前提の内容でしたが、wavファイルって一回録音して、ファイル移動して…ってひと手間が面倒ですよね。

手っ取り早くお手元の環境でSTFT解析を行い、スペクトログラム表示を試したい方のために、以下にwavファイルが無くても動作するメインの.pyファイルを示します。

コードの冒頭でscipyをimportし、chirp信号という時間で周波数が変化する信号を生成しています。

チャープ信号については「Pythonでチャープ信号!周波数スイープ正弦波の作り方」に詳しく書いてありますので、是非併せてご覧下さい!

実行すると、以下のグラフが表示されると思います。5[s]間で0[Hz]から2000[Hz]まで周波数が増加している様子がはっきりと確認できますね!

チャープ信号のスペクトログラム

補足:カラーマップの種類について

スペクトログラムは色を使ってデータを表現します。

この色のパターンは非常に多くの種類があり、imshow関数内で引数としているcmapの指定で変更することが可能です。

以下の図はcmap違いのスペクトログラムです。「参考外部サイト様:color exsample code」に良いリストがありますので、お好みの設定を見つけて下さい。

cmap違いのスペクトログラム

まとめ

Advertisements

少し長くなってしまいましたが、本記事では音声データのスペクトログラムを計算する一連の方法を書きました。

しかし新しいことはほとんど無く、これまでWATLABブログで扱ってきた内容(wavファイルの取り扱い、波形切り出し、窓関数の取り扱い、dB変換、聴感補正、多段のFFT)が多くありました。

カラーマップを変更することでスペクトログラムの見え方を大きく変えることができました。

今回は音声データの強力な分析手法であるSTFTやスペクトログラム表示を習得しました!これで色々なデータから様々な情報を得ることができるようになるね!

Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています!

  • このエントリーをはてなブックマークに追加

SNSでもご購読できます。

コメント

  1. アバター 長谷川直樹 より:

    質問させてください!

    「PythonでFFTをする前にオーバーラップ処理をしよう!」でも記述してあったのですが、
    FFTをした後に平均化処理をする理由がよく理解できません。お教えしていただけたら幸い
    です。よろしくお願いします。

    1. wat wat より:

      watです。記事を読んで頂きありがとうございます。
      平均化ですが、ここに書いてある内容は主に実験データを想定しています。
      理想波形であれば1フレームのFFTも複数フレームを平均化したFFTも変わりません。
      しかし、実験データは時間によって主要次数の振幅が安定していないことも多々あります。
      平均化処理は、不安定なデータを綺麗にするために行います。
      おそらく大抵の商用FFTアナライザー(100万円くらいするやつ)にはこのような機能が付いていると思います。

      回答が的を射てなかったらすみません…。その場合はまた質問して下さい!

      1. アバター 長谷川直樹 より:

        なるほど!!

        理解することができました!ありがとうございます!

        1. wat wat より:

          ご連絡ありがとうございます!
          また何かわからないことがありましたらコメントよろしくお願いします!

  2. アバター 山田 凌大 より:

    突然失礼します。ここにあるコードでfftを行ったのですが、 data_array[i]= data_array[i] * han # 窓関数をかける
    ValueError: operands could not be broadcast together with shapes (4096,2) (4096,) というエラーが出て困っています。このエラーの解決方法について助言をいただけたら幸いです。

    1. wat wat より:

      コメントありがとうございます。
      過去に似たような質問を頂いたことがあるのですが、ステレオ録音して作成したwavファイルをご使用していますでしょうか?
      ここに記載のコードはモノラルのwavにのみ対応していますので、まずはそこをご確認下さい。
      dataを定義した直後に、print(data.ndim)とすることでdataの次元を調査することも可能です。その次元が2とかですとステレオのwavである可能性が高いです。

      また、もしwavファイルが原因で発生しているエラーでしたら、記事内の「追記:wavファイルが無いけどすぐ試したい方は自分で時間波形を作ろう!」の項目に記載してあるコードをお試し下さい。
      このコードはwavファイルを使わないでSTFT計算をするバージョンです。
      もしこのコードがエラーなく実行できれば、他は問題ないということになります。

      もし僕のこの回答が検討はずれである場合は再度コメント頂ければと思います!よろしくお願いします!

      1. アバター 山田 凌大 より:

        本当にありがとうございます!
        ご指摘の通り、ステレオ録音で作成したwavファイルだったので、
        data, samplerate = function.wavload(path) #wavファイルを読み込む


        x = np.arange(0, len(data)) / samplerate #波形生成のための時間軸の作成

        の間に
        ♯ステレオ2chの場合、LchとRchに分割
        data_l = data[:, 0]
        data_r = data[:, 1]
        ♯ 入力をモノラル化
        data = (0.5 * data_l) + (0.5 * data_r)

        を入れて、ファイルをモノラルに変えたら無事動きました。

        1. wat wat より:

          ご連絡ありがとうございます。
          解決したようでよかったです!また何かご不明点ありましたらコメント下さい!

  3. アバター 上田 より:

    質問失礼します。
    スペクトログラムに表示されているSPLのデータの値はどこに格納されているのでしょうか?
    データを一覧でみたいのですが。。。
    fft_array?

    1. wat wat より:

      ご訪問ありがとうございます。
      スペクトログラムで表示されているデータはご推察の通り、fft_arrayに入っています。
      def fft_ave()の結果のfft_arrayをfft_array.Tと転置させてax1.imshow(fft_array,…)とmatplotlibでプロットしています。

  4. アバター 上田 より:

    ご回答ありがとうございます!
    実際にfft_arrayをprintすると以下のようなります。
    [[-139.39209144 -127.47142483 -122.92904072 … -117.88843985
    -121.33788775 -125.57795728]
    [ -35.09114161 -35.37201986 -36.05485323 … -41.30974036
    -37.42815004 -37.43793617]
    [ -17.06772899 -17.50646536 -18.35489237 … -15.32351833
    -15.61655708 -15.99538017]

    [ 1.85881599 -2.20633606 -10.41866818 … 9.30049942
    3.95455108 -3.43244636]
    [ 6.967212 6.52847562 5.68004861 … 8.71142265
    8.41838391 8.03956081]
    [ 8.08926983 7.80839159 7.12555822 … 1.87067109
    5.75226141 5.74247527]]
    これはどのように読めばいいのでしょうか?
    []は右にしたがってX軸方向、二つ目の[]は上にy軸方向など。値の間隔はあるのかどうか。
    また、カラーバーは0〜60dbまでなのになぜマイナスの値が表示されるのでしょうか?(スペクトログラム上でマウスをドラッグした場合にもマイナスの値が見受けられました。)
    至らぬ点が多いですが、よろしくお願いいたします。

    1. wat wat より:

      スペクトログラムは2次元のデータなので確かにわかりにくいですよね。
      fft_array.Tと転置する前のfft_arrayの方で説明すると、これはFFTをかけていった結果をどんどん蓄積させているデータという意味を持ちます。

      >>これはどのように読めばいいのでしょうか?
      そのため、fft_array[i]と1行だけ抽出すると、それはスペクトグラムを縦方向(任意時間)でスライスした時のFFT波形を示します。
      [[0番目のFFT波形],[1番目のFFT波形],…,[i番目のFFT波形]]です。
      FFT波形はデータが周波数順になっています。
      fft_arrayには周波数軸や時間軸の情報が無いので、fft_axisysamplerateから自分で計算する必要があります。

      >>カラーバーは0〜60dbまでなのになぜマイナスの値が表示されるのでしょうか?
      これは実際にマイナスの値があるという意味で間違いありません。
      カラーバーは見やすいように手動で調整しているだけです。指定しなければオートレンジになると思います。

      以上、よろしくお願い致します。

  5. アバター 上田 より:

    ご回答ありがとうございます!
    おかげさまで理解することができました。ありがとうございます!

    1. wat wat より:

      ご連絡ありがとうございます。
      理解して頂けてよかったです。
      今後ともよろしくお願いします!

  6. アバター can より:

    質問失礼します!
    オーバーラップする意味とかは分かったのですが、スペクトログラムで表示する際はオーバーラップによって重複したフレーム部分はどう処理しているのかが分かりません。
    教えて下さると助かります。

    1. wat wat より:

      ご訪問ありがとうございます。
      オーバーラップ率を大きくすると、平均化して1枚の周波数波形を出す時には平均化する枚数が増えます。
      スペクトログラムを計算する時は、オーバーラップによって得られた各周波数波形は平均化せず、単純に時間軸方向に並べられます。

      まずここが、通常の平均化FFTとスペクトログラムの違いです。

      そのため、スペクトログラムでオーバーラップによって重複した部分は単純に周波数波形が増える事になり、これは時間分解能が上がる事を意味します。
      オーバーラップ率0%でスペクトログラムを計算して時間分解能(dt)が0.1[s]である場合、
      オーバーラップ率50%でスペクトログラムを算出するとdt=0.05[s]と半分になります。
      (時間に対する変化を捉えやすくなります)

      是非実際に数値を変更してスペクトログラムを眺めてみて下さい。
      よろしくお願い致します。

      1. アバター ゼローク より:

        時間分解能を上げる方法としてが最もな理由ですが、もう一つ理由が(確か)合って、窓関数をかけたことでのデータ圧縮を補完するためにオーバーラップで窓関数のゲインが高いところがなるべくくっつくようにするためでもあります。

        小野測器でハニング窓を使用している場合のオーバーラップ率を66.7%(2/3)以上と推奨しているのは、ハニング窓のスペクトル振幅の窓関数補正値が2/3であるからだったと思います。
        検索してもこの記述が無く、断言出来ないところが辛いところですが、ハニング窓を66.7%で書いていけば、私の言っていることが何となく伝わるかなと思います。

        1. wat wat より:

          ゼローク様、こちらにもコメントありがとうございます。
          オーバーラップの記事で絵を描いておきながらこちらではその説明を省いてしまいました。
          ただハニング窓の補正値と紐つけた説明は非常にわかりやすかったです。
          実務では伝統のようにそのルールに則っていましたが、今頭の中がクリアになったと思います。
          ありがとうございました。

  7. アバター MIO より:

    質問失礼いたします。
    上記のメイン実行文の通りにpythonに入力したのですが、スペクトログラムの画像が表示されません…。表示されないというか、グラフは出てくるけど色がついてないという状態です。エラーの箇所を見てみると、「NameError: name ‘fft_array’ is not defined」と出たり、「NameError: name ‘im’ is not defined」と表示されています。この場合は何処を直したらよいでしょうか?初心者的な質問ですみません…。

    1. wat wat より:

      ご訪問ありがとうございます。
      コピペしたコードの中に「fft_array = 〜」というコードは入っていますか?

      「〜is not defined」というエラーは変数を定義する前に使用した時に出ますので、文字が間違っているか、コピペ忘れ等が考えられます。

      エラーの発生している行より上をみて、変なスペースや文字間違えが無いかを確認してみて下さい。
      もし解決しない場合はより詳しい情報(Python環境(PyCharm使用、Python3.7で実行…エラー行は○○、fft_array=の付近のコードをここにコピペ…)等)を記載して頂ければ、何かしら力になれると思います!

      よろしくお願いします。

      1. アバター MIO より:

        返信ありがとうございます。
        python3.9.1で実行しております。
        エラーが起こっている箇所は、行数は分からないのですが、

        # データをプロットする。
        >>> im = ax1.imshow(fft_array, \
        … vmin = 0, vmax = 60,
        … extent = [0, final_time, 0, samplerate], \
        … aspect = ‘auto’,\
        … cmap = ‘jet’)
        Traceback (most recent call last):
        File “”, line 1, in
        NameError: name ‘fft_array’ is not defined

        >>> # FFTをかける
        >>> fft_array, fft_mean, fft_axis = function.fft_ave(time_array, samplerate, Fs, N_ave, acf)
        Traceback (most recent call last):
        File “”, line 1, in
        AttributeError: module ‘function’ has no attribute ‘fft_ave’
        >>>
        >>> # スペクトログラムで縦軸周波数、横軸時間にするためにデータを転置
        >>> fft_array = fft_array.T
        Traceback (most recent call last):
        File “”, line 1, in
        NameError: name ‘fft_array’ is not defined
        という箇所です。

        1. wat wat より:

          うーん、なんでしょう?これは初めてのお問い合わせですね。

          Pythonのバージョンは僕は3.7.7を使っています。
          物によっては新しいPythonバージョンにライブラリが対応していないという事はよくあるそうですが…今は不明です。

          確認ですが、
          ①function.pyとmain.pyと二つのファイルを同じフォルダ作っている。
          ②function.pyには以下の関数が順番に書いてある。
          def wavload(path):
          def ov(data, samplerate, Fs, overlap):
          def hanning(data_array, Fs, N_ave):
          def db(x, dBref):
          def aweightings(f):
          def fft_ave(data_array, samplerate, Fs, N_ave, acf):

          という所は大丈夫でしょうか?

          エラーの発生している場所がわからない場合はprint文を1行ずつ挟んで実行して確認するという手もあります。
          できればバージョンを揃えるのが良いと思いますが、まずは上記方法等でデバッグをしてみて下さい。

コメントを残す

*