一种对真正正弦波频率进行准确估计的方法（Matlab代码实现）

 💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码、文章

---

💥1 概述

众所周知，真实正弦波的正频率和负频率分量在频谱上相互作用，导致基于周期图最大化的频率估计出现偏差。我们提出了滤除负频率分量的方法。为此，首先利用窗函数方法获得粗略频率估计，已知这一方法可以减少估计偏差，然后利用调制和离散傅立叶变换箱剔除负频率分量。接着，对经过滤波的信号使用专门为复正弦波开发的准确频率估计器进行精细估计。该方法具有O(Nlog2N)的复杂度（加法/乘法）和O(N)的复杂度（正弦/余弦运算和比较）。此外，该方法达到了克拉美尔-拉奥下界，并且对正弦波的频率和初始相位不敏感，因此表现优于现有方法。

正弦波频率估计在许多应用中都是一个基本问题，包括频谱估计、卫星/移动通信、阵列信号和雷达信号处理，以及多项式相位信号估计。在白色高斯噪声环境中，单个复正弦波的最大似然（ML）频率估计是接收信号周期图的峰值位置。用于正弦波频率估计的最流行算法基于通过在频率域中进行某种插值来最大化周期图。换句话说，周期图是通过离散傅立叶变换（DFT）计算的，由于DFT的离散性质，DFT最大位置通常不对应真实周期图的最大值。因此，需要插值来减少这两个最值之间的位移。插值是通过计算DFT最大值周围的几个额外点或者结合已计算的DFT样本来实现的。

提出的大多数正弦波频率估计方法涉及复正弦波。可以说，这个方向的研究已经达到了饱和点，因为无论在精度还是计算复杂度方面，现有方法的改进都很小。仅通过计算离散傅立叶变换网格外的几个额外点就可以达到克拉美尔-拉奥下界。然而，当处理实值（RV）正弦波时，周期图最大化方法受到正频率和负频率的复正弦波（欧拉公式）的频谱叠加，这会引入估计偏差。

📚2 运行结果

部分代码：

function pa=proposed(N,x,L,Q)%Length,signal,number of sinusoids, number of iterations
n=0:N-1;        %if you have any questions, please contact mzlaccount@163.com
s=x;
for k=1:L  %coarse estimation
    Y=fft(s);
    [~, m2]=max(abs(Y));
    if m2==1
        m2=2;
    end
    f0=(m2-1)/N;     %粗估计w0
    delta=0;S2=0;
    for Q=1:2
        b=round(1/4/f0);
        if b<=0             %in low SNRs the b is unstable
            b=1;
        end
        if b>(N+1)/2        
            b=round(N/2);
        end
        h=cos(2*pi*f0*b);   %参考信号 （8）
        h2=sin(2*pi*f0*b);  %正交分量 （8）  
        for i=1:(N-b)
            S2(i)=(s(i)*h-s(i+b))/h2;%-(z(i)*cos(2*pi*f0*b)-z(i+b))/sin(2*pi*f0*b)-z(i+b);+2*z(i+b)/sin(2*pi*f0*b);%
        end 
        for i=(N-b+1):N
            S2(i)=(s(i-b)-s(i)*h)/h2;%-(z(i-b)-z(i)*cos(2*pi*f0*b))/sin(2*pi*f0*b)+z(i-b);
        end
        SC=s+1j*S2;
        xp=sum(SC.*exp(-1i*2*pi*n*(m2-1+delta+0.5)/N));
        xp2=sum(SC.*exp(-1i*2*pi*n*(m2-1+delta-0.5)/N));
        delta=delta+0.5*(abs(xp)-abs(xp2))/(abs(xp)+abs(xp2));
        f0=(m2-1)/N+delta/N; 
    end
    A=sum(SC.*exp(-1j*f0*n*2*pi))/N;
    AA=abs(A);
    theta=angle(A);
    pa(1,:)=f0;
    pa(2,:)=AA;
    pa(3,:)=theta;
    s=s-AA*cos(2*pi*f0*n+theta);
    fatf(:,k)=pa;
end
%fine estimation
for qq=1:Q
    for k=1:L
        s=x;
        for q=1:L
            if q==k
            else
                s=s-fatf(2,q)*cos(2*pi*fatf(1,q)*n+fatf(3,q));
            end
        end     
            b=round(1/4/fatf(1,k));
            if b<=0
                b=1;
            end
            if b>(N+1)/2
                b=round(N/2);
            end
            h=cos(2*pi*fatf(1,k)*b);   %参考信号 （8）
            h2=sin(2*pi*fatf(1,k)*b);  %正交分量 （8）  
            for i=1:(N-b)
                S2(i)=(s(i)*h-s(i+b))/h2;
            end 
            for i=(N-b+1):N
                S2(i)=(s(i-b)-s(i)*h)/h2;
            end
            SC=s+1j*S2;
            xp=sum(SC.*exp(-1i*2*pi*n*(fatf(1,k)+0.5/N)));
            xp2=sum(SC.*exp(-1i*2*pi*n*(fatf(1,k)-0.5/N)));
            fatf(1,k)=fatf(1,k)+0.5*(abs(xp)-abs(xp2))/(abs(xp)+abs(xp2))/N;
        AA=sum(SC.*exp(-1j*fatf(1,k)*n*2*pi))/N;
        fatf(2,k)=abs(AA);fatf(3,k)=angle(AA);
    end
end

pa=sortrows(fatf',1)';

🎉3 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。

🌈4 Matlab代码、文章