概述
本位转载自新浪博客文章，原文地址：http://blog.sina.com.cn/s/blog_1691c09ee0102x9kd.html
空间谱估计的经典算法（MATLAB程序）
标签： matlab 空间谱估计 分类： 编程笔记
1、经典MUSIC算法
程序1
%classic music
clear all
clc
tic
%参数设定
M=10;
doa=[-10 45 60]/180pi;
P=length(doa);
f=1000;c=1500;
lambda=f/c;d=lambda/2;
snr=5;
N=400;
%阵列流型A
for i=1:P
A(:,i)=exp(-j2pid*[0:M-1]’/lambdasin(doa(i)));
end
%信源模型建立
for k=1:P
S(k,:)=sqrt(10.^(snr/10))(sind(1:N)+rand(1,N));
end
%接收信号模型建立
X=AS;
%协方差矩阵特征值分解得到噪声子空间
R=XX’/N;
[V,D]=eig®;
[Y,I]=sort(diag(D));
Un=V(:,I(1:M-P));
%谱峰搜索部分
theta=-90:0.1:90;%线阵的搜索范围为-90~90度
for i=1:length(theta)
a_theta=exp(-j*(0:M-1)‘2pidsin(pitheta(i)/180)/lambda);
Pmusic(i)=1./abs((a_theta)'UnUn’a_theta);
end
plot(theta,10log(Pmusic/max(Pmusic)),‘r’);
axis([-90 90 -90 10]);
xlabel(‘theta/degree’);
ylabel(‘归一化空间谱/dB’);
legend(‘Music Spectrum’);
title(‘经典MUSIC估计’);
grid on;
toc
程序2
%classic music
clear all
clc
tic
%参数设定
M=10;
doa=[-10 45 60]/180pi;
P=length(doa);
f=1000;c=1500;
lambda=f/c;d=lambda/2;
snr=5;
N=128;
%阵列流型A
for i=1:P
A(:,i)=exp(-j2pid[0:M-1]’/lambdasin(doa(i)));
end
%信源模型建立
for k=1:P
S(k,:)=sqrt(10.^(snr/10))(randn(1,N)+jrandn(1,N));
end
%接收信号模型建立
X=AS+1/sqrt(2)(randn(M,N)+jrandn(M,N));
%协方差矩阵特征值分解得到噪声子空间
R=XX’/N;
[V,D,U]=svd®;
US=V(:,1:P);
%谱峰搜索部分
theta=-90:0.1:90;
I=eye(M);
for ii=1:length(theta)
a=exp(-j2pid*(0:M-1)‘sin(theta(ii)/180pi)/lambda);
Pmusic(ii)=abs(1/(a’(I-USUS’)a));
end
Pmusic=10log10(Pmusic/max(Pmusic));
figure(1)
plot(theta,Pmusic,‘r’)
legend(‘classic music’)
xlabel(‘theta/degree’)
ylabel(‘归一化Pmusic/dB’)
title(‘classic music spectrum’)
grid on
toc
2、基于波束空间的MUSIC算法
clear all
clc
tic
%%%%%%%%%信号源的参数设定%%%%%%%%
M=8;
N=256;
doa=[-30 20 50]/180pi;
P=length(doa);
f=100;c=1500;
lambda=c/f;d=lambda/2;
snr=10;
B=5;
%%%%%%%%%波束变换矩阵%%%%%%%%
T=1/sqrt(M)[exp(-j*[0:M-1]'pisin(2/M)) exp(-j*[0:M-1]'pisin(4/M)) exp(-j*[0:M-1]'pisin(6/M)) exp(-j*[0:M-1]'pisin(8/M)) exp(-j*[0:M-1]‘pisin(10/M))];
T1=abs(T’T)
%%%%%%%%%导向矢量A%%%%%%%%%%%%
for i=1:P
A(:,i)=exp(-j2pidsin(doa(i))[0:M-1]/lambda);
end
%%%%%%%%接收信号%%%%%%%%%%
S=10*(randn(P,N)+jrandn(P,N));
X=awgn(AS,snr);
%%%%%%%%求协方差矩阵并分解%%%%%%%%
Y=T’X;
R=YY’/N;
[V,D]=eig®;
UN=V(:,1:B-P);
%%%%%%%%谱峰搜索%%%%%%%%
theta=[-90:0.5:90];
for ii=1:length(theta)
a=exp(-j2pi*[0:M-1]'dsin(theta(ii)/180pi)/lambda);
Pmusic(ii)=1/abs(a’TUNUN’T’a);
end
Pmusic=10log(Pmusic/max(Pmusic));
plot(theta,Pmusic)
grid on
toc
3、Root-Music算法
%root-music
clear all
clc
tic
%参数设定
M=10;
doa=[-10 45 60]/180pi;
P=length(doa);
f=1000;c=1500;
lambda=f/c;d=lambda/2;
snr=5;
N=128;
%阵列流型A
for i=1:P
A(:,i)=exp(-j2pid[0:M-1]’/lambdasin(doa(i)));
end
%信源模型建立
for k=1:P
S(k,:)=sqrt(10.^(snr/10))(randn(1,N)+jrandn(1,N));
end
%接收信号模型建立
X=AS+1/sqrt(2)(randn(M,N)+jrandn(M,N));
%协方差矩阵特征值分解得到噪声子空间
R=XX’/N;
[V,D]=eig®;
[Y,I]=sort(diag(D));
Un=V(:,I(1:M-P));
%ROOT-MUSIC算法部分
z =sym(‘z’,‘unreal’); %定义一个符号变量 z
for i=1:M
pz(i)=z^(i-1); %110 [1,z,z2,z3,z4,z5,z6,z7,z8,z9]
pzt(i)=z^(10-i); %1*10 [z9,z8,z7,z6,z5,z4,z3,z2,z,1]
end;
p=sym2poly(pztUnUn’pz.’);%多项式系数——>系数向量
zz=roots§;%求上述多项式的根，p是2(n-1)次多项式，根成对出现，相对于单位圆为镜像对，得到的根以列向量形式表示
z1=zz(abs(zz)<1);%挑出模值小于1的
z2=abs(abs(z1)-1);%计算各根离单位圆的距离
[Y,I]=sort(z2);%按升序重新排列上述各根
z3=z1(I(1:P));%取出最接近单位圆的P个根
for ii=1:P
sourcedoa(ii)=-asin(angle(z3(ii))/pi)/pi*180;%计算P个根对应的DOA
end
sourcedoa=sort(sourcedoa)%显示出按升序排列的所得DOA
stem(sourcedoa,ones(1,P));
axis([-90 90 0 2]);
grid on;
title(‘root-music算法’);
toc
4、前向平滑MUSIC算法
%FORWARD_SMOOTHNESS_MUSIC ALOGRITHM VS CLASSIC MUSIC
%DOA ESTIMATION BY FORWARD_SMOOTHNESS_MUSIC
clear all;
close all;
clc;
M=8;%阵元数
N=1024; %信号长度
snapshot_number=N;%快拍数
doa=[45,-60];%信号波达方向
P=length(doa);%信元数
f0=1000;%信号频率
fs=200;Ts=1/fs;
c=1500;%水中声速1500m/s
l=c/f0;%信号波长
d=0.5*l;%阵元间距
m=6;%每个子阵阵元数
L=M-m+1;%相互交错的子阵数
snr=0;%信噪比
A=[exp(-j*(0:M-1)d2pisin(doa(1)pi/180)/l);exp(-j(0:M-1)d2pisin(doa(2)*pi/180)/l)].’;%阵列流型
t=(0:N-1)Ts;
s1=sqrt(10.^(snr/10))exp(j(2pif0t+pi/1));
s2=sqrt(10.^(snr/10))exp(j(2pif0t+pi/2));
s=[s1;s2];%仿真信号,显然这是一个相干信源
%x=awgn(s,snr);
x=As+(1/sqrt(2))(randn(M,N)+jrandn(M,N));%加了高斯白噪声后的阵列接收信号
R=xx’/N;
Rf=zeros(m,m);
for i=1:L
y=x(i:i+m-1,:);
Rf=Rf+yy’/N;
end
Rf=Rf/L;
[U,S,V]=svd(Rf);
Un=U(:,P+1:m);
[V,D]=eig®;
UN=V(:,1:M-P);
theta=-90:0.1:90;%线阵的搜索范围为-90~90度
for ii=1:length(theta)
a=exp(-j*(0:m-1)‘2pidsin(pi*theta(ii)/180)/l);
Psmoothmusic(ii)=1./abs(a’UnUn’*a);
end
for kk=1:length(theta)
a=exp(-j*(0:M-1)‘2pidsin(pi*theta(kk)/180)/l);
Pmusic(kk)=1./abs(a’UNUN’*a);
end
plot(theta,10log(Psmoothmusic),‘R-’,theta,10log(Pmusic),‘B-’);
%axis([-90 90 -90 90]);
xlabel(‘入射角/度’);
ylabel(‘谱峰dB’);
legend(‘FORWARD-SMOOTHNESS-MUSIC’,‘CLASSIC MUSIC’);
title(‘前向空间平滑MUSIC估计与经典MUSIC的对比’);
grid on;
5、后向平滑MUSIC算法
�CKWARD_SMOOTHNESS_MUSIC ALOGRITHM VS CLASSIC MUSIC
%DOA ESTIMATION BY FORWARD_SMOOTHNESS_MUSIC
clear all;
close all;
clc;
M=8;%阵元数
N=1024; %信号长度
snapshot_number=N;%快拍数
doa=[20 40];%信号波达方向
P=length(doa);%信元数
f0=1000;%信号频率
fs=200;Ts=1/fs;
c=1500;%水中声速1500m/s
l=c/f0;%信号波长
d=0.5*l;%阵元间距
m=6;%每个子阵阵元数
L=M-m+1;%相互交错的子阵数
snr=0;%信噪比
A=[exp(-j*(0:M-1)d2pisin(doa(1)pi/180)/l);exp(-j(0:M-1)d2pisin(doa(2)*pi/180)/l)].’;%阵列流型
t=(0:N-1)Ts;
s1=sqrt(10.^(snr/10))exp(j(2pif0t+pi/1));
s2=sqrt(10.^(snr/10))exp(j(2pif0t+pi/2));
s=[s1;s2];%仿真信号,显然这是一个相干信源
%x=awgn(s,snr);
x=As+(1/sqrt(2))(randn(M,N)+jrandn(M,N));%加了高斯白噪声后的阵列接收信号
R=xx’/N;
Rf=zeros(m,m);
for i=M?M-L+1
y=x(i-m+1:i,:);
Rf=Rf+yy’/N;
end
Rf=Rf/L;
[U,S,V]=svd(Rf);
Un=U(:,P+1:m);
[V,D]=eig®;
UN=V(:,1:M-P);
theta=-90:0.1:90;%线阵的搜索范围为-90~90度
for ii=1:length(theta)
a=exp(-j*(0:m-1)‘2pidsin(pi*theta(ii)/180)/l);
Psmoothmusic(ii)=1./abs(a’UnUn’*a);
end
for kk=1:length(theta)
a=exp(-j*(0:M-1)‘2pidsin(pi*theta(kk)/180)/l);
Pmusic(kk)=1./abs(a’UNUN’*a);
end
plot(theta,10log(Psmoothmusic/max(Psmoothmusic)),‘R-’,theta,10log(Pmusic/max(Pmusic)),‘B-’);
axis([-90 90 -100 20]);
xlabel(‘入射角/度’);
ylabel(‘谱峰dB’);
legend(‘BACKWARD-SMOOTHNESS-MUSIC’,‘CLASSIC MUSIC’);
title(‘后向空间平滑MUSIC估计与经典MUSIC的对比’);
grid on;
6、双向平滑MUSIC算法
%BIOWARD_SMOOTHNESS_MUSIC ALOGRITHM VS CLASSIC MUSIC
%DOA ESTIMATION BY FORWARD_SMOOTHNESS_MUSIC
clear all;
close all;
clc;
M=16;%阵元数
N=1024; %信号长度
snapshot_number=N;%快拍数
doa=[0 5];%信号波达方向
P=length(doa);%信元数
f0=1000;%信号频率
fs=200;Ts=1/fs;
c=1500;%水中声速1500m/s
l=c/f0;%信号波长
d=0.5*l;%阵元间距
m=6;%每个子阵阵元数
L=M-m+1;%相互交错的子阵数
snr=0;%信噪比
A=[exp(-j*(0:M-1)d2pisin(doa(1)pi/180)/l);exp(-j(0:M-1)d2pisin(doa(2)*pi/180)/l)].’;%阵列流型
t=(0:N-1)Ts;
s1=sqrt(10.^(snr/10))exp(j(2pif0t+pi/1));
s2=sqrt(10.^(snr/10))exp(j(2pif0t+pi/2));
s=[s1;s2];%仿真信号,显然这是一个相干信源
%x=awgn(s,snr);
x=As+(1/sqrt(2))(randn(M,N)+jrandn(M,N));%加了高斯白噪声后的阵列接收信号
R=xx’/N;
Rb=zeros(m,m);
for i=M?M-L+1
y=x(i-m+1:i,:);
Rb=Rb+yy’/N;
end
Rb=Rb/L;
Rf=zeros(m,m);
for i=1:L
y=x(i:i+m-1,:);
Rf=Rf+y*y’/N;
end
Rf=Rf/L;
R2=(Rf+Rb)/2;
[U,S,V]=svd(R2);
Un=U(:,P+1:m);
[V,D]=eig®;
UN=V(:,1:M-P);
theta=-90:0.1:90;%线阵的搜索范围为-90~90度
for ii=1:length(theta)
a=exp(-j*(0:m-1)‘2pidsin(pi*theta(ii)/180)/l);
Psmoothmusic(ii)=1./abs(a’UnUn’*a);
end
for kk=1:length(theta)
a=exp(-j*(0:M-1)‘2pidsin(pi*theta(kk)/180)/l);
Pmusic(kk)=1./abs(a’UNUN’*a);
end
plot(theta,10log(Psmoothmusic/max(Psmoothmusic)),‘R-’,theta,10log(Pmusic/max(Pmusic)),‘B-’);
xlabel(‘入射角/度’);
ylabel(‘谱峰dB’);
legend(‘双向空间平滑MUSIC’,‘经典MUSIC’);
title(‘双向向空间平滑MUSIC估计与经典MUSIC的对比’);
grid on;
7、奇异值算法
程序1
%%%%%%针对数据处理的矢量奇异值法（DSVD）%%%%%%%%%
%%%该方法只适合线阵；必须满足 子阵阵元数m>信源个数P, 子阵数p>信源个数P；且矢量奇异值法解相干有孔径的损失
clear all
clc
%%%参数设置
M=8;
N=128;
doa=[20 45 60]/180*pi;
P=length(doa);
f=100;c=1500;
lambda=c/f;
d=lambda/2;
snr=20;
p=4;%所有阵元划分3个子阵
m=M-p+1;%每个子阵中阵元数为8-3+1=6个
%%%阵列流型A
for i=1:P
A(:,i)=exp(-j2pi*(0:M-1)'dsin(doa(i))/lambda);
end
%%%信源信号矩阵建立
Ts=1;
t=(0:N-1)Ts;
for i=1:P
s(i,:)=sqrt(10.^(snr/10))exp(j(2pift+pi/i));
end
%%%接收信号矩阵
X=As+(1/sqrt(2))(randn(M,N)+j*randn(M,N));%加了高斯白噪声后的阵列接收信号
x0=X(1,:);%取第一个阵元的接收数据位参考阵元数据矢量
X1=X*x0’/N;%得到一个包含所有信号信息的数据矢量
%%利用上述数据矢量进行矩阵重构
for k=1:m
for kk=1:p
Y(k,kk)=X1(k+kk-1);
end
end
%%对重构的矩阵进行奇异值分解
[U,S,V]=svd(Y);
Un=U(:,P+1:m);
%%%谱峰搜索
theta=-90:0.1:90;
for iii=1:length(theta)
a_theta=exp(-j2pi*(0:m-1)'dsin(theta(iii)/180*pi)/lambda);
pmusic(iii)=(abs(1/(a_theta’UnUn’a_theta)));
end
pmusic=10log10(pmusic/max(pmusic));
figure(1)
plot(theta,pmusic)
xlabel(‘theta/degree’)
ylabel(‘pmusic/dB’)
title(‘DSVD MUSIC’)
grid on
%%以求根算法给出谱峰点
syms z pz pzt;
pz=z.1.’;
pzt=z.2;
a=sym2poly(pztUnUn’pz);
zz=roots(a)
zz1=zz(abs(zz)<1)
zz2=abs(abs(zz1)-1)
[Y,I]=sort(zz2)
final=zz1(I(1:P));
sorcedoa=asin(-lambdaangle(final)/2/pi/d)180/pi;
sorcedoa=sort(sorcedoa)
程序2
%%%%%%特征矢量奇异值法（ESVD）%%%%%%%%%
%%%该方法只适合线阵；必须满足 子阵阵元数m>信源个数P, 子阵数p>信源个数P；且矢量奇异值法解相干有孔径的损失
clear all
clc
%%%参数设置
M=8;
N=128;
doa=[20 45 60]/180pi;
P=length(doa);
f=1000;c=1500;
lambda=c/f;
d=lambda/2;
snr=20;
p=4;
m=M-p+1;
%%%阵列流型A
for i=1:P
A(:,i)=exp(-j2pi*(0:M-1)'dsin(doa(i))/lambda);
end
%%%相干信源模型建立
Ts=1;
t=(0:N-1)Ts;
for i=1:P
S(i,:)=sqrt(10.^(snr/10))exp(j(2pift+pi/i));
end
%%%接收信号模型建立
X=AS+(1/sqrt(2))(randn(M,N)+jrandn(M,N));
%%%找出最大特征值对应的特征矢量e1
R=XX’/N;
[U,D]=eig®;
e1=U(:,M);
%%%重构矩阵Y
for i=1:m
for k=1:p
Y(i,k)=e1(i+k-1);
end
end
%%%对重构的矩阵进行奇异值分解,找出噪声子空间Un
[U,D,V]=svd(Y);
Un=U(:,P+1:m);
%%谱峰搜素
theta=-90:0.1:90;
for i=1:length(theta)
a=exp(-j2pi*(0:m-1)'dsin(theta(i)/180pi)/lambda);
pmusic(i)=abs(1/(a’UnUn’a));
end
pmusic=10log10(pmusic/max(pmusic));
figure(1)
plot(theta,pmusic)
xlabel(‘theta/degree’)
ylabel(‘pmusic/dB’)
title(‘ESVD SPECTRUM’)
grid on
%%%以求根算法给出谱峰点
syms z pz pzt
pz=z.3;
pzt=z.4.’;
aa=sym2poly(pzUn*Un’pzt);
zz=roots(aa);
zz1=zz(abs(zz)<1);
zz2=abs(abs(zz1)-1);
[Y,I]=sort(zz2);
final=zz1(I(1:P));
sourcedoa=asin(-lambdaangle(final)/2/pi/d)*180/pi;
sourcedoa=sort(sourcedoa)
程序3
%PSVD, 注意成功概率和信噪比有关
%DOA ESTIMATION BY VECTOR_SINGULAR_VALUE_MUSIC
clear all;
close all;
clc;
%%%参数设定
M=8;
N=128;
doa=[20 45 55]/180pi;
P=length(doa);
f=1000;c=1500;
lambda=c/f;
d=lambda/2;
snr=10;
p=4;
m=M-p+1;
doa_probable=40;%信号可能到达角度
%%%阵列流型A
for i=1:P
A(:,i)=exp(-j2pi(0:M-1)'dsin(doa(i))/lambda);
end
%%%相干信源模型建立
Ts=1;
t=(0:N-1)Ts;
for i=1:P
S(i,:)=sqrt(10.^(snr/10))exp(j(2pift+pi/i));
end
%%%接收信号模型建立
X=AS+(1/sqrt(2))(randn(M,N)+j*randn(M,N));
a_doa_probable=exp(-j*(0:M-1)‘d2pisin(doa_probable*pi/180)/lambda);%假设角的方向矢量
b=zeros(1,N);
y=zeros(8,1);
for pp=1:N
b(pp)=a_doa_probable’X(:,pp)/8;%参考阵元矢量
y=b(pp)X(:,pp)+y;
end
y1=y/N;
%进行矩阵重构
for i=1:m
for k=1:p
Y(i,k)=y1(i+k-1);
end
end
%[U,S,V]=svd(Y);
%Un=V(:,1:m-source_number);
%Gn=UnUn’;
[U,S,V]=svd(Y);
Un=U(:,P+1:m);
Gn=UnUn’;
searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
for i=1:length(searching_doa)
a_theta=exp(-j*(0:m-1)'2pidsin(pi*searching_doa(i)/180)/lambda);
Pmusic(i)=1./abs((a_theta)'Gna_theta);
end
plot(searching_doa,10log(Pmusic));
%axis([-90 90 -90 90]);
xlabel(‘入射角/度’);
ylabel(‘谱峰、dB’);
legend(‘VECTOR-SINGULAR-VALUE-MUSIC SPECTRUM’);
title(‘DSVDMUSIC估计’);
grid on;
8、线性预测算法
%%%%%%@@@@@@@@@
%%%%%%空间谱估计中的线性预测算法，包括前向预测算法，后向预测算法，双向预测算法，多阶线性预测算法◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎◎
clc;
clear;
close all;
tic
M=16;%阵元数目
N=3;%信源数目
B=8;%形成的波束个数
snap=1000;%快拍数目
C=3e8;
f0=10e6;
lamda=C/f0;
d=0.5lamda;
% k=d/lamda;
theta0=20;
theta1=30;
theta2=-10;
fs=1000;
ts=1/fs;
t=(0:snap-1)*ts;
a=[0:M-1]’;%阵列矢量
u0=5;
u1=10;
u2=20;
s0=exp(j2pi0.5u0t.^2);
s1=exp(j2pi0.5u0t.^2);
s2=exp(j2pi0.5u0*t.^2);
%阵列流行矢量
a_theta0=exp(j2pid/lamdaasin(theta0/180pi));
a_theta1=exp(j2pid/lamdaasin(theta1/180pi));
a_theta2=exp(j2pid/lamdaasin(theta2/180pi));
A=[a_theta0 a_theta1 a_theta2];
S=[s0;s1;s2];
X0=AS;
SNR=15;
randn(‘state’,0)
real_noise=randn(size(X0));
randn(‘state’,3);
imag_noise=randn(size(X0));
noise0=(real_noise+jimag_noise)/2^0.5;
noise=10^(-SNR/20)*noise0;
X=X0+noise;%%%%%%%%%%------------=========完整的基带信号
Rx=XX’/length(t);
%%%@@@@@@@@@@@@@@@@下面为一阶线性预测算法
%%%%%%%%%%前向预测算法
Xf=flipud(X(1:M-1,:));Xft=transpose(Xf);Xfh=Xf’;
Rf=XfXfh/snap;rf=XfX(M,:)’/snap;
Wflp=conj((inv(Rf)rf));
%%%%%%%%%%后向预测算法
Xb=X(2:M,:);Xbt=transpose(Xb);Xbh=Xb’;
Rb=XbXbh/snap;rb=XbX(1,:)’/snap;
Wblp=inv(Rb)rb;
%%%%%%%%%%双向预测算法
Xfbt=[Xft;conj(Xbt)];Xfb=transpose(Xfbt);Xfbh=Xfb’;
Rfb=XfbXfbh/snap;rfb=Xfb*[X(M,:)’;transpose(X(1,:))]/snap;
Wfblp=conj(inv(Rfb)*rfb);
%%%@@@@@@@@@@@@@@@@下面为多阶线性预测算法
p=10;P=M-p;%线性预测的阶数
J=fliplr(eye§);%P阶交换矩阵
RF0=zeros(P,P);%前向预测的协方差阵
Rf0=zeros(P,1);
af0=zeros(P,1);
RB=zeros(P,P);%后向预测的协方差阵
Rb=zeros(P,1);
theta=-90:0.01:90;
for ii=1:length(theta)
a_theta=exp(j2pid/lamdaasin(theta(ii)/180pi));
Pflp(ii)=10log10(1/(abs(a_theta’[1;-Wflp]))^2);
Pblp(ii)=10log10(1/(abs(a_theta’[1;-Wblp]))^2);
Pfblp(ii)=10log10(1/(abs(a_theta’[1;-Wfblp]))^2);
for jj=1:p
RF=JRx(jj:P+jj-1,jj:P+jj-1)J;
Rf=JRx(jj:P+jj-1,P+jj);
WFLP=conj((inv(RF)Rf));
PFLP(jj,ii)=10log10(1/(abs(a_theta(1:P+1)’[1;-WFLP]))^2);
RB=Rx(jj+1:P+jj,jj+1:P+jj);
Rb=Rx(jj+1:P+jj,jj);
WBLP=(inv(RB)Rb);
PBLP(jj,ii)=10log10(1/(abs(a_theta(1:P+1)’[1;-WBLP]))^2);
RFB=RF+conj(RB);
Rfb=Rf+conj(Rb);
WFBLP=conj((inv(RFB)Rfb));
PFBLP(jj,ii)=10log10(1/(abs(a_theta(1:P+1)’[1;-WFBLP]))^2);
end
end
figure(1);
plot(theta,Pflp);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘一阶前向预测’);
figure(2);
plot(theta,Pblp);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘一阶后向预测’);
figure(3);
plot(theta,Pfblp);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘一阶双向预测’);
figure(4);
plot(theta,sum(PFLP,1)/p);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘多阶前向预测’);
figure(5);
plot(theta,sum(PBLP,1)/p);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘多阶后向预测’);
figure(6);
plot(theta,sum(PFBLP,1)/p);grid on;xlabel(‘角度’);ylabel(‘峰值’);title(‘多阶双向预测’);
toc
9、旋转不变子空间算法
% ESPRIT算法
clc
clear all
close all
M=8; %阵元数
SNR=0; %信噪比
theta=[-20 20]pi/180; %信号入射角
N=256; %快拍数
f=3000; %信号频率
fs=f5; %采样频率
c=1500; %声速
lambda=c/f; %波长
d=lambda/2; %阵元间距
%%%%%%%%%%%%%%%% produce signal%%%%%%%%%%%%%%%%%
for n=1:N
S(1,n)=sqrt(210^(SNR/10))exp(j2pif/fsn+2pirand(1,1)); %独立的信号源1
S(2,n)=sqrt(210^(SNR/10))exp(j2pi2f/fsn+2pirand(1,1)); %独立的信号源2
end
SS=[S(1,:);S(2,:)];
n1=randn(M,N);
n2=randn(M,N);
Q=[exp(j2pi/lambdadsin(theta(1))),0;0,exp(j2pi/lambdadsin(theta(2)))];
A=[exp(-j2pi(0:M-1)’/lambdadsin(theta))];
X=ASS+n1; %基阵接收的信号
Y=AQ*SS+n1; %基阵接收的信号
%%%%%%%%%%%%%%% ESPRIT算法 %%%%%%%%%%%%%%%%%%%%%%
Rxx=1/NXX’; %求自协方差矩阵
Rxy=1/NXY’; %求互协方差矩阵
[eigvalxx]=eig(Rxx); %进行特征值分解
[eigvalxy]=eig(Rxy); %进行特征值分解
cxx=Rxx-min(eigvalxx)*eye(M);
z=tril(ones(M),-1);
z=triu(z,-1);
cxy=Rxy-min(eigvalxy)*z;
[eigve,eigva]=eig(cxx,cxy);
eigva=diag(eigva);
ss=abs(abs(eigva)-1);
[s,b]=sort(ss);
%计算DOA
for ii=1:2
%angle_out0(ii)=acos((1500/f0)angle(eigva0(b0(ii)))/(2pi*deta));
angle_out(ii)=asin(angle(eigva(b(ii)))/pi)*180/pi;
end
angle_out=sort(angle_out)
eig_value=eigva;
10、最大熵算法、最小模算法、最小方差法
%1、最大熵算法； 2、最小方差算法； 3、双向预测算法； 4、最小模算法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 空间谱估计中估计的线性预测算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
close all
clc
%%%%%%%%%%%%%%%%%%%%%%% generate signal %%%%%%%%%%%%%%%%%
f0=5000; %入射信号频率
fs=50f0; %采样频率
ts=1/fs; %采样间隔
% n=[1:256];
% t=nts;
M=8; %阵元数
Q=2; %信号源数
L=200; %快拍数
SNR=20; %%%%% 信噪比 %%%%%%
z1=5pi/180;
z2=10pi/180; %z1,z2为信号入射方向
c=1500;d=0.15; %c为水下声速，d为阵元间距
t=(1:L);
K=sqrt(2*10^(SNR/10));
s1=Kexp(j(2*(t-1)pif0/fs+2pirand(1,L)));
s2=Kexp(j(2*(t-1)pif0/fs+2pirand(1,L)));
s3=Kexp(j(2*(t-1)pi2f0/fs+pi/2));
s4=Kexp(j*(2*(t-1)pi2*f0/fs+pi));
s=[s1;s2]; %独立信号源矩阵
s0=[s3;s4]; %相干信号源矩阵
%%%%%%%%%%%%%%%%%%%%% 生成导向矢量矩阵 %%%%%%%%%%%%%%%%%%%
tao1=sin(z1);
tao2=sin(z2);
m=1:M;
A1=exp(-jpi(m)'tao1);
A2=exp(-jpi*(m)‘tao2);
A=[A1 A2];
%%%%%%%%%%%%%%%%%%%%% add wgn %%%%%%%%%%%%%%%%%%%%%
N=(randn(M,L)+jrandn(M,L))/sqrt(2); % 高斯白噪声
X=As+N; %阵列接收的独立信号
X0=As0+N; %阵列接收的相干信号
figure(1)
subplot(2,1,1);plot(t,real(s1));
subplot(2,1,2);plot(t,real(s2));
figure(2)
subplot(1,1,1);plot(t,real(s3),t,real(s4));
figure(3)
subplot(2,1,1);plot(t,abs(As));
subplot(2,1,2);plot(t,real(X));
figure(4)
subplot(2,1,1);plot(t,real(As0));
subplot(2,1,2);plot(t,real(X0));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 线性预测算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R=XX’/L; %采样协方差矩阵
R0=X0X0’/L;
[p,r]=eig®; %特征分解
[p0,r0]=eig(R0); %特征分解
Us(:,[1,2])=p(:,[1,2]); %对R特征分解并取出最大特征值及其对应的特征向量为信号子空间Us
Un(:,[1,2,3,4,5,6])=p(:,[3,4,5,6,7,8]); %对R特征分解并取出小特征值及其对应的特征向量为噪声子空间Un
Us0(:,[1,2])=p0(:,[1,2]); %对R特征分解并取出最大特征值及其对应的特征向量为信号子空间Us
Un0(:,[1,2,3,4,5,6])=p0(:,[3,4,5,6,7,8]); %对R特征分解并取出小特征值及其对应的特征向量为噪声子空间Un
%%%%%%%%%%% 最大商算法 %%%%%%%%%%%%
u0=[1 0 0 0 0 0 0 0]’;
i=1;
for theta=-90:0.1:90;
arfa=sin(thetapi/180)d/c;
for m=1:M;
tao0(1,m)=(m-1)arfa;
end
a=exp(-j2pif0*tao0).’; %方向向量
Pmem(i)=1./(abs(a’*inv®u0)).^2;
Pmem0(i)=1./(abs(a’inv(R0)u0)).^2;
i=i+1;
end
figure(5)
theta0=-90:0.1:90;
plot(theta0,10log10(abs(Pmem)),’-’,theta0,10log10(abs(Pmem0)),‘±’);
xlabel(‘入射角度’),ylabel(‘空间方位谱(dB)’),title(‘最大商算法’),grid on
%%%%%%%%%%%%% 最小方差算法 %%%%%%%%%%%%%%
i=1;
for theta=-90:0.1:90;
arfa=sin(thetapi/180)d/c;
for m=1:M;
tao0(1,m)=(m-1)arfa;
end
a=exp(-j2pif0*tao0).’;
Pmvm(i)=1./(a’*inv®a);
Pmvm0(i)=1./(a’inv(R0)a);
i=i+1;
end
figure(6)
theta0=-90:0.1:90;
plot(theta0,10log10(abs(Pmvm)),’-’,theta0,10log10(abs(Pmvm0)),‘±’);
xlabel(‘入射角度’),ylabel(‘空间方位谱(dB)’),title(‘最小方差算法’),grid on
%%%%%%%%%%%%%%% 双向预测算法 %%%%%%%%%%%%%%%%
%后向预测
XT=X.’;
XT0=X0.’;
for m0=2:M;
Xb(:,m0-1)=XT(:,m0);
Xb0(:,m0-1)=(XT0(:,m0));
end
%前向预测
for m1=1:M-1;
Xf(:,m1)=conj(XT(:,M-m1));
Xf0(:,m1)=conj(XT0(:,M-m1));
end
X1=XT(:,1);
X10=XT0(:,1);
XM=conj(XT(:,M));
XM0=conj(XT0(:,M));
XfbT=[Xb;Xf];
XfbT0=[Xb0;Xf0];
FF=[X1;XM];
FF0=[X10;XM0];
Wfblp=((pinv(XfbT))*FF);
Wfblp0=((pinv(XfbT0))*FF0); % 双向预测权矢量
i=1;
for theta=-90:0.1:90;
arfa=-sin(thetapi/180)d/c;
for m=1:M;
tao0(1,m)=(m-1)arfa;
end
a=exp(-j2pif0tao0).’;
Pfblp(i)=1./(abs(a’[1;-Wfblp]))^2;
Pfblp0(i)=1./(abs(a’[1;-Wfblp0]))^2;
i=i+1;
end
figure(7)
theta0=-90:0.1:90;
plot(theta0,10log10(abs(Pfblp)),’-’,theta0,10log10(abs(Pfblp0)),‘±’);
xlabel(‘入射角度’),ylabel(‘空间方位谱(dB)’),title(‘双向预测算法’),grid on
%%%%%%%%%%%%%%%%%% 最小模算法 %%%%%%%%%%%%%%%%
gT=Us(1,:); %提取出信号子空间的第一行矢量
gT0=Us0(1,:); %提取出信号子空间的第一行矢量
g=gT.’;
g0=gT0.’;
for m0=2:M; %提取出信号子空间的第二行矢量
Es(m0-1,:)=Us(m0,:);
Es0(m0-1,:)=Us0(m0,:);
end
Wmnm=[1;-Esconj(g)/(1-g.'conj(g))];
Wmnm0=[1;-Es0conj(g0)/(1-g0.'g0)];
i=1;
for theta=-90:0.1:90;
arfa=sin(thetapi/180)d/c;
for m=1:M;
tao0(1,m)=(m-1)arfa;
end
a=exp(-j2pif0tao0).’;
Pmnm(i)=1./(abs(a’Wmnm))^2;
Pmnm0(i)=1./(abs(a’Wmnm0))^2;
i=i+1;
end
figure(8)
theta0=-90:0.1:90;
plot(theta0,10log10(abs(Pmnm)),’-’,theta0,10log10(abs(Pmnm0)),‘±’);
xlabel(‘入射角度’),ylabel(‘空间方位谱(dB)’),title(‘最小模算法’),grid on
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 四种算法的空间谱比较 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(9) %% 接收独立信号源
theta0=-90:0.1:90;
plot(theta0,10log10(abs(Pmem)),’-’,theta0,10log10(abs(Pmvm)),’-b’,theta0,10log10(abs(Pfblp)),’–’,theta0,10*log10(abs(Pmnm)),‘±’)
legend(‘MEM’,‘MVM’,‘FBLP’,‘MNM’),xlabel(‘入射角度’),ylabel(‘空间方位谱(dB)’),title(‘接收独立信号源’),grid on
grid on
0:m-1 ↩︎
m-1?0 ↩︎
m-1?0 ↩︎
0:m-1 ↩︎
最后
以上就是义气春天为你收集整理的空间谱估计的经典算法（MATLAB程序）的全部内容，希望文章能够帮你解决空间谱估计的经典算法（MATLAB程序）所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复