书接上回：稀疏傅里叶变换原理说明（二）
三、流程中各部分详解
3.5 估值
在频谱中，需要知道两个信息：频率和幅值。上一步“定位”解决的是前者，也就是确定哪些频率分量是真实使用的，而这一步就解决后者，即每个频率分量的幅值是多少。不过幅值估计的想法非常简单：令筐中的幅值除以滤波器的增益得到估计值est，然后假设每个筐经过逆映射得到的频率值对应的幅值都是est。将L次循环的结果对应相加，并分别对幅值和相位取中位数，就得到了幅值估计的最终结果。
仍以稀疏傅里叶变换原理说明（二）中的例子为例。在单次定位循环中，已经知道了第25号“筐”对应着的原始频率417、1700、2983、170、1453、。。。2139、3422可能存在有用频率。（题外话，这里不是指的真实频率就是417Hz，傅里叶变换中，坐标（频点）和真实频率之间的关系是要经过换算的，所以这里用频点更合适，但之前的博文都这么写过来了，懒得改了）
假设使用的滤波器通带增益近似为1，25号“筐”的幅值是0.0218958，除以滤波器增益后的结果为0.0218152，因此原始频率417、1700、2983、170等的幅值都假设为0.0218152。同理26号“筐”的幅值约为0.01，同样的处理方式，认为26号“筐”包含的原始频率的幅值都为0.01。最终结果如图所示。
在“定位”中也提到，需要进行L次定位。同样的，也需要进行L次“估值”。将L次的结果对应相加，取出幅值最大的K个频率及其对应的幅值，分别对其幅度和相位求中位数，就得到了估计后的结果x_hat。这里给出对应的matlab代码，应该能够更好的理解。
% Z = zeros(L, N)，用来存储恢复后得到的幅值
% x_hat = zeros(N, 1)，用来存储最终结果
temp_x_hat = sum(Z ~= 0);
[~, x_hat_index] = sort(temp_x_hat, 'descend');
x_hat_index = x_hat_index(1:K); % 从估计序列中选取 K 个最大值所在的坐标，作为估计得到的坐标值
for i = 1:K % 遍历估计坐标值
x_hat_idx_i = x_hat_index(i);
re = real(Z(:, x_hat_idx_i));
im = imag(Z(:, x_hat_idx_i));
re = median(re); % 取实部的中位数
im = median(im); % 取虚部的中位数
complex_i = re + 1j * im;
x_hat(x_hat_idx_i) = complex_i;
end四、结果
这里给出一个以N=4096、K=6、无噪声的信号作为原信号，分别进行两种变换后的例子。
从结果可以看出，除了幅值不同外，频率的估计和整体趋势的估计还是非常正确的。幅值不一致，是因为我们丢掉了大部分原始序列中的点，降采样操作也会造成能量损失。
这里没给出时间上的对比结果。按照理论：当N增大时，SFT的速度会逐渐超过FFT。但在这个例子中，FFT的运算时间比SFT时间少10倍，不仅如此，当N增大，如果运算精度d和循环次数L不做出一定的改变，SFT将无法正确估计出原始信号；而增大d和L又会导致运算时间大量增加，真的是麻中麻。
推测这是因为大量的排序操作导致运算时间提高，但这一猜想没有得到太多的验证。
五、一些看法和碎碎念
关于SFT，我想提一些我自己的看法。
一维SFT主要解决的是对序列点数非常多的稀疏信号进行快速频谱分析的问题，因为在N非常大（一般而言要上万）时，SFT的理论运算速度会显著快于FFT。问题在于，这个“快”，在点数不多的时候不够明显，而且SFT不像FFT那么成熟，有非常好用的库函数，编程和调试难度也很大，程序写起来非常痛苦，至少我现在都没在matlab上编好比较满意的SFT代码。
注：论文中给出的SFT和现有的一些FFT代码库之间的运算时间对比
二维SFT听说比较好用，但我没有对其进行研究，故不予评价。
编写代码时，滤波器方面的问题不好解决，各个参数之间还相互约束，因此有一段时间我一直在考虑SFT代码性能不佳的原因是不是相关参数没弄好，调的身心俱疲也没看出什么效果。所以老板的研究方向和自己的专业不一致，真的是说不完的辛酸，读研太苦。
这里给出源码链接：sparse_fourier_transform-Matlab: 稀疏傅里叶变换（sparse fourier transform，SFT）的Matlab实现。 (gitee.com)。当然这代码也是根据稀疏傅里叶变换原理说明（一）开头提到的参考代码改的，希望能有所帮助。感谢看到最后，欢迎讨论！