使用 FFT 计算频谱值

Question

我必须计算音频的频谱值。 我在 Sources/Math/FourierTransform.cs 中使用了aForge 的 FFT，并且我使用了一个示例，使用了这个video 中使用的 16 个样本来检查结果（我在电子表格中测试了结果，就像在视频）。

FFT：

public enum Direction
{
    Forward = 1,
    Backward = -1
};

private const int minLength = 2;
private const int maxLength = 16384;
private const int minBits = 1;
private const int maxBits = 14;
private static int[][] reversedBits = new int[maxBits][];
private static Complex[,][] complexRotation = new Complex[maxBits, 2][];

static void Main(string[] args)
{
    var Data = new Complex[16];
    Data[0] = new Complex(0, 0);
    Data[1] = new Complex((float)0.998027, 0);
    Data[2] = new Complex((float)0.125333, 0);
    Data[3] = new Complex((float)-0.98229, 0);
    Data[4] = new Complex((float)-0.24869, 0);
    Data[5] = new Complex((float)0.951057, 0);
    Data[6] = new Complex((float)0.368125, 0);
    Data[7] = new Complex((float)-0.90483, 0);
    Data[8] = new Complex((float)-0.48175, 0);
    Data[9] = new Complex((float)0.844328, 0);
    Data[10] = new Complex((float)0.587785, 0);
    Data[11] = new Complex((float)-0.77051, 0);
    Data[12] = new Complex((float)-0.68455, 0);
    Data[13] = new Complex((float)0.684547, 0);
    Data[14] = new Complex((float)0.770513, 0);
    Data[15] = new Complex((float)-0.58779, 0);

    FFT(Data, Direction.Forward);

    for (int a = 0; a <= Data.Length - 1; a++)
    {
        Console.WriteLine(Data[a].Re.ToString());
    }

    Console.ReadLine();
}

public static void FFT(Complex[] data, Direction direction)
{
    int n = data.Length;
    int m = Tools.Log2(n);

    // reorder data first
    ReorderData(data);

    // compute FFT
    int tn = 1, tm;

    for (int k = 1; k <= m; k++)
    {
        Complex[] rotation = GetComplexRotation(k, direction);

        tm = tn;
        tn <<= 1;

        for (int i = 0; i < tm; i++)
        {
            Complex t = rotation[i];

            for (int even = i; even < n; even += tn)
            {
                int odd = even + tm;
                Complex ce = data[even];
                Complex co = data[odd];

                double tr = co.Re * t.Re - co.Im * t.Im;
                double ti = co.Re * t.Im + co.Im * t.Re;

                data[even].Re += tr;
                data[even].Im += ti;

                data[odd].Re = ce.Re - tr;
                data[odd].Im = ce.Im - ti;
            }
        }
    }

    if (direction == Direction.Forward)
    {
        for (int i = 0; i < n; i++)
        {
            data[i].Re /= (double)n;
            data[i].Im /= (double)n;
        }
    }
}

private static int[] GetReversedBits(int numberOfBits)
{
    if ((numberOfBits < minBits) || (numberOfBits > maxBits))
        throw new ArgumentOutOfRangeException();

    // check if the array is already calculated
    if (reversedBits[numberOfBits - 1] == null)
    {
        int n = Tools.Pow2(numberOfBits);
        int[] rBits = new int[n];

        // calculate the array
        for (int i = 0; i < n; i++)
        {
            int oldBits = i;
            int newBits = 0;

            for (int j = 0; j < numberOfBits; j++)
            {
                newBits = (newBits << 1) | (oldBits & 1);
                oldBits = (oldBits >> 1);
            }
            rBits[i] = newBits;
        }
        reversedBits[numberOfBits - 1] = rBits;
    }
    return reversedBits[numberOfBits - 1];
}

