書接上回：稀疏傅里葉變換原理說明（二）
三、流程中各部分詳解
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)。當然這代碼也是根據稀疏傅里葉變換原理說明（一）開頭提到的參考代碼改的，希望能有所幫助。感謝看到最后，歡迎討論！