はじめに

この記事では、Onset 検出や局所／大域のテンポ分析などリズム分析においてよく使われる、Novelty Function についてPythonのコードとともに紹介します。実装と理解の助けになれば幸いです。

Novelty Curve/Function (Onset_envelope)

1. Energy based
• 時間領域で、ある時間区間範囲の振幅／パワー値の総和を算出
• それぞれ先の時間フレームとの差分を算出
• 半波整流波形化：　負の値（＝エネルギー減衰）は0に置換する
2. Spectral based
• 時間周波数領域に変換
• 周波数ごとに、ある時間区間範囲の振幅／パワー値の総和を算出
• それぞれ先の時間フレームとの差分を算出
• 半波整流波形化：　負の値（＝エネルギー減衰）は0に置換する
• 周波数で合算（例：平均）

LibROSAでも採用されている 2. は Spectral Fluxの算出と似ています。（ただし、負の値を０にしており、つまり音の立ち下がりについては無視しています。）周波数ごとに差分をとることで、各楽器の立ち上げりが捉えやすくなるのだと思います。

LibROSAなど、いくつかWeb上の実装を見てみると、 - エネルギー総和を算出する時間区間の長さ - log を取る（＝人間の聴感に合わせる） - (Spectral based) 周波数重み付け といった工夫がみられます。

実装と適用例

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

# Load music singal
filepath = librosa.util.example_audio_file()
y, sr = librosa.load(filepath, offset=30)

# Energy based novelty curve
mean_square = librosa.feature.rms(y=y)
energy_novelty = np.log( mean_square[1:]/ (mean_square[:-1] + 1e-9) )
energy_novelty[energy_novelty<0.0] = 0
energy_novelty /= np.max(energy_novelty)

# Spectral based novelty curve
Y = np.log(np.abs(librosa.stft(y))+1e-9)
spectral_novelty = None
n_freq, n_tf = Y.shape
spectral_novelty = np.zeros(n_tf-1)
for f in range(0, n_freq):
tmp = Y[f,1:] - Y[f,:-1]
tmp[tmp<0.0] = 0.0
spectral_novelty += tmp
spectral_novelty /= np.max(spectral_novelty)

plt.subplot(4,1,1)
plt.plot(y)
plt.title("Original music signal")
plt.subplot(4,1,2)
plt.plot(mean_square)
plt.title("RMS")
plt.subplot(4,1,3)
plt.plot(energy_novelty)
plt.title("Energy based novelty function (Onset Envelope)")
plt.subplot(4,1,4)
plt.plot(spectral_novelty)
plt.title("Spectral based novelty function (Onset Envelope)")

plt.tight_layout()
plt.show()

Spectral based Novelty function の方が、細かいリズムを捉えやすそうな信号となっていることが分かります。

PySoundFile vs audioread ?

>>> y, sr = librosa.load(filepath, sr=sr, duration=10, mono=True)
...Python36\site-packages\librosa\core\audio.py:146: UserWarning: PySoundFile failed. Trying audioread instead.

PySoundFileでは、様々なオーディオデータを扱うことができる C言語実装のライブラリ Libsndfile を使っています。audioreadFFmpegをバックエンドとしているため、導入に一手間必要＆処理的に非効率ということなのでしょうか…。

デベロッパーの Brian McFee 氏曰く、

bmcfee commented on 12 Mar

Quick benchmark: loading a 4-minute wave file brings the average load time from 92ms down to 8ms.

Speedup is less dramatic on compressed files (eg ogg), but everything seems to work as expected.

[CR needed] Soundfile-based loader #847

とのことなので、長い or 大量の楽曲信号を扱うときに嬉しい気がします。

2. AtomicParsley のインストール

コマンドラインで、MPEG-4メタデータを編集できるソフトウェアです。

Jupyter notebooksで音楽信号処理の基礎を学べるフレームワーク：FMP Notebooks

AudioLabs @ FAU の、音楽信号処理の権威、Meinard Müller 先生が、ご自身の著書 "Fundamentals of Music Processing (FMP)."に沿った教材として、Jupyter notebooks で音楽信号処理の基礎を学べるフレームワークFMP Notebooks を公開しています。

Meinard Müller and Frank Zalkow: FMP Notebooks: Educational Material for Teaching and Learning Fundamentals of Music Processing. Proceedings of the International Conference on Music Information Retrieval (ISMIR), Delft, The Netherlands, 2019.

これほどの内容が、しかもコード付きで、Webで公開されているのは非常に驚きです…！音楽信号処理について分かりやすく体系的に書かれている文献はあまりないので、是非とも活用していきたいものです。

調波打楽器音分離とは？

そこで、分析の前処理として、 打楽器の音と非打楽器（調波楽器）の音を分離する調波打楽器音分離 (HPSS: Hermonic/Percussive Source Seperation) *1が良く使われています。

