WaveSpectraを用いた歪率の測定について

WaveSpectra（以下WS）、および姉妹ソフトWaveGene（以下WG）や外部発振器を組み合わせて歪率（THD、THD+N）を測定する場合の制限や注意点などについてまとめてみました。

（お急ぎの方は説明を飛ばしてまとめを参照してください）

WSで、ごく普通に歪率を測定すると次のようになることがほとんどではないかと思います。
(1)THDは本来の値より小さく測定される。
(2)THD+Nでは本来の値より大きく測定される。

これらの測定誤差の原因としては２つ考えられます。
窓関数による影響とFFTの際の周波数分解能です。

WGで作成した信号で調べてみることにします。
※WGで作成したWaveファイルを入力信号としてWSで直接計算するので、たとえその時サウンドデバイスで再生しながらの場合でもサウンドデバイス自体の性能は無関係です。
（サウンドデバイスへ外部信号を入力して測定する場合とは違います）

(1) THD の場合

WGで一般的な48000Hz,16bitで次のような２つのWaveファイルを作成します。（Waveファイルの長さは適当）
1つ目は、
　Wave1: 1000Hzサイン波、-3dB
　Wave2: 2000Hzサイン波、-63dB
をミックスしたもの。
および、2つ目は、
　Wave1: 1000Hzサイン波、-3dB
　Wave3: 3000Hzサイン波、-63dB
をミックスしたもの。
２つはそれぞれ 2次高調波、3次高調波のみを持つ歪率 0.1% の1kHzの信号になっているはずです。

作成した信号をWSで再生して歪率（THD）を見てみます。
※WSは測定モードにし、THD、RMS のボタンを On にしておきます。
　（また、WSのFFTサンプルデータ数は設定ダイアログで 4096 にしておきます）

THD は、THD,+N の表示部の上側の数値です。

表１

窓関数	1kHz + 2kHz	1kHz + 3kHz
Hanning
Blackman -Harris
Flat top
無し（矩形）

窓関数によって差があることがわかります。　（特に窓関数無し（矩形窓）では当然ですが酷いことになっています）
さらに、1kHz + 2kHz と 1kHz + 3kHz でも差があることがわかります。

WSのTHDの測定は、ごく普通に、基本波の成分（特に指定しない場合は最大の成分）と、その周波数の整数倍の位置の成分を集めてRMS値をとったものの比を表示しています。

窓関数の周波数応答の影響で本来とは違う多くの歪成分を計算してしまうことになる矩形窓の場合を除いて、1kHz + 2kHz ではいずれも 0.08% 前後と、本来の値 0.1% より少ない値になっています。
1kHz + 3kHz ではさらに少ない値になっていますが、こちらも勿論本来の値は 0.1% です。

ただし、Flat top窓では本来の値に近い結果が得られています。

これらは窓関数の影響に加えて、もっと根本的な次のような理由によります。

各図の最大値MaxのdB値が（Flat top窓を除いて）本来の -3dB ではなく、またさらに窓関数によって微妙に違うことと、周波数が 1000Hz ではなく 996.1Hz になっていることから想像されるように、FFTの際の周波数分解能が影響しています。

上の場合、48000Hzサンプリング、4096点FFTなので、周波数分解能は 48000 / 4096 = 11.71875Hz となります。
よって、1000Hzの信号は直流を０番目として数えた場合、85番目の 11.71875 * 85 = 996.09375Hz の周波数成分が最大値となります。
（次の86番目の 1007.8125Hz も 1000Hz に近いですが、996.09375Hz の方が 1000Hz により近いので、85番目が最大値になります）

つまり、本来の 1000Hz から約4Hz 離れた低い部分の周波数成分を表示していることになる訳で、周波数成分の値は分解能分の幅のバンドパスフィルタを通した量と考えると、上の各図のMaxのdB値が窓関数によって違うのも分かるはずです。
窓関数の中心付近の周波数応答の形の違いによる（負の）誤差が出ています。

さらに、高調波成分については、996.09375Hz の整数倍の位置の成分を集めていくことによって、さらに（負の）誤差が増加して行くのがわかるはずです。

