はじめに

f:id:Kurene:20210710233854p:plain

音信号を分析する時間周波数分析手法としてはSTFT（短時間フーリエ変換）が良く使われますが，特に音楽を対象にして音高、コード、メロディなどを分析する場合は定Q変換（CQT）という手法が便利です。

定Q変換を自分で実装するのは結構大変ですが，Pythonでは音楽信号分析モジュール LibROSA を使うことで、さくっと定Q変換を試すことができます。

はじめに
LibROSAでの定Q変換
逆定Q変換

LibROSAでの定Q変換

librosa.cqt を呼び出すことで、数行で定Q変換を実装できます。

import numpy as np
import librosa
import librosa.display
import soundfile as sf
import matplotlib.pyplot as plt


# 音源ファイル読込
sr = 44100
filepath = librosa.ex('trumpet')
y, sr = librosa.load(filepath, sr=sr)

# 定Q変換，振幅スペクトログラム
hop_length      = 128
n_octave        = 8
bins_per_octave = 12
C = librosa.cqt(y, sr=sr,
                hop_length=hop_length,
                n_bins=n_octave*bins_per_octave,
                bins_per_octave=bins_per_octave)  
M  = np.abs(C) 
librosa.display.specshow(librosa.amplitude_to_db(M, ref=np.max),
                         bins_per_octave=bins_per_octave,
                         sr=sr, x_axis='time', y_axis='cqt_note', cmap="jet")
plt.title('Constant-Q power spectrum')
plt.colorbar(format="%+2.0f dB")
plt.show()

librosa.cqt の重要な引数を列挙します．

hop_length
- 信号に定Q変換をする際にどのくらいずらしながら分析するかを指定します
- なるべく小さいサンプル数の方が時間変化を捉えやすくなりますが，計算量や出力されるCQTスペクトログラムのサイズが大きくなってしまいます
bins_per_octave
- 1オクターブを何分割するかを指定します
- 基本は12分割（12音平均律）でよいと思います
- 目的によっては12分割で解像度が足りないことがあるため、その場合は24分割にするなど、分割数を増やします
n_bins
- 全ての周波数ビンの数です。
- 分割数12の場合，n_bins/12が何オクターブ分析するかに相当します
fmax
- CQTの一番低い音高の周波数です。
- librosaの場合、fmin=librosa.note_to_hz('C2')というようにlibrosa.note_to_hz()と音高名の文字列を使って指定するのがおすすめです
tuning
- A4=440 Hz ではない基準ピッチの音楽信号を分析する場合、この引数を設定してCQT中心周波数を変えます
- librosa.estimate_tuning — librosa 0.8.1 documentation

逆定Q変換

逆STFTのように、一応逆定Q変換もあります。

しかし，CQTは、

周波数ビンによって分析する時間区間が異なる
周波数ビンの中心周波数が高くなるにつれて，分析窓長を小さくして周波数分解能を下げる

という性質があるため、逆の変換時に品質劣化が起こることに注意が必要です *1。

LibROSA では librosa.icqt として実装されています。

そこで，おそらく周波数解像度を上げれば逆変換の性能も上がると考え，bins_per_octaveを増やした時の性能変化を確認してみました．

sr = 44100
filepath = "miku.wav"
y, sr = librosa.load(filepath, sr=sr)

hop_length      = 128
n_octave        = 8
window          = "hann"
for k in range(8):
    bins_per_octave = 6 * (k+1)
    C = librosa.cqt(y, sr=sr,
                    hop_length=hop_length,
                    window=window,
                    n_bins=n_octave*bins_per_octave,
                    bins_per_octave=bins_per_octave) 
    y_hat = librosa.icqt(C=C, sr=sr, 
                         window=window,
                         hop_length=hop_length,
                         bins_per_octave=bins_per_octave)

    sf.write(f"cqt_icqt{bins_per_octave}.wav", y_hat, sr, format="WAV", subtype='PCM_16')
    error = 10 * (np.log10(np.sum(y**2)) - np.log10(np.sum((y[0:y_hat.shape[0]]-y_hat)**2)))
    print(f"{bins_per_octave},{error}")
    
sf.write("cqt_orig.wav", y,     sr, format="WAV", subtype='PCM_16')