この記事では、HPSSの概要とPython (LibROSA) のコードの解説します。

調波打楽器音分離のアイディア

フーリエ変換等を使って周波数に変換した楽曲信号を考えます。 します。*2

ただし調波／打楽器音は未知であるため、何らかの仮定をおく必要があります。

• 調波楽器音：　時間方向に連続的
• 打楽器音：　　周波数方向に連続的

であるという特徴を用いることで、調波／打楽器音を推定します。 （もしくは、調波／打楽器音について何らかの事前学習を行う手法もあります。）

アルゴリズムの設計

非負値行列因子分解＋基底クラスタリング

その基底をSVMサポートベクターマシン）で調波／打楽器音基底に2クラス分類します。そして、調波／打楽器音に分類された基底のみを用いて対応するアクティベーション行列との行列積を求めることで、調波／打楽器音のみの信号のスペクトログラムを推定します。なお、基底分類用の2クラスSVMの学習を事前学習しています。

Separation of drums from polyphonic music using non-negative matrix factorization and support vector machine - ResearchGate

最適化問題（行列因子分解）として解く

• 調波／打楽器音の時間／周波数方向の連続性に基づくコスト項
• 等式制約
• 非負値制約

Separation of a monaural audio signal into harmonic/percussive components by complementary diffusion on spectrogram - ResearchGate

メディアンフィルタベースの手法

Harmonic/Percussive Separation using Median Filtering - ResearchGate

この手法では、周波数／時間方向の1次元メディアンフィルタを元の楽曲信号にフィルタリングすることで調波／打楽器音を算出します。LibROSAの実装もこちらの手法がベースとなっています。

LibROSAにおける調波打楽器分離（HPSS）の実装

使い方

# Extract harmonic and percussive components
y, sr = librosa.load(librosa.util.example_audio_file())
y_harmonic, y_percussive = librosa.effects.hpss(y)
# Get a more isolated percussive component by widening its margin
y_harmonic, y_percussive = librosa.effects.hpss(y, margin=(1.0,5.0))

librosa.effects.hpss()

実装の詳細

LibROSAのHPSS（調波打楽器音分離）の実装を見ていきます。

librosa.effects.hpss()

y_harmonicy_percussiveがそれぞれ調波／打楽器音です。これら信号を様々な分析に利用したり、librosa.output.write_wav()で書き出したりします。

HPSSの本体はlibrosa.decompose.hpss() であり、それに時間周波数領域の楽曲信号をstft渡しています。

def hpss(y, **kwargs):

# Compute the STFT matrix
stft = core.stft(y)

# Decompose into harmonic and percussives
stft_harm, stft_perc = decompose.hpss(stft, **kwargs)

# Invert the STFTs.  Adjust length to match the input.
y_harm = util.fix_length(core.istft(stft_harm, dtype=y.dtype), len(y))
y_perc = util.fix_length(core.istft(stft_perc, dtype=y.dtype), len(y))

return y_harm, y_perc

librosa.effects.hpss()

librosa.decompose.hpss()

1. 時間周波数領域の楽曲信号Sを、振幅Sと位相phaseに分離
2. kernel_sizeより、調波／打楽器音成分用のメディアンフィルタの長さwin_harm/win_percを設定
3. marginより、調波／打楽器音成分用の閾値を設定（後述）
• margin_harm/margin_perc > 1の場合、その値に応じて閾値処理
6. 返り値
def hpss(S, kernel_size=31, power=2.0, mask=False, margin=1.0):
# (1)
if np.iscomplexobj(S):
S, phase = core.magphase(S)
else:
phase = 1

# (2)
if np.isscalar(kernel_size):
win_harm = kernel_size
win_perc = kernel_size
else:
win_harm = kernel_size
win_perc = kernel_size

# (3)
if np.isscalar(margin):
margin_harm = margin
margin_perc = margin
else:
margin_harm = margin
margin_perc = margin

# (4) Compute median filters. Pre-allocation here preserves memory layout.
harm = np.empty_like(S)
harm[:] = median_filter(S, size=(1, win_harm), mode='reflect')

perc = np.empty_like(S)
perc[:] = median_filter(S, size=(win_perc, 1), mode='reflect')

...

# (5)
power=power,
split_zeros=split_zeros)

power=power,
split_zeros=split_zeros)
# (6)

return ((S * mask_harm) * phase, (S * mask_perc) * phase)

librosa.decompose

LibROSAの実装では、メディアンフィルタによって算出した調波／打楽器音の振幅スペクトログラムをそのまま使うのではなく、ソフトマスク（ウィーナーフィルタ）を利用することで聴感的な歪みを軽減しています。

Extending Harmonic-Percussive Separation of Audio Signals - ISMIR 2014

