哈哈哈哈好久没有出现了，这次把上次比赛所用的3780fft算法写成了论文。
摘要
本文是关于3780点快速傅里叶变换算法的讨论，由于常用的快速傅里叶算法对采样点数要求为2的指数倍，而本次讨论的采样点为3780个点，并不能直接使用传统的基-2FFT算法，因此需要使用新的算法来进行FFT运算。本文采用了混合基FFT算法、素因子算法(PFA)及Winograd FFT(WFTA)算法混合进行计算，通过matlab进行仿真，速度和精度都能得到较为理想的结果。
关键词：FFT、PFA、WFTA、下标映射
简介
由于中国数字电视地面广播传输系统所采用的是时域同步正交频分复用（TDS-OFDM)分别对信号进行调制和解调，其子载波数为3780点，故对3780点FFT进行研究具有一定的现实意义。
面对3780点的FFT，通常可以比较容易想到的方法是利用内插的方法，把3780点变成4096点，再使用基-2FFT算法，得到4096个数据后减采样到3780点。由于这个过程中由于并不是整数倍的增减，必然会带来误差，而且这个操作会增大计算的复杂度，降低速度。而本文所介绍的方法，是参考专利[1]的方法，使用几种快速算法进行混合运算，底层使用WFTA算法进行7、9、3、4、5点的运算，中层采用PFA进行63、60点的运算，最顶层采用混合基FFT进行运算，如图1所示。
图 1 3780点FFT计算流程
算法解析
-
库利-图基基-2FFT
目前常用的基本FFT算法是基-2库利图基算法。要求计算的点数为2的幂次，首先定义：
（1）
则离散傅里叶变换如下：
（2）
可以通过复指数的周期性、对称性等特性在时域上抽取将此式化为：
（3）
可以观察到，内部仍为DFT，但长度缩短了一半。这样分解为了奇数和偶数两组的DFT可以继续递归，直到全部分解完。
这个算法充分利用DFT中重复使用的值，极大地减少了计算的复杂度，且由于计算的特征，使得每级计算结果可以存回原地址，虽然最终计算完地址会变成二进制编码的倒置，但只要重新排序就可以输出正确值。
-
库利-图基混合基FFT
混合基FFT或者叫多基FFT是在基-2FFT基础上的推广，由于因数并不像基-2计算那么规则，所以地址映射规则有所变化。这里引出地址映射的一般法则。
若N=N1*N2, 则分别取n1=0,1…N1-1, n2=0,1…N2-1,可以得到线性映射：
(4)
其中q1、q2均为整数此时分为两种情况：
- 第一类映射：若N1、N2互质即( N1,N2)=1,要满足一一映射，利用广义孙子定理，通过一些条件化简后可得，详见参考文献[2][3]：
(5)
k1、k2取值范围分别与n1、n2相同。
- 第二类映射：若N1、N2不互质即( N1,N2)>1，一一映射的充要条件是:
(6)
以上出现的a、b均为整数，频域下标k映射也是同理的。
所以对于任意合数的FFT运算，可以采样这种方法，按照第二类下标映射的规则设：
（7）
将（7）代入（2）中化简可得：
可以看出实际上是N1个长度为N2的DFT，再乘以旋转因子，再做N2个长度为N1的DFT。
其实基-2FFT算法中的奇偶分组的方法本质上也是下标映射，将下式代入（8）化简后即可得到（3）:
（9）
-
素因子算法
对于可以分解为两个互质因数的合数，这时可以按照第一类下标映射做运算，代入（2）化简可得：
（10）
此式和（8）十分相似，但属于不同的算法。相比于（8）少了旋转因子，乘法计算量大减，但引入了更为复杂的下标映射，且分解的因子必须是互质的。
此式可以推广到N维的分解，首先下标映射变成N维，即N= N1*N2*……*Nl，详细推导见参考文献[4]：
（11）
代入到（1）可得到如下的N维素因子算法：
（12）
-
Winograd FFT
WFTA同样可以嵌套，把做行和列DFT所需的乘法合成一个放到算法中间，形成一个更为紧凑的算法结构，使得乘法的计算次数比起PFA进一步减少，但加法次数略微增多。
前面几个算法都可以做到同址运算，而WFTA由于在计算不同点数时所使用的计算核不同并且每个计算核的算法独立，所以需要定制比较复杂的算法结构，而且没办法做到同址，在做大点数FFT时会消耗比较多的存储单元，故而通常只在计算点数比较少的DFT时才使用WFTA。
本文所用到的小N点的计算公式在参考文献[3]已经推导，此处不再重复。需要注意的是，[3]中关于9点WFTA公式的推导中S1= m0+m2+m3有误，应修改为S1= m0+2*m2。
3780点FFT的算法描述
-
第一级
通过公式（7）（8），将3780分解为N1=63，N2=60，下标映射为：
（13）
按照n的下标映射取60次63个点做同址运算，全部算完后按n的下标映射每个点乘上对应的旋转因子，之后再按n的下标映射取63次60个点做同址运算，最后最终结果按k，计算结果按n的映射读出点，k1 和n1取相同值，k2 和n2取相同值。
-
第二级
通过公式（11）（12）将60分解为N1=3，N2=4，，N3=5下标映射为：
（14）
除了无需乘旋转因子外，其余步骤与第一级相同。
-
第三级
第三级为最底层的DFT，全部直接采用WFTA所常用的小N点DFT算法。
MATLAB实现
本次测试利用matlab生成了一组时域数据，其为按下式取得的3780个点：
其中n=0,1,2……3779。
（15）
将公式（1）以及前面一节的描述转换成matlab的程序，并比较和matlab的内建函数fft()计算数据（15）的结果，得到的结果为图2。
图 2 3种方法计算数据对比
从图2可以看出，本文所采用的算法计算的结果与DFT公式计算的结果以及matlab内建函数fft计算结果高度吻合。通过使用matlab的corr2()函数计算3组数据的相关度，计算结果均为1，故而三组计算结果其实相同，可以得出本算法具有非常高精度的结论。
通过测试三种算法各100次，统计出了3组在测试电脑上各算法计算时间的数据，通过matlab绘制在图3中。
图 3 三种算法计算时间对比
从图3可以非常直观的看出本文算法对于减少DFT时间的作用，对3组数据求平均值后得到普通dft算法所使用的时间为3.9133s，本文的算法为0.046628s，内建函数fft()的时间为0.00006s。通过数据可以得到在测试环境下，本文算法相比直接使用公式（1）计算速度快了约84倍，但fft()函数仍比本算法快约745倍。但由于matlab对内建函数做了并行以及许多系统级的优化，而本函数使用了大量for循环，未针对matlab的矩阵计算特性进行优化，有如此差距是可以理解的。所以理论上在真实应用的情况下，本算法速度能够更快。
-
代码
本文的相关代码已上传到github，点击链接可以查看：代码地址
结论
本文所论述的是一种与常用算法不同的较为新颖的算法，利用了多种FFT算法混合的方法做到了能够快速准确地计算3780点DFT。而且由于前两层都可以进行同址运算，故而可以节约一定的硬件成本，在实际使用中会具有很高的实用性。
参考文献
[1]福州大学. 一种实现3780点FFT/IFFT的方法及其处理器[P]. 中国专利. 2011101385900. 2011-05-26.
[2]胡广书. 数字信号处理—理论算法与实现[M]. 北京：清华大学出版社，2003.
[3] 努斯鲍默,H.J.著, 胡光锐译. 快速傅里叶变换和卷积算法[M]. 上海: 上海科学技术文献出版社, 1982.
[4]冷建华. 傅里叶变换[M]. 北京：清华大学出版社，2004.
2018.11.23.更新
由于在GitHub上看到有人提问IFFT怎么做（快4个月过去了才看见23333），yellowko表示没有专门研究这个算法直接IFFT怎么做，但是可以直接使用FFT来做IFFT，只会增加一些共轭的运算，具体做法是这样：
x[n]=(FFT((X[k])*))*/N
推导过程很基础的啦，如下：
首先从DFT的公式开始
对IDFT做两次共轭，可以发现里面就是DFT
本作品由yellowko采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可。
咕，Github上的代码好像有错误
在ffttest.m中
第46行应该是
disp([‘fft:’,num2str(mean(tfft))])
第49行应该是
disp([‘nmix3780算法和dft算法计算结果相似度，实部：’,num2str(corr2(real(b),real(c))),…
‘虚部：’,num2str(corr2(imag(b),imag(c)))])
是的，感谢纠正