2013/02/03

DFT 基礎実験

2012年1月の記事にプログラム追加。Javaで作ったプログラムをC言語に移植し、iDFTで信号を戻す部分を追加した。

DFT(discrete Fourier transformation)とは離散フーリエ変換といって、離散信号の時間領域の波形を周波数領域に変換する手法のこと。
下の例では波形をDFTによって分析してみる。

上記波形をDFTで周波数領域に変換すると、下のように1000Hzと10000Hzの音が含まれていることがわかる。

DFTは音関連のプログラムを作る上で、とても基礎的なことなので、プログラムを作りながら確認してみることにした。

基本となる式はこれで、H(k)は最終的に欲しい各周波数の振幅。h(n)は分析する各サンプル。Nはサンプルの数。この逆の操作となる逆フーリエ変換は下のようになる。

上記の式のままでは使えないので、使いやすいように下のように式を変形する。下付き添字のrは実部を表し、iは虚部を表す。数式はあまり詳しくないので、書き方が間違っているかもしれないが、動作としては一番下のプログラムソースを参考にしてください。

逆変換は下のようになる。

フーリエ変換のプロセスは結構面白いので、そのうち子ども向けに解説しようかと思っているが、今回は流れだけ簡単に。

まず、周期性のある波形は、その整数倍のsin、cos波に分解できるという、フーリエの理論がある。具体的には、調べたい波形に整数倍の sin、cos波を掛けて積分すると、各周波数成分が、どれぐらい含まれているかが出てくるという仕掛け。含まれていなければ0で、含まれていれば、振幅のレベルが出てくる。このプロセス自体はシンプルなのだが膨大な計算量なのが欠点。より高速に処理するためのFFT(高速フーリエ変換)という手法もあるのだが、今回はシンプルなDFTをまず作ることにした。上記の式をプログラムで実行するのは簡単で、入れ子のfor文数行で完成してしまう。ただ式に虚数 i (電気ではj)があるので複素数を扱う。周波数スペクトルに使う場合は、実部と虚部の数値をそれぞれ2乗して足した平方根が、周波数スペクトルの振幅スペクトルとして利用できるようになる。式はこんなかんじ。

このXnはリニア表示なので、さらにAudacityと同じようにdB表示させるには、最大音量1.0からの比を求める。
20log(Xn/1.0)
と計算すると、最大音量0dBのdB表示となる。

実装と構成

実際の音声ファイルを扱えるように作ってみた。32bit float のwav(リトルエンディアン)を任意のサンプル数だけ読み込んで、それらをDFTする。そのデータをテキストファイルとして保存するというもの。

出力テキストファイル
Audacityと同じ内容形式のテキストファイルで出力するようにした。ファイルに書かれる内容は、以下のように周波数(Hz)とレベル(dB)。2行目以降がデータになっている。

Frequency (Hz) 
21.533203
43.066406
64.599609
86.132813
107.666016
129.199219
150.732422
172.265625
193.798828
215.332031
Level (dB)
-30.319666
-30.831787
-30.607641
-30.594364
-30.720446
-30.352894
-30.213865
-29.939169
-30.219339
-30.417690


Audacityの出力には0Hzがない。その理由が分からない・・・ 0dBは直流成分なので、結構重要だと思う。そんなDCオフセットされたデータは不必要ということなのだろうか。通常まともな音声には直流成分は0なので。試しにAudacityで無音ファイルを作成して、自作DC-offsetを使って0.5の直流ノイズにして、周波数スペクトルを見てみた。Rectangular window(窓をなし)で見ると何も表示されない。直流成分は無視されているようだ。理由がありそうだが、とりあえず自作では0dBも入れることにした。DCノイズを窓なしでDFTすると、下のように、ちゃんと0dBだけが大きなレベルになっている。

Frequency (Hz) 
0.0
44.1
88.2
132.3
176.4
220.5
264.6
Level (dB)
-3.0102999566398116
-333.07532324911443
-326.51134009780384
-324.8872661949541
-317.5628719756075
-321.77479597147874
-332.17644973225055


窓関数

通常切り出した波形は、きれいな周期関数にならない。そこで窓関数を掛けて両端を整える。そして無限に繰り返してつなぎ合わせると仮定することでフーリエ変換を適用できる。窓関数を掛けないと、つなぎ目に妙な段差ができて高域成分がおかしくなってしまうのでDFT、FFTでは必須。
やっていくと窓関数で、ちょっと疑問に思うことがあった。サンプル数が奇数なら中央値はひとつで問題ないのだが、偶数の場合どうしているのか?ということ。下は11サンプルのハミング窓の例。こんな係数が並ぶ。

0 wHamming= 0.08000000000000002
1 wHamming= 0.16785218258752427
2 wHamming= 0.39785218258752425
3 wHamming= 0.6821478174124759
4 wHamming= 0.9121478174124759
5 wHamming= 1.0
6 wHamming= 0.9121478174124759
7 wHamming= 0.6821478174124759
8 wHamming= 0.39785218258752425
9 wHamming= 0.16785218258752427
10 wHamming= 0.08000000000000002


