    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)

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