計測自動制御学会東北支部 第 295 回研究集会 (2015.6.26)
資料番号 295-5

LMS フーリエアナライザを用いた楽音の分析
Analyses of musical signals using an
LMS-based Fourier analyzer
○石橋 諭

工藤憲昌

○Satoru Ishibashi

Norimasa Kudoh

八戸高専 専攻科

八戸高専

田所嘉昭
Yoshiaki Tadokoro
豊橋技科大National Institute of Technology, Hachinohe College Advanced Engineering Course
National Institute of Technology, Hachinohe CollegeToyohashi University of Technology
キーワード：LMS アルゴリズム (LMS algorithm)，採譜 (musical transcription)
連絡先：〒039-1192

八戸市田面木上野平 16 番 1 号 八戸高専 電気情報工学科

tel:0178-27-7281， e-mail : kudohk-e@hachinohe-ct.ac.jp
ごとの正弦波の振幅を推定する．さらにその分析が正しいか、演奏のミスがあるかを判断するために，入力した信号の電力と推定した振幅の関係を導出し，入力信号電力と推定した振幅から求めた電力を比較することで分析結果を評価する．

1. はじめに日本の民俗芸能に「神楽」があり，神楽は代々弾き語りや口伝によって後世へと伝承されてきた．しかし，時代の変遷とともに後継者が減少しつつあり，やがて神楽が途絶えてしまう可能性があるという問題がある．この問題の対処法として，神楽を楽譜として残すことが有効であると考えられる．楽曲を楽譜に書き起こすことを採譜といい，採譜は音楽に精通した者でも膨大な時間を要するため，コンピュータによる自動採譜システムが研究されている．現在自動採譜ソフトはいくつか存在しているが，そのほとんどが西洋音階に従った採譜である．神楽は日本独特の音階を用いており，これを分析するには専用の自動採譜システムが必要である．
神楽では主に太鼓，鉦，笛が使われるが，本研究ではこれらのうち分析が困難とされる笛の演奏について分析を行った．分析を行うに当たって LMS (Least Mean Square：最小二乗平均) フーリエアナライザによって各音階の周波数成分

2. 原理
2.1 周波数推定
笛の各音階の周波数は振幅推定を行う際に既知である必要があり，その周波数が誤っている場合には正確な推定ができない．従って周波数適応フィルタを用いて厳密な周波数推定を行う．
Fig.1 に周波数適応フィルタのブロック図を示す．ここで𝐻𝑁 (𝑧)，𝐻𝑆 (𝑧)はそれぞれ特定の周波数を除去するノッチ特性，特定の周波数のみを通過させる BPF 特性を持つフィルタである．
入力信号を𝑥(𝑛)，ノッチフィルタの出力を𝑒(𝑛)，BPF の出力を𝑠(𝑛)，フィルタ特性の係数を𝛼̂(𝑛)とおくと，𝑥(𝑛)は(1)式のように表され，𝛼̂(𝑛)の更新式は(2)式のようになる．ここで𝜇は更新を-1-0反映する度合いを示すステップサイズパラメー
タで，𝛼̂(𝑛)の真値𝛼は𝛼 = 2 cos 𝜔であり，推定した周波数は(3)式のように表される．𝑓𝑠 はサンプリング周波数である．
𝑥(𝑛) = 𝑎(𝑛) cos 𝜔𝑛 + 𝑏(𝑛) sin 𝜔𝑛(1)(2)𝑎̂(𝑛)𝑓𝑓 = cos−1 ( 2 ) × 2𝜋𝑠(3)𝑥(𝑛)𝐻𝑁 (𝑧)𝐻𝑆 (𝑧)

amplitude[dB]𝑎̂(𝑛 + 1) = 𝑎̂(𝑛) − 𝜇𝑒(𝑛)𝑠(𝑛)𝑠(𝑛)Fig.1 frequency estimation filter

2.2 振幅推定