切り出すサンプル数

FFTでは2のべき乗しか使えないので迷うこともないのだが、DFTは自由。上の偶数奇数とかの問題もある。実際にいろいろ試してみると、必要に応じて使い分けるという感じだ。低音の分解能を上げるにはサンプル数を増やす必要があるし、高域だけ見たければ少なくてもそこそこ見れる。またサンプル数を増やすということは、それだけの時間の範囲が必要なので、分析したい部分が短ければ加工するなど、別の手段が必要になってくる。うまく使うと、サンプル数の自由度からFFTでは見れなかったものまで見れるようになる。リアルタイムに処理するならFFTではないと話にならないが、オフラインで一部だけを分析したいならDFTの方が良い場合がありそうだ。

具体的に

サンプリング周波数が44100Hzだとして、100サンプルでDFTにかけたとする。このときの解像度は44100/100で、441Hzの倍数となる。つまり周波数域で見たときには0~22050Hzまで扱えるオーディオを100/2等分していることになる。特に音の分析においては、低域は結構重要で、それなりの解像度が求められる。音楽で使う重要な基音は数10から数100Hzまで。100サンプルでは、低いほうから、0Hz、441.0Hz、882.0Hz、1323.0Hz、1764.0Hz、2205.0Hz、2646.0Hz・・・ と並んでいく。これってかなりアバウト。ベースの基音などは1段目にすべてが入っている。正確な基音など出るはずもない。ということで実際には4000サンプル以上欲しいことが多いように思う。せめて10Hz前後の解像度は欲しいと思う。さらに100サンプルで分解した周波数を並べていくと、22050Hzを超えて、22491Hz、22932Hz、23373Hz・・・44100Hzになってしまう。これって、そもそもサンプリング周波数44100Hzでは当然記録できない範囲。どういうことかというと、100等分した解像度は22050Hzを軸に対称になっている。周波数スペクトルの表示などで使うのは半分だけ。周波数領域から時間領域に逆変換するときは、すべて使う。まぁコピーして使うのもありだが。保存量を最小限にしたければ半分あれば間に合うということ。

実際にすべてをプロットしてみた。使った波形は500、1000、5000、10000、15000、20000Hzのサイン波を合成したもの。ナイキスト周波数の22025Hzを堺に対称になっているのが分かると思う。

ついでに振幅スペクトルにしていない状態の実部と虚部も出力してみた。グリーンが実部で線対称、ピンクの虚部は点対称。


各サンプルの間は?

サンプリング周波数が44100Hzで1000サンプルにした場合、0、44.1、88.2Hz・・・と刻まれていく。その間の周波数はどうなっているのか? 例えば音叉の440Hzのところを見ると、396.9000000000001と441.0000000000001の間になってしまう。音叉のピークは440Hzなので、そこに近い441Hzにそれなりのピークはあるが、完全なものではない。396.9Hzにも少し大きな数値が入っている。数値を見る限り失われているようだ。でも440Hzにピークがありそうだという予測はできる。様々な補間方法を使うことで、点のつながりから予測することが可能。そこの計算をまじめにやれば、440Hzにピークがあることは突き止められるようだ。この辺りは窓関数との関わりも大きく、窓の種類で大きな差になる。Hammingだとやや末広がりな印象。Hanningは結構ピークを立ててくれる。いくつか試したらガウス窓などはサンプル数が少なくても、かなりの精度が出ている。目的にあわせて窓も使い分けるとよさそうだ。補間方法はこれから勉強してみようかと思う。

先週DFTの逆のiDFT、逆離散フーリエ変換を作っていたので、いくつか変更を加える程度で簡単にDFTは作ることができた。時間がかかったのは正しく動いているかどうかのチェックの方。でもDFTはスピードの問題を感じる。サンプル数を1万以上にすると、かなり待たされる。FFTが必要な理由を体感できた。近々FFTも実装しようと思う。

下は自作DFTを使って、結果をテキストファイルに出力。そのファイルをLibreOfficeのCalcに読み込んで、グラフ表示したもの。Audacityの表示では隠れてしまうようなデータも、当然ながら表示することができる。ただこの手の作業ではCalcは使い勝手がイマイチなので、表示ソフトもそのうち自作しようと思う。



DFT(discrete Fourier transform)
iDFT(inverse discrete Fourier transform)
C言語ソースコード

DFT(離散フーリエ変換)、iDFT(逆離散フーリエ変換)のC言語によるプログラム。動作説明用の簡易的なもので実際に音を加工するためのプログラムではない。式を変形せずに、そのまま入れているので、分かりやすいとは思う。実用プログラムでは無駄な処理を省きもっとシンプルになる。
入力サンプルを任意の数入れておくと、まずDFTし、実部と虚部、絶対値を計算する。さらにその数値を元に、iDFTで元の入力信号に戻している。振幅スペクトルについては、資料もなく勝手な解釈なので参考にしないこと。

