Rife算法的一种Matlab实现

66 阅读 0 评论 44 点赞

我是靠谱客的博主伶俐可乐，最近开发中收集的这篇文章主要介绍Rife算法的一种Matlab实现，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

Rife算法原理及介绍

信号频谱的最大两根谱线进行插值对频率进行估计,即Rife算法。当信号频率位于离散傅里叶变换(DFT)两个相邻量化频率点的中心区域时，Rife算法精度很高,均方根误差接近克拉美罗限,且计算量小于牛顿迭代。Rife算法常用于单频信号频率估计。

加性高斯白噪声污染的正弦波信号可以表示为

$e^{j H_{0}} e^{j 2 pi f_{c} n Delta t}+w(n) quad n=0,1, ldots N-1$

式中a , $f_{c}, mathrm{H}_{0}$ 分别为振幅，频率和初相；∆t 为采样间隔，N 为样本数；w(t)为实部和虚部相互独立的高斯白噪声。
对x(n)作FFT,取其中的最大谱线值|X(K_0)|。给出利用两根谱线进行正弦波频率估计的计算公式(即Rife算法)
$widehat{f}_{c}=frac{1}{T}left(k_{0}+r frac{left|Xleft(K_{0}+rright)right|}{left|Xleft(K_{0}right)right|+left|Xleft(K_{0}+rright)right|}right)$
式中：当
$left|Xleft(K_{0}+1right)right| leqleft|Xleft(K_{0}-1right)right|$ 时，r = -1；
当 $left|Xleft(K_{0}-1right)right| leqleft|Xleft(K_{0}+1right)right|$ 时，r = 1；T = NΔt。[1]

Matlab示例程序

CSDN目前不支持Matlab的格式，建议复制到Matlab或VScode等编辑器中查看

%%
clc;clear;close all;
%%  参数设置

N = 1000;
n = 1:N;
fs = 16;
f = 4.9;        %可以在此选择模拟信号频率
T = N/fs;
dt = 1/fs;

SNR_n = 20;
SNR = -1*SNR_n+2:2:SNR_n-2;SNR = SNR+8;
%% 信号生成

signal_0 = cos(2*pi*f/fs*n)';
signal = zeros(N,SNR_n);
for k = 1:SNR_n-1
    signal(:,k+1) = awgn(signal_0,SNR(k));      %高斯白噪声
end
signal(:,1) = signal_0;
%% 窗函数

rect = rectwin(N);
hn = hann(N);
hm = hamming(N);
bk = blackman(N);

%% 加窗
signalRectWin = repmat(rect,1,SNR_n).*signal;   %将窗函数复制，与信号向量做点积
signalHann = repmat(hn,1,SNR_n).*signal;
signalHamming = repmat(hm,1,SNR_n).*signal;
signalBlackman = repmat(bk,1,SNR_n).*signal;

%%  FFT
wRectWin = real(fftshift(fft(signalRectWin),1));
wHann = real(fftshift(fft(signalHann),1));
wHamming = real(fftshift(fft(signalHamming),1));
wBlackman = real(fftshift(fft(signalBlackman),1));

x = [-N+1:2:N]/N*fs/2;
%plot(x,db(wRectWin));

%% Rife算法
%%%%%%%%%%%%%%%%%%%%%%%%%%
[valueMaxRectWin,posMaxRectWin] = max(wRectWin(N/2+1:N,:));
r = zeros(1,SNR_n);
for k = 1:SNR_n
    r(k)= 2*((wRectWin(posMaxRectWin(k)+N/2+1,k) > wRectWin(posMaxRectWin(k)+N/2-1,k))-0.5);
    rat(k) = valueMaxRectWin(k)/(valueMaxRectWin(k)+wRectWin(posMaxRectWin(k)+N/2+r(k),k));
    fcRectWin(k) = 1/T*(posMaxRectWin(k)+r(k)*rat(k)-1);
end


