C#中的快速傅里叶变换(FFT)实现

79
在哪里可以找到一份免费、非常快速、可靠的C# FFT实现?这个实现可以用于产品中吗?或者有没有任何限制?
9个回答

55
那个做AForge的家伙做得相当不错,但它不是商业质量。它很适合学习,但你可以看出他也在学习,所以他有一些非常严重的错误,比如假设图像的大小而不是使用正确的每像素位数。
我不是在抨击那个家伙,我非常尊重他学习所有这些并向我们展示如何做到这一点。我认为他现在是博士或者至少快要成为博士了,所以他非常聪明,只是这个库不能商业使用。
Math.Net库在处理傅里叶变换和复杂图像/数字时有自己的怪异之处。比如,如果我没记错的话,它会以人类可视格式输出傅里叶变换,这对于人类来说很好,如果你想查看变换的图片,但当你期望数据以某种格式(正常格式)呈现时,它就不太好了。我可能记错了,但我记得有一些奇怪的地方,所以我实际上去了他们用于傅里叶变换的原始代码,效果要好得多。(ExocortexDSP v1.2 http://www.exocortex.org/dsp/

Math.net 在处理来自 FFT 的数据时还有其他一些我不喜欢的瑕疵。我不记得那是什么了,我只知道使用 ExoCortex DSP 库更容易获得我想要的结果。虽然我不是数学家或工程师,但对于这些人来说,可能很合理。

因此!我使用从 ExoCortex 获取的 FFT 代码,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


3
我很感兴趣了解Math.NET Iridium fft实现中的怪异之处的更多细节,这样我们就可以修复它!它是否与复数的处理方式有关?不过我不明白你提到的“人类可视格式”是什么意思。示例可以在http://mathnet.opensourcedotnet.info/doc/IridiumFFT.ashx找到。 - Christoph Rüegg
7
fftw存在一种有问题的许可证;请看这个描述:“FFTW也提供了一些非自由许可证,允许不同于GPL的使用条款。” - Daniel Mošmondor
这是给Mike Bethany的一个问题。我正在尝试学习如何将数据从时域转换为频率域。您提供的exocortex链接是否是正确的方法? - T o n y
Exo Cortext在.NET 4上抛出了System Out of Range异常,但没有附加信息。不起作用。 - bh_earth0

38

AForge.net 是一个带有快速傅里叶变换支持的免费(开源)库。(使用情况请参见 Sources/Imaging/ComplexImage.cs,实现方式请参见 Sources/Math/FourierTransform.cs


Note: I have translated the content into simplified Chinese.

AForge.net链接已更改: http://www.aforgenet.com/ - SiL3NC3
这个函数的输出与Matlab的输出相匹配,而且它也很轻巧。 - Gray Programmerz

15

Math.NET的Iridium库提供了一个快速、定期更新的数学相关函数集合,包括FFT算法。它基于LGPL许可证,因此您可以在商业产品中免费使用。


2
Math.NET Iridium非常适合将使用Apache commons-math的Java代码转换为.NET,因为两者之间的类和方法非常相似。95%的时间,您只需要更改类和方法名称,一切都可以正常工作。 - finnw

13

我看到这是一个旧帖子,但是无论如何,这是一份我在2010年编写的免费(MIT许可证)1-D 2的幂长度的C# FFT实现。

我还没有将其性能与其他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;
    }


由于链接无法访问,此答案现在完全无用。 - InDieTasten
非常抱歉。我在几年前删除了我的博客,因为它吸引了太多的垃圾邮件。不幸的是,代码有点太大了,无法在这里放入评论中。请给我发送一封电子邮件至g.<mysurname>@ieee.org,我很乐意将代码发送给您。 - Gerry Beauregard
你可以更新你的答案,添加代码并删除无效链接。通过私人渠道分享你的代码将违背 Stack Overflow 的精神。 - InDieTasten
5
完成。我几天前试图将其放在评论中,但太大了。我没有意识到评论的大小限制与答案不同。 - Gerry Beauregard

7

6

5

1
仅限于少数几种变换尺寸。 - J D

2
如果您不介意手动输入,Numerical Recipes网站(http://www.nr.com/)提供了一个FFT。我正在将一个Labview程序转换为C# 2008,.NET 3.5以获取数据并查看频谱。不幸的是,Math.Net使用了最新的.NET框架,因此我无法使用该FFT。我尝试了Exocortex,虽然它能够正常工作,但结果与Labview的结果不匹配,我不知道是什么原因导致了这个问题。所以我尝试了Numerical Recipes网站上的FFT,它可以正常工作!我还能够编写Labview低旁瓣窗口(并且必须引入缩放因子)。
您可以在其网站上以访客身份阅读Numerical Recipes书的章节,但该书非常有用,我强烈建议购买它。即使您最终使用Math.NET FFT。

使用 Numerical Recipes 中的任何代码时要小心。代码没有问题,但许可证是个问题。您必须付费才能使用该代码,非商业或科学应用也不例外。请参阅此链接以获取更多信息:http://mingus.as.arizona.edu/~bjw/software/boycottnr.html。 - Bob Bryan

0

如果要针对英特尔处理器进行多线程实现,我建议使用英特尔的MKL库。虽然它不是免费的,但价格实惠(不到100美元),速度非常快 - 但您需要通过P/Invokes调用其C dll。 Exocortex项目在6年前停止了开发,因此如果这是一个重要的项目,我建议小心使用。


3
2013年6月的单用户价格为499美元。 - RickNZ
截至2015年10月,Composer版本售价为699美元。 - mcy
“社区许可计划”不是免费的吗?Intel集成性能原语库(IPP)的无费选项,自助支持,免版税 - Lati

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接