Matlab实现MUSIC算法、加权MUSIC算法和ROOT-MUSIC算法对比 上传版本.zip

clear; close all; clc; bbb=zeros(1,10); sensor_number=[3 4 6 8 10 12 14 16 18 20];%阵元数 aaa=zeros(1,10); source_number=1;%信源数? N_x=64; %信号长度 snapshot_number=N_x;%?快拍数 w=2000;%频率 l=340/w;%波长 d=0.5*l;%?阵元间距 snr=10;%信噪比 source_doa=60;%DOA角度 %WMUSIC，加权MUSIC算法 for kk=1:10 for k=1:10 A=[exp(-j*(0:sensor_number(kk)-1)*d*2*pi*sin(source_doa*pi/180)/l)].'; s=10.^((snr/2)/10)*exp(j*w*[0:N_x-1]);%信号 %x=awgn(s,snr); x=A*s+(1/sqrt(2))*(randn(sensor_number(kk),N_x)+j*randn(sensor_number(kk),N_x));%阵列接收信号加噪声 R=x*x'/N_x;%协方差矩阵 %[V,D]=eig(R); %Un=V(:,1:sensor_number(kk)-source_number); %Gn=Un*Un'; [V,D]=eig(R); %特征分解 D=diag(D); Un=V(:,1:sensor_number(kk)-source_number); Gn=Un*Un'; e=[1 zeros(1,sensor_number(kk)-1)].'; W=e*e.'; %加权值 searching_doa=-90:0.1:90;%搜索角度 for i=1:length(searching_doa) a_theta=exp(-j*(0:sensor_number(kk)-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l); Pmusic(i)=1./abs((a_theta)'*Gn*W*Gn*a_theta); end aa=diff(Pmusic); aa=sign(aa); aa=diff(aa); bb=find(aa==-2)+1; [t1 t2]=max(Pmusic(bb)); estimated_source_doa=searching_doa(bb(t2)); aaa(:,k)=estimated_source_doa; end E_source_doa=sum(aaa(1,:))/10; RMSE_source_doa=sqrt(sum((aaa(1,:)-source_doa).^2)/10);%均方误差 bbb(:,kk)=RMSE_source_doa; end hold on plot(sensor_number,bbb(1,:),'rs-'); %MUSIC算法???? for kk=1:10 for k=1:10 A=[exp(-j*(0:sensor_number(kk)-1)*d*2*pi*sin(source_doa*pi/180)/l)].'; s=10.^((snr/2)/10)*exp(j*w*[0:N_x-1]);%信号 %x=awgn(s,snr); x=A*s+(1/sqrt(2))*(randn(sensor_number(kk),N_x)+j*randn(sensor_number(kk),N_x)); R=x*x'/snapshot_number; %[V,D]=eig(R); %Un=V(:,1:sensor_number(kk)-source_number); %Gn=Un*Un'; [V,D]=eig(R); D=diag(D); Un=V(:,1:sensor_number(kk)-source_number); Gn=Un*Un'; searching_doa=-90:0.1:90;%搜索角度 for i=1:length(searching_doa) a_theta=exp(-j*(0:sensor_number(kk)-1)'*2*pi*d*sin(pi*searching_doa(i)/180)/l); Pmusic(i)=1./abs((a_theta)'*Gn*a_theta); end aa=diff(Pmusic); aa=sign(aa); aa=diff(aa); bb=find(aa==-2)+1; [t1 t2]=max(Pmusic(bb)); estimated_source_doa=searching_doa(bb(t2)); aaa(:,k)=estimated_source_doa; end E_source_doa=sum(aaa(1,:))/10; % E_source_doa=mean(aaa); RMSE_source_doa=sqrt(sum((aaa(1,:)-source_doa).^2)/10); bbb(:,kk)=RMSE_source_doa; end plot(sensor_number,bbb(1,:),'k*-'); %ROOT_MUSIC for kk=1:10 for k=1:10 A=[exp(-j*(0:sensor_number(kk)-1)*d*2*pi*sin(source_doa*pi/180)/l)].'; s=10.^((snr/2)/10)*exp(j*w*[0:N_x-1]); %x=awgn(s,snr); x=A*s+(1/sqrt(2))*(randn(sensor_number(kk),N_x)+j*randn(sensor_number(kk),N_x)); R=(x*x')/snapshot_number; [V,D]=eig(R); Un=V(:,1:sensor_number(kk)-source_number); Gn=Un*Un'; a = zeros(2*sensor_number(kk)-1,1)'; for i=-(sensor_number(kk)-1):(sensor_number(kk)-1) a(i+sensor_number(kk)) = sum( diag(Gn,i) ); end a1=roots(a); a2=a1(abs(a1)<1); [lamda,I]=sort(abs(abs(a2)-1)); f=a2(I(1:source_number)); estimated_source_doa=[asin(angle(f(1))/pi)*180/pi]; aaa(:,k)=estimated_source_doa; end E_source_doa=sum(aaa(1,:))/10; RMSE_source_doa=sqrt(sum((aaa(1,:)-source_doa).^2)/10); bbb(:,kk)=RMSE_source_doa; end hold on plot(sensor_number,bbb(1,:),'bd-'); legend('WMUSIC','MUSIC','ROOT-MUSIC'); xlabel('阵元数'); ylabel('RMSE'); grid on;