现代频谱估计的几种方法

最新推荐文章于 2022-09-28 10:23:03 发布
最新推荐文章于 2022-09-28 10:23:03 发布
阅读量1w 收藏 152
点赞数 27
文章标签： 信号处理
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/xdzhujy/article/details/105465773
版权

我要讲的几种方法

MUSIC算法

MUSIC算法简介

  • MUSIC is short for Multiple Signal Classification.
  • MUSIC算法是基于特征结构分析的空间谱估计方法，是空间谱估计技术的典型代表。
  • 其测向原理是根据矩阵特征分解的理论，对阵列输出协方差矩阵进行特征分解，将信号空间分解为噪声子空间G和信号子空间S，利用噪声子空间G与阵列的方向矩阵A的列矢量正交的性质，构造空间谱函数P(w)并进行谱峰搜索，从而估计出波达方向信息。

公式推导

在这里插入图片描述
最终结果可以表示为： p ( w ) = 1 ∥ a ⃗ ( w ) G ∥ 2 , w ∈ [ − π , π ] p(w)=\frac{1}{\left \| \vec{a}(w) G \right \|^{2}},w \in \left [ -\pi,\pi \right ] p(w)=∥a (w)G∥21​,w∈[−π,π]
流程图如下：

Created with Raphaël 2.2.0 开始 估计自相关矩阵R 分解得U和D 组成噪声子空间G 求出p(w) 结束

Matlab仿真

代码1

clc,clear all
w=2*pi*[0.05 0.40 0.42];
times=512;
Nfft=1024;
Nr=31;
figure;
	h=gca; 
	set(h,'FontSize',18); 
    x=linspace(-.5,.5,Nfft);
	ylabel('Power(dB)  ');
	xlabel('Normalized Freqency');
	title( 'Estimated Graph ');
