Реализация быстрого преобразования Фурье (БПФ) в C#

AForge.net - это бесплатная библиотека (с открытым исходным кодом) с поддержкой быстрого преобразования Фурье. (См. Sources / Imaging / ComplexImage.cs для использования, Sources / Math / FourierTransform.cs для реализации)

http://www.exocortex.org/dsp/ - это математическая библиотека C# с открытым исходным кодом с алгоритмами БПФ.

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

Парень, который сделал AForge, проделал неплохую работу, но это не коммерческое качество. На этом приятно учиться, но вы можете сказать, что он тоже учился, поэтому у него есть довольно серьезные ошибки, такие как предположение размера изображения вместо использования правильных бит на пиксель.

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

Библиотека Math.Net имеет свои странности при работе с преобразованиями Фурье и сложными изображениями / числами. Например, если я не ошибаюсь, он выводит преобразование Фурье в формате, доступном для просмотра человеком, что хорошо для людей, если вы хотите посмотреть на изображение преобразования, но это не так хорошо, когда вы ожидаете, что данные будут в определенном формате. формат (нормальный формат). Я мог ошибаться в этом, но я просто помню, что была некоторая странность, поэтому я фактически перешел к исходному коду, который они использовали для материала Фурье, и он работал намного лучше. (ExocortexDSP v1.2 http://www.exocortex.org/dsp/)

У Math.net были и другие причуды, которые мне не нравились при работе с данными из БПФ, я не могу вспомнить, что это было, я просто знаю, что было намного проще получить то, что я хотел, из библиотеки ExoCortex DSP. Но я не математик или инженер; для этих парней это могло иметь смысл.

Так! Я использую код БПФ, взятый из ExoCortex, на котором основан Math.Net, без чего-либо еще, и он отлично работает.

И, наконец, я знаю, что это не C#, но я начал рассматривать использование FFTW (http://www.fftw.org/). И этот парень уже сделал оболочку C#, так что я собирался ее проверить, но на самом деле еще не использовал. (http://www.sdss.jhu.edu/~tamas/bytes/fftwcsharp.html)

ОЙ! Я не знаю, делаете ли вы это для учебы или работы, но в любом случае есть ОТЛИЧНАЯ серия бесплатных лекций, которую читает профессор Стэнфордского университета в iTunes University.

https://podcasts.apple.com/us/podcast/the-fourier-transforms-and-its-applications/id384232849

Для многопоточной реализации, настроенной для процессоров Intel, я бы ознакомился с библиотекой Intel MKL. Это не бесплатно, но доступно (менее 100 долларов) и невероятно быстро, но вам нужно будет называть это C dll через P / Invokes. Проект Exocortex прекратил разработку 6 лет назад, поэтому я буду осторожен при его использовании, если это важный проект.

Вот еще один; порт C# для Ooura FFT. Это достаточно быстро. Пакет также включает свертку перекрытия / добавления и некоторые другие вещи DSP под лицензией MIT.

https://github.com/hughpyle/inguz-DSPUtil/blob/master/Fourier.cs

Я вижу, что это старый поток, но, чего бы то ни было, вот бесплатная (лицензия MIT) реализация C# FFT 1-D, основанная только на мощности двух, которую я написал в 2010 году.

Я не сравнивал его производительность с другими реализациями C# FFT. Я написал его в основном для сравнения производительности Flash / ActionScript и Silverlight / C#. Последнее намного быстрее, по крайней мере, для обработки чисел.

/**
 * Performs an in-place complex FFT.
 *
 * Released under the MIT License
 *
 * Copyright (c) 2010 Gerald T. Beauregard
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 */
public class FFT2
{
    // Element for linked list in which we store the
    // input/output data. We use a linked list because
    // for sequential access it's faster than array index.
    class FFTElement
    {
        public double re = 0.0;     // Real component
        public double im = 0.0;     // Imaginary component
        public FFTElement next;     // Next element in linked list
        public uint revTgt;         // Target position post bit-reversal
    }

    private uint m_logN = 0;        // log2 of FFT size
    private uint m_N = 0;           // FFT size
    private FFTElement[] m_X;       // Vector of linked list elements

    /**
     *
     */
    public FFT2()
    {
    }

    /**
     * Initialize class to perform FFT of specified size.
     *
     * @param   logN    Log2 of FFT length. e.g. for 512 pt FFT, logN = 9.
     */
    public void init(
        uint logN )
    {
        m_logN = logN;
        m_N = (uint)(1 << (int)m_logN);

        // Allocate elements for linked list of complex numbers.
        m_X = new FFTElement[m_N];
        for (uint k = 0; k < m_N; k++)
            m_X[k] = new FFTElement();

        // Set up "next" pointers.
        for (uint k = 0; k < m_N-1; k++)
            m_X[k].next = m_X[k+1];

        // Specify target for bit reversal re-ordering.
        for (uint k = 0; k < m_N; k++ )
            m_X[k].revTgt = BitReverse(k,logN);
    }

    /**
     * Performs in-place complex FFT.
     *
     * @param   xRe     Real part of input/output
     * @param   xIm     Imaginary part of input/output
     * @param   inverse If true, do an inverse FFT
     */
    public void run(
        double[] xRe,
        double[] xIm,
        bool inverse = false )
    {
        uint numFlies = m_N >> 1;   // Number of butterflies per sub-FFT
        uint span = m_N >> 1;       // Width of the butterfly
        uint spacing = m_N;         // Distance between start of sub-FFTs
        uint wIndexStep = 1;        // Increment for twiddle table index

        // Copy data into linked complex number objects
        // If it's an IFFT, we divide by N while we're at it
        FFTElement x = m_X[0];
        uint k = 0;
        double scale = inverse ? 1.0/m_N : 1.0;
        while (x != null)
        {
            x.re = scale*xRe[k];
            x.im = scale*xIm[k];
            x = x.next;
            k++;
        }

        // For each stage of the FFT
        for (uint stage = 0; stage < m_logN; stage++)
        {
            // Compute a multiplier factor for the "twiddle factors".
            // The twiddle factors are complex unit vectors spaced at
            // regular angular intervals. The angle by which the twiddle
            // factor advances depends on the FFT stage. In many FFT
            // implementations the twiddle factors are cached, but because
            // array lookup is relatively slow in C#, it's just
            // as fast to compute them on the fly.
            double wAngleInc = wIndexStep * 2.0*Math.PI/m_N;
            if (inverse == false)
                wAngleInc *= -1;
            double wMulRe = Math.Cos(wAngleInc);
            double wMulIm = Math.Sin(wAngleInc);

            for (uint start = 0; start < m_N; start += spacing)
            {
                FFTElement xTop = m_X[start];
                FFTElement xBot = m_X[start+span];

                double wRe = 1.0;
                double wIm = 0.0;

                // For each butterfly in this stage
                for (uint flyCount = 0; flyCount < numFlies; ++flyCount)
                {
                    // Get the top & bottom values
                    double xTopRe = xTop.re;
                    double xTopIm = xTop.im;
                    double xBotRe = xBot.re;
                    double xBotIm = xBot.im;

                    // Top branch of butterfly has addition
                    xTop.re = xTopRe + xBotRe;
                    xTop.im = xTopIm + xBotIm;

                    // Bottom branch of butterly has subtraction,
                    // followed by multiplication by twiddle factor
                    xBotRe = xTopRe - xBotRe;
                    xBotIm = xTopIm - xBotIm;
                    xBot.re = xBotRe*wRe - xBotIm*wIm;
                    xBot.im = xBotRe*wIm + xBotIm*wRe;

                    // Advance butterfly to next top & bottom positions
                    xTop = xTop.next;
                    xBot = xBot.next;

                    // Update the twiddle factor, via complex multiply
                    // by unit vector with the appropriate angle
                    // (wRe + j wIm) = (wRe + j wIm) x (wMulRe + j wMulIm)
                    double tRe = wRe;
                    wRe = wRe*wMulRe - wIm*wMulIm;
                    wIm = tRe*wMulIm + wIm*wMulRe;
                }
            }

            numFlies >>= 1;     // Divide by 2 by right shift
            span >>= 1;
            spacing >>= 1;
            wIndexStep <<= 1;   // Multiply by 2 by left shift
        }

        // The algorithm leaves the result in a scrambled order.
        // Unscramble while copying values from the complex
        // linked list elements back to the input/output vectors.
        x = m_X[0];
        while (x != null)
        {
            uint target = x.revTgt;
            xRe[target] = x.re;
            xIm[target] = x.im;
            x = x.next;
        }
    }

    /**
     * Do bit reversal of specified number of places of an int
     * For example, 1101 bit-reversed is 1011
     *
     * @param   x       Number to be bit-reverse.
     * @param   numBits Number of bits in the number.
     */
    private uint BitReverse(
        uint x,
        uint numBits)
    {
        uint y = 0;
        for (uint i = 0; i < numBits; i++)
        {
            y <<= 1;
            y |= x & 0x0001;
            x >>= 1;
        }
        return y;
    }

}

На веб-сайте Numerical Recipes (http://www.nr.com/) есть БПФ, если вы не против его ввести. Я работаю над проектом по преобразованию программы Labview в C# 2008, .NET 3.5 для сбора данных и затем посмотрите на частотный спектр. К сожалению, Math.Net использует последнюю версию .NET framework, поэтому я не мог использовать этот БПФ. Я попробовал Exocortex - он работал, но результаты совпадают с результатами Labview, и я недостаточно знаю теорию БПФ, чтобы понять, что вызывает проблему. Итак, я попробовал БПФ на веб-сайте числовых рецептов, и он сработал! Я также смог запрограммировать окно с низким боковым лепестком Labview (и мне пришлось ввести коэффициент масштабирования).

Вы можете прочитать главу книги «Числовые рецепты» в качестве гостя на этом сайте, но книга настолько полезна, что я настоятельно рекомендую ее приобрести. Даже если вы в конечном итоге используете Math.NET FFT.

Старый вопрос, но он все еще появляется в результатах Google ...

Очень неограниченную лицензионную библиотеку C# / .NET MIT можно найти по адресу:

https://www.codeproject.com/articles/1107480/dsplib-fft-dft-fourier-transform-library-for-net

Эта библиотека быстрая, поскольку она выполняет параллельные потоки на нескольких ядрах, и она очень полная и готова к использованию.