rmseRectWin = sqrt(abs(sum(fcRectWin-f)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%

[valueMaxHann,posMaxHann] = max(wHann(N/2+1:N,:));
r = zeros(1,SNR_n);
for k = 1:SNR_n
    r(k)= 2*((wHann(posMaxHann(k)+N/2+1,k) > wHann(posMaxHann(k)+N/2-1,k))-0.5);
    rat(k) = valueMaxHann(k)/(valueMaxHann(k)+wHann(posMaxHann(k)+N/2+r(k),k));
    fcHann(k) = 1/T*(posMaxHann(k)+r(k)*rat(k)-1);
end


rmseHann = sqrt(abs(sum(fcHann-f)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%
[valueMaxHamming,posMaxHamming] = max(wHamming(N/2+1:N,:));
r = zeros(1,SNR_n);
for k = 1:SNR_n
    r(k)= 2*((wHamming(posMaxHamming(k)+N/2+1,k) > wHamming(posMaxHamming(k)+N/2-1,k))-0.5);
    rat(k) = valueMaxHamming(k)/(valueMaxHamming(k)+wHamming(posMaxHamming(k)+N/2+r(k),k));
    fcHamming(k) = 1/T*(posMaxHamming(k)+r(k)*rat(k)-1);
end

rmseHamming = sqrt(abs(sum(fcHamming-f)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%
[valueMaxBlackman,posMaxBlackman] = max(wBlackman(N/2+1:N,:));
r = zeros(1,SNR_n);
for k = 1:SNR_n
    r(k)= 2*((wBlackman(posMaxBlackman(k)+N/2+1,k) > wBlackman(posMaxBlackman(k)+N/2-1,k))-0.5);
    rat(k) = valueMaxBlackman(k)/(valueMaxBlackman(k)+wBlackman(posMaxBlackman(k)+N/2+r(k),k));
    fcBlackman(k) = 1/T*(posMaxBlackman(k)+r(k)*rat(k)-1);
end


rmseBlackman = sqrt(abs(sum(fcBlackman-f)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 误差分析
RMSE = [rmseRectWin;rmseHann;rmseHamming;rmseBlackman];
errRW = fcRectWin-f;
errHN = fcHann-f;
errHM = fcHamming-f;
errBM = fcBlackman-f;
err = [errRW;errHN;errHM;errBM]/f*100;
plot(SNR,abs(err(:,2:SNR_n)));
grid on;
xlabel('SNR/dB');
ylabel('Relative Error/%');
leg = {'RectWin','Hann','Hamming','Blackman'};
legend(leg);
str = strcat('信号真实频率为',num2str(f),'Hz');
title(str);

为方便使用，将Rife算法独立编写成函数，代码如下：

function fc = rife(spec,fs)
%Rife算法matlab实现
%   spec:信号频谱列向量或以列向量叠加的矩阵
%   fs:信号采样率
%   输出fc为Rife算法估计出的频率
[length,SignalNum] = size(spec);
[valueMax,posMax] = max(spec(length/2+1:length,:));
T = length/fs;
for k = 1:SignalNum
    r= 2*((spec(posMax(k)+length/2+1,k) > spec(posMax(k)+length/2-1,k))-0.5);
    rat = valueMax(k)/(valueMax(k)+spec(posMax(k)+length/2+r,k));
    fc(k) = 1/T*(posMax(k)+r*rat-1);
end
end

另附试验用生成信号的函数SignalFormFFT

function [wOUT,SNR] = signalFormFFT(N,f,fs,SNR_n)
%生成模拟信号以及加窗与FFT
%   N为采样点数，f为信号频率，fs为采样率，SNR_n为信噪比个数，对应计算出的信噪比序列为
%    -1*SNR_n+2:2:SNR_n-2;SNR = SNR+8;
%   输出元胞的顺序分别是加矩形窗，汉宁窗，海明窗和布拉克曼窗
%%  参数设置

n = 1:N;
T = N/fs;
dt = 1/fs;
SNR = -1*SNR_n+2:2:SNR_n-2;SNR = SNR+8;
%% 信号生成

signal_0 = cos(2*pi*f/fs*n)';
signal = zeros(N,SNR_n);
for k = 1:SNR_n-1
    signal(:,k+1) = awgn(signal_0,SNR(k));      %高斯白噪声
end
signal(:,1) = signal_0;
%% 窗函数

rect = rectwin(N);
hn = hann(N);
hm = hamming(N);
bk = blackman(N);

%% 加窗
signalRectWin = repmat(rect,1,SNR_n).*signal;   %将窗函数复制，与信号向量做点积
signalHann = repmat(hn,1,SNR_n).*signal;
signalHamming = repmat(hm,1,SNR_n).*signal;
signalBlackman = repmat(bk,1,SNR_n).*signal;

%%  FFT
wRectWin = real(fftshift(fft(signalRectWin),1));
wHann = real(fftshift(fft(signalHann),1));
wHamming = real(fftshift(fft(signalHamming),1));
wBlackman = real(fftshift(fft(signalBlackman),1));

wOUT = {wRectWin,wHann,wHamming,wBlackman};
end

示例程序：

clc;clear;close all;
%% 原始数据输入
N = 500;
f = 4.08;
fs = 16;
SNR_n = 20;

%% 生成信号
[W,SNR] = signalFormFFT(N,f,fs,SNR_n);
wr = W{1};
whn = W{2};
whm = W{3};
wb = W{4};

%% Rife算法估算频率
fcr = rife(wr,fs);
fchn = rife(whn,fs);
fchm = rife(whm,fs);
fcb = rife(wb,fs);


fc = [fcr;fchn;fchm;fcb];
%% 分析误差
err = (fc-f)/f*100;

plot(SNR,abs(err(:,2:SNR_n)));
grid on;
xlabel('SNR/dB');
ylabel('Relative Error/%');
leg = {'RectWin','Hann','Hamming','Blackman'};
legend(leg);
str = strcat('信号真实频率为',num2str(f),'Hz');
title(str);
disp(str);

附下载链接
Rife.m