稀疏傅里叶变换（sparse FFT）
作者：桂。
时间：2018-01-06 14:00:25
链接：http://www.cnblogs.com/xingshansi/p/8214122.html
前言
对于数字接收来讲，射频域随着带宽的增加，AD、微波、FPGA资源的需求越来越高，但频域开的越宽并不意味着频谱越宽，有限信号内可认为信号在宽开频域稀疏分布，最近较为流行的稀疏FFT（SFFT）是在传统FFT的基础上，利用了信号的稀疏特性，使得计算性能优于FFT。本文简单记录自己的理解。
一、稀疏FFT
主要是12年MIT的论文：Simple and Practical Algorithm for Sparse Fourier Transform。
核心思想：
SFT 作为这样一种“输出感知”算法,其核心思路是按照一定规则 Γ ( • )将信号频点投入到一组“筐”中(数量为 B，通过滤波器实现 ) . 因频域是稀疏的,各大值点将依很高的概率在各自的“筐”中孤立存在 . 将各“筐”中频点叠加,使 N 点长序列转换为 B点的短序列并作 FFT 运算,根据计算结果,忽略所有不含大值点的“筐”,最后根据对应分“筐”规则,设计重构算法 Γ -1 ( • )恢复出 N 点原始信号频谱。
算法流程：
步骤一：选定一个sigma，实现数据频域重排。
频域重排主要是因为FFT之后，频域大值集中在一起，需要将这一伙打散，保证频域的稀疏性（sparse）。打散之后就能稀疏，依据定理（可见n越大，打散的可能性越大，从这里看，sigma与B还是有关系的，sigma体现了相邻频点的最小间隔，而B决定了每个bucket的宽度）：
这里是按一定概率（概率性）将频谱重排，MIT重排思路：
对应代码：
fs = 1000; f0 = 70; N = 256; t = [0:N-1]/fs; x = exp(1j*2*pi*t*f0) ; %%*************Step1:频谱随机重排************** sigma = 19; %inv(sigma) = 27 tao = 3; %permutation perm_num = mod([1:N]*sigma+tao,N)+1; pn = x(perm_num);
当然也有强制（确定性）将数据重排的思路，思路类似，只是实现方式不同：
sigma的选择主要利用性质：
这里sigma^-1是sigma关于模N的逆元(数论倒数)。该点主要说明：变化前后的完备性（信号可重建）。
步骤二：加窗。
这里根据筐的多少（B），平分2pi频域，即带宽2pi/B，理想窗函数为矩形窗，但时域为无限长的sinc函数，需要加窗截断，可选择gausswin。
即win = sinc.*gausswin：
这里不局限于gausswin，满足给定约束的窗函数均可：
步骤三：频域抽取并作FFT
加窗之后，保证了频谱不至于展宽严重，进一步保证了稀疏性。频域降采样：
等价于时域的混叠的傅里叶变换：
故直接在时域进行处理。处理完成后进行FFT运算。
该操作的理论基础为：
不谈加窗操作，x->p->y，可以看出x与y存在映射关系，而y的∑操作可以转化为并行，这样一来可以用多个低速率AD实现一个高速率AD的功能，前提是多个AD需要完全同步。
与上述降速率思想等价的是：如果各AD可做好同步，则多个现有AD的能力，可以做出现有AD难以完成的事。
步骤四：哈希映射
哈希映射的线索为：最终观察的数据Y(k)【即y(n)】->P(k)【即p(n)】->X(k)【即x(n)】：
第一步（y -> p）转化的理论依据：
第二步（p -> x）转化的理论依据：
序列重排后的频域变换为：
这两步解释了算法step4的参数定义：
但实际中并非第一步提到的理想情况，实际情况是：
因此存在一个频谱重建的成功概率问题：
4章节的后半部分主要在证明这个概率问题。
步骤五~六的主要目的，引用原文的说法：
Here are two versions of the inner loop: location loops and estimation loops. Location loops are given a parameter d, and output a set I ⊂ [n] of dkn/B coordinates that contains each large coefficient with“good” probability. Estimation loops are given a set I ⊂ [n] and estimate x I such that each coordinate is estimated well with “good” probability.
步骤五：定位循环
d是一个参数，存在约束
原文取值为：
该步骤的主要操作为：
将 Z ( k )中 dK 个较大值（从大到小依次排列）的坐标（ k）归入集合 J 中,通过哈希反映射得到 J 的原像：
这样便得到包含原始频谱坐标的集合，迭代L次：
迭代的目的？
步骤六：估值循环
对于每个k∈Ir，计算频率信息：
步骤5~6主要操作：
至此完成了稀疏FFT（sparse FFT）的整个运算过程，记作sfft_v01。
改进版都是依次大框架进行，不同点主要有三个方面：1）重排的实现思路；2）频谱的重建思路。不再一一展开。
原文示例：
该示例给出了步骤5~6的直观解释：
二、问题记录
仿真验证：n = 256点频信号，sigma = 19，则inv_sigma = 27，假设生成点频信号：
图1为重排的信号，图2为重排加窗(此处窗不够理想，看到长尾，论文中强调加窗，作用在于把尾巴剁掉。)，图3为原始信号频谱，图4为降采样之后的频谱。可以看出重排容易引入谐波：以sigma为周期。
原始频谱位置为k=10-1 = 9（下标从0开始），重排后对应频谱位置为mod(sigma*k,n) = 171 = 172-1，与理论分析一致。即原始k点重排后位置：
即时域信号，x(n)，重排序号perm_num = mod([1:n]*sigma,n)+1，则对应重排后的傅里叶变换，频谱顺序与原信号频谱X，经过重排序号perm_num1 = mod([1:n]*inv_sigma,n)+1，位置始终相差1（数论倒数的性质决定）:mod(sigma*inv_sigma,n)=1