for i=1:times
    smp=fix(100*rand);          %向下取整
    Z=randn(1,128) * sqrt(0.32);
    n=smp:smp+127;
    X=2*cos(w(1)*n)+3*cos(w(2)*n)+1.2*cos(w(3)*n);
    Y=(X+Z)';
    %  generate signal complete!!!
    Rline=zeros(1,255);
    R=zeros(Nr,Nr);
    Rline=xcorr(Y')/128;       %自相关
    R=toeplitz(Rline(128:128+Nr-1));    %托普利兹矩阵的特点是：除度第一行、第一列外，其他每个元知素都与它左上角的元素相同。
    %  correlation matrix complete
    [eigvct,eigval]=eig(R);
    Smusic=zeros(Nfft,1);
    for ii=1:Nr-6
        Smusic=Smusic+abs(fftshift(fft(eigvct(:,ii),Nfft))).^2;
    end
    Smusic=1./Smusic;
    hold on;
    plot(x,10*log10(abs(Smusic)) );
    %  paint complete!!!
    %  search several max value
    Shalf=abs(Smusic(Nfft/2+1:Nfft));
    for ii=1:3
        [C,Ind]=max(Shalf);
        f_estm(i,ii)=(Ind)/Nfft;
        Shalf(Ind-10:Ind+10)=0;
    end
    %   freq estimated  complete!!!
end
%   following progm for mean and var
for ii = 1:3
    mean(f_estm(:,ii));
    var(f_estm(:,ii));
end

结果1

估计的频谱

代码2

clear all;close all;clc;
j=sqrt(-1);
seita=[10,30,53]';
N=size(seita,1);    %求出信号个数
M=11;               %阵子个数
kp=1000;             
lmda=0.24;          %信号波长  
d=0.12;             %阵子的间距
jthresh=1000;       %给出一个比较的较大值
jdbeg=1;
jdend=180;
jdnum=180;
a9=20;
% for m=1:M
%     for n=1:N
%         A(m,n)=exp(j*2*pi*(m-1)*(d/lmda)*sin(seita(n)*pi/180));
%     end
% end

a1=pi/180;
for n=1:N   
    a2=j*pi*sin(seita(n)*a1);         %导向矩阵赋值
    for m=1:M
        A(m,n)=exp(a2*(m-1));
    end
end

a_jd=zeros(1,jdnum);
drt=(jdend-jdbeg)/(jdnum-1);
t1=[0:0.005:25];                      %最大快拍5000个
out1=exp(j*2*pi*10*t1);
out2=exp(j*2*pi*20*t1);               %信号模型
out3=exp(j*2*pi*50*t1); 
cd1=100;
cd2=200;
cd3=400;
S(1,:)=out1(cd1:cd1+kp-1);            %构造S矩阵  
S(2,:)=out2(cd2:cd2+kp-1);
S(3,:)=out3(cd3:cd3+kp-1);           
Ds=A*S(:,1:kp);                       %构造出 A*S
ns=norm(Ds,'fro');
cf=10;
for i1=1:jdnum
    i1;
    a1=jdbeg+(i1-1)*drt;
    a_jd(1,i1)=a1;
    jdpf(1,i1)=0;
    for m=1:M
        a(m,1)=exp(j*2*pi*(m-1)*(d/lmda)*sin(a1*pi/180));
    end
    for ig3=1:cf
        noise1=(randn(M,kp)+j*randn(M,kp));
        nn=norm(noise1,'fro');
        a9=10^(a9/20);
        xk=ns/(nn*a9);
        noise1=noise1*xk;
        Dq=Ds+noise1;                          %添加噪声

        sd=1/(kp);
        R11=sd*Dq(:,1:kp)*Dq(:,1:kp)';         %求相关函数

        [V,D] = eig(R11);       %D对角阵
        [p,k]=sort(diag(D));    %从小到大排列特征值  p特征值 k对应序号
        for i4=1:(M-N)
            Un(:,i4)=V(:,k(i4));
        end
        pf=1/abs(a'*Un*Un'*a);
        jdpf(i1)=jdpf(i1)+pf;
    end
end
        
jdpf=jdpf/cf;
plot(a_jd,jdpf);

结果2

在这里插入图片描述

其他参考

  1. https://wenku.baidu.com/view/88e19ffafab069dc5022018a.html
  2. https://blog.csdn.net/zhangziju/article/details/100730081
  3. https://wenku.baidu.com/view/9b0aa9e60740be1e650e9af5.html

ESPRIT算法

ESPRIT算法简介

  • ESPRIT is short for Estimating Signal Parameter Variational Invariance Techniques (基于旋转不变技术的信号参数估计)。
  • 类似MUSIC算法，也是基于相关矩阵特征分解的信号频率估计方法。
  • 基于ESPRIT有两种经典的方法：最小二乘（LS）法和总体最小二乘（TLS）法。

公式推导

在这里插入图片描述
流程图如下：

Created with Raphaël 2.2.0 开始 构造Rxx和Rxy 分解得最小特征值 对{Cxx,Cxy}做广义特征值分解 求K个特征值相位 结束

Matlab仿真

代码1：ESPRIT_EIG

clc,clear all,close all
Nr=7;
w=2*pi*[0.05 0.40 0.42];
N_Var=length(0.32:.01:6);
Z=zeros(Nr,Nr);
Z(2:Nr,1:Nr-1)=blkdiag(1,1,1,1,1,1);
I=blkdiag(1,1,1,1,1,1,1);
valuediag=zeros(20,3);
frq_mean=zeros(N_Var,3);
frq_var=zeros(N_Var,3);
for k=1:N_Var
    for i=1:20
        smp=fix(100*rand);
        Nise=randn(1,128)*sqrt(0.32+(k-1)*.01);
        n=smp:smp+127;
        X=2*cos(w(1)*n)+3*cos(w(2)*n)+1.2*cos(w(3)*n);
        Y=(X+Nise)';
        Rline=xcorr(Y')/128;
        Rxx=toeplitz(Rline(128:128+Nr-1));
        Rxy1=toeplitz(Rline(128:128+Nr   ));
        Rxy=Rxy1(1:Nr,2:Nr+1);
        [eigvct,eigval]=eig(Rxx);
        lamdamin=eigval(1,1);
        Cxx=Rxx-lamdamin*I;
        Cxy=Rxy-lamdamin*Z;
        [eigvct,eigval]=eig(Cxx,Cxy);
        temp=sort(unique(abs(angle(diag(eigval)')/2/pi)));
        if temp(1)==0
            temp(1:3)=temp(2:4);
        end
        valuediag(i,1:3)=temp(1:3);
    end
    frq_mean(k,:)=mean(valuediag);
    frq_var(k,:)=var(valuediag);
end

figure;
h=gca;
set(h,'FontSize',25);
x=0.32:.01:6;
subplot(2,1,1);
plot(x,frq_mean(:,1),'-.r',x,frq_mean(:,2),'b',x,frq_mean(:,3),'k.');
ylabel('Mean Estimated Frequency  ');
xlabel('Noise VAR');
title( 'Mean Frequency trends for ESPRIT in book ');
legend('0.05','0.40','0.42');
subplot(2,1,2);
plot(x,frq_var(:,1),'-.r',x,frq_var(:,2),'b',x,frq_var(:,3),'k.');
ylabel('Variance Estimated Frequency  ');
xlabel('Noise VAR');
title( 'Variance Frequency trends for ESPRIT in book');
legend('0.05','0.40','0.42');

结果1

在这里插入图片描述

代码2：ESPRIT_TLS

clc,clear all,close all
Nrow=112; %  number of row
Ncolum=17;  %number  of colum
Ntimes=20;   % iteration times
N_Var=length(0.32:.01:6);
w=2*pi*[0.05 0.40 0.42];  % normalized frq
P=length(w);          %  amount of harmonic
Xmatr=zeros(Nrow,Ncolum); % data matrix
frq_tls=zeros(Ntimes,3);   %  frq store
valid_N=0;
frq_mean=zeros(N_Var,3);
frq_var=zeros(N_Var,3);
for k=1:N_Var
    for i=1:Ntimes
        % to generate signal with Gauss Noise
        smp=fix(100*rand);
        Nise=randn(1,128) * sqrt(0.32+(k-1)*.01);
        n=smp:smp+127;
        X=2 * cos( w(1)*n )+3 * cos( w(2)*n )  + 1.2 * cos( w(3)*n );
        Y=(X + Nise);
        % to get data matrix
        for ii=1:Nrow
            Xmatr(ii,1:Ncolum)=Y(ii:ii+Ncolum-1);
        end
        % use SVD and TLS to estimate frq
        [L,S,U]=svd(Xmatr);
        Us=U(:,1:6);
        U1=Us(1:7,:);U2 = Us(2:8,:);
        [LL,SS,UU]=svd([U1 U2]);
        UU12=UU(1:6,7:12);
        UU22=UU(7:12,7:12);
        psi=-UU12/UU22;
        phi=sort(unique( abs(angle(eig(psi)))));
        if length(phi) ~= 3
            break;
        end
        valid_N=valid_N+1;
        frq_tls(valid_N,:)=phi/(2*pi);
    end
    % to get mean and Var
    frq_mean(k,:)=mean(frq_tls(1:valid_N,:));
    frq_var(k,:)=var(frq_tls(1:valid_N,:));
end

% to evaluate the performance for ESPRIT_TLS
figure;
h=gca; %
set(h,'FontSize',25); %
x=0.32:.01:6;
subplot(2,1,1);
plot(x,frq_mean(:,1),'-.r',x,frq_mean(:,2),'b',x,frq_mean(:,3),'k.');
ylabel('Mean Estimated Frequency  ');
xlabel('Noise VAR');
title( 'Mean Frequency trends ');
legend('0.05','0.40','0.42');
subplot(2,1,2);
plot(x,frq_var(:,1),'-.r',x,frq_var(:,2),'b',x,frq_var(:,3),'k.');
ylabel('Variance Estimated Frequency  ');
xlabel('Noise VAR');
title( 'Variance Frequency trends ');
legend('0.05','0.40','0.42');

结果2

在这里插入图片描述

其他参考

  1. https://blog.csdn.net/ai136172022/article/details/80963764
  2. https://wenku.baidu.com/view/a03a6c1e54270722192e453610661ed9ad515519.html

Pisarenko算法

Pisarenko算法简介

  • Pisarenko算法是由信号分析学家Pisarenko于1970s提出，是一种空间谱估计。
  • 该方法首先通过把特征多项式的系数向量和自相关矩阵特征向量结合估计频率然后利用得到的频率采用自相关算法估计振幅和相位。
  • 该方法有很多改进算法。
  • 本文介绍谐波恢复的Pisarenko。

公式推导

在这里插入图片描述
流程图如下：

Created with Raphaël 2.2.0 开始 构造Ry(k) 根据上图最后一个式子估算系数向量 通过求得的a估算特征函数中的共轭根 求得每个信号分量的频率估计值 fi=arctan[Im(zi)/Re(zi)]/2pi 结束

Matlab仿真

代码

clc,clear all,close all
w=2*pi*[0.05 0.40 0.42];
varia=0.32;
i=1;
k=1;
for varia=0.32:0.01:6
    i=1;
    f_estmt=zeros(20,3);
    for i=1:20
        Rline=zeros(1,255);
        R=zeros(7,7);
        smp=fix(100*rand);
        Z=randn(1,128)*sqrt(varia);
        a=zeros(1,4);
        n=smp:smp+127;
        X=2*cos(w(1)*n)+3*cos(w(2)*n)+1.2*cos(w(3)*n);
        Y=(X+Z)';
        Rline=xcorr(Y');
        for ii=1:7
            R(ii,:)=Rline(129-ii:135-ii);
        end
        [eigvct,eigval]=eig(R);
        [C,Ind]=min(diag(eigval));
        z=roots(eigvct(:,Ind));
        longz=length(unique(abs(angle(z))/2/pi));
        a(1:longz)=sort(unique(abs(angle(z))/2/pi))';
        maxtop3=sort(a);
        %gettop3maxvalue
        f_estmt(i,1:3)=maxtop3(2:4);
    end
    meanS(k,1:3)=mean(f_estmt);
    VarS(k,1:3)=var(f_estmt);
    k=k+1;
end

figure;
subplot(2,1,1)
h=gca;
set(h,'FontSize',16);
x=linspace(0.32,6,569);
plot(x,meanS(:,1),'-.r',x,meanS(:,2),'b',x,meanS(:,3),'k+');
ylabel('Mean Estimated Frequency');
xlabel('Noise VAR');
title('Mean Frequency trends');
legend('0.05','0.40','0.42');

subplot(2,1,2)
h=gca;%
set(h,'FontSize',16);%
x=linspace(0.32,6,569);
plot(x,VarS(:,1),'-.r',x,VarS(:,2),'b',x,VarS(:,3),'k+');
ylabel('Variance of Estimated Frequency');
xlabel('NoiseVAR');
title('Frequency Variance trends');
legend('0.05','0.40','0.42');

结果

在这里插入图片描述

其他参考

  1. http://www.doc88.com/p-6252502979147.html
  2. 基于 Hilbert 变换与 Pisarenko 谐波分解的电压闪变参数估计
  3. 一种自相关矩阵改进的 Pisarenko 算法

结论

  • 由仿真图形可以看出，MUSIC算法、ESPRIT算法等都可以实现对含噪复正弦信号的频率估计，而且能够克服 DFT 中存在能量泄漏和栅栏效应，误差较小。
  • MUSIC算法最为常规，而且能够实现超分辨，有效的克服了工程应用中由于先验信息不足而导致的分辨率降低问题，但是运算量也是很大，不利于次数较大的频率估计。
  • ESPRIT算法需要两次求特征值运算，实现较为复杂，但是有效的克服MUSIC算法需要进行谱峰搜索而带来的计算量很大的问题，计算量很小，而且随着运算次数的增大，运算时间不会明显增大，具有很好的分辨力。ESPRIT算法能够有效降低运算量，不过，在较高频率的频率估计时有少许误差。
  • 经典Pisarenko算法在对自相关矩阵进行降维的过程中存在误差，谐波恢复Pisarenko算法使用自相关函数过少也存在误差，因此存在很多改进的算法提高这一点。

其他

  • 本文缺少了对具体各个误差的分析，若有足够时间，可以在后续进行。
  • Prony、MVDR等不再进行叙述，请大家自行了解，若有错误请同学们指正！谢谢！
我亦知难而退
关注
  • 27
    点赞
  • 踩
  • 152
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 4
    评论
古典法的频谱估计MATLAB程序
12-26
这是用MATLAB编写的古典法，经典法的频谱估计；也是随机信号课程大作业之一。望对你有帮助！
现代信号处理——其他谱估计方法（最大熵谱估计）
qq_42233059的博客
09-29 3258
现代信号处理——其他谱估计方法（最大熵谱估计）
matlab仿真-谱估计_matlab_频谱估计_自相关频率_
10-01
基于matlab的自相关算法等的频率谱估计
频率估计常用方法
06-03
文中给出频率估计的常用方法，如周期图法，最大似然估计，子空间方法，加权相平均等等
jaoten.zip_esprit频谱估计
07-14
esprit算法对有干扰的信号频率进行估计,鲁棒性好，性能优越，对信号进行频谱分析及滤波。
esprit算法研究，实现空间谱估计算法，并考察算法性能。7000字word版
03-24
1） 由均匀线阵形式，确定阵列的导向矢量； 2） 由阵列导向矢量，对接收信号进行建模仿真； 3） 由ESPRIT算法实现信号DOA估计； 4） 考察算法性能与信噪比，采样率，观测时间等参数的关系。
几种常用的经典频谱估计方法原理介绍及Matlab算法实现
热门推荐
Inspire
01-27 5万+
在信号的谱分析中，由于经典的谱分析方法物理意义比较明确，能够处理的信号信噪比可以很低（相对于现代的谱估计来说），而且处理算法简单稳定；因此在大多数的情况下都是使用经典的谱分析方法对信号进行谱分析
现代信号处理——经典谱估计
qq_42233059的博客
09-28 7221
现代信号处理——经典谱估计
【经典谱和现代谱估计】
qq_45698530的博客
06-27 2020
经典谱 现代谱
matlab 最大熵谱估计,基于Burg算法的最大熵谱估计.doc
weixin_39817122的博客
03-23 507
基于Burg算法的最大熵谱估计基于Burg算法的最大熵谱估计其中，P(ejω)为功率谱密度；s2是激励白噪声的方差；ai和bi为模型参数。若ARMA 模型中bi全为0，就变成了AR模型，又称线性自回归模型，其是一个全极点模型：研究表明，ARMA 模型和MA 模型均可用无限阶的AR模型来表示。且AR模型的参数估计计算相对简单。同时，实际的物理系统通常是全极点系统。要利用AR模型进行功率谱估计，必须由...
基于Matlab的经典频谱估计周期图方法的仿真研究。
02-24
基于Matlab的经典频谱估计周期图方法的仿真研究。
浅谈频谱仪测量噪声系数方法
10-22
由于放大器本身就有噪声，输出端的信噪比和输入端信噪比是不一样的，为此，使用噪声系数来衡量放大器本身的噪声水平。该系数并不是越大越好，它的值越大，说明在传输过程中掺入的噪声也就越大，反应了器件或者信道特性的不理想。 　　用频谱仪测量噪声系数：    　　测量框图为： 　　为了能够准确测量，将频谱仪的读数置为多次平均方式，考虑到测量时间和测量准确度，一般为十五次左右。 　　基于噪声系数的定义得到的一个测量公式为： 　　NF=PNOUT-（-174dBm/Hz+20lg（BW）+Gain）     （1） 　　公式中， PNOUT是已测的总共输出噪声功率，-174dBm/Hz是290oK（室温）时
vafdvbch.zip_频谱 机器学习_频谱估计_频谱相位估计
07-14
可以得到很精确的幅值、频率、相位估计，保证准确无误，是学习通信的好帮手，使用大量的有限元法求解偏微分方程，是机器学习的例程，研究生时的现代信号处理的作业，对信号进行频谱分析及滤波。
hires_spectrum-master_Spectrum_高质量频谱估计_
10-01
Matlab代码已经过轻微修改，可以在现代版本的Matlab上运行。
一种新的认知无线电宽带盲频谱感知方法-综合文档
05-22
一种新的认知无线电宽带盲频谱感知方法
频谱校正法,频谱仪校准方法,matlab
09-10
关于雷达中频信号的频率较正方法的源程序，供大家参考
基于程序频谱的缺陷定位方法
05-06
基于此，本文提出EPStar（EP*）缺陷定位方法，该方法可以有效调整成功执行用例数的影响，以避免成功用例数量对缺陷定位效果的过度影响，从而提高缺陷定位的准确性，通过实验对比，说明了EP*方法比现有的几种缺陷定位...
快速傅里叶变化FFT频谱估计
04-18
使用快速傅里叶变换FFT实现带有高斯白噪声的信号的频谱估计
噪声中正弦信号的现代法频谱分析
12-07
噪声中正弦信号的现代法频谱分析
语音信号的频谱估计方法
最新发布
05-18
语音信号的频谱估计方法有很多种，下面列举几种常用的方法： 1. 傅里叶变换法：将语音信号分解成一系列正弦波的和，求得每个正弦波的振幅和相位，从而得到频谱信息。 2. 短时傅里叶变换法：将语音信号分成若干个短时段，对每个短时段进行傅里叶变换，得到每个时段的频谱信息。 3. 希尔伯特变换法：将语音信号转化为解析信号，即对原信号进行希尔伯特变换，得到包含原信号全部频率信息的解析信号，再对解析信号进行傅里叶变换，得到频谱信息。 4. 自相关函数法：计算语音信号的自相关函数，求得自相关函数的峰值点，从而得到频率信息。 5. LPC（线性预测编码）法：通过对语音信号进行线性预测分析，求得语音信号的LPC系数，再通过LPC系数计算得到语音信号的频谱信息。 这些方法各有优缺点，应根据具体情况选择合适的方法。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

我亦知难而退

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值