Rife算法确实不行

尝试实现了一下rife算法，参考这篇博客
Rife算法的一种Matlab实现
原博客只输出了一个相对误差分布图

看不出有没有改进。我在此基础上输出了使用rife算法前后的估计频率，发现效果确实很差

代码更改：

rife算法的函数

function [fc,fraw] = rife(spec,fs,N)
%% Rife-Jane to estimate peak frequency
% spec: input amplitude spectrum
% fs  : frequency resolution
% N   : length of signal
% fc  : estimated peak frequency
% fraw: raw estimated peak frequency without rife
[valueMax,posMax] = max(spec);
r= 2*((spec(posMax+1) > spec(posMax-1))-0.5);
spec_r=spec(posMax+r);
delta = spec_r/(valueMax+spec_r);
fc = (posMax+r*delta)*fs/N;
fraw=posMax*fs/N;
end

效果测试：

测试代码

clc;clear;close all;
%% 参数设置
N = 1000; % 数据长度
t = 1:N; % 采样时间
fs = 16; % 采样频率
f0 = 4.9; % 信号频率
nSNR = 200; % 模拟信噪比次数
SNR = (1:nSNR)*0.1;SNR=SNR-SNR(ceil(end/10)); % 模拟信噪比
%% 信号生成
signal_0 = cos(2*pi*f0/fs*t)';
signal = zeros(N,nSNR);
for k = 1:nSNR
    signal(:,k) = awgn(signal_0,SNR(k)); %加入高斯白噪声
end
%% 窗函数
rect = rectwin(N);
hn = hann(N);
hm = hamming(N);
bk = blackman(N);
%% 加窗
signalRectWin = repmat(rect,1,nSNR).*signal;
signalHann = repmat(hn,1,nSNR).*signal;
signalHamming = repmat(hm,1,nSNR).*signal;
signalBlackman = repmat(bk,1,nSNR).*signal;
%% 直接FFT获取振幅信息
wRectWin = abs(fft(signalRectWin));
wHann = abs(fft(signalHann));
wHamming = abs(fft(signalHamming));
wBlackman = abs(fft(signalBlackman));
%% Rife算法
% rectangle window
fcRectWin = zeros(1,nSNR);
fcRectWinRaw = zeros(1,nSNR);
for k = 1:nSNR
    [fcRectWin(k),fcRectWinRaw(k)]= rife(wRectWin(1:N/2,k),fs,N);
end
rmseRectWin = sqrt(abs(sum(fcRectWin-f0)));
% hanning
fcHann = zeros(1,nSNR);
fcHannRaw = zeros(1,nSNR);
for k = 1:nSNR
    [fcHann(k),fcHannRaw(k)] = rife(wHann(1:N/2,k),fs,N);
end
rmseHann = sqrt(abs(sum(fcHann-f0)));
% hamming
fcHamming = zeros(1,nSNR);
fcHammingRaw = zeros(1,nSNR);
for k = 1:nSNR
    [fcHamming(k),fcHammingRaw(k)] = rife(wHamming(1:N/2,k),fs,N);
end
rmseHamming = sqrt(abs(sum(fcHamming-f0)));
% Black man
fcBlackman = zeros(1,nSNR);
fcBlackmanRaw = zeros(1,nSNR);
for k = 1:nSNR
    [fcBlackman(k),fcBlackmanRaw(k)] = rife(wBlackman(1:N/2,k),fs,N);
end
rmseBlackman = sqrt(abs(sum(fcBlackman-f0)));
%% 误差分析
figure;hold on;box on;grid on;xlabel('SNR');ylabel('frequency');axis tight;
plot(SNR,[fcRectWin;fcHann;fcHamming;fcBlackman],'linewidth',1.5);
plot(SNR,[fcRectWinRaw;fcHannRaw;fcHammingRaw;fcBlackmanRaw],'linewidth',1.5);
title(['信号真实频率为',num2str(f0),'Hz']);
legend('RectWin','Hann','Hamming','Blackman','location','southeast');
RMSE = [rmseRectWin;rmseHann;rmseHamming;rmseBlackman];
errRW = fcRectWin-f0;
errHN = fcHann-f0;
errHM = fcHamming-f0;
errBM = fcBlackman-f0;
err = [errRW;errHN;errHM;errBM]/f0*100;
figure;semilogy(SNR,abs(err),'linewidth',1.5);grid on;box on;
xlabel('SNR/dB');ylabel('Relative Error/%');
legend('RectWin','Hann','Hamming','Blackman','location','southeast');
title(['信号真实频率为',num2str(f0),'Hz']);
%% 结论：Rife算法不行，推导过程就是有问题的
代码片

可以看到，真实频率4.9，使用rife算法前估计为4.912，使用后反而变差了。

总结

Rife算法真不行，这是原博客博主的原话

“看情况了，这就当时课设应付任务做的，算法十几年前就已经是过时货了。不加修正的irfe算法如果要使用的话受真实频率影响较大，如果真实频率落在采样两条谱线中间的话相对误差和RMSE都是比较好的，但是如果真实频率落在离谱线比较近的地方误差基本没眼看，就只能调整采样率。”