稀疏傅里叶变换原理说明（一）

前言

因为某些情况，钻研了几个月的稀疏傅里叶变换，总算从“一无所知”到“对其有了一些了解”，因此写下这篇博文以便参考。一是希望给刚开始研究稀疏傅里叶变换的网友一些经验，少走弯路；二是由于博主并非通信专业出身，对数字信号处理几乎一窍不通，导致即使我知道了稀疏傅里叶变换的大致原理，仍然无法彻底搞懂这个算法，程序中的各种细节也看不明白，因此希望和广大网友相互交流借鉴。

由于水平有限，博文之中可能存在解释不当的地方甚至是错误，欢迎各位指正！

参考代码链接：bernielampe1/sparse_fft: Implementation of sparse fast Fourier transform (sFFT). (github.com)

一、稀疏傅里叶变换简介

稀疏傅里叶变换，英文名为Sparse Fourier Transform，一般缩写为SFT或者SFFT（Sparse Fast Fourier Transform），下文中都称之为SFT。SFT是一种利用信号自身特点得到信号频谱的算法，是对信号处理中最经典的算法之一——离散傅里叶变换（缩写为DFT，下同）的一种改进。所谓“信号自身特点”，对一维信号而言，一般指的是其频率范围内只有少部分频率分量幅值较大，其他大部分频率分量的幅值都趋近或者等于0。举个例子：一个信号在0~1000Hz的范围内，只有50Hz和120Hz这两个频率上具有能量，其他频率上几乎没有能量。这一特点也被称为“稀疏性”。

注：一个只有在50Hz和120Hz上有“有用信号”的稀疏信号

稀疏性一般使用稀疏度来表示。稀疏度（用K表示）即有用频率成分的个数。以上图为例，稀疏度K=2。

SFT最先是由MIT在2012年的论文《Simple and Practical Algorithm for Sparse Fourier Transform》中提出来的，该博客也主要是基于该论文进行的研究。

二、稀疏傅里叶变换流程

SFT的理论框架如图。

这里对各个部分进行简要说明：

• 频谱重排：利用DFT中频域和时域的性质，对时域序列进行操作，使得信号的频谱重新排列
• 滤波：使用提前设计好的滤波器，将之与重排后的信号相乘，起到扩展谱线、截短序列的作用
• 降采样FFT：仍然是利用DFT的性质，通过对时域序列求和，实现对频谱的采样
• 定位：利用哈希逆变换，从降采样后得到的频谱反推出原信号谱线所在的位置
• 求值：利用滤波器的频谱和估计得到的位置，反推出原信号谱线的幅值

下面会对各个部分进行详细解释。

三、流程中各部分详解

3.1 频谱重排

频谱重排主要是利用了DFT的性质：

1. 位移性质：若时域上有：$$ x_{1}=x(n)W_{N}^{bn}$$ 则频域上有：$$ X_{1}(k)=X(k+b) $$
2. 缩放性质：若$$ x_{2}=x(\sigma n)$$ 则频域上有：$$ X_{2}(k)=X(\sigma^{-1} k) $$

其中σ^-1并非和σ互为倒数，而是一个整数，该整数满足关系：$$\sigma \sigma^{-1} mod(n) = 1$$

也就是σ^-1与σ的乘积对n取模，所得余数等于1。此时称σ^-1是σ关于模n下的逆元（博客中n和N含义相同，n论文中常用，N个人更喜欢用）。举两个例子：

7在模24下的逆元就是本身，因为$$7\times 7\div 24=2\cdots 1$$

15在模128下的逆元是111，因为$$15\times 111\div 128=13\cdots 1$$

当然并非所有数都存在逆元，不过这不在讨论范围之内。

缩放性质是SFT中重要的理论基础，其可以使用matlab进行验证。这里需要注意的是，超出了序列范围之外的下标，需要对序列长度求模进行约束。

N = 128;
F_signal = zeros(N, 1);
F_signal(60) = 5;
signal = ifft(F_signal); % 利用反变换得到信号序列

p_signal = zeros(N, 1);

% 进行频谱重排，其中重排因子为 15，其关于模 128 的倒数为 111
for i = 0:N-1
    p_signal(i+1) = signal(mod(i*15, N)+1);
end

F_p_signal = fft(p_signal); % 对新序列进行傅里叶变换，查看其频谱