例えば、
2倍の 2000Hz では、85 * 2番目の 11.71875 * 170 = 1992.1875Hz （約8Hzのズレ）
3倍の 3000Hz では、85 * 3番目の 11.71875 * 255 = 2988.28125Hz （約12Hzのズレ）
等々となり、2000Hz、3000Hz とのズレは当然ながら高次高調波ほどどんどん拡大して行きます。

つまり、本来の周波数成分のピークからズレて低下した値が集められることになるため、THDは本来の値より小さく測定されることになります。
（高次高調波まで有る場合ほど本来の値より低く測定されるのはこのような理由からです）

ちなみに上の場合、2倍の成分は、実際には 85 * 2 の 170番目ではなく 171番目 2003.90625Hz で最大になりますし、3倍の成分は 85 * 3 の 255番目ではなく 256番目のちょうど 3000Hz で最大になります。
つまり最大になる中心が基本波の整数倍からズレてしまっています。

とりあえず中心付近の周波数応答がフラットな Flat top窓を用いると低下が軽減されますが、正確な値を得るにはWGと組み合わせて測定する場合には次のようにするのが最も簡単です。

以上から、逆に言うと、測定周波数がFFTの分解能の整数倍に一致している（中心にある）場合には正確に測定できる、ということになります。
つまり、最初から分解能の整数倍に一致した周波数で測定すればよい、ということになります。
言い換えると、FFTのサンプル長が測定信号の周期のちょうど整数倍に等しくなっている場合です。
この場合、切り出した部分の最初と最後で波形が繋がり、どこで切り出してもそのための誤差が無くなります。
（窓関数無し（矩形窓）を使用すると最も良い結果が得られるケースです）

WGで次のようにしてWaveファイルを作成します。
WGにはこういう時のために、周波数を簡単に FFTの分解能の整数倍になるようにする機能、"FFT用に最適化" 機能があるので、これを使います。（詳細は WGのヘルプを参照してください）
1つ目は、
　Wave1: 996.09375Hzサイン波、-3dB
　Wave2: 1992.1875Hzサイン波、-63dB
をミックスしたもの。
2つ目は、
　Wave1: 996.09375Hzサイン波、-3dB
　Wave3: 2988.28125Hzサイン波、-63dB
をミックスしたもの。
２つは最初と同じくそれぞれ 2次高調波、3次高調波のみを持つ歪率 0.1% の1kHzの信号になっているはずです。

※Wave1: 996.09375Hzサイン波は、まず Wave1 の周波数設定コンボボックスを 1000 と設定し、その上でマウスを右クリックし、"FFT用に最適化(F)" を選択すると、自動的に変換されます。
必ず先に、同じく右クリック "FFTサンプル数(T)" で 4096 を選択しておきます。（WSでのFFT長に一致させる）

※Wave2: 1992.1875Hzサイン波は、Wave2 の周波数設定コンボボックスの上でマウスを右クリックし、"Wave1の倍数(W)" で "x 2" を選択すると、自動的に変換されます。

※Wave3: 2988.28125Hzサイン波は、Wave3 の周波数設定コンボボックスの上でマウスを右クリックし、"Wave1の倍数(W)" で "x 3" を選択すると、自動的に変換されます。

これらは自動連動するわけではないので、Wave1 の周波数を変更した場合は、Wave2、3 の設定を再度繰り返す必要があります。

作成した信号をWSで再生して歪率を見てみます。

表２

窓関数	996.09375Hz + 1992.1875kHz	996.09375Hz + 2988.28125Hz
Hanning
Blackman -Harris
Flat top
無し（矩形）

今度は、正しい値になっていることがわかります。
（今回は "FFTに最適化" されていて、FFTのサンプル長が測定信号の周期のちょうど整数倍に等しくなっているためデータ切り出しの影響が無く、本来は "窓関数無し" の時に正しい測定が行われるようになります。　←各成分のスペクトルがたった1本だけになる）

基本周波数がもっと低い場合についても見てみます。
100Hzの場合について同様に調べたものです。

表３

窓関数	100Hz + 200Hz	100Hz + 300Hz
Hanning
Blackman -Harris
Flat top
無し（矩形）

1kHz の場合より窓関数による違いが大きいことが分かります。
Flat top窓では本来の値に近い結果が得られているようです。

"FFTに最適化" した場合についても見てみます。

表４

窓関数	100Hz + 200Hz	100Hz + 300Hz
Hanning
Blackman -Harris
Flat top
無し（矩形）

