解调算法还有不少有趣的，比如DSB、LSB、USB解调器，这些以后会讲，读者也可以先把这些模块的代码对照Matlab RTL-SDR那本书看一下。另外，比如DownConverter和Oscillator也是一个经典的用软件方式实现以前硬件中实现的工作的例子。还有DCRemover在我们大多数的SDR中也很重要，Decimator在数字信号处理中也很常见。RdsDecoder是国外的在FM信号中携带文字信息的解调器。这些以后有时间再详细讲，但我认为接下来最值得讲的是SDRSharp中的FFT实现。

FFT算法可以说是数字信号处理里最重要的算法，全称是快速傅立叶变换。傅立叶变换的种类有很多，我们在通信领域会接触到的还有一些比如离散傅立叶变换、连续傅立叶变换，其实如果不考虑公式推导，只是做定性理解的话，它们的功能都差不多，都是把一个信号从时域变换到了频域，用来分析信号的频域特征的。

通过傅立叶变换，就能知道给定时间内的信号的各频率分量，比如如果是音频信号（声波信号）做傅立叶变换，就能知道是高音分量大还是低音分量大，用来判断这个声音是高音还是低音。iPhone应用商店里有不少这种APP，对着手机唱“哆 唻 咪 发 嗖 拉 西”，手机经过FFT运算后可以得出主要分量对应的音阶，有些APP还能把频谱图显示出来，你可以看到除了主要分量外的其它谐波分量，你可以看看音调升高后，频谱图上的主要分量是否在向右侧（高频率方向）移动。

无线电信号是承载在电磁波上的，电磁波很多时候与声波有相似的特性，我们也利用FFT算法来分析它的各频率分量。FFT与最原始的傅立叶变换的区别就是可以在计算机上实现，且运算速度快。

SDRSharp里其实有两个地方都实现了FFT，一个是SDRSharp/Radio/Fourier.cs，另一个是SDRSharp/DNR/Fourier.cs。前者主要用于PanView中的频谱显示，后者用于解调出的音频信号的降噪处理。其实内容差不多，但是我更喜欢DNR目录里的FFT实现，因为这个目录下的FFT更接近比较原始的碟形算法，而且还实现了IFFT傅立叶逆变换，跟教科书里的相关公式比较接近。

这里推荐一个幻灯片，把FFT算法原理讲得很好。

https://wenku.baidu.com/view/5cacb2b8bd64783e09122b9a.html

要理解FFT代码原理，首先要说一下FFT算法的前身DFT。DFT是离散傅立叶变换，公式和连续傅立叶变换很接近，只不过是把连续傅立叶变换的积分换成了累加而已。理论上，用DFT的公式（参加幻灯片第四页），只要用软件代码对一个输入的数据序列进行公式里的这些乘法和累加后，就能计算出这个输入的傅立叶变换的结果。但是从公式也可以看出来，为了计算一个频点上的傅立叶变换的值，就要做N次乘法和长度为N-1的累加，然后要计算所有N个频点上的傅立叶变换的输出值，就要再重复前面这些运算N次。也就是N*N次乘法和N*(N-1)次累加。而且这些运算还是复数运算，也就是说前面说的那些乘法和加法是复数的乘法和加法，比实数的更复杂。具体对应的实数运算数量见幻灯片第五页，这里就不再推导了。

经过上面的分析可以知道，当输入序列的长度N很大时，运算量上升很快，这样如果要求实时运算的情形，原来的方法就不能用了。这时就要考虑利用公式中的W（旋转因子）的一些特殊性质来简化运算（见幻灯片第八页），省去很多重复运算。利用这些性质简化后的DFT就是FFT了。所以FFT与DFT的结果是相同的，只是运算速度提升了。

具体的优化方法也有很多种，比较主要的两种是幻灯片里介绍的按时间抽取-基2和按频率抽取-基2。建议先看一遍时间抽取的方法，比较好理解，幻灯片里会告诉你用了这种算法可以把长度为N的DFT等价替换为两个长度为N/2的DFT，然后再经过碟形运算输出同样的结果，根据前面的分析，序列长度减少后DFT的运算量是大幅减少的，即使考虑到后面的碟形运算所增加的运算量，总体的运算量还是少了（见幻灯片第16页），所以这种分解运算的方法是有效的。

由于长度N是2的L次方，（通常是1024这样的值，幻灯片里第19页上的例子是8），它除以2以后还是可以继续再除以2，也就是可以继续分解。每分解一次，运算量又可以进一步减少。分解完以后，对更小规模的DFT（幻灯片里是长度为2的DFT），先做一级小规模的碟形运算（这个碟形运算规模小，但是要做两次），恢复出前面说的两个N/2序列的DFT，再做一次前面说的碟形运算。见幻灯片19页。

那么对于之前说的N=8的序列，现在已经简化为4个长度2的DFT和2级碟形运算了，根据幻灯片20可知，长度2的DFT本身也可以直接用碟形运算实现。也就是说，把N不停除以2，完全分解后，所有的运算都可以用碟形运算实现，不必再做DFT也可以得到原来DFT产生的一模一样的结果。

