Wizard Notes

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

トップ > LibROSA > LibROSAにおける Salience （顕著性）スペクトログラムの算出方法の解説

LibROSAにおける Salience （顕著性）スペクトログラムの算出方法の解説

LibROSA Python 楽曲解析 歌声信号処理 音響音楽信号処理

はじめに

この記事では、 以下の記事で紹介 LibROSAのSalience スペクトログラムの算出方法を解説します。

www.wizard-notes.com

実際の計算方法、プログラムはlibrosa.salience() をご参照ください。

Salience スペクトログラムとは

Salience スペクトログラムは Salience Representation (顕著性表現)とも呼ばれます。

STFTやCQTで瞬時的な周波数を求めようとすると、周波数の解像度は低い周波数ほど低くなります。

Salience スペクトログラムでは、それを改良する工夫がされています。

librosa.salience() の処理の流れ

戦略

Salience スペクトログラムは、高調波（≒倍音）のエネルギーを考慮した振幅スペクトログラムです。

端的に言えば、例えば440Hzの周波数ビンに、880, 1320, 1760Hz の成分も足してやります。

そうすることで、周波数スペクトルの構造が基音と整数倍音によって成り立っている調波楽器の信号を検出しやすい信号（＝Salience スペクトログラム）を作ることができます。

従って、librosa.salience() の処理は以下のような流れになります。

  1. 各周波数における高調波エネルギーを算出する
  2. 各高調波エネルギーの信号の重み付き平均を取る*1
  3. ピークを検出する

f:id:Kurene:20201220173651p:plain

重要な処理だけを残したlibrosa.salience()のソースコードを以下に示します。

# librosa.salience()

# 周波数方向の線形補間により、各周波数の高調波エネルギーを算出する
# どの高調波エネルギーを計算するかは、h_rangeで指定する
S_harm = interp_harmonics(S, freqs, h_range, kind="linear", axis=0)

# 重み付き平均
S_sal = np.average(S_harm, axis=0, weights=weights)

# ピーク検出
# librosa では scipy.signal.argrelmax を利用
if filter_peaks:
   S_peaks = scipy.signal.argrelmax(S, axis=0)
   S_out = np.empty(S.shape)
   S_out.fill(fill_value)
   S_out[S_peaks[0], S_peaks[1]] = S_sal[S_peaks[0], S_peaks[1]]

S_sal = S_out
return S_sal

1. 高調波エネルギーの算出

高調波エネルギーの算出は、interp_harmonics()，harmonics_1d()で行います。

#harmonics_1d()

    # Note: this only works for fixed-grid, 1d interpolation
    f_interp = scipy.interpolate.interp1d(freqs, x,
                                          kind=kind,
                                          axis=axis,
                                          copy=False,
                                          bounds_error=False,
                                          fill_value=fill_value)

    idx_out = [slice(None)] * harmonic_out.ndim

    # Compute the output index of the interpolated values
    interp_axis = 1 + (axis % x.ndim)

    # Iterate over the harmonics range
    for h_index, harmonic in enumerate(h_range):
        idx_out[0] = h_index

        # Iterate over frequencies
        for f_index, frequency in enumerate(freqs):
            # Offset the output axis by 1 to account for the harmonic index
            idx_out[interp_axis] = f_index

            # Estimate the harmonic energy at this frequency across time
            harmonic_out[tuple(idx_out)] = f_interp(harmonic * frequency)

具体的には、ある周波数 frequency において、frequency * harmonic の信号を、frequencyに対応した周波数インデックスに格納しています。

イメージとしては、harmonicの値に応じて周波数シフトしたスペクトログラムを生成する処理となっています。

ただし、入力は離散的な信号であるため、frequency * h に対応する信号を生成する必要があります。すわなち、周波数方向の補完処理が必要になります。

LibROSAでは、scipy.interpolate.interp1dを使って補完処理をしています。

以下は、interp_harmonics()のみを適応した結果となります。

各周波数において高調波エネルギーが算出されていることが確認できます。

元のスペクトログラム： f:id:Kurene:20201220180626p:plain

高調波エネルギー： h=1 は、元のスペクトログラムに相当します。 f:id:Kurene:20201220180708p:plain

ソースコード

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


# librosa.salience 用のパラメタ
h_range = [1, 2, 3, 4]
weights = [1.0, 0.5, 0.33, 0.25]

filepath = "./miku_doremi_bpm120.wav"
y, sr = librosa.load(filepath, mono=True, offset=0.3, duration=7)


# 振幅スペクトログラムでの Salience スペクトログラム
n_fft = 1024
S = np.abs(librosa.stft(y, n_fft=n_fft))
fft_freqs = librosa.fft_frequencies(sr, n_fft=n_fft) # 各周波数ビンの中心周波数

# 高調波でのエネルギー計算
S_harm = librosa.interp_harmonics(S, fft_freqs, h_range, kind="linear", axis=0)

print(f"S_harm.shape: {S_harm.shape}")

librosa.display.specshow(librosa.amplitude_to_db(S, ref=np.max),
                               sr=sr, y_axis="log", x_axis='time')
plt.colorbar(format="%+2.0f dB")
plt.tight_layout()
plt.show()


plt.subplot(2,2,1)
librosa.display.specshow(librosa.amplitude_to_db(S_harm[0], ref=np.max),
                        sr=sr, y_axis="log", x_axis='time')
plt.title("h = 1")
plt.colorbar(format="%+2.0f dB")

plt.subplot(2,2,2)
librosa.display.specshow(librosa.amplitude_to_db(S_harm[1], ref=np.max),
                               sr=sr, y_axis="log", x_axis='time')
plt.title("h = 2")
plt.colorbar(format="%+2.0f dB")

plt.subplot(2,2,3)
librosa.display.specshow(librosa.amplitude_to_db(S_harm[2], ref=np.max),
                        sr=sr, y_axis="log", x_axis='time')
plt.title("h = 3")
plt.colorbar(format="%+2.0f dB")

plt.subplot(2,2,4)
librosa.display.specshow(librosa.amplitude_to_db(S_harm[3], ref=np.max),
                               sr=sr, y_axis="log", x_axis='time')
plt.title("h = 4")
plt.colorbar(format="%+2.0f dB")

plt.tight_layout()
plt.show()

librosa.display.specshow(librosa.amplitude_to_db(np.sum(S_harm, axis=0), ref=np.max),
                               sr=sr, y_axis="log", x_axis='time')
plt.colorbar(format="%+2.0f dB")
plt.tight_layout()
plt.show()

2. 重み付き平均

weightsで指定した重みで、各高調波エネルギーの重み付き平均を取ります。

LibROSAでは、numpy.average() を使っています。

高調波エネルギーを集約することにより、スペクトル漏れや低周波数の解像度の悪さを改善したスペクトログラムを得ることができます。

ピーク検出

音高情報としては、スペクトログラムはなるべくスパースであるほうがありがたいため。ピーク検出を行っているようです。

LibROSA では、scipy.signal.argrelmax() を使うことで、複数のピークの候補を残しているようです。

まとめ

LibROSAにおける Salience （顕著性）スペクトログラムの算出方法をまとめました。

参考

f:id:Kurene:20201220145459p:plain

*1:平均だけでなく、総和や中央値などの集約関数でもよい

Kurene

読者になる
コメントを書く