%--------------------------------------------%相干信號DOA估計(jì)MD（矩陣分解）算法 %適用于均勻線陣%修正后的music算法%程序分為信源、全局變量定義、信號產(chǎn)生、預(yù)處理算法、DF算法、做圖等幾部分%盒子 于200803~200805期間研學(xué)，無版權(quán)%--------------------------------------------clc;%1：信源%--------------------------------------------Pd=2000;Fd=1;Fs=4*Fd;R=0.5;Delay=5;M=4;x1=randint(Pd,1,M); %生成一個隨機(jī)數(shù)為2000*1的矩陣，輸出范圍為0~3x2=randint(Pd,1,M);y1=modmap(x1,Fd,Fs,'qask',M);%數(shù)字信號調(diào)制模擬信號，modmap可用scatterplot代替y2=modmap(x2,Fd,Fs,'qask',M);[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);%對輸入信號進(jìn)行升余弦濾波[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);s1=amodce(rcv_a1,10,'qam');%輸出復(fù)包絡(luò)，amodce可用ammod來代替,s1為8040*1% s2=amodce(rcv_a2,10,'qam');s2=s1; %假設(shè)信號s1和s2相關(guān)save sig3 s1 s2 %把變量s1,s2保存為sig3.mat文件%--------------------------------------------%2：全局變量定義%--------------------------------------------% clearm=8; %陣元數(shù)angle1=-50;%信號角度angle2=60;th=[angle1;angle2];%th為2*1矩陣nn=1024;%采樣數(shù)SN1=10;%信號的SNRSN2=10;SN3=12;sn=[SN1;SN2];%sn為2*1矩陣degrad=pi/180;%--------------------------------------------%3：生成離散信號%--------------------------------------------load sig3 %把sig3.mat文件中的s1,s2變量裝入內(nèi)存tt=1:nn;S=[s1(tt).';s2(tt).'];%S為2*1024矩陣nr=randn(m,nn);%nr為8*1024矩陣,為正態(tài)分布隨機(jī)陣ni=randn(m,nn);U=nr+j*ni; %構(gòu)造噪聲源，U為8*1024復(fù)數(shù)矩陣Ps=S*S'/nn;%信號的協(xié)方差矩陣，Ps為2*2矩陣ps=diag(Ps);%取主對角線上值，ps為2*1矩陣refp=2*10.^(sn/10);tmp=sqrt(refp./ps);%tmp為2*1矩陣S2=diag(tmp)*S;%S2為2*1024矩陣%--------------------------------------------%4：協(xié)方差矩陣及前后向平滑%--------------------------------------------tmp2=[0:m-1]'; %陣元位置A=exp(-i*pi*tmp2*sin(th'*degrad)); %A為8*2矩陣，方向矩陣由公式得出X=A*S2+U;%X為8*1024矩陣，1024次采集后得到的數(shù)據(jù)，陣列接受信號Rxx=X*X'/nn;%Rxx為8*8矩陣，陣列接受信號協(xié)方差a=eye(8);b=a(:,8:-1:1);%構(gòu)造交換矩陣RXX=(Rxx+b*Rxx.'*b)/2;%--------------------------------------------%5：奇異值分解%--------------------------------------------[Q,SS,W]=svd(RXX); Vu=Q(:,3:8 ); %小特征值對應(yīng)的特征向量，即噪聲子空間%--------------------------------------------%6：做圖%--------------------------------------------th2=[-90:90]';tmp2=[0:m-1]';a2=-i*tmp2*pi*sin(th2'*degrad);A2=exp(a2);num=diag(A2'*A2);Ena=Vu'*A2;den=diag(Ena'*Ena);doa=num./den;semilogy(th2,abs(doa));axis([-90 90 0.1 1e5]);grid on;hold on;%--------------------------------------------