周波数推定の結果を用い，LMS フーリエアナライザによって振幅推定を行う．LMS アルゴリ
ズムは入力信号𝑥(𝑛)とその推定𝑥̂(𝑛)の誤差𝑒(𝑛)の二乗を最小化して推定値を真値に近づける方法であり，LMS 法はこのアルゴリズムを用いて対象の周波数における振幅を逐次推定するアルゴリズムである．入力信号𝑥(𝑛)が(4)式のように表されているとする．ここで 𝑎𝑖 (𝑛), 𝑏𝑖 (𝑛) は時変の振幅，𝜙(𝑛) は加法性白色雑音である．𝑥(𝑛)の 推 定 値 𝑥̂(𝑛) が (5) 式 の よ う に 表 さ れ ，𝑎𝑖 (𝑛), 𝑏𝑖 (𝑛)の推定値をそれぞれ𝑎̂𝑖 (𝑛), 𝑏̂𝑖 (𝑛)としたとき，誤差𝑒(𝑛)は(6)式のように表される．また，𝑎̂𝑖 (𝑛), 𝑏̂𝑖 (𝑛)の更新式はそれぞれ(7),(8)式のようになる．𝜔𝑖 は周波数推定で得られた既知の角周波数である．
(4)(5)

𝑒(𝑛) = 𝑥(𝑛) − 𝑥̂(𝑛)(6)𝑎̂𝑖 (𝑛 + 1) = 𝑎̂𝑖 (𝑛) + 𝜇𝑒(𝑛) cos 𝜔𝑖 𝑛(7)𝑏̂𝑖 (𝑛 + 1) = 𝑏̂𝑖 (𝑛) + 𝜇𝑒(𝑛) sin 𝜔𝑖 𝑛(8)

-40
0.1
0.2
0.3
0.4
normalized frequency

0.5

