FFTで任意区間の周波数だけを計算する 音の周波数測定（AudioFC）を作ったのですが、重いので範囲を限って計算すれば早くなると思いました。試したら全体を計算するよりずっと遅くなってしまいました。サンプル数が2のｎ乗の計算で65536サンプルを計算する時間で、単純な式通りの計算方法では何十周期分しか計算できないようです。
多くのサンプルの一部を計算するのは、分解能を維持するためです。 「１．FFTのおさらい」には、各周波数を独立に計算する場合のことを考えて見ました。 「２．サンプル数が2の整数乗の場合の計算」では、サンプル数が2の整数乗の場合の計算方法に、計算の不要箇所を除外する方法を付け加えることを考えて見ました。 結論を言うと、65536サンプル、1.486秒の計算に、410msほど要していましたが、60-1100Hzに制限すると260msで計算できました。 ただし、とにかく早くしようとしているわけではありません。C#で、素直に複素数計算、指数計算のままでできる範疇で考えています。
また、「サンプル数が2の整数乗の場合の計算」は、元々最速なので、これをどうにかしようと言うものでもありません。必要なのは一部の周波数なので、計算の不要箇所があるはずだ、と言うことです。単一の周期の計算でも、全てのサンプルに指数関数を乗じて合算することが必要なのは変わりません。計算の最初の段は、全ての計算が必要です。最終段は、単一の周期の計算なら、1つのみ計算して他は省略できます。この場合、段数が多ければ、計算の必要箇所は3角形になるでしょうが、この最良でも、1/2までしか計算量を減らせません。段数が多く（サンプル数が多く）、これに比べて、必要な周期の範囲が極端に少ない場合以外は効果がありません。
もともとの、FFTプログラムはテーブルを持っていません。今回の版は、大きなテーブルが必要になります。最大の利点は、結果を関数の戻り値としていますが、これが小さくなることだと思います。 １．FFTのおさらい Rを使って、16サンプルの例を考えます。
この16サンプルの期間に周期1,2,3のサインカーブのデータcurve1、curve2、curve3を作っておきます。 １．１．Rの例の注意- Rの配列のインデクスは、1からです。
- -1.376151e-15-4.21321e-16iのような計算結果が表示されますが、10の-15乗と言った数字で、ゼロと読んでください。
- 計算で使った変数は、変数名を入れて、「Enter」を押せば表示されます。
- plot(e1)などで、変数をグラフ表示できます。
１．２．角速度FFTの計算を考えるときは実際の時間(サンプリング時間)は出てきません。計算対象のサンプル中での周期を計算します。FFTの結果は、サンプル数と同じになっていますが、結果のインデクスが周期になります。結果[1]が周波数1(回/サンプリング時間)の強さ、結果[2]が周波数2の強さ，．．．となります。（インデクスはゼロから始まるとして）
サンプル中での時間経過を示すのに角速度を使います。16サンプルなら、2π/16がサンプルの間隔を表します。サンプル中に2周期のサイン波は角速度が2π/8*2になります。 １．３．最大周波数計算結果のインデクス（0から始まるとして）がサンプル中での周波数を表しますが、これはサンプル数の1/2までです。
この先を計算する必要はありません。計算すると、前半の計算結果の共役数が、折り返しの順番で算出されます。イメージは、「周波数の分解能は？」を参照ください。
16サンプルの例だと周波数7までが算出可能です。これが22050Hzでサンプリングしたデータなら、(22050/2) * (7/8)に相当します。 １．４．各周期の計算
16サンプルで1周期のサンプルの間隔は、r=2π/16です。
サンプリングデータがcurve[16]にあるものします。（インデクスはR流に1から）
また、e_n = exp(complex(Re=0,Im=-r*n))と書くことにします。 周期1は、e_0,e_1,...,e_15を、それぞれcurve[1],curve[2],...,curve[16]に掛け、合計します。
周期2は、2倍の速度（2*r）で進むので、e_0,e_2,...,e30を、それぞれcurve[1],curve[2],...,curve[16]に掛け、合計します。
周期3は、3倍の速度（3*r）で進むので、e_0,e_3,...,e45を、それぞれcurve[1],curve[2],...,curve[16]に掛け、合計します。 左図は、e_nをプロットしたもので、2πで一周する周期性があります。前の計算を続けるとe_225まで出てきますが、同じ値が繰り返し使われています。 １．５．Rの例 まず、周波数1のサンプルデータについて、周期3まで計算して見ます。
0-8i
-1.376151e-15-4.21321e-16i
-4.40945e-16-5.55112e-17i
と、算出されました。周期2,3の計算結果は、10の-15乗と言った数字で、ゼロと見なします。 同様に、周波数2,3のサンプルデータについて、計算して見ます。それぞれ、周期2、周期3に0-8i
が求まり、他はほぼゼロになります。 最後に、位相をずらしても、結果に影響ないことを試して見ました。 > N <- 16 #サンプル数
> r <- 2*pi/N#サンプル間の角速度
> curve1 <- sin(seq(0,by=2*pi/16,length.out=16))
> curve2 <- sin(seq(0,by=2*pi/16*2,length.out=16))
> curve3 <- sin(seq(0,by=2*pi/16*3,length.out=16))
> ##############################################
> #1周期のサンプル
> curve <- curve1
> #1周期の計算
> t1 <- seq(0,by=r,length.out=16)
> e1 <- exp(complex(re=0,im=-t))
> sum(curve*e1)
[1] 0-8i
> #2周期の計算
> t2 <- seq(0,by=r*2,length.out=16)
> e2 <- exp(complex(re=0,im=-t2))
> sum(curve*e2)
[1] -1.376151e-15-4.21321e-16i
> #3周期の計算
> t3 <- seq(0,by=r*3,length.out=16)
> e3 <- exp(complex(re=0,im=-t3))
> sum(curve*e3)
[1] -4.40945e-16-5.55112e-17i
> ##############################################
> #2周期のサンプル
> curve <- curve2
> #1周期の計算
> t1 <- seq(0,by=r,length.out=16)
> e1 <- exp(complex(re=0,im=-t))
> sum(curve*e1)
[1] -1.309283e-15-4.21321e-16i
> #2周期の計算
> t2 <- seq(0,by=r*2,length.out=16)
> e2 <- exp(complex(re=0,im=-t2))
> sum(curve*e2)
[1] 0-8i
> #3周期の計算
> t3 <- seq(0,by=r*3,length.out=16)
> e3 <- exp(complex(re=0,im=-t3))
> sum(curve*e3)
[1] -5.876376e-16-6.889021e-16i
> ##############################################
> #3周期のサンプル
> curve <- curve3
> #1周期の計算
> t1 <- seq(0,by=r,length.out=16)
> e1 <- exp(complex(re=0,im=-t))
> sum(curve*e1)
[1] 8.309325e-16-5.55112e-17i
> #2周期の計算
> t2 <- seq(0,by=r*2,length.out=16)
> e2 <- exp(complex(re=0,im=-t2))
> sum(curve*e2)
[1] 1.422202e-15-6.88902e-16i
> #3周期の計算
> t3 <- seq(0,by=r*3,length.out=16)
> e3 <- exp(complex(re=0,im=-t3))
> sum(curve*e3)
[1] 0-8i
> ##############################################
> #位相：2周期のサイン波を90度ずらして
> curve <- sin(seq(pi/2,by=2*pi/16*2,length.out=16))
> #1周期の計算
> t1 <- seq(0,by=r,length.out=16)
> e1 <- exp(complex(re=0,im=-t))
> sum(curve*e1)
[1] -1.944788e-15+1.59825e-15i
> #2周期の計算
> t2 <- seq(0,by=r*2,length.out=16)
> e2 <- exp(complex(re=0,im=-t2))
> sum(curve*e2)
[1] 8-0i
> #3周期の計算
> t3 <- seq(0,by=r*3,length.out=16)
> e3 <- exp(complex(re=0,im=-t3))
> sum(curve*e3)
[1] -1.020261e-15-8.11254e-16i
>
１．６．計算結果の値について 通常、16ビットでサンプリングしたデータは符号付の16ビット整数で、-32768から32767の範囲で振幅を表しています。例では、単にsin()の結果なので、-1から1の範囲で振幅を表しています。同じスケールにするために、32767を掛けて計算して見ます。
c<-curve1*32767
fft(c)
すると、2番目のデータが
0.000000e+00-2.621360e+05i
と、なりますた。
0+262136i
で
262136 = 32767 * 8
です。 サンプル数を32にして計算して見ます。
c <- sin(seq(0,by=2*pi/32,length.out=32))
fft(c)
すると、2番目のデータが
0.000000e+00-1.600000e+01i
となり、
0+16i
です。 FFTの計算結果は、振幅に、指数関数（-1から1）を掛けて、サンプル数分の加算をした値なので、同じ状態を計測しても、振幅の表現方法とサンプル数で値が変わります。
「440Hzの成分が多い」と言った見方をする場合の、周波数は絶対的な値ですが、「多い」の部分は、他と比較することができません。
FFTの結果の値は、同一条件でサンプリングされたデータを、同一条件で計算された場合に相対的に比較できます。
どうしても、比較しようとする場合ですが、サンプル数の影響については、入力か結果をサンプル数で割って一般化できます。
その上で、FFTの計算結果の絶対値の2乗を、0から-120dBに引き当てて表示にしているようです。FFTの計算結果の計算上の最大値を0dBとし、0を-120dB（または、-100dB）としているのだと推測します。 １．７．指数関数の周期性 x[16,16]に、計算に使った、指数関数の値をセットします。周期1の計算に使った変数の値の何番目と一致したかを、ID[16,16]に記録します。ID[]は、-1で初期化しておきます。 > r <- 2*pi/16
> x<-matrix(0,16,16)
> for(i in 0:15){ x[i+1,]<-exp(complex(re=0,im=-(seq(0,by=r*i,length.out=16))))}
>
> ID <- matrix(-1,16,16)
> for(i in 1:16){
+ for(j in 1:16){
+ for(k in 1:16){
+ if((abs(Re(x[2,k])-Re(x[i,j]))<0.0001)&&(abs(Im(x[2,k])-Im(x[i,j]))<0.0001)) {
+ ID[i,j] <- k
+ break;
+ }
+ }
+ }
+ }
> ID <- ID-1 結果は、以下のようになりました。（ゼロから始まるように最後にID<ID-1としてあります。） [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
[1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[2,] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
[3,] 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14
[4,] 0 3 6 9 12 15 2 5 8 11 14 1 4 7 10 13
[5,] 0 4 8 12 0 4 8 12 0 4 8 12 0 4 8 12
[6,] 0 5 10 15 4 9 14 3 8 13 2 7 12 1 6 11
[7,] 0 6 12 2 8 14 4 10 0 6 12 2 8 14 4 10
[8,] 0 7 14 5 12 3 10 1 8 15 6 13 4 11 2 9
[9,] 0 8 0 8 0 8 0 8 0 8 0 8 0 8 0 8
[10,] 0 9 2 11 4 13 6 15 8 1 10 3 12 5 14 7
[11,] 0 10 4 14 8 2 12 6 0 10 4 14 8 2 12 6
[12,] 0 11 6 1 12 7 2 13 8 3 14 9 4 15 10 5
[13,] 0 12 8 4 0 12 8 4 0 12 8 4 0 12 8 4
[14,] 0 13 10 7 4 1 14 11 8 5 2 15 12 9 6 3
[15,] 0 14 12 10 8 6 4 2 0 14 12 10 8 6 4 2
[16,] 0 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 指数関数は、サンプル数分計算すれば良いことが分かります。
行は、それぞれ周期0,1,...,15に相当しますが、どの行も、列方向に、周期の間隔で数字が増えていきます。周期3を見ると、
[4,] 0 3 6 9 12 15 2 5 8 11 14 1 4 7 10 13
で、15と14の2箇所の後で一周しています。15の箇所は、計算時は、exp(-2π/16*15i)、次はexp(-2π/16*18i)を計算していまが、
exp(-2π/16*2i) = exp(-2π/16*18i) (Rで計算は、exp(complex(im=-2*pi/16*18))のようにします。)
と、言うことです。(18は、16で一周し、さらに2進む。と読める。） １．８．計算速度について結局、十分な速度になりませんでした。
サンプル数が2のｎ乗の計算で65536サンプルを計算する時間は、416ms程度です。
前述の方法だと、200周期の計算で、2.8秒かかってしまいました。
指数関数の部分をテーブルにし、さらに、複素数の積和を、実数、虚数別に直接書いて見ました。（Complexの計算ではComplex型のコンストラクタが呼ばれる。）
1.1秒になりました。
複素数の積和をコメントアウトして実行すると、86msになります。 ２．サンプル数が2の整数乗の場合の計算２．１．計算方法サンプル数がNなら、サンプル全体で1周期になる角速度r=2π/Nを基に計算します。
例えば、周期1の計算は、サンプルに対して順番に
exp(-r0i)、exp(-r1i)、exp(-r2i)、．．．、exp(-r(N-1)i)を乗じて、合計します。
ここでexp(-n・r・m・i)を、En_mと書くことにします。
例えば、
周期1の計算は、サンプルに対して順番に
E1_0、E1_1、E1_2、．．．、E1_(N-1)を乗じて、合計します。
と、なります。 この表記で、8個のサンプリングデータa,b,c,d,e,f,g,hを考えて見ます。サンプル数が2のK乗の場合、K段の計算をします。各段で乗ずる値を、8=2の3乗なので、3段で書くと下表になります。 | | 0段 | 1段 | 2段 | 0 | a | a+eE4_0 | [0]+[2]E2_0 | [0]+[4]E1_0 | 1 | e | a+eE4_1 | [1]+[3]E2_1 | [1]+[5]E1_1 | 2 | c | c+gE4_0 | [0]+[2]E2_2 | [2]+[6]E1_2 | 3 | g | c+gE4_1 | [1]+[3]E2_3 | [3]+[7]E1_3 | 4 | b | b+fE4_0 | [4]+[6]E2_0 | [0]+[4]E1_4 | 5 | f | b+fE4_1 | [5]+[7]E2_1 | [1]+[5]E1_5 | 6 | d | d+hE4_0 | [4]+[6]E2_2 | [2]+[6]E1_6 | 7 | h | d+hE4_1 | [5]+[7]E2_3 | [3]+[7]E1_7 |
この表を以下のルールを使って書き換えます。 E1_0,E2_0,E4_0は、1+0i=1です。 E1_4,E2_2,E4_1は、-1+0i=-1です。 E1_5,E1_6,E1_7は、それぞれ、-E1_1,-E1_2,-E1_3 です。
| | 0段目 | 1段目 | 2段目 | 0 | a | a+e | [0]+[2] | [0]+[4] | 1 | e | a-e | [1]+[3]E2_1 | [1]+[5]E1_1 | 2 | c | c+g | [0]-[2] | [2]+[6]E1_2 | 3 | g | c-g | [1]-[3]E2_1 | [3]+[7]E1_3 | 4 | b | b+f | [4]+[6] | [0]-[4] | 5 | f | b-f | [5]+[7]E2_1 | [1]-[5]E1_1 | 6 | d | d+h | [4]-[6] | [2]-[6]E1_2 | 7 | h | d-h | [5]-[7]E2_1 | [3]-[7]E1_3 |
1行目[1]を確認して見ます。
1段目の計算結果[1]の内容
(a-e)+(c-g)E2_1
1段目の計算結果[5]の内容
(b-f)+(d-h)E2_1
を、加算して、
(a-e)+(c-g)E2_1+((b-f)+(d-h)E2_1)E1_1
a + bE1_1 + cE2_1 + dE2_1・E1_1 - e - fE1_1 - gE2_1 - hE2_1・E1_1
と、なっています。
本来は、
aE1_0 + bE1_1 + cE1_2 + dE1_3 + eE1_4 + fE1_5 + gE1_6 + hE1_7
です。
ここで、表を変換したときのルールと、exp(A+B)=exp(A)*exp(B)の公式を使うと同じ計算をしていることが分かります。たとえば、E1_4は-1なので-e=eE1_4です。E1_2は-2π/8*2の位置で、E2_1は、-2π/4*1の位置で同じ位置を示します。計算方法として、この表の見方ですが、FFTの結果を求める、サンプル数と同じ大きさの配列[0]から[7]があって、サンプルデータを1つずつ入れておきます。入れる順番にも工夫があって、最初のペアを作る順番になっています。計算は常に2個のペアで行われ、ペアの2個の間隔は、段数に応じて、1,2,4,...と変化します。最終段まで計算すると、配列[0]から[7]に周期0から7に相当するFFTの計算結果が求まっていると言うことです。
ペアの計算は、同じ値の符号の違うものの加算になっています。
基準の角速度をr=2π/N、K=log2(N)として、k(0,1,...)段目は(2の(K-k-1)乗)・ｒの角速度の計算をしています。
全周期を計算するには、全てのサンプルの、E1_1との積、E1_2との積、．．．E1_(N-1)との積、の全ての組が必要ですが、指数関数は、異なる角速度の指数関数の積で表されるので、計算の段ごとに、周波数方向に並列的に計算されています。
また、各段ごとには、各ペアのそれぞれは、1度しか使われません。格段ともN/2ペアの計算が行われます。 ２．１．計算範囲を制限するにはまず、確かなことは、どの行の計算でも全てのサンプルに係数を乗じた和を計算するので、0段目の計算は全て行う必要があると言うことです。
また、最終段は、選択的に計算可能と言うことです。 前述の表の、行方向は、サンプリング時間当たり、周期0のスペクトル、周期1のスペクトル、．．．に相当します。この周波数方向をmとします。また、計算の段数を示す、列方向をkとします。kは、ゼロから数えることにします。
0≦k＜K ただし、K=log2(サンプル数)です。
各段とも、2個の対を計算しますが、k段目で、m行を計算するときには、前段の計算結果の、ｍ行とｍ＋(2のk乗)が必要です。（サンプル数の半分まで計算することを前提として考えています。）
周期A≦m≦Bを計算することを考えて見ます。サンプル数8(K=3)以上で考えます。
最終(K-1)段の開始位置Aでは、K-2段目のA行のデータが必要ですが、ここが対の若い方なら、AがK-2段目の開始位置になります。そうでなければ、A -（２の（K-2）乗）が開始位置になります。
最終(K-1)段の終了位置Bでは、K-2段目のB行のデータが必要ですが、ここが対の若い方でないなら、BがK-2段目の終了位置になります。もし、若い方なら、B -（２の（K-2）乗）が開始位置になります。ｍが対の若い方かどうかは、mをd=(2の(k+1)乗)で割って、剰余がd/2より小さいかどうかで判断します。 ２．２．結果44100Hzのサンプリングレートで、1.486秒に相当する65536サンプルから、60-1100Hzの部分を計算してみました。分解能は0.675Hzです。サンプリング時間に対して、0から32768の周期のうち、90-1635の1546周期分を計算します。
この計算は、0から15段で計算されます。本来なら、全ての段で、32768組の計算が行われます。計算不要箇所を除くと、11から15段で計算量が減らせました。それぞれ、32768,24736,12368,6184,3092,1546でした。
時間的には、全数で460msでしたが、260msになりました。計算量は余り減っていないのに時間的には、短くなっています。これは、必要な計算をコンストラクト時に行い、テーブルを持つようになったためです。 |
インデクス