现在的问题是，这样对原序列不停地做多次奇偶分解，输入的元素的顺序比较难找，幻灯片26页给出了方法。原理就不多讲了，但是可以大概看一下怎么用程序实现，首先，按照正常的顺序0、1、2、3、4、5、6、7，可以得到它们对应的二进制值000、001、010、011、100、101、110、111等，这时候只需要把它们各自的位序倒过来，即把每一个二进制数字的左边放到右边，右边放到左边，得到：000、100、010、110、001、101、011、111。然后再把它们转换成十进制，就得到了0、4、2、6、1、5、3、7，正好就是要找的输入序列的顺序，参考幻灯片28页。

那么知道了这个这个0、4、2、6、1、5、3、7的序列顺序后，怎样从原来的按照0、1、2、3、4、5、6、7的buffer里获得新的序列顺序呢？可以参照幻灯片29页的方式，对buffer里的各元素逐个交换，这个操作在代码里也有体现。

要实现运算，还要把旋转因子W用代码实现，首先可以根据前面N=8例子中的碟形算法每一级中用到的W归纳一下，总结出旋转因子的一般形式，这样就方便用软件代码实现并做自动化的运算了。可以看到其实这个旋转因子W就是自然对数e的复指数，W上的参数变化，对应于在复数平面的旋转，所以叫旋转因子。这个旋转跟矩阵运算里的矩阵旋转没关系。

现在来看sdrsharp里的FFT运算的代码，sdrsharp/DNR/Fourier.cs中，第17～21行在求长度N对应的2的指数幂。比如N=1024时，这里算出来的结果是10，也就是说1024=2^10，如果N=8时，算出来的结果就是3。这里代码里的length变量就是N，代码里的m变量就是计算结果。原理也很好理解，i=i>>1的作用就是i=i/2，不停地把length除以2，除一次就给m加1，这样就能计算出这个2的指数幂了。

接下来，代码第25行～第46行做的是幻灯片第29页的操作。代码里的samples就是之前说的buffer，是同一个意思，samples里最早存储的是输入的时域信号序列，为了更加一般化，这个输入序列也看作复数，然后在同一个数组里交换不同位置的元素，来达到幻灯片29页的效果。第29～34行做的就是交换，只不过把复数元素的实部和虚部分开做了（应该是因为c#里不支持直接的复数运算，matlab里可以），另外tr和ti只是用于交换的临时变量。这部分代码的交换次序和幻灯片相同，只是在找哪个下标需要交换的时候没有用幻灯片28页里的倒位序的方法，但是实际效果是一样的。可以假设N=8的情况，把nd2=4代入到这段代码中，可以看到效果与幻灯片29页相同。这个循环对整个数组去头去尾，表示头和尾都不做任何操作，即数组元素0和数组元素7都保持原位，中间有几位做了交换，具体来看，数组元素0没有做任何操作，保持原位，数组元素1和数组元素4交换，数组元素2保持原位，然后数组元素3和数组元素6交换，数组元素4也没有做任何操作（因为之前已经换过一次了，换好了），然后数组元素5也保持原位，数组元素6也没有做任何操作（与数组元素4一样，换过一次了不需要再重复换），数组元素7也保持原位。长度为8时，实际有效的交换只做了2次。第37行～第45行用于找出要和当前位置i处的元素交换的位于j处的元素。之前说的倒位序方法也能达到相同的效果，即找出目标的排序，然后比较当前的排序，对当前排序的元素进行遍历，如果和目标排序相同就不交换，如果和目标排序不同就交换，如果和目标排序不同，但是之前已经做过次交换也不交换。但是为什么37行～45行也能实现这个功能还没完全理解。

幻灯片29页的变址处理（对输入数组做交换）完成后，接下来就要做碟形算法了。对应于代码第48行～第78行。要看懂这段代码，首先有几个地方要搞懂。一个是复数乘法的实现方式，可以参照幻灯片第5页，复数乘法输出的实部等于输入的两个复数的实部相乘的结果减去输入两个复数的虚部相乘的结果，而复数乘法输出的虚部等于两个输入复数的实部乘以另一个的虚部的和，这个如果不理解可以自己用两个复数推导一下复数乘法的公式。其实代码里的第65～66行以及第74～75行都在做复数乘法，前者在做的是旋转因子与输入序列间的复数乘法，后者做的是旋转因子本身的在旋转过程中要经历的复数乘法。

有几个变量可以对应看一下，第48行的l变量就对应于第33页的L，也就是第L级。（31页从左往右，最左是第一级），然后第48行里的m变量对应于第33页的M，如果N=8的话，M=3也就是总共有3级，所以这里m也等于3。le2对应于幻灯片第32页的距离，这个距离的意思是每个碟形运算的间距，碟形中间一行都没隔开，这个距离就是1，隔一行距离就是2，隔三行距离就是4。le2在这里循环里会产生的变量正好就是1、2、4。而le在这里其实也表示距离，它代表的是在一级中使用同样的旋转因子的两行之间的距离。因为从第31页碟形图中看到，尤其是中间的第二级和第三级，旋转因子本身是会变化的，比如第二级，我们可以先把使用同一个旋转因子W0N的所有的运算都做完了（从第一组运算调入第二组运算时要用这个le来找到合适的位置），再对旋转因子本身做旋转，做另一个旋转因子W2N对应的运算。