2.3 提案法
楽音はステップ関数的な立ち上がりやランプ関数的な減衰を含むことが多く，楽音の分析は追従性能を重視しなければならない．しかし実際の演奏は演奏ミスや周囲の雑音が存在し，これを正しく分析するためには十分な雑音耐性も必要である．十分な追従性能と雑音耐性を両立させるために，提案法を用いて振幅推定を行う．
提案法のブロック図を Fig.3 に示す．
提案法で−1 )−1は適応制御ループ内に(1 − 𝛾𝑧のリーク付き積分操作を導入した．ここで，𝛾 は帯域幅を制御するパラメータであり，ステップサイズパラメータ𝜇と𝛾の間には(9)式のような制約がある[1]．また，提案法における 𝑎̂𝑖 (𝑛), 𝑏̂𝑖 (𝑛)の更新式はそれぞれ(10),(11)式のようになる．

𝑧 −1𝑥̂(𝑛) = ∑𝑝𝑖=1 𝑎̂𝑖 (𝑛) cos 𝜔𝑖 𝑛 + 𝑏̂𝑖 (𝑛) sin 𝜔𝑖 𝑛-30

Fig.2 frequency characteristics of LMS method

𝑎̂(𝑛 + 1)

𝑝𝑥(𝑛) = ∑𝑖=1{𝑎𝑖 (𝑛) cos 𝜔𝑖 𝑛 + 𝑏𝑖 (𝑛) sin 𝜔𝑖 𝑛} + 𝜙(𝑛)μ=0.05μ=0.01-20-500𝑒(𝑛)−𝜇fs=8000[Hz]fi=2000[Hz]-10(1−𝛾)21+𝛾≥𝜇(9)𝑎̂(𝑛) = (1 + 𝛾)𝑎̂(𝑛 − 1) − 𝛾𝑎̂(𝑛 − 2) + 𝜇𝑒(𝑛) cos 𝜔𝑛(10)𝑏̂ (𝑛) = (1 + 𝛾)𝑏̂ (𝑛 − 1) − 𝛾𝑏̂ (𝑛 − 2) + 𝜇𝑒(𝑛) sin 𝜔𝑛(11)Fig.3 a block diagram of proposed method

LMS 法において，𝜇 の値が大きいと収束速度が向上し雑音耐性が悪くなり，𝜇 の値が小さいと雑音耐性は良くなるが収束速度が遅くなる．このことは Fig.2 に示す LMS 法の周波数特性から分かる．𝜇 = 0.05の場合は𝜇 = 0.01の場合よりも帯域幅が広がり，通過帯域から離れた周波数における振幅が大きい．これは対象周波数以外の信号を誤検出し，雑音に弱いということを示す．

提 案 法 の 定 常状 態 に おけ る 周 波 数特 性 をFig.4 に示す．𝛾 = 0.89のとき，𝛾 = 0.60よりもピーク付近が少し広くなっていることがわかる．
これは演奏のミスなどによって規定の周波数から少しずれた場合でも対応ができることを示している．また，Fig.2 の従来の LMS 法の周波数特性と比較すると，通過帯域から離れた領域でも十分な遮断が行えていることが分かる．これは提案法がより雑音に強いことを示している．

-2-

2.5 信号電力の比較
前述のような方法を用いて楽音の分析を行うが，分析結果が正しいかどうかを確かめる必要がある．そこで，入力信号電力と推定した振幅から求めた電力（以下，推定電力と呼ぶ）を比較し，分析結果の評価を行う．
入力信号を𝑥(𝑛)，角周波数𝜔𝑖 の cos 成分，sin成分の振幅推定値をそれぞれ𝑎̂𝑖 (𝑛), 𝑏̂𝑖 (𝑛)とすると，入力信号電力と振幅推定値の関係は(12)式で表される．ここで𝐸[∙]は平均操作である．(12)式の計算を行い，入力信号電力よりも推定電力が小さければ，入力信号の中に演奏のミスによる推定していない周波数が存在することになる．

0fs=8000[Hz]fi=2000[Hz]μ=0.01

amplitude[dB]-10-20γ=0.89γ=0.60-30-40-50
00.1

0.2
0.3
0.4
normalized frequency

0.5

Fig.4 frequency characteristics of proposed method

2.4 隣接音高の漏れ込み成分の除去
Fig.4 に示したように，提案法は BPF のような周波数特性を有しており，推定対象の音高と隣接する音高が存在する時，その隣接音高の影響を受けて誤検出につながる．これを漏れ込みと 呼び ，漏 れ込 みを 除去 する 必要 があ る．
Fig.5-(a) に従来法の，Fig.5-(b) に漏れ込みを除去する振幅推定のブロック図をそれぞれ示す．
𝐺𝑖 (𝑧) は音高 No. 𝑖 における LMS 法のアルゴリ
ズムの定常状態における伝達関数である．隣接音高の周波数は既知であるため，推定対象の音高の振幅を推定する際に隣接音高の振幅も同時に推定し，その振幅を破棄して推定対象の音高の振幅のみを採用することで漏れ込みを低減できる．Fig.6 に漏れ込みを除去する方法の周波数特性を示す．推定対象の音高の周波数付近のみを通過させ，隣接音高の周波数成分は除去できていることが分かる．

(a) conventional method

2𝐸[𝑥 2 (𝑛)] ≈ 𝐸 [∑𝑖{𝑎̂𝑖 (𝑛) cos 𝜔𝑖 𝑛 + 𝑏̂𝑖 (𝑛) sin 𝜔𝑖 𝑛} ]1= 2 ∑𝑖{𝑎̂𝑖2 (𝑛) + 𝑏̂𝑖2 (𝑛)}入力信号に直流成分が存在する場合，電力比較を行うと直流的な誤差が生じる．この現象の対策として，入力信号にHPFをかけてから電力計算を行う．HPFの伝達関数は(13)式で表される．
𝐻(𝑧) = 1 − 𝑧 −1
(13)

3. シミュレーション
3.1 神楽笛の各音階の周波数推定
周波数適応フィルタを用いた，神楽笛の各音階の周波数推定結果をTable.1に示す．神楽笛の音は倍音を含んでおり，原音より倍音の振幅が大きいと倍音の周波数に収束することがある．
従ってFig.1のフィルタをトリー状に接続し，4周 波推 定を 行っ た． サン プリ ング 周波 数は8[kHz]であり，各音階は音が低い順に音高No.1，No.2, … のように表記している．また，音高No.1
は笛の音階ではないが，実際の演奏の際に発せられる音の周波数である．
Table.1 frequencies of scales
推定周波数[Hz]音高No.1546.8750
No.2682.6386
No.3
846.3885
No.4
994.6145
No.5
1061.1350
No.6
1237.7460
No.7
1364.2070
No.8
1666.3798
No.9
1953.7710(b) method to remove interference

Fig.5 block diagrams of amplitude estimation0

amplitude[dB]-10μ=0.01, γ=0.8
fs=8000[Hz]-20

3.2 神楽笛の振幅推定
Table.1のデータを基に，実際の神楽笛の演奏の振幅推定を行った．入力形式はサンプリング周波数8[kHz]，モノラル，16[bit]のwavデータで，シミュレーション条件は𝜇 = 0.01, 𝛾 = 0.8である．振幅を推定する際，除去しきれなかった漏れ込み成分や倍音を取り除くために，各音高の中で一番振幅が大きいものを探し振幅が一番大きい音高以外の成分を0にする．漏れ込みや倍音-30-40-5000.1
0.2
0.3
0.4
normalized frequency

(12)

0.5

Fig.6 frequency characteristics ofinterference removal-3-

れば(12)式の関係は成り立たないため，平均するサンプル数について検討する必要がある．
40Power00150

100

00No.4No.300No.2
No.1

0

10000

20000

30000

μ=0.01, γ=0.8
fs=8000[Hz]true value
estimated value

true value
estimated value

50

50

No.5

200μ=0.01, γ=0.8
fs=8000[Hz]

150

Power

Amplitude

No.6

100
no. of blocks

Fig.9 signal power comparison 1100

No.7

2010

No.9
No.8

true value
estimated value

μ=0.01, γ=0.8
fs=8000[Hz]30Powerは原音より確実に振幅が小さいため，この方法による除去が可能である．また，それぞれの音高の振幅が視覚的にわかりやすいように（No.1の振幅）+1，（No.2の振幅）+2のようなバイアスを加えて表示している．Fig.7 は笛の音階を低い順に鳴らしたものであり，漏れ込み除去を行っているため正確な振幅推定ができていることが分かる．また，Fig.8に神楽笛の演奏を対象とした振幅推定結果を示す．はじめの音高No.5の前後においてNo.6の微小な検出があるが，これはNo.5の音が規定の周波数より少し高い周波数で演奏されたために生じた誤検出である．

20204060
no. of blocks
80804060
signal power comparison 2
no. of Fig.10
blocksno. of samples

4. まとめFig.7 amplitude estimation 1

本稿ではまず周波数適応フィルタを用いて神楽笛の各音階の周波数を推定した．次にその結果に基づき，LMS 法によって振幅推定を行った．
しかし従来の LMS 法では追従性能，雑音耐性が悪い．そこで適応制御ループ内にリーク付き積分操作を導入することによって追従性能，雑音耐性の改善を行い，これを数値実験より確認した．さらに隣接音高からの漏れ込みを除去することでより正確な振幅推定を実現した．また，
振幅推定の結果が正しいかを確認するために入力信号の電力と推定した振幅の関係を導き，これを用いて信号電力の比較の数値実験を行った．
今後は振幅推定の結果に基づいた MIDI ファイ
ルや楽譜の作成，リアルタイムな振幅推定，そして作成した楽譜が正しいかどうかを確認するシステムの構築を行うことが目標となる．

No.9

Amplitude

No.8
No.7
No.6
No.5
No.4
No.3
No.2
No.1

05000

10000

no. of samples

Fig.8 amplitude estimation 2
3.3 信号電力の比較
Fig.9 と Fig.10 に神楽笛の音階を低い方から順に演奏したものを対象としたものと，神楽笛の演奏を対象としたものの信号電力の比較結果をそれぞれ示す．シミュレーション条件は𝜇 = 0.01, 𝛾 = 0.8 であり，150 サンプルごとに平均電力を求め，これを 1block とした．Fig.9
において，50block 目までがあまり一致しておらず，それ以降はほとんど一致していることが分かる．これは一致していない部分において推定していない周波数成分が存在していることを示している．また，Fig.10 において 0～25 サンプ
ル付近で非常に大きく誤差が出ていることが分かる．これは演奏が規定の周波数に従っていないことが原因と考えられる．50block 付近においては推定値が真値を上回っているが，これは通常起こり得ないことである．考えられる原因としては，電力を計算する際の平均数が誤っていたと考えられる．各周波数の正弦波が一周期のほぼ整数倍の分だけ入るように平均を取らなけ

謝辞
本研究は JSPS 科研費 26350352 の助成を受けたものです．
参考文献[1]. N.Kudoh & Y.Tadokoro (2003) : ” PerformanceAnalysis of an LMS Based Fourier Analyzer for
Sinusoidal Signals with Time-varying Amplitude”
Proc. of IEEE TENCON’03[2]. 工藤憲昌・田所嘉昭(2013)，「周波数変動を許容した LMS フーリエアナライザの特性」 八戸工業高等専門学校紀要第 48 号

-4-