正しい値になっていることがわかります。

(2) THD+N の場合

THD+N の場合には窓関数の影響が大きいことがもっと簡単にわかるはずです。

まずは下図を見てください。

これは Hanning窓を掛けたときの 1kHz 0dB の信号のスペクトラム(これまでと同様4096点FFT)ですが、窓関数の特性による裾野が広がっていることがわかります。

これからTHD+Nを求める際に、基本波は最大値の位置の周波数成分を１つだけ採ることとしても、その他のノイズ分、つまり基本波以外の部分としてどこからどこまでを集めるか、ということが問題になることがわかるはずです。
裾野の部分を全部含めると、当然本来よりとても大きな値になることが想像されます。（実際そうです）
これはFFTサンプルの切り出し部分と窓関数による見かけの成分であり、信号自体には存在しない物だからです。

わかり易くするために上の図では THD+N の表示ではなく S/N の方で表示させているので、約96dBになっていれば良いわけです。
見ての通り、95dBほどになっており、ほぼ正しい値になっているので、逆に、あれっ、と思われるかもしれません。
しかし、これはほぼ正しくなるように基本波に近い裾野部分の一部の範囲（両側）をノイズ成分に含めないように除外しているからなのです。
（どの程度の範囲を除外するかのサジ加減は実験によって決定しました）

ためしに同じ信号をもっと裾野が広がるように 2048点FFT で表示した場合には次のようになり、約88dBとやはり増えていて、窓関数の影響が出ていることがわかります。（裾野がより広がったため除外しきれていない。2kHzにも広がっているのがわかる）

4096点以上のサンプル数にしておけば除外による影響の違いは余り無い様にしてあります。

ちなみに、1000Hz でなく、また "FFTに最適化" して 996.09375Hz とし、窓関数の裾野の影響を無くすと良い結果が出ると想像されますが、そうしたのが次の図です。
想像通り約96dBとなっています。（4096点FFT Hanning窓）

次に THD+N に戻って今度は窓関数を用いた場合の高調波成分（とノイズ）に注目すると、同様に高調波１つ１つについても窓関数による裾野の広がりがあるために、そのまま全てノイズ成分として集めると、上と同じ理由で本来の値より大きくなってしまいます。
高調波については基本波のように中心付近の除外はしていませんから全ての裾野の成分が加算されて大きな値となってしまいます。

(1)の表２において、THD より THD+N の方が大きくなっているのも同じ理由からです。
唯一、それぞれのスペクトルが一本だけで裾野が広がらない "窓関数無し" の場合だけ両者がほぼ同じ 0.1% と正しい値になっていることがわかるでしょう。

というわけで、THD+N の場合は本来の値より大きく測定されることになります。
（さらに、歪み成分の値が大きいほど、またそういう歪みの数が多いほど、裾野の影響による増加が多いこともわかるでしょう）

最後にサイン波＋ノイズの信号で確認しておきます。

WGで次のようにしてWaveファイルを作成します。
　Wave1: 1000Hzサイン波、-3dB
　Wave2: ホワイトノイズ、-55.1dB
これで、THD+N で測定した時にほぼ 0.1% の信号になっているはずです。

（下のように、２つの成分を別々にFFTしたものを見ると、RMS値の比 60dB、つまり THD+N が 0.1% となることがわかります）

1000Hz -3dB のみ　（RMS値 -6.01dB）	ホワイトノイズのみ　（RMS値 -66.01dB）

さらに比較のために、また 1000Hz を "FFTに最適化" した 996.09375Hz -3dB の信号も加えて THD+N を測定してみました。

THD+Nは、THD,+N の表示部の下側の数値です。

表５

窓関数	1kHz + ホワイトノイズ(-60dB)	996.09375Hz + ホワイトノイズ(-60dB)
Hanning
Blackman -Harris
Flat top
無し（矩形）

※測定周波数が "FFTに最適化" されていない一般的な場合には、やはりいずれも本来の THD+N の値 0.1% より大きくなっています。（Hanning窓が最も誤差が少ない？）

※測定周波数を "FFTに最適化" すると、本来の値より大きいことには違いはないですが、全てで改善され誤差が少なくなります。

