基于经典MUSIC的DOA估计Matlab仿真.doc

信息与通信工程学院 阵列信号处理实验报告 （基于经典MUSIC的DOA估计Matlab仿真）学 号：XXXXXX专 业：XXXXXX学生姓名：XXX任课教师：XXX2015年X月题目：基于经典MUSIC的DOA估计Matlab仿真1. 算法简述：基于天线阵列协方差矩阵的特征分解类DOA估计算法中，多重信号分类（MUSIC）算法具有普遍适用性，只要已知天线阵的分布形式，无论直线阵还是圆阵，不管阵元是否等间隔分布，都可以得到高分辨率的估计结果。阵列协方差矩阵R可以划分为两个空间，即。因导向矩阵的各矢量与噪声子空间正交，可得到阵列空间谱函数使变化，按照来搜索峰值来达到波达方向的估计值。2. 实验内容与结果：实验使用8阵元均匀线阵，阵元间距为信号波长的一半，输入信号为1个BPSK信号，1个非相干的单频干扰，设置载波频率10MHz、采样频率30MHz、快拍数30、信干比0dB、信号方位角、干扰方位角，分析信噪比从1到20dB估计均方误差。实验结果见下图。图1信噪比 15dB情况下的波束图图2不同信噪比测得的均方误差3. 仿真分析由仿真结果可知，在实验前提条件下，经典music算法的DOA估计可有效估计出信号的来向。当信噪比较小时，估计得均方误差较大；随着信噪比的提高，估计均方误差逐渐减小；当信噪比大于10dB以后，角度估计成功率100%。MUSIC算法在DOA估计中属于高精度算法，而算法的估计精度仍受阵元数量、信号角度、信号数量、信号间距、信噪比等影响，本实验仿真分析了一定情况下估计精度与信噪比的关系，得到了预期的实验结果。4. 程序clear all;close all;clc;%-%线阵 music DOA%-%参数设置Signal_No = 1;Interference_No = 1;S_No = Signal_No + Interference_No;sensor_No = 8;azimuth = 0*pi/180 5*pi/180;Fs = 30e6;Fc = 10e6;F1 = 8e6; %干扰功率theta = 10 45 52 93 43*pi/180;RB = 2e6;M = 2;%二进制wavelength = 3e8/Fc;d = wavelength/2;K = 300; %快拍数%-%dBPs = -30;SNR = 10;SIR = 00;M_No = 40; %码速率等于信息速率Data_No = M_No*Fs/RB;t = 1/Fs:1/Fs:Data_No/Fs;%-% for xunhuan=1:1% mean1 = 0;% mean2 = 0;% for SNR =1:20%-%功率转换Ps_l = 10(Ps/10);Pi_l = 10(Ps-SIR)/10);Pn_l = 10(Ps-SNR)/10);%-%信号生成bit = randint(1,M_No); %产生信息序列bitstream = ;for i = 1:M_No if bit(i) = 1 bitstream = bitstream, ones(1, Fs/RB); else bitstream = bitstream, -ones(1, Fs/RB); endendCarrier_R = cos(2*pi*Fc*t);S_R = Carrier_R .* bitstream;Carrier_I = sin(2*pi*Fc*t);S_I = Carrier_I .* bitstream;Signal_R = sqrt(Ps_l)*S_R/sqrt(S_R*S_R/length(S_R);Signal_I = sqrt(Ps_l)*S_I/sqrt(S_I*S_I/length(S_I);%生成BPSK复信号Signal(1,:) = complex(Signal_R, Signal_I);%-%产生干扰%I_R = sqrt(2*Pi_l)*cos(2*pi*F1*t); %I_I = sqrt(2*Pi_l)*sin(2*pi*F1*t);for i = 2:Interference_No + 1 Signal(i,:) = sqrt(2*Pi_l)*complex(sin(2*pi*F1(i-1)*t+theta(i), sin(2*pi*F1(i-1)*t+theta(i);end%-%模拟天线接收As = zeros(sensor_No,S_No);for i = 1:sensor_No for ii = 1:S_No As(i,ii) = exp(-1i*2*pi*sin(azimuth(ii)*d/wavelength*(i-1); endendx = As * Signal;%-%加噪声noise_R = randn(sensor_No,Data_No)*sqrt(Pn_l);noise_I = randn(sensor_No,Data_No)*sqrt(Pn_l);noise = complex(noise_R, noise_I);s = x + noise;x1 = s(:,(1:K);R=x1*x1/K;V,D=eig(R);D1=diag(D);Un=V(:,1:sensor_No-S_No);Gn=Un*Un;searching_doa=-90:0.1:90;%线阵的搜索范围为-9090度 for i=1:length(searching_doa)a_theta=exp(-j*(0:sensor_No-1)*2*pi*d*sin(pi*searching_doa(i)/180)/wavelength); Pmusic(i)=1./abs(a_theta)*Gn*a_theta); endplot(searching_doa,10*log(Pmusic),r);xlabel(入射角/degree);ylabel(空间谱/dB);legend(Music Spectrum);title(经典MUSIC估计);grid on;%-%查最大值% Pmax_db1=0;% for i=1:909% if Pmusic(i)Pmax_db1% Pmax_db1=Pmusic(i);% Pmax_pitching1=i;% end% end% Pmax_db2=0;% for i=910:1801% if Pmusic(i)Pmax_db2% Pmax_db2=Pmusic(i);% Pmax_pitching2=i;% end% end% err1(SNR) = (Pmax_pitching1-1)/10-90)2;% err2(SNR) = (Pmax_pitching2-1)/10-90-5)2;% end% mean1 = mean1 + err1;% mean2 = mean2 + err2;% end% mse1 = mean1/1000;% mse2 = mean2/1000;% subplot(211); plot(mse1);xl