FFTの2次補間に基づく正弦波パラメータ推定法の設計基準

社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS
信学技報
TECHNICAL REPORT OF IEICE.
FFT の 2 次補間に基づく正弦波パラメータ推定法の設計基準
∼ 擬似定常な正弦波成分の場合 ∼
安部
素嗣†
ジュリアス・スミス††
† ソニー株式会社
〒 141-0001 品川区北品川 6-7-35
†† スタンフォード大学 音楽学部 CCRMA
E-mail: †[email protected]
あらまし
オーディオ信号を構成する正弦波成分を推定する方法として，FFT スペクトルのピークを 2 次補間する方
法が広く用いられている．この方法は非常に簡便な上，実用上十分な精度が得られるのであるが，実際の精度はその
設計パラメータ，つまり窓関数形状，窓長，零詰め量などに依存する．本稿では，理論的，また数値的にこの正弦波推
定法の誤差を検証し，その設計基準を明らかにする．具体的には，許容される誤差の上限や信号成分の周波数間隔の
下限を与えた場合に，必要最小限の零詰め量や窓長などを，よく用いられる窓関数について明らかにする．また，ノ
イズを伴う信号に対する最尤性についても検討する．
キーワード
正弦波モデル，正弦波パラメータ推定，2 次補間，誤差解析，FFT
Design Criteria for Simple Sinusoidal Parameter Estimation based on
Quadratic Interpolation of FFT Magnitude Peaks
– For Quasi-Stationary Sinusoidal Components –
Mototsugu ABE† and Julius O. SMITH III††
† SONY Corporation
6-7-35 Kita-Shinagawa, Shinagawa-ku, Tokyo 141-0001
†† CCRMA, Music Department, Stanford University, CA 94305, U.S.A.
E-mail: †[email protected]
Abstract Due to its simplicity and accuracy, quadratic peak interpolation in a zero-padded Discrete Fourier
Transform (FFT) has been widely used for sinusoidal parameter estimation in audio applications. While general
criteria can guide the choice of window type, FFT length, and zero-padding factor, it is sometimes desirable in practice to know more precisely the requirements for achieving a prescribed error bound. In this paper, we theoretically
predict and numerically confirm the errors associated with various parameter choices, and provide precise criteria
for designing the estimator. Noise robustness of the estimator is also discussed.
Key words Sinusoidal Modeling, Sinusoidal Parameter Estimation, Quadratic Interpolation, Error Analysis, FFT
1. は じ め に
少なからず見られるものの — 理論的には正弦波パラメータの
最尤 (ML) 推定問題として扱われてきており，精度よく推定す
正弦波モデル [1], [2] は，オーディオ信号の直接的かつ強力な
る方法がいくつか提案されている（ [5], [6] 等）．中でも，FFT
記述モデルとして，分析，合成，変換，符号化等，様々な用途
スペクトルのピークの 2 次補間による推定法 [1] は，その簡便
（分
に広く用いられている（ [3], [4] 等）．その核となる部分に，
さと演算効率の良さに加え，大方の目的には十分な精度が得ら
析フレーム中の）信号から正弦波パラメータを推定する問題が
れることから，特によく用いられている．理論上も，よく分離
ある．この問題は— 中には漫然と直感̇的な手法を用いた例も
された正弦波成分に対して十分な零詰めを行った FFT を用い
—1—
る場合には，ほぼ ML 推定になるといってよい．
しかしながら，演算コストの観点からは零詰め量は必要最小
^
log A
True peak shape
log amplitude
限にしたく，また近傍成分からの干渉も常に無視できるとは限
らない．零詰め量の削減は，2 次補間による誤差バイアスとし
て現れる．近接成分の干渉もまたバイアスとなり，特に（窓長
と比して）ある程度以上に近接する成分は，信頼性をもっては
^
∆
Fitted parabola
推定され得ない．これらのバイアス量や周波数間隔の下限など
の正確な見積もりは，推定器の精度や効率を種々の目的に対し
kmax - 1
kmax
kmax +1
Freq. [bins]
て最適化する際の重要なポイントだと思われるが，従来は，や
図 1 スペクトルピークの 2 次補間
や直観的で漠とした十分条件に従った選択がなされてきている
Fig. 1 Quadratic interpolation of spectral peak.
ように思われる．
本稿の目的は，これらのバイアスと推定器の設計パラメータ，
すなわち窓形状，窓長および零詰め量との関係を，理論および
数値実験により考察し，この正弦波推定法の設計基準を明確化
することにある．まず第 3 節では，補間により生じる誤差バイ
アスを数値実験により検証する．これは真のピーク形状とあて
はめた 2 次関数との違いが誤差となるもので，時間領域での零
詰めにより FFT スペクトルの離散点が正確に補間されること
から，
（演算効率を度外視して）零詰め量を増やせばいくらでも
抑制できる．零詰め量の定量的な基準について，Brown [7] は，
零詰めのないハニング窓の周波数推定誤差が FFT グリッド間
隔の 5.3%以下であることを示し（注 1），Smith [1] は，矩形窓に
ついては零詰めを 5 倍以上に行うことで周波数誤差が 0.1%以
下（注 2）となると述べている．しかしながらこれらの基準も断片
的なものに過ぎず，必要十分な零詰め量を見極めるためには，
零詰め係数と誤差バイアスとの関係を，種々の窓関数について，
より詳細に明らかにする必要がある．
第 4 節では，近傍成分の干渉による誤差バイアスを検証する．
これは，概して周波数間隔が広いほど正確な推定になるのだが，
その関係は非常に非線形なものであり，特に周波数差がある程
度以下になるとバイアスは急激に増大する [8]．一般には周波数
間隔の下限を，窓関数のメインローブ幅（注 3） とする基準がよく
用いられるが，これは大まかな基準にすぎず，最適な窓形状や
必要十分な窓長を選択するという観点からは，より詳細な関係
を明らかにする必要がある．最後に第 5 節では，ノイズを伴う
信号に対するこの推定器の最尤性を検討する．
本稿では下記の表記を用いる．N : FFT サイズ，M : 窓サイ
「Quadratically Interpolated FFT 法（QIFFT 法）」は，下記
のようにまとめられる．
「適切
（ 1 ） オーディオ信号の振幅および位相スペクトルを，
に」零詰めされた FFT と「適切な」形状および「適切な」サ
イズの窓関数を用いて計算する．
（Fig. 1 の点）
（ 2 ） 振幅の最大点を見つける．(kmax )
（ 3 ） 近傍 3 点から対数振幅を 2 次補間する（点線）．
ˆ と
（kmax + ∆
（ 4 ） 振幅と周波数を補間関数から推定する．
log Â）
（ 5 ） 位相を位相スペクトルの 2 次補間から推定する．
（ 6 ） そのピークを FFT データから差引く．
（ 7 ） 残りの各ピークに対して 2–6 の処理を繰り返す．
本稿の目的は，
「適切な」窓形状，窓長，および零詰め量を与
える基準を明確化することである．なお，矩形窓で零詰め量が
1.5 倍以下の場合には，メインローブ内のサンプリング点数が
3 未満になることがあるゆえ，2 次補間は信頼性をもっては適
用され得ないことを指摘しておく．
3. 零詰め量の最小化
3. 1 数値実験のセッティング
入力信号には，離散時不変複素正弦波
x(n) = A0 ej(ω0 n+φ0 ) ,
を用いる．ただし，A0 , ω0 , φ0 は振幅，周波数，位相である．
QIFFT 法により {Â0 , ω̂0 , φ̂0 } を推定し，誤差を
ズ，Zp : 零詰め係数 (= N/M )．また，周波数の単位について，
Biasω = |ω̂0 − ω0 |/(2π/M ),
理論面では正規化された角周波数を用いるが，数値を表記する
BiasA = |Â0 − A0 |/A0 ,
際には，“窓サイズによって規格化された周波数” を用いる（注 4）．
2. QIFFT 法
スペクトルデータから正弦波パラメータを推定するための
（注 1）
：これは（対数をとらない）スペクトル振幅に対する 2 次補間を行った場
(1)
Biasφ = |φ̂0 − φ0 |/π,
(2)
により評価する．周波数バイアスを窓長の逆数周波数との比と
することで，実際の窓長によらない評価としている．
本実験では，{A0 , ω0 , φ0 } をランダムに変化させることで，
：窓長の逆数周波数に対する比率．
（注 2）
512 個 の 正 弦波 を 用 意す る ．さ ら に 8 通 りの FFT サ イズ
N = {64, 128, ..., 8192} を用い，窓長を N/Zp を超えず M >
= 31
：「（ハミング，ハニング窓の場合）窓内で 4 波以上の周波数差」と言われ
（注 3）
を満たす最大の奇数に設定する．零詰め係数 Zp を 1.0 から
合である．対数振幅に対して行うほうが一般には精度がよい [1]．
ることもある [9]．
（注 4）
：窓幅の逆数周波数を単位とする周波数．これにより実際の窓長やサンプリ
16.0 まで 0.1 ステップで変化させ，各零詰め係数に対して 4096
ングレートによらない数値となる．2π/M を乗ずれば正規化角周波数に，Fs /M
通りのテストセットの中から最大の誤差を取ることで，最大バ
を乗ずれば物理周波数に変換される．
イアス曲線を得る．
—2—
Bias [ratio to 1/M]
1
(a)
Bias [ratio to input]
W(ω - ω0 )
Rect
0.01
|X(ω)|
W(ω - ω1 )
Hann / Hamming
1e-04
Blackman
2
4
Zero padding factor
∆ω m
ω0
1e-06
1
8
16
1
Freq
ω1
(c)
∆ω r
Freq
(d)
Amplitude
Rect
Hamming
0.01
Freq
∆ω s
∆ω m /2
1e-04
Hann
1e-06
Blackman
2
4
Zero padding factor
Freq
図 3 近傍周波数成分の干渉：(a) 1 メインローブ幅（一般的な基準）
1e-08
1
8
16
(b) 2 ピークが丁度分離する周波数幅, (c) 近傍成分のサイドロー
ブの最初の零点までの幅, (d) 近傍成分のサイドローブの最初の
1
Bias [ratio to π]
ω peak
(b)
Frequency
Phase
1e-04
ピークまでの幅（本稿が提案する基準）,
Fig. 3 Interference from nearby components of various frequency
separations.
Rect
1e-08
表 3 （図 3 の）各条件に対する周波数差
1e-12
Hann / Hamming/Blackman
Table 3 Frequency separation values for the conditions in Fig. 3.
1e-16
1
2
4
Zero padding factor
8
16
Rc
図 2 零詰め係数に対する最大バイアス曲線
Fig. 2 Maximum-bias curves for various zero-padding factors
表1
Hn
Hm
Bl
∆ωm
2.00 4.00 4.00 6.00 従来の基準
∆ωr
1.37 2.00 1.84 2.35 ピーク分離条件
∆ωs
1.44 2.37 2.22 3.03 本稿の基準
所定バイアス以下に抑えるために必要な最小の零詰め係数．
Table 1 Minimum zero-padding factors for various frequency bias
bounds.
3. 3 零詰め係数の設計基準
誤差バイアスの上限を与えた際に必要な零詰め量の下限値
Frequency
Max Bias
Rc Hn Hm
Amplitude
Bl
Rc Hn Hm
を，表 1 にまとめた．例えば，周波数誤差 0.1%が許容可能で
Bl
1.00%
2.1 1.2
1.2 1.0 2.6 1.4
1.4 1.1
0.50%
2.5 1.5
1.5 1.1 3.0 1.6
1.7 1.2
0.10%
4.1 2.4
2.4 1.9 4.3 2.3
2.4 1.8
0.01%
8.7 5.0 5.1 4.0 7.5 4.0 4.2 3.1
Rc: Rectangular, Hn: Hann, Hm: Hamming, Bl: Blackman.
あれば，ハニング窓に関する零詰め量は 2.4 以上であればよい
ということである．
（ちなみに周波数誤差 0.1%は普通は聴き取
れない．
）
さらに，任意の誤差上限に対して必要な零詰め係数を得るに
は，最大バイアス曲線を簡単な関数で近似しておくと便利であ
ろう．これは指数関数でよく近似され，
表 2 最大バイアス曲線の近似係数．
Table 2 Coefficients for approximate maximum bias curves
Window
c0
ρ0
c1
ρ1
ρ0
Zp >
= c0 Biasω
(frequency),
(3)
c1 BiasAρ1
(amplitude),
(4)
Zp >
=
Rect
0.4467 −0.3218 0.8560 −0.2366
Hann
0.2436 −0.3288 0.4149 −0.2456
と表される．最小二乗法により得られた係数を表 2 に示す．例え
Hamming
0.2456 −0.3282 0.4381 −0.2451
ば，30ms のハミング窓対して 0.1Hz 以下に周波数バイアスを制
Blackman 0.1868 −0.3307 0.3156 −0.2475
限したければ，零詰め係数は 0.25 × (0.1[Hz] × 0.03[s]) −0.33 ≈
1.70 が必要である，となる．
3. 2 最大バイアス曲線
得られた最大バイアス曲線を図 2 に示す．まず，矩形窓で零
詰め係数が 1.5 以下の場合には，前述の理由により，例外的に
大きなバイアスが発生しているのが分かる．それ以外のケース
では，位相バイアスは，全ての場合において 10−12 % 以下に抑
えられている (図 2 下) ことがわかる．つまり位相は実用上は
バイアスを持たないといえる．これは対称な窓関数の線形位相
性により説明できる．また，周波数／振幅バイアスは，log-log
表記においてほぼ直線的に減少していることが分かる（図 2 上
中）．さらに，零詰め係数が同じ場合には，メインローブ幅の
大きな窓関数の方がバイアスが少ない傾向が確認される．
4. 周波数間隔の下限
4. 1 メインローブ幅による基準
複数正弦波成分の周波数間隔の下限について，一般には 1 メ
インローブ幅（零点から零点まで）の間隔があれば十分とされ
る（図 3(a)）．この場合，確かにメインローブ同士は重ならな
いため，誤差はサイドローブの影響のみに留まることが保証さ
れる．この周波数間隔を ∆ωm と表記する．代表的な窓関数に
関する具体的な値を表 3 に示す．
4. 2 詳細な基準
QIFFT 法は，パラメータ推定にピーク近傍のみを用いるの
で，上記の一般的な基準はやや十分すぎる条件と言える．より
—3—
表4
MAFS と最大誤差バイアス．
ら，QIFFT 法の周波数推定値は，振幅スペクトルの値のみで
Table 4 MAFS and maximum error bias.
はなく，その勾配の影響も受けるからであり，この周波数間隔
window
ではメインローブから継続する急な勾配が依然としてピーク推
Zp
Experiment
Theory
MAFS
Bias (%)
MAFS
Freq Amp
Rect
定に影響を与えるのである．
Pha
筆者らは，よりよい条件は，最初にサイドローブの勾配が 0
2.0
1.94
1.90
15.9
24.6 5.37
3.5
1.73
1.39
16.3
23.2 7.14
5.0
1.64
1.38
16.5
22.1 6.95
2.0
2.87
2.38
3.89
4.34 1.26
3.5
2.66
2.30
4.09
2.88 0.91
りとなるが，表 3 より，第二の条件が満たされれば，常にピー
5.0
2.57
2.28
4.15
2.74 0.87
クは分離していることが分かる．さらに，ここまで DTFT の
Hamm 2.0
2.72
2.35
1.37
2.19 0.60
連続周波数にて議論してきたが，実際には周波数軸が FFT に
3.5
2.51
2.20
1.42
0.87 0.27
より離散化されているため，最悪のサンプリング条件を考慮
5.0
2.42
2.18
1.45
0.76 0.24
2.0
3.53
3.23
0.39
0.28 0.07
3.5
3.32
3.05
0.39
0.13 0.04
5.0
3.23
3.00
0.40
0.13 0.04
Hann
Black
となる周波数間隔であると考える（図 3(d)）．これも窓関数固
有の値であり，各窓関数における数値を表 3 に示す．
これら二つの条件の大きい方が，周波数間隔の下限の見積も
しておく必要がある．ターゲット成分も妨害成分も周波数精度
は 0.5 FFT bin であるから，最終的に「最小許容周波数間隔
（Minimum Allowable Frequency Separation (MAFS) ）」の
周波数の単位は窓サイズによって規格化された周波数．バイアスの計
算は式 (2) を参照．
理論値は，
∆ωmin = ∆ωs +
詳細には，
2π
1
(2π/M ).
= ∆ωs +
N
Zp
(9)
と見積られる．零詰め係数が 2.0, 3.5, 5.0 のときの各窓関数に
（ 1 ） それぞれの成分が局所ピークをもち，
対する数値を，表 4 に示す．
（ 2 ） 他の成分のピークに対する干渉が小さい，
（必ずし
注意して欲しいことは，上記では MAFS の導出を，
ならばよいであろう．これらの条件を定量的に考察するため，2
も必要とは限らない）十分条件を用いて行ったことである．筆
つの正弦波からなる信号を考える．
者らは，MAFS の過少予測は重大な誤差を引き起こす可能性
x(n) = A0 ej(ω0 n+φ0 ) + A1 ej(ω1 n+φ1 ) ,
(5)
があるため，このやや広めの予測は，よりタイトで楽観的な予
測よりも実用的であると考えている．
ただし，{A0 , ω0 , φ0 } および {A1 , ω1 , φ1 } は，それぞれターゲッ
4. 3 数 値 実 験
ト成分と妨害成分の振幅，周波数，位相である．この信号の窓
4. 3. 1 セッティング
付き離散時間フーリエ変換 (Discrete-Time Fourier Transform
入力信号には，式 (5) に示した 2 つの複素正弦波から成る信
(DTFT)) [10] を X(ω) と表す．QIFFT 法では，最大の成分を
号を用いる．QIFFT 法により，ターゲット成分の周波数，振
最初に推定し，それを差し引いてから他の成分を推定するため，
幅，位相 (A0 , ω0 , φ0 ) を求め，式 (2) により評価する．
妨害成分からの最悪の干渉は振幅高がほぼ同じ場合であるから，
A0 = A1 の場合を考えればよい．
A0 = A1 = 1.0，そして，φ0 と φ1 は，[−π, π] の間の一様分布
第一の条件は，2 成分の中間の振幅スペクトルがピークのそ
乱数として，1024 本の信号を用意し，さらに 6 通りの FFT サ
イズ N = {256, ..., 8192} を用いる．周波数間隔 ∆ω を 0.025
れよりも低ければよいから（図 3(b)），
|X(ω)|ω=peak > |X(ω)|ω=(ω0 +ω1 )/2
本実験では，ω0 は [0, π] の間の一様分布乱数，ω1 = ω0 +∆ω ，
(6)
から 10.0 [rad/(2π/M )] まで 0.025 [rad/(2π/M )] ステップで
変化させ，各 ∆ω について 6144 通りのテストセットから最大
となるが，正確なピーク位置を知ることは難しいため，ピー
の誤差を取ることで，最大バイアス曲線を得る．零詰め係数は，
ク位置をターゲット成分の真の周波数（十分分離されたときの
2.0, 3.5 および 5.0 をテストした．
4. 3. 2 最大バイアス曲線
ピーク位置）に置き換え，
|X(ω)|ω=ω0 > |X(ω)|ω=(ω0 +ω1 )/2
(7)
零詰め係数が 5.0 のときの最大バイアス曲線を，図 4 に示す．
最大バイアス曲線には，多くの極大極小が見られるが，これは，
であれば十分である．式 (5) を用い最悪位相を考慮することで，
妨害成分のサイドローブ構造に起因するものである．周波数バ
ピークが分離するための最小周波数間隔の十分条件として
イアスの極大極小の位置が，振幅，位相バイアスのそれとは異
∆ωr = ∆ω s.t. W (0) − 2W (∆ω/2) − |W (∆ω)| = 0 (8)
を得る．これは各窓関数固有の値であり，数値解を表 3 に示す．
なることに注意して欲しい．これは，振幅，位相が主としてス
ペクトル値そのものから推定されるのに対し，周波数推定値は
スペクトルの勾配に大きく依存することによるものである．
第二の条件については，直観的には，2 成分の周波数間隔が
一般には，正確な周波数間隔は事前には分からないため，図
メインローブの半分でも，ピークに対しては互いのメインロー
4 の最大バイアス曲線の極小値を排し，各周波数間隔近傍での
ブの影響がないゆえ，十分なように思われる（図 3(c)）．しか
最大値を評価するのが有用であろう．これは，各周波数間隔に
しながら実際には，これは過小な予測となるのである．なぜな
おいて，その右側から最大値を取るようなフィルタリング演算
—4—
Hamming
Frequency
Rect
0.01
0.001
1e-04
Hann
Blackman
1e-05
0
2
4
6
8
10
Error [ratio to 2π/M]
Bias [ratio to 2π/M]
(a)
1
0.1
1
0.01
1e-04
Rect
Amplitude
Hamming
Blackman
0.001
1e-04
Hann
0
2
4
Blackman
6
8
10
Frequency Separation [window-normalized freq.]
Bias [ratio to π]
1
0.1
Rect
図 6
60
80
100
Frequency
100
1
0.01
1e-04
1e-06
Rect
Hann
Hamming
Blackman
-20
Audible limit
CRB
0
20
40
SNR [dB]
60
80
100
AWGN を 伴 う 信 号 に 対 す る 周 波 数 推 定 値 の RMS 誤 差
(M = 819, Zp = 5.0): (a) 単一正弦波信号， (b) 2 正弦波
0.001
信号．
Fig. 6 RMS errors in the frequency estimate for a signal of single
1e-04
Hann
2
4
Blackman
6
8
10
Frequency Separation [window-normalized freq.]
図4
周波数間隔に関する最大バイアス曲線 (Zp = 5.0)．
Fig. 4 Maximum-bias curves for various frequency separations.
Bias [ratio to 2π/M]
40
SNR [dB]
Hamming
0.01
0
and multiple sinusoids with AWGN.
の，依然として妨害成分のメインローブからの干渉が大きい．
急激な傾きは，メインローブ形状を反映したものである．
1
Frequency
0.1
さらに周波数差が広がると，ある点で，傾きがはっきりと緩
Rect
Hamming
0.01
やかになることが確認できる．これは，妨害成分のメインロー
Hann
0.001
ブ（傾き含む）の影響が，サイドローブの干渉量と等しくなる
Blackman
1e-04
点であり，それ以上の周波数差では誤差は十分小さいとみなし
1e-05
0
2
4
6
8
10
Frequency Separation [window-normalized freq.]
Bias [ratio to amp.]
20
Phase
1e-05
てよいであろう．徐々に減じる傾きは，窓関数のサイドローブ
高の減衰カーブを反映している．
（例えばハミング窓ではハニン
1
Amplitude
0.1
グ窓より減衰は緩やかである．
）この，傾きが変わる点を，
（実
Rect
Hamming
0.01
験的に得られた）MAFS とする．振幅や位相バイアスについて
Hann
0.001
も，概ね同様の現象が観測される．
1e-04
Blackman
他の零詰め係数に関する MAFS と，その周波数差における
1e-05
0
2
4
6
8
10
Frequency Separation [window-normalized freq.]
1
Bias [ratio to π]
0
(b)
0.01
1e-05
Hann/Hamming
Rect
Audible limit
-20
RMS Error [ratio to 2π/M]
Bias [ratio to amp.]
0.1
CRB
1e-06
Frequency Separation [window-normalized freq.]
1
Frequency, Zp=5.0
100
干零詰め係数の影響を受け，一般に零詰め係数が大きいほうが
Phase
0.1
MAFS は小さくなるが，零詰め係数が 3 を超えるあたりから
Rect
0.01
はその改善幅は非常に小さい．また，理論予測値が実験値を過
Hamming
0.001
Hann
小見積もりしていないことが確認できる．
1e-04
Blackman
1e-05
0
2
4
6
8
10
Frequency Separation [window-normalized freq.]
図5
最大バイアスを表 4 に示す．理論で述べた通りに，MAFS は若
フィルタされたの最大バイアス曲線 (Zp = 5.0)．
Fig. 5 Filtered maximum-bias curves.
4. 4 窓長に関する設計基準
ここまで，暗に窓長を固定し，信号の周波数間隔を中心に議
論してきたが，実際に制御し易いパラメータは窓長である．上
述の結果に基づくと，信号の周波数間隔が与えられたとき，必
要最小限の窓長を下記のように見積ることができる．
により得られる．図 5 にフィルタ結果を示す．この曲線は，横
M>
= 2π(MAFS)/∆ω,
(10)
軸値よりも大きな任意の周波数差に対する最大バイアスを示す
ものとなる．
MAFS は値は図 5，表 4，あるいは式 (9) により得られる.
周波数バイアス（図 5 上）を見ると，周波数差が小さいとき
例えば，零詰め係数 5.0 のハミング窓を用い，最小周波数間
は，全ての窓関数においてほぼ 100%を示している．この範囲
隔が 50Hz の場合には，窓長は 2.28/50 [Hz] = 45.6ms 以上で
では，両成分のメインローブが大きくオーバーラップするため，
ある必要がある．
ピークが分離されていないからである．周波数差がある程度を
5. QIFFT の最尤性
超えると，バイアス曲線は急激な減少を示す．この範囲では，
ターゲットのスペクトルは局所ピークとして分離しているもの
Rife [5], [8] は，1) 加算的ガウス性白色雑音 (AWGN) を伴う
—5—
単一正弦波の ML 推定は，矩形窓を用いたスペクトルピーク推
果は，文献 [13] および [15] として報告されている．
文
定と等価であることと，2) 他の窓関数を用いることで，複数成
献
分間の干渉は改善されるが推定誤差分散が上昇すること，を示
[1] J. O. Smith III and X. Serra: “PARSHL: A program for
した．QIFFT 法は，近似的にスペクトルピークを推定するも
the analysis/synthesis of inharmonic sounds based on a si-
のであるから，概ね ML 推定であることが期待される．多くの
nusoidal representation,” in Proc. ICMC’87, available at
応用では無限の精度は必要としないから，利用される範囲にお
ける ML 性が確保されれば実用上は十分であろう．
ここでは，2 次補間および近傍成分の干渉によるバイアスが，
ノイズ性正弦波の推定に与える影響を，数値実験により検証す
http://www-ccrma.stanford.edu/˜jos/parshl.
[2] R. J. McAulay and T. F. Quatieri:
“Speech Analy-
sis/Synthesis Based on a Sinusoidal Representation,” IEEE
Trans. Acoust., Speech, Signal Processing, Vol.34, No.4,
744/754 (1986).
る．QIFFT 法には 4096 点の FFT を用い，零詰め係数を 5.0
[3] K. N. Hamdy, M.Ali and A. H. Tewfik: “Low Bit Rate
に設定する．信号には，様々な S/N 比で AWGN を加えた単
High Quality Audio Coding with Combined Harmonic
一および複数の正弦波信号を用い，複数正弦波の間隔を各窓関
and Wavelet Representations,” Proc. IEEE ICASSP’96,
数における下限値（MAFS）とする．
紙面の都合で周波数誤差についての結果のみ，図 6 に示す
が，振幅，位相誤差も同様の傾向を示している．(a) は単一正
1045/1048 (1996).
[4] S. N. Levine and J. O. Smith III: “A Switched Parametric &
Transform Audio Coder,” Proc. IEEE ICASSP’99, 985/988
(1999).
弦波の，(b) は複数正弦波に対する結果である．合わせて一般
[5] D. C. Rife and R. R. Boorstyn: “Single-Tone Parameter Es-
的な可聴誤差限界（周波数誤差 0.1%を採用）と Cramer-Rao
timation from Discrete-Time Observations,” IEEE Trans.
Bound (CRB) も示した．まず，どちらのバイアスについても，
Info. Theory, 20, 5, 591/598 (1974).
概ね同様の傾向が現れていることがわかる．まず低 S/N 比の範
囲では，誤差は CRB よりも遥かに大きい．これは threshold
effect と呼ばれ [5]，ノイズ性のピークが正弦波成分のピークを
上回ることがあり，これを誤検出してしまうことに起因する．
中程度の S/N 比では，QIFFT 法は CRB に近づき，ほぼ ML
[6] D. J. Thomson: “Spectrum Estimation and Harmonic Analysis,” Proc. of the IEEE, 70, 9, 1055/1096 (1982).
[7] J. C. Brown et. al.: “A High Resolution Fundamental
Frequency Determination Based on Phase Changes of the
Fourier Transform,” J. Acoust. Soc. Am., 94, 2, 662/667
(1993).
推定器であることが確認される．厳密にいえば，矩形窓以外の
[8] D. C. Rife and R. R. Boorstyn: “Multiple Tone Parameter
場合には CRB を若干上回る誤差となるが，これは前述の通り
Estimation from Discrete-Time Observations,” Bell System
窓関数の副作用である．この範囲がオーディオ信号処理にて最
も重要な範囲であることに注意して欲しい．高 S/N 比の部分
では，誤差はほぼバイアスによって支配されている．前節まで
の議論は，このバイアス値を必要十分な精度に設計するための
指針を与えたものである．
6. お わ り に
特に正弦波周波数を推定する際，FFT スペクトルのピークの
2 次補間は，心得のある技術者ならまず当然行う工夫であろう
（補間すら行わずに FFT では精度が悪いというのは論外だろ
う）．しかしながら，それをただ漫然と行うだけでは，依然と
Technical Journal, 55, 9, 1389/1410 (1976).
[9] Ph. Depalle and T. Hélie: ”Extraction of Spectral Peak Parameters using a Short-Time Fourier Transform Modeling
and No Sidelobe Windows,” Proc. IEEE ASSP Workshop
on Applications of Signal Processing to Audio and Acoustics
(Mohonk’97), (1997).
[10] J. O. Smith III: Mathematics of the Discrete Fourier Transform (DFT), W3K Publishing (2003), web version available
at http://www-ccrma.stanford.edu/˜jos/mdft/.
[11] J. S. Marques and L. B. Almeida: “Frequency-Varying Sinusoidal Modeling of Speech,” IEEE Trans. Acoust., Speech,
Signal Processing, 37, 5, 763/765 (1989).
[12] M. Abe and J. O. Smith III: “Design Criteria for Simple
して最大 5.3%（∼半音）の誤差を覚悟しなければならない [7]．
Sinusoidal Parameter Estimation based on Quadratic Inter-
一方で，逆に精度が高いといわれる（が演算量も多い）方法を
polation of FFT Magnitude Peaks,” AES the 117th Con-
ただ盲目的に採用するのも，必要な精度は大抵有限であるゆえ，
vention, October (2004), to appear.
資源の浪費であることも多い．本稿は，簡便な 2 次補間法に対
[13] M. Abe and J. O. Smith III: “Design Criteria for the
して少々気を使う— 対数スペクトルを用い窓長や零詰め量を
Quadratically Interpolated FFT Method (I)(II)(III),” Tech-
丁寧に設計する — だけで十分高い精度が得られること，また
その精度限界に至るまではほぼ最尤推定器して働くこと，を示
したものである．
nical Reports, Dept. of Music, Stanford University, to appear.
[14] M. Abe and J. O. Smith III: “Correcting Bias in a Sinusoidal Parameter Estimator based on Quadratic Interpola-
本稿では割愛したが，筆者らは，Blackman-Harris 族以外の
tion of FFT Magnitude Peaks,” Technical Report, Dept. of
窓関数についても文献 [12] に報告している．また，本稿では定
Music, Stanford University, to appear.
[15] M. Abe and J. O. Smith III: “AM/FM Rate Estimation
and Bias Correction for Time-Varying Sinusoidal Modeling,” Technical Report, Dept. of Music, Stanford University,
to appear.
常正弦波を仮定した場合についての二つの主たるバイアス要因
について考察したが，もう一つ主要なバイアス要因として信号
の信号の時変性が挙げられる [11]．これは，一般に最大窓長に
関する限界を与えるものであるが，筆者らの現段階での考察結
—6—