フーリエ変換

はじめに

簡単に言うと、ある波形データにどの周波数成分がどのくらいあるのかを解析するための手法で、
総当たりで調べる離散フーリエ変換と、高速に計算するための高速フーリエ変換がある。

離散フーリエ変換

public static Complex[] DFT(Complex[] input)
{
    int fs = 1800; // サンプリング周波数

    var output = new Complex[input.Length];
    for (int i = 0; i < input.Length; i++)
    {
        double re = 0, im = 0;
        for (int j = 0; j < input.Length; j++)
        {
            re += input[j].Real * Math.Cos(2 * Math.PI * i * j / fs);
            im += input[j].Real * -(Math.Sin(2 * Math.PI * i * j / fs));
        }
        output[i] = new Complex(re, im);
    }
    return output;
}

高速フーリエ変換

public static Complex[] FFT(Complex[] input, int bitSize)
{
    var dataSize = 1 << bitSize;
    var reverseBitArray = BitScrollArray(dataSize);

    var output = new Complex[dataSize];
    for (int i = 0; i < dataSize; i++)
    {
        // バタフライ演算のための置き換え
        output[i] = input[reverseBitArray[i]];
    }

    // バタフライ演算
    for (int stage = 1; stage <= bitSize; stage++)
    {
        var butterflyDistance = 1 << stage;
        var numType = butterflyDistance >> 1;
        var butterflySize = butterflyDistance >> 1;

        var w = Complex.One;
        var u = new Complex(Math.Cos(Math.PI / butterflySize), -Math.Sin(Math.PI / butterflySize));

        for (int type = 0; type < numType; type++)
        {
            for (int j = type; j < dataSize; j += butterflyDistance)
            {
                int jp = j + butterflySize;
                var temp = output[jp] * w;
                output[jp] -= temp;
                output[j] += temp;
            }
            w *= u;
        }
    }

    return output;
}

// ビットを左右反転した配列を返す
private static int[] BitScrollArray(int arraySize)
{
    var result = new int[arraySize];
    var halfSize = arraySize >> 1;

    result[0] = 0;
    for (int i = 1; i < arraySize; i <<= 1)
    {
        for (int j = 0; j < i; j++) result[j + i] = result[j] + halfSize;
        halfSize >>= 1;
    }
    return result;
}

窓関数

矩形窓