定Q変換は音楽信号の音高・コード・メロディ分析に相性の良い周波数分析手法です.
この記事では,前回の定Q変換 (CQT: Constant-Q Transform) の 解説 の内容をPythonで実装する方法を解説します.
Pythonの実装ですが,C++やJavascriptなど,様々な言語での実装の参考になるように解説したいと思います.
なお、定Q変換の仕組みやパラメタについては以下の記事を参照してください。
定Q変換の実装解説
要件定義とサブルーチン、ヘルパー
リアルタイムでも動くように,
の2段階の構成とします。
cqt = MusicCQT(sr, sparse_computation=True, winlenmax=None) start_time = time.time() # 1度だけCQTする場合 spec = cqt.cqt(y[0:cqt.n_fft]) # 短時間 定Q変換 spectrogram = cqt.stcqt(y)
CQTクラスとMusicCQTクラス
MusicCQT
はCQT
クラスを音楽向けに使いやすくしたクラスです。
ピッチクラス・オクターブ表記の入力を,CQTで必要な変数に変換して渡します。
class MusicCQT(CQT): def __init__(self, sr, min_pitch="A1", n_octave=6, tuning=440, add_bins_end=0, chroma_labels = ["A", "A#", "B", "C", "C#", "D", "D#", "E", "F", "F#", "G", "G#"], **kwargs ): n_chroma = len(chroma_labels) n_bins = int(n_octave*n_chroma) # e.g. A0 (27.5Hz) is 0, A1(50Hz) is 12 (if tuning==440) calc_freq_from_pitch_index = lambda _k : tuning * 2 ** (_k/n_chroma - 4) min_pitch_class = min_pitch[0].upper() min_pitch_num = int(min_pitch[1]) min_pitch_index = chroma_labels.index(min_pitch_class) + n_chroma * min_pitch_num f1 = calc_freq_from_pitch_index(min_pitch_index) fmax = calc_freq_from_pitch_index(min_pitch_index + n_bins + add_bins_end) #print("f1, fmax", f1, fmax) super(MusicCQT, self).__init__(sr, n_octave=n_octave, f1=f1, fmax=fmax, n_bpo=n_chroma, **kwargs) self.min_pitch = min_pitch self.n_chroma = n_chroma self.labels = [ chroma_labels[k%n_chroma]+str(k//n_chroma) for k in range(min_pitch_index, min_pitch_index + n_bins)]
cqt.stcqt()
cqt.stcqt()
は、短時間フーリエ変換のように長時間の信号をずらしてcqt.cqt
を繰り返し適用します
class CQT: def stcqt(self, x): """ short-time constant q transform for batch """ n_frames = len(x) // self.hop_length x_tmp = np.r_[x, np.zeros(self.n_fft)] spec = np.zeros([n_frames, self.n_pitch], dtype=np.complex128) for k in range(n_frames): st = k * self.hop_length en = st + self.n_fft spec[k] = self.cqt(x_tmp[st:en]) return spec.T
CQTクラス
メンバ変数
CQTクラスのコンストラクタでは,CQTで利用するパラメタを計算・格納します。
class CQT(): def __init__(self, sr, n_octave=7, n_bpo=12, f1=55/2, fmax=None, window=np.hamming, sp_thresh=0.005, sparse_computation=True, hop_length=None, winlenmax=None, q_rate=1.0 # Schörkhuber, C., & Klapuri, A. , 2010. ): self.sparse_computation = sparse_computation self.sr = sr self.nyq = sr / 2 self.fmax = self.nyq if fmax is None else fmax self.n_octave = n_octave self.n_bpo = n_bpo # the number of bins per octave.
中心周波数,Q値,各周波数ビンの窓幅算出
スペクトルカーネルにおけるFFT利用に備えて,
最大の窓幅より大きい2のべき乗の値を n_fft
としています。
(self.n_pitch, n_fft)
がスペクトルカーネルの形状となります
# calc. center freqs. self.cqt_freqs = self.compute_cqt_freqs(f1) 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) self.n_fft = 2 ** (len(bin(int(self.winlens[0])))-2)
テンポラル/スペクトルカーネルの算出
テンポラルカーネルself.a
とスペクトルカーネルself.kernel
を算出します。どちらも要素は複素数であることに注意してください。
self.a = np.zeros((self.n_pitch, self.n_fft), dtype=np.complex128) self.kernel = np.zeros(self.a.shape, dtype=np.complex128) coef = 2.0 * np.pi * 1j for k in range(self.n_pitch): Nk = self.winlens[k] Fk = self.cqt_freqs[k] st = int((self.n_fft - Nk) / 2) en = Nk + st #tmp_a = np.exp(coef * Q * np.arange(0, Nk)/Nk) tmp_a = np.exp(coef * (Fk/sr) * np.arange(0, Nk)) self.a[k, st:en] = tmp_a * (window(Nk) / Nk) self.kernel[k] = np.fft.fft(self.a[k])
こちらの記事に書いた通り,tmp_a
の算出はコメントアウトしているQ
を使う式でもOKです。なお、Nk
の整数丸めの方法で結果が若干変わってきます。
疎行列計算
以下のブログを参考にさせていただきました。
適当な閾値で十分に小さい値のスロットを零とし、疎行列演算用のオブジェクトを生成しています。
【Python】 高速な Constant-Q 変換 (with FFT) - 音楽プログラミングの超入門(仮)
# prepare sparse computation self.kernel[np.abs(self.kernel) <= sp_thresh] = 0.0 if self.sparse_computation: self.kernel_sparse = csr_matrix(self.kernel) else: self.kernel_sparse = self.kernel self.kernel_sparse = self.kernel_sparse.conjugate() / self.n_fft
CQT実行部
FFTを使うので,入力信号の長さが2のべき乗になるように調整しています。
def cqt(self, x): length = len(x) if length < self.n_fft: self.x[0:length] = x self.x[length::] *= 0.0 else: self.x[:] = x[0:length] if self.sparse_computation: self.y[:] = np.fft.fft(self.x[:]) * self.kernel_sparse.T else: self.y[:] = np.dot(self.kernel_sparse, np.fft.fft(self.x[:])) return self.y
デモ
PyQtGraphを使ったデモです。