Wizard Notes

Pythonを使った音楽信号分析の技術録、作曲活動に関する雑記

定Q変換 (CQT: Constant-Q Transform) の 解説(音高・コード・メロディの分析向け)

はじめに

音信号の時間周波数分析にはFFT (高速フーリエ変換) STFT(短時間フーリエ変換がよく使われます.

しかし,FFTの性質上,音高に関わる分析をしようとすると実装や分析精度の向上に手間がかかります

そこで,この記事では音楽信号の分析によく利用される定Q変換 (CQT: Constant-Q Transform) と呼ばれる時間周波数分析の手法の仕組みを紹介します.

なお、本記事をベースにしたPython実装の解説がありますので、もしよろしければご覧ください。

www.wizard-notes.com

FFT / STFTでの音高分析の問題

音信号の時間周波数分析でよく使われるFFT (高速フーリエ変換) ) /STFT(短時間フーリエ変換)ですが,音高を分析することを考えると,

  1. FFTの各中心周波数は音高の周波数と対応が取れていない
  2. 周波数分解能の問題

という2つの問題があります.

この問題を理解するために,FFTの各周波数ビンの中心周波数を眺めてみましょう。

例えば,サンプリング周波数 8 kHz,信号長 256サンプルであれば,

0, 31.25, 62.5, 93.75, 125,156.25, 187.5, 218.75, 250, ..., 4000

というように中心周波数は等間隔に並んでいます。これを図で表すと以下のようになります.

f:id:Kurene:20210709004638p:plain:w400

図の上側は水平軸を線形に,下側は対数で表しています.

一方,12音平均律(ドレミファ…)の各音高の周波数をプロットしてみます.

f:id:Kurene:20210709005612p:plain:w400

音高では,線形軸では高周波数になるほど音高間の間隔が広がっているのが分かります.

従って,FFTのスペクトログラムを音高の観点で見ると,

  1. FFT周波数ビンの中心周波数と音高の周波数の対応がとれていない
  2. 低域では分解能が低く,高域では分解能が過剰に高い

という状態になっています.

実際の楽曲のスペクトログラムでも確認することができます.

米津玄師 - Lemon サビ

f:id:Kurene:20210709010934p:plain:w500

定Q変換のアイディア

FFTで 「低域では分解能が低く,高域では分解能が過剰に高い」現象はなぜ起こるのでしょうか。

小野測器 ー 周波数分解能はどのように決めるのか?」によると,FFTの周波数分解能\Delta fは、

\Delta f=\frac{f_{s}}{N}

となっています。ここでNは分析する時間窓のサンプル数です.

FFT各周波数で時間窓のサンプル数(長さ)が一定であるため,全ての周波数で同じ周波数分解能となっています.

定Q変換では,中心周波数に応じて時間窓のサンプル数(長さ)を変えることで周波数分解能を等比的を変えています

イメージとしては下の図のようになっています。

f:id:Kurene:20210709191110p:plain:w500

定Q変換の定式化

Constant-Q transform toolbox for music processing をベースに実装方法を解説します.

離散時間領域信号x(n)のCQTのスペクトル系列 X^{CQ}(k,n)は以下の式で求めます。

f:id:Kurene:20210710112951p:plain:w500

ここで,k = 0, 1, ... , K-1 は CQT の周波数ビン,n = 0, 1, ... , N-1 は 時間フレームを表し,\lfloor \rfloorは床関数による整数丸め)です.

a^\ast_k(n) a_k(n)複素共役であり,以下のように定義されます.

f:id:Kurene:20210710113913p:plain:w500

ここでCQT周波数ビンの窓の長さN_kは、

f:id:Kurene:20210710114952p:plain:w400

と定義されています.したがって,(2)の式は

f:id:Kurene:20210710115429p:plain:w500

とも書けます。

上記のQが quality factorと呼ばれ,CQT中心周波数f_kが高いときは窓幅を小さくし,f_kが低いときは窓幅を大きくして時間周波数分解能を一定に保つ働きをします.

なお,CQT の中心周波数は f_k

f:id:Kurene:20210710141800p:plain:w200

として等比的に決定します。

パラメタの整理

事前に決める必要があるパラメタを整理します。

  • B
    • 1オクターブを何分割するか
    • B=12 あれば12分割(12音平均律
  • f_0
    • 中心周波数の分割基準となる,最も小さい中心周波数
    • f_0 と B が決まれば,CQT中心周波数[f_{k}]が算出できます
    • 一般的な音楽データの音高分析であれば,[tex:f_0=440/2m]をよく使います
      • A4=440 Hz
  • q (0<q≦1)
    • 周波数分解能を調整するパラメタ
    • q を小さくすると周波数分解能は低下するが,窓幅 N_{k}が短くなり,計算量は少なくなる
    • 分解能としては,q=0.5, B=48とq=1, B=24は同じになる
  • w(n)
    • 窓関数 ハン窓,ハミング窓など

スパース行列演算による高速化

FFTの計算量は O(N log N) ですが,定Q変換を上記の式で実装すると計算量が多すぎるという問題があります。

具体的には,低域で窓幅 N_k が非常に大きくなってしまいます。

f:id:Kurene:20210710143303p:plain:w500

そこで,FFTを使った式変換と疎行列を使った計算によって計算量を減らす手法が提案されています.

まず,FFTを使った式変換はパーセバル定理を使います。

f:id:Kurene:20210710144251p:plain:w500

ここで,X, A^\astx, a^\astFFTした信号になります.a^\astはテンポラルカーネルA^\astはスペクトルカーネルとも呼ばれます。

最終的にCQTのスペクトル系列は以下の式で算出できます.

f:id:Kurene:20210710145107p:plain:w500

このスペクトルカーネルA^\astは事前に計算可能であり,非常にスパース(値がゼロの要素が多い)行列となっています.

f:id:Kurene:20210710120517p:plain:w400

従って、時間フレーム単位でみると,以下の図のようにスペクトルカーネル行列と信号のFFTスペクトルベクトルの積を疎行列演算すれば効率的にCQTスペクトルを算出できます.

f:id:Kurene:20210710153626p:plain:w500

定Q変換の改良手法

スパース行列演算で計算量を削減することができましたが,さらなる改良がおこなわれています.

  • 再帰的サブサンプリング方法
    • マルチレート処理風にオクターブごとに計算することで計算効率が向上
    • 実用的な逆変換を実現
    • librosa.cqt , librosa.icqt はこの方式
  • 変Q変換
    • 定Q変換では低域の窓幅が大きいことが問題
    • Q値を線形に変化させることで窓幅が大きくなるのを緩和
    • librosa.vqt

参考文献・参考ページ

関連手法

音高抽出手法としては,IIRフィルタバンクを使う手法もあります。

www.wizard-notes.com