接下来是ur和ui，它们代表旋转因子在最早期W0N的虚部和实部，从幻灯片第33页最下面的旋转因子一般形式知道，其实旋转因子就是在单位圆上用指数形式表示的复数值，当r=0时，指数幂就是0，那么这个旋转因子正好等于1，即实部是1虚部是0。所以最开头ur=1，ui=0。注意第33页最下方的公式里有2个j，其中最左侧的表示复数，另一个表示的是同一级中的j取值，就是33页上方右侧的j，第一级j只能是0，第二级j可以是0、1，第三级j可以是0、1、2、3。这个值在代码里对应于第60行的jm1。

往下看，sr和si是用来构造后期的旋转因子的。它们是cos和-sin，正好是e^(-j)的实部和虚部。然而，旋转因子在每一级中的转动幅度是不一样的，越高的级旋转幅度越小，增加一级旋转幅度缩小一倍，可以在第33页最后一行中间的公式里看到，L增加1，指数上的分母就要多乘以一个2，整个指数相当于除以2，这样每次j增加后的旋转效果也会相应缩小。这个在第55行和第66行的分母上的le2也有所体现。再往下的变量jm1前面已经说过，再往下看ip变量和i的关系就是碟形运算中两个输入值的位置的关系。只要知道碟形运算的上方的输入值的位置再加上它们间的距离le2就可以找到与它做碟形运算的另一个值的位置。然后tr和ti是碟形运算的乘法的结果，可以看到都是碟形运算下方的输入乘上了当前的旋转因子（从ur和ui中体现），乘法项做完了就做加法，根据第14页，碟形运算上方的输出是前面两个值相加，下方的输出是前面两个值相减，67行～70行就在做这个事，samples[ip]就代表下方的输出，samples[i]代表上方的输出。再往下看就是当同一级中旋转因子有好几个时，对其本身进行旋转生成下一个旋转因子的操作，u代表旋转因子的当前值，s代表旋转因子本身要转动的值，tr只是暂时保存ur的中间变量。

看完了变量可以看一下三个循环各自代表的功能，来对整个代码结果做一个整体上的理解。首先，最外侧的循环，对应于第48行的for，代表了当前运算对应第几级。第58行，对应于同一级中有几个旋转因子，所以这个循环在结束前的第73行～75行做的运算就是在计算这一级的下一个旋转因子。第62行的for循环，对应于同一级且同一个旋转因子所要做的运算。比如最左侧的第一级，一共就只有一个旋转因子，所以第58行的for循环只做了一次，而这一级中同一个旋转因子W0N用了4次，做了4个小的碟形运算，所以第62行的for循环做了4次。

另外，还可以看一下这3个循环每次变量增加的值，以及变量最大值。48行的for循环对应级的上升，显然，每个级都要做，所以每次增加1，然后级的最大值就是M所以不能超过代码里的变量m。58行的for循环，对应同一级旋转因子，每个旋转因子参与的运算都要做，所以每次增加1，由31页图可知每一级旋转因子的数量正好与这一级碟形运算的距离（碟形大小）相等，比如第二级中有2个旋转因子，而这一级的碟形大小间距就是2，所以这个循环的最大值（即旋转因子数量）就是le2。62行的循环对应于同一级中同一个旋转因子所做的所有运算，同一级同一旋转因子的距离就是用le表示的，所以这个循环的变量增加值就是le，为了遍历序列的所有输入，所以这个循环的最大值是nm1，也就是总长度-1。

看完了以上这些可以说FFT正变换就看懂了。接下来的IFFT其实是利用了傅立叶变换本身的性质，也是调用FFT实现的，并不复杂。如果要把一个复数序列做逆变换，只需要把这个输入序列的虚数部分符号反一反，然后对它做正变换，做完后再乘以一个缩放因子，然后再把它的虚部符号反一反，所得的结果正好就是直接对这个复数序列做逆变换的结果。这个可以从幻灯片第43页看出来，也可以自己进行公示推导。幻灯片第45页也再次做了证明。幻灯片里的取共轭就是实部不变，虚部符号反一反。

对应于代码可以看到第82～85行，先给输入的复数序列的各个元素的虚部加了一个负号，相当于对这个复数序列取共轭，然后87行做了FFT正变换，89行算了缩放因子。91行～95行的循环把前面的FFT正变换结果用缩放因子做了缩放，并对虚部加了负号相当于又做取了一次共轭。和幻灯片第45页的操作基本一致（只有求共轭和除以N的次序反了，但是这个没区别的）。

这样FFT的正变换和逆变换在sdrsharp里的实现就讲完了，看懂这些后再看sdrsharp里其它信号处理代码应该也没什么问题了。