subplot(2,1,1);
stem(F_signal);
subplot(2,1,2);
stem(abs(F_p_signal));

所得结果如图：

由于matlab下标不能从0开始，因此上半部分的谱线实质对应着序号59。以15作为重排因子，128作为模，所得新信号的谱线对应位置为117，满足：

$$\left\{\begin{matrix}
59\times 15\div 128=6\cdots 117
\\ 
117\times 111\div 128=101\cdots 59
\end{matrix}\right.$$

3.2 滤波

虽然叫滤波，但这个的滤波过程并非是传统意义上的滤出信号中的有用频率成分。首先需要知道的是，在时域上两个信号进行卷积操作，等于在频域上对这两个信号进行相乘操作。反过来，在时域上对两个信号的对应元素进行相乘，则对应于在频域上对这两个信号进行卷积。而这里的滤波，指的就是后一种方式。下图就显示了这两种方式的区别。

N = 256;
omega = 50;

% 生成信号和幅值
F_signal = zeros(N, 1);
F_signal(117) = 4;
F_signal(210) = 1;
signal = ifft(F_signal);

% 生成矩形窗滤波器
F_box_car = zeros(N, 1);
F_box_car((N-omega)/2 : (N+omega)/2) = 1;
box_car = ifft(F_box_car);
F_box_car = fft(box_car);

% 相乘
times_b_s = signal .* box_car;
F_times_b_s = fft(times_b_s);

% 卷积
conv_b_s = conv(signal, box_car, 'same');
F_conv_b_s = fft(conv_b_s);

% 绘图
subplot(4,1,1);
plot(abs(F_signal)); title('原信号的频谱');
subplot(4,1,2);
plot(abs(F_box_car)); title('滤波器频谱');
subplot(4,1,3);
plot(abs(F_times_b_s)); title('时域上相乘后信号的频谱');
subplot(4,1,4);
plot(abs(F_conv_b_s)); title('时域上卷积后信号的频谱');

 

第四个子图中展示的行为就是我们常见的滤波器的操作。在第三个子图中，原信号的频率分量被滤波器展开并移动，而从幅值的大小可以大致猜出发生了怎样的变化。

使用时域相乘这种方法的要点有两个：

1. 同一频率分量经滤波器扩展后，幅值相同（假如滤波器理想），因此在频域上采样得到的幅值可以代表原信号频率成分的幅值；
2. 降低了运算量，比如原信号序列长为4096，但滤波器序列可以仅为256，只需取出信号序列中的256个点与滤波器进行相乘即可。

值得一提的是，滤波器的性能直接决定了SFT的可靠性，而滤波器问题也是困扰我的问题之一。理想情况下，直接使用矩形窗滤波器就能完成所有事情，但理想的矩形窗滤波器对应着一个时域上无限长的sinc序列，这显然无法实现。现实中，我们常用过渡带平滑的滤波器与理想矩形窗相乘或卷积来改善矩形窗性能，使之在一定序列长度下，幅频特性尽量接近理想矩形窗的幅频特性。处理后的矩形窗被称为“平坦窗”。

如何实现平坦窗的自动化生成，一直是我在考虑的问题。论文中关于这一段的描述如下：

Let B be a parameter that divides n, to be determined later. Let G be a (1/B, 1/(2B), δ, ω) flat window function for some δ and ω= O(Blog(n/δ)). We will have δ≈1/n^c, so one can think of it as negligibly small.

B是分筐数，后面会提。B一般满足关系

$$B\approx \sqrt{nk}$$

以序列长度n=8192、稀疏度K=7为例，B=239，由于B应该能整除n，因此取与之最接近的256。按照上面的公式，取δ≈1/n^c≈1/n，则平坦窗宽度ω=256*log(8192²)≈4614，显然与时域相乘中的要点（2）相冲突。即使将ω进行缩小，缩小到何种程度也是一个问题。另一方面，平坦窗由两个窗卷积得到，一般选择矩形窗和高斯窗。如何利用已知参数（分筐数B、稀疏度K、序列总点数N等）计算出合适的参数来生成矩形窗和高斯窗，两者的卷积能满足(1/B, 1/(2B), δ, ω)这一条件，也是一个非常麻烦的问题。

---

 

由于篇幅限制，后续内容放在了“稀疏傅里叶变换原理说明（二）”和“稀疏傅里叶变换原理说明（三）”当中。