private static Complex[] GetComplexRotation(int numberOfBits, Direction direction)
{
    int directionIndex = (direction == Direction.Forward) ? 0 : 1;

    // check if the array is already calculated
    if (complexRotation[numberOfBits - 1, directionIndex] == null)
    {
        int n = 1 << (numberOfBits - 1);
        double uR = 1.0;
        double uI = 0.0;
        double angle = System.Math.PI / n * (int)direction;
        double wR = System.Math.Cos(angle);
        double wI = System.Math.Sin(angle);
        double t;
        Complex[] rotation = new Complex[n];

        for (int i = 0; i < n; i++)
        {
            rotation[i] = new Complex(uR, uI);
            t = uR * wI + uI * wR;
            uR = uR * wR - uI * wI;
            uI = t;
        }

        complexRotation[numberOfBits - 1, directionIndex] = rotation;
    }
    return complexRotation[numberOfBits - 1, directionIndex];
}

// Reorder data for FFT using
private static void ReorderData(Complex[] data)
{
    int len = data.Length;

    // check data length
    if ((len < minLength) || (len > maxLength) || (!Tools.IsPowerOf2(len)))
            throw new ArgumentException("Incorrect data length.");

    int[] rBits = GetReversedBits(Tools.Log2(len));

    for (int i = 0; i < len; i++)
    {
        int s = rBits[i];

        if (s > i)
        {
            Complex t = data[i];
            data[i] = data[s];
            data[s] = t;
        }
    }
}

这些是转换后的结果：

Output FFT results:             Excel FFT results:
0,0418315622955561              0,669305
0,0533257974328085              0,716163407
0,137615673627316               0,908647001         
0,114642731070279               1,673453043
0,234673940537634               7,474988602
0,0811255020953362              0,880988382          
0,138088891589122               0,406276784
0,0623766891658306              0,248854492
0,0272978749126196              0,204227
0,0124250144575261              0,248854492
0,053787064184711               0,406276784
0,00783331226557493             0,880988382
0,0884368745610118              7,474988602
0,0155431246384978              1,673453043
0,0301093757152557              0,908647001
0                               0,716163407

结果完全不相似。哪里错了？ 是复杂（Data）的实现错误还是FFT方法错误或其他？

提前致谢！

Answer 1

首先，生成的 FFT 通常是一个复数函数。您只显示代码中的真实部分，但您要比较的是显示幅度，所以它们当然会有所不同：您是在比较苹果和橘子。

当您使用大小并将苹果与苹果进行比较时，您应该会得到：

for (int a = 0; a <= Data.Length - 1; a++)
{
    Console.WriteLine(Data[a].Magnitude.ToString());
}

...

0.0418315622955561
0.0447602132472683
0.0567904388057513
0.104590813761862
0.46718679147454
0.0550617784710375
0.025392294285886
0.0155534081359397
0.0127641875296831
0.0155534081359397
0.025392294285886
0.0550617784710375
0.46718679147454
0.104590813761862
0.0567904388057513
0.0447602132472683

这看起来好一些——它具有与 Excel 输出相同的对称属性，并且在相同的位置似乎有峰值。

看起来天平好像关闭了。如果我将每个元素除以 Excel 输出中的相应元素，我会得到：

16
16
16
16
16
16
16
16
16
16
16
16
16
16
16
16

所以你的结果是非常正确的，只是一个比例因子。

您在FFT 的最后一步将所有内容除以n：

if (direction == Direction.Forward)
{
    for (int i = 0; i < n; i++)
    {
        data[i].Re /= (double)n;
        data[i].Im /= (double)n;
    }
}

这通常是为inverse transform 完成的，而不是前向变换。

总而言之，将输出从Data[a].Re 更改为Data[a].Magnitude 并将FFT 末尾的条件从if (direction == Direction.Forward) 更改为if (direction == Direction.Backward)，我得到了这个输出：

0.669304996728897
0.716163411956293
0.908647020892022
1.67345302018979
7.47498866359264
0.880988455536601
0.406276708574176
0.248854530175035
0.20422700047493
0.248854530175035
0.406276708574176
0.880988455536601
7.47498866359264
1.67345302018979
0.908647020892022
0.716163411956293

与 Excel 输出相匹配。