月別アーカイブ: 2010年1月

離散フーリエ変換 その1

MacとかiPhoneとかあまり関係なく基礎を固めようと思っていろいろ勉強モードに入っていまして、すこしずつまとめていこうと思います。まずはフーリエ変換についてです。

あくまでプログラムで使う事を前提に書いていきますので、いろいろおかしかったりするかもしれませんがご了承ください。もし、明らかに変だったら突っ込んでいただけるとうれしいです。基本的に小難しそうな数式とかはできるだけ出さずにコード重視、ビジュアル重視で書いていくつもりです。

プログラミングでフーリエ変換というと、高速フーリエ変換(FFT)を使うという事になると思うのですが、FFTについては特に詳しく書きません。プログラム化されたコードを見た場合、FFTだと高速化されたアルゴリズムだけで本質的な部分がわからないので、高速化していないノーマルな離散フーリエ変換(DFT)のプログラムを見て、フーリエ変換の仕組みを調べていきます。

オーディオ信号に対してDFTを使うと何ができるかというと、時間単位で並んでいる時間領域のオーディオデータの一部分を切り出して、含まれている周波数ごとの成分に分けた周波数領域のデータに変換する事ができます。逆に、その周波数領域のデータを時間領域のオーディオデータに戻す事もでき、それが逆離散フーリエ変換(IDFT)になります。

DFT01.jpg

離散フーリエ変換というのは、デジタルでサンプリングされたオーディオ信号のような飛び飛びのデータに対して行うフーリエ変換です。いきなり「離散」と出てくると何の事かわからないかもしれませんが、「離散」=「デジタル」と考えればよいと思います。アナログの「連続」に対して、デジタルの「離散」です。そのDFTを高速化したのがFFTという関係になります。

まずは何より使ってみない事にはしょうがないということで、サンプルソースです。

#include <stdio.h>
#include <math.h>
void DFT(int n, double *real, double *imag)
{
    double tmpReal[n], tmpImag[n];
    
    for (int i = 0; i < n; i++) {
        
        tmpReal[i] = 0.0;
        tmpImag[i] = 0.0;
        
        double d = 2.0 * M_PI * i / n;
        
        for (int j = 0; j < n; j++) {
            
            double phase = d * j;
            
            tmpReal[i] += real[j] * cos(phase);
            tmpImag[i] -= real[j] * sin(phase);
        }
    }
    
    for (int i = 0; i < n; i++) {
        real[i] = tmpReal[i];
        imag[i] = tmpImag[i];
    }
}
int main (int argc, const char * argv[]) {
    
    int n = 16;
    double real[n], imag[n];
    
    double d = 2.0 * M_PI / n;
    
    for (int i = 0; i < n; i++) {
        real[i] = sin(1.0 * i * d); //1Hzのサイン波
        real[i] += sin(3.0 * i * d + M_PI_4); //3Hzのサイン波(1/4πずらし)
        real[i] += sin(5.0 * i * d + M_PI_2); //5Hzのサイン波(1/2πずらし)
        imag[i] = 0.0;
    }
    
    //フーリエ変換
    DFT(n, real, imag);
    
    for (int i = 0; i < n; i++) {
        printf("%dHz %f\n", i, sqrt(real[i] * real[i] + imag[i] * imag[i]));
    }
    
    return 0;
}

Xcodeで実行する場合はCのCommand Line Toolでプロジェクトを作って、main.cを書き換えてください。

今回のコードでは、1Hzと3Hzと5Hzのサイン波をいろいろ位相を変えてミックスして元のデータを作ってから変換しています。結果は、

0Hz 0.000000
1Hz 8.000000
2Hz 0.000000
3Hz 8.000000
4Hz 0.000000
5Hz 8.000000
6Hz 0.000000
7Hz 0.000000
8Hz 0.000000
9Hz 0.000000
10Hz 0.000000
11Hz 8.000000
12Hz 0.000000
13Hz 8.000000
14Hz 0.000000
15Hz 8.000000

と表示されて、1Hzと3Hzと5Hzが含まれているということがわかります。上の周波数の11Hzと13Hzと15Hzにもデータが現れていますが、これは、ナイキスト周波数以上の周波数の音は、ナイキスト周波数を対称に折り返した低い周波数へ現れるということで出てきているものですので、単純に周波数特性を調べるだけなら無視してください(けっして無意味なデータというわけではありませんが...)。

コードの説明に入りまして、DFT関数がまさに離散フーリエ変換する関数になります。realに変換したいオーディオデータの配列、imagに同じサイズの空の配列、nに配列のサイズをそれぞれ渡すと、realとimagにフーリエ変換された値が返ってきます。

DFTをすると1つの周波数ごとにデータが実部と虚部という2つに分かれて返ってきます。今回のコードではrealが実部、imagが虚部のデータです。実部と虚部は2つでひとつの意味を持ったデータですので、どちらかだけを使うことはあまり無いと思います。main関数の中でやっていますが、周波数特性の振幅を得るという場合はrealとimagをそれぞれ2乗した値の平方根を求めます。

DFTして得られる周波数は、DFTに渡すサンプル数の長さを1Hzとして、0からサンプル数-1までの整数倍の周波数です。16サンプルの場合は、0・1・2・3・4・5・6・7・8・9・10・11・12・13・14・15Hzの16個です。ただこれが1
秒を1Hzと考えるときには、たとえば44.1kHzのオーディオの中の64サンプルを変換したという状況だと、0・689・1378...というふうにもっと広い間隔の周波数の並びになります。これはどこを基準にするかで違ってくるというだけの話ですので、このDFTの説明ページでは変換するサンプル数を1Hzとして表現します。

とりあえず、DFTができましたというところで1回目は終わりです。次へつづきます。