如何利用matlab修改已有程序，然后绘制图中的所给图像？ |
您所在的位置：网站首页 › Matlab绘制玫瑰图 › 如何利用matlab修改已有程序，然后绘制图中的所给图像？ |
如何利用matlab修改已有程序，然后绘制图中的所给图像？
|
问题：如何利用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; |
【本文地址】
今日新闻 |
推荐新闻 |