Рассчитать значения спектра с помощью БПФ - PullRequest
Новая
       14

Рассчитать значения спектра с помощью БПФ

0 голосов
/ 29 апреля 2019

Я должен рассчитать значения спектра звука.Я использовал aForge FFT в Sources / Math / FourierTransform.cs, и я использовал пример выборки из 16 образцов, как в этом видео , чтобы проверить результаты с помощью Excel (я тестировалрезультаты в электронной таблице, как в видео).

БПФ: 

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

Результаты не ввсе похоже.Где это не так?Это неправильная реализация complex (Data) или метод FFT или другой?

Заранее спасибо!

1 Ответ

1 голос
/ 29 апреля 2019

Во-первых, в результате БПФ представляет собой комплексную функцию в целом.Вы отображаете только реальные части в своем коде, но то, с чем вы сравниваете, отображает величины , поэтому, конечно, они будут другими: вы сравниваете яблоки с апельсинами.

Когда вы используете амплитуды и сравниваете яблоки с яблоками, вы должны получить следующее: 

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.

...