More than 3 years have passed since last update.

音質評価技術鈴虫とコオロギの音色の差を検知する変動音解析

Python

Last updated at 2020-04-26Posted at 2020-04-26

はじめに

騒音評価は多くの場合、フーリエ変換などによって求めた周波数特性に基づいて行われます。この周波数特性は周波数成分ごとの強さを示すパワースペクトルが頻繁に用いられ、多くの場合位相は無視されます。伊藤毅(1955), 「音響工学原論上巻」，コロナ社にも位相は聴覚上重要ではないと記載されていた記憶があります(手元にないので記載ページはわかりません)。しかし、極端な話、無限の周波数成分が均一に含まれている信号の位相について、全てがランダムならホワイトノイズ、全て揃っていればパルス波になるように、位相によっては似ても似つかない音になります。近年では位相の重要性が見直されているようで、騒音計で有名な小野測器では従来の騒音評価とは異なる様々な手法が提案されており、音質評価と称しているようです。
今回はその中でも「変動音」と呼ばれる音の解析方法について学んだことを書きたいと思います。知覚騒音は門外漢なので間違っていることが多々あると思います。そのときはコメントで教えていただけるとありがたいです。

参考

こちらを読んで勉強しました。ここでは解析手順の記載はありますが、具体的にどのようなアルゴリズムで解析するかはわかりません。
小野測器:変動音解析とは

変動音解析に使うアルゴリズムはおそらくヒルベルト変換して包絡線を求め、それをFFTするのだと思います。
小野測器:計測コラム180号用基礎からの周波数分析(29)-「ヒルベルト変換と解析信号」

それで合っているのか、このコラムの内容と照らし合わせて確かめました。
小野測器:身近な計測虫の音

今回は簡単のため変動成分の単位を無視して2つの音の相対的な評価に留めました。更に詳細な内容は以下に記載があります。
小野測器:音質評価とは(page3)

変動音解析

変動音

変動音とは、例えば400Hzと404Hzの音叉を同時に鳴らすとうねりが聴こえるように、周期的に振幅が変動する騒音のことです。
ところで耳で聴くと4Hzのうねりが聴こえますが、この音をフーリエ変換しても400Hzと404Hzしか出てきません。このことからも音質評価はフーリエ変換して周波数特性を見るだけでは不足だということがわかると思います。

これは100Hzの信号に10HzでAM変調をかけたものです。

解析信号

中学で習う三角関数の和積の公式から上記のことを考えてみます。

\begin{align}
f(t) &=& A \left( \sin{(2 \pi 404 t)} + \sin{(2 \pi 400 t)} \right) \\
&=& 2A \sin{\left( 2 \pi \frac{404-400}{2} t \right)} \sin{\left( 2 \pi \frac{404+400}{2} t \right)}
\end{align}

和を積の形に直すと、4Hzのうねりの成分が出てきました。

フーリエ変換はある関数を三角関数の線形結合に変換するものですが、その周波数成分の強さを表すフーリエ係数は時間変化しません。一方で、和積の公式からわかるように、三角関数の和の形式で表される信号は積の形式にも直せます。すると、振幅と位相が時間変化する信号とも捉えることができます。ここで解析信号という考え方を導入します。実部に元の信号、虚部に位相が90度遅れた信号を持つ関数です。
$$f(t) = A(t) \left( \cos{(\theta (t))} + j \sin{(\theta (t))} \right)$$
この信号はヒルベルト変換して位相を90度遅らせて虚部につけるか、フーリエ変換して負の周波数を0にして逆変換することで得られます。

今回知りたいのは$A(t)$成分です。$\cos^2{\theta} + \sin^2{\theta}=1$の関係を使えば実部と虚部の2乗和平方根より、
$$A(t) = A(t) \sqrt{\cos^2{(\theta (t))} + \sin^2{(\theta (t))}}$$
と得られます。これは包絡線(Envelope)と呼ばれています。これをFFTすることで変動成分の周波数と強度を得られます。

import numpy as np

fft = np.fft.fft(sig)
fft[N//2:] = 0.0
envelope = np.abs(np.fft.ifft(2.0*fft))

fluctuation_spectrum = np.abs(np.fft.fft(envelope))