定Q変換 (CQT)は歌声・音楽の音高／メロディ／和音を分析するのに便利な周波数分析方法です。
例えば，CQTスペクトログラムを算出すれば各音高（ドレミファ…）がどのくらいの強さで鳴っているかを時系列で観察することができます。
以下の記事では、定Q変換のリアルタイム向け実装としてパーセバル定理と疎行列計算を使った定Q変換の仕組みと実装方法を紹介しました。
しかし、この方法は疎行列計算ライブラリが使えないような環境では非常に低速になってしまいます。
そこで、以下の論文で紹介されている、疎行列計算を利用しなくても高速に計算できる再帰的ダウンサンプリング法による定Q変換を実装してみました。
実は、Pythonの音楽信号分析モジュール LibROSA の定Q変換メソッド librosa.cqtもこの手法に基づいています。
この記事では，再帰的ダウンサンプリング法による定Q変換の簡易実装とその原理を説明します．
再帰的ダウンサンプリング法の原理
従来手法の課題
以下が再帰的ダウンサンプリング法が提案された論文です。
"Constant-Q transform toolbox for music processing." , 2010.
定Q変換の基本的な仕組みはこちらの記事をご覧ください。
上記記事の実装では，最終的に以下の行列計算していました。
この式でCQTスペクトルを算出する時の問題点として、スペクトルカーネル行列の横幅が大きいため，スペクトルカーネル行列と分析信号のFFTスペクトルベクトルの積を計算するのに時間がかかります。
それに対して，前回の実装ではスペクトラムカーネル行列が疎（ほとんどの要素が零に近い値）であることから、疎行列計算ライブラリを使うことで実用的な計算量で処理できました。
しかし、利便性を考えると疎行列計算ライブラリを使わなくても実用的な計算量で処理したいところです。
そこで提案されたのが再帰的ダウンサンプリングによる定Q変換です。
再帰的ダウンサンプリング
この手法の肝は、CQTスペクトル計算の行列-ベクトルの積を1オクターブごとに処理するところです。
ブロック図で見ると、イメージしやすいと思います。
従来手法の計算量がネックなのはスペクトルカーネル行列の横幅＝分析信号が長いことが原因でした。
分析信号が長い理由は、低域の音高の周波数分解能を高域の音高と同じにする，定Q変換の性質のためです。
結局のところ、実は高域の音高を算出するのに使う分析信号の長さはかなり短いので、従来手法は周波数が高い音高ほど無駄な行列計算をしていることになります。
そこで、低域の音高に対する計算では分析信号をデシメーションする（間引く）ことで、スペクトルカーネル行列の横幅を小さくすることができます。
しかも、高域から低域に向かって1オクターブずつ再帰的に処理を行うことで，スペクトルカーネルを使い回すことができます。
この計算は以下の式で表すことができます。ここで、dがどのオクターブを処理しているかを示します。
イメージとしては、に最初は一番高域の１オクターブを処理します。
分析信号 x の一部だけを使っていることに注意してください。
以下は次のオクターブを処理する時のイメージです。
この処理では、最初の１オクターブ処理に使った分析信号 x の一部の，2倍の長さの信号を切り出して使います。
その切り出した信号に対してローパスフィルタリングと間引きを行うことで，最初の処理に用いた分析信号 Xd と同じ長さになるので、スペクトルカーネルを使い回すことができます。
3つ目以降の1オクターブも同様な処理となっていますが、2つ目の1オクターブに対する処理でローパスフィルタ→間引きを行った信号に対して再度ローパスフィルタ→間引きを行うことに注意してください（おそらくこの部分が最も分かりにくい部分だと思います。。）
実装解説
再帰的サブサンプリング法の実装は，基礎的な部分はパーセバル定理＋疎行列計算の定Q変換の実装方法と同じなので、その解説をした以下の記事も合わせてお読みください。
パラメタ準備 ～ CQT中心周波数ごとの窓長Nkの算出
上記記事とほとんど同じですが、lpf_orderでローパスフィルタの次数を与えています。
class RDCQT(): def __init__(self, sr, n_octave=8, # num. of octaves n_bpo=12, # num. of bins per octave. f0=55/2, # lowest center frequency window=np.hanning, # window func. lpf_order=6, # order of lowpass filter using recursive downsampling hop_length=128, winlenmax=None, q_rate=1.0, sp_thresh=0.0, ): # Params self.sr = sr self.nyq = sr / 2 self.n_octave = n_octave self.n_bpo = n_bpo self.hop_length = hop_length # calc. center freqs. self.cqt_freqs = self.compute_cqt_freqs(f0) self.n_pitch = len(self.cqt_freqs) # calc. Q and Nk(winlens) Q = q_rate / (2 **(1.0 / self.n_bpo) - 1.0) self.winlens = np.ceil(Q * sr / self.cqt_freqs).astype(np.int64)
低域通過フィルタ係数算出
低域通過フィルタは同じものを使い回すのでフィルタ係数を算出しておきます
正規化周波数の半分で高域をカットしていきたいので0.5をセットしています。
IIRフィルタを使用する場合はbutter()バターワースフィルタやbessel()ベッセルフィルタなどの利用が考えられます。
self.lpf_order = lpf_order
#self.lpf_b, self.lpf_a = butter(self.lpf_order, 0.5, btype='low')
self.lpf_b, self.lpf_a = bessel(self.lpf_order, 0.5, btype='low')
1オクターブ分のテンポラル／スペクトラムカーネル算出
まず、分析信号全体の長さn_fftとカーネルのサイズn_fft_1octを算出します。
このn_fft_1octは一番高域の1オクターブ内で最も長い窓幅に相当します。
パーセバル定理＋疎行列計算の定Q変換ではn_fftがカーネルの列幅でしたが、再帰的ダウンサンプリング法ではn_fft_1octになるのでだいぶ小さくなります。
なお、n_bpoは1オクターブを何分割するかを決める変数で、12音平均律なら12です。
# calc. kernel size self.n_fft_exp = len(bin(int(self.winlens[0])))-2 self.n_fft = 2 ** (self.n_fft_exp) self.n_fft_1oct = 2 ** (self.n_fft_exp - (self.n_octave - 1)) # 1-octave temporal/spectral kernel self.temporal_kernel = np.zeros((self.n_bpo, self.n_fft_1oct), dtype=np.complex128) self.spectral_kernel = self.temporal_kernel.copy()
次に、1オクターブ分のテンポラル／スペクトラムカーネルを算出します。
パーセバル定理＋疎行列計算の定Q変換でのカーネル算出と同じ計算ですが、一番高域の1オクターブ分しか計算していません。
# calc. 1-octave kernel coef = 2.0 * np.pi * 1j for m, k in enumerate(range(self.n_pitch-self.n_bpo, self.n_pitch)): Nk = self.winlens[k] Fk = self.cqt_freqs[k] #print(m, k, int(Nk), int(Fk), self.n_fft, self.n_fft_1oct) st = int((self.n_fft_1oct - Nk) / 2) en = Nk + st tmp_kernel = np.exp(coef * (Fk/sr) * np.arange(0, Nk)) self.temporal_kernel[m, st:en] = tmp_kernel * (window(Nk) / Nk) self.spectral_kernel[m] = np.fft.fft(self.temporal_kernel[m]) self.spectral_kernel = self.spectral_kernel.conjugate() / self.n_fft
信号に対するCQT処理
まず、全体のソースコードです。
_rdcqtが再帰的サブサンプリング法のサブルーチンになります。
def cqt(self, x): # zero clear self.x_list[:, :] = 0.0 self.tmp_x[:] = 0.0 self.y[:] = 0.0 # set x length = len(x) if length < self.n_fft: self.x_list[0, 0:length] = x self.x_list[0, length::] = 0.0 else: self.x_list[0, :] = x[0:self.n_fft] # subroutine _rdcqt(self.spectral_kernel, self.xd_fft, self.x_list, self.tmp_x, self.y, self.n_fft, self.n_fft_1oct, self.n_bpo, self.n_octave, self.lpf_b, self.lpf_a) return self.y[:]
再帰的ダウンサンプリング法
再帰的ダウンサンプリング法では、ここの処理が少し複雑です。
私の実装では実装しやすさを重視し、
に処理を分けました。
この実装では、ローパスフィルタリングはscipy.signal.lfilterを使っています。
もしオフライン処理で位相歪みが気になる場合は，ゼロ位相フィルタ (scipy.signal.filtfilt)を使うこともできます．
なお，間引き係数2の間引き処理は numpy のスライス機能で[::2]としています。
def _rdcqt(spectral_kernel, xd_fft, x_list, tmp_x, y, n_fft, n_fft_1oct, n_bpo, n_octave, lpf_b, lpf_a): # 1. recursive downsampling for k in range(1, n_octave): idx1 = n_fft // 2 ** (k-1) idx2 = idx1 // 2 tmp_x[0:idx1] = lfilter(lpf_b, lpf_a, x_list[k-1, 0:idx1]) #filtfilt #tmp_x[0:idx1] = x_list[k-1, 0:idx1] # w/o lpf x_list[k, 0:idx2] = tmp_x[0:idx1][::2] # 2. calc. cqt spectrum using Parseval's theorem for k in range(0, n_octave): center = n_fft // 2 ** (k + 1) st = center - n_fft_1oct//2 en = center + n_fft_1oct//2 st2 = n_bpo * (n_octave-k-1) en2 = st2 + n_bpo xd_fft[:] = np.fft.fft(x_list[k, st:en]) y[st2:en2] = np.dot(spectral_kernel, xd_fft)
実行結果
左がパーセバル定理＋疎行列計算の定Q変換，右が今回の実装（再帰的サブサンプリング法）です。
私の計算機環境とPython実装の場合、再帰的サブサンプリング法のほうが1～1.5倍くらい早いようです。
もう少し早くなるかと思っていましたが、lfilterの計算に時間がかかっていることが確認できています。
また、以下の正規化したスペクトログラムの結果は目検ではほぼ一緒ですが，スペクトログラム全体やオクターブごとのスケールが異なるため、結果を完全に一致させたい場合にはスケール調整が必要です。
チャープ信号
各音高のサイン波
米津玄師 - Lemon サビ冒頭
ソースコード
まとめ
再帰的ダウンサンプリング法による定Q変換 (CQT) のPython 実装と仕組みを解説しました。
パーセバル定理＋疎行列計算の定Q変換と比較すると、処理速度は若干良くなり、疎行列計算ライブラリが必要ありません。
今回は同じく論文で提案されている逆定Q変換は実装していませんが、分析用途であればこの実装で十分だと思います。
補足
カーネルの比較
ナイーブな定Q変換のカーネル
再帰的サブサンプリング法 定Q変換のカーネル
ローパスフィルタ→間引きが上手くできてない場合のスペクトログラム
ローパスフィルタリング周りがうまく実装出来ていない場合は、図のような折り返しひずみが出てしまいます。