LibROSA: 調波打楽器分離（HPSS）の適用例

コード

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

# Main
sr = 16000
margin = 1.5
offset = 116
duration = 20
filepath = "02 Retrograde Amnesia.mp3"

y, sr = librosa.load(filepath, sr=sr, offset=offset, duration=duration, mono=True)
y_harm, y_perc = librosa.effects.hpss(y, margin=margin)

# Dump signals
librosa.output.write_wav('hpss_org.wav', y, sr)
librosa.output.write_wav('hpss_harm.wav', y_harm, sr)
librosa.output.write_wav('hpss_perc.wav', y_perc, sr)

# Plot
plt.subplot(3,1,1)
plt.plot(y)
plt.title("Original signal")

plt.subplot(3,1,2)
plt.plot(y_harm)
plt.title("Harmonic signal")

plt.subplot(3,1,3)
plt.plot(y_perc)
plt.title("Percussive signal")

plt.tight_layout()
plt.show()

実際の楽曲への適用

• 調波楽器（エレキギター、ボーカル）のアタック音が打楽器音の方に含まれてしまう
• （特に）打楽器音の音質が悪い

が挙げられます。

コード進行やリズムの分析などの用途としては十分かもしれませんが、調波／打楽器音単体を聞きたいというような用途に対してはまだ余地があるように思います。そのような場合、パラメタ（例：フィルタサイズ、マージン）をチューニングしてみるか、最適化ベースの手法を検討してみるのもよいかもしれません。

また、調波音を取り出すメディアンフィルタが時間方向に大きな幅を持つことが、リアルタイム処理の場合は問題となります。従って、フィルタの形状等に何らかの工夫が必要になります。

まとめ

*1:HPS: Hermonic/Percussive Seperationと呼ばれることもあります

*2:時間周波数領域における楽曲信号の振幅／パワースペクトログラムを入力として与えます

はじめに

この記事では、Pythonの音楽分析モジュールである LibROSA

BPM自動算出の概要・設計方針については、以下の記事をご参考ください。

LibROSAのBPM自動算出の詳細

LibROSA によるBPM算出

# Estimate a static tempo
y, sr = librosa.load(librosa.util.example_audio_file())
onset_env = librosa.onset.onset_strength(y, sr=sr)
tempo = librosa.beat.tempo(onset_envelope=onset_env, sr=sr)

>>> tempo
array([129.199])

librosa.beat.tempo Examplesより

2. librosa.onset.onset_strength関数により、onset strength envelope（発音強度包絡）onset_envを得る
3. librosa.beat.tempo()onset_envを与え、BPM tempoが算出される

librosa.onset.onset_strength() について

（SuperFluxでは、さらに Peak-picking と呼ばれる識別処理により、Onsetの時間フレームか否かの2クラスに分類します。）

Maximum filter vibrato suppression for onset detection - DAFx-13

コアとなる処理部のコードを見てみましょう。
※デフォルトでは、ref = Slag = 1です。

# Compute difference to the reference, spaced by lag
onset_env = S[:, lag:] - ref[:, :-lag]

# Discard negatives (decreasing amplitude)
onset_env = np.maximum(0.0, onset_env)
...

if aggregate is None:
aggregate = np.mean
...

