Wizard Notes

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

定Q変換のリアルタイム向けPython実装の解説(音高・コード・メロディの分析向け)

f:id:Kurene:20210710162748p:plain:w500

定Q変換音楽信号の音高・コード・メロディ分析に相性の良い周波数分析手法です.

この記事では,前回の定Q変換 (CQT: Constant-Q Transform) の 解説 の内容をPythonで実装する方法を解説します.

Pythonの実装ですが,C++Javascriptなど,様々な言語での実装の参考になるように解説したいと思います.

なお、定Q変換の仕組みやパラメタについては以下の記事を参照してください。

www.wizard-notes.com

ソースコードGithubにおきました。

github.com

定Q変換の実装解説

要件定義とサブルーチン、ヘルパー

リアルタイムでも動くように,

  1. 事前に定Q変換の行列(スペクトルカーネル)を算出する
  2. 入力信号に対して,スペクトルカーネルを適用する
    1. 1回だけCQTする
    2. STFTのように短時間の信号に繰り返し適用する

の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クラス

MusicCQTCQTクラスを音楽向けに使いやすくしたクラスです。

ピッチクラス・オクターブ表記の入力を,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)

www.wizard-notes.com

テンポラル/スペクトルカーネルの算出

テンポラルカーネル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を使ったデモです。

www.youtube.com

www.wizard-notes.com