问题：如何利用matlab修改已有程序，然后绘制图中的所给图像？
问题描述：已有程序比类似代码多了一个优化矩阵，结果是绘制信噪比与均方根差的关系图。类似代码就是提供一些参照代码以及求解均方根差的代码，类似代码和已有代码维度可能不一样，可以修改的。
结果类似图：
已有程序：
clear;
close all;
%%%%%%%%%%%参数设定%%%%%%%%%%%
R = pi/180; %角度->弧度
M = 3; % 信源数目
N = 8; % 阵元个数
theta = [-30 0 60]; % 待估计角度
snr = 10; % 信噪比
Sample =500; % 快拍数(抽样密度)
element_spacing = 0.5; % 阵元间距
d=0:element_spacing:(N-1)*element_spacing;
A=exp(-1i*2*pi*d.'*sin(theta*R)); %方向矢量
%%%%%%%%%%构建信号模型%%%%%%%%
S=randn(M,Sample); %信源信号，入射信号
X=A*S; %构造接收信号
X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中
% 计算协方差矩阵
Rxx=X1*X1'/Sample;
% 特征值分解
[Gn,D]=eig(Rxx); %特征值分解
EVA=diag(D)'; %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA); %将特征值排序 从小到大
Gn=Gn(:,I); % 对应特征矢量排序
%%%%%%%%%%% 优化矩阵U %%%%%%%%%
q = EVA(1);
SUM = 0;
for k = N-M+1:N
u = Gn(:,k);
SUM = SUM+(EVA(k)*u*u')/(q-EVA(k))^2;
end
U=q*SUM;
%%%%%%%遍历每个角度，计算空间谱%%%%%%
for A = 1:361
angle(A)=(A-181)/2;
phim=R*angle(A);
a=exp(-1i*2*pi*d*sin(phim)).';
G=Gn(:,1:N-M); % 取矩阵的第1到N-M列组成噪声子空间
Pmusic(A)=(a'*U*a)/(a'*G*G'*a);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic);
h=plot(angle,Pmusic);
set(h,'Linewidth',2);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;
绘制结果图像类似代码：
%MUSIC ALOGRITHM
%DOA ESTIMATION BY CLASSICAL_MUSIC
clear all;
%close all;
clc;
bbb=zeros(1,11);
for kk=1:11
snr=[-10 -8 -6 -4 -2 0 2 4 6 8 10];%信噪比(dB)
aaa=zeros(1,300);
for k=1:300
iwave=1;%信元数
kelm=8;%阵元数
N_x=1024; %信号长度
snapshot_number=N_x;%快拍数
w=pi/4;%信号频率
L=2*pi*3e8/w;%信号波长
d=0.5*L;%阵元间距
source_doa=50;%信号的入射角度
A=[exp(-j*(0:kelm-1)*d*2*pi*sin(source_doa*pi/180)/L)].';
s=10.^((snr(kk)/2)/10)*exp(j*w*[0:N_x-1]);%仿真信号
%x=awgn(s,snr);
x=A*s+(1/sqrt(2))*(randn(kelm,N_x)+j*randn(kelm,N_x));%加了高斯白噪声后的阵列接收信号
R=x*x'/N_x;
%[V,D]=eig(R);
%Un=V(:,1:sensor_number-source_number);
%Gn=Un*Un';
[V,D]=eig(R);
D=diag(D);
disp(D);
Un=V(:,1:kelm-iwave);
Gn=Un*Un';
searching_doa=-90:0.1:90;%线阵的搜索范围为-90~90度
for i=1:length(searching_doa)
a_theta=exp(-j*(0:kelm-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
disp(aaa);
%求解均方根误差和标准偏差
E_source_doa=sum(aaa(1,:))/300;%做300次试验的均值
disp('E_source_doa');
RMSE_source_doa=sqrt(sum((aaa(1,:)-source_doa).^2)/300);%做300次试验的均方根误差
disp('RMSE_source_doa');
disp(RMSE_source_doa);
bbb(:,kk)=RMSE_source_doa;
end
disp(bbb);
plot(snr,bbb(1,:),'k*-');
save CLASSICAL_MUSIC_snr_rmse.mat;
legend('经典MUSIC');
xlabel('信噪比（snr）/dB');
ylabel('估计均方根误差');
grid on;
