mrife.m

M-rife估频算法的实现

%%%%%%% -------------------------------------------------------------------
%%%%% -----修正Rife算法用于频率估计
%
function [fe mi] = mrife(x,fs,N)
%
% x:输入数据
% fs:数据采样率
% N:FFT的点数
% fe:估计出的频率
% mi:最长谱线的频域采样系数

% N = length(x);

delta_f = fs/N;
n = 0:N-1;   %%% ?

Z = fft(x,N);
P = abs(Z);
[Pm mi] = max(P);
fo = fs*(mi-1)/N;

if mi == 1 || P(mi+1) >= P(mi-1)
    r = 1;
else
    r = -1;
end

fe = ((mi-1)+r*P(mi+r)/(P(mi)+P(mi+r)))*fs/N;

if abs(fe-fo) > delta_f/3 && abs(fe-fo) <= delta_f/2    %%% ---落在两采样频率的中心区域？
    fe = fe;
else
    shift = 1/2 - P(mi+r)/(P(mi)+P(mi+r));
%     shift = 1/3;
%     x1 = x.*exp(j*2*pi*n*r*shift/N);
    Pr = abs(x.*exp(-j*2*pi*n*(mi+r-1-r*shift)/N));
    Pm = abs(x.*exp(-j*2*pi*n*(mi-1-r*shift)/N));
    fe = ((mi-1)+r*min(Pm,Pr)/(Pm+Pr))*fs/N-r*shift*delta_f;
end