周波数の分解能は?
パソコンで、音の周波数を測定するとき、どこまで低い周波数を検出できるのかという疑問をもっていました。
パソコンで音を扱うのは歴史が古く、ほぼ終わったテーマのようです。しかし、音楽向けに周波数の観点から音を計測するプログラムと言うのは見たことがありません。声紋や音声認識の例はありますが、楽器の音合わせといったものは聞いたことがありません。
- 基準音に使われる440Hz（A、ラ）の音とA#(466.１6Hz)の間隔は,26.16Hzで、これぐらい違うと、私にも音の違いがわかります。
- 理科年表にあった音階の周波数の表は、小数点以下2桁まで書かれています。
- 人の可聴範囲は、20Hz-20KHz程度らしい。
たとえば、「小数点以下2桁まで」を考えて見ます。0.01Hzというのは、1波長が100秒（1分40秒）と言うことなので、測定時間が1分40秒以上必要な気がします。これは実用的ではなく、同じ状態が継続するのはせいぜい数秒程度ではないかと思います。もし、1波長が3秒なら1/3Hzとなります。
音声認識では、数十ms単位で周波数を分析しているようですが、40msとすると25Hz以下は使われていないことになるのでしょうか
こんな疑問を持っていました。PCMデータから周波数を知るには、FFTを使うものらしいので、FFTの算出結果の意味を考えてみたいと思います。
前提として、データは一定の周期でサンプリングされた、偶数個のサンプルとします。
ここで、作図に使っているのは、「R言語」「サンプル」のスクリプトです。
１．式から
tは時間。Nは、サンプル数。ωは、２π/N。
こんな感じの式のようです。
１．１．exp(ｘ・i)の性質 
t=0 から 2π/16 刻みで、
exp(0・i)、exp(2π/16・i)、．．．
を計算する。
結果のそれぞれの実部、虚部を、座標（実部、虚部）としてプロットすると、左図ができます。
exp(x・i)は、絶対値が1で、2πで一周する周期性があります。
また、
z(n)=exp(n・2π/16・i) nは、0,1,...,15)
とあらわすと、実軸に対象な点は、
z(-1)=z(15)
z(-2)=z(14)
：
z(-7)=z(9)
と、なっています。
たとえば、
z(-7)=exp(-14π/16・i)=exp(18π/16・i)
z(9)=exp(18π/16・i)
となります。
サンプル数/2で折り返す。n > N/2 の箇所は、N-n 番目の値の共役数になっています。
また、
exp(ix) = cos(x) + i・sin(x)
exp(-ix) = cos(x) - i・sin(x)
exp(ix)・exp(iy) = exp(ix・iy)
2cos(x) = exp(ix) + exp(-ix)
-2i・sin(x) = exp(ix) - exp(-ix)
で、あるらしい。
１．２．式の意味
分からないなりに考えて見ました。
この式は、元の信号を多項式で表現しようとするもので、サンプル数の項でできています。
exp(x・i)の性質から、各項はそれぞれ周波数を示しているのだと思います。どんな周波数を示しているかを考えて見ます。
nは、0,1,2,...,N-1です。n=0,1を除くと、各項は、順次前の項の2倍の速さで回転します。
n=1のときの周波数はいくつでしょうか？ ω=2π/NなのでN回で一周することが分かりますが、N回が何秒なのかは表されていないので、Hzで周波数を示すことはできません。サンプリングに要した時間（サンプリング時間 = サンプリング数 / サンプリングレート）を基準として、n=1,2,3,...,N-1は、周波数1,2,3,...,N-1を表しているのだと解釈しました。
たとえば、16サンプルのデータがあるとして、
サンプリングに要した時間が1秒（または、サンプリングレートが16Hz）なら、2π=1秒で、最長1Hzです。
サンプリングに要した時間が0.01秒（または、サンプリングレートが1600Hz）なら、2π=0.01秒で、最長100Hz です。
各項の係数a(n)は、振幅を表すことになります。サンプリングデータに該当する周波数が含まれないときは、a(n)はゼロになります。
１．３．FFTの計算
目的は、サンプリングされた値から含まれる周波数を知りたいと言うことです。各項の係数a(n)を計算するのがFFTの計算になります。
※振幅と書いたが、振幅から算出される値で、その周波数の強さを示しているのだと思う。また、周波数(N-1)/tと書いたが、サンプリングレートの半分の周波数以上は、本来検出不能なのだと思う。
３．プログラムから
シンプルなFFTのプログラムをC#で書いてみました。前述の式の各項の係数a(n)を求めるので、出力は、サンプル数と同じ個数の複素数になります。
2重のループからなり、外側のループは、前述の式の項に対応します。内側のループは、サンプルを列挙するものです。
ループ内では、全サンプルにexp(-n・ω・i)を掛けて積算します。結果の配列には、全サンプルにexp(-n・ω・i)を掛けた値の和が求められますが、
d[1]には、nを0,1,2,3,...
d[2]には、ｎを0,2,4,6,,...
のようにずらしながら計算したexp()の値がつかわれます。これは、exp()の周期性を利用して、1周期から順に強さを求めることになるのだと思います。
public class FFT
{
public static Complex[] basic3(double[] s)
{
Complex[] d = new Complex[s.Length];
double pip2 = -2.0 * Math.PI / s.Length;//-ixの定数部分を先に計算しておく
for (int i = 0; i < d.Length; i++)
{
d[i] = new Complex();//0 + 0i に初期化
for (int j = 0; j < s.Length; j++)
{
Complex ix = new Complex(0, pip2 * i * j);
d[i] += s[j] * Complex.Exp(ix);
}
}
return d;
}
}
３．１．値の大きさについて
同じ入力を同じ時間、サンプリングレートを変えてサンプリングした場合を考えます。たとえば、ｎHzと2nHzでサンプリングした場合、サンプル数が1:2になりますが、FFTの結果の総和も概ね1:2になります。（当然、周波数成分のあることを示す、絶対値の大きな値のある位置は変わりません。サンプル数回加算しているので絶対値がサンプル数にほぼ比例して大きくなります。）異なるサンプル数で求めた値を比較する必要があるなら、サンプリングデータか、FFTの結果のいずれかをサンプル数で割った値に正規化する必要があります。
３．２．高速化について
前述のように対象性がるので、サンプル数の半分を計算すれば、残りは符号の反転で求まります。また、exp()の部分もサンプル数分のテーブルに計算しておけます。
４．数値を入れて確認
４．１．16サンプルの例
下図のようなデータをFFTしてみます。 このデータは、「１Hzのサイン波を16Hzで1秒間サンプリングした」とも、「100Hzのサイン波を1600Hzで0,01秒間サンプリングした」とも解釈できます。
このデータをRのfft()で、演算すると以下のような結果が得られます。結果の値はサンプル数16で割ってあります。
上が次数部の値、下が虚数部をプロットしたものです。
虚数部の1と15に有意な値があり
f(t) = -0.5i・exp(ti) + 0.5i・exp(15ti)
= -0.5i・exp(ti) + 0.5i・exp(-ti) --- N/2で折り返し（共役）
= -0.5i(exp(ti) - exp(-ti))
= -0.5i・2・i・sin(t) ---- 2i・sin(x) = exp(ix) - exp(-ix)
= sin(t)
これは、サンプリング時間で一周するようなサイン波を示します。秒を単位とした時間sで表すと、
16Hzで1秒間サンプリングしたものなら、f(s) = sin(2π/1・s)
1600Hzで0.01秒間サンプリングしたものなら、f(s) = sin(2π/0.01・s) = sin(200π・s)
と、解釈できます。
４．２．位相をずらすと
おなじ、周期のサイン波とコサイン波を足して、下図、黒線の信号を作り16サンプルを作りました。
このデータをRのfft()で、演算すると以下のような結果が得られます。結果の値はサンプル数16で割ってあります。
上が次数部の値、下が虚数部をプロットしたものです。
実部、虚数部の1と15に有意な値があり
f(t) = (0.25-0.25i)・exp(ti) + (0.25+0.25i)・exp(15ti)
= 0.25・exp(ti) - 0.25i・exp(ti) + 0.25・exp(-ti) + 0.25i・exp(-ti) --- N/2で折り返し（共役）
= 0.25(exp(ti)+exp(-ti)) - 0.25i(exp(ti) - exp(-ti))
= 0.25・2・cos() - 0.25i・2・i・sin(t) ---- 2・cos(x) = exp(ix) + exp(-ix) 、2i・sin(x) = exp(ix) - exp(-ix)
= 0.5・cos(t) + 0.5・sin(t)
これは、サンプリング時間で一周するようなサイン波を示します。秒を単位とした時間sで表すと、
16Hzで1秒間サンプリングしたものなら、f(s) = 0.5cos(2π/1・s) + 0.5sin(2π/1・s)
1600Hzで0.01秒間サンプリングしたものなら、f(s) = 0.5cos(200π・s) + 0.5sin(200π・s)
と、解釈できます。
４．３．サンプリングレート1000Hzの例
１００Hz、３００Hz、８００Hzのサイン波を合成した波形を、1000Hzで1秒サンプリングしたデータを作って見ます。
最初の10サンプル(10ms)をプロットして見る。黒が合成した波形。他は、元のサイン波。
100Hz(赤）、300Hzは、概ねもとの信号がわかりますが、800Hz（青）については、本来8周期であることが分かりません。
800Hzだけを、元の信号を入れてプロットして見ます。
fft()の結果を見て見る。
有意そうな値のある箇所の値は以下の通り。(図では実部の200,800に離れたマークがあるが-14乗なので、ゼロと表示される。)
100 0-0.1666667i
200 0+0.1666667i
300 0-0.1666667i
700 0+0.1666667i
800 0-0.1666667i
900 0+0.1666667i
700,800,900の位置は、500で折り返して、それぞれ300,200,100に対応する周波数成分の強さを表します。サンプリングレート1000Hzで、1000サンプル(1秒）なので、１００Hz、２００Hz,300Hzに周波数成分があることになります。FFTの結果として求まる値を、スペクトルと言うようです。この、絶対値の２乗をパワースペクトルと言うようです。
５．まとめ
FFTの話は、３つの数値に依存しています。サンプリングレート(Hz)：S、サンプリング時間(秒):T、サンプル数：Nとすると
サンプリング時間：T=N/S
なので、2つの数値があれば良いことになります。
■分解能
FFTの結果の値の間の間隔（刻み）は、何Hzに相当するか
分解能 d = 1/T または、d = S/N
1秒間サンプリングすると１Hz、100msサンプリングすると10Hz。
■最小周波数
分解能に同じ。
■最大周波数
S/2-d または、N/(2・T)-d
と、言うことのようです。
可聴範囲の上限の20KHzまでを、小数点以下1桁で求めるとすると、サンプリングレート40KHz超、サンプリング時間10秒超、データ数400K個超が必要と言うことのようです。
試しに、Cと刻んである音叉を44.1KHzで10秒間サンプリングして見ました。この音叉は、C5(523.25Hz)なのだろうと思います。
|
インデクス
