阵列信号处理MUSIC法扩展复现

function [error] = softdomain_perform(SNR,snapshot,theta) %{ 函数从信噪比，快拍数，信源间隔三个方面对MUSIC法进行性能测试,参数定义为 SNR--------信噪比 snapshot---快拍数 theta------信源间隔，输入为行矢量 error------来向角估计平均错误角度 %} %% ===========================参数设置================================ error=zeros(2,1); %初始化 f0=1e10; %信号频率为X波段10GHz fs=3*f0; %采样频率 c=3e8; %光速 lamda=c/f0; %求波长lamda d=lamda/2; %阵元间隔 M=16; %阵元数目 L = 4; %空间平滑次数 K = M-L+1; %子阵阵元个数 theta=theta'; %方向来向角theta=30度 theta_p=-90:0.01:90; %方向角度扫描从-90到90 %% ============================信号源生成============================== atheta=exp(-j*2*pi*d*sin(theta/180*pi)*(0:M-1)/lamda);%a(theta)为空域导向矢量 % s1=sin(2*pi*f0*(0:1/fs:(snapshot-1)/fs)); %假设发射信号为sin(t)函数，对其采样1000次 % s2=cos(2*pi*f0*(0:1/fs:(snapshot-1)/fs)); %两个信源另一个可设为cos函数 s1 = exp(j*2*pi*f0*(0:1/fs:(snapshot-1)/fs)); s2 = exp(j*2*pi*f0*(0:1/fs:(snapshot-1)/fs))*exp(j*(pi/6)); S=[s1;s2]; %% ============================接收信号============================== X=atheta'*S;%采样信号为X=a(theta)*singal+n(t) X=awgn(X,SNR);%添加白高斯噪声，snr=20 %% ========================MUSIC法求空间谱============================== Rhat=X*X'/snapshot;%R^=X*X'/N [sigma,lamdaM_temp] = eig(Rhat); [lamdaM,index] = sort(diag(lamdaM_temp),'descend'); sigma = sigma(:,index); Us = sigma(:,1:2);%Us为信源特征向量 Un = sigma(:,3:M);%Un为噪声特征向量 %% =======================空间平滑法================================= allx=zeros(K,snapshot,L); Rhat2=zeros(K,K,L); RL=zeros(K,K); for i=1:L allx(:,:,i)=X(i:i+K-1,:); Rhat2(:,:,i)=allx(:,:,i)*allx(:,:,i)'/snapshot; RL=Rhat2(:,:,i)+RL; end RL=RL/L; %求RL的特征值和特征向量，sigma为从大到小排列的特征向量，lamdaK为K个特征值 [sigma_soft,lamdaK_temp]=eig(RL); [lamdaK,index2]=sort(diag(lamdaK_temp),'descend'); sigma_soft=sigma_soft(:,index2); UN_soft=sigma_soft(:,1:3);%将信号子空间提取出来 %% ==========================求空间谱================================= for i = 1:length(theta_p)%从-90到90度扫描，每一度都对信号相位补偿求功率 %直接MUSIC法 atheta_p1=exp(-j*2*pi*d*sin(theta_p(i)/180*pi)*(0:M-1)/lamda);%每一度的补偿相位 p(i)=abs((atheta_p1*atheta_p1')/(atheta_p1*(eye(M)-Us*Us')*atheta_p1'));%求功率 %空间平滑MUSIC法 atheta_p2=exp(-j*2*pi*d*sin(theta_p(i)/180*pi)*(0:K-1)/lamda);%每一度的补偿相位 p_soft(i)=abs((atheta_p2*atheta_p2')/(atheta_p2*(eye(K)-UN_soft*UN_soft')*atheta_p2'));%求功率 end %% ==========================来向角估计=============================== P3=10*log10(p/max(p));%归一化功率 P_soft = 10*log10(p_soft/max(p_soft)); [P3peak,peakpos3]=findpeaks(P3);%寻找空间谱峰值 [Ppeak_soft,peakpos_soft]=findpeaks(P_soft); [~,pos3]=sort(P3peak,'descend');%对峰值进行排序 [~,pos_soft]=sort(Ppeak_soft,'descend'); theta3=pos3(1:2);%取最大的两个谱峰进行来向角估计，记录索引 theta_soft=pos_soft(1:2); theta3=peakpos3(theta3);%找到空间谱中最大峰值的位置 theta_soft=peakpos_soft(theta_soft); theta3=theta_p(theta3);%得到来向角估计 [theta3,~]=sort(theta3,'ascend'); theta_soft=theta_p(theta_soft); [theta_soft,~]=sort(theta_soft,'ascend'); error(1) = sqrt(((theta(1)-min(theta3))^2 + (theta(2)-max(theta3))^2)/length(theta));%求误差 error(2) = sqrt(((theta_soft(1)-min(theta_soft))^2 + (theta_soft(2)-max(theta_soft))^2)/length(theta_soft)); end