【笔记】Cooley–Tukey FFT Algorithm - Iterative Edition

DFT

IDFT

FFT

1 当输入向量长度为2^k时

参考《Introduction to Algorithms》chapter 30.3

每次迭代的元素赋值及辅助方法reorder

/// <summary>
/// Radix-2 Step Helper Method
/// </summary>
/// <param name="data">Sample vector.</param>
/// <param name="exponentSign">Fourier series exponent sign.</param>
/// <param name="levelSize">Level Group Size.</param>
/// <param name="k">Index inside of the level.</param>
private static void Radix2Step(Complex[] data, int exponentSign, int levelSize, int k)
{
    // Twiddle Factor
    var exponent = (exponentSign * k) * Math.PI / levelSize;
    var w = new Complex(Math.Cos(exponent), Math.Sin(exponent));

    var step = levelSize << 1;
    for (var i = k; i < data.Length; i += step)
    {
        var ai = data[i];
        var t = w * data[i + levelSize];
        data[i] = ai + t;
        data[i + levelSize] = ai - t;
    }
}

/// <summary>
/// Radix-2 Reorder Helper Method
/// </summary>
private static void Radix2Reorder<T>(T[] data)
{
    var j = 0;
    for (var i = 0; i < data.Length - 1; i++)
    {
        if (i < j)
        {
            var temp = data[i];
            data[i] = data[j];
            data[j] = temp;
        }

        var m = data.Length;

        do
        {
            m >>= 1;
            j ^= m;
        }
        while ((j & m) == 0);
    }
}

当exponentSign = -1时计算的是DFT, 1时计算的是IDFT

/// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
private static void Radix2Forward(Complex[] data)
{
    Radix2Reorder(data);
    for (var levelSize = 1; levelSize < data.Length; levelSize <<= 1)
    {
        for (var k = 0; k < levelSize; k++)
        {
            Radix2Step(data, -1, levelSize, k);
        }
    }
}

/// <summary>
/// Radix-2 generic FFT for power-of-two sized sample vectors.
/// </summary>
private static void Radix2Inverse(Complex[] data)
{
    Radix2Reorder(data);
    for (var levelSize = 1; levelSize < data.Length; levelSize <<= 1)
    {
        for (var k = 0; k < levelSize; k++)
        {
            Radix2Step(data, 1, levelSize, k);
        }
    }
}

2 当输入向量长度!=2^k时

计算exp(I*Pi*n^2/N)

/// <summary>
/// Sequences with length greater than Math.Sqrt(Int32.MaxValue) + 1
/// will cause k*k in the Bluestein sequence to overflow (GH-286).
/// </summary>
private static readonly int BluesteinSequenceLengthThreshold = 46341;
/// <summary>
/// Generate the bluestein sequence for the provided problem size.
/// </summary>
/// <returns>Bluestein sequence exp(I*Pi*k^2/N)</returns>
private static Complex[] BluesteinSequence(int n)
{
    var s = Math.PI / n;
    var sequence = new Complex[n];

    // TODO: benchmark whether the second variation is significantly
    // faster than the former one. If not just use the former one always.
    if (n > BluesteinSequenceLengthThreshold)
    {
        for (var k = 0; k < sequence.Length; k++)
        {
            var t = (s * k) * k;
            sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
        }
    }
    else
    {
        for (var k = 0; k < sequence.Length; k++)
        {
            var t = s * (k * k);
            sequence[k] = new Complex(Math.Cos(t), Math.Sin(t));
        }
    }

    return sequence;
}

时域的卷积就是频域的乘积，所以卷积计算就是IDFT(DFT(a)*DFT(b))

/// <summary>
/// Convolution with the bluestein sequence.
/// </summary>
private static void BluesteinConvolution(Complex[] data)
{
    var n = data.Length;
    var sequence = BluesteinSequence(n);

    // Padding to power of two >= 2N–1 so we can apply Radix-2 FFT.
    var m = (int)CeilingPowerOf2(((uint)n << 1) - 1);
    var b = new Complex[m];
    var a = new Complex[m];

    // Build and transform padded sequence b_k = exp(I*Pi*k^2/N)
    for (var i = 0; i < n; i++)
    {
        b[i] = sequence[i];
    }

    for (var i = m - n + 1; i < b.Length; i++)
    {
        b[i] = sequence[m - i];
    }

    Radix2Forward(b);


    // Build and transform padded sequence a_k = x_k * exp(-I*Pi*k^2/N)
    for (var i = 0; i < data.Length; i++)
    {
        a[i] = Complex.Conjugate(sequence[i]) * data[i];
    }

    Radix2Forward(a);

    for (var i = 0; i < a.Length; i++)
    {
        a[i] *= b[i];
    }

    Radix2Inverse(a);

    var nbinv = 1.0 / m;
    for (var i = 0; i < data.Length; i++)
    {
        data[i] = nbinv * Complex.Conjugate(sequence[i]) * a[i];
    }
}

/// <summary>
/// least power of 2 greater than or equal to <paramref name="n"/>
/// </summary>
public static uint CeilingPowerOf2(uint n)
{
    --n;
    n |= (n >> 1);
    n |= (n >> 2);
    n |= (n >> 4);
    n |= (n >> 8);
    n |= (n >> 16);

    return n + 1;
}

IDFT先虚实交换再调用上述方法，最后再对结果虚实交换。

/// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary>
private static void BluesteinForward(Complex[] data)
    => BluesteinConvolution(data);

/// <summary>
/// Bluestein generic FFT for arbitrary sized sample vectors.
/// </summary>
private static void BluesteinInverse(Complex[] data)
{
    SwapRealImaginary(data);
    BluesteinConvolution(data);
    SwapRealImaginary(data);
}

/// <summary>
/// Swap the real and imaginary parts of each sample.
/// </summary>
private static void SwapRealImaginary(Complex[] data)
{
    for (var i = 0; i < data.Length; i++)
    {
        data[i] = new Complex(data[i].Imaginary, data[i].Real);
    }
}

结合上诉两种情况，得出最终DFT和IDFT代码

/// <summary>
/// Cooley–Tukey FFT algorithm,
/// </summary>
public static void DFT(Complex[] data)
{
    var len = (uint)data.Length;
    if (IsPowerOf2(len))
    {
        Radix2Forward(data);
    }
    else
    {
        BluesteinForward(data);
    }
}

/// <summary>
/// Cooley–Tukey FFT algorithm,
/// </summary>
public static void IDFT(Complex[] data)
{
    var len = (uint)data.Length;
    if (IsPowerOf2(len))
    {
        Radix2Inverse(data);
    }
    else
    {
        BluesteinInverse(data);
    }

    var factor = 1.0 / len;
    for (var i = 0; i < len; i++)
    {
        data[i] *= factor;
    }
}

/// <summary>
/// return true if n is power of 2
/// </summary>
public static bool IsPowerOf2(uint n)
    => n < 1u ? false : (n & (n - 1u)) == 0;