/* dft.c
DFT iDFT 実験用プログラム */
#include <stdio.h>
#include <math.h>
/***********************************************************/
int main(void){
/* 入力信号 実部 */
  double sr[] ={
    1,
    1,
    1,
    1,
    0,
    0,
    0,
    0
    };

  /* 入力信号 虚部 ここではすべて0を入れておく。
     sizeofは1byteなのでdoubleを得るには8で割る */
  int N = sizeof(sr)/8, i, j;
  double si[N];
  for(i=0; i<sizeof(sr)/8; i++){
     si[i]=0;
  }
/***********************************************************/
/* フーリエ変換後の数値を入れる配列を3つ用意 */
  double hr[N];/* 実部 */
  double hi[N];/* 虚部 */
  double p[N];/* 振幅スペクトル */
/***********************************************************/
/* DFT処理 フーリエ変換 */
 double t = 2*M_PI/N;
 for(i=0; i<N; i++){ 
   hr[i] = 0;
   hi[i] = 0;
   for(j=0; j<N; j++){ 
     hr[i] += (sr[j]*cos(t*i*j))+(si[j]*sin(t*i*j));/* 実部 */
     hi[i] += (si[j]*cos(t*i*j))-(sr[j]*sin(t*i*j));/* 虚部 */
   }
 }
 /* 振幅を求めつつ結果を表示 */
 printf("\nN:\tDFT Real\t\tDFT Imaginary\t\tDFT Power\n");
 for(i=0; i<N; i++){
     /* 振幅スペクトル */
     if(i==0){
       p[i] = sqrt(hr[i]*hr[i] + hi[i]*hi[i])/N;
     }else{
       p[i] = sqrt(hr[i]*hr[i] + hi[i]*hi[i])/N*2;
     }
     /* 実部 虚部 振幅スペクトル 表示 */
     printf("%d:\t%19.18f\t%19.18f\t%19.18f\n"
                         ,i,hr[i],hi[i],p[i]);
 }
/***********************************************************/
/* iDFT処理 逆フーリエ変換 */
 for(i=0; i<N; i++){
   sr[i] = 0;/*初期化*/
   si[i] = 0;/*初期化*/
   for(j=0; j<N; j++){ 
     sr[i] += (hr[j]*cos(t*i*j))-(hi[j]*sin(t*i*j));/* 実部 */
     si[i] += (hi[j]*cos(t*i*j))+(hr[j]*sin(t*i*j));/* 虚部 */
   }
   sr[i] /= N;
   si[i] /= N;
 }
 printf("\nN:\tiDFT Real\t\tiDFT Imaginary\n");
 /* 信号を戻す 精度はdouble */
 for(i=0; i<N; i++){
     printf("%d:\t%19.18f\t%19.18f\n",i,sr[i],si[i]);
 }
/***********************************************************/
 return 0;
}

上記のままコンパイルして実行すると以下のような内容が表示される。精度がdoubleなので15桁ぐらいまで有効。ただ1.0を入力すると0.999999のようになってしまうが、丸めたり、1 > x > -1 であれば問題ないはず。またナイキスト周波数は妙な値になるので気にしない。
N: DFT Real               DFT Imaginary          DFT Power
0: 4.000000000000000000   0.000000000000000000   0.500000000000000000
1: 1.000000000000000000   -2.414213562373094923  0.653281482438188288
2: -0.000000000000000184  -0.000000000000000222  0.000000000000000072
3: 1.000000000000000000   -0.414213562373094923  0.270598050073098451
4: 0.000000000000000000   -0.000000000000000245  0.000000000000000061
5: 0.999999999999999223   0.414213562373095923   0.270598050073098395
6: 0.000000000000000329   -0.000000000000000333  0.000000000000000117
7: 0.999999999999999667   2.414213562373095368   0.653281482438188288

N: iDFT Real              iDFT Imaginary
0: 0.999999999999999889   0.000000000000000056
1: 1.000000000000000222   -0.000000000000000025
2: 0.999999999999999889   -0.000000000000000143
3: 1.000000000000000000   0.000000000000000083
4: -0.000000000000000078  -0.000000000000000059
5: -0.000000000000000593  0.000000000000000485
6: -0.000000000000000339  -0.000000000000000176
7: 0.000000000000001255   0.000000000000000748

サイン波のサンプル毎の値を出すJavaScriptを下記に作ってみた。サンプル数N(64以内ぐらい)と周期f(Nにもよるけど1,2,3,4ぐらい)を入れてボタンを押すとサンプル毎の数値が表示される。それを上記プログラムの配列部にコピーして、いろいろなパターンを実行してみればDFTの雰囲気が分かると思う。
N:
f:
sound programming 目次はこちら