Есть ли ошибки в реализации fft?

Question

[eliasum]

Реализован расчет спектра сигнала через fft, файл fft.cs:

namespace FFTW
{
    using NAudio.Wave;
    using System;
    using System.Numerics;

    internal static class fft
    {
        internal static byte[] OpenWAVFile(string file)
        {
            using (WaveFileReader reader = new WaveFileReader(file))
            {
                byte[] buffer = new byte[reader.Length];
                reader.Read(buffer, 0, buffer.Length);
                return buffer;
            }
        }

        internal static Complex[] Convert(int[] value)
        {
            Complex[] buffer = new Complex[value.Length];   
            for (int i = 0; i < value.Length; i++)
            {
                // Re = value[i], Im = 0
                buffer[i] = new Complex(value[i], 0);
                buffer[i] *= 1;                            
            }
            return buffer;
        }

        internal static class FFT_V1
        {
            private static Complex w(int k, int n)
            {
                if (k % n == 0) return 1;
                double arg = -2 * Math.PI * k / n;
                return new Complex(Math.Cos(arg), Math.Sin(arg));
            }

            public static Complex[] Calculate(Complex[] value)
            {
                // Check if it is splitted enough
                if (value != null && value.Length <= 1) { return value; }

                int n = value.Length >> 1;  

                // Split even and odd
                Complex[] odd = new Complex[n];     
                Complex[] even = new Complex[n];    

                for (int i = 0; i < n; i++)
                {
                    even[i] = value[i * 2];       
                    odd[i] = value[i * 2 + 1];     
                }

                // Split on tasks
                even = Calculate(even); 
                odd = Calculate(odd);   

                // Calculate DFT
                for (int k = 0; k < n; k++)
                {
                    value[k] = even[k] + w(k, value.Length) * odd[k];
                    value[n + k] = even[k] - w(k, value.Length) * odd[k];
                }
                return value;
            }
        }
    }
}

Файл Program.cs:

using System.Numerics;

namespace FFTW
{
    internal class Program
    {
        private static void Main()
        {
            byte[] source = fft.OpenWAVFile(@".\source.wav");

            int[] buffer = new int[source.Length >> 1];

            for (int i = 0, k = 0; k < buffer.Length; i += 2, k++)
                buffer[k] = System.Convert.ToUInt16(source[i]) |
                    System.Convert.ToUInt16(source[i + 1] << 8);

            Complex[] spectrum1 = fft.FFT_V1.Calculate(fft.Convert(buffer));
        }
    }
}

Вопрос к знатокам алгоритмов, есть ли ошибки или недочёты в реализации fft?

Comments

[Василий Банников]

ну напиши тесты и узнаешь.

Answer 1

[Wataru]

Зачем buffer[i] *= 1;?

Есть ошибка с тем, что у вас код работает только с длиной равной степеням двойки, но нигде этого не проверяется. Обычно, чтобы не возиться с частыми случаями, расширяют буфер до степени двойки.

Если же делать с нечетными n, то эти стандартные формулы не работают.

Еще возможно по коду разбросаны всякие off by one error и перепутанные знаки (почему в w() аргумент у синуса/косинуса берется со знаком минус?)

Еще недочет - вы 2 раза вызываете w(), когда как можно было бы и запомнить результат первого вызова.

И еще, обычно рекурсию разворачивают в циклы. В английской вики приведен псевдокод.

Answer 2

[Армянское Радио]

Ну супир, синусоиду на вход давали? Что на выходе?

Недостатки - память перевыделяется на динамике при каждом вызове, взрослые библиотеки так не делают (они перекладывают на пользователя) - это влияет на общую скорость работы.