快速傅立葉變換（FFT）的C#代碼

這個代碼是從《快速傅立葉變換（FFT）的C++實現與Matlab實驗》這篇文章裏的源代碼轉換而來，請注意查看原文。

在這裏自己轉換成了C#代碼，並作了一些改動，主要是對N值的確定，原文的N值為常量1024，自己通過對輸入的數組的長度來確定N值，N值的確定符合2的n次方，函數返回N值。通過作者提供的測試變量進行測試，能得到相應的結果。代碼如下：

FFT代碼：

using System;
using System.Collections.Generic;
using System.Text;

namespace ConsoleApplication1
{
    /// <summary>
    /// 快速傅立葉變換(Fast Fourier Transform)。
    /// </summary>
    public class TWFFT
    {
        private TWFFT()
        {
        }

        private static void bitrp(float[] xreal, float[] ximag, int n)
        {
            // 位反轉置換 Bit-reversal Permutation
            int i, j, a, b, p;

            for (i = 1, p = 0; i < n; i *= 2)
            {
                p++;
            }
            for (i = 0; i < n; i++)
            {
                a = i;
                b = 0;
                for (j = 0; j < p; j++)
                {
                    b = b * 2 + a % 2;
                    a = a / 2;
                }
                if (b > i)
                {
                    float t = xreal[i];
                    xreal[i] = xreal[b];
                    xreal[b] = t;

                    t = ximag[i];
                    ximag[i] = ximag[b];
                    ximag[b] = t;
                }
            }
        }

        public static int FFT(float[] xreal, float[] ximag)
        {
            //n值為2的N次方
            int n = 2;
            while (n <= xreal.Length)
            {
                n *= 2;
            }
            n /= 2;

            // 快速傅立葉變換，將複數 x 變換後仍保存在 x 中，xreal, ximag 分別是 x 的實部和虛部
            float[] wreal = new float[n / 2];
            float[] wimag = new float[n / 2];
            float treal, timag, ureal, uimag, arg;
            int m, k, j, t, index1, index2;

bitrp(xreal, ximag, n);

            // 計算 1 的前 n / 2 個 n 次方根的共軛複數 W'j = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
            arg = (float)(-2 * Math.PI / n);
            treal = (float)Math.Cos(arg);
            timag = (float)Math.Sin(arg);
            wreal[0] = 1.0f;
            wimag[0] = 0.0f;
            for (j = 1; j < n / 2; j++)
            {
                wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
                wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
            }

            for (m = 2; m <= n; m *= 2)
            {
                for (k = 0; k < n; k += m)
                {
                    for (j = 0; j < m / 2; j++)
                    {
                        index1 = k + j;
                        index2 = index1 + m / 2;
                        t = n * j / m;    // 旋轉因子 w 的實部在 wreal [] 中的下標為 t
                        treal = wreal[t] * xreal[index2] - wimag[t] * ximag[index2];
                        timag = wreal[t] * ximag[index2] + wimag[t] * xreal[index2];
                        ureal = xreal[index1];
                        uimag = ximag[index1];
                        xreal[index1] = ureal + treal;
                        ximag[index1] = uimag + timag;
                        xreal[index2] = ureal - treal;
                        ximag[index2] = uimag - timag;
                    }
                }
            }

return n;
}

        public static int IFFT(float[] xreal, float[] ximag)
        {
            //n值為2的N次方
            int n = 2;
            while (n <= xreal.Length)
            {
                n *= 2;
            }
            n /= 2;

            // 快速傅立葉逆變換
            float[] wreal = new float[n / 2];
            float[] wimag = new float[n / 2];
            float treal, timag, ureal, uimag, arg;
            int m, k, j, t, index1, index2;

bitrp(xreal, ximag, n);

            // 計算 1 的前 n / 2 個 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, ... , n / 2 - 1
            arg = (float)(2 * Math.PI / n);
            treal = (float)(Math.Cos(arg));
            timag = (float)(Math.Sin(arg));
            wreal[0] = 1.0f;
            wimag[0] = 0.0f;
            for (j = 1; j < n / 2; j++)
            {
                wreal[j] = wreal[j - 1] * treal - wimag[j - 1] * timag;
                wimag[j] = wreal[j - 1] * timag + wimag[j - 1] * treal;
            }

            for (j = 0; j < n; j++)
            {
                xreal[j] /= n;
                ximag[j] /= n;
            }

            return n;
        }
    }
}

這裏隻給出了函數FFT()的，函數IFFT()可以進行相應測試。測試代碼：

using System;
using System.Collections.Generic;
using System.Text;

namespace ConsoleApplication1
{
    class Program
    {
        static void Main(string[] args)
        {
            float[] a = {
                0.5751f,0.4514f,0.0439f,0.0272f,0.3127f,0.0129f,0.3840f,0.6831f,
                0.0928f,0.0353f,0.6124f ,0.6085f,0.0158f,0.0164f,0.1901f,0.5869f};

            float[] b = {
                0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,
                0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f};

int n = TWFFT.FFT(a,b);

Console.WriteLine("FFT:");

Console.WriteLine("FFT的返回值N = {0}",n);
Console.WriteLine();

Console.WriteLine("i\ta\tb");

            Console.WriteLine();
            for (int i = 0; i < n; i++)
            {
                Console.WriteLine("{0}\t{1}\t{2}",i,a[i],b[i]);
            }

            Console.ReadLine();
        }
    }
}

最後更新：2017-04-02 15:15:00

快速傅立葉變換（FFT）的C#代碼

上一篇：黑客能否黑掉好奇號火星車？

下一篇：卡巴斯基躋身全球IT安全軟件廠商前三強

相關內容

熱門內容

最新內容

快速傅立葉變換（FFT）的C#代碼

上一篇： 黑客能否黑掉好奇號火星車？

下一篇： 卡巴斯基躋身全球IT安全軟件廠商前三強

相關內容

熱門內容

最新內容

上一篇：黑客能否黑掉好奇號火星車？

下一篇：卡巴斯基躋身全球IT安全軟件廠商前三強