さて、前回はProcessingを使ってwaveファイルのデータ読みこみ・波形表示まで実装した。今回が山場、スペアナの実装である。みんな大好きFFT(高速フーリエ変換)でスペクトラムを計算する。
スペアナのプログラムよりむしろFFTについてかなりの量を割いて説明しているので、ぜひ読んで欲しい(できるだけわかりやすく書いたつもりなのだが...)。
・フーリエ変換とは
まずは軽くフーリエ変換について触れておく。フーリエ変換とは、「時間についての非周期関数を、周波数の異なる正弦波の重ね合わせで表す」ことである。そしてどの周波数がどれだけ含まれているかを表すのが振幅スペクトルである。
フーリエ変換は(1)式のように定義される。ここで時間領域の関数\(f(t)\),周波数領域の関数\( F(\omega)\)はいずれも連続な複素関数である。
$$F(\omega)=\int_{-\infty}^{\infty} f(t)e^{-i\omega t}\mathrm{d}t\quad\cdots(1)$$
・離散時間フーリエ変換
ところが、実際の音声データをフーリエ変換にかけようと思った時、そのデータは一定時間おきにサンプリングされたもの、つまり離散的なもので、連続的ではない。
そこで、まず離散時間でのフーリエ変換を試みる。この音声データが整数\(n\)を用いて\( f[n]\)と表されるとき、その離散時間フーリエ変換は(2)式のようになる。
$$F(\omega)=\sum_{n=-\infty}^{\infty} f[n]e^{-i\omega n}\quad\cdots(2)$$
まあ積分が総和になっただけである。そして、この\( F(\omega)\)は周期\(2π\)の周期関数になる*1*2。
・離散フーリエ変換(DFT: Discrete Fourier Transform)
ところがまだ問題があって、このままだと周波数領域が連続的なのでコンピューターでは扱えない。そこで、総和の範囲を有限にしてしまって、その区間が延々と繰り返される周期関数にしてしまおう*3ということになったのが(3)式である。
$$F(\omega)=\sum_{n=0}^{N-1} f[n]e^{-i\omega n} \quad \cdots(3)$$
ただ範囲が\(0\)から\(N-1\)になっただけである。総和の範囲が\(0\)から\(N-1\)なので、この\( f[n]\)は周期\(N\)であり、基本角周波数は\(\displaystyle \frac{2\pi}{N}\)になる。すなわち、\(\omega\)は\( \displaystyle \frac{2\pi}{N}\)の整数倍の値を取ることになるから、整数\(k\)を用いて \( \displaystyle \omega=\frac{2\pi k}{N}\)と置き換えると、(3)式は
$$F\left(\frac{2\pi k}{N}\right)= F[k] = \sum_{n=0}^{N-1} f[n]e^{-i\frac{2\pi k}{N} n} \quad \cdots(4)$$
と書き換えることができる。これが離散フーリエ変換である。
また、\( F(\omega)\)が周期\(2\pi\)であったことを考えれば、\( F[k]\)が周期\(N\)であることもわかる。つまり\( f[n], F[k]\)はともに離散的で、かつ周期\(N\)で周期的なのである。
・高速フーリエ変換(FFT: Fast Fourier Transform)
とりあえず(4)式さえあれば計算はできるのだが、\(0\leq n < N, 0\leq k < N\) なので真面目にやると計算量が\( O(N^2)\)必要になる。
\(N\)が小さいうちはそこそこ早く終わるが、\(N\)が大きくなってくると時間がかかりすぎるということで高速化したのがFFTである。
実はFFTを使うには制限があって、\(N\)が2の累乗でなければならない。その理由はすぐわかるので、まずはFFTの原理を説明していく。
FFTでは、(5)式のようにまず\(n\)の偶奇で項を分ける。
$$F[k] = \sum_{n=0}^{\frac{N}{2}-1} f[2n]e^{-i\frac{2\pi k}{N}\cdot 2n}+\sum_{n=0}^{\frac{N}{2}-1} f[2n+1]e^{-i\frac{2\pi k}{N}\cdot (2n+1)}\quad\cdots(5)$$
\(e\)の右肩がゴチャゴチャするので、\( \displaystyle W=e^{-i\frac{2\pi}{N}}\)と置き換える。すると(5)式は(6)式のようになる。
$$F[k] = \sum_{n=0}^{\frac{N}{2}-1} f[2n]W^{k\cdot 2n}+W^k\sum_{n=0}^{\frac{N}{2}-1} f[2n+1]W^{k\cdot 2n}\quad\cdots(6)$$
つまり奇数の項には\( W^k\)がかかる。
そして、それぞれの項についてまた\(n\)の偶奇が分けられて(7)式のようになる。
今度は奇数の項に\( W^{2k}\)がかかっている。
$$\begin{alignat*}{3}
F[k] &=& && &\sum_{n=0}^{\frac{N}{4}-1} f[4n]W^{k\cdot 4n} &\;+\; W^{2k}&\sum_{n=0}^{\frac{N}{4}-1} f[4n+2]W^{k\cdot 4n}\\
&& &+ W^k& &\sum_{n=0}^{\frac{N}{4}-1} f[4n+1]W^{k\cdot 4n} &\;+\; W^{k}W^{2k}&\sum_{n=0}^{\frac{N}{4}-1} f[4n+3]W^{k\cdot 4n}\quad\cdots(7)
\end{alignat*}$$
そしてこれを\(\log N\)回繰り返すと、(8)図のようになってしまうのである。\(N\)が2の累乗でなければならない理由がおわかりいただけただろうか。
(図は\(N=8\)のとき。やっつけなのは勘弁してほしい。)
これで、各々の\(k\)について\(F[k]\)が\(O(N\log N)\)で求められる..........
........???????????
それぞれの\(k\)について\(O(N\log N)\)で計算するのだから、最終的な計算量が
\( O(N^2 \log N)\)になってむしろ悪化している。あれれ~おかしいぞ~。(棒読み)
・\( e^{i\theta} = \cos\,\theta + i\sin\,\theta\) の使いどころ
この謎を解く鍵は、「\(W\)が\(k\)についての複素指数関数、すなわち周期的である」というところにある。\( \displaystyle W=e^{-i\frac{2\pi}{N}}\)なのだから、
$$\begin{alignat*}{2}
W^{k+N}&=&&e^{-i\frac{2\pi}{N}(k+N)}&&\\
&=&&e^{-i\left(\frac{2\pi k}{N}+2\pi\right)}&&\\
&=&&e^{-i\frac{2\pi k}{N}}&&\\
&=&&W^k &\cdots(9)&
\end{alignat*}$$
が成り立つのである。
試しに(8)図に\( F[4]\)を書き加えてみると(10)図のようになる。
赤枠の中は\( F[0]\)の計算結果を流用できるのがわかるだろうか。
そしてそれを\( F[0]\)から\( F[N]\)まで縦に並べて重ねたのがいわゆるバタフライダイアグラムだ...図(11)。
赤破線とその右端の数字は、\(W\)をその数字の分だけ冪乗して足すという意味である。
これを、黒破線で区切られた部分を一つのユニットとして左から計算していけば、ユニット当たり\(O(N)\)、それが\(\log N\)ユニットあるので、晴れて\(O(N\log N)\)のアルゴリズムの完成である。
・ビット逆転
ここで、左側の\( f[n]\)の数字の並びに戸惑った人もいるのではないだろうか。結論から言ってしまうと、ビットが逆転しているのである。反転ではない、逆転である。
例えば\(N=8\)のとき、上から2進数で000->100->010->110->001->101->011->111
となっている。このそれぞれを逆から読んでみると、000->001->010->011->100->101->110->111
となってちゃんと順番になっている。
考えてみれば特に不思議でもなんでもなく、偶奇すなわち下位ビットが0か1かで項を分けているのだからそりゃ逆転するわな、という感じである。
ここまでで理論は終わりで、あとは実際に実装していく。
・やっとProcessingの出番
前述のとおりフーリエ変換は複素関数について定義されているので、簡単に自前の複素数クラスを用意してみた。
class Complex{
private double re;
private double im;
public Complex(double real,double imaginary){
re=real;
im=imaginary;
}
public Complex(double radius,float arg){
re=radius*cos(arg);
im=radius*sin(arg);
}
public double Re(){
return re;
}
public double Im(){
return im;
}
public void Add(Complex c){
re+=c.Re();
im+=c.Im();
}
public Complex Product(Complex c){
return new Complex(re*c.Re()-im*c.Im(),re*c.Im()+im*c.Re());
}
public double Abs(){
double re_sq=re*re;
double im_sq=im*im;
return Math.pow(re_sq+im_sq,0.5);
}
public String toString(){
return String.format("%f+%fi",re,im);
}
}
使うのは和、積、大きさくらいなので最小限しか用意していない。
コンストラクタには\(a+bi\)の形と、\(W\)の計算で便利なので極形式\( r({\rm cos}\,\theta +i{\rm sin}\,\theta)\)の両方を用意してみた(ド・モアブル万歳)。toString()
はデバッグ用につけてみたら案外便利だった。
\(W\)の累乗は事前に計算しておいて配列に格納しておく。そこまで軽量化にこだわってはいないので\(N\)通り全部計算している。
W=new Complex[N];
for(int i=0;i<N;i++){
W[i]=new Complex(1,-2*i*P/N); //極座標表記が活躍
}
それと、ビット逆転の配列も作っておく。これは某所で見つけたビット逆転インクリメントのコードをそのまま写した。
誰が考えたのか知らないが、よくもまあこんなのが思いつくものだ。
void ReverseIncrement(){
revBit=new int[N];
int p=0;
int m=1<<power;
for(int i=1;i<N;i++){
revBit[i-1]=p;
p^=m-m/2/(i&(-i&m-1));
}
revBit[N-1]=p;
}
窓関数も忘れてはいけない*4。今回は矩形窓、サイン窓、ハン(ハニング)窓、ハミング窓、ブラックマン窓、Vorbis窓などを自由に切り替えられるようにしてみた。
float Window(int i){
switch(window){
case 0:return 1; //Rect
case 1:return sin(P*i/N); //Sine
case 2:return 0.5-0.5*cos(2*P*i/N); //Hann
case 3:return 0.54-0.46*cos(2*P*i/N); //Hamming
case 4:return 0.42-0.5*cos(2*P*i/N)+0.08*cos(4*P*i/N);//Blackman
case 5:return sin(P/2*pow(sin(P*i/N),2)); //Vorbis
default:return 0;
}
}
ここまできたら準備完了。あとはFFTを組んでいく。
void FFT(){
if(timePos+N>sampleNum)return;
Complex[] prev=new Complex[N];
Complex[] next=new Complex[N];
for(int i=0;i<N;i++){
prev[i]=new Complex(Window(revBit[i])*waveData[0][timePos+revBit[i]],0);
next[i]=new Complex(0,0);
}
for(int i=0;i<power;i++){
int k=1<<i;
for(int j=0;j<N;j++){
if(j%(k*2)<k){
next[j].Add(prev[j]);
next[j+k].Add(prev[j]); //...(12)
}else{
int w=N/k/2;
next[j].Add(prev[j].Product(W[w*(j%(k*2))]));
next[j-k].Add(prev[j].Product(W[w*((j-k)%(k*2))]));
//...(13)
}
}
System.arraycopy(next,0,prev,0,N); //...(14)
next=new Complex[N];
for(int a=0;a<N;a++)next[a]=new Complex(0,0);
}
for(int i=0;i<N;i++){
freqData[i]=prev[i].Abs(); //...(15)
}
}
細かい説明は省くが、偶数側のときはそのまま足し...(12)、奇数側のときは\(W\)の冪乗をかけてから足している...(13)。ユニットひとつの計算が終わったら配列の中身を移し替えて次のユニットの計算をする...(14)。
全部終わったら大きさを取れば振幅スペクトルが得られる...(15)。
これでFFTが終わったので、あとはプロットすればスペアナが完成する。
void PlotFreq(){
translate(0,height-25);
textAlign(CENTER,BOTTOM);
text(String.format("Amp:%s, Freq:%s, Sample:%.0f, Window:%s",
logAmp?"Log":"Linear",logFreq?"Log":"Linear", //...(18)
pow(2,power),windowName),width/2,-5);
nyqFreq=logFreq?log(sampleRate/2):sampleRate/2;
fundFreq=logFreq?log(sampleRate/N):sampleRate/N;
freqWidth=nyqFreq-fundFreq;
for(int i=1;i<=N/2;i++){
float freq=(float)sampleRate/N*i;
float freqPos=logFreq?log(freq):freq; //...(16)
float pos=20+(freqPos-fundFreq)/freqWidth*localWidth;
//...(20)
float val=(float)freqData[i];
float fLevel=Float.NaN;
if(logAmp){ //...(17)
float v=log(val)+thres;
fLevel=50*fAmpZoom*(v>0?v:0);
}else{
fLevel=ampRate[power-10]*fAmpZoom*val;
}
if(fLevel>height){
fLevel=height;
stroke(255,0,0); //...(20)
}else{
stroke(255);
}
line(pos,-30,pos,-fLevel-30);
textAlign(i==1?LEFT:RIGHT,BOTTOM);
if(i==1||i==N/2)text(String.format("%.0f",freq),pos,-5);
//...(19)
}
stroke(255);
line(20,-30,width-20,-30);
}
大きなポイントとしては両軸それぞれリニア/対数を切り替えられるようにしたり...(16,17)、軸のスケールやFFTの\(N\)の値や窓関数の種類、基本周波数とナイキスト周波数を下部に表示したりしている...(18,19)。また、スケールを切り替えても全体の幅が一定になるように計算している...(20)。
細かいところとしては、振幅の値が画面外に出るレベルのときはバーが赤くなるようにしてみたりしている...(21)。
実際に動いているのがこちら。
Fastと謳うだけあって\(N=16384\)でも30fpsは堅い。むしろプロットするほうが重いまである。
以上でスペアナが完成した。しかし、ここまででは波形とスペアナの映像しかなくて物足りないので次回は音声をつけていく。
こんなに長い記事を書いたのは初めてでメチャクチャ時間がかかったが、いい復習になった。
余談:「2π/N」の「π/」の部分をπスラッシュと誤認されているようで大変遺憾です。