なお、この場合には、本当は窓関数無しの場合に正しい測定ができるはずですが、約0.096% と本来の値 0.1% より逆に小さくなっています。
これは少し上で説明したように、基本波近傍の一部の範囲をノイズ成分から除外しているため、その分だけ逆にノイズ分が小さくなっているからです。
"FFTに最適化" して窓関数無しの場合には、もともと裾野が無く（スペクトルがただ一本）、除外する必要がないにもかかわらず除外したためです。
これについては将来、窓関数無しの場合には除外しないようにするかもしれません。
（窓関数無しの場合には、測定周波数を "FFTに最適化" した測定周波数だけが窓に同期していて、他のノイズ成分は同期していないため、矩形窓による切り出しのための誤差が出るはずです。　しかし、ノイズ成分があまり大きくない一般的な場合には、特に問題にならない程度だと考えられます）

(3) FFTサンプルデータ数について

FFTサンプルデータ数による違いについてまとめてみました。

これまでのように、基本周波数が1kHzの場合は以下のようにほとんど違いはありません。
窓関数は Hanning窓です。 THD 0.1% が正しい値です。

表６

サンプル数	1kHz + 2kHz	1kHz + 3kHz
4096
16384
65536

ただし、基本周波数がもっと低い場合は違って来るので注意が必要です。
下は 100Hz の場合について、同様に調べたものです。
サンプルデータ数を増やさないと高調波が分離できていないことがわかります。
ただ、増やしても、(1)の理由からTHDの値はそれほど正しくないこともわかります。（THD 0.1% が正しい値です）

表７

サンプル数	100Hz + 200Hz	100Hz + 300Hz
4096
16384
65536

"FFTに最適化" した場合は、4096点でも正しくなります。

表８

サンプル数	105.46875Hz + 210.9375Hz	105.46875Hz + 316.40625Hz
4096
16384
65536

まとめ

WSでの歪率測定の際に注意しておいたほうが良いと思われることをまとめました。

測定周波数の設定に注意する。

具体的には、なるべく外部発振器ではなくWGを信号源として用い、測定信号の周波数をWSのFFTサンプルデータ長に合わせてWG側で "FFTに最適化" するようにして、周波数がFFTの分解能の整数倍に一致するようにします。
このようにすると測定誤差が少なくなります。
用いる窓関数にも注意する。

THD のみの場合には 1. のように "FFTに最適化" した場合は、どの窓関数でもあまり変わらないはずですが、本来は窓関数無し（矩形窓）が良いはずです。
"FFTに最適化" しない場合（外部発振器を用いる場合等、最適化できない場合）には Flat top窓を用いるのが良いでしょう。
また、基本周波数が低い場合にも Flat top窓を用いるのが良いようです。

THD+N の場合は、同様に "FFTに最適化" しておいて、Hanning窓か窓関数無し（矩形窓）。
Hanning窓の場合は若干多い目、窓関数無し（矩形窓）の場合は現状では若干少ない目の値になることに注意。
FFTサンプルデータ数にも注意する。

"FFTに最適化" した場合には、あまり影響しませんが、"FFTに最適化" しない場合（外部発振器を用いる場合等、最適化できない場合）には、2次高調波が確実に分離できる程度のデータ長になるようにします。
無闇に多くしすぎると(1)の理由で誤差が大きくなります。
WGとWSを動作させるサウンドデバイスは、録音再生が同時に出来る（全２重動作が可能な）物を用いる。

再生/録音専用に別のPCを用いる場合など、別々のサウンドデバイスを用いると、２者間のクロックが完全に同期しないため、1.の "FFTに最適化" を行っても正確に測定できません。
（水晶発振器であってもこういう用途では僅かな差があってもダメです。
　やってみると分かりますが、スペクトルの幅がふらふらと変わります）
サウンドデバイス自体の歪特性を前もって自身のループ接続などで充分確認しておく。

当然ですが、測定精度は使用するサウンドデバイスの性能に依存します。
測定結果はデバイス自身の歪やノイズが加わったものになるわけですから、デバイスがどのような歪やノイズ特性なのか、事前によく把握しておいてください。

デバイスによっては、0dB近辺まで入力すると急激に歪み、-2～3dBまでの入力に制限しておく必要があるものがあるので注意してください。　（サウンドブラスター系カードの一部など）