用FFT计算频谱值

Question

我必须计算音频的频谱值。我在Sources/Math/FourierTransform.cs中使用了aForge的快速傅立叶变换，并使用视频中使用的16个样本进行采样，用excel检查结果(我在像视频一样的电子表格中测试了结果)。

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

结果一点也不相似。哪里出问题了？复杂(数据)的实现是错误的，还是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

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

在您的n的最后一步，您将所有事情除以FFT

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

这通常用于逆变换，而不是正向转换。

总之，将输出从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输出相匹配。