傅里叶变化-Math.Net中Fourier变换使用详解

 最近公司做一个实验与仿真协同平台的项目，采用的是BS/CS的一种混合式架构方案。最初我提案用Electron技术实现客户端，后来经过几轮讨论限于人力将该方案枪决了。最终，客户端采用基于开源项目CSharpBrowser实现。但项目迭代到第二轮，产品想要加入一些计算资源框架，来解决实验数据处理中的高等运算问题。你懂的，小项目组的一贯原则是首选免费开源框架，在千挑万选之下，Math.Net渐渐浮出水面，虽然这套框架的API和Demo目前尚缺，但是不失为一套有力的数值计算开源框架，限于晚上对于该框架的Fourier介绍较少，这里仅做一点补充。（简书的排版能力还是不行啊→_→!!!）

1. Fourier变换简介

        网络博客中关于连续/离散Fourier变换的文章已经非常详实，本无需赘述。但毕竟下文要用到，所以这里还要简明扼要的说一下。
        简单说，Fourier变换就是将周期信号沿正交基分解，而一组良好的正交基就是正弦/余弦函数。不过套用伟大的欧拉公式后，我们更多是把

作为正交基。基于此连续域上的Fourier变换及其逆变换可以写为

不过，对于归一化参数可以略作调整，从而将Fourier变换对写为

但是，对于计算机是无法处理连续变量的，从而在上述工作基础之上发展了离散Fourier变换(DFT)，将其变换对写为

形式1

此即是Matlab做快速Fourier变换(FFT)的基础。但是同样的，我们可以对变换对的归一化参数，继而将DFT写为

形式2

以上只是扼要的介绍了Fourier变换。

2. Math.Net——C#开源数值计算框架

        说道数值计算，脑海中第一浮现的当然是Matlab/Octave这种高级语言，如果嫌弃它们安装麻烦，那么我们还有伟大的Python可用，就算再不济C/C++中也有很多相当不错的开源项目来支撑高效率的数值计算环境。所以C#当中的计算框架就显得弱鸡了一些。事实也证明，大多时候C#更适合客户端的表现和业务处理，鲜有用它去做计算类的东西。但是当真的遇到这种需求的时候，我其实还是更多的维护一个原则，那就是原生框架的效率优先于差异语言编译。
        C#当中Math.Net框架是一个相当不错开源工具包，但是相关的资料却不甚丰富，也缺乏深度。Math.Net能够支撑大部分数值计算处理，例如微分，积分，积分变换，线性代数计算，统计，信号处理，复变函数以及大规模稀疏矩阵的存储和并行等等问题。在语言方面能够兼顾C#和F#，同时支持一定程度的符号运算，并且符号运算的结果可以是Matlab展示形式，也可以直接Format成Latex样式。所以说Math.Net的综合能力还是令人为之一振的，敲起代码有一种Matlab/Mathematica/Maple任你摆布的感觉。

注意：如果Math.Net库版本与.netframework版本兼容异常需要统一版本否则程序会报错，在nuget包管理器中下载安装正确的版本库才可以用。

3. Math.Net——从复数(Complex)开始

        老实说，类似于C#/Java这些老牌业务型语言，在兼顾计算资源的时候，真的很难像Matlab等专业软件表现的那样优雅。例如创建一个复数

  Complex32 c = new Complex32(1f,2f);
  Console.WriteLine(c);//print (1,2)

如此可想而知，要创建一个复数序列是多麻烦的事情。但是在Matlab这种软件当中大可以优雅的多

>>c = 1+2j;

当然，优雅当不了饭吃，能用就是好滴。Math.Net中构造的复数，其实部和虚部只能是float类型，当然他也给出了许多复数的常规计算，例如

  Complex32 c = new Complex32(1f,2f);
  c.CommonLogarithm();\\取对数
  c.Conjugate();\\共轭
  var img = c.Imaginary;\\return 2
  var real = c.Real;\\return 1
  var magnitude = c.Magnitude;\\return 幅值
  var phase = c.Phase;\\return 相位

4. Fourier变换

        正常来说，Fourier变换是指对一个复数序列进行变换，如下例子

Complex32[] series = new Complex32[] 
{ 
  new Complex32(1, 2), 
  new Complex32(2, 1), 
  new Complex32(3, -2), 
  new Complex32(-1, 4) 
};
//print series  未经过FFT的samples序列
Fourier.Forward(series);//Fast Fourier变换
//Fourier.Forward(series, FourierOptions.Default);上一行等同于该行
//print series  已被FFT的spectrum序列
Fourier.Inverse(series);//Fast Fourier逆变换
//print series  已被逆变换的spectrum序列

Math.Net中采用的是内部运算，所以当执行Fourier变换之后，结果装入原有序列，其中值得关注的是，当FourierOptions采用Default值时，FFT遵循的是我们在第一部分讨论的形式2FFT，而如果想得到与Matlab运算相同的结果可以将FourierOptions修改为FourierOptions.Matlab，即

Fourier.Forward(series, FourierOptions.Matlab)

但是现实中，我们的实际信号采用往往是一个实数序列，而不是上述的复数序列。其实Forward方法还有很多其他的重载，这里我们不妨提供一个实现实序列Fourier变换的思路

double[] RealSamples = new double[100];//构造采样信号
double[] ImgSamples = new double[RealSamples .Length];//辅助虚部
Random r = new Random();
for(int i = 0; i<RealSamples.Length;i++)
{
  RealSamples[i] = r.NextDouble()*10-5;//随机噪声信号
  ImgSamples[i] = 0;
}
Fourier.Forward(RealSamples , ImgSamples, FourierOptions.Matlab );
double[] result = new double[RealSamples.Length];//结果序列
for(int i = 0; i<RealSamples.Length;i++)
{
  Complex32 temp = new Complex32(RealSamples[i],ImgSamples [i]);
  result[i] = temp.Magnitude();//从计算结果中检出结果序列
}

大家看到了，上面过程中使用了一个全部为0的虚数列来补足信号的采样序列，从而来完成Fourier变换过程。但是显然我们为这样一个API付出的太多了，有没有更简洁的办法呢，答案当然是有的

注意："虚数单位不可以等于零。
在复数a+bi中，a称为复数的实部，b称为复数的虚部，i称为虚数单位。当虚部等于零时，这个复数就是实数；当虚部不等于零时，这个复数称为虚数，虚数的实部如果等于零，则称为纯虚数。由上可知，复数集包含了实数集，因而是实数集的扩张。
在计算中规定i的平方=-1，并且i可以与实数在一起按照同样的运算律进行四则运算，i叫做虚数单位。
所以，虚数单位不等于零。"

int N = 512;
 Fourier.ForwardReal(NewRealSamples,N);
//这里NewRealSamples是对RealSamples重新构造的一个实序列

是不是瞬间简洁多了，不过这里面API对RealSamples要求是如果RealSamples.Length是偶数，那么NewRealSample.Length就要是RealSamples.Length+2的，而如果RealSamples.Length是奇数，那么NewRealSample.Length就等于RealSamples.Length+1；那么你可能会问为什么是这个样子的？另外，这里为什么无缘无故多出来一个整数N？其实这些问题就牵扯出了另外一个话题——DFT计算流，限于篇幅和主题，本人计划以后再去讨论。这里值得一提的是，FFT 的结果数据长度等于采样信号序列的长度，但由于FFT计算结果的对称性，其实真正的结果是结果数组中的前N/2或(N+1)/2个元素。