clear;
close all;
clc;
Ns=512;%信号长度??????????
snapshots=512;%?快拍数?????
snrdB=0;snr=10^(snrdB/20); %??信噪比???
c=34000;%波速
f0=1500;%中心频率,频率范围100-3100??
lambda=c/f0;%波长?????
d=0.8*lambda;%?阵元间距??
M=5;%?阵元数??
angle=[60 80];%?信号方向?
P=length(angle);%?信号数目
frequency_band=3000;%频率带宽?????
fs=3000;%?采样率
Ts=1/fs;
t=Ts:Ts:Ns*Ts;
f1=100;%??频率? ???
ffd=3000/20;
for ffk=1:21
fff(ffk)=f1+ffd*(ffk-1);
end
% % for ffk=12:21
% % fff(ffk)=f0+ffd*(ffk-10-1);
% % end ?
% % ff=c/(2*d);
% % for k=1:21
% % fmin(k)=abs(fff(k)-ff);??
% % end
% % Min=min(fmin);
% % for k=1:21
% % if fmin(k)==Min;
% % fy=fff(k);?
% % end
% % end
% % fy=repmat(fy,[1 P]);
fy=repmat(f0,[1 P]);
for k=1:M
for sig=1:P
A0(k,sig)=exp(j*2*pi*fy(sig)*(k-1)*d*cos(angle(sig)*pi/180)/c);
end;
end;
s=zeros(P,snapshots);
Y=zeros(M,1024);
Am=[3,5];%????
phi=[20,30];%????
for jj=1:21
f=100+(jj-1)*150;
for ii=1:P
s(ii,:)=s(ii,:)+Am(ii)*cos(2*pi*f*t+phi(ii)*pi/180);
%s(ii,:)=s(ii,:)+cos(2*pi*f*t+pi/180);
end
end
clear jj f;
for jj=1:21
f=100+(jj-1)*150;
for i=1:P
a0(:,i)=exp(j*[0:M-1]'*d*2*pi*cos(angle(i)*pi/180)*f/c);
end
n=exp(j*2*pi*randn(M,snapshots))/snr;%????
X=a0*s+n;
XX11=conj(X');
XX1=fft(XX11,1024);
XX2=conj(XX1');
for sig=1:P
%A(:,sig)=exp(j*2*pi*fff(jj)*(0:M-1)'*d*cos(angle(sig)*pi/180)/c);
A(:,sig)=exp(j*2*pi*f*[0:M-1]'*d*cos(angle(sig)*pi/180)/c);
end
AA=A0*A';
[U S V]=svd(AA);
T=U*V';
Y=Y+T*XX2;
end
Y=Y/jj;
clear i;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%MUSIC算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ryy=(Y*Y')/snapshots;
[EigenVectors,EigenValues]=eig(Ryy);
Lemda=diag(EigenValues);
[SortedLemda,Index]=sort(Lemda);
Index=flipud(Index);
NoiseSubspace(1:M,1:M-P)=EigenVectors(1:M,Index(P+1:M));
for i=1:1:180
a=exp(j*2*pi*f0*d*[0:M-1]'*cos(i*pi/180)/c);
MUSIC_Spec(i)=a'*a/(a'*NoiseSubspace*(NoiseSubspace)'*a);
end
M_MUSIC_Spec=max(MUSIC_Spec);
abslog_MUSIC_Spec=20*log10(abs(MUSIC_Spec));
ang=[];
index=find(abslog_MUSIC_Spec<60);
abslog_MUSIC_Spec(index)=0;
for i=1:1:180
if(i==1)
left=180;
right=i+1;
elseif(i==180)
left=i-1;
right=1;
else
left=i-1;
right=i+1;
end;
if(abslog_MUSIC_Spec(i)>0)
if(abslog_MUSIC_Spec(i)>abslog_MUSIC_Spec(left)&&abslog_MUSIC_Spec(i)>abslog_MUSIC_Spec(right))
ang=[ang,i];
fprintf('\tpeak_angle= %d\n',i);
end;
end
end
figure
MUSIC=abs(MUSIC_Spec)/max(abs(MUSIC_Spec));
Delta=[1:180];
plot(Delta,20*log10(abs(MUSIC_Spec)),'r'); %NonNormative Spatial Spectrum
hold on;
stem(ang,abslog_MUSIC_Spec(ang))
xlabel('angle/\circ'),ylabel('Power Spectrum(dB)');