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