C#实现快速傅里叶变换（FFT）

1、FFT类

using System;
using System.Collections.Generic;
using System.Linq;
using System.Runtime.InteropServices.ComTypes;
using System.Text;
using System.Threading.Tasks;

namespace DFT_FFTApp.Utils
{
    public class FFT
    {
        /// <summary>
        /// FFT
        /// </summary>
        /// <param name="xL"></param>
        /// <returns></returns>
        public static double[] ButterflyFFT(List<double> xL)
        {
            int len = xL.Count;
            if (len == 0)
            {
                return null;
            }
            double[] data_r = new double[len];
            double[] data_i = new double[len];
            for (int i = 0; i < len; i++)
            {
                data_r[i] = xL[i];
            }
            double[] fft_r = new double[data_r.Length];
            double[] fft_i = new double[data_r.Length];
            double[] result = new double[data_r.Length];
            BF_FFT(ref data_r, ref data_i, ref fft_r, ref fft_i);
            GetMod(ref fft_r, ref fft_i, ref result);
            return result;
        }

        /// <summary>
        /// FFT算法
        /// </summary>
        /// <param name="data_r"></param>
        /// <param name="data_i"></param>
        /// <param name="result_r"></param>
        /// <param name="result_i"></param>
        public static void BF_FFT(ref double[] data_r, ref double[] data_i, ref double[] result_r, ref double[] result_i)
        {
            if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length)
                return;
            int len = data_r.Length;
            double[] X_r = new double[len];
            double[] X_i = new double[len];
            for (int i = 0; i < len; i++)//将源数据复制副本，避免影响源数据的安全性
            {
                X_r[i] = data_r[i];
                X_i[i] = data_i[i];
            }
            DataSort(ref X_r, ref X_i);//位置重排
            double WN_r, WN_i;//旋转因子
            int M = (int)(Math.Log(len) / Math.Log(2));//蝶形图级数
            for (int l = 0; l < M; l++)
            {
                int space = (int)Math.Pow(2, l);
                int num = space;//旋转因子个数
                double temp1_r, temp1_i, temp2_r, temp2_i;
                for (int i = 0; i < num; i++)
                {
                    int p = (int)Math.Pow(2, M - 1 - l);//同一旋转因子有p个蝶
                    WN_r = Math.Cos(2 * Math.PI / len * p * i);
                    WN_i = -Math.Sin(2 * Math.PI / len * p * i);
                    for (int j = 0, n = i; j < p; j++, n += (int)Math.Pow(2, l + 1))
                    {
                        temp1_r = X_r[n];
                        temp1_i = X_i[n];
                        temp2_r = X_r[n + space];
                        temp2_i = X_i[n + space];//为蝶形的两个输入数据作副本，对副本进行计算，避免数据被修改后参加下一次计算
                        X_r[n] = temp1_r + temp2_r * WN_r - temp2_i * WN_i;
                        X_i[n] = temp1_i + temp2_i * WN_r + temp2_r * WN_i;
                        X_r[n + space] = temp1_r - temp2_r * WN_r + temp2_i * WN_i;
                        X_i[n + space] = temp1_i - temp2_i * WN_r - temp2_r * WN_i;
                    }
                }
            }
            result_r = X_r;
            result_i = X_i;
        }

        /// <summary>
        /// 对原数据组进行重排
        /// </summary>
        /// <param name="data_r"></param>
        /// <param name="data_i"></param>
        public static void DataSort(ref double[] data_r, ref double[] data_i)
        {
            if (data_r.Length == 0 || data_i.Length == 0 || data_r.Length != data_i.Length)
                return;
            int len = data_r.Length;
            int[] count = new int[len];
            int M = (int)(Math.Log(len) / Math.Log(2));
            double[] temp_r = new double[len];
            double[] temp_i = new double[len];

            for (int i = 0; i < len; i++)
            {
                temp_r[i] = data_r[i];
                temp_i[i] = data_i[i];
            }
            for (int l = 0; l < M; l++)
            {
                int space = (int)Math.Pow(2, l);
                int add = (int)Math.Pow(2, M - l - 1);
                for (int i = 0; i < len; i++)
                {
                    if ((i / space) % 2 != 0)
                        count[i] += add;
                }
            }
            for (int i = 0; i < len; i++)
            {
                data_r[i] = temp_r[count[i]];
                data_i[i] = temp_i[count[i]];
            }
        }

        /// <summary>
        /// 求模
        /// </summary>
        /// <param name="complex_r"></param>
        /// <param name="complex_i"></param>
        /// <param name="mod"></param>
        public static void GetMod(ref double[] complex_r, ref double[] complex_i, ref double[] mod)
        {
            if (complex_r.Length == 0 || complex_i.Length == 0 || complex_r.Length != complex_i.Length)
                return;
            for (int i = 0; i < complex_r.Length; i++)
                mod[i] = Math.Sqrt(complex_r[i] * complex_r[i] + complex_i[i] * complex_i[i]);
        }
    }
}

2、应用

private void DFTInitialize()
        {
            List<double> list = GetCaseData(128);
            double[] res = FFT.ButterflyFFT(list);
            foreach (double item in res)
            {
                Console.WriteLine(item);
            }
            easyChartX1.Plot(res, 0, 1);
        }
        
        /// <summary>
        /// 获取时域数据
        /// </summary>
        /// <param name="points"></param>
        /// <returns></returns>
        private List<double> GetCaseData(int points)
        {
            List<double> data = new List<double>();
            double w = 2*Math.PI / points;
            for(int i = 0; i < points;i++)
            {
                double d=Math.Sin(w * i);
                //Console.WriteLine($"i={i},d={d.ToString("G2")}");
                data.Add(d);
            }
            return data;
        }