if aggregate:
onset_env = util.sync(onset_env, channels,
aggregate=aggregate,

librosa.beat.onset_strength()

これは音楽分析でよく使われる、 Spectral flux と呼ばれる特徴量です。ただし、Spectral fluxの負の成分は0としています。

LibROSAのBPM検出では、このSpectral flux を元にBPMを算出します。この周辺技術については、以下のドキュメントが参考になります。

librosa.beat.tempo()について

librosa.feature.tempogram で算出されるテンポグラムを元に（グローバルな）BPMが算出されます。

テンポグラム行列tgの形状は(win_length, len(onset_envelop))であり、1次元目がBPMに対応する周波数（以降、BPM周波数）となっています。テンポグラムは自己相関関数により算出された局所的なテンポ変化を捉えるためのリズム特徴量であり、端的に言えば、各時間フレームでどのBPM周波数っぽいかを表しています。

（テンポグラムの詳細については、別途まとめたいと思います。）

1. onset_envelopeより、tempogram関数よりテンポグラムtgを算出
2. tgを時間方向に集約（平均）する
3. tgBPM周波数のインデックスに対応するBPMが格納されたベクトルbpmsを取得
4. 事前確立（重み付け）ベクトルの負の対数であるlogpriorを作成
5. logpriorを考慮して、最も楽曲のBPMらしい、BPM周波数インデックスbest_periodを算出
6. best_periodbpmsインデックスとして与え、推定されたBPMを返り値とする
...
# (1)
tg = tempogram(y=y, sr=sr,
onset_envelope=onset_envelope,
hop_length=hop_length,
win_length=win_length)

# (2) Eventually, we want this to work for time-varying tempo

if aggregate is not None:
tg = aggregate(tg, axis=1, keepdims=True)

# (3) Get the BPM values for each bin, skipping the 0-lag bin
bpms = core.tempo_frequencies(tg.shape, hop_length=hop_length, sr=sr)

# (4) Weight the autocorrelation by a log-normal distribution
if prior is None:
logprior = -0.5 * ((np.log2(bpms) - np.log2(start_bpm)) / std_bpm)**2
else:
logprior = prior.logpdf(bpms)

...

# (5) Get the maximum, weighted by the prior
# Using log1p here for numerical stability ...(5)
best_period = np.argmax(np.log1p(1e6 * tg) + logprior[:, np.newaxis], axis=0)

# (6)
return bpms[best_period]

librosa.beat

最後に

LibROSAのBPM算出は処理が少し複雑ですが、BPM自動算出の設計方針にも書いたように、拍検出がしやすい信号（Onset envelope）に対して周波数分析（テンポグラム）を行うことでBPMを推定しています。

• 周波数帯域への重み付け
• 邪魔な成分の除去／抑圧
• 時間フレームへの重み付け、取捨選択
• ノイジーな（BPM算出が難しそうな）区間を捨てる／信頼性を重みづける
• BPMの値を制限
• 最小／最大のBPMを限定する
• 解像度

補足

Onset-envelope, Tempogramのプロット

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

y, sr = librosa.load(librosa.util.example_audio_file(), offset=15, duration=15)
onset_env = librosa.onset.onset_strength(y, sr=sr)
#tempo = librosa.beat.tempo(onset_envelope=onset_env, sr=sr)

ac_size=8.0
hop_length=512
win_length = librosa.core.time_to_frames(ac_size, sr=sr, hop_length=hop_length).item()

tg = librosa.feature.tempogram(y=y, sr=sr, onset_envelope=onset_env,hop_length=hop_length, win_length=win_length)

plt.subplot(3,1,1)
plt.plot(y)
plt.title("Audio signal")
plt.subplot(3,1,2)
plt.plot(onset_env)
plt.title("Onset envelope: onset_env")
plt.subplot(3,1,3)
plt.imshow(tg, aspect="auto")
plt.title("Tempogram: tg")
plt.tight_layout()
plt.show()

はじめに

どこかで聞いたことがあるのではないでしょうか。

あまり楽器が出回っておらず、演奏者も比較的少ないせいか、 音源の種類も少ないようです。

無料音源

SamsaraHurdy Gurdy FREE

32-bit版と64-bit版がありますが、自分の環境（FL Studio）では32-bit版のみ動作を確認しました。

有料音源

Hurdy Gurdy - Sonokinetic

1オクターブのドローン弦と、 2オクターブのチャンター弦を備えています。
ハーディガーディの画像の中央部にある、鍵盤設定パネルでは、

• Flip:
• ドローン弦とチャンター弦のキーボード割り当てを左右反転します
• Always Trans.:
• OFFの場合、チャンター弦をレガート（次の音へと滑らかに繋がる）で演奏します
• ONの場合、次の打鍵の際の打鍵音が含まれます
• IR
• 反響・残響の制御
• 3種類の空間プリセット（CASTLE HALL, ROOM, INN）
• ２つのパラメータ：　Room size, DRY/WET
• Aftertouch Response
• アフタータッチ（打鍵後の鍵盤の押し込みによる音色変化）を制御できます
• Volumes

などの設定が可能であり、細かい音作りができます。 音色が高品質であり、ハーディガーディを使った本格的な民族音楽を作るのに

• €49,90
• 48 kHz, 24bit mono aif format
• 1.04GB uncompressed
• Kontakt or Kontakt player 5.7.1 以上

Hurdy Gurdy - Rhythmic Robot

デモ音源を聞くと、民族音楽らしい空気感が作りやすい音源のように感じました。

また、上記サイトにはハーディガーディに関する詳しい説明が記載されています。
ハーディガーディの音作りの参考になるかもしれません。

また、Rhythmic Robot Audioではハーディガーディにインスパイアされたシンセ音源 Crank があります。
とてもアナログな温かみが感じられる音色です。
お得なHurdy Gurdyとのバンドルもあるのでご興味がある方はチェックしてみてください。

Hurdy Gurdy - Talesweaver Orchestra

こちらも非常に高音質なリアルな音源です。

この音色で比較的安価なのも魅力的です。

ERA II Medieval Legends - Best Service

\$259 と少し値が張りますが、他の中世ヨーロッパ系の民族楽器音源を揃えられるのがメリットです。

また、同社のフリーのサンプラー"Engine"向けのライブラリであるため、
Kontaktが必要ないのも注目すべきところです。