我要讲的几种方法
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∈[−π,π]
流程图如下:
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
其他参考
- https://wenku.baidu.com/view/88e19ffafab069dc5022018a.html
- https://blog.csdn.net/zhangziju/article/details/100730081
- https://wenku.baidu.com/view/9b0aa9e60740be1e650e9af5.html
ESPRIT算法
ESPRIT算法简介
- ESPRIT is short for Estimating Signal Parameter Variational Invariance Techniques (基于旋转不变技术的信号参数估计)。
- 类似MUSIC算法,也是基于相关矩阵特征分解的信号频率估计方法。
- 基于ESPRIT有两种经典的方法:最小二乘(LS)法和总体最小二乘(TLS)法。
公式推导
流程图如下:
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
其他参考
- https://blog.csdn.net/ai136172022/article/details/80963764
- https://wenku.baidu.com/view/a03a6c1e54270722192e453610661ed9ad515519.html
Pisarenko算法
Pisarenko算法简介
- Pisarenko算法是由信号分析学家Pisarenko于1970s提出,是一种空间谱估计。
- 该方法首先通过把特征多项式的系数向量和自相关矩阵特征向量结合估计频率然后利用得到的频率采用自相关算法估计振幅和相位。
- 该方法有很多改进算法。
- 本文介绍谐波恢复的Pisarenko。
公式推导
流程图如下:
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');
结果
其他参考
- http://www.doc88.com/p-6252502979147.html
- 基于 Hilbert 变换与 Pisarenko 谐波分解的电压闪变参数估计
- 一种自相关矩阵改进的 Pisarenko 算法
结论
- 由仿真图形可以看出,MUSIC算法、ESPRIT算法等都可以实现对含噪复正弦信号的频率估计,而且能够克服 DFT 中存在能量泄漏和栅栏效应,误差较小。
- MUSIC算法最为常规,而且能够实现超分辨,有效的克服了工程应用中由于先验信息不足而导致的分辨率降低问题,但是运算量也是很大,不利于次数较大的频率估计。
- ESPRIT算法需要两次求特征值运算,实现较为复杂,但是有效的克服MUSIC算法需要进行谱峰搜索而带来的计算量很大的问题,计算量很小,而且随着运算次数的增大,运算时间不会明显增大,具有很好的分辨力。ESPRIT算法能够有效降低运算量,不过,在较高频率的频率估计时有少许误差。
- 经典Pisarenko算法在对自相关矩阵进行降维的过程中存在误差,谐波恢复Pisarenko算法使用自相关函数过少也存在误差,因此存在很多改进的算法提高这一点。
其他
- 本文缺少了对具体各个误差的分析,若有足够时间,可以在后续进行。
- Prony、MVDR等不再进行叙述,请大家自行了解,若有错误请同学们指正!谢谢!