C# Csharp调用FFTW进行傅里叶变换
为了对一个星期以来工作做一个总结，也可能为大家提供一些参考
完整代码可以从官网下载：链接
我也把代码放到百度网盘了：链接 提取码：0h1x
视频教程：链接（自己熬夜录的，一看就懂）
代码部分：
首先看一下这个目录结构：
一共有3个代码文件，我们主要关注的是FFTWtest.cs这个文件，因为interop.cs实际上是对依赖文件的一个引入，Managed.cs是一些方法的封装
以下是FFTWtest.cs的源码（含我的标注），配合视频教程食用口味更佳
using System;
using System.Runtime.InteropServices;
using FFTWSharp;
namespace FFTWSharp_test
{
public class FFTWtest
{
const int sampleSize = 16384;
const int repeatPlan = 10000;
const bool forgetWisdom = false; //Wisdom 说明http://www.fftw.org/fftw-wisdom.1.html
//pointers to unmanaged arrays 指向数组的指针
IntPtr pin, pout; //这个现在理解成指针
//managed arrays 管理数组
float[] fin, fout;
double[] din, dout;
//handles to managed arrays, keeps them pinned in memory
//托管数组的句柄，把他们留在数组中 ？这个GCHandle怎么理解
//为了防止托管代码把对象回收，需要把那个对象“钉”住(pin)，让它的内存地址固定，
//而不被垃圾回收掉，然后最后我们自己管理，自己释放内存,这时候就需要GCHandle
GCHandle hin, hout, hdin, hdout;
//pointers to the FFTW plan objects FFTW plan 对象的指针
IntPtr fplan1, fplan2, fplan3, fplan4, fplan5;
// and an example of the managed interface
fftwf_plan mplan1, mplan2, mplan3, mplan4; //plan对象型变量
fftwf_complexarray mfin, mfout; //对象型变量
fftw_plan mplan5;
fftw_complexarray mdin, mdout; //复数数组
// length of arrays
int fftLength = 0; //fft的长度
// Initializes FFTW and all arrays 初始化所有的fftw和所有的数组
// n: Logical size of the transform
public FFTWtest(int n) //构造函数
{
Console.WriteLine($"Start testing with n = {n.ToString("#,0")} complex numbers. All plans will be executed {repeatPlan.ToString("#,0")} times on a single thread.");
Console.WriteLine("Please wait, creating plans...");
fftLength = n;
// create two unmanaged arrays, properly aligned
//pin和pout都是IntPtr类型的
pin = fftwf.malloc(n * 8); //创建空间，用pin指向空间的首地址
pout = fftwf.malloc(n * 8);
// create two managed arrays, possibly misalinged
// n*2 because we are dealing with complex numbers
fin = new float[n * 2]; //给数组分配存储空间，float类型
fout = new float[n * 2];
// and two more for double FFTW
din = new double[n * 2]; //给数组分配存储空间，float类型
dout = new double[n * 2];
// get handles and pin arrays so the GC doesn't move them
//为了防止托管代码把对象回收，需要把那个对象“钉”住(pin)，让它的内存地址固定，
//而不被垃圾回收掉，然后最后我们自己管理，自己释放内存,这时候就需要GCHandle
//C#的垃圾处理机制的触发条件：
//当处理newobj(cil 的new指令)指令时，
//如果CLR判定托管堆没有足够的空间来分配所请求的类型,它会执行一次垃圾回收来尝试释放内存
//意思是fin可能会被自动回收，但是hin不会回收
// 类型说明：GCHandle hin, hout, hdin, hdout;
hin = GCHandle.Alloc(fin, GCHandleType.Pinned); //意思是fin可能会被自动回收，但是hin不会回收
hout = GCHandle.Alloc(fout, GCHandleType.Pinned);
hdin = GCHandle.Alloc(din, GCHandleType.Pinned);
hdout = GCHandle.Alloc(dout, GCHandleType.Pinned);
// create a few test transforms
/// <dft_1d summary>
/// 为一维的复数到复数的离散傅里叶变换创建一个plan
/// Creates a plan for a 1-dimensional complex-to-complex DFT
/// </summary>
/// <param name="n">The logical size of the transform</param> 逻辑空间大小
/// <param name="direction">Specifies the direction of the transform</param>指定转换的方向
/// <param name="input">Pointer to an array of 8-byte complex numbers</param>指向复试数组的指针
/// <param name="output">Pointer to an array of 8-byte complex numbers</param>
/// <param name="flags">Flags that specify the behavior of the planner</param> 指定计划器行为的标志
///创建各个plan用来做性能比较
//pin->pout
//pin和pout都是IntPtr类型
//这为什么用的pin，不怕对象被回收么？这里的pin为什么能直接用了
//回收的是fin？因为fin是float[]类型？
fplan1 = fftwf.dft_1d(n, pin, pout, fftw_direction.Forward, fftw_flags.Estimate);
fplan2 = fftwf.dft_1d(n, hin.AddrOfPinnedObject(), hout.AddrOfPinnedObject(), fftw_direction.Forward, fftw_flags.Estimate);
fplan3 = fftwf.dft_1d(n, hout.AddrOfPinnedObject(), pin, fftw_direction.Backward, fftw_flags.Measure);
// end with transforming back to original array 以转换回原始数组结束
fplan4 = fftwf.dft_1d(n, hout.AddrOfPinnedObject(), hin.AddrOfPinnedObject(), fftw_direction.Backward, fftw_flags.Estimate);
// and check a quick one with doubles, just to be sure
fplan5 = fftw.dft_1d(n, hdin.AddrOfPinnedObject(), hdout.AddrOfPinnedObject(), fftw_direction.Backward, fftw_flags.Measure);
// create a managed plan as well 同时创建一个管理计划
mfin = new fftwf_complexarray(fin); // 从浮点数组创建与FFTW兼容的数组，仅初始化为单精度
mfout = new fftwf_complexarray(fout);
mdin = new fftw_complexarray(din);
mdout = new fftw_complexarray(dout);
Complex<->Complex transforms
mplan1 = fftwf_plan.dft_1d(n, mfin, mfout, fftw_direction.Forward, fftw_flags.Estimate);
mplan2 = fftwf_plan.dft_1d(n, mfin, mfout, fftw_direction.Forward, fftw_flags.Measure);
mplan3 = fftwf_plan.dft_1d(n, mfin, mfout, fftw_direction.Forward, fftw_flags.Patient);
mplan4 = fftwf_plan.dft_1d(n, mfout, mfin, fftw_direction.Backward, fftw_flags.Measure);
mplan5 = fftw_plan.dft_1d(n, mdin, mdout, fftw_direction.Forward, fftw_flags.Measure);
// fill our arrays with an arbitrary complex sawtooth-like signal
//用任意复杂的锯齿状信号填充阵列
for (int i = 0; i < n * 2; i++)
fin[i] = i % 50;
for (int i = 0; i < n * 2; i++)
fout[i] = i % 50;
for (int i = 0; i < n * 2; i++)
din[i] = i % 50;
//for(int i = 0; i < n * 2; i++)
//{
// Console.WriteLine(fin[i]);
//}
// copy managed arrays to unmanaged arrays
//将托管数组复制到非托管数组即本地执行
// 参数:
// source:
// 从中进行复制的一维数组。
//
// startIndex:
// 源数组中从零开始的索引，在此处开始复制。
//
// destination:
// 要复制到的内存指针。
//
// length:
// 要复制的数组元素的数目,从fin复制到pin。
Marshal.Copy(fin, 0, pin, n * 2);
Marshal.Copy(fout, 0, pout, n * 2);
Console.WriteLine();
}
public void TestAll()
{
TestPlan(fplan1, "single (malloc) | pin, pout, Forward, Estimate");
TestPlan(fplan2, "single (pinned) | hin, hout, Forward, Estimate");
TestPlan(fplan3, "single (pinned) | hout, pin, Backward, Measure ");
// set fin to 0, and try to refill it from a backwards fft from fout (aka hin/hout)
for (int i = 0; i < fftLength * 2; i++) fin[i] = 0;
TestPlan(fplan4, "single (pinned) | hout, hin, Backward, Estimate");
// check and see how we did, don't say anyt
for (int i = 0; i < fftLength * 2; i++)
{
// check against original values
// note that we need to scale down by length, due to FFTW scaling by N
if (System.Math.Abs(fin[i]/fftLength - (i % 50)) > 1e-3)
{
System.Console.WriteLine("FFTW consistency error!");
return;
}
}
TestPlan(fplan5, "double (pinned) | hdin, hdout, Backward, Measure ");
Console.WriteLine();
TestPlan(mplan1, "#1 single (managed) | mfin, mfout, Forward, Estimate");
TestPlan(mplan2, "#2 single (managed) | mfin, mfout, Forward, Measure ");
TestPlan(mplan3, "#3 single (managed) | mfin, mfout, Forward, Patient ");
Console.WriteLine();
// fill our input array with an arbitrary complex sawtooth-like signal
for (int i = 0; i < fftLength * 2; i++) fin[i] = i % 50;
for (int i = 0; i < fftLength * 2; i++) fout[i] = 0;
mfin.SetData(fin);
mfout.SetData(fout);
TestPlan(mplan2, "#2 single (managed) | mfin, mfout, Forward, Measure ");
fout = mfout.GetData_Float(); // let's see what's in mfout
// at this point mfout contains the FFT'd mfin
TestPlan(mplan4, "#4 single (managed) | mfout, mfin, Backward, Measure ");
// at this point we have transfarred backwards the mfout into mfin, so mfin should be very close to the original complex sawtooth-like signal
fin = mfin.GetData_Float();
for (int i = 0; i < fftLength * 2; i++) fin[i] /= fftLength;
// at this point fin should be very close to the original (sawtooth-like) signal
// check and see how we did, don't say anyt
for (int i = 0; i < fftLength * 2; i++)
{
// check against original values
if (System.Math.Abs(fin[i] - (i % 50)) > 1e-3)
{
System.Console.WriteLine("FFTW consistency error!");
return;
}
}
Console.WriteLine();
TestPlan(mplan2, "#2 single (managed) | mfin, mfout, Forward, Measure ");
TestPlan(mplan5, "#5 double (managed) | mdin, mdout, Forward, Measure ");
Console.WriteLine();
}
// Tests a single plan, displaying results
// plan: Pointer to plan to test
public void TestPlan(object plan, string planName)
// fplan1, "single (malloc) | pin, pout, Forward, Estimate
{
// a: adds, b: muls, c: fmas
double a = 0, b = 0, c = 0;
int start = System.Environment.TickCount; //plan的开始时间
if (plan is IntPtr)
{
IntPtr umplan = (IntPtr)plan; //这里为什么要强转一下
for (int i = 0; i < repeatPlan; i++)
{
fftwf.execute(umplan); //为什么不能把IntPtr直接传进去
Console.WriteLine(Marshal.ReadIntPtr(pin));
//Console.WriteLine();
//Console.WriteLine(Marshal.ReadIntPtr(pout));
}
// Console.WriteLine(Marshal.ReadIntPtr(pin));
//Console.WriteLine(Marshal.ReadIntPtr(pout));
/// <summary>
/// Returns (approximately) the number of flops used by a certain plan
/// </summary>
/// <param name="plan">The plan to measure</param>
/// <param name="add">Reference to double to hold number of adds</param> 引用双精度保存加法数
/// <param name="mul">Reference to double to hold number of muls</param>
/// <param name="fma">Reference to double to hold number of fmas (fused multiply-add)</param>
/// <remarks>Total flops ~= add+mul+2*fma or add+mul+fma if fma is supported</remarks>
fftwf.flops(umplan, ref a, ref b, ref c);
}
if (plan is fftw_plan)
{
fftw_plan mplan = (fftw_plan)plan;
for (int i = 0; i < repeatPlan; i++)
{
mplan.Execute();
}
fftw.flops(mplan.Handle, ref a, ref b, ref c);
}
if (plan is fftwf_plan) //这个if跟上边的if有什么区别
{
fftwf_plan mplan = (fftwf_plan)plan;
for (int i = 0; i < repeatPlan; i++)
{
mplan.Execute();
}
fftwf.flops(mplan.Handle, ref a, ref b, ref c);
}
double mflops = (((a + b + 2 * c)) * repeatPlan) / (1024 * 1024); //这个评价指标是浮点运算次数，但是为什么这么算
long ticks = (System.Environment.TickCount - start); //plan所花费的时间
//Console.WriteLine($"Plan '{planName}': {ticks.ToString("#,0")} us | mflops: {FormatNumber(mflops)} | mflops/s: {(1000*mflops/ticks).ToString("#,0.0")}");
Console.WriteLine("Plan '{0}': {1,8:N0} us | mflops: {2,8:N0} | mflops/s: {3,8:N0}", planName, ticks, mflops, (1000 * mflops / ticks));
//Console.WriteLine();
///pin是IntPtr类型的，理解为指向数组的指针，这样的数组怎么遍历
/*
for(int i = 0; i < 5; i++)
{
// Console.WriteLine(pin);
//这里应该用Marshal的某个方法读取IntPtr的数，类似Marshal.readInt16(Pin),是不是这个方法还得查
//Marshal.ReadIntPtr(pin);
Console.WriteLine(Marshal.ReadIntPtr(pin));
Console.WriteLine(Marshal.ReadIntPtr(pout));
}
*/
//Console.WriteLine();
//Console.WriteLine(pout);
}
// Releases all memory used by FFTW/C# 析构函数
~FFTWtest()
{
fftwf.export_wisdom_to_filename("wisdom.wsd");
// it is essential that you call these after finishing
// that's why they're in the destructor. See Also: RAII
fftwf.free(pin);
fftwf.free(pout);
fftwf.destroy_plan(fplan1);
fftwf.destroy_plan(fplan2);
fftwf.destroy_plan(fplan3);
fftwf.destroy_plan(fplan4);
fftwf.destroy_plan(fplan5);
hin.Free();
hout.Free();
}
static void Main(string[] args)
{
if (forgetWisdom)
{
fftwf.fftwf_forget_wisdom();
}
else
{
Console.WriteLine("Importing wisdom (wisdom speeds up the plan creation process, if that plan was previously created at least once)");
fftwf.import_wisdom_from_filename("wisdom.wsd");
}
// initialize our FFTW test class
FFTWtest test = new FFTWtest(sampleSize); //调用有参构造函数
// run the tests, print debug output
test.TestAll();
// pause for user input, then quit
System.Console.WriteLine("\nDone. Press any key to exit.");
String str = System.Console.ReadLine();
}
}
}
扫一扫
1395
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
