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

離散フーリエ変換 その2

離散フーリエ変換の第2回目という事で、実際にどういうことをやっているのかをコードで見ていきたいと思います。前回のDFTのコードの変換を行っているところが以下の部分です。

// iは抜き出す周波数
for (int i = 0; i < n; i++) {
    
    tmpReal[i] = 0.0;
    tmpImag[i] = 0.0;
    
    // 1サンプルの位相の差分
    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);
    }
}

変換するオーディオデータに対して、抜き出す周波数のcosをかけ算して全て足し合わせた値が実部、-sinをかけて足し合わせた値が虚部になります。具体的に説明はしませんが、このようにする事で、それぞれの周波数のcosやsinがオーディオデータにどれくらい含まれているかを得る事ができます。(ちゃんと理解したい方は、なにかしらフーリエ変換の本などを参考にしてください)

抜き出す周波数というのは、0Hzから始まって1Hz、2Hz...と続いて、変換するサンプル数-1Hzまでです。つまり、変換するサンプル数と同じ数の周波数成分に分けられます。DFTは、変換するオーディオデータが延々とループしていると仮定して周波数成分を取り出しますので、1Hz以上の周波数はぴったり整数の周波数になります。

DFT02ri.jpg

上の図を見ると、cosもsinも0Hzの時は横一直線になっています。0HzはDC成分といわれたりしてちょっと特別です。cos0Hzでは全体に1がかけられますので、波形全体がプラスとマイナスのどちらにバランスが偏っているかという要素になります。sinは0をかけているので、元がどんな波形だろうと結果は0になります。

ちなみに、下の図の1.23Hzみたいな中途半端な周波数のサイン波などというのは途中で急激に途切れた状態で繰り返されることになりますので、DFTをすると1Hzを中心に高い周波数まで全体的に成分が含まれることがわかります。

DFT02s123.jpg

また、ちょっと説明が後回しになりましたが、なぜDFTでcosとsinを抜き出すかといえば、同じ周波数のcosとsinを足し合わせることで、360度どの位相のサイン波でも表現できるからです。

本当にcosとsinの足し算でサイン波が作れるのか、検証してみたのが以下のコードです。stRadを0〜2πの間で変えてみて実行してみてください。sinはDFTと同じようにマイナスにしています。もちろん計算の誤差はあると思いますので、近ければOKという基準で判断すれば、どの位相のサイン波でも、スタート位置がゼロのcosとsinの足し合わせで作り出せる事がわかると思います。(てきとうに書いたので、なんかおかしかったらすみません。0〜2πの範囲外では正確な結果が出るようになってません)いちおうDFTと同じように、sinはマイナスsinでやってます。

#include <stdio.h>
#include <math.h>
int main (int argc, const char * argv[]) {
    
    double stRad = 0.0; //sinの開始位置。0〜2πの間で!
    int n = 16;
    
    double xcos = 0;
    double xsin = 0;
    double co = fabs(sin(stRad));
    double si = sqrt(1.0 - co * co);
    
    if (stRad < M_PI_2) {
        //第一象限
        xcos = co;
        xsin = -si;
    } else if (stRad < M_PI) {
        //第二象限
        xcos = co;
        xsin = si;
    } else if (stRad < M_PI_2 * 3) {
        //第三象限
        xcos = -co;
        xsin = si;
    } else {
        //第四象限
        xcos = -co;
        xsin = -si;
    }
    
    for (int i = 0; i < n; i++) {
        
        double phase = (double)i / n * 2.0 * M_PI;
        double sinVal = sin(stRad + phase);
        double mixVal = cos(phase) * xcos - sin(phase) * xsin;
        
        printf("sin %03d %f / sin+cos %f\n", i, sinVal, mixVal);
        if (fabs(sinVal - mixVal) > 0.000001) {
            printf("error\n");
            return 0;
        }
    }
    
    printf("success\n");
    
    return 0;
}