An implementation of the fast Fourier transform (FFT) in C# [closed]

Question

Where can I find a free, very quick, and reliable implementation of FFT in C#?

That can be used in a product? Or are there any restrictions?

Kobylak

Answer 1

AForge.net is a free (open-source) library with Fast Fourier Transform support. (See Sources/Imaging/ComplexImage.cs for usage, Sources/Math/FourierTransform.cs for implemenation)

Needs

AForge.net link changed: aforgenet.com – Schonfield

This function output matches matlab output and also is lite weight. – Ogawa

downloaded AForge from code.google.com/archive/p/aforge/downloads, third link is a zip with the code. – Candice

Answer 2

The guy that did AForge did a fairly good job but it's not commercial quality. It's great to learn from but you can tell he was learning too so he has some pretty serious mistakes like assuming the size of an image instead of using the correct bits per pixel.

I'm not knocking the guy, I respect the heck out of him for learning all that and show us how to do it. I think he's a Ph.D now or at least he's about to be so he's really smart it's just not a commercially usable library.

The Math.Net library has its own weirdness when working with Fourier transforms and complex images/numbers. Like, if I'm not mistaken, it outputs the Fourier transform in human viewable format which is nice for humans if you want to look at a picture of the transform but it's not so good when you are expecting the data to be in a certain format (the normal format). I could be mistaken about that but I just remember there was some weirdness so I actually went to the original code they used for the Fourier stuff and it worked much better. (ExocortexDSP v1.2 http://www.exocortex.org/dsp/)

Math.net also had some other funkyness I didn't like when dealing with the data from the FFT, I can't remember what it was I just know it was much easier to get what I wanted out of the ExoCortex DSP library. I'm not a mathematician or engineer though; to those guys it might make perfect sense.

So! I use the FFT code yanked from ExoCortex, which Math.Net is based on, without anything else and it works great.

And finally, I know it's not C#, but I've started looking at using FFTW (http://www.fftw.org/). And this guy already made a C# wrapper so I was going to check it out but haven't actually used it yet. (http://www.sdss.jhu.edu/~tamas/bytes/fftwcsharp.html)

OH! I don't know if you are doing this for school or work but either way there is a GREAT free lecture series given by a Stanford professor on iTunes University.

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

Motorway

I'd be interested in more details about the weirdness in the Math.NET Iridium fft implementation - so we can fix it! ;). Is it related to how complex numbers are handled? No idea what you mean with the "human viewable format" though. Samples: mathnet.opensourcedotnet.info/doc/IridiumFFT.ashx – Arak

fftw has some kind of problematic license; check this out: "Non-free licenses for FFTW are also available that permit different terms of use than the GPL." – Rives

This is a question to Mike Bethany. I am trying to learn how to convert data from time domain to frequency domain. Is your exocortex link the correct way to do so? – Athallia

exo cortext throws system out of range exception without additional info on.net4 . not working. – Heppman

Answer 3

Math.NET's Iridium library provides a fast, regularly updated collection of math-related functions, including the FFT. It's licensed under the LGPL so you are free to use it in commercial products.

Briannebriano

+1. Math.NET Iridium is great for translating Java code (that uses Apache commons-math) to .NET thanks to the close correspondence between classes and methods of each. 95% of the time all you have to do is change class and method names and everything will work. – Fichu

FFT on Mathnet.Numerics needs native providers (MKL) which are problematic to install on linux. The useManaged() alternative to native acceleration throws unimplemented. – Candice

Answer 4

I see this is an old thread, but for what it's worth, here's a free (MIT License) 1-D power-of-2-length-only C# FFT implementation I wrote in 2010.

I haven't compared its performance to other C# FFT implementations. I wrote it mainly to compare the performance of Flash/ActionScript and Silverlight/C#. The latter is much faster, at least for number crunching.

/**
 * 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;
    }

}

Sharlenesharline

This answer is now completely useless due to the link leading nowhere... – Mazurka

Sorry about that. I removed my blog a couple of years ago as it was attracting too much spam. Unfortunately the code is a bit too big to put in a comment here. Drop me a note at g.<mysurname>@ieee.org, and I'd be happy to send the code to you. – Sharlenesharline

You can update your answer, add the code and remove your dead link. Sharing your code via private channels would be against the spirit of Stack Overflow. – Mazurka

Done. I tried to put it in a comment a couple of days ago, and but it was too big. It didn't occur to me that the size limit would be different for comments than for answers. – Sharlenesharline

Answer 5

An old question but it still shows up in Google results...

A very un-restrictive MIT Licensed C# / .NET library can be found at,

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

This library is fast as it parallel threads on multiple cores and is very complete and ready to use.

Edacity

Answer 6

Here's another; a C# port of the Ooura FFT. It's reasonably fast. The package also includes overlap/add convolution and some other DSP stuff, under the MIT license.

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

Blasphemous

Answer 7

http://www.exocortex.org/dsp/ is an open-source C# mathematics library with FFT algorithms.

Ultramontane

Limited to only a few transform sizes. – Cassady

Answer 8

The Numerical Recipes website (http://www.nr.com/) has an FFT if you don't mind typing it in. I am working on a project converting a Labview program to C# 2008, .NET 3.5 to acquire data and then look at the frequency spectrum. Unfortunately the Math.Net uses the latest .NET framework, so I couldn't use that FFT. I tried the Exocortex one - it worked but the results to match the Labview results and I don't know enough FFT theory to know what is causing the problem. So I tried the FFT on the numerical recipes website and it worked! I was also able to program the Labview low sidelobe window (and had to introduce a scaling factor).

You can read the chapter of the Numerical Recipes book as a guest on thier site, but the book is so useful that I highly recomend purchasing it. Even if you do end up using the Math.NET FFT.

Trepang

Be careful with any code that you use from Numerical Recipes. Nothing wrong with the code, it is the license that is the problem. You have to pay to use the code, and no exceptions for non-commercial or scientific applications. See this link for more info. – Haubergeon

Answer 9

For a multi-threaded implementation tuned for Intel processors I'd check out Intel's MKL library. It's not free, but it's afforable (less than $100) and blazing fast - but you'd need to call it's C dll's via P/Invokes. The Exocortex project stopped development 6 years ago, so I'd be careful using it if this is an important project.

Glooming

Single-user price as of June 2013 is US$499. – Lambent

As of October 2015, composer edition is $699 – Famagusta

Isn't it free with 'Community Licensing' program? No Cost Options for Intel Integrated Performance Primitives Library (IPP), Support Yourself, Royalty-Free – Degrading