音楽が好きな人であればShazamというアプリはご存知かと思います。
街中で気になった曲をスマートフォンに聞かせると、 相当マイナーな曲でない限り大抵の場合10秒程度聞かせると曲名やアーティスト名を教えてくれて、 しかも多少雑音が混じってしまっていても大丈夫という非常に便利なアプリです。
仕組みが気になったので、調べて一部をJuliaで実装してみました。
Juliaのバージョンは1.5.0です。
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.5.0 (2020-08-01)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
意外と単純な基本原理
最初、ディープラーニング的なテクニックを使っているのかと思ったのですが、調べてみると意外と手法自体はシンプルなものでした。
曲を認識させるにあたり、原曲の曲データの特徴量を以下の方法でフィンガープリント (Audio fingerprint/Acoustic fingerprint) します。
1. 曲データのスペクトログラムを作成する (スペクトログラムに関しては、前 にやりましたのでそちらを参照してください)
2. スペクトログラムのピークを探す。これは、適当な大きさの最大値フィルターをかけて元の画像と一致したらピークとします。
3. スペクトログラムのピークからペアを作る。と言っても全ての組み合わせのペアを考えるのではなく、二点が近いペアのみを考えます。
4. ペアからハッシュを生成する。ピークを(時間, 周波数)のペアで表した時に、前項で見つけたペアが $((t_1, f_1), (t_2, f_2)), t_1 < t_2$ であるとします。
この時、ハッシュテーブルに登録するキーと値はそれぞれ $[f_1:f_2:\Delta t]$ と $t_1$ になります。ここで、$\Delta t = t_2 - t_1$です。
という方法で曲から特徴量を抽出して曲のフィンガープリンティングを行います。
曲データにノイズが加わっていたとしても（全く同じではないですが、一致していることが確認するには十分に）似たようなハッシュデータが得られるようになっているのがポイントです。
マッチングの方法
Shazamを使ったことがある人ならわかると思うのですが、Shazamを使って曲名を認識させる際には、 曲を全部聞かせる必要はなく、途中から一部分を聞かせるだけで十分です。
ハッシュが Hash:time = $[f_1:f_2:\Delta t]:t_1$ であることがキーポイントで、 原曲と聞かせた曲の一部分では曲のスタート位置が当然違いますが、時間差は曲のスタート位置に依存しないので、ハッシュのキー $[f_1:f_2:\Delta t]$ は双方に出現します。 それぞれ対応する時間を$t_1$と$t^\prime_1$として、時間のオフセット$t_1 - t^\prime_1$を計算すると、 正しくマッチングできていれば原曲と聞かせた曲のスタート位置の差である定数になるはずです。
ハッシュが一致した時に、原曲と聞かせた曲で見つかった時間のプロットを確認してみましょう。
曲が一致している場合は上のようにマッチしたハッシュが傾き1の直線の上に並びます。 線から外れている点は、間違ってマッチングしてしまった点です。
逆に曲が一致していない場合であれば、間違ってハッシュがマッチングしてしまったとしても時間の差はランダムに分布することが予想されます。
あとで具体例を見せますが、時間のオフセットをヒストグラムにすると、マッチングしている場合は単一のピークが出現します。 理屈としては単純ですが鮮やかな手法で、最初に考えた人は頭いいですね。
更に詳しく知りたい方は原論文 A.WangのAn Industrial-Strength Audio Search Algorithm を読んでください。
Juliaで実装
本格的にShazamのようなマッチングを実現しようとすれば、 大量の楽曲データをfingerprintして、ハッシュデータをデータベースに格納して、 曲の認識の際はデータベースからの検索を行ってスコアリングし、最も類似度が高い曲を提示する、という処理が必要になると思います。
今回は、基本原理を調べて実装する段階で力尽きてしまったので、核となるfingerprintとマッチングできた場合・出来なかった場合のシグナルの可視化に絞って行いたいと思います。
コードは(Pythonですが)、worldveil/dejavuのfingerprint.pyを参考にしました。
WAVファイル
WAVファイルの読み込みは WAV.jl を使います。 詳細はここ を見てください。
julia> using WAV
julia> ys, fs, _, _ = wavread("test/data/original.wav")
([-0.44181646168401134 -0.40189825128940704; -0.43061616870632036 -0.3761101107821894; … ; 0.17578661458174383 -0.04297006134220405; 0.1151768547624134 -0.1497848445081942], 44100.0f0, 0x0010, WAVChunk[WAVChunk(Symbol("fmt "), UInt8[0x10, 0x00, 0x00, 0x00, 0x01, 0x00, 0x02, 0x00, 0x44, 0xac, 0x00, 0x00, 0x10, 0xb1, 0x02, 0x00, 0x04, 0x00, 0x10, 0x00])])
スペクトログラム
今回は簡単のため、スペクトログラムは DSP.jl を使って作りましょう。
DSP.Periodograms.spectrogramのドキュメンテーション はかなり簡素なので matplotlib.mlab.specgram や scipy.signal.spectrogram も参照したりしていました。
spectrogram が返す Spectrogram オブジェクトからスペクトログラムの2次元配列を取得するには .power でアクセスすれば良いです。
using DSP
function songspectrogram(samples, n, fs)
spec = spectrogram(samples, n; fs = fs, window = DSP.Windows.hanning).power
normspec = log10.(spec)
return normspec
end
julia> n = 4096
4096
julia> songspectrogram(ys[:, 1], n, fs)
2049×284 Array{Float64,2}:
-9.70616 -9.07873 -9.17428 -7.92614 -7.63062 -8.47055 … -9.40998 -8.13383 -7.47217 -9.79269 -9.69248 -9.65601
-8.32769 -8.22531 -6.85743 -6.56021 -6.97262 -7.67239 -7.66505 -6.78739 -6.7776 -7.62074 -7.43576 -8.01717
-8.03 -7.558 -5.53623 -5.9818 -6.17798 -7.68783 -6.94568 -5.43333 -5.61844 -6.74417 -8.42553 -7.95556
-5.43267 -5.78585 -4.39696 -4.4776 -5.94032 -5.69941 -6.9579 -4.4222 -4.48091 -7.12551 -6.49895 -6.33804
-4.77574 -4.66916 -3.62902 -4.67299 -4.92199 -4.97365 -5.47771 -3.91721 -3.9714 -5.75402 -5.86911 -6.84826
-4.4678 -4.45818 -3.10731 -3.27021 -4.21703 -4.90298 … -4.07139 -4.17006 -4.64172 -5.11589 -5.22632 -4.74766
-3.58357 -3.90762 -2.7872 -3.43402 -4.40083 -5.44588 -3.23751 -3.09109 -3.55516 -3.60761 -3.69366 -3.61307
-3.80939 -4.70587 -2.69206 -4.20379 -4.66742 -6.45696 -3.87312 -3.02429 -3.4179 -3.71929 -3.72232 -3.55507
-3.26219 -3.08835 -2.81918 -3.87016 -4.02507 -4.6202 -3.44302 -3.35965 -4.34106 -4.80432 -5.04157 -4.25322
⋮ ⋮ ⋱ ⋮
-12.6353 -12.9023 -10.7465 -12.6273 -15.0526 -11.7402 -11.5892 -12.8856 -11.8646 -11.332 -12.3234 -12.8574
-13.0681 -12.2647 -10.8669 -12.1085 -12.4588 -11.9597 -13.332 -13.1384 -11.9654 -12.3453 -11.8526 -14.436
-12.4502 -12.5558 -10.7107 -12.4894 -11.7158 -11.8924 -11.7527 -12.745 -12.8113 -12.007 -13.0835 -12.283
-11.849 -11.9355 -10.6198 -11.6698 -11.7177 -11.8629 -11.4374 -12.7903 -12.5631 -12.7018 -11.7313 -12.4594
-11.8144 -12.0432 -10.7412 -11.8539 -13.1211 -12.3507 … -11.5013 -11.6661 -12.4602 -12.5109 -11.4371 -12.5457
-11.9763 -11.1635 -10.687 -13.7323 -12.827 -11.7099 -11.7457 -12.372 -12.088 -12.1428 -11.8306 -12.9932
-11.5124 -11.445 -10.5719 -11.9994 -11.9061 -11.9427 -11.5007 -12.151 -13.2719 -11.2705 -12.2111 -12.879
-12.8636 -14.4831 -10.7474 -12.3807 -12.008 -13.2068 -11.5433 -12.7996 -13.0267 -11.3213 -12.6664 -13.6277
最大値フィルター
Scipyのscipy.ndimage.maximum_filter
に該当するフィルターが見つからなかったので、下のように最大値フィルターを構成します。
function maxfilter(matrix, filtersize)
temp, result = zero(matrix), zero(matrix)
n1, n2 = size(matrix)
for i in 1:n1
imin = max(i - filtersize, 1)
imax = min(i + filtersize, n1)
temp[i, :] = maximum(view(matrix, imin:imax, :), dims=1)
end
for j in 1:n2
jmin = max(j - filtersize, 1)
jmax = min(j + filtersize, n2)
result[:, j] = maximum(view(temp, :, jmin:jmax), dims=2)
end
return result
end
julia> x = [0.94 0.54 0.67 0.48 0.10;
0.77 0.98 0.54 0.63 0.32;
0.83 0.54 0.57 0.96 0.50;
0.56 0.40 0.79 0.69 0.43;
0.47 0.61 0.10 0.18 0.88]
5×5 Array{Float64,2}:
0.94 0.54 0.67 0.48 0.1
0.77 0.98 0.54 0.63 0.32
0.83 0.54 0.57 0.96 0.5
0.56 0.4 0.79 0.69 0.43
0.47 0.61 0.1 0.18 0.88
julia> maxfilter(x, 1)
5×5 Array{Float64,2}:
0.98 0.98 0.98 0.67 0.63
0.98 0.98 0.98 0.96 0.96
0.98 0.98 0.98 0.96 0.96
0.83 0.83 0.96 0.96 0.96
0.61 0.79 0.79 0.88 0.88
julia> maxfilter(x, 2)
5×5 Array{Float64,2}:
0.98 0.98 0.98 0.98 0.96
0.98 0.98 0.98 0.98 0.96
0.98 0.98 0.98 0.98 0.96
0.98 0.98 0.98 0.98 0.96
0.83 0.96 0.96 0.96 0.96
ピークを探す
次にmaxfilterをかけて得られた行列からピーク(実際には局所的なピークなわけですが)を取得します。
ピークを求めるときに背景部分は排除しています。
function getmaskindex(mask)
maskindex = getindex.(findall(mask), [1 2])
freqs = maskindex[:, 1]
times = maskindex[:, 2]
return collect(zip(times, freqs))
end
function findpeaks(matrix, filtersize)
maxmatrix = maxfilter(matrix, filtersize)
maxmask = maxmatrix .== matrix
nobackground = matrix .!= minimum(matrix)
return getmaskindex(maxmask .* nobackground)
end
julia> findpeaks(x, 1)
3-element Array{Tuple{Int64,Int64},1}:
(2, 2)
(4, 3)
(5, 5)
julia> y = zeros(5, 6); y[2, 2] = 1; y
5×6 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0
julia> findpeaks(y, 1)
1-element Array{Tuple{Int64,Int64},1}:
(2, 2)
dejavuではピークに出現させたくない箇所のマスクを作成する際に、背景を抽出した後 scipy.ndimage.binary_erosion を使って領域の内部を取得して作成、と言う手順を取っているのですが、このようにする意味はよくわかりませんでした。
余談ですがこのようなMathematical morphologyの操作は JuliaではImageMorphology.jlを使うとできます。 Scipyのscipy.ndimage.generate_binary_structureみたいに近傍の構造を指定する機能はなさそうです。
ハッシュを生成
次にピークの集合のペアリングを行った後、ハッシュを作成します。
全ての組み合わせてペアを作ると大変なので、各ピークに対してfanvalue個先までのピークをチェックし、
時間、周波数の位置の差分がそれぞれtimerange, freqrange(これらはPair{Int64}の元とします)に入っていたらペアリングすることにしましょう。
ハッシュは標準ライブラリSHAで生成出来ます。
using SHA
function paringpeaks(peaks, fanvalue, timerange, freqrange)
data = Vector{NTuple{4, Int64}}()
ntimes = Base.length(peaks)
# println(peaks)
mintdelta, maxtdelta = timerange
minfdelta, maxfdelta = freqrange
for (i1, (t1, f1)) in pairs(IndexLinear(), peaks)
for i in 1:fanvalue
i2 = i1 + i
i2 > ntimes && break
t2, f2 = peaks[i2]
dt = t2 - t1
df = f2 - f1
(mintdelta <= dt && dt <= maxtdelta && minfdelta <= df && df <= maxfdelta) || continue
push!(data, (f1, f2, dt, t1))
end
end
return data
end
function hashpeaks(peaks, fanvalue, timerange, freqrange)
hashdict = Dict{String, Int64}()
pairs = paringpeaks(peaks, fanvalue, timerange, freqrange)
for (f1, f2, dt, t1) in pairs
info = "$f1|$f2|$dt"
hash = bytes2hex(sha256(info))
hashdict[hash] = t1
# println("($t1, $f1) - ($t2, $f2), $info [$hash] -> $t1")
end
return hashdict
end
julia> peaks = [(1, 5), (2, 3), (2, 4), (3, 2), (4, 2), (4, 4), (6, 4)]
7-element Array{Tuple{Int64,Int64},1}:
(1, 5)
(2, 3)
(2, 4)
(3, 2)
(4, 2)
(4, 4)
(6, 4)
julia> fanvalue = 2
2
julia> timerange = 0 => 1
0 => 1
julia> freqrange = -3 => 3
-3 => 3
julia> hashdict = hashpeaks(peaks, fanvalue, timerange, freqrange)
Dict{String,Int64} with 8 entries:
"db724d0a500003163dce50a08d4cb5199d837df32ff9bea778229f6f89e0ec49" => 2
"9fc7745cc33e507d9ad28f16e9bd8d717b0de72ed078424da70292feb19248e4" => 1
"3533f2977e2cb6bb57e7135baff39dbb15e418fa6e7841216ebb6979110a5da4" => 2
"fe812d99d40bccea0f739ed5716d6f77af15b68fc73190b56bd11d579b7ed5d7" => 2
"798fb293c5fce0ad4b6c4405d361322ada948757accfbfb43c79b094702c419f" => 1
"8b0bbed14dafb6086bf675ee4fe2ab1e3a33a79bcc45f990918ba5cb24f12089" => 3
"1c22c9113f4f50934e03ac63bf365ddebfb220f6af556f6dd31ee9c8be4eb619" => 3
"c241ae9f10dd86f58a0c97363f4de1f08d10e3b0dae3cf139efd6027f2c75482" => 4
曲をフィンガープリント
今までの準備で部品は揃ったので、あとはつなげていくだけです。
今回は左右それぞれのチャンネルで独立にフィンガープリントして合体させています。
フィンガープリントの結果を可視化したくなったので、fingerprint_songに関連するコードが入っています。
長い時間の曲データのフィンガープリントの結果を可視化するのには時間がかかるのでご注意ください。
using Plots
function fingerprint(samples, n, fs, filtersize, fanvalue, timerange, freqrange)
spec = songspectrogram(samples, n, fs)
peaks = findpeaks(spec, filtersize)
hashdict = hashpeaks(peaks, fanvalue, timerange, freqrange)
return hashdict
end
function fingerprint_song(ys, fs, n, filtersize, fanvalue, timerange, freqrange; path = nothing)
hashdict = Dict{String, Int64}()
for i in 1:size(ys, 2)
samples = view(ys, :, i)
if !isnothing(path)
spec = songspectrogram(samples, n, fs)
peaks = findpeaks(spec, filtersize)
pairs = Hanauta.paringpeaks(peaks, fanvalue, timerange, freqrange)
heatmap(spec)
for (f1, f2, dt, t1) in pairs
plot!([t1, t1 + dt], [f1, f2], label="", linecolor=:blue)
end
scatter!(peaks, label="")
output = path * "_ch$i.png"
savefig(output)
end
newdict = fingerprint(samples, n, fs, filtersize, fanvalue, timerange, freqrange)
merge!(hashdict, newdict)
end
return hashdict
end
マッチングさせてみる
以上でWAVファイルからフィンガープリントを作成することができるようになったので、実際にマッチングできることを確かめてみます。 元のWAVにホワイトノイズを加えたデータや、スピーカーで曲を再生したものをiPhoneのボイスメモで録音してWAVファイルに変換したものでマッチングできることを確認します。
サンプル曲として、Windows 95のCD-ROMにビデオが収録されていたことでも知られる WeezerのBlue Albumに収録されている Buddy Hollyを使うことにします。 原曲データは本当はCDからWAVファイルとしてリッピングした方が良いと思うのですが、 再度リッピングするのが面倒だったので以前リッピングしたAACファイルをiTunesでWAV形式に変換して原曲データとしました。まあ大丈夫でしょう。
「聞かせた」時のデータですが、今回は2種類を試してみました。一つは原曲データにホワイトノイズを加えたもの(編集にはAudacityを使いました)と、 もう一つはより実践的なサンプルとして、原曲データをスピーカーで再生し、iPhoneのボイスメモ機能を使って録音し、 コンピューターに転送してMP4からWAVファイルに変換したものです。時間は13秒程度です。
ノイズを加える際は、Amplitudeが0.4と普通に聞いてかなり邪魔に感じるレベルのノイズを加えました。
以下のコードで原曲とオリジナル曲のハッシュのマッチング結果の図示を行います。
using Plots
using Measures
function plot_fingerprint(input, n, filtersize, fanvalue, timerange, freqrange)
ys, fs, _, _ = wavread(input)
# ysview = view(ys, 100000:300000, :)
ysview = ys
path, _ = splitext(input)
return fingerprint_song(ysview, fs, n, filtersize, fanvalue, timerange, freqrange; path = path)
end
function hashmatching(hash1, hash2)
i = 0
ts1 = Vector{Int64}()
ts2 = Vector{Int64}()
for h in keys(hash2)
if haskey(hash1, h)
t1, t2 = hash1[h], hash2[h]
push!(ts1, t1)
push!(ts2, t2)
i += 1
end
end
println("Match: $i")
if i > 0
scatter(ts1, ts2, label="", margin=5mm)
savefig(joinpath(@__DIR__, "output/scatter.png"))
tsdiff = ts1 .- ts2
histogram(tsdiff, label="", bins=50, margin=5mm)
savefig(joinpath(@__DIR__, "output/hist.png"))
end
end
n = 4096
filtersize = 10
fanvalue = 50
timerange = 1 => 20
freqrange = -200 => 200
input1 = joinpath(@__DIR__, "data/04 Buddy Holly.wav")
input2 = joinpath(@__DIR__, "data/noise_added.wav")
hash1 = plot_fingerprint(input1, n, filtersize, fanvalue, timerange, freqrange)
hash2 = plot_fingerprint(input2, n, filtersize, fanvalue, timerange, freqrange)
hashmatching(hash1, hash2)
まずはWAVファイルにノイズを加えたものの比較結果を見ましょう。横軸は原曲ファイルの時間、縦軸は比較対象ファイルの時間です。
点を画面に納める都合上アスペクト比率は1ではないです。
時間のオフセットのヒストグラムも見てみましょう。
ノイズの影響は受けつつもヒストグラムに単一のピークが認められました。
次にスピーカーから流した音をiPhoneで録音して試したところ、こちらはなぜかうまくいきませんでした。
何が原因なのか分からずかなり時間を使ってしまったのですが、ピークをプロットした点の画像を重ねてやっと原因がわかりました。
原因はサンプリング周波数の差で、音楽CDのサンプリング周波数は44.1kHzであるため、 CDからリッピングしたWAVファイルのサンプリング周波数も同じく44.1kHzになりますが、 iPhoneのボイスメモで録音したサウンドのサンプリング周波数は48kHzだったのが原因でした。
録音した音声ファイルを、Audacityでサンプリング周波数を変えてもう一度試してみます。
今度は上手く出来ました。
念のため全然関係ない曲だとマッチングできないことも確かめておきましょう。 本当は大量の曲を予め解析しておいて、実際に検索できることを確かめる必要があると思いますので手抜きです。
なお、曲は適当に15秒程度にトリミングしています。
- Blue Album一曲目のMy Name Is Jonas
誤ってマッチングしているハッシュが思いの外多いですが、ヒストグラムを見ると時間のオフセットはばらつきがあり特徴的な単一のピークはありません。
ハイパーパラメーターを調整して無駄なマッチ数を減らす、バックグラウンドを上手に除去するなど、 改良の余地はまだまだありそうですが、今回はここまでとします。
あまり綺麗ではないですが、コードは下に置いておきます。
https://github.com/matsueushi/AudioFingerprinting
参考
- Wang, A. (2003, October). An industrial strength audio search algorithm. In Ismir (Vol. 2003, pp. 7-13).
- Acoustic fingerprint (Wikipedia, The Free Encyclopedia.)
- The Basics of Audio Fingerprinting
- Audio Fingerprinting with Python and Numpy
- Shazamのしくみをちょっと理解してみる
- worldveil/dejavu
- JuliaDSP/DSP.jl
- 【Audacity】サンプリングレートを変換する方法