短時間フーリエ解析では、時間方向と周波数方向の分解能が固定であるため、高周波領域と低周波領域の両方に十分な解像度を得るのが難しくなっています。 この点を解決できるのが、ウェーブレット変換です。 基本的な考え方は、解析に用いる波形(ウェーブレット)を低周波を解析する場合には時間方向に引き延ばし、高周波領域を解析する場合には縮めることにより、各周波数領域において最適な時間分解能を得ることができるというものです。 離散信号に拡張した離散ウェーブレット変換は、効率的な符号化が得られることから画像圧縮などに用いられています。 しかし、音声波形の周波数解析という面からいうと、この離散ウエーブレット変換は解析に用いるのが困難です。なぜならば、解析に用いるウエーブレットの周波数成分が単一でないため、ある解像度のシグナルがどの周波数に一致するのかがわかりにくいからです。
可変窓を用いた高速再帰的スペクトル解析(著:中辻 秀人) amazon link で述べられてるアルゴリズムは、短時間フーリエ解析の欠点をウェーブレット変換と同様の手法で解決しつつ、周波数解析に適したものとなっています。
詳細なアルゴリズムや数式は、上記文献を参照していただくとして、ここでは概略を述べます。 基本的なアイデアは、ウェーブレット変換と同様に切り出し波を解析したい周波数に応じて変化させるというものです。 具体的には、周期数T=16回くり返すsin波とcos波を用意して、信号波形との内積を求めることによりスペクトルを計算します。 正直に計算をすると、各解析周波数の数と波形の長さと切り出し波の長さの積の計算が必要となり大変ですが、切り出し波の長さ及び周波数をうまく調整することにより、計算を減らすことができます。 切り出し波T=16周期の長さが、サンプリングの整数倍になる解析周波数を選ぶと、解析波形と切り出し波の内積計算の際に切り出し波形を使い回すことができるようになります。 新しく内積の計算に入ってくるサンプルと、区間から外れるサンプルのデータのみを計算対象にできるため、正直に計算した場合に比べ大幅に計算量を削減できます。 また、各周波数ユニットは独立であるため、分散計算を行うことにより高速化できます。
上記文献の巻末に、高速計算バージョンのプログラムが掲載されています。 このBasicコードを参考に、C++で書いたものを https://github.com/lithium0003/FMRS に置いています。
並列実装が可能とのことなので、cudaを用いて並列計算を行うように改良し、GTKで図が出るようにしたプログラムを https://github.com/lithium0003/sound_split に置きました。
https://github.com/lithium0003/sound_split からソースコードをダウンロードした上で、コンパイルには次の環境が必要です。
開発機は以下の条件です。Ubuntu 18.04.1 LTS / Intel Core i7-6850K / 64GB memory / GeForce GTX 1080Ti 11GB
linuxであればこれ以外の環境でも動くと思います。計算結果が大きくなるので、メモリに展開しています。ある程度のメモリがない場合はプログラムが落ちる場合があります。この場合は、バッファを減らすソースコードの修正をしてコンパイルするか、解析対象を短くしてください。
短時間フーリエ解析で用いた音源と同じものをこの手法で解析しました。
300Hz, 440Hz, 880Hzを混合しています。
110Hzから1760Hzまでスイープする信号です。
ようこそジャパリパークへ(どうぶつビスケッツ×PPP) 作詞/作曲:大石昌良 (amazon link) の冒頭部分です。
短時間フーリエ変換に比べて、時間分解能と周波数分解能共によい結果が得られています。 単純なsin波やスイープ信号では違いが見えにくいですが、楽曲などの複雑に音が混ざっているものに対してはよい分解ができていると感じられます。 20Hz-22000Hzまで0.5%増分で解析した場合1353unitsの解析チャンネルが得られます。半音の周波数差は6%なので十分に楽曲の音程が分解できます。 時間方向には、オリジナルのサンプリング周波数と同じサンプル数の解析結果が得られます。しかし、T=16の場合該当周波数における16波長分の音が継続しないと、解析応答が安定に達しません。これは440Hzにおいて36msであり、44.1kHzサンプリング周波数で1600samplesに相当します。
解析のサイドローブについて、T=8-64まで変化させた場合について見てみます。緑の実線が各周波数ユニットのスペクトル出力で、緑の影はT=16の場合で近似したサイドローブマスクです。 T=8の場合、時間分解能は非常によいのですが、メインローブの幅が太く周波数分解能があまりよくありません。 T=16の場合、メインローブの幅は細くなり周波数分解能はよくなります。サイドローブはかなりあります。 T=32の場合、メインローブの幅はさらに急峻になり、サイドローブも抑えられます。しかし、時間分解能が悪いです。 T=64の場合、さらにピークが急峻になりますが、時間分解能が非常に悪くなります。