基于cyclic-MUSIC的DOA估计Matlab仿真

信息与通信工程学院阵列信号处理实验报告（基于cyclic-MUSIC的DOA估计Matlab仿真）学号：XXXXXX专业：XXXXXX学生姓名：XXX任课教师：XXX2015年X月题目：基于cyclic-MUSIC的DOA估计Matlab仿真1.算法简述：DOA估计算法一般使用与平稳信号，而生活中测量的信号往往是非平稳的。在非平稳的信号中，有一类信号属于周期平稳信号，利用周期平稳性得到的算法称循环（共轭循环）平稳算法。循环平稳的循环自相关函数表示为1*j2πn1()()()0NxxnxnxneNR　　其中为循环角频率。共轭循环平稳的循环自相关函数表示为1j2πn1()()()0NxxnxnxneNR　　基于天线阵列协方差矩阵的特征分类DOA估计算法中，MUSIC算法具有普遍的实用性。只要已知天线阵的分布形式，无论直线阵还是圆阵，不管阵元是否等间隔分布，都可以得到高分辨率的估计结果。2.实验内容与结果：实验先仿真分析BPSK信号的循环平稳和共轭循环平稳性与循环频率和延时的的关系，然后使用12阵元均匀线阵，阵元间距为信号波长的一半，输入信号为1个BPSK信号、设置载波频率10MHz、采样频率50MHz、快拍数1000、信号方位角-40、信噪比为25dB情况下，得到DOA估计结果。020406080-10-50510012延时循环频率循环频率模图1BPSK的循环平稳性020406080-30-20-100102030012延时循环频率循环频率模图2BPSK共轭循环平稳性-100-80-60-40-20020406080100-40-30-20-100102030入射角/degree空间谱/dBcyclicMUSIC估计MusicSpectrum图3DOA仿真结果3.仿真分析由仿真结果可知，BPSK信号具有循环平稳特性，当循环频率为bkf，信号具有循环平稳性，当循环频率为bckfmf，信号具有共轭循环平稳性。由图一图二可得到合适的循环频率和延时，例如图一，循环频率为bf、延时为40Ts时其循环自相关函数值最大。实验最后针对cyclic-music算法的DOA估计，使用共轭循环平稳性设计算法，并在给定实验条件下仿真，有效估计出了信号的来向。4.程序1.程序1：clearall;closeall;clc;%---------------------------------%线阵波束形成%---------------------------------%参数设置S_No=1;sensor_No=12;azimuth=[40*pi/180];Fs=40e6;Fc=10e6;%干扰功率RB=1e6;M=2;%二进制wavelength=3e8/Fc;d=wavelength/2;K=1000;%快拍数%---------------%dBPs=-0;SNR=20;M_No=40;%码速率等于信息速率Data_No=M_No*Fs/RB;t=1/Fs:1/Fs:Data_No/Fs;%---------------------------------%功率转换Ps_l=10^(Ps/10);Pn_l=10^((Ps-SNR)/10);%---------------------------------%信号生成bit=randint(1,M_No);%产生信息序列bitstream=[];fori=1:M_Noifbit(i)==1bitstream=[bitstream,ones(1,Fs/RB)];elsebitstream=[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);%---------------------------------%{%---------------------------------%模拟天线接收As=zeros(sensor_No,S_No);fori=1:sensor_Noforii=1:S_NoAs(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;%}X=Signal(1,(1:K));Nm=length(X(1,:));%产生共轭循环矩阵flagg=1;Num1=100;Num2=5;forn=1:Num1forKm=0:39Fcj=(n-Num1/2)/Num2*1e6;%选择循环频率x1=X(1:Nm-Km)-mean(X(1:Nm-Km));x2=(X(1+Km:Nm)-mean(X(1+Km:Nm))).*exp(-1i*2*pi*Fcj/Fs*(1:Nm-Km));ifflagg==1xx=x1*x2';%循环平稳算法xn=xx/(Nm-Km);elseifflagg==2xx=x1*conj(x2)';%共轭循环平算法xn=xx/(Nm-Km);endphas=exp(1i*2*pi*Fcj/Fs*Km);R(n,Km+1)=xn*phas;endend%Max=0;%form1=1:Num1%form2=1:50%ifR(m1,m2)=Max;%m3=m1;m4=m2;%Max=R(m1,m2);%end%end%end%m3=m3-Num1/2%m4aaa=-(Num1/2-1)/Num2:1/Num2:Num1/2/Num2;bbb=1:2:80;Rr=abs(R);meshc(bbb,aaa,Rr);ylabel('循环频率');xlabel('延时');zlabel('循环频率模');2.程序2：clearall;closeall;clc;%---------------------------------%线阵musicDOA%---------------------------------%参数设置Signal_No=1;Interference_No=1;S_No=Signal_No+Interference_No;sensor_No=8;azimuth=[0*pi/1805*pi/180];Fs=30e6;Fc=10e6;F1=[8e6];%干扰功率theta=[1045529343]*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;%---------------------------------%forxunhuan=1:1%mean1=0;%mean2=0;%forSNR=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=[];fori=1:M_Noifbit(i)==1bitstream=[bitstream,ones(1,Fs/RB)];elsebitstream=[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);fori=2:Interference_No+1Signal(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);fori=1:sensor_Noforii=1:S_NoAs(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;%线阵的搜索范围为-90~90度fori=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('MusicSpectrum');title('经典MUSIC估计');gridon;%%---------------------------------------------%%查最大值%%---------------------------------------------%%Pmax_db1=0;%fori=1:909%ifPmusic(i)Pmax_db1%Pmax_db1=Pmusic(i);%Pmax_pitching1=i;%end%end%Pmax_db2=0;%fori=910:1801%ifP