关于MUSIC算法的改进算法，不需要已知信源数，科研取得较...

admin

1. %%%利用盖式圆半径法进行信源数目的估计
   
2. % close  all;
   
3. clear all;clc  
   
4. f=[420,420]*10^6;
   
5. fs=2*420*10^6;
   
6. c=3e8;
   
7. numda=c/max(f); %波长
   
8. d=numda/(2);
   
9. N=516;%%快拍数
   
10. M=5;%%阵元数
    
11. P=2; %%信源数
    
12. doa=[15,30]*pi/180;%来波方向
    
13. SNR=10;
    
14. coef=10^(SNR/20);
    
15. % A=rand(N,P);
    
16. n=[1:M]';
    
17. A=exp(-j*2*pi*d*(n-1)*sin(doa)/numda);
    
18. X=rand(M,N);
    
19. n0=1:N;
    
20. X=A*coef*exp(-j*(2*pi*f'*n0/fs+randn(P,N)))+randn(M,N);
    
21. % X=A*coef*sin(-j*2*pi*f'*n0/fs+randn(P,M))+randn(N,M);
    
22. R=1/N*X*X';
    
23. % load E:\xinyuan_guji\终板\x_ref\极低\x_ref5.dat;
    
24. % load E:\xinyuan_guji\终板\x_ref\极低\x_ref6.dat;
    
25. % load E:\xinyuan_guji\终板\x_ref\极低\x_ref7.dat;
    
26. % load E:\xinyuan_guji\终板\x_ref\极低\x_ref8.dat;
    
27. %
    
28. % X=[x_ref5.';x_ref6.';x_ref7.';x_ref8.'];
    
29. % M=4;
    
30. % data=zeros(4,(length(x_ref5))/2);
    
31. %  for n=1:4
    
32. %     for k=1:(length(X(1,:)))/2
    
33. %             Idata(k)=X(1,2*k-1);
    
34. %             Qdata(k)=data(2*k);
    
35. %         data(n,k)=X(n,2*k)*sqrt(-1)+X(n,2*k-1);
    
36. %     end
    
37. %  end
    
38. %
    
39. % data_length=length(data(1,:));
    
40. % R=data*data'/data_length;
    
41.   
    
42. %%%协方差矩阵分块
    
43. R1=R(1:M-1,1:M-1);
    
44. [U1,D1]=eig(R1);%%R1特征分解
    
45. DD=diag(D1);
    
46. oo=zeros(length(U1(:,1)),1);
    
47. U=[U1 oo;oo' 1];
    
48. S=U'*R*U;
    
49. Cmm=sum(S(1:M-1,M));
    
50. 
    
51. dd=[sqrt(Cmm.*DD);1];
    
52. G=diag(dd);
    
53. SS=inv(G)*S*G;
    
54. 
    
55. r=zeros(M,1);
    
56. r1=zeros(M,1);
    
57. 
    
58. r(1:M-1,1)=S(1:M-1,M);
    
59. r1(1:M-1,1)=SS(1:M-1,M);
    
60. 
    
61. %r(M)=sum(S(1:M-1,M));
    
62. r=abs(r);
    
63. r1=abs(r1);
    
64. % D_N=0.2;
    
65. D_N=4./log(512);
    
66. K=0;
    
67. K1=0;
    
68. K2=0;
    
69. for k=1:M-1
    
70.      GDE=r(k)-(D_N/(M-1))*sum(r(1:M-1));
    
71.      if GDE<=0
    
72.       K=k-1; break;
    
73.      end   
    
74. end
    
75. K
    
76. 
    
77. 
    
78.   clear k;
    
79. for k=1:M-1
    
80.      GDE1=r1(k)-(D_N/(M-1))*sum(r1(1:M-1));
    
81.      if GDE1<=0
    
82.       K1=k-1; break;
    
83.      end   
    
84.   end
    
85. K1
    
86. 
    
87. L=16;
    
88. for nn=1:L
    
89.     for k=1:M-1
    
90.         GDE=r(k)-(D_N/(M-1))*sum(r(1:M-1));
    
91.         if GDE<=0
    
92.             K(nn)=k-1;  break;           
    
93.         end
    
94.         
    
95.      end
    
96. end
    
97. 
    
98. K2=mean(K)
    
99. 
    
100. 
     
101. %%%各阵元输出信号的协方差矩阵的特征值分解
     
102.   [V,D]=eig(R);
     
103.   d1=diag(D);
     
104.   lambtal=d1;
     
105. %   [lambtal,p]=sort(d1);%%lambtal按升序排列   k输出的是每个值所在的行后或列
     
106. %   lambtal=(fliplr(lambtal'))';
     
107. %   V=V(:,p);
     
108. % %   z=fliplr(z);
     
109.          
     
110.       
     
111.       %%%产生空间谱函数
     
112.     VV=V(:,K+1:M);
     
113. %       en=en/norm(en);%%归一化
     
114. angel=[0:0.1:90];
     
115. angel1=[0:0.1:90]*pi/180;
     
116. alfa=exp(-1i*2*pi/numda*d*(n-1)*sin(angel1));
     
117. pum=1./diag((alfa'*VV*VV'*alfa));
     
118. figure(1);
     
119. plot(angel,10*log10(abs(pum)));
     
120. grid on
     
121. title('music方法')
     
122. xlabel('DOA(deg)')
     
123. ylabel('amplitude(dB)')
     
124. %%源数目的估计
     
125. 
     
126. 
     
127. %%%产生空间谱函数
     
128.     V1=V(:,K1+1:M);
     
129. %       en=en/norm(en);%%归一化
     
130. angel=[0:0.1:90];
     
131. angel1=[0:0.1:90]*pi/180;
     
132. alfa=exp(-1i*2*pi/numda*d*(n-1)*sin(angel1));
     
133. pum1=1./diag((alfa'*V1*V1'*alfa));
     
134. figure(2);
     
135. plot(angel,10*log10(abs(pum1)));
     
136. grid on
     
137. title('music方法')
     
138. xlabel('DOA(deg)')
     
139. ylabel('amplitude(dB)')
     
140. %%源数目的估计
     
141. 
     
142. 
     
143. 
     
144. 
     
145. 
     
146. 
     
147. 
     
148. 
     
149. 
     
150. 
     
151. 
     
152. 
     
153. 
     
154. 
     
155. 
     
156. 
     

复制代码